%% Table 8
%% Date: 30/10/2025
%% Relation between lagged rr and bf to future intervention; But now, we decompose into left and right

cd(ini.dirs.code)
%% Initialize all data matrix
load ccydates;
T=size(ccydates,1);
%% intervention
load intervention; load interventiondates;
idx=ismember(ccydates,interventiondates);
int=zeros(T,1);
int(idx)=intervention;
int(isnan(int))=0;
lagint=lagmatrix(int,0:10);

load lagliborchange;load lagtelborchange;


%% RR & BF
load('lagrr'); load('lagbf1'); load('lagbf2');
load('bf25lag');load('bf10lag');load('rr25lag');load('rr10lag');
%% CDS
load('lagcds')

%% ccys
load ilsnyccy;
load ilslondonccy;
load NEERny;

%% Change in NEER
load lagNEER1
load lagNEER2
load lagNEER5
load lagNEER60
%load lagNEER3

%% Flows
load laginstflows
load lagforeignflows

%% vol smile

volsmile1w=NaN(T,5); %% 5 IV'S (-0.1,-0.25,0,.25,.1)

volsmile1w(:,1)=usdilsbf10(:,1)+usdilsatmv(:,1)-0.5*usdilsrr10(:,1);
volsmile1w(:,2)=usdilsbf25(:,1)+usdilsatmv(:,1)-0.5*usdilsrr25(:,1);
volsmile1w(:,3)=usdilsatmv(:,1);
volsmile1w(:,4)=usdilsbf25(:,1)+usdilsatmv(:,1)+0.5*usdilsrr25(:,1);
volsmile1w(:,5)=usdilsbf10(:,1)+usdilsatmv(:,1)+0.5*usdilsrr10(:,1);


%%% Create vol. smile for 1 month
volsmile1m=NaN(T,5); %% 5 IV'S (-0.1,-0.25,0,.25,.1)
volsmile1m(:,1)=usdilsbf10(:,2)+usdilsatmv(:,2)-0.5*usdilsrr10(:,2);
volsmile1m(:,2)=usdilsbf25(:,2)+usdilsatmv(:,2)-0.5*usdilsrr25(:,2);
volsmile1m(:,3)=usdilsatmv(:,2);
volsmile1m(:,4)=usdilsbf25(:,2)+usdilsatmv(:,2)+0.5*usdilsrr25(:,2);
volsmile1m(:,5)=usdilsbf10(:,2)+usdilsatmv(:,2)+0.5*usdilsrr10(:,2);

%%% Create vol. smile for 3 months
volsmile3m=NaN(T,5); %% 5 IV'S (-0.1,-0.25,0,.25,.1)
volsmile3m(:,1)=usdilsbf10(:,3)+usdilsatmv(:,3)-0.5*usdilsrr10(:,3);
volsmile3m(:,2)=usdilsbf25(:,3)+usdilsatmv(:,3)-0.5*usdilsrr25(:,3);
volsmile3m(:,3)=usdilsatmv(:,3);
volsmile3m(:,4)=usdilsbf25(:,3)+usdilsatmv(:,3)+0.5*usdilsrr25(:,3);
volsmile3m(:,5)=usdilsbf10(:,3)+usdilsatmv(:,3)+0.5*usdilsrr10(:,3);

%%% Create vol. smile for 6 month
volsmile6m=NaN(T,5); %% 5 IV'S (-0.1,-0.25,0,.25,.1)
volsmile6m(:,1)=usdilsbf10(:,4)+usdilsatmv(:,4)-0.5*usdilsrr10(:,4);
volsmile6m(:,2)=usdilsbf25(:,4)+usdilsatmv(:,4)-0.5*usdilsrr25(:,4);
volsmile6m(:,3)=usdilsatmv(:,4);
volsmile6m(:,4)=usdilsbf25(:,4)+usdilsatmv(:,4)+0.5*usdilsrr25(:,4);
volsmile6m(:,5)=usdilsbf10(:,4)+usdilsatmv(:,4)+0.5*usdilsrr10(:,4);

%%% Create vol. smile for 9 month
volsmile9m=NaN(T,5); %% 5 IV'S (-0.1,-0.25,0,.25,.1)
volsmile9m(:,1)=usdilsbf10(:,5)+usdilsatmv(:,5)-0.5*usdilsrr10(:,5);
volsmile9m(:,2)=usdilsbf25(:,5)+usdilsatmv(:,5)-0.5*usdilsrr25(:,5);
volsmile9m(:,3)=usdilsatmv(:,5);
volsmile9m(:,4)=usdilsbf25(:,5)+usdilsatmv(:,5)+0.5*usdilsrr25(:,5);
volsmile9m(:,5)=usdilsbf10(:,5)+usdilsatmv(:,5)+0.5*usdilsrr10(:,5);

%%% Create vol. smile for 12 month
volsmile12m=NaN(T,5); %% 5 IV'S (-0.1,-0.25,0,.25,.1)
volsmile12m(:,1)=usdilsbf10(:,6)+usdilsatmv(:,6)-0.5*usdilsrr10(:,6);
volsmile12m(:,2)=usdilsbf25(:,6)+usdilsatmv(:,6)-0.5*usdilsrr25(:,6);
volsmile12m(:,3)=usdilsatmv(:,6);
volsmile12m(:,4)=usdilsbf25(:,6)+usdilsatmv(:,6)+0.5*usdilsrr25(:,6);
volsmile12m(:,5)=usdilsbf10(:,6)+usdilsatmv(:,6)+0.5*usdilsrr10(:,6);

