%% Table 9
%% Date: 11/10/2022
%% Relation between lagged prob. and intervention.
cd(ini.dirs.code)
%% Initialize all data matrix
load ccydates; 
T=size(ccydates,1);
%% intervention
load intervention; load interventiondates;
idx=ismember(ccydates,interventiondates);
int=NaN(T,1);
int(idx)=intervention;
int(isnan(int))=0;
intlag=lagmatrix(int,0:100);


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

%% RR & BF
load('lagrr'); load('lagbf'); load('lagatmv');

%% CDS
load('lagcds')

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

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


%% Flows
load laginstflows
load lagforeignflows


%%%% Exogenous
%% lagged interest rate change
%% 5-day change in interest rate
load lagtelborchange;load lagliborchange;
%% vix
load vix; load vixdates;
vixfull=NaN(T,1);
idxtilde=ismember(vixdates,ccydates);
idxtilde1=ismember(ccydates,vixdates);
vixfull(idxtilde1)=vix(idxtilde);
vixfull=fillmissing(vixfull,'previous');
p=22*12+10;
lagvix=lagmatrix(vixfull,0:p);
save('lagvix','lagvix');
%%% dilsny

load optionsdates;load usdilsrr10;load usdilsrr25;


load optionsdates;load usdilsrr10;load usdilsrr25;
%% 10 day change RR and BF
lead=10;
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,0:10);
rr10lag=lagmatrix(drr1w10,0:10);

load usdilsbf10; load usdilsatmv;load usdilsrr10;load usdilsrr25;load usdilsbf25;
startdate='01/01/2013';
enddate='31/12/2019';
load ccydates
idx1=find(ccydates==datenum(startdate,'dd/mm/yyyy'));
idx2=find(ccydates==datenum(enddate,'dd/mm/yyyy'));



%%% RND %%%%

load('probl1_all');
load('probr1_all');

load('probl3_all');
load('probr3_all');

load('probl5_all');
load('probr5_all');

load('probl8_all');
load('probr8_all');

load('probl10_all');
load('probr10_all');


load('probl13_all');
load('probr13_all');

load('probl15_all');
load('probr15_all');

load('probl20_all');
load('probr20_all');

%%% RW probabilities %%%
load('probl1_all_p');
load('probr1_all_p');

load('probl3_all_p');
load('probr3_all_p');

load('probl5_all_p');
load('probr5_all_p');

load('probl8_all_p');
load('probr8_all_p');

load('probl10_all_p');
load('probr10_all_p');


load('probl13_all_p');
load('probr13_all_p');

load('probl15_all_p');
load('probr15_all_p');

load('probl20_all_p');
load('probr20_all_p');
%% Initialize all data matrix


%% Prob 1+3+5+8+10+13+15+20 change
probl1_all=probl1_all*100;
probr1_all=probr1_all*100;
probl3_all=probl3_all*100;
probr3_all=probr3_all*100;
probl5_all=probl5_all*100;
probr5_all=probr5_all*100;
probl8_all=probl10_all*100;
probr8_all=probr10_all*100;
probl10_all=probl10_all*100;
probr10_all=probr10_all*100;
probl13_all=probl13_all*100;
probr13_all=probr13_all*100;
probl15_all=probl15_all*100;
probr15_all=probr15_all*100;
probl20_all=probl20_all*100;
probr20_all=probr20_all*100;
%%% rwd
probl1_all_p=probl1_all_p*100;
probr1_all_p=probr1_all_p*100;
probl3_all_p=probl3_all_p*100;
probr3_all_p=probr3_all_p*100;
probl5_all_p=probl5_all_p*100;
probr5_all_p=probr5_all_p*100;
probl8_all_p=probl10_all_p*100;
probr8_all_p=probr10_all_p*100;
probl10_all_p=probl10_all_p*100;
probr10_all_p=probr10_all_p*100;
probl13_all_p=probl13_all_p*100;
probr13_all_p=probr13_all_p*100;
probl15_all_p=probl15_all_p*100;
probr15_all_p=probr15_all_p*100;
probl20_all_p=probl20_all_p*100;
probr20_all_p=probr20_all_p*100;


lead=1;
dprobl1=[repmat(NaN,lead,6); probl1_all(1+lead:end,:)-probl1_all(1:end-lead,:)];
dprobr1=[repmat(NaN,lead,6); probr1_all(1+lead:end,:)-probr1_all(1:end-lead,:)];

dprobl3=[repmat(NaN,lead,6); probl3_all(1+lead:end,:)-probl3_all(1:end-lead,:)];
dprobr3=[repmat(NaN,lead,6); probr3_all(1+lead:end,:)-probr3_all(1:end-lead,:)];

dprobl5=[repmat(NaN,lead,6);probl5_all(1+lead:end,:)-probl5_all(1:end-lead,:)];
dprobr5=[repmat(NaN,lead,6);probr5_all(1+lead:end,:)-probr5_all(1:end-lead,:)];


