% Copyright (c) 2001-2009 by Ales Cerny

%************************************************************************%
% chapter13sect3c.m - supplementary program to                            %
% Ales Cerny (2009) Mathematical Techniques in Finance (2nd ed.)         % 
% Princeton University Press http://press.princeton.edu/titles/9079.html %
%************************************************************************%
    
% This code is provided 'as-is', without any express or implied warranty.  
% 
% Permission is granted to anyone to use this code for any purpose, 
% subject to the following restrictions:
% 
% 1. The origin of this code must not be misrepresented; you must not
%    claim that you wrote the original code.
% 2. Modified code versions must be plainly marked as such, and must not 
%    be misrepresented as being the original code.
% 3. This notice may not be removed from any source distribution.

% NOTICE TO STUDENTS: To avoid accusations of plagiarism, if you use this 
% code or its modifications in assessed work you should prepend it with a 
% note stating: 
%   "This is the original/modified version of the code chapter13sect3C.m by 
%    Ales Cerny (2009), Mathematical Techniques in Finance (2nd ed.),  
%    Princeton University Press. The original version is available from 
%    http://www.martingales.info/mtfweb2".
% A similar acknowledgement should appear prominently inside your written 
% report.

% This program uses backward recursion to evaluate squared hedging error
% to maturity in a multinomial model calibrated using empirical FTSE100 
% equity return distribution.

clear;
clc;

%***************************%
%      trading time         %
%       parameters          %
%***************************%
Minute = 1;
Hour = 60;
HoursInDay  =  8;
DaysInWeek = 5;
DaysInMonth = 21;
Day  = HoursInDay*Hour;
Week = DaysInWeek*Day;
Month = DaysInMonth*Day;
Year = 12*Month;

%***************************%
%                           %
%     Hedging Parameters    %
%                           %
%***************************%

T = 33*Day;                                       % Time to maturity    
dt=60*Minute;                                      % Rehedging frequency 
S0=3567;                                          % 3576 Starting stock price, value 26/02/2003
Tidx=ceil(T/dt)+1;                                % Number of rebalancing dates
strike=3625;                                      % Strike 
Rsafe=1.04^(dt/Year);
mu = 0.0*dt/Year;

%***************************%
%                           %
%        Market             %
%                           %
%***************************%
%*********************************************************%
%             empirical distribution density              %
%*********************************************************%

lnR = [  0.0361547  0.0346407  0.0331266  0.0316125  0.0300985  0.0285844  0.0270704  0.0255563  0.0240423  0.0225282... 
         0.0210142  0.0195001  0.0179861  0.0164720  0.0149580  0.0134439  0.0119299  0.0104158  0.0089018  0.0073877... 
         0.0058737  0.0043596  0.0028455  0.0013315 -0.0001826 -0.0016966 -0.0032107 -0.0047247 -0.0062388 -0.0077528...
        -0.0092669 -0.0107809 -0.0122950 -0.0138090 -0.0153231 -0.0168371 -0.0183512 -0.0198652 -0.0213793 -0.0228934...
        -0.0244074 -0.0259215 -0.0274355 -0.0289496 -0.0304636 -0.0319777 -0.0334917 -0.0350058 -0.0365198 -0.0380339]'; 
 
PDistr=[ 0.0000136  0.0000000  0.0000068  0.0000000  0.0000339  0.0000136  0.0000136  0.0000136  0.0000136  0.0000272...
         0.0000136  0.0000679  0.0000950  0.0000475  0.0000883  0.0001358  0.0002648  0.0002037  0.0004956  0.0009098...
         0.0018331  0.0071627  0.0341976  0.2181871  0.5790646  0.1221527  0.0245500  0.0063140  0.0019349  0.0008962...
         0.0004888  0.0001833  0.0001222  0.0000883  0.0000950  0.0000883  0.0000543  0.0000339  0.0000136  0.0000068...
         0.0000000  0.0000272  0.0000000  0.0000000  0.0000068  0.0000068  0.0000000  0.0000272  0.0000000  0.0000068]'; 
