format long;
clear;

load ECON_variable_a.mat 
load NONFUND_variable_new.mat 
    

    Mt=aligned_Mt(469:546,1); %combined non-fundamental predictor in whole sample
    mu_t=aligned_Yt(469:546,1); %combined fundamental predictor in whole sample
    
input_file='data20200113.xls';
input_sheet='fears';
fears=xlsread(input_file,input_sheet,'e2:e79');
aaa = sort(fears, 'descend');
Ih2=zeros(78,1);
for i=1:78
    if fears(i)<aaa(39)
        Ih2(i,1)=0;
    else
        Ih2(i,1)=1;
    end
end
th2=find(Ih2==1);
tl2=find(Ih2==0);
MtH=zscore(Mt(th2,1));
MtL=zscore(Mt(tl2,1));
muH=zscore(mu_t(th2,1));
muL=zscore(mu_t(tl2,1));
Rt=Rt(469:end,1);
RtH=Rt(th2,1);
RtL=Rt(tl2,1);
DataSize=78;
Number_High=size(th2,1);

%% whole sample


A1=[ones(DataSize,1) Mt(1:DataSize) mu_t(1:DataSize)];
[X1,STDX1]=lscov(A1,Rt(1:DataSize));
R_squared1=1-sum((Rt(1:DataSize)-X1(1)-X1(2)*Mt(1:DataSize)-X1(3)*mu_t(1:DataSize)).^2)/sum((Rt(1:DataSize)-mean(Rt(1:DataSize))).^2);
nlag=12; 
RegressionCoef1=nwest(Rt(1:DataSize),A1,nlag);
Newey_West_t=RegressionCoef1.tstat;

%% high regime


    NumberH=Number_High;
    A2=[ones(NumberH,1) MtH muH];
    [X2,STDX2]=lscov(A2,RtH);
    R_squared2=1-sum((RtH-X2(1)-X2(2)*MtH-X2(3)*muH).^2)/sum((RtH-mean(RtH)).^2);
    RegressionCoef2=nwest(RtH,A2,nlag);
    Newey_West_tH=RegressionCoef2.tstat;
%% low regime    
    NumberL=DataSize-Number_High;
    A3=[ones(NumberL,1) MtL muL];
    [X3,STDX3]=lscov(A3,RtL);
    R_squared3=1-sum((RtL-X3(1)-X3(2)*MtL-X3(3)*muL).^2)/sum((RtL-mean(RtL)).^2);
    RegressionCoef3=nwest(RtL,A3,nlag);
    Newey_West_tL=RegressionCoef3.tstat;  

X1=X1';
X2=X2';
X3=X3';
TX1=Newey_West_t';
TX2=Newey_West_tH';
TX3=Newey_West_tL';

%% Bootstrap p-values
NW_t=zeros(3,2000);
NW_tH=zeros(3,2000);
NW_tL=zeros(3,2000);
for j=1:2000
    seeds=j;
% full sample
    epsilon_t=Rt(1:DataSize)-X1(1)-X1(2)*Mt(1:DataSize)-X1(3)*mu_t(1:DataSize);
    X_m=lscov([ones(DataSize-1,1) Mt(1:DataSize-1)],Mt(2:DataSize));
    rho1_m=X_m(2)+((1+3*X_m(2))/(DataSize))+(3*(1+3*X_m(2))/(DataSize)^2); 
    rho0_m=mean(Mt(2:DataSize)-rho1_m*Mt(1:DataSize-1));
    phi_m=Mt(2:DataSize)-rho0_m-rho1_m*Mt(1:DataSize-1);
    X_mu=lscov([ones(DataSize-1,1) mu_t(1:DataSize-1)],mu_t(2:DataSize));
    rho1_mu=X_mu(2)+((1+3*X_mu(2))/(DataSize))+(3*(1+3*X_mu(2))/(DataSize)^2); 
    rho0_mu=mean(mu_t(2:DataSize)-rho1_mu*mu_t(1:DataSize-1));
    phi_mu=mu_t(2:DataSize)-rho0_mu-rho1_mu*mu_t(1:DataSize-1);
    randn('state',seeds)
    Wt=normrnd(0,1,DataSize,1); 
    pseudo_R=mean(Rt(1:DataSize))+epsilon_t.*Wt;
    pseudo_m=zeros(DataSize,1);
    pseudo_m(1)=Mt(1);
    pseudo_mu=zeros(DataSize,1);
    pseudo_mu(1)=mu_t(1);
    for i=2:DataSize
        pseudo_m(i)=rho0_m+rho1_m*pseudo_m(i-1)+phi_m(i-1)*Wt(i-1);
        pseudo_mu(i)=rho0_mu+rho1_mu*pseudo_mu(i-1)+phi_mu(i-1)*Wt(i-1);
    end
    RegressionCoef4=nwest(pseudo_R,[ones(DataSize,1) pseudo_m pseudo_mu],nlag);
    NW_t(:,j)=RegressionCoef4.tstat;
