Analysing Geometric Brownian Motion

Question

For uni I'm doing this exercise on brownian motion, specifically geometric brownian motion. For this part of the exercise I'm required to compute certain values: sample mean, sample variance, and the quadratic variation for m = 10, 100, 1000. The output I get is posted below but I'm not quite happy with the results. The distribution of the ends is a sort of shifted normal distribution do the variance and mean being off a little seems reasonable, but the thing I'm worried about is the quadratic variation, as this should be equal to T. T is set to 5 but the QV is nothing close to that. I cant see what I'm doing wrong, can anyone help me? any help is appreciated!

import yfinance as yf
import numpy as np
import math
from matplotlib import pyplot as plt
import scipy.stats as stats


class BMG:  #Geometric Brownian Motion
    def __init__(self, drift, variance_term, m, T):
        self.dt = T/m
        self.Z = np.random.standard_normal(m)
        self.Y = np.zeros(m+1)
        for i in range(1, m+1):
            self.Y[i] = math.e**(variance_term*math.sqrt(self.dt)*self.Z[i-1]) + drift*self.dt
        self.arr = np.ones(m+1)
        for i in range(m):
            self.arr[i+1] = self.arr[i]*self.Y[i+1]

        self.end = self.arr[-1]


fig, axs = plt.subplots(2)

N = 500                 #Number of BM simulations
drift = -0.1
variance_term = 0.25
m = 500                 #Number of time-steps 
T = 5                   #Time in years

ends = np.zeros(N)

for i in range(N):              
    bmg = BMG(drift, variance_term, m, T)
    ends[i]=bmg.end
    axs[0].plot(np.arange(0,m+1), bmg.arr)          #Plotting BM trajectories
axs[1].hist(ends, density = True, bins=int(N/5))    #Plotting distribution of end points

plt.show()



#Sample Mean
sample_mu = 0
for i in range(N):
    sample_mu += ends[i]
sample_mu = sample_mu/N
print("sample mean:  "+ str(sample_mu) + "  analytic mean:  " + str(bmg.arr[0]*math.e**(drift*T)))

#Sample Variance
sample_var = 0 
for i in range(N):
    sample_var += (ends[i] - sample_mu)**2
sample_var = sample_var/(N-1)
print("sample variance:  "+ str(sample_var) + 
"  analytic variance:  " + str( (bmg.arr[0]**2) * (math.e**(2*drift*T)) * ((math.e**(variance_term**2))-1)))

#Quadratic variation
steps = [10, 100, 1000]
for n in steps:
    bmg1 = BMG(-0.1, 0.25, n, T)
    QV = 0

    for i in range(n):
        step = bmg1.arr[i+1] - bmg1.arr[i]
        QV += step**2

    print("for m = " + str(n) +":   QV = " + str(QV))

/* Output:
sample mean:  0.6942405094555087  analytic mean:  0.6065306597126334
sample variance:  0.16157387382343052  analytic variance:  0.023726185505356632
for m = 10:   QV = 0.12398361878757773
for m = 100:   QV = 0.4954051008996326
for m = 1000:   QV = 0.10524633607091849
*/

Answer 1

I will provide here analytic results with which one can compare a (good) simulation. The SDE of a GBM is $dX_t=\mu X_tdt+\sigma X_t dW_t,\,X_0>0$. The solution is found by applying Ito to $f(x)=\ln x$ and we find $$\ln X_t=\underbrace{\ln X_0 + \bigg(\mu-\frac{1}{2}\sigma^2\bigg)t}_{:=d_t}+\sigma W_t \implies X_t=X_0e^{(\mu-\frac{1}{2}\sigma^2)t+\sigma W_t}$$ Since $X_t=e^{\ln X_t}$ and $\ln X_t \sim \mathcal{N}(d_t,\sigma^2t)$ we have that $X_t$ is lognormally distributed: for $x>0$ $$\begin{aligned}P(X_t\leq x)&=P(\ln X_t \leq \ln x)= \\&=\int_{(-\infty,\ln x]}\frac{1}{\sqrt{2\pi \sigma^2t}}\exp\bigg\{-\frac{(y-d_t)^2}{2\sigma^2 t}\bigg\}dy\end{aligned}$$ So by Leibniz integral rule and the fact that the Gaussian density decays at $y\downarrow -\infty$, $$f_{X_t}(x)=\frac{d}{dx}P(X_t\leq x)=\frac{1}{x}\frac{1}{\sqrt{2\pi \sigma^2t}}\exp\bigg\{-\frac{(\ln x-d_t)^2}{2\sigma^2 t}\bigg\}$$ which is not Gaussian and not a shifted Gaussian distribution. The mean and second moment are, by using the Gaussian MGF, $$E[X_t]=X_0e^{(\mu-\frac{1}{2}\sigma^2)t}E[e^{\sigma W_t}]=X_0e^{(\mu-\frac{1}{2}\sigma^2)t}e^{\frac{\sigma^2t}{2}}=X_0e^{\mu t}$$ $$E[X_t^2]=X_0^2e^{2(\mu-\frac{1}{2}\sigma^2)t}E[e^{2\sigma W_t}]=X_0^2e^{2(\mu-\frac{1}{2}\sigma^2)t}e^{2\sigma^2 t}=X_0^2e^{2\mu t+\sigma^2 t}$$ The variance is $$\textrm{Var}[X_t]=E[X_t^2]-E[X_t]^2=X_0^2e^{2\mu t}(e^{\sigma^2 t}-1)$$ and one can use the usual estimators for fixed $t$. Since $X_t$ is a Ito process, the quadratic variation is $$d[X,X]_t=\sigma^2X_t^2dt\implies [X,X]_t=\sigma^2\int_{[0,t]}X_s^2ds$$ which is not the QV of a Brownian motion, which is $t$. This value can be approximated pathwise in a simulation.