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


%parameters 

alphamH    = 1;    % matching efficiency
alphamL    = 1;
alpha      = 0.5;  % matching elasticity
vartheta   = 0.5;  % bargaining power
eps_p      = 6;    % within sector elasticity of substitution
eps_s      = 2;    % elasticity of substitution across sectors
theta      = 0.2;  % price rigidity
psi        = 0.5;    % Taylor rule, inflation coefficient
gamma      = 0.5;  % high sector's share in output

rho_H_ss   = 0.984; %0.984
rho_L_ss   = 0.967; %0.967
rho_rH     = 0.8437; %0.8437
sigma_rH   = 0.0008;  
rho_rL     = 0.3389; %0.3389
sigma_rL   = 0.0028; %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;
sigma_tauH  = 0.01;
sigma_tauL  = 0.01;

% SOLVE FOR STEADY STATE 

z_ss=1;
i_bar=0.003;
pi_H=1;
pi_L=1;
pi=1;

%prices
pHstar_P=1;
pLstar_P=1;
pHstar_PH=1;
pLstar_PL=1;

q_H_ss=(eps_p-1)/eps_p;
q_L_ss=(eps_p-1)/eps_p;
chi_bar_H   = (1/2)*q_H_ss;
chi_bar_L   = (1/2)*q_L_ss;

% Labor market 

u_H_ss=0.023; %0.023
u_L_ss=0.068; %0.068

e_H_ss=1-u_H_ss;
e_L_ss=1-u_L_ss;
s_H_ss=1-rho_H_ss*e_H_ss;
s_L_ss=1-rho_L_ss*e_L_ss;

fs_H_ss=(e_H_ss*(1-rho_H_ss))/(1-rho_H_ss*e_H_ss); 
fs_L_ss=(e_L_ss*(1-rho_L_ss))/(1-rho_L_ss*e_L_ss); 
tightnessH=(fs_H_ss/alphamH)^(1/(1-alpha));
tightnessL=(fs_L_ss/alphamL)^(1/(1-alpha));
v_H_ss=tightnessH*s_H_ss;
v_L_ss=tightnessL*s_L_ss;
fv_H_ss=alphamH*(tightnessH^(-alpha));  
fv_L_ss=alphamL*(tightnessL^(-alpha));


%discount factor
betaH=0.72/((1+i_bar)*(0.72*e_H_ss+u_H_ss));
betaL=0.72/((1+i_bar)*(0.72*e_L_ss+u_L_ss));
LambdaH=betaH;
LambdaL=betaL;

%output

Y_H_ss=e_H_ss;
Y_L_ss=e_L_ss;
Y_ss=(1/gamma)*Y_H_ss;
%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));


% defining the constant terms to simplify wage equation
constH=0.28+log(0.72)+chi_bar_H;
constL=0.28+log(0.72)+chi_bar_L;
const2H=(1-q_H_ss)*Y_ss+q_H_ss*e_H_ss;
const2L=(1-q_L_ss)*Y_ss+q_L_ss*e_L_ss;

x0_tau = 0.05;
x0_wH = 0.8; 
x0_wL = 0.8;
x0=[x0_tau, x0_wH, x0_wL];
[varw,fval,exitflag] = fsolve(@(x)Findvarw(x, betaH, betaL, vartheta, rho_H_ss, rho_L_ss, u_H_ss, u_L_ss, e_H_ss, e_L_ss, q_H_ss, q_L_ss, v_H_ss, v_L_ss, fv_H_ss, fv_L_ss, fs_H_ss, fs_L_ss, constH, constL, const2H, const2L),x0);
tau_ss = varw(1);
w_H_ss = varw(2);
w_L_ss = varw(3);

%unemployment benefits
tau_uH_ss=0.3*w_H_ss; 
tau_uL_ss=0.3*w_L_ss; 

% vacancy costs
k_H=(fv_H_ss*(q_H_ss-w_H_ss))/(1-betaH*rho_H_ss); % from optimal hiring condition
k_L=(fv_L_ss*(q_L_ss-w_L_ss))/(1-betaL*rho_L_ss); % from optimal hiring condition

% dividends
dw_H_ss=(q_H_ss-w_H_ss)*e_H_ss-k_H*v_H_ss;
dw_L_ss=(q_L_ss-w_L_ss)*e_L_ss-k_L*v_L_ss;
D_H_ss=Y_ss-q_H_ss*Y_ss+dw_H_ss;
D_L_ss=Y_ss-q_L_ss*Y_ss+dw_L_ss;


