function xxx = lect03ivtree
clear;
format compact;
format short;
%format long;
pts = load('c:/VolSurf0.dat');
K = pts(:,1);
F = pts(:,2);
C = pts(:,3);
T = pts(:,4)/252;
IV = pts(:,5);
r = 0.03;
VM = mean(IV)
% fit volatility surface
D = mness(K, F, T);
pars = [VM, 5., -1., 0., 0., 0.];
opts = optimset('Display','iter','TolFun',0.00001);
x = lsqnonlin(@ivpolyerr,pars,[], [], opts, D, T, IV);
x
figure(1);
plot3(D,T,IV,'o'); hold on;
[DM,TM] = meshgrid(-0.2:.02:0.6,0:.02:1.2);
ZM = ivpoly(x, DM, TM);
mesh(DM, TM, ZM); hold off;
% 
n = 10;
r = 0.03;
T = 0.2;
[S,p,L] = ivtree(@volfun, n, r, T, x);
Dt = T/n;
SP = [];
for k = 1:n
    SP = [SP ; k*Dt*ones(1+k,1) S(k,1:(1+k))' L(k,1:(1+k))'];
end
figure(3);
stem3(SP(:,1),SP(:,2),SP(:,3))

function y = ivpoly(pars, D0, T0)
D = D0;
T = sqrt(T0);
y = pars(1) + pars(2)*D + pars(3)*T + pars(4)*D.^2 + pars(5)*T.^2 + pars(6)*T.*D;

function y = ivpolyerr(pars, D, T, IV)
y = ivpoly(pars, D, T) - IV;

function y = bsc(K, sig, r, T)
% The simple Black-Scholes European call option Price
d1 = ( -log(K)+(r+.5*sig^2)*T ) / sig / sqrt(T);
d2 = d1 - sig*sqrt(T);
n1 = normcdf(d1);
n2 = normcdf(d2);
y =  n1 - K*exp(-r*T)*n2;

function y = bsp(K, sig, r, T)
% The simple Black-Scholes European put option Price
d1 = ( -log(K)+(r+.5*sig^2)*T ) / sig / sqrt(T);
d2 = d1 - sig*sqrt(T);
n1 = normcdf(-d1);
n2 = normcdf(-d2);
y =  -n1 + K*exp(-r*T)*n2;

function y = mness(K, F, T)
% The moneyness function
y = log(F./K) ./ sqrt(T);

function y = imness(D, F, T)
% The inverse moneyness function
y = F .* exp( - D .* sqrt(T) );

function y = volfun(K, T, pars)
%y = 0.1;
y = ivpoly(pars, mness(K, 1, T), T);
%y = 0.3 + 0.1 * (K-1);

function [S,p,L] = ivtree(SFH, n, r, T, pars)
Dt = T/n;
S = zeros(n+1, n); p = zeros(n+1, n);
L = zeros(n+1, n); F = zeros(n+1);
% n = 1
even = 1;
F(1) = exp(r*Dt); 
sig = feval(SFH, F(1), Dt, pars);
C = bsc(F(1), sig, r, Dt); 
D = exp(r*Dt) * C;
S(1,1) = F(1) * (F(1)-D) / (F(1)+D); 
S(2,1) = F(1) * (F(1)+D) / (F(1)-D);
p(1,1) = (F(1)-S(1,1)) / (S(2,1)-S(1,1));
L(1,1) = exp(r*Dt) * (1-p(1,1)); 
L(2,1) = exp(r*Dt) * p(1,1);
for k=2:n % loop time k
    for j = 1:k % set forwards
        F(j) = S(j,k-1) * exp(r*Dt);
    end
    if (even) % even step
        k0 = k/2; delta = 0; even = 0;
        S(k0+1,k) = exp(r*k*Dt);
    else % odd step
        k0 = (k+1)/2; delta = 1; even = 1;
        sig = feval(SFH, F(k0), k*Dt, pars);
        C = bsc(F(k0), sig, r, k*Dt); % set middle call
        D = exp(r*Dt) * C; % set middle delta
        for i = (k0+1):k
            D = D - L(i,k-1) * (F(i)-F(k0));
        end
        S(k0  ,k) = F(k0) * (L(k0,k-1)*F(k0)-D) / (L(k0,k-1)*F(k0)+D); % set middle prices
        S(k0+1,k) = F(k0) * (L(k0,k-1)*F(k0)+D) / (L(k0,k-1)*F(k0)-D);
    end
    % compute upper part
    for j = (k0-delta):-1:1 % loop states
        % set puts
        sig = feval(SFH, F(j), k*Dt, pars);
        P = bsp(F(j), sig, r, k*Dt);
        % set delta
        D = exp(r*Dt) * P;
        for i = 1:(j-1)
            D = D - L(i,k-1) * (F(j)-F(i));
        end
        % set price
        S(j,k) = ( L(j,k-1) * F(j) * (S(j+1,k)-F(j)) - D * S(j+1,k) ) / ...
                 ( L(j,k-1)        * (S(j+1,k)-F(j)) - D            );
        % check bounds
        if ( j==1 & S(1,k)>F(1) )         
            S(1,k) = 0.95*F(1);
            %S(1,k) = max( F(1) + (S(2,k)-F(1))*(F(2)-S(2,k))/(S(3,k)-F(2)) , 0.9*F(1) );
        end
        if ( j~=1 & ( S(j,k)>F(j) | S(j,k)<F(j-1)  ) )
            %S(j,k) = .5 * (F(j-1)+F(j));
            S(j,k) = 0.75*F(j-1) + 0.25*F(j);
        end
    end
    % compute lower part
    for j = (k0+1):k
        % set calls
        sig = feval(SFH, F(j), k*Dt, pars);
        C = bsc(F(j), sig, r, k*Dt);
        % set delta
        D = exp(r*Dt) * C;
        for i = (j+1):k
            D = D - L(i,k-1) * (F(i)-F(j));
        end
        % set price
        S(j+1,k) = ( D * S(j,k) - L(j,k-1) * F(j) * (F(j)-S(j,k)) ) / ...
                   ( D          - L(j,k-1)        * (F(j)-S(j,k)) );
        % check bounds
        if ( j==k & S(k+1,k)<F(k) )
            S(k+1,k) = 1.05 * F(k);
            %S(k+1,k) = F(k) + (F(k)-S(k,k))*(S(k,k)-F(k-1))/(F(k-1)-S(k-1,k));
        end
        if ( j~=k & ( S(j+1,k)>F(j+1) | S(j+1,k)<F(j)  ) )
            %S(j+1,k) = .5 * (F(j)+F(j+1));
            S(j+1,k) = 0.75*F(j) + 0.25*F(j+1);
        end
    end
    % set transition probabilities
    for j = 1:k
        p(j,k) = (F(j)-S(j,k)) / (S(j+1,k)-S(j,k));
    end;
    % Arrow-Debreau prices
    L(1,k) = exp(-r*Dt) * L(1,k-1) * (1-p(1,k));
    for j = 2:k
        L(j,k) = exp(-r*Dt) * ( L(j,k-1)*(1-p(j,k)) + L(j-1,k-1)*p(j-1,k) );
    end
    L(k+1,k) = exp(-r*Dt) * L(k,k-1) * p(k,k);
end
S = S';
p = p';
L = L';