Simulation of Geometric Brownian Motion

Question

I generate 10000 random binomial paths for a stock whose price is from S(0) = 10 out to S(t) where t = 1 year. Assume geometric Brownian motion for the stock price with a drift of 15% per year and a volatility of 20%. I use 10000 equally spaced time steps of length along each path and $\Delta S = \mu S \Delta t \pm \sigma S \sqrt{\Delta t}$, where the + or – moves at each time step are generated at random with equal probability.

S = 10
S1 = []
mean = 0.15
sd = 0.2

for i in range(10000):
    S = 10

    for j in range(10000):
        roll = np.random.rand()
    
        if roll < 0.5:
            S = S + mean*S/10000 + sd*S/100
        else:
            S = S + mean*S/10000 - sd*S/100

    S1.append(S)

I use the above codes to simulate the GBM and my result is the same as the theoretical values. Mean should be $S(0)e^{\mu t}$ and the variance should be $S^2(0)e^{2\mu t}(e^{\sigma^2 t} - 1)$. But there is one question want me to calculate the inverse, which is $G(t) = 1/S(t)$. I consider that mean should be $\frac{1}{S(0)}e^{\mu t}$ and variance should be $\frac{1}{S^2(0)}e^{2\mu t}(e^{\sigma^2 t} - 1)$. Am I correct? But I'm not sure how to change my codes. I have tried to write down the following codes. But the result is far from the theoretical values mean = 0.116 and variance = 0.00055. Could someone explain where I'm doing wrong? Thank you very much.

S = 10
G1 = [1/10]
mean = 0.15
sd = 0.2

for i in range(10000):
    S = 10

    for j in range(10000):
        roll = np.random.rand()
    
        if roll < 0.5:
            G = 1/S
            S = S + mean*S/10000 + sd*S/100
            G = G + mean*(1/S)/10000 + sd*(1/S)/100
        else:
            G = 1/S
            S = S + mean*S/10000 - sd*S/100
            G = G + mean*(1/S)/10000 - sd*(1/S)/100

    G1.append(G)

Answer 1

Just to explain the formulae for the mean and variance! We can start with the solution of the GBM SDE:

$S_t=S_0\,e^{\left(\mu-\frac{1}{2}\sigma^2\right)t+\sigma W_t}$

then,

$\frac{1}{S_t}=\frac{1}{S_0}\,e^{-\left(\mu-\frac{1}{2}\sigma^2\right)t-\sigma W_t}$

The mean (and variance) are easily obtained via this identity for a Gaussian random variable Y:

$E\left[e^Y\right]=e^{E\left[Y\right]+\frac{1}{2}V\left[Y\right]}$

Applying this to S, and noting that:

$E\left[\left(\mu-\frac{1}{2}\sigma^2\right)t+\sigma W_t\right]=\left(\mu-\frac{1}{2}\sigma^2\right)t$

$V\left[\left(\mu-\frac{1}{2}\sigma^2\right)t+\sigma W_t\right]=\sigma^2t$

we get the familiar expression for the mean of S:

$E\left[S_t\right]=S_0\,e^{\left(\mu-\frac{1}{2}\sigma^2\right)t+\frac{1}{2}\sigma^2 t}=S_0\,e^{\mu t}$

And applying the identity to the expression for 1/S, we get its mean as follows:

$E\left[\frac{1}{S_t}\right]=\frac{1}{S_0}\,e^{-\left(\mu-\frac{1}{2}\sigma^2\right)t+\frac{1}{2}\sigma^2 t}=\frac{1}{S_0}\,e^{-\mu t+\sigma^2t}$

Variance can be determined using the same identity.