format long;
clear;
load discount_rate_CRSP.mat
load NONFUND_variable_new.mat 
load ECON_variable_a.mat 


%% regimes H and L
load Para_sentiment.mat
    STATES=(sign(pi1(2:end)-0.5)+1)/2; %STATES=1 if at state 1;  STATES=0 if at state 2
    Number_High=high_m;
DP_tMINUS1H=zeros(Number_High,1); % dividend price ratio as independent variable
DPH=zeros(Number_High,1);      % dividend price ratio as dependent variable
DGH=zeros(Number_High,1);
DP_tMINUS1L=zeros(DataSize-Number_High,1);
DPL=zeros(DataSize-Number_High,1);
DGL=zeros(DataSize-Number_High,1);
NumberH=0;
NumberL=0;
for i=1:DataSize
    if STATES(i)==1
        NumberH=NumberH+1;
        DP_tMINUS1H(NumberH)=DP(i);
        DPH(NumberH)=DP(i+1);
        DGH(NumberH)=DG(i);
    elseif STATES(i)==0
        NumberL=NumberL+1;
        DP_tMINUS1L(NumberL)=DP(i);
        DPL(NumberL)=DP(i+1);
        DGL(NumberL)=DG(i);
    end
end
%% regressions
Mt=aligned_Mt;
mu_t=aligned_Yt; 

Switch_dependent=1; %dependent variable as 1: DP, 2: DG
if Switch_dependent==1
    dependent_v=DP(2:end);
    dependent_vH=DPH;
    dependent_vL=DPL;
elseif Switch_dependent==2
    dependent_v=DG;
    dependent_vH=DGH;
    dependent_vL=DGL;
end
%% whole sample
A1=[ones(DataSize,1) Mt(1:DataSize) mu_t(1:DataSize) DP(1:end-1)];
[X1,STDX1]=lscov(A1,dependent_v); % independent variables:  DP(2:end) DG 
R_squared1=1-sum((dependent_v-X1(1)-X1(2)*Mt(1:DataSize)-X1(3)*mu_t(1:DataSize)-X1(4)*DP(1:end-1)).^2)/sum((dependent_v-mean(dependent_v)).^2);
R_squared_adjusted1=1-(DataSize-1)/(DataSize-4)*sum((dependent_v-X1(1)-X1(2)*Mt(1:DataSize)-X1(3)*mu_t(1:DataSize)-X1(4)*DP(1:end-1)).^2)/sum((dependent_v-mean(dependent_v)).^2);
nlag=12;
RegressionCoef1=nwest(dependent_v,A1,nlag);
Newey_West_t=RegressionCoef1.tstat;


%% high regime
MtH=aligned_Mt_H;
muH=aligned_Yt_H;
MtL=aligned_Mt_L;
muL=aligned_Yt_L;
    NumberH=Number_High;
    A2=[ones(NumberH,1) MtH muH DP_tMINUS1H];
    [X2,STDX2]=lscov(A2,dependent_vH); % independent variables:  DPH DGH 
    R_squared2=1-sum((dependent_vH-X2(1)-X2(2)*MtH-X2(3)*muH-X2(4)*DP_tMINUS1H).^2)/sum((dependent_vH-mean(dependent_vH)).^2);
    R_squared_adjusted2=1-(NumberH-1)/(NumberH-4)*sum((dependent_vH-X2(1)-X2(2)*MtH-X2(3)*muH-X2(4)*DP_tMINUS1H).^2)/sum((dependent_vH-mean(dependent_vH)).^2);
    RegressionCoef2=nwest(dependent_vH,A2,nlag);
    Newey_West_tH=RegressionCoef2.tstat;
%% low regime
    NumberL=DataSize-Number_High;
    A3=[ones(NumberL,1) MtL muL DP_tMINUS1L];
    [X3,STDX3]=lscov(A3,dependent_vL); % independent variables:  DPL DGL 
    R_squared3=1-sum((dependent_vL-X3(1)-X3(2)*MtL-X3(3)*muL-X3(4)*DP_tMINUS1L).^2)/sum((dependent_vL-mean(dependent_vL)).^2);
    R_squared_adjusted3=1-(NumberL-1)/(NumberL-4)*sum((dependent_vL-X3(1)-X3(2)*MtL-X3(3)*muL-X3(4)*DP_tMINUS1L).^2)/sum((dependent_vL-mean(dependent_vL)).^2);
    RegressionCoef3=nwest(dependent_vL,A3,nlag);
    Newey_West_tL=RegressionCoef3.tstat;  



%% Bootstrap p-values
NW_t=zeros(4,2000);
NW_tH=zeros(4,2000);
NW_tL=zeros(4,2000);

for j=1:2000
    seeds=j;
