% THIS VERSION: A first version at trying to replicate LWZ. It calls for
% the mod file. It first sets the steady state values, taken from the
% paper. The parameters are taken from the paper. 





%% STEP 1. CALIBRATION


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I take the parameters from Table 1 and Table 2, Posterior Mode. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Here are the parameters to calibrate: 
% ggammap Omegahp Omegaep lambdaabarp lambdaqbarp gammahp gammaep betap 
% Omegap lambdakp mubp muep thetabarp thetap alphap phip deltap Chp Yp Cep 
% Ip Lhp Lep qlp Bp Rp Lbarp rhozp rhonuzp rhoqp rhonuqp rhozp rhovarphip rhopsip rhothetap  

% Group I 

gammahp=0.4655;                             % Habit persistence of households 
gammaep=0.6050;                             % Habit persistence of entrepreneurs 
Omegap=0.1881;                              % Capital accumulation frictions
ggammap=1+0.3682/100;                       % The aggregate growth factor of the trend Big Gamma
lambdaqbarp=1+0.12530/100;                    % The growth factor of Q 



% Group II 

betap=0.9870;                               % The discount factor 
lambdaabarp=0.0068;                         % The growth rate of A - a preference shock of Hhs
varphibarp=0.0497;                          % The housing preference mean in Hhs problem
phip=0.0697;                                % The weight of land in the production function 
deltap=0.0369;                              % The depreciation rate 
alphap=0.3;                                 % Calibrated non-labor income share in production function
thetabarp=0.75;                             % The credit thightness mean
psibarp=0.5; %missing !! determiens hours worked, from the preference for leisure but doesn't affect the dynamics





% Group 3 (shocks: persistence and variance) 

rhoap=0.9108;                                
rhozp=0.4743;
rhonuzp=0.0074;
rhoqp=0.6078;
rhonuqp=0.2920;
rhovarphip=0.9998;
rhopsip=0.9799;
rhothetap=0.9790;

sigmaap=0.1387;
sigmazp=0.0036;
sigmanuzp=0.0038;
sigmaqp=0.0037;
sigmanuqp=0.0025;
sigmavarphip=0.0543;
sigmapsip=0.0073;
sigmathetap=0.0126;


% NEW PARAMETERS 

% betacp lambdacbarp gammacp  phicp  rhocp
betacp=0.90; %discount factor for constrained Hhs
lambdacbarp=lambdaabarp;  % The growth rate of A
gammacp=gammahp; % habit persistence
phicp=0.5; %thetabarp; % credit tightness, for now the same as entrepreneurs
rhocp=rhoap; 




%some are defined in text as function of the above params: 

lambdakp=ggammap*lambdaqbarp;
Omegahp=(ggammap-betap*(1+lambdaabarp)*gammahp)*(ggammap-gammahp);
Omegaep=(ggammap-betap*gammaep)*(ggammap-gammaep);

%% STEP 2. STEADY STATES 

R_ss=ggammap/(betap*(1+lambdaabarp)); 
mub_mue_ss=betap*lambdaabarp/ggammap;
ql_Le_Y_ss=betap*alphap*phip/(1-betap-betap*lambdaabarp*thetabarp);
I_K_ss=1-((1-deltap)/lambdakp); 
K_Y_ss=betap*alphap*(1-phip)/(1-(betap/lambdakp)*(lambdaabarp*thetabarp+1-deltap));
Y_K_ss=1/K_Y_ss;
I_Y_ss=I_K_ss*K_Y_ss; 
B_Y_ss=thetabarp*ggammap*ql_Le_Y_ss+(thetabarp/lambdaqbarp)*K_Y_ss;
Ce_Y_ss=alphap-I_Y_ss-(B_Y_ss)*(1-betap*(1+lambdaabarp))/(ggammap);

mucb_muc_ss=(betap-betacp)*(1+lambdaabarp)/ggammap;


Ch_Y_ss=(1-Ce_Y_ss-I_Y_ss)*(ggammap-betap*gammahp*(1+lambdaabarp))/(2*ggammap-(1+lambdaabarp)*(betap+betacp));
Cc_Y_ss=(1-Ce_Y_ss-I_Y_ss)*(ggammap-betacp*gammahp*(1+lambdacbarp))/(2*ggammap-(1+lambdacbarp)*(betap+betacp));

