%% Contemporaenous sub-sample analysis
%% Date: 05/09/2021
%% Intervention effect on usd/ils and NEER using IVs and GMM
%% 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(isnan(int))=0;

%%%% First sub-sample

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

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

%% ccys
load ilsnyccy; load dilsnyccy;
load NEERny;

%% Change in NEER
load lagNEER1
load lagNEER2
load lagNEER5
load lagNEER60
load lagNEER10
load newforwards
load newbasis
%load lagNEER3

%% Flows
load laginstflows
load lagforeignflows


%%%% Exogenous
%% lagged interest rate change
%% 1-day change in interest rate
load uslibor; load telbor;
lead=5;
dtelbor=[repmat(NaN,lead,6);telbor(lead+1:end,:)-telbor(1:end-lead,:)];
dlibor=[repmat(NaN,lead,6);uslibor(lead+1:end,:)-uslibor(1:end-lead,:)];
%% Take the three-month change
lagtelborchange=lagmatrix(dtelbor(:,2),0:p);
lagliborchange=lagmatrix(dlibor(:,2),0:p);
save('lagtelborchange','lagtelborchange');
save('lagliborchange','lagliborchange');

%1-day libor change 
lead=1;
dlibor=[repmat(NaN,lead,6);uslibor(lead+1:end,:)-uslibor(1:end-lead,:)];
dtelbor=[repmat(NaN,lead,6);telbor(lead+1:end,:)-telbor(1:end-lead,:)];
lagtelborchange1=lagmatrix(dtelbor(:,2),0:p);
%% Take the three-month change
lagliborchange1=lagmatrix(dlibor(:,2),0:p);
save('lagliborchange1','lagliborchange1');
save('lagtelborchange1','lagtelborchange1');


load lageurusd;
load intweeklag;
%% 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');
%%% dforward
dforwards=[repmat(NaN,1,5);100*log(newforwards(2:end,:)./newforwards(1:end-1,:))]; %% in %
save('dforwards','dforwards');
load newbasis;

%% Panel A:
%%  Regular OLS %%%
%% there use to be an error here
data=[dilsnyccy,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,2),lagliborchange1(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%
cd(ini.dirs.code)
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');

intlag=lagmatrix(int,0:10);
load lagatmvshort; load lagatmvsimple; load lagrr; load lagbf1; load lagbf2; load lagusdils5;load lagNEER1;load lagNEER5; load lagusdils60;
data=[dilsny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');
options.MaxFunEvals=10000;
%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
%startval=[x1;x2;x3;x4;x5]; 
%startval=[0;0;0;0;0];
startval=[0;x2;x3;x4;x5]; 
[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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];


%%%  Panel B: NEER estimation results
%%  Regular OLS %%%
cd(ini.dirs.code)

data=[NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,2),lagliborchange1(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,1);
x1=coeff(2),x2=coeff(3),x3=coeff(4);x4=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%
cd(ini.dirs.code)
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
data=[NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange1(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,2),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');
options.MaxFunEvals=1000;

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));

optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[x1;x2;x3;x4;x5];
[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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];

%%%% is there a difference effect between NEER and USD/ILS?
cd(ini.dirs.code)

data=[dilsny-NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,5);
x1=coeff(2),x2=coeff(3),x3=coeff(4);x4=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
cd(ini.dirs.code)

options=optimset('fminunc');
data=[dilsny-NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));

optset('cugmmest','bandw',maxlag);

options.MaxFunEvals=1000;
startval=[0;0;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]

%%% Are the two series statistically different? Answer is NO: t-stat of
%%% -0.0465 with point estimate of -.2633

cd(ini.dirs.code);
%% Panel C: 3-month forward:

data=[dforwards(:,2),ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,2),lagliborchange1(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);



