%% Date: 26/12/2021
%% Predictive regressinos using GMM+OLS
%% 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;

%% 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');
%% CDS
load('lagcds')

%% ccys
load ilsnyccy; load dilsnyccy;
dilsny=dilsnyccy;
%load ilslondonccy;
load NEERny;

%% Change in NEER
load lagNEER1
load lagNEER2
load lagNEER3
load lagNEER5
load lagNEER10
load lagNEER60
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;
p=100;
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;
load lagint;
%% 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
dforward=[repmat(NaN,1,5);100*log(newforwards(2:end,:)./newforwards(1:end-1,:))]; %% in %
save('dforward','dforward');
load newbasis;
ilsnyccy=ilsnyccy';

results=NaN(8,3); %% t=2, t=3, t=6, t=11;

%% usd
d1usdils=lagmatrix(dilsny,0:p);

[coeff,t,r,n,resid1,pred]=fastnw(d1usdils(idx:idx1,1),lagint(idx:idx1,2),1,0); %% u
[coeff,t,r,n,resid2,pred]=fastnw(lagint(idx:idx1,1),lagint(idx:idx1,2),1,0); %% v
rho_usd=coeff(2);
sigmauv_usd=nancov(resid1,resid2);
sigmauv_usd=sigmauv_usd(1,1);
sigmav_usd=nanvar(resid2);
%% neer
d1NEER=lagmatrix(NEERny,0:p);

[coeff,t,r,n,resid1,pred]=fastnw(d1NEER(idx:idx1,1),lagint(idx:idx1,2),1,0); %% u
[coeff,t,r,n,resid2,pred]=fastnw(lagint(idx:idx1,1),lagint(idx:idx1,2),1,0); %% v
rho_neer=coeff(2);
sigmauv_neer=nancov(resid1,resid2);
sigmauv_neer=sigmauv_neer(1,1);
sigmav_neer=nanvar(resid2);
%% forward
[coeff,t,r,n,resid1,pred]=fastnw(dforward(idx:idx1,1),lagint(idx:idx1,2),1,0); %% u
[coeff,t,r,n,resid2,pred]=fastnw(lagint(idx:idx1,1),lagint(idx:idx1,2),1,0); %% v
rho_frwrd=coeff(2);
sigmauv_frwrd=nancov(resid1,resid2);
sigmauv_frwrd=sigmauv_frwrd(1,1);
sigmav_frwrd=nanvar(resid2);
T=idx1-idx;

%% Panel A:

