Ndata=141;
S=200;

datamomvec    = zeros(10,1);
datamomvec(1) = 3.08465;
datamomvec(2) = 1.35315;
datamomvec(3) = 1.25521;
datamomvec(4) = 2.54147;
datamomvec(5) = 0.370018;
datamomvec(6) = 3.07079;
datamomvec(7) = 2.962625;
datamomvec(8) = 1.224991;
datamomvec(9) = 3.02269;
datamomvec(10)= -0.784446;


datcovmat = importdata('datcovmat_DE.txt');
datcovmat = datcovmat(1:10,1:10);
datcovmat = datcovmat*Ndata;

nmom = length(datamomvec);
npar = 7;


opt = {'CommentStyle','%', 'CollectOutput',true};
[fid,msg] = fopen('param_est_DE_uncon.txt','rt');
assert(fid>=3,msg)
C = textscan(fid, '%f', opt{:});
fclose(fid);
params = C{1};
PARAMHAT_uncon = params;

%read in model moments at estimated parameters
opt = {'CommentStyle','%', 'CollectOutput',true};
[fid,msg] = fopen('model_mom_minimized_DE_uncon.txt','rt');
assert(fid>=3,msg)
C = textscan(fid, '%f', opt{:});
fclose(fid);
model_mom = C{1};
baseline_mom_uncon = model_mom;

opt = {'CommentStyle','%', 'CollectOutput',true};
[fid,msg] = fopen('param_est_DE_con.txt','rt');
assert(fid>=3,msg)
C = textscan(fid, '%f', opt{:});
fclose(fid);
params = C{1};
PARAMHAT_con = params;

%read in model moments at estimated parameters
opt = {'CommentStyle','%', 'CollectOutput',true};
[fid,msg] = fopen('model_mom_minimized_DE_con.txt','rt');
assert(fid>=3,msg)
C = textscan(fid, '%f', opt{:});
fclose(fid);
model_mom = C{1};
baseline_mom_con = model_mom;


%% Do the SMM SE calculation
diff_val=-0.01;

%% Unconstrained
mom_mat = zeros(npar,nmom);
for paramct=1:npar
  xdiff = PARAMHAT_uncon;
  xdiff(paramct) = PARAMHAT_uncon(paramct)*(1+diff_val);
  mom_mat(paramct,:) = SMM_model_moments_DE_uncon(xdiff,datamomvec);
end

mom_diff_mat = mom_mat - repmat(baseline_mom_uncon',npar,1);

mom_jac = zeros(npar,nmom);

for paramct=1:npar
   mom_jac(paramct,:) = mom_diff_mat(paramct,:)/(PARAMHAT_uncon(paramct)*diff_val);
end

%now, construct the asymptotic std errors
OMEGA = (1+(1/S))* eye(npar)/( mom_jac * inv(datcovmat) * mom_jac' );
PARAMSE_uncon = sqrt(diag(OMEGA/Ndata));

%% Constrained
mom_mat = zeros(npar-1,nmom);
for paramct=1:npar-1
  xdiff = PARAMHAT_con;
  xdiff(paramct) = PARAMHAT_con(paramct)*(1+diff_val);
  mom_mat(paramct,:) = SMM_model_moments_DE_con(xdiff,datamomvec);
end

mom_diff_mat = mom_mat - repmat(baseline_mom_con',npar-1,1);

mom_jac = zeros(npar-1,nmom);

for paramct=1:npar-1
   mom_jac(paramct,:) = mom_diff_mat(paramct,:)/(PARAMHAT_con(paramct)*diff_val);
end

%now, construct the asymptotic std errors
OMEGA = (1+(1/S))* eye(npar-1)/( mom_jac * inv(datcovmat) * mom_jac' );
PARAMSE_con = sqrt(diag(OMEGA/Ndata));

%% Table D2
Description = ["Persistence of latent state","State innovation dispersion","Public signal noise",...
    "Private signal noise","Probability of quarterly update","Probability of annual update","Diagnosticiy"]';

disp('%%%%%%%%%%%%%%%%%%')
disp('Table D2')
disp(num2str([PARAMHAT_uncon PARAMSE_uncon [PARAMHAT_con;0] [PARAMSE_con;NaN]]))
disp(' ')


%% Simulating bgms coefficient 

Nsim=200;
N=141;
T=244;
burn=50;

% Draw shocks
shocks=[];
rng(5301);
v=randn(Nsim,N,T);%private signal noise
e=randn(Nsim,T+burn);%public signal noise
w=randn(Nsim,T+burn);%state innovations
%for probabilties, draw U[0,1] and compare to (p,q)
pdraw=rand(Nsim,N,T);
qdraw=rand(Nsim,N,T);


% Counterfactual: Unconstrained 
nontarget_uncon  = NaN(Nsim,3);
xval_uncon       = PARAMHAT_uncon;
for simct=1:Nsim
    shocks.v               = squeeze(v(simct,:,:));
    shocks.e               = squeeze(e(simct,:));
    shocks.w               = squeeze(w(simct,:));
    shocks.pdraw           = squeeze(pdraw(simct,:,:));
    shocks.qdraw           = squeeze(qdraw(simct,:,:));
    nontarget_uncon(simct,:) = SMM_model_DE_uncon_nontarget(xval_uncon,shocks,N,T,burn);
end

% Counterfactual: Constrained
nontarget_con  = NaN(Nsim,3);
xval_con       = PARAMHAT_con;
xval_con(7)    = 0;
for simct=1:Nsim
    shocks.v               = squeeze(v(simct,:,:));
    shocks.e               = squeeze(e(simct,:));
    shocks.w               = squeeze(w(simct,:));
    shocks.pdraw           = squeeze(pdraw(simct,:,:));
    shocks.qdraw           = squeeze(qdraw(simct,:,:));
    nontarget_con(simct,:) = SMM_model_DE_uncon_nontarget(xval_con,shocks,N,T,burn);
end

% - 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;

X = categorical({'Errors on revisions','Revision autocorrelation', 'Error on outcomes'});
X = reordercats(X,{'Errors on revisions','Revision autocorrelation', 'Error on outcomes'});

Y_uncon      = mean(nontarget_uncon);
Y_con        = mean(nontarget_con);

%% Figure D1
Y_p = Y_uncon - Y_con;
Y   = [Y_con;Y_p];
figure('units','centimeter','Position',[10 10 1.75*height*aspect 1.75*height],'color','w')
hold on;
xlim([0 0.22]);
b1=bar(X,Y,'stacked','BarWidth',0.9);
set(gca,'FontSize',16,'TickLabelInterpreter', 'latex','ycolor','k');
set(b1,{'FaceColor'},{[0.64,0.08,0.18];[0.00,0.45,0.74]},{'FaceAlpha'},{0.3});
ylabel('Percent','interpreter','latex')
legend(b1,'Annual smoothing','Diagnostic expectations', 'interpreter','latex','Location','southeast');
legend boxoff;
saveas(gcf,fullfile(figure_folder,'figureD1.eps'),'epsc');
