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?

• Should your T and ttm both be 3? Why are they different? – Gordon Nov 26 '19 at 21:29
• @Gordon $T$ is the maturity time, but $T-t$ is the time-to-maturity (right?). As I understand it, I have to discount by the time-to-maturity $T-t=10-7=3$, but integrate the geometric Brownian motion up to $T=10$ (?) I get this from Feynman-Kac: $C(S, t=7) = E[e^{-r(T-t)}C(T=10, S)]$ – Tyler D Nov 26 '19 at 21:43
• I do not completely understand. But try setting both to 3 and see what will happen – Gordon Nov 26 '19 at 21:49
• It is a bad coding practice to have exp(-r*3). What is 3, is it the time to maturity? If so it should be ttm or some other symbolic reference. Same for X_T(r, S, sigma, 10, 1000), what is 10? It seems inconsistent with 3. – Alex C Nov 27 '19 at 2:47

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.