%% Figure 4
%% Date: 30/10/2025

%% Change in USD ILS. I choose NY as that's what I halagatmve data for RR
%% Initialize table
cd(ini.dirs.code)
%% Initialize all data matrix
load ccydates; 
T=size(ccydates,1);
%% interlagatmvention
load intervention; load interventiondates;
load lagint;
idx=ismember(ccydates,interventiondates);
int=zeros(T,1);
int(idx)=intervention;
int(isnan(int))=0;
lagint=lagmatrix(int,0:100);
%% dates
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
%idx=find(ccydates==datenum('01/01/2016','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/12/2019','dd/mm/yyyy'));
%idx1=find(ccydates==datenum('30/12/2016','dd/mm/yyyy'));
load CITIindex.mat;load nsindex1.mat;load nsindex2.mat;
load ILshockindex; load CPIindex.mat; load lageurusd.mat;
load lagvix; load lagNEER10.mat;
load intweeklag;
load lagint;

%% RR & BF
load('lagrr'); load('lagbf');
%% CDS
load('lagcds')

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

%% Change in NEER
load lagNEER1
load lagNEER2
load lagNEER5
load lagNEER60
%load lagNEER3
%% Change in USD/ILS
load lagusdils1;

%% Flows
load laginstflows
load lagforeignflows
load lageurusd
load lagNEER10;
load intweeklag;

%%%% Exogenous
%% lagged interest rate change
%% 5-day change in interest rate
load lagtelborchange;load lagliborchange;
%% vix
load vix; load vixdates;
vixfull=NaN(T,1);

% Find matching dates using intersect (more reliable)
[~, idx_vix, idx_ccy] = intersect(vixdates, ccydates);

% Assign VIX values to matching dates
vixfull(idx_ccy) = vix(idx_vix);
vixfull = fillmissing(vixfull,'previous');

p=22*12+10;
lagvix=lagmatrix(vixfull,0:p);
save('lagvix','lagvix');
%%% dilsny

load optionsdates;load usdilsrr10;load usdilsrr25;
%% 1 day change RR and BF
lead=1;
drr1w10=[repmat(NaN,lead,6);usdilsrr10(1+lead:end,:)-usdilsrr10(1:end-lead,:)];
drr1w25=[repmat(NaN,lead,6);usdilsrr25(1+lead:end,:)-usdilsrr25(1:end-lead,:)];



rr25lag=lagmatrix(drr1w25,-150:0);
rr10lag=lagmatrix(drr1w10,-150:0);

[L,M]=size(rr25lag);
eyemat=zeros(M,M);
for i=1:M
    eyemat(M-i+1,i)=1;
end

rr25lag=rr25lag*eyemat;
rr10lag=rr10lag*eyemat;

load usdilsbf10;
load usdilsbf25;
load usdilsatmv;


%% 1-day
lead=1;
dbf1w10=[repmat(NaN,lead,6);usdilsbf10(1+lead:end,:)-usdilsbf10(1:end-lead,:)];
dbf1w25=[repmat(NaN,lead,6);usdilsbf25(1+lead:end,:)-usdilsbf25(1:end-lead,:)];
bf10lag=lagmatrix(dbf1w10,-150:0);
bf25lag=lagmatrix(dbf1w25,-150:0);
bf10lag=bf10lag*eyemat;
bf25lag=bf25lag*eyemat;
%% 1- day atmv
lead=1;
datmv=[repmat(NaN,lead,6);usdilsatmv(1+lead:end,:)-usdilsatmv(1:end-lead,:)];
lagatmv=lagmatrix(datmv,-150:0);
lagatrm=lagatmv*eyemat;

%% Zhang probabilities %%

%% 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=1;
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,:);

rise10dlag=lagmatrix(rise10d,-150:0);
rise25dlag=lagmatrix(rise25d,-150:0);
fall10dlag=lagmatrix(fall10d,-150:0);
fall25dlag=lagmatrix(fall25d,-150:0);


[L,M]=size(rise10dlag);
eyemat=zeros(M,M);
for i=1:M
    eyemat(M-i+1,i)=1;
end

rise10dlag=rise10dlag*eyemat;
rise25dlag=rise25dlag*eyemat;

fall10dlag=fall10dlag*eyemat;
fall25dlag=fall25dlag*eyemat;

%% dates
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/12/2019','dd/mm/yyyy'));

res=[];
lowbound10=[];
highbound10=[];

lowbound5=[];
highbound5=[];