n = length(lnR);
dlnR = lnR(1)-lnR(2);

%*************************************************************%
%   population moments from the empirical density function    %
%*************************************************************%
R = exp(lnR);                                       % risky return 
ER=PDistr'* R;                                      % mean return E[R] 
ElnR=PDistr'* lnR;                                  % mean logreturn E[lnR]    
sig=sqrt(PDistr'* (R.^2) - ER^2);                   % standard deviation of risky return
siglnR=sqrt(PDistr'*((lnR-ElnR).^2));               % standard deviation of log risky return
skew=((R-ER).^3)'* PDistr/(sig^3);                  % skewness of risky return
skewlnR=((lnR-ElnR).^3)'* PDistr/(siglnR^3);        % skewness of log risky return
kurt=((R-ER).^4)'* PDistr/(sig^4);                  % kurtosis of risky return
kurtlnR=((lnR-ElnR).^4)'* PDistr/(siglnR^4);        % kurtosis of log risky return

X = R - Rsafe;                                      % excess return 
EX= PDistr'* X;                                     % E[X] 
EX2=PDistr'*(X.^2);                                 % E[X^2]
a=EX/EX2;                                           % coefficient a
b=1-EX^2/EX2;                                       % coefficient b
QDistr=PDistr.*((1-a*X)/b);                         % risk-neutral probabilities as a vector
EQR=R'* QDistr;                                     % moments of X and R under V-O risk-neutral measure Q %  
EQlnR=lnR'* QDistr;
sigQ=sqrt(((R-EQR).^2)'* QDistr);
sigQlnR=sqrt(((lnR-EQlnR).^2)'*QDistr);
skewQ=(((R-EQR).^3)'*QDistr)/sigQ^3;               
skewQlnR=(((lnR-EQlnR).^3)'* QDistr)/sigQlnR^3;
kurtQ=(((R-EQR).^4)'* QDistr)/sigQ^4;
kurtQlnR=(((lnR-EQlnR).^4)'* QDistr)/sigQlnR^4;

% format /rd 16,5;  
disp( '_______________________________________________________________________');
disp(' ');
disp('                             empirical density                          ');
disp(' ');
disp('   population moment            under P                     under Q ');
disp(' ');
disp(sprintf('   E[R]               %16.5f            %16.5f ', ER , EQR));
disp(sprintf('   StDev[R]           %16.5f            %16.5f ', sig, sigQ));
disp(sprintf('   Skew[R]            %16.5f            %16.5f ', skew,skewQ));
disp(sprintf('   Kurt[R]            %16.5f            %16.5f ', kurt,kurtQ));
disp(' ');
disp(sprintf('   E[lnR]             %16.5f            %16.5f ', ElnR , EQlnR));
disp(sprintf('   StDev[lnR]         %16.5f            %16.5f ', siglnR, sigQlnR));
disp(sprintf('   Skew[lnR]          %16.5f            %16.5f ', skewlnR,skewQlnR));
disp(sprintf('   Kurt[lnR]          %16.5f            %16.5f ', kurtlnR,kurtQlnR));
disp(' ');
disp(sprintf('   annualized E[lnR]  %16.5f            %16.5f', [ElnR EQlnR]*Year/dt));
disp(sprintf('annualized StDev[lnR] %16.5f            %16.5f', [siglnR sigQlnR]*sqrt(Year/dt)));
disp(' ');
disp(' ');

starttime = clock;
%************************%
%   grid indexation      %
%************************%
highlnR = lnR(1);
dlnS = dlnR;

% there are 1+(n-1)*(tt-1) live cells at time tt, highest stock price at the top 
% log price at cell  1 at time tt is ln(S0)+(tt-1)*highlnRdt
% log price at cell ii at time tt is ln(S0)+(tt-1)*highlnRdt-(ii-1)*dlnS

%************************%
%      option payoff     %
%************************%
MaxIdx = 1+(n-1)*(Tidx-1);
S_T = log(S0)+(Tidx-1)*highlnR-(0:(MaxIdx-1))*dlnS;                     % log price at maturity                              
S_T = exp(S_T');                                                        % stock price at maturity  
strikeidx = floor((log(S0)-log(strike)+(Tidx-1)*highlnR)/dlnS)+1;
H = min([100*S0*ones(1,MaxIdx);max([(S_T-strike)';zeros(1,MaxIdx)])])'; % option payoff at maturity  
k_L(1) = 1; 
k_D(1) = 1;
for i = 2:Tidx
    k_L(i) = k_L(i-1)*Rsafe^2;
    k_D(i) = k_D(i-1)*b*Rsafe^2;
end
k_L=fliplr(k_L)';                                        % process k for loc. optimal strategy 
k_D=fliplr(k_D)';                                        % process k for dyn. optimal strategy
eps2_D=zeros(MaxIdx,1);
%************************%
%       main loop        %
%************************%
for tt = Tidx-1 :-1 : 1
    Hnext=H;
    epsnext=eps2_D;                                                    %mean value in next period  
    for ii = 1 : 1+(n-1)*(tt-1)
        focus=Hnext(ii:ii+n-1);
        H(ii)=(QDistr'*focus)./Rsafe;                                    % risk-neutral pricing     
        hedge=(PDistr.*X)'*(focus-Rsafe*H(ii))...
               /EX2;                                                   % optimal hedge           
        HedgeError=focus-hedge*X-Rsafe*H(ii);                          % hedging error              
        ESRE=(HedgeError.*PDistr)'*HedgeError;                          % ESRE                        
        eps2_D(ii)=PDistr'*epsnext(ii:ii+n-1)...
                 +k_D(tt+1)*ESRE;
    end
    disp(sprintf('period     %6.0f', tt))
end
finishtime = clock;
lattime=etime(finishtime, starttime);

%*********************************%
% computation of Toft's function  %
%    E_0[(S_t^2 * gamma_t)^2]     %
%      for t=0,1,2,...,Tidx-1     %
%*********************************%

ttime = (0:1:Tidx-2)';                    % time points from 0 to Tidx-1 

lambda = log(S0/strike)+log(Rsafe)*(Tidx-1-ttime)+mu*ttime;
ToftCoef = S0^2*exp(2*mu*ttime)./(2*pi*sig^2*sqrt((Tidx-1)^2-ttime.^2)).*...
           exp(-((lambda+sig^2/2*(Tidx-1)).^2+2*sig^2*ttime.*lambda)./(sig^2*(Tidx-1+ttime)));

ToftError = sig^4*(kurt-1+4*mu/sig*skew+4*(mu/sig)^2+(mu/sig)^4)*(ToftCoef'*k_D(2:Tidx))/4;
sig1=siglnR*sqrt(Year/dt);
BSPrice=BSCall(S0,strike,log(1.04),T/Year,sig1);

%*******************%
%      output       %
%*******************%
% format /rd 6,2;
disp(' ');
disp(' ');
disp(sprintf('squared error                                 %6.2f    ', eps2_D(1)));
disp(sprintf('runtime in seconds                            %6.2f    ', lattime));
disp(sprintf('Toft formula with empirical distr. parameters %6.2f    ', ToftError));
disp(sprintf('mean value based on empirical distribution    %6.2f    ', H(1)));
disp(sprintf('Black-Scholes price                           %6.2f    ', BSPrice ));
disp(sprintf('delta based on empirical distribution         %6.2f    ', hedge/S0));
disp(sprintf('Black-Scholes delta                           %6.2f    ', normcdf((log(S0/strike)+(sig1^2/2+log(1.04))*T/Year)/(sig1*sqrt(T/Year)))));