Verification of order of convergence of Implicit Euler Method to numerically solve Black-Scholes PDE

I'm trying to verify the order of convergence for implicit Euler method to numerically solve Black-Scholes PDE. Theory says that it should be $$O(\Delta t + \Delta S^2).$$ My code is working absolutely fine. I'm getting almost zero error when I compare my solution to the one obtained by Black-Scholes formula. However, I'm getting incorrect order of convergence. Does it have something to do with how I choose parameters? Could someone have a look at my code and tell me where I'm going wrong if there is a mistake somewhere? Or perhaps convey to me what option parameters need to be chosen so that I get the correct order of convergence?

Here is the MATLAB code:

function [V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS)
% [V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS)
% Prices European put or call option using implicit finite difference method.
% -- Input arguments --
%  pc = 'put' or 'call'
%  K  = strike price
%  T  = expiry time
%  r  = risk-free interest rate
%  sigma = volatility
%  Smax = max value of S used in finite difference grid
%  Nt = number of time steps
%  NS = number of grid points along S axis
% -- Output arguments --
%  V  = values of option at asset values in S
%  t  = vector of time points
%  S  = vector of asset prices

% Time discretization
Dt = T / Nt;
t = linspace(0, T, Nt+1);

% Asset price discretization
DS = Smax / (NS);
S = linspace(0, Smax, NS+1)';

% Storage for option values
V = zeros(length(S), length(t));

% Set initial condition at expiry and boundary values
switch lower(pc(1))
case 'p' % put option
V(1,:) = K * exp(-r*(T-t));
V(end,:) = 0;
V(1:NS-1,end) = max( K - S(2:end-1), 0);
case 'c' % call option
V(1,:) = 0;
V(end,:) = Smax - K * exp(-r*(T-t));
V(2:end-1,end) = max( S(2:end-1) - K, 0);
otherwise
error(['unknown option type ' pc])
end

% Define index vector of interior asset prices
J = 2:NS;

c1 = sigma^2 * S(J).^2 * Dt / ( 2 * DS^2 );
c2 = r * S(J) * Dt / ( 2 * DS );

alpha = -c1 + c2;
beta  = 1 + r*Dt + 2*c1;
gamma = -c1 - c2;

% Create and factor the sparse matrix for the linear system.
A = spdiags([[alpha(2:NS-1); 0], beta, [0;gamma(1:NS-2)]], [-1, 0, 1], NS-1, NS-1);
[L,U,p] = lu(A,'vector');

% Time step backwards from expiry
for n = Nt:-1:1
% Set up RHS vector from values at time step n+1
b = V(J,n+1);
% Adjust for boundary values at S = 0 and S = Smax
b(1) = b(1) - alpha(1)*V(1,n);
b(NS-1) = b(NS-1) - gamma(NS-1)*V(NS+1,n);

% Solve linear system A * V(J,n) = b using existing LU factorization
% V(J,n) = A \ b; would recalculate factorization at each time step
V(J,n) = U \ ( L \ b(p) );
end
function [c, dcds] = blackscholes(S, K, r, sigma, Tmt)
% [c, dcds] = blackscholes(S, K, r, sigma, Tmt)
% Black and Scholes formula for the value of a call option
% and its derivative with respect to volatility sigma
% S   = underlying asset price
% K   = strike price
% r   = risk-free interest rate
% Tmt = time to maturity = T - t where T = expiry
% If sigma is a vector of volatilities, then both the
% call value and its derivatives are vectors of the same size.
%
% Uses normpdf and normcdf from Statistics toolbox.

s = sigma * sqrt(Tmt);

d1 = ( log(S/K) + ( r + sigma.^2/2)*(Tmt) ) ./ s;
d2 = d1 -  s;

% Use normpdf and normcdf from Statistics toolbox
c = S .* normcdf(d1) - K * exp(-r*Tmt) * normcdf(d2);

% Derivative of call value w.r.t. volatility sigma
dcds = S .* normpdf(d1) * sqrt(Tmt);
% Test script for the functions to solve Black - Scholes PDE
% for the value of a call or put option.
% Define option parameters
r = 0.1;
sigma = 0.4;
T = 5/12;
K = 50;
pc = 'call';

% Define discretization parameters
% Asset price S goes from 0 to Smax
Smax = 100;

% Compare with Black and Scholes formula for call
Tmt = T;
[c, dcds] = blackscholes(S, K, r, sigma, Tmt);
Err = c - V(:,1);

fprintf("\nConvergence with respect to Dt\n")

