% I. This code generates the steady state for the baseline model

function flag = MinWageModel
x0_wH = 6;
finalWh = fsolve(@RunMinWageModel,x0_wH);
w_H_ss = finalWh;
[test, ~,vecSS] = RunMinWageModel(finalWh);
flag = 'done';
 
function [test, par, vecSS] = RunMinWageModel(GuessWh)

%set key parameters values

alpham_H   = 0.1875;
alpham_L   = 0.2297;
alpha      = 0.5;
vartheta   = 0.5;
varrho     = 1;   
chi_bar_H  = 1;
chi_bar_L  = 1;
iota       = 1;
phi        = 1;
eps_p      = 7;
eps_s      = 3;
eta        = 1;
theta      = 0.2;
psi        = 2;
mu         = 0.9388;
gamma      = 0.615;

rho_H_ss   = 0.984;  %0.984
rho_L_ss   = 0.967;  %0.967
rho_rh     = 0.8; %0.8437
sigma_rh   = 0.0008;  
rho_rl     = 0.3389;
sigma_rl   = 0.0028;
rho_b      = 0.9448;
sigma_b    = 0.0038;
rho_z      = 0.9752;
sigma_z    = 0.0031;
rho_i      = 0.8527;
sigma_i    = 0.0013;

% make a guess of high-sector wage
w_H_ss = GuessWh; 

% SOLVE FOR STEADY STATE 

z_ss=1;
pi_ss=1;
w_L=0.3;            % low sector wage is exogeneous
tau_u_ss=0.5*w_L; 
varsigma_H=1;
varsigma_L=1;

beta=0.96;
i_bar=0.003;
Lambda=beta;

% job filling/finding probabilities:
fs_H_ss=0.1718;
fs_L_ss=0.1608;
fv_H_ss=alpham_H*(fs_H_ss/alpham_H)^(-alpha/(1-alpha));
fv_L_ss=alpham_L*(fs_L_ss/alpham_L)^(-alpha/(1-alpha));

% Labor force participation =0.6465 which is the average over
% 2001-2019 period; I take labor participation of low-skill to correspond to their share in population ;

l_ss=0.645;
l_L_ss=0.06*0.6465;
l_H_ss=0.6465-l_L_ss;
e_H_ss=(fs_H_ss/(1-(1-fs_H_ss)*rho_H_ss))*l_H_ss;
e_L_ss=(fs_L_ss/(1-(1-fs_L_ss)*rho_L_ss))*l_L_ss;
e_ss=e_H_ss+e_L_ss;
u_ss=l_ss-e_ss;
u_H_ss=l_H_ss-e_H_ss;
u_L_ss=l_L_ss-e_L_ss;

% searchers and vacancies equal:
s_H_ss=l_H_ss-rho_H_ss*e_H_ss;
s_L_ss=l_L_ss-rho_L_ss*e_L_ss;
v_H_ss=s_H_ss*(fs_H_ss/fv_H_ss);
v_L_ss=s_L_ss*(fs_L_ss/fv_L_ss);


tau_ss=0.2; 

% Consumption part:

% start with a guess for b_bar and then compute it so that c_u/c_eH=0.72
% x0_bbar = 0.36; 
% [vecbbar,fval, exitflag] = fsolve(@(j)Findbbar(j,beta,i_bar, w_L, tau_u_ss, u_ss, l_ss, e_H_ss, e_L_ss, tau_ss),x0_bbar);
% b_bar_ss=vecbbar;

b_bar_ss=0.4;

x_ss=b_bar_ss;
c_n_ss=x_ss;
c_u_ss=x_ss+tau_u_ss;
c_eL_ss=x_ss+(1-tau_ss)*w_L; 
c_eH_ss=(1-beta*(1+i_bar)*e_H_ss)/beta*(1+i_bar)*((e_L_ss/c_eL_ss)+(u_ss/c_u_ss)+((1-l_ss)/c_n_ss));

chi_H_ss=chi_bar_H*(e_H_ss^(1+iota))/(1+iota); % disutility cost of labor
chi_L_ss=chi_bar_L*(e_L_ss^(1+iota))/(1+iota);
u_cu=log(c_u_ss);
u_cn=log(c_n_ss);
u_ceL=log(c_eL_ss);
u_ceH=log(c_eH_ss);

