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

Question

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?

Thanks in advance.

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);

Answer 1

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.