% American put valuation by Crank-Nicolson Method.
% Example taken from Clewlow and Strickland's book.
% "Implementing Derivatives Models", Figure 3.17 (page 76).

clc; clear;

% Option features
T = 1;            % Maturity
Spot = 100;       % Spot price
K = 100;          % Strike price
v = 0.20;         % Volatility
r = 0.06;         % Risk free rate
q = 0.03 ;        % Dividend yield

% PDE features.
N = 3;                 % Number of time steps
M = 3;                 % Number of stock price steps
dt = T/N;              % Time increment
mu = r - q - v^2/2;    % Drift for stock process
dx = v*sqrt(3*dt);     % Increment for stock price

% Probabilities.
pu = -1/4*dt*(v^2/dx^2 + mu/dx);        % Up probability
pm =  1 + dt*v^2/2/dx^2 + r*dt/2;       % Middle probability
pd = -1/4*dt*(v^2/dx^2 - mu/dx);        % Down probability

% Initialize stock price and option value.
S = zeros(2*M+1, N+1);
V = zeros(2*M+1, N+1);

% Indices for stock price step.
J = [M:-1:-M]';

% Stock price at maturity
S = Spot*exp(J.*dx);

% Option price at maturity
V(:,end) = max(K - S, 0);

% Initialize the required matrices
pmp = zeros(2*M+1,N+1);
pp  = zeros(2*M+1,N+1);
C   = zeros(2*M+1,N+1);

% Create the P' matrix
for j=N+1:-1:2
	pmp(2*M,j) = pd+pm;
 	for i=2*M-1:-1:2
 		pmp(i,j) = pm - pd/pmp(i+1,j)*pu;
 	end
end

% Upper and lower bounds for the American put
lambda_L = S(2*M+1) - S(2*M);
lambda_U = 0;

% Work backwards and obtain the matrix of values, V
for j=N+1:-1:2;
	for i=2*M:-1:2
		if i==2*M
			pp(i,j) = -pu*V(i-1,j) - (pm-2)*V(i,j) - pd*V(i+1,j) + pd*lambda_L;
		else
			pp(i,j) = -pu*V(i-1,j) - (pm-2)*V(i,j) - pd*V(i+1,j) - pd/pmp(i+1,j)*pp(i+1,j);
		end
	end
	j=j-1;
	for i=1:2*M+1
		if i==1
			C(i,j) = (pp(i+1,j+1) + pmp(i+1,j+1)*lambda_U) / (pmp(i+1,j+1) + pu);
			V(i,j) = max(K - S(i), C(i,j));
		elseif i<2*M+1
			C(i,j) = (pp(i,j+1) - pu*C(i-1,j))/pmp(i,j+1);
			V(i,j) = max(K - S(i), C(i,j));
		else
			C(i,j) = C(i-1,j) - lambda_L;
			V(i,j) = max(K - S(i), C(i,j));
		end
	end
	j=j+1;
end

% Output the value matrix and the price of the American put
V
AmerPutPrice = V(M+1,1)