% full sample

    epsilon_t=dependent_v-X1(1)-X1(2)*Mt(1:DataSize)-X1(3)*mu_t(1:DataSize)-X1(4)*DP(1:end-1);
    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);
    X_DP=lscov([ones(DataSize-1,1) DP(1:DataSize-1)],DP(2:DataSize));
    rho1_DP=X_DP(2)+((1+3*X_DP(2))/(DataSize))+(3*(1+3*X_DP(2))/(DataSize)^2);
    rho0_DP=mean(DP(2:DataSize)-rho1_DP*DP(1:DataSize-1));
    phi_DP=DP(2:DataSize)-rho0_DP-rho1_DP*DP(1:DataSize-1);
    randn('state',seeds)
    Wt=normrnd(0,1,DataSize,1); 
    pseudo_R=mean(dependent_v)+epsilon_t.*Wt;
    pseudo_m=zeros(DataSize,1);
    pseudo_m(1)=Mt(1);
    pseudo_mu=zeros(DataSize,1);
    pseudo_mu(1)=mu_t(1);
    pseudo_DP=zeros(DataSize,1);
    pseudo_DP(1)=DP(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);
        pseudo_DP(i)=rho0_DP+rho1_DP*pseudo_DP(i-1)+phi_DP(i-1)*Wt(i-1);
    end
    RegressionCoef4=nwest(pseudo_R,[ones(DataSize,1) pseudo_m pseudo_mu pseudo_DP],nlag);
    NW_t(:,j)=RegressionCoef4.tstat;
% regime H

    epsilon_tH=dependent_vH-X2(1)-X2(2)*MtH-X2(3)*muH-X2(4)*DP_tMINUS1H;
    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);
    X_DPH=lscov([ones(NumberH-1,1) DP_tMINUS1H(1:NumberH-1)],DP_tMINUS1H(2:NumberH));
    rho1_DPH=X_DPH(2)+((1+3*X_DPH(2))/(NumberH))+(3*(1+3*X_DPH(2))/(NumberH)^2);
    rho0_DPH=mean(DP_tMINUS1H(2:NumberH)-rho1_DPH*DP_tMINUS1H(1:NumberH-1));
    phi_DPH=DP_tMINUS1H(2:NumberH)-rho0_DPH-rho1_DPH*DP_tMINUS1H(1:NumberH-1);
    pseudo_RH=mean(dependent_vH(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);
    pseudo_DPH=zeros(NumberH,1);
    pseudo_DPH(1)=DP_tMINUS1H(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);
        pseudo_DPH(i)=rho0_DPH+rho1_DPH*pseudo_DPH(i-1)+phi_DPH(i-1)*Wt(i-1);
    end
    RegressionCoef5=nwest(pseudo_RH,[ones(NumberH,1) pseudo_mH pseudo_muH pseudo_DPH],nlag);
    NW_tH(:,j)=RegressionCoef5.tstat;

% regime L

    epsilon_tL=RtL-X3(1)-X3(2)*MtL-X3(3)*muL-X3(4)*DP_tMINUS1L;
    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);
    X_DPL=lscov([ones(NumberL-1,1) DP_tMINUS1L(1:NumberL-1)],DP_tMINUS1L(2:NumberL));
    rho1_DPL=X_DPL(2)+((1+3*X_DPL(2))/(NumberL))+(3*(1+3*X_DPL(2))/(NumberL)^2);
    rho0_DPL=mean(DP_tMINUS1L(2:NumberL)-rho1_DPL*DP_tMINUS1L(1:NumberL-1));
    phi_DPL=DP_tMINUS1L(2:NumberL)-rho0_DPL-rho1_DPL*DP_tMINUS1L(1:NumberL-1);
    pseudo_RL=mean(dependent_vL(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);
    pseudo_DPL=zeros(NumberL,1);
    pseudo_DPL(1)=DP_tMINUS1L(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);
        pseudo_DPL(i)=rho0_DPL+rho1_DPL*pseudo_DPL(i-1)+phi_DPL(i-1)*Wt(i-1);
    end
    RegressionCoef6=nwest(pseudo_RL,[ones(NumberL,1) pseudo_mL pseudo_muL pseudo_DPL],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_t(4)>0
	     p_full_DP=mean(NW_t(4,:)>=Newey_West_t(4));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	elseif Newey_West_t(4)<0
	     p_full_DP=mean(NW_t(4,:)<=Newey_West_t(4));
	elseif Newey_West_t(4)==0
	     p_full_DP=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_tH(4)>0
	     p_H_DP=mean(NW_tH(4,:)>=Newey_West_tH(4));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	 elseif Newey_West_tH(4)<0
	     p_H_DP=mean(NW_tH(4,:)<=Newey_West_tH(4));
	 elseif Newey_West_tH(4)==0
	     p_H_DP=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
    if Newey_West_tL(4)>0
	     p_L_DP=mean(NW_tL(4,:)>=Newey_West_tL(4));% the percentage of NW_t(:,j) is greater than Newey_West_t(2)
	elseif Newey_West_tL(4)<0
	     p_L_DP=mean(NW_tL(4,:)<=Newey_West_tL(4));
	elseif Newey_West_tL(4)==0
	     p_L_DP=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_full_DP=mean(abs(NW_t(4,:))>=abs(Newey_West_t(4)));
    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_H_DP=mean(abs(NW_tH(4,:))>=abs(Newey_West_tH(4)));
    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)));
    p_L_DP=mean(abs(NW_tL(4,:))>=abs(Newey_West_tL(4)));
end

TX1=Newey_West_t;
TX2=Newey_West_tH;
TX3=Newey_West_tL;
%% save results
Table1=[X1(3)*100 X1(2)*100 R_squared1*100;
    p_full_mu  p_full_m 0;
    TX1(3) TX1(2) 0;
    X2(3)*100 X2(2)*100 R_squared2*100;
    p_H_mu  p_H_m 0;
    TX2(3) TX2(2) 0;
    X3(3)*100 X3(2)*100 R_squared3*100;
    p_L_mu  p_L_m 0;
    TX3(3) TX3(2) 0];