dprobl8=[repmat(NaN,lead,6);probl8_all(1+lead:end,:)-probl8_all(1:end-lead,:)];
dprobr8=[repmat(NaN,lead,6);probr8_all(1+lead:end,:)-probr8_all(1:end-lead,:)];

dprobl10=[repmat(NaN,lead,6);probl10_all(1+lead:end,:)-probl10_all(1:end-lead,:)];
dprobr10=[repmat(NaN,lead,6);probr10_all(1+lead:end,:)-probr10_all(1:end-lead,:)];


dprobl13=[repmat(NaN,lead,6);probl13_all(1+lead:end,:)-probl13_all(1:end-lead,:)];
dprobr13=[repmat(NaN,lead,6);probr13_all(1+lead:end,:)-probr13_all(1:end-lead,:)];


dprobl15=[repmat(NaN,lead,6);probl15_all(1+lead:end,:)-probl15_all(1:end-lead,:)];
dprobr15=[repmat(NaN,lead,6);probr15_all(1+lead:end,:)-probr15_all(1:end-lead,:)];


dprobl20=[repmat(NaN,lead,6);probl20_all(1+lead:end,:)-probl20_all(1:end-lead,:)];
dprobr20=[repmat(NaN,lead,6);probr20_all(1+lead:end,:)-probr20_all(1:end-lead,:)];

%% one day


%% one week
oneweekprobleft1_lag=lagmatrix(dprobl1(:,1),0:10);
oneweekprobright1_lag=lagmatrix(dprobr1(:,1),0:10);


oneweekprobleft3_lag=lagmatrix(dprobl3(:,1),0:10);
oneweekprobright3_lag=lagmatrix(dprobr3(:,1),0:10);

oneweekprobleft5_lag=lagmatrix(dprobl5(:,1),0:10);
oneweekprobright5_lag=lagmatrix(dprobr5(:,1),0:10);

oneweekprobleft10_lag=lagmatrix(dprobl10(:,1),0:10);
oneweekprobright10_lag=lagmatrix(dprobr10(:,1),0:10);

oneweekprobleft20_lag=lagmatrix(dprobl20(:,1),0:10);
oneweekprobright20_lag=lagmatrix(dprobr20(:,1),0:10);


%% one month

onemonthprobleft1_lag=lagmatrix(dprobl1(:,2),0:10);
onemonthprobright1_lag=lagmatrix(dprobr1(:,2),0:10);

onemonthprobleft3_lag=lagmatrix(dprobl3(:,2),0:10);
onemonthprobright3_lag=lagmatrix(dprobr3(:,2),0:10);

onemonthprobleft5_lag=lagmatrix(dprobl5(:,2),0:10);
onemonthprobright5_lag=lagmatrix(dprobr5(:,2),0:10);

onemonthprobleft8_lag=lagmatrix(dprobl8(:,2),0:10);
onemonthprobright8_lag=lagmatrix(dprobr8(:,2),0:10);


onemonthprobleft10_lag=lagmatrix(dprobl10(:,2),0:10);
onemonthprobright10_lag=lagmatrix(dprobr10(:,2),0:10);


onemonthprobleft13_lag=lagmatrix(dprobl13(:,2),0:10);
onemonthprobright13_lag=lagmatrix(dprobr13(:,2),0:10);

onemonthprobleft15_lag=lagmatrix(dprobl15(:,2),0:10);
onemonthprobright15_lag=lagmatrix(dprobr15(:,2),0:10);

onemonthprobleft20_lag=lagmatrix(dprobl20(:,2),0:10);
onemonthprobright20_lag=lagmatrix(dprobr20(:,2),0:10);

%% three month
threemonthprobleft1_lag=lagmatrix(dprobl1(:,3),0:10);
threemonthprobright1_lag=lagmatrix(dprobr1(:,3),0:10);

threemonthprobleft3_lag=lagmatrix(dprobl3(:,3),0:10);
threemonthprobright3_lag=lagmatrix(dprobr3(:,3),0:10);

threemonthprobleft5_lag=lagmatrix(dprobl5(:,3),0:10);
threemonthprobright5_lag=lagmatrix(dprobr5(:,3),0:10);

threemonthprobleft8_lag=lagmatrix(dprobl8(:,3),0:10);
threemonthprobright8_lag=lagmatrix(dprobr8(:,3),0:10);

threemonthprobleft10_lag=lagmatrix(dprobl10(:,3),0:10);
threemonthprobright10_lag=lagmatrix(dprobr10(:,3),0:10);

threemonthprobleft13_lag=lagmatrix(dprobl13(:,3),0:10);
threemonthprobright13_lag=lagmatrix(dprobr13(:,3),0:10);

threemonthprobleft15_lag=lagmatrix(dprobl15(:,3),0:10);
threemonthprobright15_lag=lagmatrix(dprobr15(:,3),0:10);


threemonthprobleft20_lag=lagmatrix(dprobl20(:,3),0:10);
threemonthprobright20_lag=lagmatrix(dprobr20(:,3),0:10);


