% This wrapper file calls either the SMM estimation procedure or a single function solution call
% given a set of parameter values. 

clc
clear all

%% Set up the parallel environment
npool=7;
p = gcp('nocreate'); 
if isempty(p)
    poolsize = 0;
else
    poolsize = p.NumWorkers;
end
if poolsize <npool
    delete(gcp('nocreate'))
    parpool('local',npool)
end
warning('off','MATLAB:mir_warning_maybe_uninitialized_temporary');

%% Set if code runs for estimation (estimate=1) or for single function call (estimate=0)
estimate=0;

%% Set if saving data at each run
savedata= 0;    % Choose 1 to save the entire output for the function solution

%% Set if run IR as part of the solution
irflag= 1; % Choose 1 to run the IR which is used for the calculation of elasticities (and save the results), and 0 otherwise

%% Set non-estimated parameters according to MRS estimation and parametrization
PSI(1) = 0.211;   % psi_l1
PSI(2) = 0.162;   % psi_l2
PSI(3) = 0.535;   % rho_l
PSI(4) = 0.107;   % psi_t1
PSI(5) = 0.411;   % psi_t2
PSI(6) = -0.288;  % rho_t
PSI(7) = 0.903;   % eta_cp

chi = [2.2 2.2];    % level parameter in the tax function
taumu=[0.1 0.1];    % slope parameter in the tax function
tau_b=[12*(526+476)/4160 12*367/4160];    % minimum pre-tax transfers that the household gets regardless of work. Apply the scaling by 4160

%% If chose to estimate  
if estimate==1   
%% Theta 0
THETA_0(1) = 0.000133034;          % mu_l1_k
THETA_0(2) = 0.000461992;        % mu_l1_nk
THETA_0(3) = 0.000083;          % mu_l2_k
THETA_0(4) = 0.000148335;            % mu_l2_nk
THETA_0(5) = 0.000000000011407;     % mu_t1
THETA_0(6) = 0.0078413;            % mu_t2
THETA_0(7) = 1.1412;              % mu_0_k (nk normalized to 6)
THETA_0(8) = 4790/10^9;         % tau_gamma_nk/10^9
THETA_0(9) = 1.4763;              % sd_eps_k
THETA_0(10) = 0.6987;              % sd_eps_nk

lb = [THETA_0(1)/10,...
      THETA_0(2)/10,...
      THETA_0(3)/10,...
      THETA_0(4)/10,...
      THETA_0(5)/1000,...
      THETA_0(6)/10,...
      1.01,...
      100/10^9,...
      0.1,...
      0.1,...
      ];
ub = [THETA_0(1)*10,...
      THETA_0(2)*10,...
      THETA_0(3)*10,...
      THETA_0(4)*10,...
      THETA_0(5)+0.0002,...
      THETA_0(6)*10,...
      1.5,...
      12000/10^9,...
      4,...
      4,...
      ];

%% call estimation with Simulated Annealing   
fun = @(THETA)main_estimation(THETA,PSI,chi,taumu,tau_b,savedata,irflag);
options = saoptimset('PlotFcns',{@saplotbestf,@saplottemperature,@saplotf,@saplotstopping,@saplotbestx},'HybridFcn',{@fmincon},'HybridInterval',100,'InitialTemperature',10^4,'TemperatureFcn',@mytemp);
[THETA,fval,exitFlag,output] = simulannealbnd(fun,THETA_0,lb,ub,options);
save estimation_final_output

else
%% call single run with THETA_0
% load theta from file
THETA_0=csvread('THETA.csv');   % Can plug here the theta for the baseline solution
main_estimation(THETA_0,PSI,chi,taumu,tau_b,savedata,irflag);
end
