How to estimate the parameters of a geometric Brownian motion (GBM)? It seems rather simple but actually took me quite some time to solve it. The most intuitive way is by using the method of moments.

Estimation of ABM

First let us consider a simpler case, an arithmetic Brownian motion (ABM). The evolution is given by

$$ dS = \mu dt + \sigma dW. $$

By integrating both sides over $(t,t+T]$ we have

$$ \Delta \equiv S(t+T) - S(t) = \left(\mu - \frac{\sigma^2}{2}\right) T + \sigma W(T) $$

which follows a normal distribution with mean $(\mu - \sigma^2/2)T$ and variance $\sigma^2 T$. That is to say, given $T$ and i.i.d. observations $\\{\Delta_1,\Delta_2,\ldots,\Delta_n\\}$ for different $t$ values¹, with sample mean

and modified sample variance

we have unbiased estimator for $\mu$

and for $\sigma^2$ we have

$$ \hat{\sigma}^2 = \frac{\hat{\sigma}_{\Delta}^2}{T}. $$

Now we prove the consistency. First we consider the variance of $\hat{\mu}_{\Delta}$

and the variance of $\hat{\sigma}_{\Delta}^2$

The variance of $\hat{\mu}$ is therefore given by

and the variance of $\hat{\sigma}^2$ is given by

$$ Var(\hat{\sigma}^2) = \frac{Var(\hat{\sigma}_{\Delta}^2)}{T^2} = \frac{2\sigma^4}{n}. $$

So the two estimators are also both consistent. It should be noticed that there exists certain “trade-off” between the efficiency of $\hat{\mu}_{\Delta}$ and $\hat{\sigma}_{\Delta}^2$ by varying the value of $T$.

Estimation of GBM

For a general GBM with drift $\mu$ and diffusion $\sigma$, we have PDE

$$ \frac{dS}{S} = \mu dt + \sigma dW, $$

so we can integrate² the both sides within $(t,t+T]$ for any $t$ and get

$$ \Delta \equiv \ln S(t+T) - \ln S(t) = \left(\mu - \frac{\sigma^2}{2}\right) T + \sigma W(T). $$

The rest derivation is exactly the same.

Validation

Now we numerically validate this against monte Carlo simulation.

import numpy as np
import pandas as pd
from scipy.stats import ttest_1samp as ttest

np.random.seed(1)

timespan = 10000
mu = 0.002
sigma = 0.06
T = 50
S0 = 10
n_sim = 20000
sigma2 = sigma**2
mu_hat_list = []
sigma2_hat_list = []

for i in range(n_sim):
    # simulate log_S
    log_S0 = np.log(S0)
    log_St = mu - sigma2 / 2 + sigma * np.random.normal(size=timespan)
    log_S = np.cumsum(np.append(log_S0, log_St))
    # calculate delta
    delta = log_S[T::T] - log_S[:-T:T]
    n = len(delta)
    # estimate mu and sigma2
    mu_delta = delta.mean()
    sigma2_delta = ((delta - mu_delta)**2).sum() / (n - 1)
    mu_hat = (2 * mu_delta + sigma2_delta) / T / 2
    sigma2_hat = sigma2_delta / T
    mu_hat_list.append(mu_hat)
    sigma2_hat_list.append(sigma2_hat)

E_mu_hat_mc = np.mean(mu_hat_list)
Var_mu_hat_mc = np.var(mu_hat_list)
Var_mu_hat_mm = sigma2 * (2 + sigma2 * T) / n / T / 2
E_sigma2_hat_mc = np.mean(sigma2_hat_list)
Var_sigma2_hat_mc = np.var(sigma2_hat_list)
Var_sigma2_hat_mm = 2 * sigma2**2 / n

result = pd.DataFrame([
    [E_mu_hat_mc, mu, ttest(mu_hat_list, mu)[1]],
    [Var_mu_hat_mc, Var_mu_hat_mm, '-'],
    [E_sigma2_hat_mc, sigma2, ttest(sigma2_hat_list, sigma2)[1]],
    [Var_sigma2_hat_mc, Var_sigma2_hat_mm, '-']
], columns = ['monte Carlo', 'Method of moment', 'P Value'],
   index = ['E(mu_hat)', 'Var(mu_hat)', 'E(sigma2_hat)', 'Var(sigma2_hat)'])
result

Now we may safely apply this estimation in application.