% u_cu=(c_u_ss^(1-varrho))/(1-varrho);
% u_cn=(c_n_ss^(1-varrho))/(1-varrho);
% u_ceL=(c_eL_ss^(1-varrho))/(1-varrho);
% u_ceH=(c_eH_ss^(1-varrho))/(1-varrho);

GammaH_ss=(c_eH_ss-x_ss)+(u_cu-(u_ceH-(2+iota)*chi_H_ss))*c_eH_ss;
GammaL_ss=(c_eL_ss-x_ss)+(u_cu-(u_ceL-(2+iota)*chi_L_ss))*c_eH_ss;

% from participation condition:
xi_bar_H=((fs_H_ss/(1-beta*rho_H_ss*(1-fs_H_ss)))*((1-tau_ss)*w_H_ss-GammaH_ss)+c_eH_ss*(log(c_u_ss)-log(c_n_ss)))/c_eH_ss*(1+log(mu-l_H_ss));
xi_bar_L=((fs_L_ss/(1-beta*rho_L_ss*(1-fs_L_ss)))*((1-tau_ss)*w_L-GammaL_ss)+c_eH_ss*(log(c_u_ss)-log(c_n_ss)))/c_eH_ss*(1+log(1-mu-l_L_ss));

xi_H_ss = xi_bar_H*log(mu-l_H_ss); % utility from home goods
xi_L_ss = xi_bar_L*log(1-mu-l_L_ss);
MRS_NC_Hss=xi_bar_H*c_eH_ss;
MRS_NC_Lss=xi_bar_L*c_eH_ss;

% xi_H_ss = xi_bar_H*((mu-l_H_ss)^(1-phi))/(1-phi);
% xi_L_ss = xi_bar_L*((1-mu-l_L_ss)^(1-phi))/(1-phi);
% MRS_NC_Hss=xi_bar_H*((mu-l_H_ss)^(1-phi))*(c_eH_ss)^(varrho);
% MRS_NC_Lss=xi_bar_L*((1-mu-l_L_ss)^(1-phi))*(c_eH_ss)^(varrho);

OmegaH_ss=(u_cu-(u_cn+xi_H_ss))*c_eH_ss;
OmegaL_ss=(u_cu-(u_cn+xi_L_ss))*c_eH_ss;

% prices  

revper=0.01;
q_H_ss=w_H_ss*fv_H_ss*(1-rho_H_ss)/(fv_H_ss*(1-rho_H_ss)-(1-beta*rho_H_ss)*revper);
q_L_ss=w_L*fv_L_ss*(1-rho_L_ss)/(fv_L_ss*(1-rho_L_ss)-(1-beta*rho_L_ss)*revper);

% B_x is calibrated so that hiring costs take 1% of revenues;
B_H = revper*q_H_ss/(1-rho_H_ss)*(fs_H_ss^eta);
B_L = revper*q_L_ss/(1-rho_L_ss)*(fs_L_ss^eta);
k_H_ss =z_ss*B_H*(fs_H_ss)^eta; 
k_L_ss =z_ss*B_L*(fs_L_ss)^eta;

PstarH_P_ss=q_H_ss*eps_p/(eps_p-1);
PstarL_P_ss=q_L_ss*eps_p/(eps_p-1);
PH_P_ss=PstarH_P_ss;
PL_P_ss=PstarL_P_ss;

Y_H_ss=z_ss*e_H_ss/varsigma_H;  
Y_L_ss=z_ss*e_L_ss/varsigma_L;

Y_ss=(1/gamma)*Y_H_ss*PH_P_ss^eps_s;

%Y_ss=(gamma^(1/eps_s)*Y_H_ss^((eps_s-1)/eps_s)+(1-gamma)^(1/eps_s)*Y_L_ss^((eps_s-1)/eps_s))^(eps_s/(eps_s-1));

varpsi_H_ss=(PH_P_ss^(eps_p-eps_s))*Y_ss;
varpsi_L_ss=(PL_P_ss^(eps_p-eps_s))*Y_ss;