fprintf("\n%6s %10s %10s %8s\n\n","NS", "Nt", "max error", "rate")

nrows = 8;
NS = 100;
Nt = 200;
err = zeros(nrows,1);
for row = 1:nrows
Nt = 2 * Nt;
[V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS);
for n = 1:Nt+1
err_n = norm( V(:,n) - blackscholes(S, K, r, sigma, T - t(n)), Inf);
if err_n > err(row)
err(row) = err_n;
end
end
if row == 1
fprintf("%6d %10d %10.2e\n", NS, Nt, err(row))
else
ratio = err(row-1) / err(row);
rate = log2(ratio);
fprintf("%6d %10d %10.2e %8.3f\n", NS, Nt, err(row), rate)
end
end

fprintf("\nConvergence with respect to Ds\n")

fprintf("\n%6s %10s %10s %8s\n\n","Ns", "Nt", "max error", "rate")

nrows = 5;
NS = 50;
Nt = 2000;
err = zeros(nrows,1);
for row = 1:nrows
NS = 2 * NS;
%[U, x, t] = iEuler(a, L, T, f, gamma0, gammaL, u0, NS, Nt);
[V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS);
for n = 1:Nt+1
err_n = norm( V(:,n) - blackscholes(S, K, r, sigma, T - t(n)), Inf);
if err_n > err(row)
err(row) = err_n;
end
end
if row == 1
fprintf("%6d %10d %10.2e\n", NS, Nt, err(row))
else
ratio = err(row-1) / err(row);
rate = log2(ratio);
fprintf("%6d %10d %10.2e %8.3f\n", NS, Nt, err(row), rate)
end
end
$$`$$

You forgot to multiply the norm of the difference between the numerical solution and the exact one by the discretization step. In your case it is enough to divide the err_n by Nt when computing the order of time discretization, and by NS when computing the order of space discretization. I got a perfect first order accuracy for the time discretization, in the case of space discretization the order seems to approach the first one from above as the initial function is not differentiable in the space.

Here is the output of the corrected code:

Convergence with respect to Dt

NS         Nt  max error     rate

100        400   2.59e-04
100        800   1.19e-04    1.120
100       1600   5.79e-05    1.042
100       3200   2.83e-05    1.032
100       6400   1.40e-05    1.017
100      12800   6.95e-06    1.009
100      25600   3.46e-06    1.005
100      51200   1.73e-06    1.002

Convergence with respect to Ds

Ns         Nt  max error     rate

100       2000   9.17e-04
200       2000   2.53e-04    1.858
400       2000   7.03e-05    1.848
800       2000   2.20e-05    1.678
1600       2000   8.92e-06    1.300
• Thanks!! It works now but I'm not sure as to why we need to divide err_n by Nt when computing the order of time discretization, and so forth by NS . We know that order of convergence $\rho = \frac{log(E_1/E_2)}{log(h_1/h_2)}$ where $E_1 = u_{h_1} - u*$ and $E_2 = u_{h_2} - u*.$ My code calculates rate = log2(ratio) which is the same as the formula. So I don't understand why we need to do the divisions that you have suggested? Also, not sure why you say multiply in the first line when you are actually dividing and how did you deduce that initial function is not differentiable in the space
– user36184
Jan 24 '21 at 20:01
• I am answering on a phone. Have a look to definitions of errors, your one is wrong. You have to e.g. multiply the difference with e.g. DeltaS that is same as divided by NS. Because DeltaS=Smax/NS, it is the same. The initial condition is not differentiable, because it contains the function max that is not differentiable in zero. Jan 25 '21 at 6:16
• Thanks!!! Could you please refer me to some relevant literature where the definition of error is given. I'm having a hard time finding out relevant resources. Also, in my code I'm using infinity norm to compute error , what will happen if I use L2 norm and do I have to make any changes?
– user36184
Jan 25 '21 at 9:53
• A simple explanation for a "correct way" to compute the error is that you should do a fair comparison for different discretizations steps (i.e. different NS and Nt). Then an "averaged error", e.g. the sum of absolute differences in each grid node divided by the number of nodes, can do the job. More theoretically based approach is that you should place (extend) all numerical solutions to one functional space, e.g. using piecewise linear functions by interpolating values in nodes, and compare errors with a norm of that space (e.g. an integral norm) using enough precise numerical quadrature. Jan 26 '21 at 15:29
• A good reference to be found also on the web is the book of Randy LeVeque on "Finite Difference Methods for Ordinary and Partial Differential Equations", where there is a Chapter (Appendix) on Measuring errors. Jan 26 '21 at 15:39