% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
% 
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license 
%
% This code is being provided solely for information and general 
% illustrative purposes. The authors will not be responsible for the 
% consequences of reliance upon using the code or for numbers produced 
% from using the code. 




% get calibrated data
load 'MODEL_Hist';

% Discretization parameters
NTime = 12;                    % number of time steps
NSim = 10000;                    % number of simulations
NBatches = 1;                  % number of batches

NHist = 251;                     % number of historic scenarios

NModels = 10;                  % number of models

% Options
%K = 6800;                                % Strike
%C = 1;                                  % call (1) or put (0)    
LF = 0; LC = 0.04; GF = 0; GC = 0.15;   % cliquet params
L = 80; U = 120;
T= 1;                                   % maturity of option

%option1 = 'Arithmetic Asian Call';
%option2 = 'Arithmetic Asian Put';
%exotic_price_c = @(x,y,r) Price_ArithmeticAsian(x,y,1,r,T);
%exotic_price_p = @(x,y,r) Price_ArithmeticAsian(x,y,0,r,T);

% exotic_price = @(x,r) Price_ArithmeticCliquet(x,LF,LC,GF,GC,r,T);
option1 = 'FixedStrike Lookback Call';
option2 = 'FixedStrike Lookback Put';
exotic_price_c = @(x,y,r) Price_LookBackFixedStrike(x,y,1,r,T);
exotic_price_p = @(x,y,r) Price_LookBackFixedStrike(x,y,0,r,T);

%option1 = 'FloatingStrike Lookback Call';
%option2 = 'FloatingStrike Lookback Put';
%exotic_price_c = @(x,y,r) Price_LookBackFloatingStrike(x,y,1,r,T);
%exotic_price_p = @(x,y,r) Price_LookBackFloatingStrike(x,y,0,r,T);

%exotic_price = @(x,r) Price_KnockOut(x,x,L,U,C,r,T);
vanilla_price_c = @(x,y,r) Price_CallPut(x,y,1,r,T);
vanilla_price_p = @(x,y,r) Price_CallPut(x,y,0,r,T);

% output
ep1c = zeros(NHist,NModels); %ep2c = ep1c; ep3c = ep1c;
ep1p = ep1c; %ep2p = ep1c; ep3p = ep1c;

vp_c = ep1c; vp_p = ep1c;
qp_c = ep1c; qp_p = ep1c;

vc = zeros(NHist,1); vp = vc;

% Black Scholes
for l = 1: NHist
    pathS = MC_B(Spot(l),r1y(l),0,T,BS_Par_Hist(l),NTime,NSim,NBatches);
    ep1c(l,1) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,1) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,1) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,1) = vanilla_price_p(pathS,Spot(l),r1y(l));
    [vc(l), vp(l)] = blsprice(Spot(l), Spot(l), r1y(l), T, BS_Par_Hist(l),0);
    qp_c(l,1) = ep1c(l,1)/vp_c(l,1);
    qp_p(l,1) = ep1p(l,1)/vp_p(l,1);
end

% Merton
for l = 1: NHist
    pathS = MC_M(Spot(l),r1y(l),0,T,...
        MJ_Par_Hist(l,1),MJ_Par_Hist(l,2),MJ_Par_Hist(l,3),MJ_Par_Hist(l,4),...
        NTime,NSim,NBatches);
    ep1c(l,2) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,2) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,2) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,2) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,2) = ep1c(l,2)/vp_c(l,2);
    qp_p(l,2) = ep1p(l,2)/vp_p(l,2);
end

% Heston
for l = 1: NHist
    [pathS, pathV] = MC_QE(Spot(l),r1y(l),0,T,...
        HESTON_Par_Hist(l,1),HESTON_Par_Hist(l,2),HESTON_Par_Hist(l,3),...
        HESTON_Par_Hist(l,4),HESTON_Par_Hist(l,5),...
        NTime,NSim,NBatches);
    ep1c(l,3) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,3) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,3) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,3) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,3) = ep1c(l,3)/vp_c(l,3);
    qp_p(l,3) = ep1p(l,3)/vp_p(l,3);
end

% Bates
for l = 1: NHist
    [pathS, pathV] = MC_QE_j(Spot(l),r1y(l),0,T,...
        BATES_Par_Hist(l,1),BATES_Par_Hist(l,2),BATES_Par_Hist(l,4),BATES_Par_Hist(l,5),...
        BATES_Par_Hist(l,3),BATES_Par_Hist(l,7),BATES_Par_Hist(l,7),BATES_Par_Hist(l,8),...
        NTime,NSim,NBatches);
    ep1c(l,4) = exotic_price_c(pathS,Spot(l),r1y(l));
    ep1p(l,4) = exotic_price_p(pathS,Spot(l),r1y(l));
    vp_c(l,4) = vanilla_price_c(pathS,Spot(l),r1y(l));
    vp_p(l,4) = vanilla_price_p(pathS,Spot(l),r1y(l));
    qp_c(l,4) = ep1c(l,4)/vp_c(l,4);
    qp_p(l,4) = ep1p(l,4)/vp_p(l,4);
