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

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

%% Comparative Statics
stepvec =[-0.1 -0.05 -0.025 0.00 0.025 0.05 0.1];
stepvec2=[-0.25 -0.20 -0.15 -0.10 -0.05 -0.025 0.0];
matmom=zeros(length(baseline_mom),length(PARAMHAT),length(stepvec));
for stepct=1:length(stepvec)
    stepval=stepvec(stepct);
    for thetact=1:length(PARAMHAT)
        if thetact==5
            theta=PARAMHAT;
            theta(thetact) = PARAMHAT(thetact)*(1+stepval);
            matmom(:,thetact,stepct)=SMM_model_moments(theta,datamomvec);
        else
            theta=PARAMHAT;
            theta(thetact) = PARAMHAT(thetact)*(1+stepval);
            matmom(:,thetact,stepct)=SMM_model_moments(theta,datamomvec);
        end
    end
end 

% Transform matmom to std and corr
momtrans        = matmom;
momtrans(1,:,:) = matmom(1,:,:).^0.5;
momtrans(2,:,:) = matmom(2,:,:)./sqrt(matmom(3,:,:).*matmom(1,:,:));
momtrans(3,:,:) = matmom(3,:,:).^0.5;
momtrans(4,:,:) = matmom(4,:,:).^0.5;
momtrans(5,:,:) = matmom(5,:,:)./sqrt(matmom(4,:,:).*matmom(6,:,:));
momtrans(6,:,:) = matmom(6,:,:).^0.5;
momtrans(7,:,:) = matmom(7,:,:).^0.5;
momtrans(8,:,:) = matmom(8,:,:).^0.5;

%Figure
f1a=figure;
subplot(2,2,1)
    hold on;
    X=stepvec*100;
    y=squeeze(momtrans(2,1,:)); 
	y=((y-y(4))/(abs(y(4))))*100;
    h1=plot(X,y,'LineWidth',3,'Color','k');
    set(gca,'FontSize',14,'TickLabelInterpreter', 'latex')
    ylabel('Percent','interpreter','latex')
    xlabel('','interpreter','latex')
    title('Panel 1: $\rho$','interpreter','latex')
    legend('Corr(f0q,f0y)','interpreter','latex','Location','northwest');
    legend boxoff;
    hold off;
 subplot(2,2,2)  
    hold on;
    X=stepvec*100;
    y1=squeeze(momtrans(3,2,:)); 
	y1=((y1-y1(4))/(abs(y1(4))))*100;
    h2=plot(X,y1,'LineWidth',3,'Color','k');
    set(gca,'FontSize',14,'TickLabelInterpreter', 'latex')
    ylabel('','interpreter','latex')
    xlabel('','interpreter','latex')
    title('Panel 2: $\sigma_w$','interpreter','latex')
    legend('Std(f0y)','interpreter','latex','Location','northwest');
    legend boxoff;
    hold off;
  subplot(2,2,3) 
    hold on;
    X=stepvec*100;
    y=squeeze(momtrans(6,3,:)); 
	y=((y-y(4))/(abs(y(4))))*100;
    h3=plot(X,y,'LineWidth',3,'Color','k');
    set(gca,'FontSize',14,'TickLabelInterpreter', 'latex')
    ylabel('Percent','interpreter','latex')
    xlabel('','interpreter','latex')
    title('Panel 3: $\sigma_e$','interpreter','latex')
    legend('Std(fe0q$_{-1}$)','interpreter','latex','Location','northwest');
    legend boxoff;
    hold off;
  subplot(2,2,4) 
    hold on;
    X=stepvec*100;
    y1=squeeze(momtrans(5,4,:)); 
	y1=((y1-y1(4))/(abs(y1(4))))*100;
    h4=plot(X,y1,'LineWidth',3,'Color','k');
    set(gca,'FontSize',14,'TickLabelInterpreter', 'latex')
    ylabel('','interpreter','latex')
    xlabel('','interpreter','latex')
    title('Panel 4: $\sigma_v$','interpreter','latex')
    legend('Corr(fr0q,fe0q$_{-1}$)','interpreter','latex','Location','northwest');
    legend boxoff;
    hold on;
saveas(gcf,'tables_and_figures/FigureC1_PanelA','epsc')
    
 f1b=figure;  
 subplot(2,2,1) 
    hold on;
    X=stepvec2*100;
    y=squeeze(momtrans(7,5,:)); 
	y=((y-y(4))/(abs(y(4))))*100;
    h5=plot(X,y,'LineWidth',3,'Color','k');
    set(gca,'FontSize',14,'TickLabelInterpreter', 'latex')
    ylabel('Percent','interpreter','latex')
    xlabel('Percent of estimated value','interpreter','latex') 
    title('Panel 5: $q$','interpreter','latex')
    legend('MSEQ','interpreter','latex','Location','northeast');
    legend boxoff;
    hold off;
subplot(2,2,2) 
    hold on;
    X=stepvec*100;
    y=squeeze(momtrans(8,6,:)); 
	y=((y-y(4))/(abs(y(4))))*100;
    h6=plot(X,y,'LineWidth',3,'Color','k');
    set(gca,'FontSize',14,'TickLabelInterpreter', 'latex')
    xlabel('Percent of estimated value','interpreter','latex') 
    title('Panel 6: $p$','interpreter','latex')
    legend('MSEY','interpreter','latex','Location','northeast');
    legend boxoff;
    hold off;    
saveas(gcf,'tables_and_figures/FigureC1_PanelB','epsc')

