% Reproduces the implied Trinomial tree from Clewlow & Strickland's book
% "Implementing Derivatives Models, Chapter 5.
% Reproduces Figures 5.3, 5.6, and 5.9.

clc; clear;

% Inline function for the Black Scholes call.
C = inline('s.*normcdf((log(s./K) + (r+v.^2./2).*T)./v./sqrt(T)) - K.*exp(-r.*T).*normcdf((log(s./K) + (r+v.^2./2).*T)./v./sqrt(T) - v.*sqrt(T))',...
	's','K','r','v','T');

% Inline function for the Black Scholes put.
P = inline('K.*exp(-r.*T).*normcdf(-(log(s./K) + (r+v.^2./2).*T)./v./sqrt(T) + v.*sqrt(T)) - s.*normcdf(-(log(s./K) + (r+v.^2./2).*T)./v./sqrt(T))',...
	's','K','r','v','T');

% Input option features.
T = 1;                  % Maturity.
Spot = 100;             % Spot price.
v = 0.15;               % Volatility skew.
r = 0.06;               % Risk free rate.
N = 4;                  % Number of steps.
maxSig = 0.2914;        % Maximum sigma.
dt = T/N;               % Time step.
dx = maxSig*sqrt(3*dt); % Stock price step.

% Create the strikes.
K = Spot*exp([N:-1:-N]'.*dx);

% Create the time steps.
time = [0:dt:T];

% Create the implied volatility surface, 
% and the calls and puts based on this surface.
% Note the functional form of the implied volatilty used
% by Clewlow and Strickland: IV(K,T) = a + b(K-Spot)^2*exp(-cT).
a = .2;
b = .00001;
c = 2;
for t=1:length(time);
	for k=1:length(K);
		IV(k,t) = a + b*(K(k) - Spot)^2*exp(-c*time(t));
		if t>1
			Call(k,t) = C(Spot,K(k),r,IV(k,t),time(t));
			Put(k,t)  = P(Spot,K(k),r,IV(k,t),time(t));
		else
			Call(k,t) = max(Spot - K(k), 0);
			Put(k,t)  = max(K(k) - Spot, 0);
		end
	end
end

% Override with calls and puts with Clewlow and Strickland values.
% This ensures that all results are exactly in-line with theirs.
% These overrides are only used to reproduce their example,
% and must be removed for general use.
  CallOverride = [ ...
   0.0000000000000000	 0.0000001592388074	 0.0000865556968692	 0.0011125848992856	 0.0055197870895514;
   0.0000000000000000	 0.0000225006924577	 0.0049681379449013	 0.0412422862265776	 0.1398623494408960;
   0.0000000000000000	 0.0403082950717086	 0.4194965534523710	 1.1158727257204500	 1.9965170268512400;
   0.0000000000000142	 4.7468781888800400	 7.1558872433358600	 9.1754456282794400	10.9895467387105000;
  22.3034866353223000	23.4718876538093000	24.7109893309608000	25.9912065293228000	27.2682274946354000;
  39.6325181097246000	40.5312733542971000	41.4169658955222000	42.2921635514483000	43.1598301216955000;
  53.0965713652027000	53.7948724435416000	54.4827772803334000	55.1604490030943000	55.8281067499040000];
  Call(2:8,:) = CallOverride;

  PutOverride = [...
  113.2040298772760000  110.0298355628310000 106.9029850846420000  103.8236282642440000 100.7935135427450000;
   65.6520975676292000	 63.1858816348555000  60.7613063453924000	58.4042304212511000	 56.1451327524717000;
   28.7059041255020000	 26.8300311465108000  25.3215663533683000	24.1583929667528000	 23.2071727951437000; 
    0.0000000000000000	  3.2580721491862900   4.2004405981866700	 4.7751938115894200	  5.1660000971353900;
    0.0000000000000000	  0.0116506348822343   0.1112236980042950	 0.2688776531664620	  0.4400481646434410;
    0.0000000000000000	  0.0000005281791151   0.0003190675104799	 0.0033242231568523	  0.0117835477462800;
   0.0000000000000000	  0.0000000000005577   0.0000001070836316	 0.0000086672989891   0.0000923416559747];
  Put(2:8,:) = PutOverride;
  clear CallOveride PutOverride
% End of overrides.

% Set the prices and initialize the Arrow-Debreu state prices, Q.
S = K;
Q = zeros(2*N+1,N+1);

% Middle node, time zero.
Q(N+1,1) = Call(N+2,1)/(S(N+1) - S(N+2));

% Figure 5.3 of Clewlow & Strickland.
% Construct the tree for the Arrow-Debreu prices, Q.
% Upper Nodes of the tree for Q(i,j).
for j=2:N+1
	for i=N-j+2:N+1;
		SumUp = 0;
		for k=1:i-1
			SumUp = SumUp + Q(k,j)*(S(k)-S(i+1));
		end
		Q(i,j) = (Call(i+1,j) - SumUp) / (S(i) - S(i+1));
	end
end

% Lower Nodes of the tree for Q(i,j).
for j=2:N+1
	for i=N+j:-1:N+2
		SumLow = 0;
		for k=N+j:-1:i+1
			SumLow = SumLow + Q(k,j)*(S(i-1) - S(k));
		end
		Q(i,j) = (Put(i-1,j) - SumLow) / (S(i-1) - S(i));
	end
end

fprintf('The Arrow Debreu prices -------------------------------')
Q

% Figure 5.6 of Clewlow & Strickland.
% Construct the trees for the lower (pd), middle (pm),
% and upper (pu) transition probabilities.
infl = exp(r*dt);          % Inflation factor.
pu = zeros(2*N+1,N+1);     % Up probabilities.
pm = zeros(2*N+1,N+1);     % Middle probabilities.
pd = zeros(2*N+1,N+1);     % Down probabilities.
st = zeros(2*N+1,N+1);     % Stability check.

% Upper and middle nodes.
for j=N:-1:1;
	for i=N-j+2:N+1
		if i==N-j+2
			pu(i,j) =  infl*Q(i-1,j+1)/Q(i,j); 
			pm(i,j) = (infl*S(i) - S(i+1) - pu(i,j)*(S(i-1)-S(i+1))) / (S(i) - S(i+1));
			pd(i,j) = 1 - pu(i,j) - pm(i,j);			
		elseif i==N-j+3
			pu(i,j) = (infl*Q(i-1,j+1) - pm(i-1,j)*Q(i-1,j)) / Q(i,j);
			pm(i,j) = (infl*S(i) - S(i+1) - pu(i,j)*(S(i-1)-S(i+1))) / (S(i) - S(i+1));
			pd(i,j) = 1 - pu(i,j) - pm(i,j);
		else
 			pu(i,j) = (infl*Q(i-1,j+1) - pd(i-2,j)*Q(i-2,j) - pm(i-1,j)*Q(i-1,j)) / Q(i,j);
 			pm(i,j) = (infl*S(i) - S(i+1) - pu(i,j)*(S(i-1)-S(i+1))) / (S(i) - S(i+1));
			pd(i,j) = 1 - pu(i,j) - pm(i,j);

		end
		st(i,j) = sqrt(3*(pu(i,j) + pd(i,j) - pu(i,j)^2 + 2*pu(i,j)*pd(i,j) - pd(i,j)^2));
	end
end

% Lower nodes.
for j=N:-1:2
	for i=N+j:-1:N+2
		if i==N+j
			pd(i,j) = infl*Q(i+1,j+1)/Q(i,j);
			pm(i,j) = (infl*S(i) - S(i-1) - pd(i,j)*(S(i+1) - S(i-1))) / (S(i) - S(i-1));
			pu(i,j) = 1 - pd(i,j) - pm(i,j);
		elseif i==N+j-1
			pd(i,j) = (infl*Q(i+1,j+1) - pm(i+1,j)*Q(i+1,j)) / Q(i,j);
			pm(i,j) = (infl*S(i) - S(i-1) - pd(i,j)*(S(i+1) - S(i-1))) / (S(i) - S(i-1));
			pu(i,j) = 1 - pd(i,j) - pm(i,j);
		else
			pd(i,j) = (infl*Q(i+1,j+1) - pu(i+2,j)*Q(i+2,j) - pm(i+1,j)*Q(i+1,j)) / Q(i,j);
			pm(i,j) = (infl*S(i) - S(i-1) - pd(i,j)*(S(i+1) - S(i-1))) / (S(i) - S(i-1));
			pu(i,j) = 1 - pd(i,j) - pm(i,j);
		end
		st(i,j) = sqrt(3*(pu(i,j) + pd(i,j) - pu(i,j)^2 + 2*pu(i,j)*pd(i,j) - pd(i,j)^2));
	end
end

fprintf('The "up" probabilities -------------------------------')
pu
fprintf('The "middle" probabilities ---------------------------')
pm
fprintf('The "down" probabilities -----------------------------')
pd

% Figure 5.9 of Clewlow & Strickland.
% Value the American Up and Out Put barrier option via the implied tree.
% Use the pu, pm, and pd probabilities from the implied tree.
Rebate = 1;             % Rebate when the put gets knocked out
Strike = 100;           % Strike price
H = 110;                % Barrier
V = zeros(2*N+1,N+1);   % Value of the option

% At maturity, the UO put is worth max(0,K-S(T)) if barrier is not breached, 
% and it is worth the rebate otherwise.
for i=1:2*N+1;
	if S(i) < H
		V(i,N+1) = max(Strike - S(i));
	else
		V(i,N+1) = Rebate;
	end
end

% Remaining nodes of the tree
for j=N:-1:1
	for i=N-j+2:N+j
		if S(i) < H
			V(i,j) = max(exp(-r*dt)*(pu(i,j)*V(i-1,j+1) + pm(i,j)*V(i,j+1) + pd(i,j)*V(i+1,j+1)), Strike - S(i));
		else
			V(i,j) = Rebate;
		end
	end
end

fprintf('The UO American put ----------------------------------')

UOPut = V(N+1,1)