%   Copyright (c) 2001-2009 by Ales Cerny

%************************************************************************%
% chapter7sect4.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 chapter7sect4.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 binomial option pricing model using Discrete
% Fourier Transform.

clear;                    % clear the workspace 
clc;                      % clear screen        

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

%***************************%
%     Hedging Parameters    %
%***************************%
T = 3*Month;                                               % Time to maturity    
RehedgeInterval = 1*Minute;                                % Trading period      
S0 = 5100;                                                 % Initial stock price 
strike = 5355;

%***************************%
%    Transformation of      %
%       log returns         %
%***************************%
UnitTime = Month;
R1safe = 1.0033;                            % monthly safe return                
R1 = [1.053 0.965];                         % monthly return                     
PDistr = [0.5 0.5];                         % prob. density of monthly  returns  
lnR1 = log(R1);                             % monthly log return                 
mu1 = lnR1 * PDistr';                       % expected monthly log return     
sig1= sqrt(((lnR1-mu1).^2)*PDistr');        % volatility of monthly log return

dt = RehedgeInterval/UnitTime;
lnRdt = mu1*dt+(lnR1-mu1)*sqrt(dt);
Rdt=exp(lnRdt);
Rdtsafe=R1safe^dt; 
%************************%
%      Risk-neutral      %
%      probabilities     %
%************************%
QDistr=[Rdtsafe-Rdt(2) Rdt(1)-Rdtsafe]./(Rdt(1)-Rdt(2));

%************************%
%   grid indexation      %
%************************%
Tidx=ceil(T/RehedgeInterval)+1;           % Number of trading dates            
dlnS=lnRdt(1)-lnRdt(2);                   % increment on log price grid        
highlnRdt=lnRdt(1);                       % the highest return 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)*dlnS

%************************%
%      option payoff     %
%************************%
S_T = log(S0)+(Tidx-1)*highlnRdt-(0:(Tidx-1))*dlnS;  % log price at maturity        
S_T = exp(S_T');                                     % stock price at maturity
C_T = max([(S_T-strike)';zeros(1,length(S_T))]);     % option payoff at maturity    
C_T = min([C_T ; 2*strike*ones(1,length(S_T))]);     % truncation of high payoffs to avoid numerical instability of dft

%***********************%
%         pricing       %
%***********************%
tic;
b=[(QDistr'/Rdtsafe) ; zeros(Tidx-2,1)]';  % pricing kernel
fftsize=nextn(length(b));
CharFction=fft(b,fftsize);                 %fftsize must be chosen judiciously so that the transform is fast
Cimage = ifft(C_T,fftsize);
C = fft(Cimage.*(CharFction.^(Tidx-1)));
Cnext = fft(Cimage.*(CharFction.^(Tidx-2)));
elapsedtime = toc;

disp( '_____________________________________________________');
disp(' ');
disp(sprintf('trading frequency in minutes       %12.0f  ', RehedgeInterval));
disp(' '); 
disp(sprintf('binomial price                     %12.6f  ',real(C(1))));
disp(sprintf('Black-Scholes price                %12.6f  ',BSCall(S0,strike,log(R1safe),T/UnitTime,sig1)));
disp(' '); 
disp(sprintf('binomial delta                     %12.9f  ',real((Cnext(1)-Cnext(2))/S0/(Rdt(1)-Rdt(2)))));
disp(sprintf('Black-Scholes delta                %12.9f  ',normcdf((log(S0/strike)+(sig1^2/2+log(R1safe))*T/UnitTime)/(sig1*sqrt(T/UnitTime)))));
disp(' ');
%format /rd 12,3; 
disp(sprintf('required time in seconds           %12.3f  ',elapsedtime));
disp( '_____________________________________________________');
disp(' ');