sixmonthprobleft1_lag=lagmatrix(dprobl1(:,4),0:10);
sixmonthprobright1_lag=lagmatrix(dprobr1(:,4),0:10);

sixmonthprobleft3_lag=lagmatrix(dprobl3(:,4),0:10);
sixmonthprobright3_lag=lagmatrix(dprobr3(:,4),0:10);

sixmonthprobleft5_lag=lagmatrix(dprobl5(:,4),0:10);
sixmonthprobright5_lag=lagmatrix(dprobr5(:,3),0:10);

sixmonthprobleft8_lag=lagmatrix(dprobl8(:,4),0:10);
sixmonthprobright8_lag=lagmatrix(dprobr8(:,4),0:10);

sixmonthprobleft10_lag=lagmatrix(dprobl10(:,4),0:10);
sixmonthprobright10_lag=lagmatrix(dprobr10(:,4),0:10);


sixmonthprobleft13_lag=lagmatrix(dprobl13(:,4),0:10);
sixmonthprobright13_lag=lagmatrix(dprobr13(:,4),0:10);


sixmonthprobleft15_lag=lagmatrix(dprobl15(:,4),0:10);
sixmonthprobright15_lag=lagmatrix(dprobr15(:,4),0:10);

sixmonthprobleft20_lag=lagmatrix(dprobl20(:,4),0:10);
sixmonthprobright20_lag=lagmatrix(dprobr20(:,4),0:10);

ninemonthprobleft3_lag=lagmatrix(dprobl3(:,5),0:10);
ninemonthprobright3_lag=lagmatrix(dprobr3(:,5),0:10);


ninemonthprobleft5_lag=lagmatrix(dprobl1(:,5),0:10);
ninemonthprobright5_lag=lagmatrix(dprobr1(:,5),0:10);


ninemonthprobleft8_lag=lagmatrix(dprobl8(:,5),0:10);
ninemonthprobright8_lag=lagmatrix(dprobr8(:,5),0:10);


ninemonthprobleft10_lag=lagmatrix(dprobl10(:,5),0:10);
ninemonthprobright10_lag=lagmatrix(dprobr10(:,5),0:10);


ninemonthprobleft13_lag=lagmatrix(dprobl13(:,5),0:10);
ninemonthprobright13_lag=lagmatrix(dprobr13(:,5),0:10);

ninemonthprobleft15_lag=lagmatrix(dprobl15(:,5),0:10);
ninemonthprobright15_lag=lagmatrix(dprobr15(:,5),0:10);

ninemonthprobleft20_lag=lagmatrix(dprobl20(:,5),0:10);
ninemonthprobright20_lag=lagmatrix(dprobr20(:,5),0:10);


twelvemonthprobleft1_lag=lagmatrix(dprobl1(:,6),0:10);
twelvemonthprobright1_lag=lagmatrix(dprobr1(:,6),0:10);


twelvemonthprobleft3_lag=lagmatrix(dprobl3(:,6),0:10);
twelvemonthprobright3_lag=lagmatrix(dprobr3(:,6),0:10);


twelvemonthprobleft5_lag=lagmatrix(dprobl5(:,6),0:10);
twelvemonthprobright5_lag=lagmatrix(dprobr5(:,6),0:10);


twelvemonthprobleft8_lag=lagmatrix(dprobl8(:,6),0:10);
twelvemonthprobright8_lag=lagmatrix(dprobr8(:,6),0:10);


twelvemonthprobleft10_lag=lagmatrix(dprobl10(:,6),0:10);
twelvemonthprobright10_lag=lagmatrix(dprobr10(:,6),0:10);


twelvemonthprobleft13_lag=lagmatrix(dprobl13(:,6),0:10);
twelvemonthprobright13_lag=lagmatrix(dprobr13(:,6),0:10);


twelvemonthprobleft15_lag=lagmatrix(dprobl15(:,6),0:10);
twelvemonthprobright15_lag=lagmatrix(dprobr15(:,6),0:10);

twelvemonthprobleft20_lag=lagmatrix(dprobl20(:,6),0:10);
twelvemonthprobright20_lag=lagmatrix(dprobr20(:,6),0:10);

%%% RWD


lead=5;
dprobl1_p=[repmat(NaN,lead,6); probl1_all_p(1+lead:end,:)-probl1_all_p(1:end-lead,:)];
dprobr1_p=[repmat(NaN,lead,6); probr1_all_p(1+lead:end,:)-probr1_all_p(1:end-lead,:)];

dprobl3_p=[repmat(NaN,lead,6); probl3_all_p(1+lead:end,:)-probl3_all_p(1:end-lead,:)];
dprobr3_p=[repmat(NaN,lead,6); probr3_all_p(1+lead:end,:)-probr3_all_p(1:end-lead,:)];

dprobl5_p=[repmat(NaN,lead,6);probl5_all_p(1+lead:end,:)-probl5_all_p(1:end-lead,:)];
dprobr5_p=[repmat(NaN,lead,6);probr5_all_p(1+lead:end,:)-probr5_all_p(1:end-lead,:)];


