%% Revisiting Local Expansions for 0DTE Option Pricing 
% Author: Matteo Molta
clear, close all, clc, warning off

%% Data for TTM and RF rate
% Underlying of interest: SPX, DAX or EuroStoxx50
underlying = 'SPX';
%underlying = 'DAX';
%underlying = 'EuroStoxx50';
str = strcat('Data',underlying,'\Maturities.xlsx');
TTM = readtable(str);
tau = TTM.TTM_days_/365;
N = length(tau);
r = TTM.RF/100;

%% Model calibration and performance evaluation
% RMSE_t for both the expansion-based and the BSM models
RMSE_EXP_t = zeros(N,1);
RMSE_BSM_t = zeros(N,1);

% Number of options and volumes
opt_number = zeros(N,1);
Volumes = zeros(14,1);

% Calibrated parameters
Params_EXP = zeros(N,8);
Params_BSM = zeros(N,4);

% Moneyness range
Moneyness_range = zeros(N,2);

% Skewness and kurtosis metrics
Skewness.Mkt = zeros(N,1);
Skewness.EXP = zeros(N,1);
Kurtosis.Mkt = zeros(N,1);
Kurtosis.EXP = zeros(N,1);

% RMSEs for different moneyness levels, first column for BSM, second for EDG/GC
m_RMSE_DOTMP = [];
m_RMSE_OTMP = [];
m_RMSE_ATM = [];
m_RMSE_OTMC = [];
m_RMSE_DOTMC = [];

% Definition of a commod grid to compute the average IV values across days
% For the same reason, values of Market, Edgeworth/Gram-Charlier and BSM 
% implied volatilities are interpolated on that grid 
xx = -8:1/64:5;
Npoints = length(xx);
MatrixIV.Mkt = zeros(N,Npoints);
MatrixIV.EXP = zeros(N,Npoints);
MatrixIV.BSM = zeros(N,Npoints);

% Expansion considered: Edgeworth or Gram-Charlier
%expansion = 'Edgeworth';
expansion = 'Gram-Charlier';

% Jump types: 1 for Merton, 2 for Kou, 3 for VG, 4 for NIG
JT = 4;

% Flag  variable: 1 for cleaning on the tails, 0 for no cleaning
% Cleaning is performed only for SPX. DAX and Euro Stoxx 50 tails are almost fine
tail_filter = 1;

for i = 1:N
    % Data reading
    [Data, Spot, NewVolumes] = DataCleaning(r(i),tau(i),i,tail_filter,underlying);
    opt_number(i) = length(Data(:,1));
    Volumes = Volumes + NewVolumes;
    Moneyness_range(i,:) = [Data(1,4), Data(end,4)];

    % Calibration
    [IV, Params, RMSE, residuals] = Calibration(Data, Spot, r(i), tau(i), TTM.Date(i), JT, underlying, expansion);

    % Collection of RMSEs for different moneyness levels
    m_RMSE_DOTMP = [m_RMSE_DOTMP; residuals(Data(:,4)<-4,:)];
    m_RMSE_OTMP = [m_RMSE_OTMP; residuals(Data(:,4)>=-4 & Data(:,4)<-2,:)];
    m_RMSE_ATM = [m_RMSE_ATM; residuals(Data(:,4)>=-2 & Data(:,4)<=2,:)];
    m_RMSE_OTMC = [m_RMSE_OTMC; residuals(Data(:,4)>2 & Data(:,4)<=4,:)];
    m_RMSE_DOTMC = [m_RMSE_DOTMC; residuals(Data(:,4)>4,:)];

    % Implied volatility interpolation just for graphical purposes
    MatrixIV.Mkt(i,:) = interp1(Data(:,4),Data(:,5),xx);
    MatrixIV.EXP(i,:) = interp1(Data(:,4),IV.EXP,xx);
    MatrixIV.BSM(i,:) = interp1(Data(:,4),IV.BSM,xx); 

    % Skewness and kurtosis metrics update
    [~,index_atm] = min(abs(Data(:,4)));
    Skewness.Mkt(i) = interp1(Data(:,4),Data(:,5),1,'spline')-interp1(Data(:,4),Data(:,5),-1,'spline'); 
    Skewness.EXP(i) = interp1(Data(:,4),IV.EXP,1,'spline')-interp1(Data(:,4),IV.EXP,-1,'spline');
    Kurtosis.Mkt(i) = interp1(Data(:,4),Data(:,5),1,'spline')+interp1(Data(:,4),Data(:,5),-1,'spline')-2*Data(index_atm,5); 
    Kurtosis.EXP(i) = interp1(Data(:,4),IV.EXP,1,'spline')+interp1(Data(:,4),IV.EXP,-1,'spline')-2*IV.EXP(index_atm);

    % Calibrated parameters update
    Params_EXP(i,:) = Params.EXP;
    Params_BSM(i,:) = Params.BSM;

    % RMSE_t update
    RMSE_EXP_t(i) = RMSE.EXP;
    RMSE_BSM_t(i) = RMSE.BSM;
end

%% Model RMSEs and plots
RMSE_EXP = mean(RMSE_EXP_t)
RMSE_BSM = mean(RMSE_BSM_t)

% Plot of the ratio of the RMSEs over time
figure()
plot(categorical(TTM.Date),RMSE_EXP_t./RMSE_BSM_t,'LineWidth',1.2)
hold on
yline(RMSE_EXP/RMSE_BSM,'--','LineWidth',1.2)
title(strcat('Pricing error: RMSE(',expansion,') / RMSE(BSM)'))
set(gca, 'XTick', categorical(TTM.Date(3:5:30)))
set(findall(gca,'Type','Line'),'Linewidth',2)
set(gca,'Fontsize',16)

% RMSEs for different moneyness levels
m_RMSE_BSM = [median(m_RMSE_DOTMP(:,1)), median(m_RMSE_OTMP(:,1)), median(m_RMSE_ATM(:,1)), ...
    median(m_RMSE_OTMC(:,1)), median(m_RMSE_DOTMC(:,1))];
m_RMSE_EXP = [median(m_RMSE_DOTMP(:,2)), median(m_RMSE_OTMP(:,2)), median(m_RMSE_ATM(:,2)), ...
    median(m_RMSE_OTMC(:,2)), median(m_RMSE_DOTMC(:,2))];
TableData = [m_RMSE_BSM; m_RMSE_EXP; m_RMSE_EXP./m_RMSE_BSM];
T = array2table(TableData, 'VariableNames', {'DOTMP', 'OTMP', 'ATM', 'OTMC', 'DOTMC'}, 'RowNames', {'BSM', expansion, strcat(expansion,'/BSM')});
disp('Table: Pricing performance on different moneyness levels')
disp(T)

% Other plots for the analysis
PlotVisual(TTM.Date, opt_number, Params_EXP, Moneyness_range, Skewness, Kurtosis, MatrixIV, Volumes, expansion);
