dynare simpleSMM

set_param_value('std_u',1.0)
%% simulation -- figure out log-normal distribution 
aalpha = M_.params(1);
dis = M_.params(2);
beta = M_.params(3);
delta = M_.params(4);
mu = M_.params(5);
mean_u = M_.params(6);
gamma = M_.params(7);
xi = M_.params(8);
theta = M_.params(9);
psi = M_.params(10);
C1 = M_.params(11);
tauT = M_.params(12);
tauNT = M_.params(13);
pii = M_.params(14);
std_eT = M_.params(15);
std_eNT = M_.params(16);
std_u = M_.params(17);
kappa = M_.params(18);
std_i = M_.params(19);
phi = M_.params(20);
phipi = M_.params(21);
phiL = M_.params(22);
eta = M_.params(23);
q = M_.params(24);
p = M_.params(25);
mH = M_.params(26);
ibar = M_.params(27);
rho_c = M_.params(28);

fprintf("std_u = %.3f \n", std_u)


%% simulate  ergodic wealth distribution
tic 
% stage 2 is about calibration
Nsample2 = 100;
periods2 = 400;

% stage 1 is about ergodic

Nsample = 1 ;
periods = 1000;
pn = 10000;
perc = 0.01;
starting = 2;
effective_periods = periods-starting+2;
uhh_matrix = zeros(Nsample,effective_periods);
pitaylor_matrix = zeros(Nsample,effective_periods);
piT_matrix = zeros(Nsample,effective_periods);
piNTh0_matrix = zeros(Nsample,effective_periods);
piNTf0_matrix = zeros(Nsample,effective_periods);
moment_cal = zeros(Nsample,effective_periods);
mNTh_matrix = zeros(Nsample, effective_periods);
mNTf_matrix = zeros(Nsample, effective_periods);
dyTh_matrix = zeros(Nsample,effective_periods);
dyNTh_matrix = zeros(Nsample,effective_periods);
dybarh_matrix = zeros(Nsample,effective_periods);

dch_matrix = zeros(Nsample,effective_periods);
dcf_matrix = zeros(Nsample,effective_periods);

%dybarh_matrix = zeros(Nsample,periods);


dyTf_matrix = zeros(Nsample,effective_periods);
dyNTf_matrix = zeros(Nsample,effective_periods);
dybarf_matrix = zeros(Nsample,effective_periods);
deNT_matrix = zeros(Nsample,effective_periods);
% document moments.
correlation_y_uhh = zeros(Nsample,1);

correlation_pi_y = zeros(Nsample,1);
correlation_pi_u = zeros(Nsample,1);


corrh = zeros(Nsample,1);
corrf = zeros(Nsample,1);

beta_hash = zeros(Nsample,4);
betai_hash = zeros(Nsample,4);

paircorr_dc = zeros(Nsample,1);
paircorr_dy = zeros(Nsample,1);
paircorr_dret = zeros(Nsample,1);

correlation_pi_HML = zeros(Nsample,1);
correlation_u_HML = zeros(Nsample,1);
correlation_pi_income= zeros(Nsample,1);

correlation_pi_ineq = zeros(Nsample,1);
correlation_u_ineq = zeros(Nsample,1);

correlation_pi_GDPEmp = zeros(Nsample,1);



uni_reg_u = zeros(Nsample,2);
uni_reg_u_pval = zeros(Nsample,2);
uni_reg_u_se = zeros(Nsample,2);

uni_reg_prod = zeros(Nsample,2);
uni_reg_prod_pval = zeros(Nsample,2);
uni_reg_prod_se = zeros(Nsample,2);

uni_reg_prod_resid = zeros(Nsample,2);
uni_reg_prod_pval_resid = zeros(Nsample,2);
uni_reg_prod_se_resid = zeros(Nsample,2);

bi_reg_prod_income = zeros(Nsample,2);
bi_reg_prod_income_pval = zeros(Nsample,2);
bi_reg_prod_income_se = zeros(Nsample,2);

uni_reg_ineq = zeros(Nsample,2);
uni_reg_ineq_pval = zeros(Nsample,2);
uni_reg_ineq_se = zeros(Nsample,2);


bi_reg_u_ineq = zeros(Nsample,2);
bi_reg_u_ineq_pval = zeros(Nsample,2);
bi_reg_u_ineq_se = zeros(Nsample,2);


bi_reg_u_income = zeros(Nsample,2);
bi_reg_u_income_pval = zeros(Nsample,2);
bi_reg_u_income_se = zeros(Nsample,2);