rightsmile=[.5*(volsmile1w(:,5)+volsmile1w(:,4)),.5*(volsmile1m(:,5)+volsmile1m(:,4)),.5*(volsmile3m(:,5)+volsmile3m(:,4)),.5*(volsmile6m(:,5)+volsmile6m(:,4)),.5*(volsmile9m(:,5)+volsmile9m(:,4)),.5*(volsmile12m(:,5)+volsmile12m(:,4))];
leftsmile=[.5*(volsmile1w(:,2)+volsmile1w(:,1)),.5*(volsmile1m(:,2)+volsmile1m(:,1)),.5*(volsmile3m(:,2)+volsmile3m(:,1)),.5*(volsmile6m(:,2)+volsmile6m(:,1)),.5*(volsmile9m(:,1)+volsmile9m(:,2)),.5*(volsmile12m(:,1)+volsmile12m(:,2))];

lead=10;
drightsmile=[repmat(NaN,5,6);rightsmile(lead+1:end,:)-rightsmile(1:end-lead,:)];
dleftsmile=[repmat(NaN,5,6);leftsmile(lead+1:end,:)-leftsmile(1:end-lead,:)];
drightsmilelag=lagmatrix(drightsmile,1);
dleftsmilelag=lagmatrix(dleftsmile,1);

res=NaN(8,6);
save('drightsmilelag','drightsmilelag');save('dleftsmilelag','dleftsmilelag');

%% Load crash probabilities
tend=readtable('..\_data\crash probabilities\tail prob (Zhang, Xiong 10D) with graphs.xlsx','Sheet','data');
tfd=readtable('..\_data\crash probabilities\tail prob (Zhang, Xiong 25D) with graphs','Sheet','data');

rise10d=[tend.X2Rise1M,tend.X3Rise3M,tend.X6Rise6M,tend.X9Drop9M,tend.X10Rise12M];
fall10d=[tend.X2Drop1M,tend.X3Drop3M,tend.X6Drop6M,tend.X9Drop9M,tend.X10Drop12M];

rise25d=[tfd.X2Rise1M,tfd.X3Rise3M,tfd.X6Rise6M,tfd.X9Rise9M,tfd.X10Rise12M];
fall25d=[tfd.X2Drop1M,tfd.X3Drop3M,tfd.X6Drop6M,tfd.X9Drop9M,tfd.X10Drop12M];

lead=10;
drise10d=[repmat(NaN,lead,5);rise10d(lead+1:end,:)-rise10d(1:end-lead,:)];
drise10dlag=lagmatrix(drise10d,1);
dfall10d=[repmat(NaN,lead,5);fall10d(lead+1:end,:)-fall10d(1:end-lead,:)];
dfall10dlag=lagmatrix(dfall10d,1);


drise25d=[repmat(NaN,lead,5);rise25d(lead+1:end,:)-rise25d(1:end-lead,:)];
drise25dlag=lagmatrix(drise25d,1);
dfall25d=[repmat(NaN,lead,5);fall25d(lead+1:end,:)-fall25d(1:end-lead,:)];
dfall25dlag=lagmatrix(drise25d,1);

%% Two pct 1 month
alldates=tend.dates;
idx=ismember(ccydates,datenum(alldates));
idx1=ismember(datenum(alldates),ccydates);
rise10d=[];
fall10d=[];
rise25d=[];
fall25d=[];

rise10d(idx,:)=drise10d(idx1,:);
fall10d(idx,:)=dfall10d(idx1,:);
rise25d(idx,:)=drise25d(idx1,:);
fall25d(idx,:)=dfall25d(idx1,:);

idx=ismember(ccydates,datenum(alldates));
idx1=ismember(datenum(alldates),ccydates);



load CITIindex.mat;load nsindex1.mat;load nsindex2.mat;
load ILshockindex; load CPIindex.mat; load lageurusd.mat;
load lagvix; load lagNEER10.mat;
load intweeklag;
%% dates
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/12/2019','dd/mm/yyyy'));

for j=2:6 %% run on all BF and all RR

    [coeff,tstat,r,n]=fastnw(lagint(idx:idx1,1)/1000,[CITIindex(idx:idx1),NSindex1(idx:idx1),ILshockindex(idx:idx1),CPIindex(idx:idx1),NSindex2(idx:idx1),lagint(idx:idx1,2)/1000,lageurusd(idx:idx1,2),lagvix(idx:idx1,1)-lagvix(idx:idx1,5),intweeklag(idx:idx1,2)/1000,lagNEER2(idx:idx1,1),lagNEER10(idx:idx1,1),0.5*fall10d(idx:idx1,j-1)+0.5*fall25d(idx:idx1,j-1),0.5*rise10d(idx:idx1,j-1)+0.5*rise25d(idx:idx1,j-1)],1); %% divide int by 1000.
    %% intercept
    res(1,j)=coeff(1);
    res(2,j)=tstat(1);

    %% left
    res(3,j)=coeff(end-1);
    res(4,j)=tstat(end-1);
    %% right
    res(5,j)=coeff(end);
    res(6,j)=tstat(end);
    res(end,j)=r;
end

