format long;
clear;
clc;
input_file='investor_Sentiment_Data_v23_POST.xlsx';% load Baker-Wurgler data 
input_sheet='POST monthly';
date=xlsread(input_file,input_sheet,'a102:a647');
%% pdnd
pdnd=xlsread(input_file,input_sheet,'f90:f635');
%% nipo
nipo=xlsread(input_file,input_sheet,'g102:g647');
%% ripo
ripo=xlsread(input_file,input_sheet,'aa90:aa635');
%% turn
%turn=xlsread(input_file,input_sheet,'i90:i635');
turnover=xlsread(input_file,input_sheet,'i30:i635');
turn=nan(606,1);
for t=61:606
turn(t)=log(turnover(t))-sum(log(turnover(t-60:t-1)))/60;
end
turn=turn(61:end);
%% cefd
cefd=xlsread(input_file,input_sheet,'j102:j647');
%% eqti
se=xlsread(input_file,input_sheet,'k102:k647');
sd=xlsread(input_file,input_sheet,'l102:l647');
eqti=se./(se+sd);
IP=xlsread(input_file,input_sheet,'n102:n647');
lagIP=xlsread(input_file,input_sheet,'n101:n646');
gIP=log(IP)-log(lagIP);
Cons=xlsread(input_file,input_sheet,'o102:o647');
lagCons=xlsread(input_file,input_sheet,'o101:o646');
gCons=log(Cons)-log(lagCons);
employ=xlsread(input_file,input_sheet,'t102:t647');
lagemploy=xlsread(input_file,input_sheet,'t101:t646');
gemploy=log(employ)-log(lagemploy);
recess=xlsread(input_file,input_sheet,'s102:s647');
ECON=[gIP gCons gemploy recess];
IP2=xlsread(input_file,input_sheet,'n90:n635');
lagIP2=xlsread(input_file,input_sheet,'n89:n634');
gIP2=log(IP2)-log(lagIP2);
Cons2=xlsread(input_file,input_sheet,'o90:o635');
lagCons2=xlsread(input_file,input_sheet,'o89:o634');
gCons2=log(Cons2)-log(lagCons2);
employ2=xlsread(input_file,input_sheet,'t90:t635');
lagemploy2=xlsread(input_file,input_sheet,'t89:t634');
gemploy2=log(employ2)-log(lagemploy2);
recess2=xlsread(input_file,input_sheet,'s90:s635');
ECON2=[gIP2 gCons2 gemploy2 recess2];
%% sentiment component
X_SENT6=[cefd turn nipo ripo pdnd eqti];
T=546;
X_SENT6=X_SENT6(1:T,:);
%%
load SP500_sentment_new.mat;
load PCA_ECON.mat;
load NONFUND_variable.mat;

TrainingSample=273; 
R=TrainingSample; 
P=T-R; % forecast evaluation period 

% Computing recursive forecasts

    for t=1:P+1;
        results_SENT_1=ols(X_SENT6(1:R+t-1,1),[ones(R+t-1,1) ECON(1:R+(t-1),:)]);
        X_SENT6hat(1:R+t-1,1)=results_SENT_1.yhat;
        X_SENT6res(1:R+t-1,1)=results_SENT_1.resid;
        results_SENT_2=ols(X_SENT6(1:R+t-1,2),[ones(R+t-1,1) ECON2(1:R+(t-1),:)]);
        X_SENT6hat(1:R+t-1,2)=results_SENT_2.yhat;
        X_SENT6res(1:R+t-1,2)=results_SENT_2.resid;
        results_SENT_3=ols(X_SENT6(1:R+t-1,3),[ones(R+t-1,1) ECON(1:R+(t-1),:)]);
        X_SENT6hat(1:R+t-1,3)=results_SENT_3.yhat;
        X_SENT6res(1:R+t-1,3)=results_SENT_3.resid;
        results_SENT_4=ols(X_SENT6(1:R+t-1,4),[ones(R+t-1,1) ECON2(1:R+(t-1),:)]);
        X_SENT6hat(1:R+t-1,4)=results_SENT_4.yhat;
        X_SENT6res(1:R+t-1,4)=results_SENT_4.resid;
        results_SENT_5=ols(X_SENT6(1:R+t-1,5),[ones(R+t-1,1) ECON2(1:R+(t-1),:)]);
        X_SENT6hat(1:R+t-1,5)=results_SENT_5.yhat;
        X_SENT6res(1:R+t-1,5)=results_SENT_5.resid;
        results_SENT_6=ols(X_SENT6(1:R+t-1,6),[ones(R+t-1,1) ECON(1:R+(t-1),:)]);
        X_SENT6hat(1:R+t-1,6)=results_SENT_6.yhat;
        X_SENT6res(1:R+t-1,6)=results_SENT_6.resid;

        X_SENT6res(1:R+t-1,[1 5])=(-1).*X_SENT6res(1:R+t-1,[1 5]);
        X_SENT6res(1:R+(t-1),:)=mvavg(X_SENT6res(1:R+(t-1),:),6);
        X_SENT6res(1:R+(t-1),:)=zscore(X_SENT6res(1:R+(t-1),:));
    
        [Lambda_hat,F_hat(1:R+(t-1),:)]=princomp(zscore(X_SENT6res(1:R+(t-1),:)));
        PCres(1:R+(t-1),t)=F_hat(1:R+t-1,1);
        PCres(1:R+(t-1),t)=zscore(PCres(1:R+(t-1),t));%real time BW index
    end;
    Constraint_level=0.05;
    [ R_squared_OS,R_squared_OSH,R_squared_OSL,p_value_all,p_value_all_h,p_value_all_l] = OOS_A(PCres,TrainingSample,Mt2,y1,DataSize,Rt,Constraint_level);
    [ R_squared_OS_m,R_squared_OSH_m,R_squared_OSL_m,p_value_m,p_value_m_h,p_value_m_l] = OOS_m_A(PCres,TrainingSample,Mt2,y1,DataSize,Rt,Constraint_level);
    [ R_squared_OS_mu,R_squared_OSH_mu,R_squared_OSL_mu,p_value_mu,p_value_mu_h,p_value_mu_l] = OOS_mu_A(PCres,TrainingSample,Mt2,y1,DataSize,Rt,Constraint_level);
    [ R_squared_OS_diff,p_value_diff] = dif_vari_dif_regime_A(PCres,TrainingSample,Mt2,y1,DataSize,Rt,Constraint_level);


Table1=[R_squared_OS_mu*100  R_squared_OS_m*100  R_squared_OS*100  R_squared_OS_diff*100;
    p_value_mu p_value_m p_value_all p_value_diff;
    R_squared_OSH_mu*100    R_squared_OSH_m*100   R_squared_OSH*100   NaN;
    p_value_mu_h p_value_m_h p_value_all_h NaN;
    R_squared_OSL_mu*100  R_squared_OSL_m*100   R_squared_OSL*100   NaN;
    p_value_mu_l p_value_m_l p_value_all_l NaN;];

