Pricing European call with Feynman-Kac

Question

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 1

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.