% I. This code generates the steady state for the baseline model with no SI
% and no LFP


%parameters 

alpham     = 1;    % matching efficiency
alpha      = 0.5;  % matching elasticity
vartheta   = 0.5;  % bargaining power
varrho     = 1;    % risk aversion
eps_p      = 6;    % elasticity of substitution
theta      = 0.2;  % price rigidity
psi        = 2;    % Taylor rule, inflation coeff.

rho_ss     = 0.967; 
rho_r      = 0.3389;
sigma_r    = 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_tau  = 0.01;

% SOLVE FOR STEADY STATE 

z_ss=1;
pi_ss=1;
varsigma=1;
i_bar=0.003;

q_ss=(eps_p-1)/eps_p;
chi_bar    = (1/2)*q_ss;

% Labor market 

u_ss=0.062;
e_ss=1-u_ss;
s_ss=1-rho_ss*e_ss;

fs_ss=e_ss*(1-rho_ss)/(1-rho_ss*e_ss); 
tightness=(fs_ss/alpham)^(-alpha/(1-alpha));
v_ss=tightness*s_ss;
fv_ss=alpham*tightness^(-alpha);  

beta=0.72/(1+i_bar)*(0.72*e_ss+u_ss);
Lambda=beta;

% defining the constant terms to simplify wage equation
const=0.28+log(0.72)+chi_bar;
A1=0.3+(const*v_ss*fv_ss/(1-beta*rho_ss)*(e_ss+0.72*u_ss));
A2=(const*e_ss/(e_ss+0.72*u_ss))-(const*v_ss*fv_ss*q_ss/(e_ss+0.72*u_ss)*(1-beta*rho_ss));

x0_w = 0.8; 
[varw,fval,exitflag] = fsolve(@Findvarw,x0_w);
w_ss = varw;

k=fv_ss*(q_ss-w_ss)/(1-beta*rho_ss);
dw_ss=(q_ss-w_ss)*e_ss-k*v_ss;
D_ss=(1-w_ss)*e_ss-k*v_ss;
c_e_ss=(e_ss*w_ss+D_ss)/(e_ss+0.72*u_ss);
c_u_ss=0.72*c_e_ss;
tau_u_ss=0.3*w_ss;
b_bar_ss=0.72*c_e_ss-tau_u_ss;
Gamma_ss=w_ss*A1+A2;

tau_ss=u_ss*tau_u_ss/(e_ss*w_ss+D_ss);
x_ss=b_bar_ss;
b_ss=b_bar_ss;
a_ss=b_ss;
Y_ss=e_ss;
p_starP_ss=1;
p_a_ss=(eps_p/(eps_p-1))*Y_ss*q_ss/(1-beta*(1-theta)) ;
p_b_ss= Y_ss/(1-beta*(1-theta));
C=e_ss*c_e_ss+u_ss*c_u_ss;

par = [beta alpham alpha vartheta chi_bar eps_p ...
                  theta psi z_ss pi_ss tau_u_ss i_bar fs_ss fv_ss b_bar_ss ...
                  rho_ss rho_r sigma_r rho_b sigma_b rho_z sigma_z rho_i sigma_i sigma_tau];

vecSS=[e_ss u_ss s_ss k v_ss x_ss b_ss a_ss c_e_ss c_u_ss Gamma_ss ... 
                  p_starP_ss w_ss p_a_ss p_b_ss q_ss Lambda tau_ss Y_ss D_ss dw_ss];


save vecSS_test beta alpham alpha vartheta chi_bar eps_p ...
                  theta psi z_ss pi_ss tau_u_ss i_bar fs_ss fv_ss b_bar_ss ...
                  rho_ss rho_r sigma_r rho_b sigma_b rho_z sigma_z rho_i sigma_i sigma_tau ...
                  e_ss u_ss s_ss k v_ss x_ss b_ss a_ss c_e_ss c_u_ss Gamma_ss ... 
                  p_starP_ss w_ss p_a_ss p_b_ss q_ss Lambda tau_ss Y_ss D_ss dw_ss fv_ss


function y = Findvarw(x)
global  beta vartheta rho_ss fs_ss q_ss A1 A2 u_ss e_ss v_ss fv_ss
y = x*(1+(vartheta*beta*rho_ss*fs_ss/(1-beta*rho_ss)))-vartheta*q_ss*(1+(beta*rho_ss*fs_ss/(1-beta*rho_ss)))-(1-vartheta)*(A1*x+A2)/(1-(0.3*u_ss*x*(1-beta*rho_ss)/((1-beta*rho_ss)*e_ss-v_ss*fv_ss*q_ss+v_ss*fv_ss*x)));
end