T=length(intervention(idx:idx1));

z=[ones(T,1),lageurusd(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),NSindex1(idx:idx1),NSindex2(idx:idx1),ILshockindex(idx:idx1),CPIindex(idx:idx1),CITIindex(idx:idx1),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3)];
P_z=z/(z'*z)*z'; %% projection matrix

x=[ones(T,1),lagint(idx:idx1,1)/1000,lageurusd(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),NSindex1(idx:idx1),NSindex2(idx:idx1),ILshockindex(idx:idx1),CPIindex(idx:idx1),CITIindex(idx:idx1)];
xtag=P_z*x;

load ilsnyccy;


%% What happens to probabilities
%% 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];


alldates=tend.dates;
idx=ismember(ccydates,datenum(alldates));
idx1=ismember(datenum(alldates),ccydates);



%% dates
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/12/2019','dd/mm/yyyy'));


res=[];
lowbound10=[];
highbound10=[];

lowbound5=[];
highbound5=[];

T=length(intervention(idx:idx1));

z=[ones(T,1),lageurusd(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),NSindex1(idx:idx1),NSindex2(idx:idx1),ILshockindex(idx:idx1),CPIindex(idx:idx1),CITIindex(idx:idx1),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3)];
P_z=z/(z'*z)*z'; %% projection matrix

x=[ones(T,1),lagint(idx:idx1,1)/1000,lageurusd(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),NSindex1(idx:idx1),NSindex2(idx:idx1),ILshockindex(idx:idx1),CPIindex(idx:idx1),CITIindex(idx:idx1)];
xtag=P_z*x;

load ilsnyccy;
%% LP-IV: FALL
for j=2:6
    for k=1:101
        ret=100*(ilsnyccy(k+1:end)./ilsnyccy(1:end-k)-1);
        ret=ret';
        ret=[NaN(k,1);ret];
        totalfall=[0.5*(fall10d(k+1:end,j-1)-fall10d(1:end-k,j-1)+fall25d(k+1:end,j-1)-fall25d(1:end-k,j-1))];
        totalrise=[0.5*(rise10d(k+1:end,j-1)-rise10d(1:end-k,j-1)+rise25d(k+1:end,j-1)-rise25d(1:end-k,j-1))];
        totalfall=[NaN(k,1);totalfall];
        totalrise=[NaN(k,1);totalrise];
        lead=k;
        xtaglag=lagmatrix(xtag,k-1);
        [coeff,t,r,n]=fastnw(totalfall(idx:idx1),[xtaglag,ret(idx:idx1)],0,k-1);
        res(j-1,k)=coeff(2);
        lowbound10(j-1,k)=coeff(2)-1.65*coeff(2)/t(2);
        highbound10(j-1,k)=coeff(2)+1.65*coeff(2)/t(2);

        lowbound5(j-1,k)=coeff(2)-1.96*coeff(2)/t(2);
        highbound5(j-1,k)=coeff(2)+1.96*coeff(2)/t(2);
    end

end


tiledlayout(3,2);
nexttile
plot(0:100,res(1,1:101),'k','LineWidth',2); hold on;plot(0:100,lowbound5(1,1:101),'b--','LineWidth',1); hold on;plot(0:100,highbound5(1,1:101),'b--','LineWidth',1)
hold on; plot(0:100,lowbound10(1,1:101),'k--','LineWidth',1); plot(0:100,highbound10(1,1:101),'k--','LineWidth',1)
yline(0)
title('2% appreciation')
set(gca,'FontSize',17);
ylim([-0.6,0.1])

nexttile
plot(0:100,res(2,1:101),'k','LineWidth',2); hold on;plot(0:100,lowbound5(2,1:101),'b--','LineWidth',1); hold on;plot(0:100,highbound5(2,1:101),'b--','LineWidth',1)
hold on; plot(0:100,lowbound10(2,1:101),'k--','LineWidth',1); plot(0:100,highbound10(2,1:101),'k--','LineWidth',1)
yline(0)
title('3% appreciation')
set(gca,'FontSize',17);
ylim([-0.6,0.1])

nexttile
plot(0:100,res(3,1:101),'k','LineWidth',2); hold on;plot(0:100,lowbound5(3,1:101),'b--','LineWidth',1); hold on;plot(0:100,highbound5(3,1:101),'b--','LineWidth',1)
hold on; plot(0:100,lowbound10(3,1:101),'k--','LineWidth',1); plot(0:100,highbound10(3,1:101),'k--','LineWidth',1)
yline(0)
title('6% appreciation')
set(gca,'FontSize',17);
ylim([-0.6,0.1])

nexttile;
plot(0:100,res(4,1:101),'k','LineWidth',2); hold on;plot(0:100,lowbound5(4,1:101),'b--','LineWidth',1); hold on;plot(0:100,highbound5(4,1:101),'b--','LineWidth',1)
hold on; plot(0:100,lowbound10(4,1:101),'k--','LineWidth',1); plot(0:100,highbound10(4,1:101),'k--','LineWidth',1)
yline(0)
title('9% appreciation')
set(gca,'FontSize',17);
ylim([-0.6,0.1])

nexttile;
plot(0:100,res(5,1:101),'k','LineWidth',2); hold on;plot(0:100,lowbound5(5,1:101),'b--','LineWidth',1); hold on;plot(0:100,highbound5(5,1:101),'b--','LineWidth',1)
hold on; plot(0:100,lowbound10(5,1:101),'k--','LineWidth',1); plot(0:100,highbound10(5,1:101),'k--','LineWidth',1)
yline(0)
title('10% appreciation')
set(gca,'FontSize',17);
ylim([-0.6,0.1])


legend('Response of Appreciation Probabilities','95\% CI','','90\% CI','Box','off','Location','southoutside')
lgd = legend;
lgd.Layout.Tile = 6;

set(gca,'FontSize',17);


% Define common x vector for fill
x = [0:100, fliplr(0:100)];

% Define color for the confidence band
band_color = [0.6 0.6 1];

% Define x-axis tick spacing
xtick_spacing = 20;

% Tiled layout with custom padding and tile spacing
tiledlayout(3,2, 'TileSpacing', 'normal', 'Padding', 'normal');

% Plotting for 2% appreciation
nexttile
fill(x, [lowbound10(1,:), fliplr(highbound10(1,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(1,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('2% appreciation')
xticks(0:xtick_spacing:100);
set(gca,'FontSize',17);
ylim([-0.6,0.1]);

% Plotting for 3% appreciation
nexttile
fill(x, [lowbound10(2,:), fliplr(highbound10(2,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(2,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('3% appreciation')
xticks(0:xtick_spacing:100);
set(gca,'FontSize',17);
ylim([-0.6,0.1]);

% Plotting for 6% appreciation
nexttile
fill(x, [lowbound10(3,:), fliplr(highbound10(3,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(3,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('6% appreciation')
xticks(0:xtick_spacing:100);
set(gca,'FontSize',17);
ylim([-0.6,0.1]);

% Plotting for 9% appreciation
nexttile
fill(x, [lowbound10(4,:), fliplr(highbound10(4,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(4,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('9% appreciation')
xticks(0:xtick_spacing:100);
set(gca,'FontSize',17);
ylim([-0.6,0.1]);

% Plotting for 10% appreciation
nexttile
fill(x, [lowbound10(5,:), fliplr(highbound10(5,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(5,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('10% appreciation')
xticks(0:xtick_spacing:100);
set(gca,'FontSize',17);
ylim([-0.6,0.1]);

legend('90\% CI', 'Response of Appreciation Probabilities','Box','off','Location','southoutside','interpreter','latex')
lgd = legend;
lgd.Layout.Tile = 6;
set(gca,'FontSize',17);




%%% depreciation
%% dates
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/12/2019','dd/mm/yyyy'));

res=[];
lowbound10=[];
highbound10=[];
lowbound5=[];
highbound5=[];

% CORRECT: T should be the length of the sample period, not full ccydates
T = length(intervention(idx:idx1));

% Build z matrix with correct dimensions
z = [ones(T,1), ...
     lageurusd(idx:idx1,1), ...
     lagvix(idx:idx1,1)-lagvix(idx:idx1,6), ...
     NSindex1(idx:idx1), ...
     NSindex2(idx:idx1), ...
     ILshockindex(idx:idx1), ...
     CPIindex(idx:idx1), ...
     CITIindex(idx:idx1), ...
     lagint(idx:idx1,2)/1000, ...
     intweeklag(idx:idx1,2), ...
     lagNEER2(idx:idx1,2), ...
     lagNEER10(idx:idx1,3)];

P_z = z/(z'*z)*z'; %% projection matrix

x = [ones(T,1), ...
     lagint(idx:idx1,1)/1000, ...
     lageurusd(idx:idx1,1), ...
     lagvix(idx:idx1,1)-lagvix(idx:idx1,6), ...
     NSindex1(idx:idx1), ...
     NSindex2(idx:idx1), ...
     ILshockindex(idx:idx1), ...
     CPIindex(idx:idx1), ...
     CITIindex(idx:idx1)];

xtag = P_z*x;

load ilsnyccy;

%% LP-IV: RISE (or FALL)
for j=2:6
    for k=1:101
        ret=100*(ilsnyccy(k+1:end)./ilsnyccy(1:end-k)-1);
        ret=ret';
        ret=[NaN(k,1);ret];
        
        totalfall=[0.5*(fall10d(k+1:end,j-1)-fall10d(1:end-k,j-1)+fall25d(k+1:end,j-1)-fall25d(1:end-k,j-1))];
        totalrise=[0.5*(rise10d(k+1:end,j-1)-rise10d(1:end-k,j-1)+rise25d(k+1:end,j-1)-rise25d(1:end-k,j-1))];
        totalfall=[NaN(k,1);totalfall];
        totalrise=[NaN(k,1);totalrise];
        
        lead=k;
        xtaglag=lagmatrix(xtag,k-1);
        
        [coeff,t,r,n]=fastnw(totalrise(idx:idx1),[xtaglag,ret(idx:idx1)],0,k-1);
        
        res(j-1,k)=coeff(2);
        lowbound10(j-1,k)=coeff(2)-1.65*coeff(2)/t(2);
        highbound10(j-1,k)=coeff(2)+1.65*coeff(2)/t(2);
        lowbound5(j-1,k)=coeff(2)-1.96*coeff(2)/t(2);
        highbound5(j-1,k)=coeff(2)+1.96*coeff(2)/t(2);
    end
end

% Define common x vector for fill
x = [0:100, fliplr(0:100)];

% Define color for the confidence band
band_color = [0.6 0.6 1];

% Define x-axis tick spacing
xtick_spacing = 20;
ytick_spacing = 0.1;

% Tiled layout with custom padding and tile spacing
tiledlayout(3,2, 'TileSpacing', 'normal', 'Padding', 'normal');

% Plotting for 2% depreciation
nexttile
fill(x, [lowbound10(1,:), fliplr(highbound10(1,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(1,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('2% depreciation')
xticks(0:xtick_spacing:100);
yticks(-0.3:0.1:0.1);
set(gca,'FontSize',17);
ylim([-0.3,0.15]);

% Plotting for 3% depreciation
nexttile
fill(x, [lowbound10(2,:), fliplr(highbound10(2,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(2,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('3% depreciation')
xticks(0:xtick_spacing:100);
yticks(-0.3:0.1:0.1);
set(gca,'FontSize',17);
ylim([-0.3,0.15]);

% Plotting for 6% depreciation
nexttile
fill(x, [lowbound10(3,:), fliplr(highbound10(3,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(3,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('6% depreciation')
xticks(0:xtick_spacing:100);
yticks(-0.3:0.1:0.1);
set(gca,'FontSize',17);
ylim([-0.3,0.15]);

% Plotting for 9% depreciation
nexttile
fill(x, [lowbound10(4,:), fliplr(highbound10(4,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(4,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('9% depreciation')
xticks(0:xtick_spacing:100);
yticks(-0.3:0.1:0.1);
set(gca,'FontSize',17);
ylim([-0.3,0.15]);

% Plotting for 10% depreciation
nexttile
fill(x, [lowbound10(5,:), fliplr(highbound10(5,:))], band_color, 'LineStyle', 'none'); hold on;
plot(0:100,res(5,1:101),'b','LineWidth',2); 
yline(0, 'k', 'LineWidth', 2)
title('10% depreciation')
xticks(0:xtick_spacing:100);
yticks(-0.3:0.1:0.1);
set(gca,'FontSize',17);
ylim([-0.3,0.15]);

legend('90\% CI', 'Response of depreciation Probabilities','Box','off','Location','southoutside','interpreter','latex')
lgd = legend;
lgd.Layout.Tile = 6;
set(gca,'FontSize',17);


% Save the current figure
fig = gcf;  % Get current figure

output_dir = '..\_results\Figure 4';
if ~exist(output_dir, 'dir')
    mkdir(output_dir);
end

% Save as PDF
exportgraphics(fig, fullfile(output_dir, 'depreciation.pdf'), 'ContentType', 'vector', 'Resolution', 300);
fprintf('Figure saved to: %s\n', fullfile(output_dir, 'depreciations.pdf'));