bi_reg_u_GDPEmp = zeros(Nsample,2);
bi_reg_u_GDPEmp_pval = zeros(Nsample,2);
bi_reg_u_GDPEmp_se = zeros(Nsample,2);

bi_reg_ineq_GDPEmp = zeros(Nsample,2);
bi_reg_ineq_GDPEmp_pval = zeros(Nsample,2);
bi_reg_ineq_GDPEmp_se = zeros(Nsample,2);

bi_reg_income_ineq = zeros(Nsample,2);
bi_reg_income_ineq_pval = zeros(Nsample,2);
bi_reg_income_ineq_se = zeros(Nsample,2);

% HML regression 
inflation_reg_unemp_HML = zeros(Nsample,2);
inflation_reg_unemp_HML_pval = zeros(Nsample,2);
inflation_reg_unemp_HML_se = zeros(Nsample,2);

inflation_reg_unemp_HML_MKT = zeros(Nsample,3);
inflation_reg_unemp_HML_MKT_pval = zeros(Nsample,3);
inflation_reg_unemp_HML_MKT_se = zeros(Nsample,3);

inflation_reg_HML = zeros(Nsample,1);
inflation_reg_HML_pval = zeros(Nsample,1);
inflation_reg_HML_se = zeros(Nsample,1);

inflation_reg_unemp_HMLt0 = zeros(Nsample,2);
inflation_reg_unemp_HMLt0_pval = zeros(Nsample,2);
inflation_reg_unemp_HMLt0_se = zeros(Nsample,2);


inflation_reg_unemp_HMLt0_MKT = zeros(Nsample,3);
inflation_reg_unemp_HMLt0_MKT_pval = zeros(Nsample,3);
inflation_reg_unemp_HMLt0_MKT_se = zeros(Nsample,3);

inflation_reg_HMLt0 = zeros(Nsample,1);
inflation_reg_HMLt0_pval = zeros(Nsample,1);
inflation_reg_HMLt0_se = zeros(Nsample,1);

% NMT regression 

inflation_reg_unemp_NMT = zeros(Nsample,2);
inflation_reg_unemp_NMT_pval = zeros(Nsample,2);
inflation_reg_unemp_NMT_se = zeros(Nsample,2);

inflation_reg_unemp_NMT_MKT = zeros(Nsample,3);
inflation_reg_unemp_NMT_MKT_pval = zeros(Nsample,3);
inflation_reg_unemp_NMT_MKT_se = zeros(Nsample,3);

inflation_reg_NMT = zeros(Nsample,1);
inflation_reg_NMT_pval = zeros(Nsample,1);
inflation_reg_NMT_se = zeros(Nsample,1);

inflation_reg_unemp_NMTt0 = zeros(Nsample,2);
inflation_reg_unemp_NMTt0_pval = zeros(Nsample,2);
inflation_reg_unemp_NMTt0_se = zeros(Nsample,2);

inflation_reg_unemp_NMTt0_MKT = zeros(Nsample,3);
inflation_reg_unemp_NMTt0_MKT_pval = zeros(Nsample,3);
inflation_reg_unemp_NMTt0_MKT_se = zeros(Nsample,3);

inflation_reg_NMTt0 = zeros(Nsample,1);
inflation_reg_NMTt0_pval = zeros(Nsample,1);
inflation_reg_NMTt0_se = zeros(Nsample,1);

tri_reg_u_income_ineq  = zeros(Nsample,3);
tri_reg_u_income_ineq_pval  = zeros(Nsample,3);
tri_reg_u_income_ineq_se = zeros(Nsample,3);

%-----------------erogidc------------------------
%-----------------erogidc------------------------
%-----------------erogidc------------------------
%-----------------erogidc------------------------
%-----------------erogidc------------------------
%-----------------erogidc------------------------
k=1;

    rng(k,'twister');
    simul_periods = periods;
    
    % generate shocks
    
    uh = normrnd(0,1,[simul_periods,1]);
    uf = normrnd(0,1,[simul_periods,1]);

    eTh = normrnd(0,1,[simul_periods,1]);
    eTf = normrnd(0,1,[simul_periods,1]);
    eNTh = normrnd(0,1,[simul_periods,1]);
    eNTf = normrnd(0,1,[simul_periods,1]);
    ei = normrnd(0,1,[simul_periods,1]);