% regime H
    epsilon_tH=RtH-X2(1)-X2(2)*MtH-X2(3)*muH;
    X_mH=lscov([ones(NumberH-1,1) MtH(1:NumberH-1)],MtH(2:NumberH));
    rho1_mH=X_mH(2)+((1+3*X_mH(2))/(NumberH))+(3*(1+3*X_mH(2))/(NumberH)^2); 
    rho0_mH=mean(MtH(2:NumberH)-rho1_mH*MtH(1:NumberH-1));
    phi_mH=MtH(2:NumberH)-rho0_mH-rho1_mH*MtH(1:NumberH-1);
    X_muH=lscov([ones(NumberH-1,1) muH(1:NumberH-1)],muH(2:NumberH));
    rho1_muH=X_muH(2)+((1+3*X_muH(2))/(NumberH))+(3*(1+3*X_muH(2))/(NumberH)^2); 
    rho0_muH=mean(muH(2:NumberH)-rho1_muH*muH(1:NumberH-1));
    phi_muH=muH(2:NumberH)-rho0_muH-rho1_muH*muH(1:NumberH-1);
    pseudo_RH=mean(RtH(1:NumberH))+epsilon_tH.*Wt(1:NumberH);
    pseudo_mH=zeros(NumberH,1);
    pseudo_mH(1)=MtH(1);
    pseudo_muH=zeros(NumberH,1);
    pseudo_muH(1)=muH(1);
    for i=2:NumberH
        pseudo_mH(i)=rho0_mH+rho1_mH*pseudo_mH(i-1)+phi_mH(i-1)*Wt(i-1);
        pseudo_muH(i)=rho0_muH+rho1_muH*pseudo_muH(i-1)+phi_muH(i-1)*Wt(i-1);
    end
    RegressionCoef5=nwest(pseudo_RH,[ones(NumberH,1) pseudo_mH pseudo_muH],nlag);
    NW_tH(:,j)=RegressionCoef5.tstat;
% regime L
    epsilon_tL=RtL-X3(1)-X3(2)*MtL-X3(3)*muL;
    X_mL=lscov([ones(NumberL-1,1) MtL(1:NumberL-1)],MtL(2:NumberL));
    rho1_mL=X_mL(2)+((1+3*X_mL(2))/(NumberL))+(3*(1+3*X_mL(2))/(NumberL)^2); 
    rho0_mL=mean(MtL(2:NumberL)-rho1_mL*MtL(1:NumberL-1));
    phi_mL=MtL(2:NumberL)-rho0_mL-rho1_mL*MtL(1:NumberL-1);
    X_muL=lscov([ones(NumberL-1,1) muL(1:NumberL-1)],muL(2:NumberL));
    rho1_muL=X_muL(2)+((1+3*X_muL(2))/(NumberL))+(3*(1+3*X_muL(2))/(NumberL)^2); 
    rho0_muL=mean(muL(2:NumberL)-rho1_muL*muL(1:NumberL-1));
    phi_muL=muL(2:NumberL)-rho0_muL-rho1_muL*muL(1:NumberL-1);
    pseudo_RL=mean(RtL(1:NumberL))+epsilon_tL.*Wt(1:NumberL);
    pseudo_mL=zeros(NumberL,1);
    pseudo_mL(1)=MtL(1);
    pseudo_muL=zeros(NumberL,1);
    pseudo_muL(1)=muL(1);
    for i=2:NumberL
        pseudo_mL(i)=rho0_mL+rho1_mL*pseudo_mL(i-1)+phi_mL(i-1)*Wt(i-1);
        pseudo_muL(i)=rho0_muL+rho1_muL*pseudo_muL(i-1)+phi_muL(i-1)*Wt(i-1);
    end
    RegressionCoef6=nwest(pseudo_RL,[ones(NumberL,1) pseudo_mL pseudo_muL],nlag);
    NW_tL(:,j)=RegressionCoef6.tstat;