ql_Lh_Ch_ss=(varphibarp*(ggammap-gammahp))/(ggammap*(1-ggammap/R_ss)*(1-gammahp/R_ss));

ql_Lh_Y_ss=ql_Lh_Ch_ss*Ch_Y_ss;

ql_Lc_Cc_ss=(varphibarp*(ggammap-gammahp))/((ggammap-betap*gammacp*(1+lambdacbarp))*(1-betacp*phicp*(betap-betacp)*(1+lambdaabarp)-betacp*(1+lambdacbarp)));

ql_Lc_Y_ss=ql_Lc_Cc_ss*Cc_Y_ss;

Lc_Le_ss=1/ql_Le_Y_ss*Cc_Y_ss*ql_Lc_Cc_ss;
Lh_Le_ss=1/ql_Le_Y_ss*Ch_Y_ss*ql_Lh_Ch_ss;
Le_Lbar_ss=1/(1+Lc_Le_ss+Lh_Le_ss);
Lc_Lbar_ss=Le_Lbar_ss*Lc_Le_ss; 
Lh_Lbar_ss=Le_Lbar_ss*Lh_Le_ss;

Lbar_lc_ss=1/Lc_Lbar_ss;
Lbar_lh_ss=1/Lh_Lbar_ss;
Lbar_le_ss=1/Le_Lbar_ss;

N_ss=(1/Ch_Y_ss)*((1-alphap)*ggammap*(1-gammahp/R_ss))/(psibarp*(ggammap-gammahp));

Bc_Y_ss=phicp*ql_Lc_Y_ss;
w_Nc_Y_ss=Cc_Y_ss+Bc_Y_ss*(1/R_ss-ggammap);

ql_Le_B_ss=ql_Le_Y_ss/B_Y_ss;
S_Y_ss=Bc_Y_ss+B_Y_ss;
Bc_S_ss=(1-B_Y_ss*1/S_Y_ss);
B_S_ss=1-Bc_S_ss;

S_Y_ss=B_Y_ss+Bc_Y_ss;
w_Nh_Y_ss=Ch_Y_ss+S_Y_ss*(1/R_ss-1/ggammap);

%% Save Params  




save LWZ_params.mat  gammahp gammaep Omegap ggammap lambdaqbarp betap lambdaabarp varphibarp phip deltap alphap thetabarp rhoap rhozp ...; 
     rhonuzp rhoqp rhonuqp rhovarphip rhopsip rhothetap sigmaap sigmazp sigmanuzp sigmaqp sigmanuqp sigmavarphip sigmapsip sigmathetap lambdakp ...;
     Omegahp Omegaep betacp lambdacbarp gammacp phicp rhocp;
 


%% Save Steady State 

save LWZ_ss.mat R_ss mub_mue_ss ql_Le_Y_ss I_K_ss K_Y_ss Y_K_ss I_Y_ss B_Y_ss Ce_Y_ss Ch_Y_ss ql_Lh_Ch_ss Lh_Le_ss ... ; 
    N_ss Le_Lbar_ss Lh_Lbar_ss ql_Le_B_ss mucb_muc_ss Cc_Y_ss ql_Lc_Cc_ss ql_Lc_Y_ss Le_Lbar_ss Lh_Lbar_ss Lc_Lbar_ss ...;
    Bc_Y_ss B_S_ss w_Nc_Y_ss Bc_S_ss S_Y_ss ql_Lh_Y_ss w_Nh_Y_ss   ;





% ggammap Omegahp Omegaep lambdaabarp lambdaqbarp gammahp gammaep betap Omegap lambdakp thetabarp alphap phip deltap  rhozp rhoap rhonuzp 
% rhoqp rhonuqp  rhovarphip rhopsip rhothetap R_ss mub_mue_ss ql_Le_Y_ss I_K_ss K_Y_ss Y_K_ss I_Y_ss B_Y_ss Ce_Y_ss Ch_Y_ss ql_Lh_Ch_ss Lh_Le_ss N_ss 
% Lbar_Le_ss Le_Lbar_ss Lh_Lbar_ss ql_Le_B_ss; 


%% Run Dynare
dynare liu_wang_zha_modif1.mod