clear; clc; close all;
cd('/Users/julioortiz/Dropbox/ME/replication/code/Appendix_D');

%% Define estimated parameters
xvec_RE = [2.556 0.434 -0.006 1.663 1.5304];
xvec_OC = [2.556 0.434 -0.006 1.663 1.5304 0.7224];
xvec_DE = [2.556 0.434 -0.006 1.663 1.5304 0.2328];
xvec_ME = [2.556 0.434 -0.006 1.663 1.5304 0.5641];

%% Simulate each model 
N=100;
T=51;
burn=200;
Nsim=2000;

rng(2501)
shocks.v=randn(N,T+burn,Nsim);
shocks.w=randn(1,T+burn,Nsim);

% First run no shock
results_RE_noshock=aggshocksim_fnc_RE(N,T,burn,Nsim,shocks,xvec_RE,0);
results_OC_noshock=aggshocksim_fnc_OC(N,T,burn,Nsim,shocks,xvec_OC,0);
results_DE_noshock=aggshocksim_fnc_DE(N,T,burn,Nsim,shocks,xvec_DE,0);
results_ME_noshock=aggshocksim_fnc_ME(N,T,burn,Nsim,shocks,xvec_ME,0);

% Then run shock
results_RE_shock=aggshocksim_fnc_RE(N,T,burn,Nsim,shocks,xvec_RE,1);
results_OC_shock=aggshocksim_fnc_OC(N,T,burn,Nsim,shocks,xvec_OC,1);
results_DE_shock=aggshocksim_fnc_DE(N,T,burn,Nsim,shocks,xvec_DE,1);
results_ME_shock=aggshocksim_fnc_ME(N,T,burn,Nsim,shocks,xvec_ME,1);


%% Plot IRFs of average one-step ahead FE
results_RE.fe=results_RE_shock.fe - results_RE_noshock.fe;

results_OC.fe=results_OC_shock.fe - results_OC_noshock.fe;

results_DE.fe=results_DE_shock.fe - results_DE_noshock.fe;

results_ME.fe=results_ME_shock.fe - results_ME_noshock.fe;

% now, load the empirical IRF
dat = readmatrix('/Users/julioortiz/Dropbox/ME/replication/results/ushock.csv');
conf = 1; %68 percent CI
beta = dat(1,2:end);
tstat = dat(2,2:end);
beta_se = beta./tstat;
beta_lb = beta - conf*beta_se;
beta_ub = beta + conf*beta_se;

x0=10;
y0=10;
width=750;
height = 11;
aspect = 16/9;
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];

%% Figure D3
w1=23; w2=30; period=1:(w2-w1+1); period=period-2;

figure('units','centimeter','Position',[10 10 1.75*height*aspect 1.75*height],'color','w');
hold on;
patch([[0:6]';flip(0:6)'],[beta_ub(1:7)';flip(beta_lb(1:7)')],[1,1,1]*0.9,'EdgeColor','none');
h5=plot(period(2:end),beta(1:(w2-w1)),'-','Color','k','LineWidth',3);
h6=plot(period(2:end),beta_ub(1:(w2-w1)),'--','Color','k','LineWidth',2);
h7=plot(period(2:end),beta_lb(1:(w2-w1)),'--','Color','k','LineWidth',2);
yline(0,'--');
xlabel('Quarter','interpreter','latex');
ylabel('Forecast error','interpreter','latex');
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex');
hold off
saveas(gcf,'/Users/julioortiz/Dropbox/ME/results/results/FigureD3','epsc')


%% Figure D4
w1=23; w2=30; period=1:(w2-w1+1); period=period-2;
figure('units','centimeter','Position',[10 10 1.75*height*aspect 1.75*height],'color','w');
hold on;
h1=plot(period(2:end),results_RE.fe(w1+1:w2),'-o','Color',myblue,'MarkerSize',10,'LineWidth',3,'MarkerFaceColor',myblue);
h2=plot(period(2:end),results_OC.fe(w1+1:w2),'-x','Color',myred,'MarkerSize',10,'LineWidth',3,'MarkerFaceColor',myred);
h3=plot(period(2:end),results_DE.fe(w1+1:w2),'-s','Color',green,'MarkerSize',10,'LineWidth',3,'MarkerFaceColor',green);
h4=plot(period(2:end),results_ME.fe(w1+1:w2),'-*','Color','m','MarkerSize',10,'LineWidth',3,'MarkerFaceColor','m');
yline(0,'--');
xlabel('Quarter','interpreter','latex');
ylabel('Forecast error','interpreter','latex');
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex');
legend('Rational Expectations','Overconfidence','Diagnostic Expectations','Misspecified Expectations','Location', 'northeast','interpreter','latex'); 
legend('boxoff')
hold off
saveas(gcf,'/Users/julioortiz/Dropbox/ME/replication/results/FigureD4','epsc')



