clear all
clc

warning off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setting Empirical Moments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MOMS_hp_data =...
[2.32;...
 1.31;...
 3.79;...
 0.78;...
 0.81;...
 0.92;...
 0.95;...
 -0.82;...
 0.42;...
 0.54;...
 0.56;...
 1.07;...
 0.74;...
 -0.45];

MOMS_hp_data_corrs_yF_hI = ...
[-0.2009	0.0707	0.3322;...
 -0.3095	-0.0481	0.22;...
 -0.4064	-0.1607	0.1067;...
 -0.5245	-0.3057	-0.049;...
 -0.6392	-0.4562	-0.2243;...
 -0.6987	-0.5356	-0.3196;...
 -0.669	    -0.4927	-0.264;...
 -0.66	    -0.4784	-0.2441;...
 -0.6368	-0.4449	-0.2011];

% Setting entire vector of moments for quadratic function
MOMS_hp_data_all = [MOMS_hp_data;...
                    MOMS_hp_data_corrs_yF_hI(1:4,2);...
                    MOMS_hp_data_corrs_yF_hI(6:9,2)];
    
% Grid for omega between (0,1) with equally spaced grids (except few cases
% to avoid numerical indeterminacy of the solution)
omega_grid = [0.05,  0.10,  0.15,  0.201, 0.251, 0.30,  0.349, 0.399, 0.449,  0.506, 0.551, 0.601, 0.651, 0.70,  0.752, 0.802, 0.853,  0.901, 0.95];
[x,size_omega] = size(omega_grid);

% Grid for sigma_a
sigma_a_first = 0.0001; sigma_a_last  = 0.01;
sigma_a_num_grids = 1; 
sigma_a_increase = (sigma_a_last-sigma_a_first)/sigma_a_num_grids;
sigma_a_grid = sigma_a_first:sigma_a_increase:sigma_a_last;
[x,size_sigma_a] = size(sigma_a_grid);

% Grid for sigma_g
sigma_g_first = 0.0001; sigma_g_last  = 0.01;
sigma_g_num_grids = 1; 
sigma_g_increase = (sigma_g_last-sigma_g_first)/sigma_g_num_grids;
sigma_g_grid = sigma_g_first:sigma_g_increase:sigma_g_last;
[x,size_sigma_g] = size(sigma_g_grid);


 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calibration Benchmark Case 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calibrating some of the other key parameter values. Same calibration as
% in the benchmark case
DELTA_F = 0.05;
DELTA_I = 0.05;
TFPF_to_TFPI = 2.1901;
ALPHA_I = 0.2;
MU = 1.006;
elasticity_F_I = 8; 
E     = 1-(1/elasticity_F_I); 
GAMMA = (MU^(-0.35+ALPHA_I)/TFPF_to_TFPI)^(1/ALPHA_I); % This calibration changes relative to the benchmark case. See appendix.
TAUN  = .1142;  
TAUY  = .07223; 
RHOA  = 0.94;
RHOG  = 0.72;

CHIF = 0;
CHII=0;


%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid - moment matching
%%%%%%%%%%%%%%%%%%%%%%%%%

Mat_diffs = 100*ones(size_omega,size_sigma_a,size_sigma_g);
num_cells = size_omega*size_sigma_a*size_sigma_g;

for cont_omega   = 1:size_omega
    cont_omega
    size_omega
end
for cont_sigma_a = 1:size_sigma_a    
for cont_sigma_g = 1:size_sigma_g  
    OMEGA_a = omega_grid(1,cont_omega);
    OMEGA_g = OMEGA_a; 
    SIGMA_a = sigma_a_grid(1,cont_sigma_a); 
    SIGMA_g = sigma_g_grid(1,cont_sigma_g); 

    
%%%%%%%%%%%%%%%%%%%
% Solving Model
%%%%%%%%%%%%%%%%%%%

[fx,fxp,fy,fyp] = benchmark_model(RHOA,SIGMA_a,RHOG,SIGMA_g,TAUN,TAUY,E,GAMMA,OMEGA_a,OMEGA_g,DELTA_F,DELTA_I,ALPHA_I,CHIF,CHII);
 [gx,hx] = gx_hx(fy,fx,fyp,fxp); 


