clear all; clc
% This script estimates the discrete planned retirement age model 
%mex mexCombi.cpp COMPFLAGS="/openmp $COMPFLAGS"
LOAD = 0;
SAVE = 1;
do_no_cov = 0;

NumStart = 4;
do_particleswarm = 1;

str_cov = '';
if do_no_cov==1
    str_cov = 'no_cov';
end

diary off
if SAVE==1
    diary(['logs\' char(datetime('now','Format','yyyy_MM_dd''T''HHmmss')) str_cov '.txt'])
end
tic

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Specify which parameters to estimate (in excess of the explanatory variables)
if do_no_cov==1
    EstPar  = {'joint_h','SPA_w','beta_h','beta_w',...
           'RetAge_cons','RetAge_trend','RetAge55','RetAge60','RetAge65','var_w'};
else
EstPar  = {'joint_h','SPA_w','beta_h','beta_w',...
           'RetAge_cons','RetAge_trend','RetAge55','RetAge60','RetAge65','var_w','cov'};
end
ExplVar.h.Own = {'Highskilled','PoorHealth','ExpWorseHealth','Income','pensPPP','pensEPS','cohort'}; % these are all individual-specific
ExplVar.w.Own = {'Highskilled','PoorHealth','ExpWorseHealth','Income','pensPPP','pensEPS','cohort'}; % these are all individual-specific
ExplVar.h.Spouse = {'Income','pensPPP','pensEPS'}; 
ExplVar.w.Spouse = {'Income','pensPPP','pensEPS'};

% exclusively in the Auxiliary model:
ExplVar_excl.w = {'SPA_5155','SPA_56'};  % these variables are only entering the wife's utility. Nkids: this is a simple way to enusre that in the auxiliary models household-specific variables enter only once!
ExplVar_excl.h = {}; % these variables are only entering the wife's utility. Nkids: this is a simple way to enusre that in the auxiliary models household-specific variables enter only once!

%%%%%%%%%%%%%%
% 2. Load data
seed = 9210; % seed used to draw the normal errors for both spouses
S    = 2000;  % number of simulated shocks - used to evaluate the integral
% we use only 100 in the initial multi-start run and then switch to this in
% a final estimation
data = fun.load_data(seed,S);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Specify the explanatory variables (in excess of age dummies and, thus, no constant)
data = fun.explanatory_vars(data,ExplVar,ExplVar_excl,EstPar);
 % replace expectations with NaN if one of the explanatory variables are nan.
remove = any(isnan([data.h.X data.w.X data.w.SPA])')';
data.w.ExpRetAge(remove) = NaN;
data.h.ExpRetAge(remove) = NaN;

%%%%%%%%%%%%%%%%%%%%%
% 4. Setup parameters
par = fun.setup(data);
if do_no_cov==1
    par.cov = 0;
end
par.do_particleswarm = do_particleswarm;
par.NumStart = NumStart; % reduce to get some results
par.names    = data.names;

par.EstPar  = EstPar;
par.ExplVar = ExplVar;
par.ExplVar_excl = ExplVar_excl;

if LOAD==0
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 5. Calculate moments and auxiliary parameters and weight matrix using data
    b_aux                    = fun.aux_est_data(data,par);
    [mom,mom_var,names_mom]  = fun.aux_moments(b_aux,data,par,'data'); % N-by-J vector, where the age-dummies (cuf-offs in the oprobit) are collapsed into one element related to the relevant age of the individual

    [bsCov, bsDiagVar] = fun.bootstrap_cov(b_aux,data,par);
    W         = inv(bsDiagVar);
    par.CoVar = bsCov*data.Nobs; % multiply by N because the bootstrapped covariance is based on the sample-mean of the moments
    par.W     = W;
    par.b_aux = b_aux;
    par.names_mom = names_mom;
    par.mom   = mom;
    par.S     = S;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 6. Estimate parameters using indirect inference

    par = fun.estimate(EstPar,data,par,W,b_aux);
    disp DONE!
else
    load(['Results\est' str_cov]);
end

toc
if SAVE==1
    save(['Results\est' str_cov],'par')
end

[par.SE ,par.grad,par.sens,par.Avar] = fun.sensitivity(par.est,par.EstPar,data,par,par.W,par.b_aux);

[par.est' par.SE]
diary off

if SAVE==1
save(['Results\est' str_cov],'par')
end
