% Copyright (c) 2001-2009 by Ales Cerny

%************************************************************************%
% chapter13sect3.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 chapter13sect3.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 implements multinomial option pricing model. 
% The log returns are calibrated to give Brownian motion in the limit.
% The hedging errors are evaluated recursively in the stock price lattice

clear;                    % clear the workspace 
clc;                      % clear screen     

%***************************%
%      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 = 6*Week;                                                % Time to maturity    
RehedgeInterval = Hour;                                    % Trading period      
S0 = 5100;                                                 % Initial stock price 
strike = 5355;

%***************************%
%    Transformation of      %
%       log returns         %
%***************************%
UnitTime = Week;
R1safe   = 1.04 ^(1/52);                                    % weekly safe return
lnR1     = [0.060 0.040 0.02  0.000 -0.02 -0.04 -0.06 ];    % weekly return
PDistr   = [0.013 0.067 0.273 0.384  0.199 0.050 0.014];    % prob. density of weekly  returns
mu1      = lnR1*PDistr';                                    % expected weekly log return 
sig1     = sqrt(((lnR1-mu1).^2)*PDistr');                   % volatility of weekly log return 

dt = RehedgeInterval/UnitTime;
lnRdt=mu1*dt+(lnR1-mu1)*sqrt(dt);                           % log return over rehedging interval
Rdt=exp(lnRdt);
Rdtsafe=R1safe^dt; 

%************************%
%      Risk-neutral      %
%      probabilities     %
%************************%
X = Rdt - Rdtsafe;                               % excess return          
EX= X*PDistr';                                   % E[X]                   
EX2=(X.^2)*PDistr';                              % E[X^2]                 
sigX=sqrt(EX2-(EX)^2);                           % st. dev. of returns    
kurt=((X-EX).^4)*PDistr'/(sigX^4);               % kurtosis of returns    
a=EX/EX2;                                        % coefficient a          
b=1-EX^2/EX2;                                    % coefficient b          
QDistr=PDistr.*((1-a*X)/b);                      % r-n prob-s as a vector 

%************************%
%   grid indexation      %
%************************%
Tidx = ceil(T/RehedgeInterval)+1;      % Number of trading dates     
dlnR =lnRdt(1)-lnRdt(2);               % increment on log price grid
highlnRdt=lnRdt(1);                    % the highest return over one period %
n = length(lnRdt);                     % number of branches over one period %

% there are tt 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)*dlnR

%************************%
%      option payoff     %
%************************%
MaxDim=1+(n-1)*(Tidx-1);                     % no. of cells at time T 
S_T = log(S0)+(Tidx-1)*highlnRdt...
        -(0:MaxDim-1)*dlnR;                  % log price at maturity        
S_T=exp(S_T);                                % stock price at maturity
H=max([(S_T-strike); zeros(1,MaxDim)]);      % option payoff at maturity  
eps2_D=zeros(length(S_T),1);                 % time T value of sq. error 
k_D = (b*Rdtsafe^2).^(Tidx-1:-1:0);          % dynamically-optimal k  
%************************%
%       main loop        %
%************************%
tic;                                                     % start of computation    
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)/Rdtsafe;                    % risk-neutral pricing
        hedge=(PDistr.*X)*(focus-Rdtsafe*H(ii))/EX2;     % optimal hedge 
        HedgeError=focus'-hedge*X-Rdtsafe*H(ii);         % hedging error  
        ESRE=(HedgeError.*PDistr)*HedgeError';           % ESRE  
        eps2_D(ii)=PDistr*epsnext(ii:ii+n-1)...
                 +k_D(tt+1)*ESRE;
    end             
end    

disp( '_______________________________________________________________________');
disp(' ');
% format /rd 12,4; 
disp(sprintf('option''s mean value at t=0        %12.4f   ', H(1)));
disp(' '); 
disp(sprintf('sq. hedging error to maturity     %12.4f   ', eps2_D(1)));
disp(' ');
%format /rd 12,2; 
disp(sprintf('kurtosis of one-period return     %12.2f   ', kurt));
disp(' ');
disp(sprintf('number of trading dates           %12.2f   ', Tidx));
disp(' ');
%format /rd 12,2; 
disp(sprintf('required time in seconds          %12.1f    ', toc));
disp(' ');
disp( '_______________________________________________________________________');