% consumption target: c_u/c_e=0.72 
c_eH_ss=((1-tau_ss)*(e_H_ss*w_H_ss+D_H_ss)+tau_uH_ss*u_H_ss)/(e_H_ss+0.72*u_H_ss);
c_eL_ss=((1-tau_ss)*(e_L_ss*w_L_ss+D_L_ss)+tau_uL_ss*u_L_ss)/(e_L_ss+0.72*u_L_ss);

b_barH_ss=0.72*c_eH_ss-tau_uH_ss;
b_barL_ss=0.72*c_eL_ss-tau_uL_ss;
c_uH_ss=0.72*c_eH_ss;
c_uL_ss=0.72*c_eL_ss;

GammaH_ss=c_eH_ss-b_barH_ss+ c_eH_ss*(log(c_uH_ss/c_eH_ss)+chi_bar_H);
GammaL_ss=c_eL_ss-b_barL_ss+ c_eL_ss*(log(c_uL_ss/c_eL_ss)+chi_bar_L);
xH_ss=b_barH_ss;
bH_ss=b_barH_ss;
aH_ss=b_barH_ss;

xL_ss=b_barL_ss;
bL_ss=b_barL_ss;
aL_ss=b_barL_ss;

C=e_H_ss*c_eH_ss+e_L_ss*c_eL_ss+u_H_ss*c_uH_ss+u_L_ss*c_uL_ss;

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


save vecSS_2sectorsA betaH betaL alphamH alphamL alpha eps_p eps_s gamma vartheta chi_bar_H chi_bar_L ...
                  theta psi z_ss tau_uH_ss tau_uL_ss i_bar fs_H_ss fs_L_ss fv_H_ss fv_L_ss ...
                  rho_H_ss rho_rH sigma_rH rho_L_ss rho_rL sigma_rL rho_b sigma_b rho_z sigma_z rho_i sigma_i sigma_tauH sigma_tauL ...
                  w_H_ss w_L_ss b_barH_ss b_barL_ss e_H_ss e_L_ss u_H_ss u_L_ss s_H_ss s_L_ss k_H k_L v_H_ss v_L_ss D_H_ss D_L_ss dw_H_ss dw_L_ss...
                  xH_ss xL_ss aH_ss aL_ss bH_ss bL_ss c_eH_ss c_eL_ss c_uH_ss c_uL_ss p_aL_ss p_bL_ss p_aH_ss p_bH_ss varpsi_H_ss varpsi_L_ss...
                  GammaH_ss GammaL_ss q_H_ss q_L_ss tau_ss Y_ss Y_H_ss Y_L_ss 


function y = Findvarw(x, betaH, betaL, vartheta, rho_H_ss, rho_L_ss, u_H_ss, u_L_ss, e_H_ss, e_L_ss, q_H_ss, q_L_ss, v_H_ss, v_L_ss, fv_H_ss, fv_L_ss, fs_H_ss, fs_L_ss, constH, constL, const2H, const2L)
y =[x(1)-((u_H_ss*0.3*x(2)+u_L_ss*0.3*x(3))/(const2H+const2L-v_H_ss*((fv_H_ss*(q_H_ss-x(2)))/(1-betaH*rho_H_ss)) - v_L_ss*((fv_L_ss*(q_L_ss-x(3)))/(1-betaL*rho_L_ss))));
    x(2)*(1+((vartheta*betaH*rho_H_ss*fs_H_ss)/(1-betaH*rho_H_ss)))-(vartheta*q_H_ss)*(1+((betaH*rho_H_ss*fs_H_ss)/(1-betaH*rho_H_ss)))-((1-vartheta)/(1-x(1)))*((constH/(e_H_ss+0.72*u_H_ss))*((1-x(1))*(const2H-v_H_ss*((fv_H_ss*(q_H_ss-x(2)))/(1-betaH*rho_H_ss)))+0.3*x(2)*u_H_ss)+0.3*x(2));
    x(3)*(1+((vartheta*betaL*rho_L_ss*fs_L_ss)/(1-betaL*rho_L_ss)))-(vartheta*q_L_ss)*(1+((betaL*rho_L_ss*fs_L_ss)/(1-betaL*rho_L_ss)))-((1-vartheta)/(1-x(1)))*((constL/(e_L_ss+0.72*u_L_ss))*((1-x(1))*(const2L-v_L_ss*((fv_L_ss*(q_L_ss-x(3)))/(1-betaL*rho_L_ss)))+0.3*x(3)*u_L_ss)+0.3*x(3))];
end