0
$\begingroup$

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
```
Improve this question
$\endgroup$
1
$\begingroup$

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
Improve this answer
$\endgroup$
5
  • $\begingroup$ 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 $\endgroup$ – Rhombus Jan 24 at 20:01
  • $\begingroup$ 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. $\endgroup$ – Peter Frolkovič Jan 25 at 6:16
  • $\begingroup$ 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? $\endgroup$ – Rhombus Jan 25 at 9:53
  • 1
    $\begingroup$ 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. $\endgroup$ – Peter Frolkovič Jan 26 at 15:29
  • 1
    $\begingroup$ 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. $\endgroup$ – Peter Frolkovič Jan 26 at 15:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.