format long;
load CRSP_SENT_value_weighted_new.mat  % orthogonal sentiment index


P120=0.01; P210=0.1; psi10=1; psi20=-1; sigma_X0=0.5;
A0=[P120  P210  psi10 psi20 sigma_X0]; %give the initial values of parameters for search

L_bound=[  0;  0;-2;-2;0];
U_bound=[0.5;0.5; 2; 2;3];
options = optimset('Display','iter','TolX',1e-8,'TolFun',1e-6,'MaxIter',10000, 'MaxFunEvals', 10000,'Algorithm','interior-point');% this is for MATLAB R2013a; for latest MATLAB version, can use options = optimset('Display','iter') 
%options = optimset('Display','iter');%,'TolX',1e-8,'TolFun',1e-6,'MaxIter',10000, 'MaxFunEvals', 10000,'Algorithm','interior-point');
ff = @(u) Log_likelihood_fun_Sentiment(u,Xt,DataSize); 
[A_text,F_value]=fmincon(ff,A0,[],[],[],[],L_bound,U_bound,[],options);

    P12=A_text(1);
    P21=A_text(2);
    psi1=A_text(3);
    psi2=A_text(4);
    sigma_X=A_text(5);
    P11=1-P12;
    P22=1-P21;
    
 %% pi(t)
 
    f1=zeros(DataSize,1);
    f2=zeros(DataSize,1);
    pi1=zeros(DataSize+1,1); % pi1(t) stands for \pi(t-\Delta t|t-\Delta t), so pi1(1) stands for \pi(0|0)
    pi2=zeros(DataSize,1); % pi2(t) stands for \pi(t|t-\Delta t)
    pi1(1)=P21/(P12+P21);
    
    for j=1:DataSize
        f1(j)=1/sqrt(2*pi)/sigma_X*exp(-(Xt(j)-psi1)^2/2/sigma_X^2);
        f2(j)=1/sqrt(2*pi)/sigma_X*exp(-(Xt(j)-psi2)^2/2/sigma_X^2);
        pi2(j)=pi1(j)*P11+(1-pi1(j))*P21;
        pi1(j+1)=f1(j)*pi2(j)/(f1(j)*pi2(j)+f2(j)*(1-pi2(j)));
    end
    
    high_m=0;
    for i=2:547
        if pi1(i)>=0.5
            high_m=high_m+1;
        end
    end
    

%% the high-low sentiment month defined in Stambaugh_Yu_Yuan_JFE_12

aaa = sort(Xt, 'descend');
RegimeSYY=zeros(546,1);
for i=1:546
    if Xt(i)<aaa(273);
        RegimeSYY(i)=-9;
    else
        RegimeSYY(i)=3;
    end
end

%save('Para_sentiment.mat','P12','P21','psi1','psi2','sigma_X','pi1','RegimeSYY','high_m');

%% plot of sentiment regime

figure(1) % Markov regime model
plot(pi1(2:end),'k-','LineWidth',3)
%xlabel('t','FontSize',12,'FontWeight','bold')
ylabel('High regime probability','FontSize',12,'FontWeight','bold','Rotation',90)%,'Position',[-50 0.88])

figure(2) % Stambaugh_Yu_Yuan_JFE_12
plot((sign(Xt-mean(Xt))+1)/2,'b-','LineWidth',1.5)
xlabel('t','FontSize',12,'FontWeight','bold')
ylabel('regime','FontSize',12,'FontWeight','bold','Rotation',90)%,'Position',[-50 0.88])

figure(3)
plot((sign(Xt-mean(Xt))+1)/2,'r*','LineWidth',1)
hold on
plot(pi1(2:end),'b-','LineWidth',1.5)
hold off
xlabel('t','FontSize',12,'FontWeight','bold')
legend('SYY','\pi_t','Location','Northeast')

x=1:546;
y1=ones(1,546)*(-3);
y2=ones(1,546);     % regime switching model
for i=1:546
    if pi1(i+1)<0.5;
        y2(i)=-9;
    else
        y2(i)=3;
    end
end
X=[x,fliplr(x)];                %#create continuous x value array for plotting
Y=[y1,fliplr(y2)];              %#create y values for out and then back
figure(4)
fill(X,Y,[1 1 1] * 0.8); 
hold on
plot(Xt,'k-','LineWidth',3)
hold off
%xlabel('t','FontSize',12,'FontWeight','bold')
ylabel('Sentiment','FontSize',12,'FontWeight','bold','Rotation',90)%,'Position',[-50 0.88])
ylim([-3 3])

load Recession_Expansion.mat
figure(7)
fill(X,Y,[1 1 1] * 0.8);
hold on
plot(REt*2-1,'k*','LineWidth',0.5)
hold off
ylim([0 1])




Y=[y1,fliplr(RegimeSYY')];              %SYY model
figure(5)
fill(X,Y,[1 1 1] * 0.8); 
hold on
plot(Xt,'k-','LineWidth',3)
hold off
%xlabel('t','FontSize',12,'FontWeight','bold')
ylabel('Sentiment','FontSize',12,'FontWeight','bold','Rotation',90)
ylim([-3 3])

correlation=corr(RegimeSYY,(sign(pi1(2:end)-mean(pi1(2:end)))+1)/2);

