计算科学 - 验证隐式欧拉法数值求解 Black-Scholes PDE 的收敛阶 - 吾爱随笔录

我正在尝试验证隐式欧拉方法的收敛顺序以数值求解 Black-Scholes PDE。理论说应该 $O(\Delta t + \Delta S^2).$ 我的代码工作得很好。当我将我的解决方案与 Black-Scholes 公式获得的解决方案进行比较时，我得到的错误几乎为零。但是，我的收敛顺序不正确。它与我如何选择参数有关吗？如果某处有错误，有人可以看看我的代码并告诉我哪里出错了吗？或者也许向我传达需要选择哪些选项参数才能获得正确的收敛顺序？

这是MATLAB代码：

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
```