end

% % Variance Gamma CGM
% for l = 1:NHist
%     pathS  = MC_VG_CGM(Spot(l),r1y(l),0,T,...
%         VG_Par_Hist(l,1),VG_Par_Hist(l,2),VG_Par_Hist(l,3),...
%         NTime,NSim,NBatches);
%     ep1c(l,5) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,5) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,5) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,5) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,5) = ep1c(l,5)/vp_c(l,1);
%     qp_p(l,5) = ep1p(l,5)/vp_p(l,1);
% end
% 
% % Normal Inverse Gaussian
% for l = 1:NHist
%     pathS  = MC_NIG(Spot(l),r1y(l),0,T,...
%         NIG_Par_Hist(l,1),NIG_Par_Hist(l,2),NIG_Par_Hist(l,3),...
%         NTime,NSim,NBatches);
%     ep1c(l,6) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,6) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,6) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,6) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,6) = ep1c(l,1)/vp_c(l,1);
%     qp_p(l,6) = ep1p(l,6)/vp_p(l,1);
% end
% 
% % Variance Gamma - GOU
% for l = 1:NHist
%     pathS  = MC_VGGOU(Spot(l),r1y(l),0,T,...
%         VGGOU_Par_Hist(l,1), VGGOU_Par_Hist(l,2), VGGOU_Par_Hist(l,3),...
%         VGGOU_Par_Hist(l,6), VGGOU_Par_Hist(l,4), VGGOU_Par_Hist(l,5),...
%         NTime,NSim,NBatches);
%     ep1c(l,7) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,7) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,1) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,7) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,7) = ep1c(l,7)/vp_c(l,1);
%     qp_p(l,7) = ep1p(l,7)/vp_p(l,1);
% end
% 
% % % Variance Gamma - CIR
% for l = 1: NHist
%     pathS =  MC_VGCIR(Spot(l),r1y(l),0,T,...
%         VGCIR_Par_Hist(l,1),VGCIR_Par_Hist(l,2),VGCIR_Par_Hist(l,3), ...
%         VGCIR_Par_Hist(l,4),VGCIR_Par_Hist(l,5),VGCIR_Par_Hist(l,6),...
%         NTime,NSim,NBatches);
%     ep1c(l,8) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,8) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,8) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,8) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,8) = ep1c(l,8)/vp_c(l,1);
%     qp_p(l,8) = ep1p(l,8)/vp_p(l,1);
% end
% 
% % NIG - GOU
% for l = 1: NHist
%     pathS =  MC_NIGGOU(Spot(l),r1y(l),0,T,...
%         NIGGOU_Par_Hist(l,1),NIGGOU_Par_Hist(l,2),NIGGOU_Par_Hist(l,3), ...
%         NIGGOU_Par_Hist(l,4),NIGGOU_Par_Hist(l,5),NIGGOU_Par_Hist(l,6),NTime,NSim,NBatches);
%     ep1c(l,9) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,9) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,9) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,9) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,9) = ep1c(l,9)/vp_c(l,1);
%     qp_p(l,9) = ep1p(l,9)/vp_p(l,1);
% end
% 
% % NIG - CIR
% for l = 1: NHist
%     pathS =  MC_NIGCIR(Spot(l),r1y(l),0,T,...
%         NIGCIR_Par_Hist(l,1),NIGCIR_Par_Hist(l,2),NIGCIR_Par_Hist(l,3), ...
%         NIGCIR_Par_Hist(l,5),NIGCIR_Par_Hist(l,6),NIGCIR_Par_Hist(l,4),NTime,NSim,NBatches);
%     ep1c(l,10) = exotic_price_c(pathS,Spot(l),r1y(l));
%     ep1p(l,10) = exotic_price_p(pathS,Spot(l),r1y(l));
%     vp_c(l,10) = vanilla_price_c(pathS,Spot(l),r1y(l));
%     vp_p(l,10) = vanilla_price_p(pathS,Spot(l),r1y(l));
%     qp_c(l,10) = ep1c(l,10)/vp_c(l,1);
%     qp_p(l,10) = ep1p(l,10)/vp_p(l,1);
% end
m1 = 'BS'; m2 = 'MJ'; m3 = 'Heston'; m4 = 'Bates';

plotit3(ep1c,option1,m1, m2, m3, m4);
plotit3(ep1p,option2,m1, m2, m3, m4);
plotit3(qp_c,strcat('Normalized ',option1),m1, m2, m3, m4);
plotit3(qp_p,strcat('Normalized ',option2),m1, m2, m3, m4);

modelriskc = repmat(max(ep1c,[],2),1,NModels)-ep1c;
modelriskp = -repmat(max(ep1p,[],2),1,NModels)+ep1p;