% Analytical-Model-Based HP Second Moments
        varshock      = zeros(size(hx,1));
        varshock(4,4) = SIGMA_a^2; 
        varshock(5,5) = SIGMA_g^2; 
        [var_y,var_x] = mom(gx,hx,varshock,0); %Compute Variance/covariance matrix
        [var_y1,var_x1] = mom(gx,hx,varshock,1); %First-order autocovariance
        table_modelmoms(1:13,1) = (diag(var_y(1:13,1:13)).^(1/2))*100;
        table_modelmoms(1:13,2) = diag(var_y1(1:13,1:13)) ./ diag(var_y(1:13,1:13));
        table_modelmoms(1:13,3) = var_y(1:13,1) ./ (diag(var_y(1:13,1:13)).^(1/2)) /var_y(1,1)^(1/2);

        % Table of serial correlation Labor and Output
        corr_yFt_hIt = var_y(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        corr_yFt_hIt_plus1 = var_y1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        [var_y2,var_x2] = mom(gx,hx,varshock,2); 
        corr_yFt_hIt_plus2 = var_y2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        [var_y3,var_x3] = mom(gx,hx,varshock,3); 
        corr_yFt_hIt_plus3 = var_y3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        [var_y4,var_x4] = mom(gx,hx,varshock,4); 
        corr_yFt_hIt_plus4 = var_y4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        [var_ym1,var_xm1] = mom(gx,hx,varshock,-1); 
        corr_yFt_hIt_minus1 = var_ym1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        [var_ym2,var_xm2] = mom(gx,hx,varshock,-2); 
        corr_yFt_hIt_minus2 = var_ym2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        [var_ym3,var_xm3] = mom(gx,hx,varshock,-3); 
        corr_yFt_hIt_minus3 = var_ym3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
        [var_ym4,var_xm4] = mom(gx,hx,varshock,-4); 
        corr_yFt_hIt_minus4 = var_ym4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));

        MOMS_hp_model = ...
         [table_modelmoms(1,1);                      %  std(yF)*100       %1
          table_modelmoms(2,1)/table_modelmoms(1,1); %  std(cF)/std(yF)   %2
          table_modelmoms(3,1)/table_modelmoms(1,1); %  std(ivF)/std(yF); %3
          table_modelmoms(7,1)/table_modelmoms(1,1); %  std(tby)/std(yF); %4
          table_modelmoms(1,2);                      %  corr(yF_t,yF_t-1);%5
          table_modelmoms(2,3);                      %  corr(yF,cF);      %6  
          table_modelmoms(3,3);                      %  corr(yF,ivF);     %7
          table_modelmoms(7,3);                      %  corr(yF,tby);     %8
          table_modelmoms(6,1)/table_modelmoms(1,1); %  std(h)/std(yF);   %9
          table_modelmoms(6,3);                      %  corr(yF,h);       %10      
          table_modelmoms(4,1)/table_modelmoms(1,1); %  std(hF)/std(yF);  %11
          table_modelmoms(5,1)/table_modelmoms(1,1); %  std(hI)/std(yF);  %12
          table_modelmoms(4,3);                      %  corr(yF,hF);      %13        
          table_modelmoms(5,3);                      %  corr(yF,hI);      %14
          corr_yFt_hIt_minus4;                       %  corr(yF_t,hI_t-4);%15
          corr_yFt_hIt_minus3;                       %  corr(yF_t,hI_t-3);%16
          corr_yFt_hIt_minus2;                       %  corr(yF_t,hI_t-2);%17
          corr_yFt_hIt_minus1;                       %  corr(yF_t,hI_t-1);%18
          corr_yFt_hIt_plus1;                        %  corr(yF_t,hI_t+1);%19
          corr_yFt_hIt_plus2;                        %  corr(yF_t,hI_t+2);%20
          corr_yFt_hIt_plus3;                        %  corr(yF_t,hI_t+3);%21
          corr_yFt_hIt_plus4];                       %  corr(yF_t,hI_t+4);%22

         diff = ((MOMS_hp_model-MOMS_hp_data_all)./MOMS_hp_data_all);
         diff_1 = diff'*diff;

        Mat_diffs(cont_omega,cont_sigma_a,cont_sigma_g) = diff_1;
        end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computing moments of the lowest value in the grid.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[mini,posi] = min(Mat_diffs(:));
[loc_omega,loc_sigma_a,loc_sigma_g] = ind2sub([size_omega size_sigma_a size_sigma_g],posi);
omega_a_est = omega_grid(loc_omega);     % Calibrated omega presented in upper panel of Table 6, Column labelled "extension 2".
omega_g_est = omega_a_est;
sigma_a_est = sigma_a_grid(loc_sigma_a); % Calibrated sigma_a presented in upper panel of Table 6, Column labelled "extension 2".
sigma_g_est = sigma_g_grid(loc_sigma_g); % Calibrated sigma_g presented in upper panel of Table 6, Column labelled "extension 2".
[fx,fxp,fy,fyp] = benchmark_model(RHOA,SIGMA_a,RHOG,SIGMA_g,TAUN,TAUY,E,GAMMA,OMEGA_a,OMEGA_g,DELTA_F,DELTA_I,ALPHA_I,CHIF,CHII);
 [gx,hx] = gx_hx(fy,fx,fyp,fxp); 




