 
%Chatterjee, Morley, and Singh (CMS) "Estimating Household Consumption Insurance",
%forthcoming, Journal of Applied Econometrics
%Working paper version available at https://ssrn.com/abstract=2933226
%
%If you use this code, please cite the paper
% 
%This version of the qmle_CMS_FD.m code: 11 November 2020
%Output file is qmle_CMS_M_FD.txt with no bias correction and qmle_CMS_bbc_M_FD.txt with bias correction 
%
% For a new dataset, it is a good idea to check for different starting
% values on line 61 to get global max. Modify assumptions for different shock
% variances over time for different datasets. Modify the state - space form 
% for different models than BPP. Set missing observations in dataset to 999.

clear, clc;
warning('off');

skip_estimation = 0; %set to 1 if estimation already done and just conducting estimation of biases given estimates and 0 otherwise
bootstrap = 0; %set to 1 if using semi-parametric bootstrap to bias-adjust estimates and 0 otherwise
bsim = 50; %set maximum number of bootstrap simulations to estimate biases
number_sim_start = 1; %update this after any crash with estimation for bootstrap simulations to the number_sim at which it crashed

% Select data code here. See below for descriptions.
data_code = 1;

if data_code == 1 % Select for whole sample
    N = 1765;
    T = 15;
    data_im = dlmread('cohA_CMS_added1978.txt');

elseif data_code == 2 % Select for low education
    N = 883;
    T = 15; 
    data_im = dlmread('cohA_low_CMS_added1978.txt');

elseif data_code == 3 % Select for high education
    N = 882;
    T = 15; 
    data_im = dlmread('cohA_high_CMS_added1978.txt');

elseif data_code == 4 % Select for age 1 (30-47 years)
    N = 1413;
    T = 15; 
    data_im = dlmread('cohA_age1_CMS_added1978.txt');

elseif data_code == 5 % Select for age 2 (48-65 years)
    N = 708;
    T = 15; 
    data_im = dlmread('cohA_age2_CMS_added1978.txt');
end

data_sample = [data_im(1:N*T, 1), data_im(1:N*T, 2)]; % Residual income is in column 1 and residual consumption is in column 2. 

start = 2; % Evaluate likelihood in second quarter of levels data, as in Morley, Nelson, and Zivot (2003).

y = data_sample;

%% MLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prmtr_in = -ones(36, 1); % Initial parameter values

if (skip_estimation==1)
estimates = dlmread('qmle_CMS_M_FD_store.txt');
prm_fnl = estimates(1:36,1);
sd_fnl = estimates(1:36,2);
xout = estimates(1:36,3); 
fout = estimates(1,4);   
    
else

% Optimization routine.
fprintf('Optimization in progress...\n');
options=optimoptions('fminunc','MaxFunctionEvaluations',15000,'MaxIterations',10000,'FiniteDifferenceType','central');
[xout,fout,cout,out,gout,hout] = fminunc(@(prmtr)lik_fcn_FD(prmtr,N,T,y,start),prmtr_in,options);
% Returns paramter estimates, -LL value, code, gradient, hessian

% Final parameter values
prm_fnl = trans(xout);

% Use Hessian to find parameter standard errors
hessn0 = -hout;
cov0_h = inv(-hessn0);

grdn_fnl = numeric_jacobian(@(xout)trans(xout),xout);

grad0 = numeric_jacobian(@(xout)scr_fcn_FD(xout,N,T,y),xout);
ggg=grad0'*grad0;
cov0_op=inv(ggg); 

ihess=inv(hessn0);    
cov0_sw=ihess*ggg*ihess;    

cov_h=grdn_fnl*cov0_h*grdn_fnl';
sd_fnl_h =sqrt(diag(cov_h)); %Standard errors of the estimated coefficients based on Hessian
    
cov_op=grdn_fnl*cov0_op*grdn_fnl';
sd_fnl_op =sqrt(diag(cov_op)); %Standard errors of the estimated coefficients based on outer product of gradients

cov_sw=grdn_fnl*cov0_sw*grdn_fnl';
sd_fnl =sqrt(diag(cov_sw)); %Standard errors of the estimated coefficients based on Huber-White sandwich estimator--use for QMLE@


store_results=[prm_fnl, sd_fnl, xout, fout.*ones(size(prm_fnl,1),1)];
dlmwrite('qmle_CMS_M_FD_store.txt',store_results);

end

results = [prm_fnl, sd_fnl];

% Creates output file to store results
output = fopen('qmle_CMS_M_FD.txt','w');