t=2;
%%  Regular OLS %%%
lead=t;
dilsny=[repmat(NaN,lead,1);log(ilsnyccy(lead+1:end,:)./ilsnyccy(1:end-lead,:) )]*100; %% change in pct.
[coeff,tstat,r,n]=fastnw(dilsny(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,0);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
[coeff(2),tstat(2)/sqrt(lead)]
results(1,1)=coeff(2)+(1/T)*(lead*(1+rho_usd)+2*rho_usd*( (1-rho_usd^lead)/(1-rho_usd) ) )*sigmauv_usd/sigmav_usd;
results(2,1)=tstat(2)/sqrt(lead);

t=3;
%%  Regular OLS %%%

lead=t;
dilsny=[repmat(NaN,lead,1);log(ilsnyccy(lead+1:end,:)./ilsnyccy(1:end-lead,:) )]*100; %% change in pct.
[coeff,tstat,r,n]=fastnw(dilsny(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,-1);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
[coeff(2),tstat(2)/sqrt(lead)]
results(3,1)=coeff(2)+(1/T)*(lead*(1+rho_usd)+2*rho_usd*( (1-rho_usd^lead)/(1-rho_usd) ) )*sigmauv_usd/sigmav_usd;;
results(4,1)=tstat(2)/sqrt(lead);



t=6;
%%  Regular OLS %%%

lead=t;
dilsny=[repmat(NaN,lead,1);log(ilsnyccy(lead+1:end,:)./ilsnyccy(1:end-lead,:) )]*100; %% change in pct.
[coeff,tstat,r,n]=fastnw(dilsny(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,0);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
[coeff(2),tstat(2)/sqrt(lead)]
results(5,1)=coeff(2)+(1/T)*(lead*(1+rho_usd)+2*rho_usd*( (1-rho_usd^lead)/(1-rho_usd) ) )*sigmauv_usd/sigmav_usd;;
results(6,1)=tstat(2)/sqrt(lead);

t=11;
%%  Regular OLS %%%

lead=t;
dilsny=[repmat(NaN,lead,1);log(ilsnyccy(lead+1:end,:)./ilsnyccy(1:end-lead,:) )]*100; %% change in pct.
[coeff,tstat,r,n]=fastnw(dilsny(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,0);
res=[coeff,tstat];
res=res'; [m,n]=size(res);
res=reshape(res,m*n,1);
[coeff(2),tstat(2)/sqrt(lead)]
results(7,1)=coeff(2)+(1/T)*(lead*(1+rho_usd)+2*rho_usd*( (1-rho_usd^lead)/(1-rho_usd) ) )*sigmauv_usd/sigmav_usd;;
results(8,1)=tstat(2)/sqrt(lead);


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

t=2;
lead=t;
[coeff,tstat,r,n]=fastnw(lagNEER2(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,0);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(1,2)=coeff(2)+(1/T)*(lead*(1+rho_neer)+2*rho_usd*( (1-rho_neer^lead)/(1-rho_neer) ) )*sigmauv_neer/sigmav_neer;
results(2,2)=tstat(2)/sqrt(lead);


%%%%%%% 

t=3;
lead=t;
[coeff,tstat,r,n]=fastnw(lagNEER3(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,0);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(3,2)=coeff(2)+(1/T)*(lead*(1+rho_neer)+2*rho_usd*( (1-rho_neer^lead)/(1-rho_neer) ) )*sigmauv_neer/sigmav_neer;
results(4,2)=tstat(2)/sqrt(lead);



t=5;
lead=t;
[coeff,tstat,r,n]=fastnw(lagNEER5(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,0);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(5,2)=coeff(2)+(1/T)*(lead*(1+rho_neer)+2*rho_usd*( (1-rho_neer^lead)/(1-rho_neer) ) )*sigmauv_neer/sigmav_neer;
results(6,2)=tstat(2)/sqrt(lead);



t=10;
lead=t;
[coeff,tstat,r,n]=fastnw(lagNEER10(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,0);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(7,2)=coeff(2)+(1/T)*(lead*(1+rho_neer)+2*rho_usd*( (1-rho_neer^lead)/(1-rho_neer) ) )*sigmauv_neer/sigmav_neer;
results(8,2)=tstat(2)/sqrt(lead);

%% Panel C: 3-month forward:

t=2;
lead=t;
dforward=[repmat(NaN,lead,1);100*log(newforwards(1+lead:end,2)./newforwards(1:end-lead,2))]; %% in %
[coeff,tstat,r,n]=fastnw(dforward(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,-1);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(1,3)=coeff(2)+(1/T)*(lead*(1+rho_frwrd)+2*rho_usd*( (1-rho_frwrd^lead)/(1-rho_frwrd) ) )*sigmauv_frwrd/sigmav_frwrd;
results(2,3)=tstat(2)/sqrt(lead);


t=3;
lead=t;
dforward=[repmat(NaN,lead,1);100*log(newforwards(1+lead:end,2)./newforwards(1:end-lead,2))]; %% in %
[coeff,tstat,r,n]=fastnw(dforward(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,-1);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(3,3)=coeff(2)+(1/T)*(lead*(1+rho_frwrd)+2*rho_usd*( (1-rho_frwrd^lead)/(1-rho_frwrd) ) )*sigmauv_frwrd/sigmav_frwrd;
results(4,3)=tstat(2)/sqrt(lead);


t=5;
lead=t;
dforward=[repmat(NaN,lead,1);100*log(newforwards(1+lead:end,2)./newforwards(1:end-lead,2))]; %% in %
[coeff,tstat,r,n]=fastnw(dforward(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,-1);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(5,3)=coeff(2)+(1/T)*(lead*(1+rho_frwrd)+2*rho_usd*( (1-rho_frwrd^lead)/(1-rho_frwrd) ) )*sigmauv_frwrd/sigmav_frwrd;
results(6,3)=tstat(2)/sqrt(lead);


t=10;
lead=t;
dforward=[repmat(NaN,lead,1);100*log(newforwards(1+lead:end,2)./newforwards(1:end-lead,2))]; %% in %
[coeff,tstat,r,n]=fastnw(dforward(idx:idx1,1),[lagint(idx:idx1,1+lead-1)/1000,lageurusd(idx:idx1,1+lead-1),lagvix(idx:idx1,1+lead-1)-lagvix(idx:idx1,6+lead-1),lagliborchange(idx:idx1,1+lead-1)],1,-1);
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);
[coeff(2),tstat(2)/sqrt(lead)]
results(7,3)=coeff(2)+(1/T)*(lead*(1+rho_frwrd)+2*rho_usd*( (1-rho_frwrd^lead)/(1-rho_frwrd) ) )*sigmauv_frwrd/sigmav_frwrd;
results(8,3)=tstat(2)/sqrt(lead);


e=y-x*tsls_beta;
n=size(data,1);
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;

res=[tsls_beta,tstat_tsls];

%% GMM
%% for some reason in Bruce Hansen's program the z and the x are confused?
z=x; x=inst;
cu_gmm; %% output at cuout.txt
