Monte Carlo Pricing of a European Call with Exact GBM Sampling

# Pricing European call with Feynman-Kac


# Pricing European call with Feynman-Kac












I am trying to calculate the solution to the Black-Scholes (BS) equation using the Feynman-Kac (FK) formula for a simple European call. According to FK, the solution to BS is the discounted average of the process `S(T)` where $S$ follows a geometric Brownian motion.

The analytical solution to the BS-equation I have implemented in Python:

```
from math import *

#first define these 2 functions
def d1(S,X,T,r,sigma):
    return (log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T))

def d2(S,X,T,r,sigma):
    return d1(S,X,T,r,sigma)-sigma*sqrt(T)

#define the call option price function
def bs_call(S,X,T,r,sigma):
     return S*CND(d1(S,X,T,r,sigma))-X*exp(-r*T)*CND(d2(S,X,T,r,sigma))

#define cumulative standard normal distribution
def CND(X):
     (a1,a2,a3,a4,a5)=(0.31938153,-0.356563782,1.781477937,-1.821255978,1.330274429)
     L = abs(X)
     K=1.0/(1.0+0.2316419*L)
     w=1.0-1.0/sqrt(2*pi)*exp(-L*L/2.)*(a1*K+a2*K*K+a3*pow(K,3)+a4*pow(K,4)+a5*pow(K,5))
     if X<0:
        w=1.0-w
     return w
```

I evaluate it using these parameters (risk-free return $r$, volatility $\sigma$, maturity $K$, time-to-maturity `ttm` and underlying stock price $S$:

```
r = 0.05
sigma = 0.003
K = 25
S = 30
ttm = 3
bs_call(S, K, ttm, r, sigma)
```

which gives `8.48`. I now try with FK:

```
def X_T(r, S0, sigma, T, n):
    dt = T / n
    X = S0
    for i in range(n):
        dx = r*X*dt + sigma*X*np.random.normal(0, np.sqrt(dt))
        X = X + dx
    return X

ave = 0
N = 1000
for i in range(N):
    ave += np.max([X_T(r, S, sigma, 10, 1000) - K, 0])
print(ave*np.exp(-r*3) / N)
```

which gives `21.05`, which doesn't give the BS-result (I have taken $T=10$).

Where does the discrepancy arise?

## Answer by siou0107 (score 3, accepted)

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

First, you'd rather simulate $\log(X)$ rather than $X$; thus, there is no level dependency in your discretisation scheme, making it more accurate. $$Z_t = \log(S_t)$$ $$dZ_t = \left(r - \frac{\sigma^2}{2}\right)dt + \sigma dW_t$$ You can even run one single time step, and the distribution of your final price will still be as accurate!

Second, the different time to maturity (3 vs 10) certainly explains a large part of the difference :)

Third, I would advise you to take advantage of vectorisation rather than using a loop.

My suggestions are summarised in the following code:

```
Z_T = np.log(S0) + (r - 0.5 * sigma**2)*ttm + sigma*np.random.normal(0, np.sqrt(ttm), N)
S_T = np.exp(Z_T)
payoff = np.maximum(S_T - K, 0)
price = np.exp(-r*ttm)*np.average(payoff)
price
```

That yields 8.48 as well.