Euler Explict Method for solving the PDE for option prices under the Schwartz mean reverting model. Numerical Finance

I have to solve the following PDE for a Call option :

$$\partial_tV + \{ \alpha - (\mu - \lambda/ \alpha -log(S))\}S\partial_SV + 1/2 \sigma^{2}S^{2}\partial_{S}^{2}V - rV = 0$$

Where V(S,t) is the price of the option as a function of stock price S and time t, $$\alpha$$ the strength of mean reversion, $$\mu$$ mean rate of return, $$r$$ is the risk-free interest rate, and $$\sigma$$ is the volatility of the stock, $$\lambda$$ Market price of Risk. I have solved the PDE with finite difference for S and, as required for the exercise, Euler Explict for the time (I upload the Matlab code here with all the data). My question is about the boundary condictions. For the Black-Scholes equation is easy:

• Left end S = 0: So payoff(S) = payoff(0) = max(0-K,0) = 0, hence $$V(0,t) = 0 \forall t\ge 0$$.
• Right end S = M: Where M is 'large' with respect to K, payoff(S) = max(S-K,0) = S - K, so the value of an option with payoff equal to S - K is just $$V(S,t) = S - Ke^{-rT}$$, $$\forall t\ge 0$$.

The left end condiction remains the same, but i don't understand how to change the right end.

% Generate a Path for the Schwartz single factor model.
% Call
% We start from : dS = ...., then X = log(S).
clear all
close all

%% Parameters
T = 1;           % End time
K = 1.5;         % Strike price
alpha = 10;     % strength of mean reversion
mu = 1;         % long term mean
sigma = 0.7;    % volatility
r = 0.1;        % risk-free interest rate
lambda = mu - r;% Market Price of Risk
lambdaHat = lambda/alpha;
Nt = 5000;      % number of time steps

%% Eu for PDE
ds = 0.01;
dt = T/Nt;
t = 0:dt:T;
Ns = round(3*(K/ds));   % Number of "price" step
S = ds*(0:(Ns));
S1 = S(2:(end-1));       % all but endpoints
V = max(S-K,0);
hold on
plot(S,V,'o-b')

%% PDE solution
Ws = alpha*(mu - lambdaHat - log(S1)).*S1;
Zs = 0.5*sigma^2*S1.^2;
Vnew = zeros(Ns+1,1);
Vold = V;
for k = 2:Nt
F = zeros(Ns+1,1)';
dV = (Vold(3:end) - Vold(1:(end-2)))/(2*ds);
ddV = (Vold(3:end) - 2*Vold(2:(end-1)) + Vold(1:(end-2)))/(ds^2);
F(2:(end-1)) = Ws.*dV + Zs.*ddV - r*Vold(2:(end-1));
Vnew = Vold + dt*F;
%% Boundary cond
Vnew(1) = 0; % left end condition
Vnew(end) = S(end) - K*exp(-r*k*dt); % right end condition
% Update solution
Vold = Vnew;
end
hold on
plot(S,Vnew,'-r')

• Did you find out the right end condition? – sound wave Sep 23 '20 at 18:29