# Monte Carlo and PDE results are different for a Call Option!

Okay so this might be a fairly trivial question but I'm having an issue with valuing a call option using both a Monte Carlo method and a PDE method.

When I started I first used the parameters:

Spot = 0 to 20

Strike = 10

Interest rate = 0

Volatility = 0.25

Time = 1

And both the Monte Carlo and PDE methods came out identical. However, when I changed the spot price range from 90 to 110 and the strike to 100, the Monte Carlo and PDE methods now give different results! For instance, with a spot price of 110, the Monte Carlo method gives (approximately) 16.19 as the option price. However, the PDE method instead gives 10 as the price. And the thing is, this happens with every Monte Carlo and PDE code I try (I found a few off the internet and I've been trying them out). Can someone tell me why this happens?

Edit: Just as a small favor, can someone tell me what needs to be modified in this code in order to give the true result? I got the code from the link I mentioned below

r=0.0; % Interest rate
sigma=0.25; % Volatility of the underlying
M=1600; % Number of time points
N=160; % Number of share price points
Smax=110; % Maximum share price considered
Smin=90; % Minimum share price considered
T=1.; % Maturation (expiry)of contract
E=100; % Exercise price of the underlying

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=(T/M); % Time step
ds=(Smax-Smin)/N; % Price step
% Initializing the matrix of the option value
v(1:N,1:M) = 0.0;
% Initial conditions prescribed by the European Call payoff at         expiry:V(S,T)=max(S-E,0);
v(1:N,1)=max((Smin+(0:N-1)*ds-E),zeros(size(1:N)))';
%v(1:N,1)=max((Smin+(0:N-1)*ds-E),(5))';

% Boundary conditions prescribed by the European Call:
v(1,2:M)=zeros(M-1,1)'; % V(0,t)=0
v(N,2:M)=((N-1)*ds+Smin)-E*exp(-r*(1:M-1)*dt); % V(S,t)=S-Eexp[-r(T-t)] as S ->infininty.

% Determining the matrix coeficients of the explicit algorithm
aa=0.5*dt*(sigma*sigma*(1:N-2).*(1:N-2)-r*(1:N-2))';
bb=1-dt*(sigma*sigma*(1:N-2).*(1:N-2)+r)';
cc=0.5*dt*(sigma*sigma*(1:N-2).*(1:N-2)+r*(1:N-2))';

% Implementing the explicit algorithm
for i=2:M,
v(2:N-1,i)=bb.*v(2:N-1,i-1)+cc.*v(3:N,i-1)+aa.*v(1:N-2,i-1);
end
% Reversal of the time components in the matrix as the solution of the     BlackScholes
% equation was performed backwards
v=fliplr(v);

• You must be doing something wrong. Are you sure you use the same kind of data in input e.g. annualised rates and volatilities, rates expressed in % and not in absolute values etc. Also, is it the MC which fails or the PDE, I mean did you compare to analytic BS price? Jun 21, 2016 at 10:52
• @ThePlowKing Do you have Matlab ?
– user16651
Jun 21, 2016 at 10:58
• If $S_0=110$ ,$K=100$ , $\tau=1$ ,$r=0$ and $\sigma=0.25$ and $q=0$ then call option price=$16.19$ . Check blsprice$(110,100,0,1,0.25,0)$
– user16651
Jun 21, 2016 at 11:01
• stability of your numerical method is not enough.
– user16651
Jun 21, 2016 at 11:03
• @Behrouz Maleki sorry about that, yes I actually got about 16.19 for the Monte Carlo method, I accidentally used my old strike price of 98 which is why I put 17.4174 so I'll edit my post now. And yes I do have MATLAB, I actually downloaded two PDE codes off that (in addition to the one I wrote in C++) and they all seem to give me a value of 12 for the option price for some reason. For instance, I used the MATLAB PDE code from this link: ehu.eus/aitor/irakas/fin/apuntes/pde.pdf (as on page 47) and when I plugged in the values above I instead got 12 as the option price Jun 21, 2016 at 11:41

I would definitely recommend Volopta as a reliable source of self-contained and commented financial engineering source codes (useful for prototyping/understanding but clearly not production code). I have for instance copy-pasted, the explicit PDE solver you are looking for (centred in space, backward in time) below (+ edited for clarity + improved performance, see [Edits]).

Running the code in your configuration: $S_0=110, K=100, T=1, r=q=0, \sigma=0.25$ gives $16.1948$ (you can still refine the spot/time meshes, beware of the fact that explicit schemes are only conditionally stable though), while the Black-Scholes price is $16.1904$

% European Call Option by Explicit Finite Differences
% Example taken from Clewlow and Strickland's book
% "Implementing Derivatives Models", Figure 3.8 (page 61)

clc; clear;

% Config
T = 1;                                  % Maturity
Spot = 110;                             % Spot price
K = 100;                                % Strike price
v = 0.25;                               % Volatility
r = 0.0;                                % Risk free rate
q = 0.0;                                % Dividend yield
N = 250;                                % Number of time steps
M = 250;                                % Number of stock price steps

% Time and spot meshes
J = 0:1:N;                              % Indices for time step
dt = T/(length(J)-1);                   % Time increment
I = M:-1:-M;                            % Indices for stock price step
dx = v*sqrt(3*dt);                      % Increment for stock price
S = Spot*exp(I.*dx);                    % Initialize spot mesh (non-unif)

% PDE discretisation
mu = r - q - v^2/2;                     % Drift for stock process
pu = dt*(v^2/2/dx^2 + mu/2/dx);         % Up probability
pm = 1 - dt*v^2/dx^2 - r*dt;            % Middle probability
pd = dt*(v^2/2/dx^2 - mu/2/dx);         % Down probability
V = zeros(length(I), length(J));        % Initialize option price

% Terminal condition (absence of arbitrage)
V(:,end) = max(S - K, 0);

% Work backwards through the lattice
for j=N:-1:1

% Inner spatial nodes at time t_j = (j-1)*dt
I = 2:2*M;
V(I,j) = pu*V(I-1,j+1) + pm*V(I,j+1) + pd*V(I+1,j+1);

% Lower boundary condition @ S_{min} = S(2*M+1)
V(2*M+1,j) = V(2*M,j);

% Upper boundary condition @ S_{max} = S(1)
V(1,j) = V(2,j) + (S(1) - S(2));

end

% Finite difference price
fdPrice = V(M+1,1);

% Analytic BS price
DF0T = exp(-r*T);
F0T = Spot*exp((r-q)*T);
bsPrice = bsPrice(v, DF0T, F0T, K, T, 'C');

fprintf('Finite Difference: %0.4f\n', fdPrice)
fprintf('Analytic: %0.4f\n', bsPrice)


[Edits]

Apologies for I didn't read the code that I provided in my answer: I only made sure it produced the right output. AMOF this code was pretty messy (it mixed the $i$ and $j$ indices used once for space/time respectively then later for time/space), also it had an additional set of lines inside the main backward loop

for i=1:2*M
V(i,j) = max(K - S(j), V(i,j));
end


that made no sense. I've added my own corrected version above. It should run faster also since I've removed useless "for" loops.

As far as your own implementation is concerned, there are so many problems that I would advise you start from the one I gave you and ask questions. If you need proper code review post it on Code Review Stack Exchange. Examples of problems with your implementation include:

• Either embed your MATLAB code inside a proper function, or make sure you clear all previously used variables before you run it. If you don't, your results can be wrong (especially if you change the number of time/space nodes in your case, from big M and N, to small ones).

• You pick $S_{\text{min}}=90$ and $S_{\text{max}}=110$ to price an option struck at $E=100$... $S_{\text{min}}$ and $S_{\text{max}}$ should correspond to the boundaries of the spatial domain i.e., ideally, $S_{\text{min}}\approx 0$ and $S_{\text{max}}\approx +\infty$. By over-restricting the spatial domain, you can't expect to reasonably approach the real solution (especially if your boundary conditions are wrong).

• You do have a problem with your boundary conditions, but also with your matrix coefficients IMHO.

• Thanks heaps for that, I think my issue might actually be a boundary/initial condition problem in the PDE, would it be okay if you could help me out with it? (This is actually the first PDE code I've written so I'm not sure how things work properly yet). For the initial condition I currently have the following: $\max((Smin + j*ds) - K, 0)$ where $(Smin + j*ds)$ is the spot (and in particular $j$ goes from 0 to N and $ds$ is the price increment. And for the boundary condition I have: $(N*ds + Smin) - K*exp(-r*j*dt)$ where $(N*ds + Smin)$ is the max spot price and $dt$ is time increment. Thanks! Jun 21, 2016 at 22:56
• Thanks for all the help, your comment made me realize the issue I had - once I changed $S_{min}$ to 0 and $S_{max}$ to 200 the PDE gave the correct prices. I appreciate all the help! (Also, the MATLAB code I gave isn't the one I used, it's just one I copied off the internet since I'm quite familiar with it and since my code is written in C++) Jun 23, 2016 at 8:05
• Glad I could help :) Jun 23, 2016 at 8:14