p_aH_ss=(eps_p/(eps_p-1))*q_H_ss*varpsi_H_ss/(1-beta*(1-theta)) ;
p_aL_ss=(eps_p/(eps_p-1))*q_L_ss*varpsi_L_ss/(1-beta*(1-theta)) ;
p_bH_ss= varpsi_H_ss/(1-beta*(1-theta));
p_bL_ss= varpsi_L_ss/(1-beta*(1-theta));

dw_H_ss=q_H_ss*z_ss*e_H_ss-w_H_ss*e_H_ss-k_H_ss*v_H_ss;
dw_L_ss=q_L_ss*z_ss*e_L_ss-w_L*e_L_ss-k_L_ss*v_L_ss;
d_H_ss=Y_ss*(PH_P_ss^(1-eps_s))-(1/gamma)*q_H_ss*z_ss*e_H_ss+dw_H_ss;
d_L_ss=Y_ss*(PL_P_ss^(1-eps_s))-(1/(1-gamma))*q_L_ss*z_ss*e_L_ss+dw_L_ss;
D_ss=d_H_ss+d_L_ss;

a_ss=b_bar_ss;
b_ss=b_bar_ss;


test=w_H_ss-vartheta*(q_H_ss*z_ss+beta*rho_H_ss*fs_H_ss*k_H_ss/fv_H_ss)-(1-vartheta)*GammaH_ss/(1-tau_ss);

par = [beta alpham_H alpham_L alpha vartheta chi_bar_H chi_bar_L xi_bar_H xi_bar_L iota eta eps_p eps_s theta ...
                   gamma psi mu rho_H_ss rho_L_ss z_ss i_bar B_H B_L b_bar_ss rho_z sigma_z rho_rh sigma_rh rho_rl sigma_rl rho_b sigma_b rho_i sigma_i ];

vecSS = [l_ss e_H_ss e_L_ss u_H_ss u_L_ss u_ss e_ss l_H_ss l_L_ss s_H_ss s_L_ss k_H_ss k_L_ss v_H_ss v_L_ss ...
         x_ss c_eH_ss c_eL_ss c_u_ss c_n_ss OmegaH_ss OmegaL_ss GammaL_ss GammaH_ss  w_H_ss w_L tau_u_ss tau_ss...
         fs_H_ss fs_L_ss fv_H_ss fv_L_ss PstarH_P_ss PstarL_P_ss PH_P_ss PL_P_ss p_aH_ss p_bH_ss varpsi_H_ss p_aL_ss p_bL_ss varpsi_L_ss q_H_ss q_L_ss Y_ss Y_H_ss Y_L_ss D_ss d_H_ss d_L_ss dw_H_ss dw_L_ss]; 
          
save param beta alpham_H alpham_L alpha vartheta chi_bar_H chi_bar_L xi_bar_H xi_bar_L iota eta eps_p eps_s theta ...
                   gamma psi mu rho_H_ss rho_L_ss z_ss i_bar B_H B_L b_bar_ss rho_z sigma_z rho_rh sigma_rh rho_rl sigma_rl rho_b sigma_b rho_i sigma_i
          
save vecSS l_ss e_H_ss e_L_ss u_H_ss u_L_ss u_ss e_ss l_H_ss l_L_ss s_H_ss s_L_ss k_H_ss k_L_ss v_H_ss v_L_ss ...
                  x_ss c_eH_ss c_eL_ss c_u_ss c_n_ss OmegaH_ss OmegaL_ss GammaL_ss GammaH_ss  w_H_ss w_L tau_u_ss tau_ss...
                  fs_H_ss fs_L_ss fv_H_ss fv_L_ss PstarH_P_ss PstarL_P_ss PH_P_ss PL_P_ss p_aH_ss p_bH_ss varpsi_H_ss p_aL_ss p_bL_ss varpsi_L_ss q_H_ss q_L_ss Y_ss Y_H_ss Y_L_ss D_ss d_H_ss d_L_ss dw_H_ss dw_L_ss 
          
          
% function y = Findbbar(j, beta, i_bar, w_L, tau_u_ss, u_ss, l_ss, e_H_ss, e_L_ss, tau_ss)
% y=0.72-(j+tau_u_ss)*beta*(1+i_bar)*((e_L_ss/(j+(1-tau_ss)*w_L))+(u_ss/(j+tau_u_ss))+((1-l_ss)/j))/(1-beta*(1+i_bar)*e_H_ss);

