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

function flag = Model
x0_P = 1;
finalP = fsolve(@RunModel,x0_P);
P_ss = finalP;
[test, ~,vecSS] = RunModel(finalP);
flag = 'done';
 
function [test, par, vecSS] = RunModel(GuessP)
%parameters 

alpham     = 0.6;
alpha      = 0.5;
vartheta   = 0.5;
varrho     = 1;   
chi_bar    = 0.42;
%xi_bar    =1;
iota       = 1;
phi        = 1;
eps_p      = 6;
eta        = 1;
theta      = 0.2;
psi        = 2;


rho_ss     = 0.9644; 
rho_r      = 0.8437;
sigma_r    = 0.0008;
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 P
P_ss = GuessP; 

% SOLVE FOR STEADY STATE 

z_ss=1;
pi_ss=1;
varsigma=1;
w=5.15;
tau_u=0.5*w;

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

% Labor market 

l_ss=0.645;

fs_ss=0.3;
u_ss=((1-fs_ss)*(1-rho_ss)/(1-rho_ss*(1-fs_ss)))*l_ss;

% or assume average unemployment rate of 6% and determine job finding probability fs
% ur=0.06;
% fs_ss=(1-rho_ss)*(ur-1)/(rho_ss-1-rho_ss*ur);
% u_ss=ur*l_ss;

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

fv_ss=alpham*(fs_ss/alpham)^(-alpha/(1-alpha));

% Consumption:

% start with a guess for b_bar and then compute it so that c_u/c_e=0.72
 
x0_bbar = 0.4; 
[vecbbar,fval, exitflag] = fsolve(@Findbbar,x0_bbar);
b_bar_ss=vecbbar;

x_ss=b_bar_ss;
c_n_ss=x_ss;
c_u_ss=x_ss+(tau_u/P_ss);
% c_e_ss is obtained from Euler equation:
c_e_ss= ((beta*(1+i_bar)/(1-beta*(1+i_bar)*e_ss))*(u_ss*(c_u_ss^(-varrho))+(1-l_ss)*(c_n_ss^(-varrho))))^(-1/varrho);
% tau_ss is obtained combining end-of-period asset constraint and govt budget
% constraint
tau_ss=(tau_u*u_ss)/P_ss*(e_ss*c_e_ss+u_ss*c_u_ss+(1-l_ss)*c_n_ss);

um_ce=1/c_e_ss; %marginal utility

% B is calibrated so that hiring costs are 1% of GDP
revper=0.01;
B=revper/(1-rho_ss)*(fs_ss)^eta;
k_ss=z_ss*B*(fs_ss)^eta;

u_cu=log(c_u_ss);
u_cn=log(c_n_ss);
u_ce=log(c_e_ss);

chi_ss=chi_bar*((e_ss)^(1+iota))/(1+iota);

% wage
%w=(vartheta/(1-vartheta))*(1-beta*rho_ss+beta*rho_ss*fs_ss)*(k_ss/fv_ss)+(1/(1-tau_ss))*(Gamma_ss+chi_bar*(e_ss^(1+iota))*(c_e_ss)^(-varrho));

Gamma_ss=c_e_ss-x_ss+(u_cu-u_ce+chi_ss)/um_ce;

%xi_bar is obtained from
xi_bar=((fs_ss/(1-beta*rho_ss*(1-fs_ss)*um_ce))*((1-tau_ss)*(w/P_ss)-chi_bar*(e_ss^(1+iota))*(c_e_ss)^(varrho)-Gamma_ss)+(log(c_u_ss)-log(c_n_ss))*(c_e_ss)^(varrho))*um_ce/(1+log(1-l_ss));

MRS_NC_ss=xi_bar/um_ce;
xi_ss=xi_bar*log(1-l_ss);
Omega_ss=(u_cu-(u_cn+xi_ss))/um_ce;

% prices  

q_ss=P_ss*(1-beta*rho_ss)*(k_ss/fv_ss)+w;

p_star_ss=q_ss*(eps_p-1)/eps_p;

Y_ss=z_ss*e_ss;

p_a_ss=((eps_p-1)/eps_p)*q_ss*Y_ss/(1-beta*(1-theta)) ;
p_b_ss= Y_ss/(1-beta*(1-theta));

%dividends
dw_ss=(q_ss/P_ss)*z_ss*e_ss-(w/P_ss)*e_ss-k_ss*v_ss;
D_ss=Y_ss-(q_ss/P_ss)*z_ss*e_ss+dw_ss;

b_ss=b_bar_ss*P_ss;
a_ss=b_ss;

test=p_star_ss-P_ss;

par = [beta alpham alpha vartheta varrho chi_bar xi_bar iota phi eps_p eta ...
                  theta psi z_ss pi_ss tau_u i_bar fs_ss fv_ss b_bar_ss B ...
                  rho_ss rho_r sigma_r rho_b sigma_b rho_z sigma_z rho_i sigma_i];

vecSS=[l_ss e_ss u_ss s_ss k_ss v_ss x_ss b_ss a_ss c_e_ss c_u_ss c_n_ss xi_ss chi_ss MRS_NC_ss Omega_ss Gamma_ss ... 
                  P_ss w p_a_ss p_b_ss q_ss Lambda tau_ss p_star_ss Y_ss D_ss dw_ss];

save param beta alpham alpha vartheta varrho chi_bar xi_bar iota phi eps_p eta ...
                  theta psi z_ss pi_ss tau_u i_bar fs_ss fv_ss b_bar_ss B ...
                  rho_ss rho_r sigma_r rho_b sigma_b rho_z sigma_z rho_i sigma_i 

save vecSS l_ss e_ss u_ss s_ss k_ss v_ss x_ss b_ss a_ss c_e_ss c_u_ss c_n_ss xi_ss chi_ss MRS_NC_ss Omega_ss Gamma_ss ... 
                  P_ss w p_a_ss p_b_ss q_ss Lambda tau_ss p_star_ss Y_ss D_ss dw_ss


function y = Findbbar(x)
global  beta w x_ss b_bar_ss tau_u_ss u_ss l_ss c_e_ss c_u_ss tau_u  i_bar P_ss e_ss varrho
y =0.72- c_u_ss/c_e_ss;
