%% Sub-sample analysis of intervention on NEER and USD ILS
%% Date: 05/09/02021
%% Change in USD ILS. I choose NY as that's what I have data for RR
%% Initialize table
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(int==NaN)=0;

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

%% 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

%% 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
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('01/01/2020','dd/mm/yyyy'));


%% reaction function
[coeff,tstat,r,n]=fastnw(int(idx:idx1)/1000,[intlag(idx:idx1,2)/1000,lagusdils60(idx:idx1,2),lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagatmvshort(idx:idx1,3),lagatmvlong(idx:idx1,3),lagrr(idx:idx1,4),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2)],1,10); %% divide int by 1000.
res=NaN(size(coeff,1)*2+1,1);
for i=1:size(coeff,1)
    res(2*i-1,1)=coeff(i)
    res(2*i,1)=tstat(i);
end
res(end)=r;



%%  Regular OLS %%%
PanelA=NaN(8,2);
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('30/12/2016','dd/mm/yyyy'));


data=[dilsny,ones(T,1),int/1000,lagvix(:,1)-lagvix(:,2)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,5);
x1=coeff(1),x2=coeff(2),x3=coeff(3);

coeffs=NaN(size(coeff,1)*2+1,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
%coeffs(end)=r;
PanelA(1:6,1)=coeffs(1:6);
    
%%% GMM %%%

%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('30/12/2016','dd/mm/yyyy'));


inst=[lagvix(idx:idx1,1)-lagvix(idx:idx1,2),intlag(idx:idx1,2)/1000,lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagatmvshort(idx:idx1,3),lagatmvlong(idx:idx1,3),lagrr(idx:idx1,4),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2)]; %% instruments
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4]);
inst=temp(:,[2 5:end]); %% col. 2 is the constant
t=size(data,1);

cd('gmm');
optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','hacc_b');
%optset('cugmmest','method','serunc');
%% 1st stage regression for Kleinbergen-Paap. Used in Table 3.
[coeff,tstat,r,n,~,~,~,~,v]=fastnw(data(:,3),[inst],0,5); %% we exclude the constant
wald=coeff(3:end)'*inv(v(3:end,3:end))*coeff(3:end);%% test for joint significance without the coefficient
dof=size(coeff(4:end),1); 
pwald=1-chi2cdf(wald,dof);
robustF=wald/size(coeff(3:end),1);

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));
%maxlag=1;
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[x1;x2;x3];
[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]

coeffs=NaN(size(coeff,1)*2,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
PanelA(1:6,2)=coeffs(1:6);
PanelA(7 ,2)=J_test;
PanelA(8 ,2)=probJ;


%%%  NEER estimation results
%%  Regular OLS %%%
PanelB=NaN(8,2);
data=[NEERny,ones(T,1),int/1000,lagvix(:,1)-lagvix(:,2)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,5);
x1=coeff(1),x2=coeff(2),x3=coeff(3);

coeffs=NaN(size(coeff,1)*2+1,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
coeffs(end)=r;
PanelB([1:6],1)=coeffs(1:6);
    
%%% GMM %%%
cd(ini.dirs.code)
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');

inst=[lagvix(idx:idx1,1)-lagvix(idx:idx1,2),intlag(idx:idx1,2)/1000,lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagatmvshort(idx:idx1,3),lagatmvlong(idx:idx1,3),lagrr(idx:idx1,5),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2)]; %% instruments
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4]);
inst=temp(:,[2 5:end]); %% col. 2 is the constant
t=size(data,1);

cd('gmm');
optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','hacc_b');
%optset('cugmmest','method','serunc');
%% 1st stage regression for Kleinbergen-Paap. Used in Table 3.
[coeff,tstat,r,n,~,~,~,~,v]=fastnw(data(:,3),[inst],0,5); %% we exclude the constant
wald=coeff(3:end)'*inv(v(3:end,3:end))*coeff(3:end);%% test for joint significance without the coefficient
dof=size(coeff(4:end),1); 
pwald=1-chi2cdf(wald,dof);
robustF=wald/size(coeff(3:end),1);
%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));
%maxlag=1;
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[x1;x2;x3];
[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]

coeffs=NaN(size(coeff,1)*2,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
PanelB([1:6],2)=coeffs;
PanelB(7 ,2)=J_test;
PanelB(8 ,2)=probJ;


%%%%%%%%%%%%%%%% second sub-sample %%%%%%%%%%%%%%%%%%%
%% Change in USD ILS. I choose NY as that's what I have data for RR
%% Initialize table
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;


%% dates
idx=find(ccydates==datenum('02/01/2017','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/01/2020','dd/mm/yyyy'));

%% 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

%% 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');
idx=find(ccydates==datenum('02/01/2017','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/01/2020','dd/mm/yyyy'));

%% reaction function
[coeff,tstat,r,n]=fastnw(int(idx:idx1)/1000,[intlag(idx:idx1,2)/1000,lagusdils60(idx:idx1,2),lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagatmvshort(idx:idx1,3),lagatmvlong(idx:idx1,3),lagrr(idx:idx1,4),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2)],1,10); %% divide int by 1000.
res=NaN(size(coeff,1)*2+1,1);
for i=1:size(coeff,1)
    res(2*i-1)=coeff(i)
    res(2*i)=tstat(i);
