%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WARGING: Because of the double recursive estimation and the use of Bootstrap p-values 
%          this code may take some time to conclude the computations.
%
%          Make sure they files are in your working directory of MATLAB.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clc;
clear all;
close all;
warning('off')
cv01=chi2inv(0.99,1); 
cv05=chi2inv(0.95,1); 
cv10=chi2inv(0.90,1);
%% Change if necessary
tau1=0.25;          % As in BLS -  recursive                      
tau2=1/4;           % As in BLS - rolling
tau3=1/4;           % Double recursive

DREC=0;             % if 1 compute double recursive; if 0 do not compute double recursive
B=999;              % Bootstrap Replications 

%% DATA
input_file='Predictive_Variables';
input_sheet='Predictive Reg variables';
Data=xlsread(input_file,input_sheet,'b276:q1092');  

%% 
T=length(Data);
T1=round(tau1*T);
T2=round(tau2*T);

rng('default');

Tests=zeros(14,10);
P_VAL=zeros(14,10);

for kk=1:2

    y0=Data(:,1);
    x0=Data(:,kk+1);

%% IVX Instrument    
 R_nz=1-1/((T)^(0.95));
 d_X=[0; (x0(2:end)-x0(1:end-1))];
 Z_tilde1=zeros(T,1); % instrument matrix

 for t=2:length(x0)
     Z_tilde1(t,:)=R_nz*Z_tilde1(t-1,:)+d_X(t,1);
 end
 