end

Switch_sided=1; %1: one-sided bootstrap p-value? 2: two-sided bootstrap p-value
if Switch_sided==1 
	if Newey_West_t(2)>0 
	     p_full_m=mean(NW_t(2,:)>=Newey_West_t(2));% the percentage of NW_t(2,:) is greater than Newey_West_t(2)
	elseif Newey_West_t(2)<0
	     p_full_m=mean(NW_t(2,:)<=Newey_West_t(2));
	elseif Newey_West_t(2)==0
	     p_full_m=1;
	end
	 if Newey_West_t(3)>0 
	     p_full_mu=mean(NW_t(3,:)>=Newey_West_t(3));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	 elseif Newey_West_t(3)<0
	     p_full_mu=mean(NW_t(3,:)<=Newey_West_t(3));
	 elseif Newey_West_t(3)==0
	     p_full_mu=1;
	 end
	 if Newey_West_tH(2)>0 
	     p_H_m=mean(NW_tH(2,:)>=Newey_West_tH(2));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	 elseif Newey_West_tH(2)<0
	     p_H_m=mean(NW_tH(2,:)<=Newey_West_tH(2));
	 elseif Newey_West_tH(2)==0
	     p_H_m=1;
	 end
	 if Newey_West_tH(3)>0 
	     p_H_mu=mean(NW_tH(3,:)>=Newey_West_tH(3));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	 elseif Newey_West_tH(3)<0
	     p_H_mu=mean(NW_tH(3,:)<=Newey_West_tH(3));
	 elseif Newey_West_tH(3)==0
	     p_H_mu=1;
	 end
	 if Newey_West_tL(2)>0 
	     p_L_m=mean(NW_tL(2,:)>=Newey_West_tL(2));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	 elseif Newey_West_tL(2)<0
	     p_L_m=mean(NW_tL(2,:)<=Newey_West_tL(2));
	 elseif Newey_West_tL(2)==0
	     p_L_m=1;
	 end
	 if Newey_West_tL(3)>0 
	     p_L_mu=mean(NW_tL(3,:)>=Newey_West_tL(3));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	 elseif Newey_West_tL(3)<0
	     p_L_mu=mean(NW_tL(3,:)<=Newey_West_tL(3));
	 elseif Newey_West_tL(3)==0
	     p_L_mu=1;
	end
elseif Switch_sided==2
    p_full_m=mean(abs(NW_t(2,:))>=abs(Newey_West_t(2)));
    p_full_mu=mean(abs(NW_t(3,:))>=abs(Newey_West_t(3)));
    p_H_m=mean(abs(NW_tH(2,:))>=abs(Newey_West_tH(2)));
    p_H_mu=mean(abs(NW_tH(3,:))>=abs(Newey_West_tH(3)));
    p_L_m=mean(abs(NW_tL(2,:))>=abs(Newey_West_tL(2)));
    p_L_mu=mean(abs(NW_tL(3,:))>=abs(Newey_West_tL(3)));
end
%% save results
Table1=[X1(3)*100 X2(3)*100 X3(3)*100;
    p_full_mu p_H_mu p_L_mu;
    TX1(3) TX2(3) TX3(3);
    X1(2)*100 X2(2)*100 X3(2)*100;
    p_full_m  p_H_m p_L_m;
    TX1(2) TX2(2) TX3(2);
    R_squared1*100 R_squared2*100 R_squared3*100];  