Fix the GBM model and the price simulation

finance/models.ts

@@ -1,16 +1,30 @@
 /**
- * μ = sample mean
- * σ = sample volatility
- * Δt = 1 (1 day)
- * φ = normally distributed random number
+ * μ = sample mean (drift per unit time)
+ * σ = sample volatility (per √(unit time))
+ * Δt = time step
+ * φ = standard normal random number
+ *
+ * Log return over Δt: (μ - σ² / 2)Δt + σ√(Δt)φ
+ * so that S(t + Δt) = S(t) * exp(log return) and E[S(t + Δt)] = S(t) * exp(μΔt).
 */
 export const geometric_brownian = (estimators, time, random) => {
     let miu = estimators.mean;
     let sigma = estimators.std;
 
-    let drift = miu * time;
-    let shock = sigma * random * time;
-    console.log(`miu = ${miu}, sigma = ${sigma}, drift = ${drift}, shock = ${shock}`);
+    let drift = (miu - (sigma * sigma) / 2) * time;
+    let shock = sigma * Math.sqrt(time) * random;
     let result = drift + shock;
     return result;
-};
+};
+
+/**
+ * Box-Muller transform: turns two uniforms on (0, 1) into a standard normal.
+ * Math.random alone is uniform on [0, 1) and would only ever shock the price up.
+*/
+export const standard_normal = () => {
+    let u = 0;
+    let v = 0;
+    while (u === 0) { u = Math.random(); }
+    while (v === 0) { v = Math.random(); }
+    return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
+};

finance/stock.ts

@@ -1,7 +1,7 @@
 import { historical, snapshot  } from 'yahoo-finance'
 import { mean, std  } from 'mathjs'
 import { run  } from './../src/montecarlo/index'
-import { geometric_brownian  } from './models'
+import { geometric_brownian, standard_normal  } from './models'
 
 const one = 1
 const thousand = 1000
@@ -41,15 +41,17 @@ export class StockHelper {
     return quotes;
   }
 
-  static calculate_price(sample_size)
+  // run() calls the estimator with a single argument, so the price, the
+  // estimators and the horizon are closed over instead of being passed through.
+  static calculate_price(price, estimators, time, sample_size?)
  {
     if (sample_size == undefined)
     {
       sample_size = thousand;
     }
 
-    let generator = Math.random;
-    let estimator_function = geometric_brownian;
+    let generator = standard_normal;
+    let estimator_function = (random) => price * Math.exp(geometric_brownian(estimators, time, random));
     return run(sample_size, generator, estimator_function);
   }
 

finance/test/models.test.ts

@@ -0,0 +1,82 @@
+import { geometric_brownian, standard_normal } from '../models'
+import { StockHelper } from '../stock'
+
+const sample = (generator, size) => {
+    let samples = [];
+    for (let i = 0; i < size; i++) {
+        samples.push(generator());
+    }
+
+    return samples;
+};
+
+describe('standard_normal: ', () =>
+{
+    it('has mean zero and unit variance', () =>
+    {
+        let samples = sample(standard_normal, 200000);
+        let stats = StockHelper.getStats(samples);
+
+        expect(stats['mean']).toBeCloseTo(0, 1);
+        expect(stats['std']).toBeCloseTo(1, 1);
+    });
+
+    // Math.random is uniform on [0, 1), so the shock could only ever be positive.
+    it('produces negative draws', () =>
+    {
+        let samples = sample(standard_normal, 1000);
+
+        expect(samples.filter((draw) => draw < 0).length).toBeGreaterThan(300);
+    });
+});
+
+describe('geometric_brownian: ', () =>
+{
+    let estimators = { mean: 0.1, std: 0.4 };
+
+    it('applies the ito correction to the drift', () =>
+    {
+        // With no shock the log return is (mu - sigma^2 / 2) * time.
+        let time = 2;
+        let expected = (0.1 - 0.4 * 0.4 / 2) * time;
+
+        expect(geometric_brownian(estimators, time, 0)).toBeCloseTo(expected, 10);
+    });
+
+    it('scales the shock by the square root of time', () =>
+    {
+        // Doubling the horizon must scale the shock by sqrt(2), not by 2.
+        let drift_one = geometric_brownian(estimators, 1, 0);
+        let drift_two = geometric_brownian(estimators, 2, 0);
+
+        let shock_one = geometric_brownian(estimators, 1, 1) - drift_one;
+        let shock_two = geometric_brownian(estimators, 2, 1) - drift_two;
+
+        expect(shock_two / shock_one).toBeCloseTo(Math.sqrt(2), 10);
+    });
+});
+
+describe('StockHelper.calculate_price: ', () =>
+{
+    let price = 100;
+    let estimators = { mean: 0.08, std: 0.25 };
+    let time = 2;
+
+    it('converges to the analytical mean of the lognormal price', () =>
+    {
+        // E[S_T] = S_0 * exp(mu * T)
+        let expected = price * Math.exp(estimators.mean * time);
+        let estimated = StockHelper.calculate_price(price, estimators, time, 400000);
+
+        expect(estimated).not.toBeNaN();
+        expect(estimated / expected).toBeCloseTo(1, 1);
+    });
+
+    it('defaults the sample size and stays finite', () =>
+    {
+        let estimated = StockHelper.calculate_price(price, estimators, time);
+
+        expect(Number.isFinite(estimated)).toBe(true);
+        expect(estimated).toBeGreaterThan(0);
+    });
+});