%% Sine Instrument  
 Z_tilde2=sin(pi*(0:T-1)'/(2*(T)));

%% Full Sample
    [t_beta_comb_rob, t_beta_comb_homosc]=Compute_2SLS_Comb_New_demeaned(y0,x0, Z_tilde1, Z_tilde2,1);    
        
%% Recursive Estimation
    [Comb_het, Comb_hom, Comb_het_r, Comb_hom_r, Rec_f, Rec_r]=Recursive_for_back(y0, x0, T1, Z_tilde1, Z_tilde2);
       
%% Rolling Estimation
     [Roll_het, Roll_hom, Beta_est, Var_est,results_Comb_a, results_Comb_b]=rolling_test(y0,x0,T2, Z_tilde1, Z_tilde2);
    
%% Double Recursive
if DREC==1
 [Double_test_Comba, Double_test_Combb]=double_rec(y0,x0,Z_tilde1, Z_tilde2, tau3);
else
    Double_test_Comba=0;
    Double_test_Combb=0;
end

%% Bootstrap   
% Full
Boot_full_rob=zeros(B,1);
Boot_full=zeros(B,1);

% Recursive
Boot_res_Comba_f=zeros(B,1); 
Boot_res_Comba_r=zeros(B,1);

Boot_res_Combb_f=zeros(B,1);
Boot_res_Combb_r=zeros(B,1);

% Rolling
results_Comba_d_B=zeros(B,1);
results_Combb_d_B=zeros(B,1);

% Double_rec Bootstrap
Boot_full_DR_rob=zeros(B,1);
Boot_full_DR=zeros(B,1);

y_star=Fixed_Reg_WildBoot(y0,B);

parfor boot=1:B
%% Full Estimation
[t_beta_comb_rob_B, t_beta_comb_homosc_B]=Compute_2SLS_Comb_New_demeaned(y_star(:,boot), x0,Z_tilde1, Z_tilde2,1);    

Boot_full_rob(boot,1)=t_beta_comb_rob_B>t_beta_comb_rob;
Boot_full(boot,1)=t_beta_comb_homosc_B>t_beta_comb_homosc;
     
%% Recursive Estimation
[Comb_het_B, Comb_hom_B, Comb_het_r_B, Comb_hom_r_B, Rec_f_B, Rec_r_B]=Recursive_for_back(y_star(:,boot), x0, T1,Z_tilde1, Z_tilde2);

    Boot_res_Comba_f(boot,1)=Comb_het_B>Comb_het;
    Boot_res_Comba_r(boot,1)=Comb_het_r_B>Comb_het_r;

    Boot_res_Combb_f(boot,1)=Comb_hom_B>Comb_hom;
    Boot_res_Combb_r(boot,1)=Comb_hom_r_B>Comb_hom_r;

    
%% Rolling Estimation
[Roll_het_B, Roll_hom_B, Beta_est_B, Var_est_B, results_Comb_a_B, results_Comb_b_B]=rolling_test(y_star(:,boot),x0,T2,Z_tilde1, Z_tilde2);   
 
   results_Comba_d_B(boot,1)=Roll_het_B>Roll_het;
   results_Combb_d_B(boot,1)=Roll_hom_B>Roll_hom; 
    
% %% Double Recursive
if DREC==1
  [t_beta_DR_rob_B, t_beta_DR_homosc_B]=double_rec(y_star(:,boot), x0, Z_tilde1, Z_tilde2, tau3);    
      Boot_full_DR_rob(boot,1)=t_beta_DR_rob_B>Double_test_Comba;
      Boot_full_DR(boot,1)=t_beta_DR_homosc_B>Double_test_Combb;
else
    Boot_full_DR_rob(boot,1)=0;
    Boot_full_DR(boot,1)=0;
end
    
end
%% Tests OUTPUT and p-values
Tests(kk,:)=[t_beta_comb_rob, t_beta_comb_homosc, Comb_het, Comb_hom, Comb_het_r, Comb_hom_r, Roll_het, Roll_hom, Double_test_Comba, Double_test_Combb];
P_VAL(kk,:)=[mean(Boot_full_rob), mean(Boot_full) mean(Boot_res_Comba_f), mean(Boot_res_Combb_f), mean(Boot_res_Comba_r), mean(Boot_res_Combb_r), mean(results_Comba_d_B), mean(results_Combb_d_B), mean(Boot_full_DR_rob), mean(Boot_full_DR)];


kk


end
   

[Tests(1,:);P_VAL(1,:);Tests(2,:);P_VAL(2,:);Tests(3,:);P_VAL(3,:);Tests(4,:);P_VAL(4,:);
 Tests(5,:);P_VAL(5,:);Tests(6,:);P_VAL(6,:);Tests(7,:);P_VAL(7,:);Tests(8,:);P_VAL(8,:);
 Tests(9,:);P_VAL(9,:);Tests(10,:);P_VAL(10,:);Tests(11,:);P_VAL(11,:);Tests(12,:);P_VAL(12,:);
 Tests(13,:);P_VAL(13,:);Tests(14,:);P_VAL(14,:)]

%% Functions used for estimation
function results=ols(y,x)
if (nargin ~= 2); error('Wrong # of arguments to ols'); 
else
 [nobs nvar] = size(x); [nobs2 junk] = size(y);
 if (nobs ~= nobs2); error('x and y must have same # obs in ols'); 
 end
end

results.meth = 'ols';
results.y = y;
results.nobs = nobs;
results.nvar = nvar;

if nobs < 10000
  [q r] = qr(x,0);
  xpxi = (r'*r)\eye(nvar);
else % use Cholesky for very large problems
  xpxi = (x'*x)\eye(nvar);
end

results.beta = xpxi*(x'*y);
results.yhat = x*results.beta;
results.resid = y - results.yhat;
end
function [t_beta_robust, t_beta_homoscedastic]=Compute_2SLS_Comb_New_demeaned(y,X,Z_tilde1,Z_tilde2, CC)
%Residuals
epsilon_hat=y(2:end)-[ones(length(y)-1,1) X(1:end-1,:)]*([ones(length(y)-1,1) X(1:end-1,:)]\y(2:end));
U_hat=X(2:end,1)-[ones(length(y)-1,1) X(1:end-1,1)]*([ones(length(y)-1,1) X(1:end-1,1)]\X(2:end,1));

% Compute variances and covariances;

Sigma_hat_ee=(1/length(epsilon_hat))*(epsilon_hat'*epsilon_hat); % Predictive Regression variance

Sigma_hat_uu=(1/length(epsilon_hat))*(U_hat'*U_hat);             % Variance predictor AR 
Omega_hat_uu=Sigma_hat_uu;
Sigma_hat_eu=(1/length(epsilon_hat))*(epsilon_hat'*U_hat);       % Covariance predictive and predictor errors
Omega_hat_eu=Sigma_hat_eu;

% OLS estimation of predictive regression system
n_K=length(y)-1;
X_1=X(1:end-1,:);
y_hat=y(2:end,:)-sum(y(2:end,:))/length(y(2:end,:));

%% Demeaning
Reg_X=X_1-sum(X_1)/length(X_1);
Reg_Z1=Z_tilde1(1:end-1); 
Reg_Z2=Z_tilde2(1:end-1)-sum(Z_tilde2(1:end-1))/length(Z_tilde2);

%%
  INST=[Reg_Z1 Reg_Z2];
  INST1=[Reg_Z1.*y_hat Reg_Z2.*y_hat];
  INST_M=[mean(Reg_Z1) 0];


%% 2SLS
AT=(Reg_X'*INST);
BT=(INST'*INST);
CT=(INST'*y_hat);
DT=(INST1'*INST1);

Sigma_hat_yy2=y_hat'*y_hat/n_K;    

Omega_hat_FM=Sigma_hat_yy2-Omega_hat_eu*inv(Omega_hat_uu)*Omega_hat_eu';
if CC == 1
CORRECT=diag([n_K, 1])*(INST_M'*INST_M)*Omega_hat_FM;
else
    CORRECT=0;
end

t_beta_rob=(AT*inv(BT)*CT)/sqrt(AT*inv(BT)*(DT-CORRECT)*inv(BT)*AT');
t_beta_robust=t_beta_rob*t_beta_rob;
 
t_beta_homosc=(AT*inv(BT)*CT)/sqrt(AT*inv(BT)*(Sigma_hat_yy2*BT-CORRECT)*inv(BT)*AT');
t_beta_homoscedastic=t_beta_homosc*t_beta_homosc;

end 
function [t_beta_robust, t_beta_homoscedastic, beta_roll, Var_beta]=Compute_2SLS_Comb_New_demeaned_roll(y,X,Z_tilde1,Z_tilde2)

%Residuals
epsilon_hat=y(2:end)-[ones(length(y)-1,1) X(1:end-1,:)]*([ones(length(y)-1,1) X(1:end-1,:)]\y(2:end));
U_hat=X(2:end,1)-[ones(length(y)-1,1) X(1:end-1,1)]*([ones(length(y)-1,1) X(1:end-1,1)]\X(2:end,1));

% Compute variances and covariances;

Sigma_hat_ee=(1/length(epsilon_hat))*(epsilon_hat'*epsilon_hat); % Predictive Regression variance

Sigma_hat_uu=(1/length(epsilon_hat))*(U_hat'*U_hat);             % Variance predictor AR 
Omega_hat_uu=Sigma_hat_uu;
Sigma_hat_eu=(1/length(epsilon_hat))*(epsilon_hat'*U_hat);       % Covariance predictive and predictor errors
Omega_hat_eu=Sigma_hat_eu;

% OLS estimation of predictive regression system
n_K=length(y)-1;
X_1=X(1:end-1,:);
y_hat=y(2:end,:)-sum(y(2:end,:))/length(y(2:end,:));

%% Demeaning
Reg_X=X_1-sum(X_1)/length(X_1);
Reg_Z1=Z_tilde1(1:end-1); 
Reg_Z2=Z_tilde2(1:end-1)-sum(Z_tilde2(1:end-1))/length(Z_tilde2);

%%
  INST=[Reg_Z1 Reg_Z2];
  INST1=[Reg_Z1.*y_hat Reg_Z2.*y_hat];
  INST_M=[mean(Reg_Z1) 0];


%% 2SLS
AT=(Reg_X'*INST);
BT=(INST'*INST);
CT=(INST'*y_hat);
DT=(INST1'*INST1);

Sigma_hat_yy2=y_hat'*y_hat/n_K;    
Omega_hat_FM=Sigma_hat_yy2-Omega_hat_eu*inv(Omega_hat_uu)*Omega_hat_eu';

CORRECT=diag([n_K, 1])*(INST_M'*INST_M)*Omega_hat_FM;

beta_roll=(AT*inv(BT)*CT)/(AT*inv(BT)*AT');
Var_beta=inv(AT*inv(BT)*AT')*(AT*inv(BT)*(Sigma_hat_yy2*BT-CORRECT)*inv(BT)*AT')*inv(AT*inv(BT)*AT');

 t_beta_rob=(AT*inv(BT)*CT)/sqrt(AT*inv(BT)*(DT-CORRECT)*inv(BT)*AT');
 t_beta_robust=t_beta_rob*t_beta_rob;
 
 t_beta_homosc=(AT*inv(BT)*CT)/sqrt(AT*inv(BT)*(Sigma_hat_yy2*BT-CORRECT)*inv(BT)*AT');
 t_beta_homoscedastic=t_beta_homosc*t_beta_homosc;
end
function y_star=Fixed_Reg_WildBoot(y,B)
T=size(y,1);
y_star=zeros(T,B);
y_tilde=y-ones(length(y),1)'*y/length(y);
w=randn(T,B); % Normal
    for b=1:B
        y_star(:,b)=w(:,b).*y_tilde;
    end
end
function [Comb_het, Comb_hom, Comb_het_r, Comb_hom_r, Rec_f, Rec_r, Beta_est_f, Var_est_f]=Recursive_for_back(y0, x0, T1,Z_tilde1,Z_tilde2, tau1)
T=length(y0);
y0r=flip(y0);  X0r=flip(x0);   
Z_tilde1r=flip(Z_tilde1); % Reversing the Sample
Z_tilde2r=flip(Z_tilde2);

Comba=zeros(T-T1+1,1); Combb=zeros(T-T1+1,1);
Comba_r=zeros(T-T1+1,1); Combb_r=zeros(T-T1+1,1);
Beta_est_f=zeros(T-T1+1,1);
Var_est_f=zeros(T-T1+1,1);

%% Recursive Estimation
for TL=T1+1:T
y=y0(1:TL,1);
X=x0(1:TL,1);
Z1=Z_tilde1(1:TL,1);
Z2=Z_tilde2(1:TL,1);

yr=flip(y0r(1:TL,1));  
Xr=flip(X0r(1:TL,1));
Z1r=flip(Z_tilde1r(1:TL,1));
Z2r=flip(Z_tilde2r(1:TL,1));

%[t_beta_comb_rob, t_beta_comb_homosc]=Compute_2SLS_Comb_New_demeaned(y,X,Z1,Z2);
[t_beta_comb_rob, t_beta_comb_homosc, beta_rec, Var_beta_rec]=Compute_2SLS_Comb_New_demeaned_roll(y,X,Z1,Z2);
Comba(TL-T1,1)=t_beta_comb_rob; 
Combb(TL-T1,1)=t_beta_comb_homosc;

Beta_est_f(TL-T1,1)=beta_rec;
Var_est_f(TL-T1,1)=Var_beta_rec;


[t_beta_comb_rob_r, t_beta_comb_homosc_r]=Compute_2SLS_Comb_New_demeaned(yr,Xr,Z1r,Z2r,1);
Comba_r(TL-T1,1)=t_beta_comb_rob_r; 
Combb_r(TL-T1,1)=t_beta_comb_homosc_r;
end
Rec_f=Comba;
Rec_r=Combb_r;
Comb_het=max(Comba); 
Comb_hom=max(Combb); 
Comb_het_r=max(Comba_r); 
Comb_hom_r=max(Combb_r);
end
function [Test_het, Test_hom,Beta_est,Var_est, results_Comb_a, results_Comb_b]=rolling_test(y0,x0,T1,Z_tilde1, Z_tilde2)

T=length(y0);
results_Comb_a=zeros(T-T1+1,1);
results_Comb_b=zeros(T-T1+1,1);
Beta_est=zeros(T-T1+1,1);
Var_est=zeros(T-T1+1,1);

   for vvv=0:1:T-T1
        y_roll=y0(1+vvv:T1+vvv,:);
        x_roll=x0(1+vvv:T1+vvv,:);
        Z1_roll=Z_tilde1(1+vvv:T1+vvv,:);
        Z2_roll=Z_tilde2(1+vvv:T1+vvv,:);
        [t_beta_comb_rob, t_beta_comb_homosc, beta_roll, Var_beta]=Compute_2SLS_Comb_New_demeaned_roll(y_roll,x_roll,Z1_roll,Z2_roll);
        %[t_beta_comb_rob, t_beta_comb_homosc]=Compute_2SLS_Comb_New_demeaned(y_roll,x_roll,Z1_roll,Z2_roll);
        results_Comb_a(vvv+1)=t_beta_comb_rob;
        results_Comb_b(vvv+1)=t_beta_comb_homosc;
        Beta_est(vvv+1)=beta_roll;
        Var_est(vvv+1)=Var_beta;
    end
 Test_het=max(results_Comb_a);
 Test_hom=max(results_Comb_b);
 
end
function [Double_test_Comba, Double_test_Combb]=double_rec(y0,x0,Z_tilde1, Z_tilde2, tau3)

T=length(y0);
T1=round(tau3*T);
dim=T-T1+1;
results_Comba_d=zeros(1,dim);
results_Combb_d=zeros(1,dim);

for r2=T1:1:T
    dim0=r2-T1+1;
    results_Comb_a=zeros(dim0,1);
    results_Comb_b=zeros(dim0,1);
    
    for r1=1:1:dim0
        y_dbl=y0(r1:r2);
        x_dbl=x0(r1:r2);
        Z1_dbl=Z_tilde1(r1:r2);
        Z2_dbl=Z_tilde2(r1:r2);
 if tau3==1/3 
     CC=0; 
 else
     CC=1;
 end       
        [t_beta_comb_rob, t_beta_comb_homosc]=Compute_2SLS_Comb_New_demeaned(y_dbl,x_dbl,Z1_dbl,Z2_dbl,CC);
        results_Comb_a(r1)=t_beta_comb_rob;
        results_Comb_b(r1)=t_beta_comb_homosc;
    end
    results_Comba_d(1,r2)=max(results_Comb_a);
    results_Combb_d(1,r2)=max(results_Comb_b);
   
end

Double_test_Comba=max(results_Comba_d(1,:),[],2);
Double_test_Combb=max(results_Combb_d(1,:),[],2);

end