clear; clc; close all;
cd('/Users/julioortiz/Dropbox/ME/replication');
addpath '/Users/julioortiz/Dropbox/ME/replication/processed_data';

%Loop through subsamples
Nsamp   = 10;
xvec    = NaN(Nsamp,4);
xvec_RE = NaN(Nsamp,1);
xvec_OC = NaN(Nsamp,1);
xvec_DE = NaN(Nsamp,1);
xvec_ME = NaN(Nsamp,1);

fval    = NaN(Nsamp,1);
fval_RE = NaN(Nsamp,1);
fval_OC = NaN(Nsamp,1);
fval_DE = NaN(Nsamp,1);
fval_ME = NaN(Nsamp,1);

AIC_RE = NaN(Nsamp,1);
AIC_OC = NaN(Nsamp,1);
AIC_DE = NaN(Nsamp,1);
AIC_ME = NaN(Nsamp,1);

N         = NaN(Nsamp,1);
subsample = 1:Nsamp;

for j=1:length(subsample)

%% First, solve actuals
dat=readmatrix(strcat('processed_data/rgdp_actuals_subsample',num2str(subsample(j)),'.csv'));

act_rt=dat(:,3);
tStart=tic;
[xvec(j,:),fval(j)]=filterfnc_ar2_H(act_rt);
tEnd=toc(tStart);
fprintf('Fundamental parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

%% Next, solve panel MLE 
dat=readmatrix(strcat('processed_data/spf9219_rgdp','_fefrlfe_subsample',num2str(subsample(j)),'.csv'));
year    = dat(:,2);
quarter = dat(:,3);
fe = dat(:,4);
fr = dat(:,5);
lfe = dat(:,6);
f0data=dat(:,7);
f1data=dat(:,8);

if subsample(j)==10
    T=96;
else
    T=100;
end

N(j) = size(dat(:,1),1)/T;

fe     = reshape(fe,[T,N(j)])';
fr     = reshape(fr,[T,N(j)])';
lfe    = reshape(lfe,[T,N(j)])';
f0data = reshape(f0data,[T,N(j)])';
f1data = reshape(f1data,[T,N(j)])';

clear Y
Y(:,:,1)=fe;
Y(:,:,2)=fr;
Y(:,:,3)=lfe;

tStart=tic;
[xvec_RE(j),fval_RE(j)]=filtervecfnc_RE(Y,N(j),T,xvec(j,2:4)); %RE model
[xvec_OC(j),fval_OC(j)]=filtervecfnc_OC(Y,N(j),T,xvec(j,2:4),xvec_RE(j)); %OC model
[xvec_DE(j),fval_DE(j)]=filtervecfnc_DE(Y,N(j),T,xvec(j,2:4),xvec_RE(j)); %DE model
[xvec_ME(j),fval_ME(j)]=filtervecfnc_ME(Y,N(j),T,xvec(j,2:4),xvec_RE(j)); %ME model
tEnd=toc(tStart);
fprintf('Info/bias parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

n_RGDP=sum(~isnan(fe),'all');

% Model selection exercise
AIC_RE(j) = 2*5 - 2*(-fval_RE(j));
AIC_DE(j) = 2*6 - 2*(-fval_DE(j));
AIC_OC(j) = 2*6 - 2*(-fval_OC(j));
AIC_ME(j) = 2*6 - 2*(-fval_ME(j));

xvec_RE_all = [xvec(j,2:4) xvec_RE(j)];
xvec_OC_all = [xvec(j,2:4) xvec_RE(j) xvec_OC(j)];
xvec_DE_all = [xvec(j,2:4) xvec_RE(j) xvec_DE(j)];
xvec_ME_all = [xvec(j,2:4) xvec_RE(j) xvec_ME(j)];

[Yact_RE,Yhat_RE,Yf1_RE,Yf0_RE,Ycq_RE] = predictfnc_RE(Y,xvec_RE_all); 
[Yact_OC,Yhat_OC,Yf1_OC,Yf0_OC,Ycq_OC] = predictfnc_OC(Y,xvec_OC_all); 
[Yact_DE,Yhat_DE,Yf1_DE,Yf0_DE,Ycq_DE] = predictfnc_DE(Y,xvec_DE_all); 
[Yact_ME,Yhat_ME,Yf1_ME,Yf0_ME,Ycq_ME] = predictfnc_ME(Y,xvec_ME_all); 

f1_RE = Yf1_RE{1}(1:3:end);
f1_OC = Yf1_OC{1}(1:3:end);
f1_DE = Yf1_DE{1}(1:3:end);
f1_ME = Yf1_ME{1}(1:3:end);

for t=2:T    
    f1_RE = [f1_RE ; Yf1_RE{t}(1:3:end)];
    f1_OC = [f1_OC ; Yf1_OC{t}(1:3:end)];
    f1_DE = [f1_DE ; Yf1_DE{t}(1:3:end)];
    f1_ME = [f1_ME ; Yf1_ME{t}(1:3:end)];
end

f1tmp = f1data(:,1);
f1tmp(isnan(f1tmp))=[];
f1vec = f1tmp;
for t=2:T
    f1tmp = f1data(:,t);
    f1tmp(isnan(f1tmp))=[];
    f1vec = [f1vec;f1tmp];
end

RE_o  = f1_RE;
OC_o  = f1_OC;
DE_o  = f1_DE;
ME_o  = f1_ME;
Act_o = f1vec;

lb  = [-Inf zeros(1,4)];
ub  = [Inf ones(1,4)];
Aeq = [0 ones(1,4)];
beq =  1;
[coef,resnorm,residual,exitflag,output,lambda] = lsqlin([ones(size(RE_o,1),1) RE_o OC_o DE_o ME_o],Act_o,[],[],Aeq,beq,lb,ub);
beta(:,j) = coef;

end

%% FIGURE
% - Colors -
green = [0.0 0.6 0.0];
blue  = [0.0 0.0 0.6];
red   = [0.6 0.0 0.0];
myred= [0.6350 0.0780 0.1840];
myblue = [0 0.4470 0.7410];

% - Line width -
heavy = 3;
light = 0.8;

height = 11;
aspect = 16/9;

dAIC = ([AIC_RE AIC_OC AIC_DE AIC_ME] - min([AIC_RE AIC_OC AIC_DE AIC_ME],[],2));
AICw=exp(-0.5.*dAIC)./sum(exp(-0.5.*dAIC),2);


figure('units','centimeter','Position',[10 10 1.75*height*aspect 2.25.*height],'color','w')
subplot(4,1,1)
    plot(beta(2,:),'-x','MarkerSize',10,'LineWidth',heavy,'Color','k')
    ylim([0 1])
    yline(0.50,'Color','k','LineStyle','--')
    title('Panel A: Rational expectations','Interpreter','latex');
    hold off;
    set(gca,'FontSize',16,'TickLabelInterpreter','latex','XtickLabel',{'1993Q4','1996Q4','1999Q4'...
        ,'2002Q4','2005Q4','2008Q4','2011Q4','2014Q4','2017Q4','2019Q4'})
    ylabel('$w_k$','interpreter','latex')
subplot(4,1,2)
    plot(beta(3,:),'-x','MarkerSize',10,'LineWidth',heavy,'Color','k')
    ylim([0 1])
    yline(0.50,'Color','k','LineStyle','--')
    title('Panel B: Overconfident Expectations','Interpreter','latex');
    set(gca,'FontSize',16,'TickLabelInterpreter','latex','XtickLabel',{'1993Q4','1996Q4','1999Q4'...
        ,'2002Q4','2005Q4','2008Q4','2011Q4','2014Q4','2017Q4','2019Q4'})
    ylabel('$w_k$','interpreter','latex')
subplot(4,1,3)
    plot(beta(4,:),'-x','MarkerSize',10,'LineWidth',heavy,'Color','k')
    ylim([0 1])
    yline(0.50,'Color','k','LineStyle','--')
    title('Panel C: Diagnostic Expectations','Interpreter','latex');
    set(gca,'FontSize',16,'TickLabelInterpreter','latex','XtickLabel',{'1993Q4','1996Q4','1999Q4'...
        ,'2002Q4','2005Q4','2008Q4','2011Q4','2014Q4','2017Q4','2019Q4'})
    ylabel('$w_k$','interpreter','latex')
subplot(4,1,4)
    plot(beta(5,:),'-x','MarkerSize',10,'LineWidth',heavy,'Color','k')
    ylim([0 1])
    yline(0.50,'Color','k','LineStyle','--')
    title('Panel D: Misspecified Expectations','Interpreter','latex');
    set(gca,'FontSize',16,'TickLabelInterpreter','latex','XtickLabel',{'1993Q4','1996Q4','1999Q4'...
        ,'2002Q4','2005Q4','2008Q4','2011Q4','2014Q4','2017Q4','2019Q4'})
    ylabel('$w_k$','interpreter','latex')
    xlabel('Date','interpreter','latex','Fontsize',16)
saveas(gcf,'/Users/julioortiz/Dropbox/ME/replication/results/Figure3','epsc')


