
Ndata=141;
S=200;

%Data Moments
datamomvec    = zeros(8,1);
datamomvec(1) = 3.04638;
datamomvec(2) = 1.33218;
datamomvec(3) = 1.24311;
datamomvec(4) = 2.52651;
datamomvec(5) = 0.384526;
datamomvec(6) = 3.05899;
datamomvec(7) = 2.946679;
datamomvec(8) = 1.228861;

%Read in covariance matrix of moments
datcovmat = importdata('datcovmat.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.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.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.01;

mom_mat               = zeros(npar,nmom);
for paramct=1:npar
  xdiff               = PARAMHAT;
  xdiff(paramct)      = PARAMHAT(paramct)*(1+diff_val);
  mom_mat(paramct,:)  = SMM_model_moments(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 5, 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 5, Panel B')
disp('Model v Data Mom')
disp([ num2str(round([baseline_mom_scaled data_mom_scaled (datamomvec-baseline_mom)./DATAMOMSEs],3))])
disp(' ')


%% TABLE 6
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);

nontarget=NaN(Nsim,23);
kvec=NaN(Nsim,2);
xval=PARAMHAT;
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(simct,:), kvec(simct,:)] = SMM_nontarget(xval,shocks,N,T,burn);
end

mean_nontarget = mean(nontarget);
std_nontarget  = std(nontarget);

mom_nontarget_description = ["b(FECQ,FRCQ)","b(FE1Q,FR1Q)","b(FE2Q,FR2Q)","b(FE3Q,FR3Q)","b(FRCQ,FR1Q_{-1})" ...
    "b(FR1Q,FR2Q_{-1})","b(FR2Q,FR3Q_{-1})","b(FEYY,FRYY)","b(FEYY,Outcome)","b(FECQ,FECQ_{-1})"]';

nontarget_table_c1 = round(mean_nontarget(1:10),3)';
nontarget_table_c2 = round(std_nontarget(1:10),3)';
nontarget_table_c3 = [-0.266;-0.145;-0.334;-0.657;-0.131;-0.302;-0.424;-0.236;-0.134;0.147];
nontarget_table_c4 = [0.059;0.073;0.066;0.087;0.057;0.055;0.050;0.137;0.070;0.054];


disp('%%%%%%%%%%%%%%%%%%')
disp('Table 6')
disp('Model v Data Mom')
disp([mom_nontarget_description nontarget_table_c1 nontarget_table_c2 nontarget_table_c3 nontarget_table_c4])
disp(' ')


%% Table 8
pr_update                = [0.999;0.581;1.000;0.552];
implied_friction         = [(1-pr_update(1,1))+pr_update(1,1)*(1-0.800-0.026); (1-pr_update(2,1))+pr_update(2,1)*(1-0.800-0.026);...
    (1-pr_update(3,1))+pr_update(3,1)*(1-0.783-0.028);(1-pr_update(4,1))+pr_update(4,1)*(1-0.783-0.028)];
sticky_info_contribution = [(1-pr_update(1,1))/implied_friction(1,1);(1-pr_update(2,1))/implied_friction(2,1);...
    (1-pr_update(3,1))/implied_friction(3,1);(1-pr_update(4,1))/implied_friction(4,1)];
noisy_info_contribution  = [(pr_update(1,1)*(1-0.800-0.026))/implied_friction(1,1);(pr_update(2,1)*(1-0.800-0.026))/implied_friction(2,1);...
    (pr_update(3,1)*(1-0.783-0.028))/implied_friction(3,1);(pr_update(4,1)*(1-0.783-0.028))/implied_friction(4,1)];


table8_c1=["Quarterly";"Annual";"Quarterly";"Annual"];

disp('%%%%%%%%%%%%%%%%%%')
disp('Table 8')
disp('')
disp(["" "Probabilty of updating" "Implied friction" "Sticky info contribution" "Noisy info contribution";table8_c1 pr_update implied_friction sticky_info_contribution noisy_info_contribution])
disp(' ')