end
res(end)=r;

%%  Regular OLS %%%
PanelA=NaN(9,2);
%data=[NEERny,ones(T,1),int/1000,lagvix(:,1)-lagvix(:,2),lagliborchange(:,1)]; %% dataset used for the estimation
idx=find(ccydates==datenum('02/01/2017','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/01/2020','dd/mm/yyyy'));

data=[dilsny,ones(T,1),int/1000,lagvix(:,1)-lagvix(:,2)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,5);
x1=coeff(1),x2=coeff(2),x3=coeff(3);

coeffs=NaN(size(coeff,1)*2+1,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
coeffs(end)=r;
PanelA([1:6],1)=coeffs(1:6);
    

%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
idx=find(ccydates==datenum('02/01/2017','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/01/2020','dd/mm/yyyy'));

inst=[lagvix(idx:idx1,1)-lagvix(idx:idx1,2),intlag(idx:idx1,2)/1000,lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagusdils60(idx:idx1,2),lagcds(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagatmvshort(idx:idx1,3),lagatmvlong(idx:idx1,3),lagrr(idx:idx1,3),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2)]; %% instruments

temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4]);
inst=temp(:,[2 5:end]); %% col. 2 is the constant
t=size(data,1);

cd('gmm');
%optset('cugmmest','method','hacc_b');
%optset('cugmmest','method','hacc_p');
optset('cugmmest','method','serunc');
%% 1st stage regression for Kleinbergen-Paap. Used in Table 3.
[coeff,tstat,r,n,~,~,~,~,v]=fastnw(data(:,3),[inst],0,5); %% we exclude the constant
wald=coeff(3:end)'*inv(v(3:end,3:end))*coeff(3:end);%% test for joint significance without the coefficient
dof=size(coeff(3:end),1); 
pwald=1-chi2cdf(wald,dof);
robustF=wald/size(coeff(3:end),1);

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));
%maxlag=60;
%maxlag=10;
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[0;1;0];
[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]

coeffs=NaN(size(coeff,1)*2,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
PanelA([1:6],2)=coeffs;
PanelA(7 ,2)=J_test;
PanelA(8 ,2)=probJ;

%%%  NEER estimation results
PanelB=NaN(11,2);
idx=find(ccydates==datenum('02/01/2017','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/01/2020','dd/mm/yyyy'));


data=[NEERny,ones(T,1),int/1000,lagvix(:,1)-lagvix(:,2)]; %% dataset used for the estimation
%data=[dilsny,ones(T,1),int/1000,lagvix(:,1)-lagvix(:,2),lagliborchange(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,5);
x0=coeff(1);x1=coeff(2),x2=coeff(3)

coeffs=NaN(size(coeff,1)*2+1,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
coeffs(end)=r;
PanelB([1:6],1)=coeffs(1:6);
    
%%% GMM %%%
cd(ini.dirs.code)
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');

inst=[lagvix(idx:idx1,1)-lagvix(idx:idx1,2),intlag(idx:idx1,2)/1000,lagNEER1(idx:idx1,2),lagNEER5(idx:idx1,3),lagcds(idx:idx1,2),lagusdils60(idx:idx1,2),lagforeignflows(idx:idx1,2)+laginstflows(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagatmvshort(idx:idx1,3),lagatmvlong(idx:idx1,3),lagrr(idx:idx1,4),lagbf1(idx:idx1,3),lagbf2(idx:idx1,3)]; %% instruments
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4]);
inst=temp(:,[2 5:end]); %% col. 2 is the constant
t=size(data,1);

%cd('gmm');
optset('cugmmest','method','hacc_p');
%optset('cugmmest','method','hacc_b');
%optset('cugmmest','method','serunc');
%% 1st stage regression for Kleinbergen-Paap. Used in Table 3.
[coeff,tstat,r,n,~,~,~,~,v]=fastnw(data(:,3),[inst],0,5); %% we exclude the constant
wald=coeff(3:end)'*inv(v(3:end,3:end))*coeff(3:end);%% test for joint significance without the coefficient
dof=size(coeff(3:end),1); 
pwald=1-chi2cdf(wald,dof);
robustF=wald/size(coeff(3:end),1);

%% GMM estimation
%maxlag=floor(4*(t/100)^(2/9));
maxlag=1;
%maxlag=5;
%maxlag=3;
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[0;0;0];
[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]

coeffs=NaN(size(coeff,1)*2,1);
for i=1:size(coeff,1)
    coeffs(2*i-1)=coeff(i);
    coeffs(2*i)=tstat(i);
end
PanelB([1:6],2)=coeffs;
PanelB(7 ,2)=J_test;
PanelB(8 ,2)=probJ;

