Ndata=94;
S=100;

datamomvec=zeros(8,1);
datamomvec(1) = 2.36651;
datamomvec(2) = 1.24554;
datamomvec(3) = 1.12321;
datamomvec(4) = 1.50027;
datamomvec(5) = 0.405598;
datamomvec(6) = 2.30464;
datamomvec(7) = 2.276894;
datamomvec(8) = 0.9388027;

%Read in covariance matrix of moments
datcovmat = importdata('datcovmat0019.txt');
datcovmat = datcovmat*Ndata;
DATAMOMSEs = sqrt(diag(datcovmat/Ndata));

%% Read in results
%read in estimated parameters
opt = {'CommentStyle','%', 'CollectOutput',true};
[fid,msg] = fopen('param_est_0019.txt','rt');
assert(fid>=3,msg)
C = textscan(fid, '%f', opt{:});
fclose(fid);
params = C{1};
PARAMHAT = params;

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

nmom = length(datamomvec);
npar = length(PARAMHAT);

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

mom_mat               = zeros(npar,nmom);
for paramct=1:npar
  xdiff               = PARAMHAT;
  xdiff(paramct)      = PARAMHAT(paramct)*(1+diff_val);
  mom_mat(paramct,:)  = SMM_model_moments_0019(xdiff,datamomvec);
end

mom_diff_mat          = mom_mat - repmat(baseline_mom',npar,1);
mom_jac               = zeros(npar,nmom);

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

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

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

disp('%%%%%%%%%%%%%%%%%%')
disp('Table D6, Panel A')
disp('Parameter Estimates')
disp([Description num2str([round(PARAMHAT,3) round(PARAMSE,3)])])
disp(' ')

%now, scale the model & data moments
baseline_mom_scaled    = baseline_mom;
baseline_mom_scaled(1) = baseline_mom(1)^0.5; 
baseline_mom_scaled(2) = baseline_mom(2)/sqrt(baseline_mom(3)*baseline_mom(1)); 
baseline_mom_scaled(3) = baseline_mom(3)^0.5; 
baseline_mom_scaled(4) = baseline_mom(4)^0.5; 
baseline_mom_scaled(5) = baseline_mom(5)/sqrt(baseline_mom(4)*baseline_mom(6)); 
baseline_mom_scaled(6) = baseline_mom(6)^0.5; 
baseline_mom_scaled(7) = baseline_mom(7)^0.5;
baseline_mom_scaled(8) = baseline_mom(8)^0.5;

data_mom_scaled    = datamomvec;
data_mom_scaled(1) = datamomvec(1)^0.5; 
data_mom_scaled(2) = datamomvec(2)/sqrt(datamomvec(3)*datamomvec(1)); 
data_mom_scaled(3) = datamomvec(3)^0.5; 
data_mom_scaled(4) = datamomvec(4)^0.5; 
data_mom_scaled(5) = datamomvec(5)/sqrt(datamomvec(4)*datamomvec(6)); 
data_mom_scaled(6) = datamomvec(6)^0.5; 
data_mom_scaled(7) = datamomvec(7)^0.5;
data_mom_scaled(8) = datamomvec(8)^0.5;


disp('%%%%%%%%%%%%%%%%%%')
disp('Table D6, Panel B')
disp('Model v Data Mom')
disp([ num2str(round([baseline_mom_scaled data_mom_scaled (datamomvec-baseline_mom)./DATAMOMSEs],3))])
disp(' ')