dprobl8_p=[repmat(NaN,lead,6);probl8_all_p(1+lead:end,:)-probl8_all_p(1:end-lead,:)];
dprobr8_p=[repmat(NaN,lead,6);probr8_all_p(1+lead:end,:)-probr8_all_p(1:end-lead,:)];

dprobl10_p=[repmat(NaN,lead,6);probl10_all_p(1+lead:end,:)-probl10_all_p(1:end-lead,:)];
dprobr10_p=[repmat(NaN,lead,6);probr10_all_p(1+lead:end,:)-probr10_all_p(1:end-lead,:)];


dprobl13_p=[repmat(NaN,lead,6);probl13_all_p(1+lead:end,:)-probl13_all_p(1:end-lead,:)];
dprobr13_p=[repmat(NaN,lead,6);probr13_all_p(1+lead:end,:)-probr13_all_p(1:end-lead,:)];


dprobl15_p=[repmat(NaN,lead,6);probl15_all_p(1+lead:end,:)-probl15_all_p(1:end-lead,:)];
dprobr15_p=[repmat(NaN,lead,6);probr15_all_p(1+lead:end,:)-probr15_all_p(1:end-lead,:)];


dprobl20_p=[repmat(NaN,lead,6);probl20_all_p(1+lead:end,:)-probl20_all_p(1:end-lead,:)];
dprobr20_p=[repmat(NaN,lead,6);probr20_all_p(1+lead:end,:)-probr20_all_p(1:end-lead,:)];


%% one week
oneweekprobleft1_lag_p=lagmatrix(dprobl1_p(:,1),0:10);
oneweekprobright1_lag_p=lagmatrix(dprobr1_p(:,1),0:10);


oneweekprobleft3_lag_p=lagmatrix(dprobl3_p(:,1),0:10);
oneweekprobright3_lag_p=lagmatrix(dprobr3_p(:,1),0:10);

oneweekprobleft5_lag_p=lagmatrix(dprobl5_p(:,1),0:10);
oneweekprobright5_lag_p=lagmatrix(dprobr5_p(:,1),0:10);

oneweekprobleft10_lag_p=lagmatrix(dprobl10_p(:,1),0:10);
oneweekprobright10_lag_p=lagmatrix(dprobr10_p(:,1),0:10);

oneweekprobleft20_lag_p=lagmatrix(dprobl20_p(:,1),0:10);
oneweekprobright20_lag_p=lagmatrix(dprobr20_p(:,1),0:10);


%% one month

onemonthprobleft1_lag_p=lagmatrix(dprobl1_p(:,2),0:10);
onemonthprobright1_lag_p=lagmatrix(dprobr1_p(:,2),0:10);

onemonthprobleft3_lag_p=lagmatrix(dprobl3_p(:,2),0:10);
onemonthprobright3_lag_p=lagmatrix(dprobr3_p(:,2),0:10);

onemonthprobleft5_lag_p=lagmatrix(dprobl5_p(:,2),0:10);
onemonthprobright5_lag_p=lagmatrix(dprobr5_p(:,2),0:10);

onemonthprobleft8_lag_p=lagmatrix(dprobl8_p(:,2),0:10);
onemonthprobright8_lag_p=lagmatrix(dprobr8_p(:,2),0:10);


onemonthprobleft10_lag_p=lagmatrix(dprobl10_p(:,2),0:10);
onemonthprobright10_lag_p=lagmatrix(dprobr10_p(:,2),0:10);


onemonthprobleft13_lag_p=lagmatrix(dprobl13_p(:,2),0:10);
onemonthprobright13_lag_p=lagmatrix(dprobr13_p(:,2),0:10);

onemonthprobleft15_lag_p=lagmatrix(dprobl15_p(:,2),0:10);
onemonthprobright15_lag_p=lagmatrix(dprobr15_p(:,2),0:10);

onemonthprobleft20_lag_p=lagmatrix(dprobl20_p(:,2),0:10);
onemonthprobright20_lag_p=lagmatrix(dprobr20_p(:,2),0:10);

%% three month
threemonthprobleft1_lag_p=lagmatrix(dprobl1_p(:,3),0:10);
threemonthprobright1_lag_p=lagmatrix(dprobr1_p(:,3),0:10);

threemonthprobleft3_lag_p=lagmatrix(dprobl3_p(:,3),0:10);
threemonthprobright3_lag_p=lagmatrix(dprobr3_p(:,3),0:10);

threemonthprobleft5_lag_p=lagmatrix(dprobl5_p(:,3),0:10);
threemonthprobright5_lag_p=lagmatrix(dprobr5_p(:,3),0:10);

threemonthprobleft8_lag_p=lagmatrix(dprobl8_p(:,3),0:10);
threemonthprobright8_lag_p=lagmatrix(dprobr8_p(:,3),0:10);