% Print table of results into 'qmle_CMS_M_FD.txt'
fprintf(output, 'Quasi log likelihood value is %1.8f\n', -fout);
fprintf(output, 'Parameter, QMLE parameter Estimates, Standard errors\n');
fprintf(output, 'theta        %0.8f              %0.8f\n', results(1, :));
fprintf(output, 'gamma_ctau   %0.8f              %0.8f\n', results(2, :));
fprintf(output, 'gamma_ceps   %0.8f              %0.8f\n', results(3, :));
fprintf(output, 'sig_u        %0.8f              %0.8f\n', results(4, :));
fprintf(output, 'sig_eta79-81 %0.8f              %0.8f\n', results(5, :));
fprintf(output, 'sig_eta82    %0.8f              %0.8f\n', results(6, :));
fprintf(output, 'sig_eta83    %0.8f              %0.8f\n', results(7, :));
fprintf(output, 'sig_eta84    %0.8f              %0.8f\n', results(8, :));
fprintf(output, 'sig_eta85    %0.8f              %0.8f\n', results(9, :)); 
fprintf(output, 'sig_eta86    %0.8f              %0.8f\n', results(10, :));
fprintf(output, 'sig_eta87    %0.8f              %0.8f\n', results(11, :));
fprintf(output, 'sig_eta88    %0.8f              %0.8f\n', results(12, :));
fprintf(output, 'sig_eta89    %0.8f              %0.8f\n', results(13, :));
fprintf(output, 'sig_eta90-92 %0.8f              %0.8f\n', results(14, :));
fprintf(output, 'sig_eps79    %0.8f              %0.8f\n', results(15, :));
fprintf(output, 'sig_eps80    %0.8f              %0.8f\n', results(16, :));
fprintf(output, 'sig_eps81    %0.8f              %0.8f\n', results(17, :));
fprintf(output, 'sig_eps82    %0.8f              %0.8f\n', results(18, :));
fprintf(output, 'sig_eps83    %0.8f              %0.8f\n', results(19, :));
fprintf(output, 'sig_eps84    %0.8f              %0.8f\n', results(20, :));
fprintf(output, 'sig_eps85    %0.8f              %0.8f\n', results(21, :));
fprintf(output, 'sig_eps86    %0.8f              %0.8f\n', results(22, :));
fprintf(output, 'sig_eps87    %0.8f              %0.8f\n', results(23, :));
fprintf(output, 'sig_eps88    %0.8f              %0.8f\n', results(24, :));
fprintf(output, 'sig_eps89    %0.8f              %0.8f\n', results(25, :));
fprintf(output, 'sig_eps90-92 %0.8f              %0.8f\n', results(26, :));
fprintf(output, 'sig_v79      %0.8f              %0.8f\n', results(27, :));
fprintf(output, 'sig_v80      %0.8f              %0.8f\n', results(28, :));
fprintf(output, 'sig_v81      %0.8f              %0.8f\n', results(29, :));
fprintf(output, 'sig_v82      %0.8f              %0.8f\n', results(30, :));
fprintf(output, 'sig_v83      %0.8f              %0.8f\n', results(31, :));
fprintf(output, 'sig_v84      %0.8f              %0.8f\n', results(32, :));
fprintf(output, 'sig_v85      %0.8f              %0.8f\n', results(33, :));
fprintf(output, 'sig_v86      %0.8f              %0.8f\n', results(34, :));
fprintf(output, 'sig_v89      %0.8f              %0.8f\n', results(35, :));
fprintf(output, 'sig_v90-92   %0.8f              %0.8f\n', results(36, :));
fclose(output);

if (bootstrap==1)
    
smoother = 1; %set to 1 if using Kalman smoother to estimate shocks for bootstrap and 0 for Kalman filter


%%Estimate empirical shocks  
[shocks_v,shocks_u,shocks_eps,shocks_eta] = shocks_FD(xout,N,T,y,smoother);

bias_store = []; 
prmtr_index_store = [];    
    
if (number_sim_start==1)
output = fopen('qmle_bias_bs_FD.txt','w');
fclose('all');
store = [];
elseif (number_sim_start>1)    
store = dlmread('qmle_bias_bs_FD.txt');
bias_store = store(1:number_sim_start-1,2:size(prm_fnl,1)+1);
prmtr_index_store = store(1:number_sim_start-1,size(prm_fnl,1)+2:2*size(prm_fnl,1)+1);    
end

number_sim = 1;
while number_sim <= bsim 
    
   if number_sim >= number_sim_start 
    fprintf('Simulation number...\n');
    disp(number_sim);
   end 
 
%% Draw sample from bootstrap DGP
y = dgp_FD(prm_fnl,N,T,shocks_v,shocks_u,shocks_eps,shocks_eta);     

col_iter = 1;
while col_iter <= size(y,2)
    
    row_iter = 1;
    while row_iter <= size(y,1)
        if (data_sample(row_iter,col_iter)==999)
            y(row_iter,col_iter) = 999; 
        end
    row_iter = row_iter + 1;
    end

col_iter = col_iter + 1;
end
           
if (number_sim<number_sim_start) 
    number_sim=number_sim+1; 
else 

%Estimation for bootstrap sample@
% Optimization routine.
fprintf('Optimization in progress...\n');
options=optimoptions('fminunc','MaxFunctionEvaluations',15000,'MaxIterations',10000,'FiniteDifferenceType','central');
[xout_bs,fout_bs,cout_bs,out_bs] = fminunc(@(prmtr)lik_fcn_FD(prmtr,N,T,y,start),prmtr_in,options);
% Returns paramter estimates, -LL value, code, gradient, hessian