% Analytical-Model-Based HP Second Moments
% Table of second moments
varshock      = zeros(size(hx,1));
varshock(4,4) = sigma_a_est^2; 
varshock(5,5) = sigma_g_est^2; 
[var_y,var_x] = mom(gx,hx,varshock,0); %Compute Variance/covariance matrix
[var_y1,var_x1] = mom(gx,hx,varshock,1); %First-order autocovariance
table_modelmoms(1:13,1) = (diag(var_y(1:13,1:13)).^(1/2))*100;
table_modelmoms(1:13,2) = diag(var_y1(1:13,1:13)) ./ diag(var_y(1:13,1:13));
table_modelmoms(1:13,3) = var_y(1:13,1) ./ (diag(var_y(1:13,1:13)).^(1/2)) /var_y(1,1)^(1/2);

% Table of serial correlation
corr_yFt_hIt = var_y(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
corr_yFt_hIt_plus1 = var_y1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y2,var_x2] = mom(gx,hx,varshock,2); %Second-order autocovariance
corr_yFt_hIt_plus2 = var_y2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y3,var_x3] = mom(gx,hx,varshock,3); %Third-order autocovariance
corr_yFt_hIt_plus3 = var_y3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_y4,var_x4] = mom(gx,hx,varshock,4); %Fourth-order autocovariance
corr_yFt_hIt_plus4 = var_y4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym1,var_xm1] = mom(gx,hx,varshock,-1); 
corr_yFt_hIt_minus1 = var_ym1(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym2,var_xm2] = mom(gx,hx,varshock,-2); 
corr_yFt_hIt_minus2 = var_ym2(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym3,var_xm3] = mom(gx,hx,varshock,-3); 
corr_yFt_hIt_minus3 = var_ym3(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));
[var_ym4,var_xm4] = mom(gx,hx,varshock,-4); 
corr_yFt_hIt_minus4 = var_ym4(1,5)/((var_y(1,1)^(1/2)*var_y(5,5)^(1/2)));

MOMS_hp_model_est = ...
 [table_modelmoms(1,1);                      %  std(yF)*100       %1
  table_modelmoms(2,1)/table_modelmoms(1,1); %  std(cF)/std(yF)   %2
  table_modelmoms(3,1)/table_modelmoms(1,1); %  std(ivF)/std(yF); %3
  table_modelmoms(7,1)/table_modelmoms(1,1); %  std(tby)/std(yF); %4
  table_modelmoms(1,2);                      %  corr(yF_t,yF_t-1);%5
  table_modelmoms(2,3);                      %  corr(yF,cF);      %6  
  table_modelmoms(3,3);                      %  corr(yF,ivF);     %7
  table_modelmoms(7,3);                      %  corr(yF,tby);     %8
  table_modelmoms(6,1)/table_modelmoms(1,1); %  std(h)/std(yF);   %9
  table_modelmoms(6,3);                      %  corr(yF,h);       %10      
  table_modelmoms(4,1)/table_modelmoms(1,1); %  std(hF)/std(yF);  %11
  table_modelmoms(5,1)/table_modelmoms(1,1); %  std(hI)/std(yF);  %12
  table_modelmoms(4,3);                      %  corr(yF,hF);      %13        
  table_modelmoms(5,3);                      %  corr(yF,hI);      %14
  corr_yFt_hIt_minus4;                       %  corr(yF_t,hI_t-4);%15
  corr_yFt_hIt_minus3;                       %  corr(yF_t,hI_t-3);%16
  corr_yFt_hIt_minus2;                       %  corr(yF_t,hI_t-2);%17
  corr_yFt_hIt_minus1;                       %  corr(yF_t,hI_t-1);%18
  corr_yFt_hIt_plus1;                        %  corr(yF_t,hI_t+1);%19
  corr_yFt_hIt_plus2;                        %  corr(yF_t,hI_t+2);%20
  corr_yFt_hIt_plus3;                        %  corr(yF_t,hI_t+3);%21
  corr_yFt_hIt_plus4];                       %  corr(yF_t,hI_t+4);%22

MOMS_comparison_2moms = [MOMS_hp_data_all(1:14,1),MOMS_hp_model_est(1:14,1)]; % Final results reported in Table 3 / Column 3

MOMS_comparison_serial = [MOMS_hp_data_corrs_yF_hI,...
                          [corr_yFt_hIt_minus4;...
                           corr_yFt_hIt_minus3;...
                           corr_yFt_hIt_minus2;...
                           corr_yFt_hIt_minus1;...
                           corr_yFt_hIt;...
                           corr_yFt_hIt_plus1;...
                           corr_yFt_hIt_plus2;...
                           corr_yFt_hIt_plus3;...
                           corr_yFt_hIt_plus4;...
                           ]]; 
                       % Final results reported in Figure 8. In workspace we
                       % call this as serial_hI_AG
time_serial = [-4,-3,-2,-1,0,1,2,3,4]';
figure(1)
plot(time_serial,MOMS_comparison_serial(:,1),'--r',time_serial,MOMS_comparison_serial(:,2),'b',...
     time_serial,MOMS_comparison_serial(:,3),'--r',time_serial,MOMS_comparison_serial(:,4))