threemonthprobleft10_lag_p=lagmatrix(dprobl10_p(:,3),0:10);
threemonthprobright10_lag_p=lagmatrix(dprobr10_p(:,3),0:10);

threemonthprobleft13_lag_p=lagmatrix(dprobl13_p(:,3),0:10);
threemonthprobright13_lag_p=lagmatrix(dprobr13_p(:,3),0:10);

threemonthprobleft15_lag_p=lagmatrix(dprobl15_p(:,3),0:10);
threemonthprobright15_lag_p=lagmatrix(dprobr15_p(:,3),0:10);


threemonthprobleft20_lag_p=lagmatrix(dprobl20_p(:,3),0:10);
threemonthprobright20_lag_p=lagmatrix(dprobr20_p(:,3),0:10);


sixmonthprobleft1_lag_p=lagmatrix(dprobl1_p(:,4),0:10);
sixmonthprobright1_lag_p=lagmatrix(dprobr1_p(:,4),0:10);

sixmonthprobleft3_lag_p=lagmatrix(dprobl3_p(:,4),0:10);
sixmonthprobright3_lag_p=lagmatrix(dprobr3_p(:,4),0:10);

sixmonthprobleft5_lag_p=lagmatrix(dprobl5_p(:,4),0:10);
sixmonthprobright5_lag_p=lagmatrix(dprobr5_p(:,3),0:10);

sixmonthprobleft8_lag_p=lagmatrix(dprobl8_p(:,4),0:10);
sixmonthprobright8_lag_p=lagmatrix(dprobr8_p(:,4),0:10);

sixmonthprobleft10_lag_p=lagmatrix(dprobl10_p(:,4),0:10);
sixmonthprobright10_lag_p=lagmatrix(dprobr10_p(:,4),0:10);


sixmonthprobleft13_lag_p=lagmatrix(dprobl13_p(:,4),0:10);
sixmonthprobright13_lag_p=lagmatrix(dprobr13_p(:,4),0:10);


sixmonthprobleft15_lag_p=lagmatrix(dprobl15_p(:,4),0:10);
sixmonthprobright15_lag_p=lagmatrix(dprobr15_p(:,4),0:10);

sixmonthprobleft20_lag_p=lagmatrix(dprobl20_p(:,4),0:10);
sixmonthprobright20_lag_p=lagmatrix(dprobr20_p(:,4),0:10);

ninemonthprobleft3_lag_p=lagmatrix(dprobl3_p(:,5),0:10);
ninemonthprobright3_lag_p=lagmatrix(dprobr3_p(:,5),0:10);


ninemonthprobleft5_lag_p=lagmatrix(dprobl1_p(:,5),0:10);
ninemonthprobright5_lag_p=lagmatrix(dprobr1_p(:,5),0:10);


ninemonthprobleft8_lag_p=lagmatrix(dprobl8_p(:,5),0:10);
ninemonthprobright8_lag_p=lagmatrix(dprobr8_p(:,5),0:10);


ninemonthprobleft10_lag_p=lagmatrix(dprobl10_p(:,5),0:10);
ninemonthprobright10_lag_p=lagmatrix(dprobr10_p(:,5),0:10);


ninemonthprobleft13_lag_p=lagmatrix(dprobl13_p(:,5),0:10);
ninemonthprobright13_lag_p=lagmatrix(dprobr13_p(:,5),0:10);

ninemonthprobleft15_lag_p=lagmatrix(dprobl15_p(:,5),0:10);
ninemonthprobright15_lag_p=lagmatrix(dprobr15_p(:,5),0:10);

ninemonthprobleft20_lag_p=lagmatrix(dprobl20_p(:,5),0:10);
ninemonthprobright20_lag_p=lagmatrix(dprobr20_p(:,5),0:10);


twelvemonthprobleft1_lag_p=lagmatrix(dprobl1_p(:,6),0:10);
twelvemonthprobright1_lag_p=lagmatrix(dprobr1_p(:,6),0:10);


twelvemonthprobleft3_lag_p=lagmatrix(dprobl3_p(:,6),0:10);
twelvemonthprobright3_lag_p=lagmatrix(dprobr3_p(:,6),0:10);


twelvemonthprobleft5_lag_p=lagmatrix(dprobl5_p(:,6),0:10);
twelvemonthprobright5_lag_p=lagmatrix(dprobr5_p(:,6),0:10);


twelvemonthprobleft8_lag_p=lagmatrix(dprobl8_p(:,6),0:10);
twelvemonthprobright8_lag_p=lagmatrix(dprobr8_p(:,6),0:10);


twelvemonthprobleft10_lag_p=lagmatrix(dprobl10_p(:,6),0:10);
twelvemonthprobright10_lag_p=lagmatrix(dprobr10_p(:,6),0:10);


twelvemonthprobleft13_lag_p=lagmatrix(dprobl13_p(:,6),0:10);
twelvemonthprobright13_lag_p=lagmatrix(dprobr13_p(:,6),0:10);


