Skip to content
All library documents

Estimating GBM Volatility from Simulated Log Returns

Article Quant Q&A · Author: Arthur

Summary

The document addresses why volatility estimated from simulated geometric Brownian motion prices may not match the sigma parameter used to generate them. It points to three issues in the example: the calculated log returns are per time step, only a subset of simulated paths is processed, and the standard deviation uses a population divisor rather than a sample divisor.

For a GBM with annualized volatility σ and step length Δt, log-return standard deviation per step is σ√Δt; dividing the estimate by √Δt annualizes it. The answer recommends processing every path and using the sample standard deviation when estimating from a finite sample. Its code illustrates those adjustments, but the written answer reverses the per-step scaling relationship by stating σ/√Δt before recommending division by √Δt. That inconsistency matters: the latter operation is appropriate for annualizing per-step returns. Estimates from finite simulations will still vary around the input parameter.

Key ideas

  • GBM log returns over a time step of length Δt have standard deviation σ√Δt.
  • To estimate annualized volatility from per-step returns, divide their standard deviation by √Δt.
  • The example processes only a few paths despite allocating storage for many, leaving other results at zero.
  • A sample standard deviation uses n−1 in the denominator, though finite samples remain noisy.
  • The document’s prose gives an inconsistent per-step scaling formula, while its code applies the annualization adjustment.

Tags

Full text
# Volatility is not becoming sigma in simulated GBM


# Volatility is not becoming sigma in simulated GBM












I am trying to simulate GBM using values for sigma and then after calculating the price, calculate the log of returns and the associated volatility. I was expecting the calculated volatility of log returns to be equal to sigma used in the simulation, but its not the case. I would like to know if this result was expected and if so, what is wrong with my code that is preventing me from getting to the right answer.

```
import random
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

M=1000
n=126
T=0.5
r=0.1175
sigma=0.2
S0=100
dt = T/n
np.random.seed(100)

mudt = r - 1/2*sigma**2
St=np.zeros((M,(n+1)))
St[:,0]=S0
for i in range(M):
    for j in range(n):
        St[i,j+1] = St[i,j]*np.exp(mudt*dt+sigma*np.random.normal()*dt**(0.5))

for i in range(M):
    plt.plot(St[i,:])

Ravg=np.zeros(M)
R=np.zeros((M,(n)))
for i in range(5):
    RR = pd.Series(St[i])
    RR = RR.pct_change()
    RR = np.log(RR+1)
    Ravg[i] = RR.mean()
    RR=RR.dropna()
    R[i] = RR.values
    
   
var=np.zeros(M)
for i in range(M):
    vari=0
    for j in range(n):
        vari += (R[i,j]-Ravg[i])**2
    var[i]=(vari/n)**(1/2)

```
```

## Answer by Achrbot (score 3)

https://quant.stackexchange.com/a/77819

Under GBM, the log returns $\log S_{t+dt} - \log S_t$ have a standard deviation of $\frac{\sigma}{\sqrt{dt}}$. So you need to scale your computed standard deviation by $1/\sqrt{dt}$.

You only compute the returns for range(5), and leave the rest as zeros - So you get the wrong result.

Finally, you are using a biased estimator for the standard deviation. The divisor should be $1/(n-1)$, not $1/n$.

I have amended your code, and I get the expected result.

```
import random
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

M=1000
n=126
T=0.5
r=0.1175
sigma=0.2
S0=100
dt = T/n
np.random.seed(100)

mudt = r - 1/2*sigma**2
St=np.zeros((M,(n+1)))
St[:,0]=S0
for i in range(M):
    for j in range(n):
        St[i,j+1] = St[i,j]*np.exp(mudt*dt+sigma*np.random.normal()*dt**(0.5))

Ravg=np.zeros(M)
R=np.zeros((M,(n)))
for i in range(M):
    RR = pd.Series(St[i])
    RR = RR.pct_change()
    RR = np.log(RR+1)
    Ravg[i] = RR.mean()
    RR=RR.dropna()
    R[i] = RR.values
    
   
var=np.zeros(M)
for i in range(M):
    vari=0
    for j in range(n):
        vari += (R[i,j]-Ravg[i])**2
    var[i]=(vari/(n-1))

std = var**0.5
np.mean(std)/(dt**0.5)
```
```

Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)

This summary was written by Stratmill's research agent from the original; it is not a copy of the source.