% Final parameter values
prm_fnl_bs = trans(xout_bs);

bias_store=[bias_store;
           (prm_fnl_bs-prm_fnl)'];
    
prmtr_index = ones(size(prm_fnl,1),1);
dbb0 = ones(size(prm_fnl,1),1);
dbb1 = ones(size(prm_fnl,1),1);
    
if (size(bias_store,1)>2) 
bb0 = mean(bias_store)';
bb1 = mean(bias_store(1:size(bias_store,1)-1,1:size(prm_fnl,1)))';
bb2 = mean(bias_store(1:size(bias_store,1)-2,1:size(prm_fnl,1)))';  
dbb0 = abs(bb1-bb0);
dbb1 = abs(bb2-bb1);
end

if (sum(selif(prmtr_index,dbb0 < 0.001 & dbb1 < 0.001)) ~= size(prmtr_index,1))
prmtr_index(selif((1:size(prmtr_index,1))',dbb0 >= 0.001 | dbb1 >= 0.001)) = zeros(size(selif((1:size(prmtr_index,1))',dbb0 >= 0.001 | dbb1 >= 0.001),1),1);    
end

if (sum(prmtr_index)>0)
prmtr_index(selif((1:size(prmtr_index,1))',prmtr_index == 1 | prmtr_index_store(size(prmtr_index_store,1),:)' == 1)) = ones(size(selif((1:size(prmtr_index,1))',prmtr_index == 1 | prmtr_index_store(size(prmtr_index_store,1),:)' == 1),1),1);
end    

prmtr_index_store=[prmtr_index_store;
                   prmtr_index'];    

if (number_sim>1)    
store = dlmread('qmle_bias_bs_FD.txt');
end
bs_results=[store;
            number_sim,(prm_fnl_bs-prm_fnl)',prmtr_index'];               
output = fopen('qmle_bias_bs_FD.txt','w');
dlmwrite('qmle_bias_bs_FD.txt',bs_results);

      
if (sum(prmtr_index)==size(prm_fnl,1)) 
    number_sim = bsim; 
end

number_sim = number_sim + 1;
end

end

prm_fnl_bs = prm_fnl-mean(bias_store)'; %adjusts the estimates for the bootstrap estimate of bias
sd_fnl_bs = std(bias_store)';

results = [prm_fnl, prm_fnl_bs, sd_fnl, sd_fnl_bs, mean(bias_store)'];

% Creates output file to store results
output = fopen('qmle_CMS_bbc_M_FD.txt','w');

% Print table of results into 'qmle_CMS_M_FD.txt'
fprintf(output, 'Quasi log likelihood value is %1.8f\n', -fout);
fprintf(output, 'Parameter,    QMLE,     Bias-corrected QMLE,     SEs,     Bootstrap SEs,    Biases\n');
fprintf(output, 'theta        %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(1, :));
fprintf(output, 'gamma_ctau   %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(2, :));
fprintf(output, 'gamma_ceps   %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(3, :));
fprintf(output, 'sig_u        %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(4, :));
fprintf(output, 'sig_eta79-81 %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(5, :));
fprintf(output, 'sig_eta82    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(6, :));
fprintf(output, 'sig_eta83    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(7, :));
fprintf(output, 'sig_eta84    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(8, :));
fprintf(output, 'sig_eta85    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(9, :)); 
fprintf(output, 'sig_eta86    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(10, :));
fprintf(output, 'sig_eta87    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(11, :));
fprintf(output, 'sig_eta88    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(12, :));
fprintf(output, 'sig_eta89    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(13, :));
fprintf(output, 'sig_eta90-92 %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(14, :));
fprintf(output, 'sig_eps79    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(15, :));
fprintf(output, 'sig_eps80    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(16, :));
fprintf(output, 'sig_eps81    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(17, :));
fprintf(output, 'sig_eps82    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(18, :));
fprintf(output, 'sig_eps83    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(19, :));
fprintf(output, 'sig_eps84    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(20, :));
fprintf(output, 'sig_eps85    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(21, :));
fprintf(output, 'sig_eps86    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(22, :));
fprintf(output, 'sig_eps87    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(23, :));
fprintf(output, 'sig_eps88    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(24, :));
fprintf(output, 'sig_eps89    %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(25, :));
fprintf(output, 'sig_eps90-92 %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(26, :));
fprintf(output, 'sig_v79      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(27, :));
fprintf(output, 'sig_v80      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(28, :));
fprintf(output, 'sig_v81      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(29, :));
fprintf(output, 'sig_v82      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(30, :));
fprintf(output, 'sig_v83      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(31, :));
fprintf(output, 'sig_v84      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(32, :));
fprintf(output, 'sig_v85      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(33, :));
fprintf(output, 'sig_v86      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(34, :));
fprintf(output, 'sig_v89      %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(35, :));
fprintf(output, 'sig_v90-92   %0.8f     %0.8f   %0.8f   %0.8f   %0.8f\n', results(36, :));
fprintf(output, 'Number of bootstrap sims %1.8f\n', size(bias_store,1));
fclose(output);
end