twelvemonthprobleft15_lag_p=lagmatrix(dprobl15_p(:,6),0:10);
twelvemonthprobright15_lag_p=lagmatrix(dprobr15_p(:,6),0:10);

twelvemonthprobleft20_lag_p=lagmatrix(dprobl20_p(:,6),0:10);
twelvemonthprobright20_lag_p=lagmatrix(dprobr20_p(:,6),0:10);

%% average atmv per horizon
[nanmean(usdilsatmv(idx1:idx2,1))*sqrt(7/365),nanmean(usdilsatmv(idx1:idx2,2))*sqrt(30/365),nanmean(usdilsatmv(idx1:idx2,3))*sqrt(30*3/365),nanmean(usdilsatmv(idx1:idx2,4))*sqrt(30*6/365),nanmean(usdilsatmv(idx1:idx2,5))*sqrt(30*9/365),nanmean(usdilsatmv(idx1:idx2,6))*sqrt(30*12/365)]

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

%% 1 Week

%% Probabilities
%[coeff,tstat,r,n]=fastnw(intervention(idx:idx1)/1000,[onemonthprobleft5_lag(:,2),onemonthprobright5_lag(:,2),lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagNEER60(idx:idx1,2),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagtelborchange(idx:idx1,2),lagliborchange(idx:idx1,2)],1,60); %% divide int by 1000.
%% We'll use GMM as the regressions show that there is noise.

load dleftsmilelag;
load drightsmilelag;
res=NaN(9,6);
%% GMM

%{
%% first stage
cd('gmm');
optset('cugmmest','method','hacc_p');
optset('cugmmest','method','serunc');

options=optimset('fminunc');

n=size(int(idx:idx1),1);
inst=[lagint(idx:idx1,2)/1000,lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagNEER60(idx:idx1,2),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagtelborchange(idx:idx1,2),lagliborchange(idx:idx1,2),drightsmilelag(idx:idx1,7),dleftsmilelag(idx:idx1,7)]; %% instruments
data=[int(idx:idx1)/1000,ones(n,1),oneweekprobleft3_lag(1:end-23,2),oneweekprobright3_lag(1:end-23,2),lagint(idx:idx1,2)/1000,lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagNEER60(idx:idx1,2),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagtelborchange(idx:idx1,2),lagliborchange(idx:idx1,2)]; %% dataset used for the estimation
temp=[data,inst];
temp=packr(temp);
data=temp(:,[1:12]);
inst=temp(:,[2 13:end]); %% we need the constant as another instrument
t=size(data,1);
maxlag=floor(4*(t/100)^(2/9));
%maxlag=60;
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
options=optimset('MaxFunEval',10^5);

[coeff,tstat,r,n]=fastnw(int(idx:idx1)/1000,[oneweekprobleft3_lag(1:end-23,2),oneweekprobright3_lag(1:end-23,2),lagint(idx:idx1,2)/1000,lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagNEER60(idx:idx1,2),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagtelborchange(idx:idx1,2),lagliborchange(idx:idx1,2)],1,60); %% divide int by 1000.
%startval=coeff;
startval=zeros(size(coeff,1),1);

[theta_final, S_final, J_test, probJ, bandw, var_theta, std_theta, conf_inter] = cugmmest(options, data, @popmomnew, startval, inst);
coeff=theta_final;
tstat=theta_final./std_theta;
[coeff,tstat]
j=1;
res(1,j)=coeff(1);
res(2,j)=tstat(1);
res(3,j)=coeff(2);
res(4,j)=tstat(2);
res(5,j)=coeff(3);
res(6,j)=tstat(3);
res(end-1,j)=J_test;
res(end,j)=probJ;
%}

%% 1 month

%% Probabilities

%% We'll use GMM as the regressions show that there is noise.

%% GMM

%% first stage
%[coeff,t,r,n]=fastnw(onemonthprobleft5_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)
%[coeff,t,r,n]=fastnw(onemonthprobright5_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)

load lagusdils60; load lagrr10; load lagbf1; load lagbf2; load lagNEER10; load lagint; load intweeklag;

n=size(int(idx:idx1),1);
inst=[lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)]; %% instruments
data=[int(idx:idx1)/1000,ones(n,1),onemonthprobleft5_lag(1:end-23,2),onemonthprobright5_lag(1:end-23,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2)]; %% dataset used for the estimation
[coeff,t]=fastnw(data(:,1),[data(:,2:end)],0);
temp=[data,inst];
temp=packr(temp);
data=temp(:,[1:11]);
inst=temp(:,[2 12:end]); %% we need the constant as another instrument
t=size(data,1);
maxlag=floor(4*(t/100)^(2/9));
maxlag=60;
%optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.

%options=optimset('MaxFunEval',10^5);
%optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','serunc');


z=inst;
x=data(:,2:end);
y=data(:,1);