%%
    M = [uh,uf,eTh,eTf,eNTh,eNTf];
    R = [1 0 0 0 0 0 ;...
        0 1 0 0 0 0; ...
        0,0,1,0 0 0; ...
        0,0,0,1, 0, 0; ...
        0, 0, 0,0, 1,0; ...
        0,0,0,0,0,1];

    L = chol(R);
    M = M*L;
    uh = M(:,1);
    uf = M(:,2);
    eTh = M(:,3);
    eTf = M(:,4);
    eNTh = M(:,5);
    eNTf = M(:,6);
    shocks = zeros(7,simul_periods); 
    shocks(1,1:end) = uh;
    shocks(2,1:end) = uf;
    shocks(3,1:end) = eTh;
    shocks(4,1:end) = eTf;
    shocks(5,1:end) = eNTh;
    shocks(6,1:end) = eNTf;
    shocks(7,1:end) = ei;

    %% simulate 
    initial_condition_states = repmat(oo_.dr.ys,1,M_.maximum_lag);
    options_.irf = simul_periods;
    shock_matrix = zeros(options_.irf,M_.exo_nbr); %create shock matrix with number of time periods in columns
    order = 3;
    r = simult_(M_,options_,initial_condition_states,oo_.dr,shocks',options_.order);
    periods = simul_periods;

    %% extract simulation series 
    % burn in here 
    

    uhh = (r(find(strcmp(M_.endo_names,'uhh')),starting:end));
    uff = (r(find(strcmp(M_.endo_names,'uff')),starting:end));


    uh_basic = exp(uh*std_u+mean_u);
    uf_basic = exp(uf*std_u+mean_u);

    uhh = rho_c*uf_basic + (1-rho_c)*uh_basic;
    uff = rho_c*uh_basic + (1-rho_c)*uf_basic;
    

    uhh5 = movsum(uhh,5);
    uff5 = movsum(uff,5);

    pdTh = r(find(strcmp(M_.endo_names,'pdTh')),starting:end);
    pdTf = r(find(strcmp(M_.endo_names,'pdTf')),starting:end);
    pdNTh = r(find(strcmp(M_.endo_names,'pdNTh')),starting:end);
    pdNTf = r(find(strcmp(M_.endo_names,'pdNTf')),starting:end);
    
    wbar = r(find(strcmp(M_.endo_names,'wbar')),starting:end);

    w = r(find(strcmp(M_.endo_names,'w')),starting:end);
    lambda = r(find(strcmp(M_.endo_names,'lambda')),starting:end);
    
    sTh = r(find(strcmp(M_.endo_names,'sTh')),starting:end);
    sNTh = r(find(strcmp(M_.endo_names,'sNTh')),starting:end);
    sTf = r(find(strcmp(M_.endo_names,'sTf')),starting:end);
    sNTf = r(find(strcmp(M_.endo_names,'sNTf')),starting:end);

    logyTr = r(find(strcmp(M_.endo_names,'logyTr')),starting:end);

    Lh = r(find(strcmp(M_.endo_names,'Lh')),starting:end);
    Lf = r(find(strcmp(M_.endo_names,'Lf')),starting:end);

    wHTh = exp(r(find(strcmp(M_.endo_names,'wHTh')),starting:end));
    wHNTh = exp(r(find(strcmp(M_.endo_names,'wHNTh')),starting:end));
    wHTf = exp(r(find(strcmp(M_.endo_names,'wHTf')),starting:end));
    wHNTf = exp(r(find(strcmp(M_.endo_names,'wHNTf')),starting:end));

    wSTh = exp(r(find(strcmp(M_.endo_names,'wSTh')),starting:end));
    wSNTh = exp(r(find(strcmp(M_.endo_names,'wSNTh')),starting:end));
    wSTf = exp(r(find(strcmp(M_.endo_names,'wSTf')),starting:end));
    wSNTf = exp(r(find(strcmp(M_.endo_names,'wSNTf')),starting:end));

    HS_ratio = (wHTh+wHNTh+wHTf+wHNTf)./((wHTh)+(wHNTh)+(wHTf)+(wHNTf) + (wSTh)+(wSNTh)+(wSTf)+(wSNTf));

    wh = exp(wbar)*1./(1+exp(w));
    wf = exp(wbar).*exp(w)./(1+exp(w));
    
    whprime = wh - (1-exp(-uhh))'*eta.*(1./(1+exp(logyTr)).*exp(pdTh) + (1-aalpha)/aalpha./(1+exp(lambda)).*exp(pdNTh)) - xi.*(wHTh+wHNTh); 
    wfprime = wf - (1-exp(-uff))'*eta.*(exp(logyTr)./(1+exp(logyTr)).*exp(pdTf) + (1-aalpha)/aalpha.*exp(lambda)./(1+exp(lambda)).*exp(pdNTf))- xi.*(wHTf+wHNTf); 
    
    Dprime = (1-(1-exp(-uhh)).*eta)'.*(1./(1+exp(logyTr)).*(1-sTh) +1./(1+exp(lambda)).*(1-sNTh)) + (1-(1-exp(-uff)).*eta)'.*(exp(logyTr)./(1+exp(logyTr)).*(1-sTf) +exp(lambda)./(1+exp(lambda)).*(1-sNTf));
    wageprime = (1-xi).*(1./(1+exp(logyTr)).*(sTh)+1./(1+exp(lambda)).*(sNTh)+1./(1+exp(-logyTr)).*(sTh)+1./(1+exp(-lambda)).*(sNTh)) ;

    incomeprime = Dprime+wageprime;
    Dprimeh = incomeprime.*whprime./(whprime+wfprime)./wh;
    Dprimef = incomeprime.*wfprime./(whprime+wfprime)./wf;
    
    bh = r(find(strcmp(M_.endo_names,'bh')),starting:end);
    bf = r(find(strcmp(M_.endo_names,'bf')),starting:end);
    Omegah = r(find(strcmp(M_.endo_names,'Omegah')),starting:end);
    Omegaf = r(find(strcmp(M_.endo_names,'Omegaf')),starting:end);

    leisure_exp_h = 1./(1+exp(lambda)).*(1-psi)/psi./wh;
    leisure_exp_f = exp(lambda)./(1+exp(lambda)).*(1-psi)/psi./wf;
    
    leisure_h = r(find(strcmp(M_.endo_names,'Nh')),starting:end);
    leisure_f = r(find(strcmp(M_.endo_names,'Nf')),starting:end);
    pitaylor = r(find(strcmp(M_.endo_names,'pitaylor')),starting:end);
    piT = r(find(strcmp(M_.endo_names,'piT')),starting:end);
    piNTf0 = r(find(strcmp(M_.endo_names,'piNTf0')),starting:end);
    piNTh0 = r(find(strcmp(M_.endo_names,'piNTh0')),starting:end);


    dch = r(find(strcmp(M_.endo_names,'dch')),starting:end);
    dcf = r(find(strcmp(M_.endo_names,'dcf')),starting:end);


    dyTh_matrix(k,:) = r(find(strcmp(M_.endo_names,'dyTh')),starting:end);
    dyNTh_matrix(k,:) = r(find(strcmp(M_.endo_names,'dyNTh')),starting:end);
    dybarh = aalpha*exp(r(find(strcmp(M_.endo_names,'dyTh')),starting:end))+(1-aalpha)*exp(r(find(strcmp(M_.endo_names,'dyNTh')),starting:end));
    dybarh = exp(r(find(strcmp(M_.endo_names,'dGDPh')),starting:end));


    dyTf_matrix(k,:) = r(find(strcmp(M_.endo_names,'dyTf')),starting:end);
    dyNTf_matrix(k,:) = r(find(strcmp(M_.endo_names,'dyNTf')),starting:end);
    dybarf = aalpha*exp(r(find(strcmp(M_.endo_names,'dyTf')),starting:end))+(1-aalpha)*exp(r(find(strcmp(M_.endo_names,'dyNTf')),starting:end));

    
    
    
    dch_matrix(k,:) = r(find(strcmp(M_.endo_names,'dch')),starting:end);
    dcf_matrix(k,:) = r(find(strcmp(M_.endo_names,'dcf')),starting:end);
    

    %%
    %response_noshock(k,:) = bh;

    %disp_h_path = exp(-abs(bh));
    disp_h_path = 1- (wSTh+wSNTh).*(1-exp(-uhh))'*eta./wh;
    %%
    %disp_f_path = exp(-abs(bf));
    disp_f_path = 1- (wSTf+wSNTf).*(1-exp(-uff))'*eta./wf;

    LaborLoss_H_path = (1-xi.*HS_ratio);
    LaborLoss_F_path = (1-xi.*HS_ratio);
    
    
    % remember this needs to be divided by wh, wf respectively, since we will
    % calculate in terms of wealth share.
    
    allo_h_path = Dprimeh;
    allo_f_path = Dprimef;
    
    % initial distribution
    h_initial_dis = ones(pn,1)/pn;
    f_initial_dis = ones(pn,1)/pn;

toc

%%

mean(uhh),mean(uff);
pdTh(1:5)+pdTf(1:5)