%%% GMM %%%
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
intlag=lagmatrix(int,0:10);
load lagatmvshort; load lagatmvsimple; load lagrr; load lagbf1; load lagbf2; load lagusdils5;load lagNEER1;load lagNEER5; load lagusdils60;
data=[dforwards(:,2),ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');

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

%% are forwards and spots really different?

data=[dforwards(:,2)-dilsny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
intlag=lagmatrix(int,0:10);
load lagatmvshort; load lagatmvsimple; load lagrr; load lagbf1; load lagbf2; load lagusdils5;load lagNEER1;load lagNEER5; load lagusdils60;
data=[dforwards(:,2)-dilsny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[0;x2;x3;x4;x5];
[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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];


%%%% Basis
cd(ini.dirs.code);

%% why is newbasis messed up?
%rho=(360/(30*3))*(dforwards(:,2)-dilsny); %% yearly change in forward premium

dnewbasis3m=[NaN;newbasis(2:end,2)-newbasis(1:end-1,2)];
data=[dnewbasis3m,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];

[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%

%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
intlag=lagmatrix(int,0:10);

data=[dnewbasis3m,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7:end]); %% col. 2 is the constant
t=size(data,1);

cd('gmm');
optset('cugmmest','method','hacc_b');

%% 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=[0;0;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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];


%%%%% Second sub-sample %%%

%% dates
idx=find(ccydates==datenum('01/01/2016','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 dilsnyccy;
load NEERny;

%% Change in NEER
load lagNEER1
load lagNEER2
load lagNEER5
load lagNEER60
load lagNEER10
load newforwards
load newbasis
%load lagNEER3

%% Flows
load laginstflows
load lagforeignflows


%%%% Exogenous
%% lagged interest rate change
%% 1-day change in interest rate
load uslibor; load telbor;
lead=5;
dtelbor=[repmat(NaN,lead,6);telbor(lead+1:end,:)-telbor(1:end-lead,:)];
dlibor=[repmat(NaN,lead,6);uslibor(lead+1:end,:)-uslibor(1:end-lead,:)];
%% Take the three-month change
lagtelborchange=lagmatrix(dtelbor(:,2),0:p);
lagliborchange=lagmatrix(dlibor(:,2),0:p);
save('lagtelborchange','lagtelborchange');
save('lagliborchange','lagliborchange');

%1-day libor change 
lead=1;
dlibor=[repmat(NaN,lead,6);uslibor(lead+1:end,:)-uslibor(1:end-lead,:)];
dtelbor=[repmat(NaN,lead,6);telbor(lead+1:end,:)-telbor(1:end-lead,:)];
lagtelborchange1=lagmatrix(dtelbor(:,2),0:p);
%% Take the three-month change
lagliborchange1=lagmatrix(dlibor(:,2),0:p);
save('lagliborchange1','lagliborchange1');
save('lagtelborchange1','lagtelborchange1');


load lageurusd;
load intweeklag;
%% 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');
%%% dforward
dforwards=[repmat(NaN,1,5);100*log(newforwards(2:end,:)./newforwards(1:end-1,:))]; %% in %
save('dforwards','dforwards');
load newbasis;

%% Panel A:
%%  Regular OLS %%%
%% there use to be an error here
data=[dilsnyccy,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,2),lagliborchange1(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%
cd(ini.dirs.code)
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');

intlag=lagmatrix(int,0:10);
load lagatmvshort; load lagatmvsimple; load lagrr; load lagbf1; load lagbf2; load lagusdils5;load lagNEER1;load lagNEER5; load lagusdils60;
data=[dilsny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,2),lagliborchange1(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange1(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,2),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');
options.MaxFunEvals=10000;
%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[x1;x2;x3;x4;x5]; 
%startval=[0;0;0;0;0];
%startval=[0;x2;x3;x4;x5]; 
[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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];


%%%  Panel B: NEER estimation results
%%  Regular OLS %%%
cd(ini.dirs.code)

data=[NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,2),lagliborchange1(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,1);
x1=coeff(2),x2=coeff(3),x3=coeff(4);x4=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%
cd(ini.dirs.code)
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
data=[NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange1(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,2),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');
options.MaxFunEvals=1000;

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));

optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[x1;x2;x3;x4;x5];
[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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];

%%%% is there a difference effect between NEER and USD/ILS?
cd(ini.dirs.code)

data=[dilsny-NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0,5);
x1=coeff(2),x2=coeff(3),x3=coeff(4);x4=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
cd(ini.dirs.code)

options=optimset('fminunc');
data=[dilsny-NEERny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));

optset('cugmmest','bandw',maxlag);

options.MaxFunEvals=1000;
startval=[0;0;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]

%%% Are the two series statistically different? Answer is NO: t-stat of
%%% -0.0465 with point estimate of -.2633

cd(ini.dirs.code);
%% Panel C: 3-month forward:

data=[dforwards(:,2),ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,2),lagliborchange1(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);



%%% GMM %%%
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
intlag=lagmatrix(int,0:10);
load lagatmvshort; load lagatmvsimple; load lagrr; load lagbf1; load lagbf2; load lagusdils5;load lagNEER1;load lagNEER5; load lagusdils60;
data=[dforwards(:,2),ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');

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

%% are forwards and spots really different?

data=[dforwards(:,2)-dilsny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%
%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
intlag=lagmatrix(int,0:10);
load lagatmvshort; load lagatmvsimple; load lagrr; load lagbf1; load lagbf2; load lagusdils5;load lagNEER1;load lagNEER5; load lagusdils60;
data=[dforwards(:,2)-dilsny,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7: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');

%% GMM estimation
maxlag=floor(4*(t/100)^(2/9));
optset('cugmmest','bandw',maxlag); %% still need to understand how to properly adjust this number.
startval=[0;x2;x3;x4;x5];
[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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];


%%%% Basis
cd(ini.dirs.code);

%% why is newbasis messed up?
%rho=(360/(30*3))*(dforwards(:,2)-dilsny); %% yearly change in forward premium

dnewbasis3m=[NaN;newbasis(2:end,2)-newbasis(1:end-1,2)];
data=[dnewbasis3m,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];

[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);

%%% GMM %%%

%% we built the moment condition function in momgrad.
%% usd/ils
options=optimset('fminunc');
intlag=lagmatrix(int,0:10);

data=[dnewbasis3m,ones(T,1),int/1000,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),lagliborchange(:,1)]; %% dataset used for the estimation
inst=[lageurusd(idx:idx1,1),lagliborchange(idx:idx1,1),lagvix(idx:idx1,1)-lagvix(idx:idx1,6),intlag(idx:idx1,2),lagusdils60(idx:idx1,2),lagNEER2(idx:idx1,2),lagNEER10(idx:idx1,3),lagcds(idx:idx1,2),lagvix(idx:idx1,2)-lagvix(idx:idx1,20),lagrr10(idx:idx1,2),lagbf1(idx:idx1,2),lagbf2(idx:idx1,2),intweeklag(idx:idx1,2)];
temp=[data(idx:idx1,:),inst];
temp=packr(temp);
data=temp(:,[1 2 3 4 5 6]);
inst=temp(:,[2 7:end]); %% col. 2 is the constant
t=size(data,1);

cd('gmm');
optset('cugmmest','method','hacc_b');

%% 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=[0;0;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]
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
res=[res;J_test;probJ];