P_z=z/(z'*z)*z'; %% projection matrix
tsls_beta=(x'*P_z*x)\(x'*P_z*y);
e=y-x*tsls_beta;
n=size(data,1);
varcovar=inv(x'*P_z*x)*x'*e*(inv(x'*P_z*x)*x'*e)'/n;
std=sqrt(diag(varcovar));
e=y-x*tsls_beta;
vc=diag( diag(e*e') );
varcovar=(inv(x'*x)*x')*vc*(inv(x'*x)*x')';
std=sqrt(diag(varcovar));
tstat_tsls=tsls_beta./std;

z=x; x=inst;
cu_gmm; %% output at cuout.txt

startval=coeff;

[theta_final, S_final, J_test, probJ, bandw, var_theta, std_theta, conf_inter] = cugmmest(options, data, @popmomnew, startval, inst);
coeff=theta_final;
tstat=theta_final./std_theta;
[coeff,tstat]
j=2;
res(1,j)=coeff(1);
res(2,j)=tstat(1);
res(3,j)=coeff(2);
res(4,j)=tstat(2);
res(5,j)=coeff(3);
res(6,j)=tstat(3);
res(end-1,j)=J_test;
res(end,j)=probJ;

%% three month

%% first stage
[coeff,t,r,n]=fastnw(threemonthprobleft8_lag(1:end-23,2),[lagrr10(idx:idx1,3),lagbf1(idx:idx1,3),lagbf2(idx:idx1,3),lagatmv(idx:idx1,3)],1,0)
[coeff,t,r,n]=fastnw(threemonthprobright8_lag(1:end-23,2),[lagrr10(idx:idx1,3),lagbf1(idx:idx1,3),lagbf2(idx:idx1,3),lagatmv(idx:idx1,3)],1,0)

n=size(int(idx:idx1),1);

inst=[lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)]; %% instruments
data=[int(idx:idx1)/1000,ones(n,1),threemonthprobleft8_lag(1:end-23,2),threemonthprobright8_lag(1:end-23,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2)]; %% dataset used for the estimation
[coeff,t]=fastnw(data(:,1),[data(:,2:end)],0);
temp=[data,inst];
temp=packr(temp);
data=temp(:,[1:11]);
inst=temp(:,[2 12:end]); %% we need the constant as another instrument
t=size(data,1);
maxlag=floor(4*(t/100)^(2/9));
maxlag=60;
%optset('cugmmest','method','serunc');
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
options=optimset('MaxFunEval',10^5);
optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','serunc');


z=inst;
x=data(:,2:end);
y=data(:,1);
P_z=z/(z'*z)*z'; %% projection matrix
xtag=P_z*x;
tsls_beta=(xtag'*xtag)\(xtag'*y);
e=y-x*tsls_beta;
n=size(data,1);
e=y-x*tsls_beta;
vc=diag( diag(e*e') );
varcovar=(inv(xtag'*xtag)*xtag')*vc*(inv(xtag'*xtag)*xtag')';
std=sqrt(diag(varcovar));
tstat_tsls=tsls_beta./std;


z=x; x=inst;
cu_gmm; %% output at cuout.txt

startval=coeff;

[theta_final, S_final, J_test, probJ, bandw, var_theta, std_theta, conf_inter] = cugmmest(options, data, @popmomnew, startval, inst);
coeff=theta_final;
tstat=theta_final./std_theta;
[coeff,tstat]
j=3;
res(1,j)=coeff(1);
res(2,j)=tstat(1);
res(3,j)=coeff(2);
res(4,j)=tstat(2);
res(5,j)=coeff(3);
res(6,j)=tstat(3);
res(end-1,j)=J_test;
res(end,j)=probJ;

%% 6 month

%% first stage
[coeff,t,r,n]=fastnw(sixmonthprobleft10_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)
[coeff,t,r,n]=fastnw(sixmonthprobright10_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)

n=size(int(idx:idx1),1);

inst=[lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)]; %% instruments
data=[int(idx:idx1)/1000,ones(n,1),sixmonthprobleft10_lag(1:end-23,2),sixmonthprobright10_lag(1:end-23,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2)]; %% dataset used for the estimation

[coeff,t]=fastnw(data(:,1),[data(:,2:end)],0,-1);

temp=[data,inst];
temp=packr(temp);
data=temp(:,[1:11]);
inst=temp(:,[2 12:end]); %% we need the constant as another instrument
t=size(data,1);
maxlag=floor(4*(t/100)^(2/9));
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
options=optimset('MaxFunEval',10^5);
optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','serunc');


z=inst;
x=data(:,2:end);
y=data(:,1);
P_z=z/(z'*z)*z'; %% projection matrix
xtag=P_z*x;
tsls_beta=(xtag'*xtag)\(xtag'*y);
e=y-x*tsls_beta;
n=size(data,1);
e=y-x*tsls_beta;
vc=diag( diag(e*e') );
varcovar=(inv(xtag'*xtag)*xtag')*vc*(inv(xtag'*xtag)*xtag')';
std=sqrt(diag(varcovar));
tstat_tsls=tsls_beta./std;

z=x; x=inst;
cu_gmm; %% output at cuout.txt

startval=coeff;

[theta_final, S_final, J_test, probJ, bandw, var_theta, std_theta, conf_inter] = cugmmest(options, data, @popmomnew, startval, inst);
coeff=theta_final;
tstat=theta_final./std_theta;
[coeff,tstat]

j=4;
res(1,j)=coeff(1);
res(2,j)=tstat(1);
res(3,j)=coeff(2);
res(4,j)=tstat(2);
res(5,j)=coeff(3);
res(6,j)=tstat(3);
res(end-1,j)=J_test;
res(end,j)=probJ;

%% 9 month
%% first stage
[coeff,t,r,n]=fastnw(ninemonthprobleft13_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)
[coeff,t,r,n]=fastnw(ninemonthprobright13_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)


n=size(int(idx:idx1),1);
inst=[lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)]; %% instruments
data=[int(idx:idx1)/1000,ones(n,1),ninemonthprobleft13_lag(1:end-23,2),ninemonthprobright13_lag(1:end-23,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2)]; %% dataset used for the estimation
[coeff,t]=fastnw(data(:,1),[data(:,2:end)],0);
temp=[data,inst];
temp=packr(temp);
data=temp(:,[1:11]);
inst=temp(:,[2 12:end]); %% we need the constant as another instrument
t=size(data,1);
maxlag=floor(4*(t/100)^(2/9));

optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
options=optimset('MaxFunEval',10^5);
optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','serunc');




z=inst;
x=data(:,2:end);
y=data(:,1);

P_z=z/(z'*z)*z'; %% projection matrix
tsls_beta=(x'*P_z*x)\(x'*P_z*y);
e=y-x*tsls_beta;
n=size(data,1);
varcovar=inv(x'*P_z*x)*x'*e*(inv(x'*P_z*x)*x'*e)'/n;
std=sqrt(diag(varcovar));
e=y-x*tsls_beta;
vc=diag( diag(e*e') );
varcovar=(inv(x'*x)*x')*vc*(inv(x'*x)*x')';
std=sqrt(diag(varcovar));
tstat_tsls=tsls_beta./std;

z=x; x=inst;
cu_gmm; %% output at cuout.txt



startval=coeff;

[theta_final, S_final, J_test, probJ, bandw, var_theta, std_theta, conf_inter] = cugmmest(options, data, @popmomnew, startval, inst);
coeff=theta_final;
tstat=theta_final./std_theta;
[coeff,tstat]

j=5;
res(1,j)=coeff(1);
res(2,j)=tstat(1);
res(3,j)=coeff(2);
res(4,j)=tstat(2);
res(5,j)=coeff(3);
res(6,j)=tstat(3);
res(end-1,j)=J_test;
res(end,j)=probJ;

%% 12 month 
[coeff,t,r,n]=fastnw(twelvemonthprobleft20_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)
[coeff,t,r,n]=fastnw(twelvemonthprobright20_lag(1:end-23,2),[lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)],1,0)


n=size(int(idx:idx1),1);
inst=[lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),lagatmv(idx:idx1,2)]; %% instruments
data=[int(idx:idx1)/1000,ones(n,1),twelvemonthprobleft20_lag(1:end-23,2),twelvemonthprobright20_lag(1:end-23,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagint(idx:idx1,2)/1000,intweeklag(idx:idx1,2)]; %% dataset used for the estimation
[coeff,t]=fastnw(data(:,1),[data(:,2:end)],0);
temp=[data,inst];
temp=packr(temp);
data=temp(:,[1:11]);
inst=temp(:,[2 12:end]); %% we need the constant as another instrument
t=size(data,1);
maxlag=floor(4*(t/100)^(2/9));
%maxlag=60;
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
options=optimset('MaxFunEval',10^5);
optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','serunc');

z=inst;
x=data(:,2:end);
y=data(:,1);
P_z=z/(z'*z)*z'; %% projection matrix
xtag=P_z*x;
tsls_beta=(xtag'*xtag)\(xtag'*y);
e=y-x*tsls_beta;
n=size(data,1);
e=y-x*tsls_beta;
vc=diag( diag(e*e') );
varcovar=(inv(xtag'*xtag)*xtag')*vc*(inv(xtag'*xtag)*xtag')';
std=sqrt(diag(varcovar));
tstat_tsls=tsls_beta./std;

z=x; x=inst;
cu_gmm; %% output at cuout.txt



startval=coeff;

[theta_final, S_final, J_test, probJ, bandw, var_theta, std_theta, conf_inter] = cugmmest(options, data, @popmomnew, startval, inst);
coeff=theta_final;
tstat=theta_final./std_theta;
[coeff,tstat]

j=6;
res(1,j)=coeff(1);
res(2,j)=tstat(1);
res(3,j)=coeff(2);
res(4,j)=tstat(2);
res(5,j)=coeff(3);
res(6,j)=tstat(3);
res(end-1,j)=J_test;
res(end,j)=probJ;


%%% Panel B

