%% Steady State
clc;clear;
   
%----------------------------------------------------------------------------------------------------------
% Targets (2019 values)
%----------------------------------------------------------------------------------------------------------
Emp_US = 157536000;                                                     % U.S. Employment level from FRED  
Vac_US = 7161750;                                                       % U.S. Unfilled job vacancy from FRED
LF_US = 163537000;                                                      % U.S. Labor force from FRED

Nss = Emp_US/LF_US;                                                     % Implied employment rate
%Vss = Vac_US/LF_US;                                                     % Implied vacancy rate
THss = 0.5595;                                                          % Slack

Pi_ss = 1;                                                              % Inflation
Pi_star_ss = 1;                                                         % Optimal reset price
P_ss = 1;                                                               % Price level
Ass = 1;                                                                % Steady state TFP
v_pss = 1;                                                              % Price dispersion

Qss = 8.98;                                                             % Steady state LMFI


%----------------------------------------------------------------------------------------------------------
% Picked parameters
%----------------------------------------------------------------------------------------------------------
beta = (1.035)^(-1/4);                                                  % Discount rate
vth = 1.10^(1/4) - 1;                                                   % Depreciation rate
eta = 0.34;                                                             % Capital share
delta = 0.036;                                                          % Separation rate
omg = 0.5;                                                              % LMFI parameter
nu = 1;                                                                 % Adjustment cost paramter
chi = 3/4;                                                              % Price stickiness parameter
eps = 9;                                                                % Elasticity of substitution
psi = 0.541;                                                            % Disutility of labor
alpha_r = 0.5;                                                          % Index of real wage rigidity
phi_q = 0.25;                                                           % Match efficiency response to LMFI


%----------------------------------------------------------------------------------------------------------
% Estimated parameters
%----------------------------------------------------------------------------------------------------------
gam = 0.8833;                                                           % Elasticity of match
alpha = 0.4938;                                                         % Workers bargaining power


%----------------------------------------------------------------------------------------------------------
% Implied parameters: Rss r_kss vp_ss Kss Iss Yss Mss Uss THss Zss...
% Tauss Muss Kapss Css x_aux_1_ss x_aux_2_ss W_real_ss W_nash_ss

% Steady state implied parameter: kap_c
%----------------------------------------------------------------------------------------------------------
Rss = Pi_ss/beta;  
r_kss = 1/beta - 1 + vth;
vp_ss = (1-beta*chi*Pi_ss^(eps))*(eps-1)*Pi_star_ss /((1-beta*chi*Pi_ss^(eps-1))*eps*Pi_ss);
Kss = Nss*(r_kss/(eta*Ass*vp_ss))^(1/(eta-1));
Iss = vth*Kss;
Yss = Ass*Kss^(eta)*Nss^(1-eta)/v_pss;
Mss = Nss*delta;
Uss = 1 -(1-delta)*Nss;
Vss = THss*Uss;
Zss = Mss/(Uss^(gam)*Vss^(1-gam));
Tauss = Mss/Uss;
Muss = Mss/Vss;
kap_c = 2.844346501429049;
%kap_c = (1/Qss^(-omg))*(((1-alpha)*((1-eta)*vp_ss*Ass*Kss^eta*Nss^(-eta) - psi*(Yss-Iss)))/(((1-beta*(1-delta))/Muss) + alpha*beta*(1-delta)*THss -(1-alpha)*psi*Vss));
Kapss = kap_c*Qss^(-omg);
Css = Yss - Iss - Kapss*Vss;
W_nash_ss = alpha*((1-eta)*vp_ss*Ass*Kss^(eta)*Nss^(-eta) + beta*(1-delta)*Kapss*THss) + (1-alpha)*psi*Css;
W_real_ss = W_nash_ss;
x_aux_1_ss = vp_ss*Css^(-1)*Yss/(1-beta*chi*Pi_ss^(eps));
x_aux_2_ss= Css^(-1)*Yss/(1-beta*chi*Pi_ss^(eps-1));


baseline = [Yss Css Uss Nss Vss Kss Rss r_kss vp_ss Iss Mss Tauss Muss Kapss x_aux_1_ss x_aux_2_ss W_real_ss W_nash_ss Pi_ss Pi_star_ss v_pss];

func = @(x)[1-beta*x(7)/x(19); % (1) Consumption Euler equation
                1-beta*(x(8)+1-vth); % (2) Capital Euler equation
                (x(14)/x(13))-(x(14)/x(13))*beta*(1-delta)-(1-eta)*x(9)*Ass*(x(6)/x(4))^(eta) + x(17); % (3) Job creation
                x(8)-eta*x(9)*Ass*(x(6)^(eta-1))*x(4)^(1-eta); % (4) Rental rate
                x(18)-alpha*((1-eta)*x(9)*Ass*(x(6)/x(4))^(eta)+beta*(1-delta)*x(14)*THss)-(1-alpha)*psi*x(2); % (5) Nash wage
                x(1)-(Ass*x(6)^(eta)*x(4)^(1-eta))/x(21); % (6) Production
                x(10)-vth*x(6); % (7) Capital adjustment
                x(1)-x(2)-x(10)-x(14)*x(5); % (8) ARC
                x(11)-delta*x(4); % (9) Labor transition
                x(3)-1+(1-delta)*x(4); % (10) Unemployment
                x(11)-Zss*(x(3)^(gam))*x(5)^(1-gam); % (11) Matching
                x(12)-x(11)/x(3); % (12) Job Finding
                x(13)-x(11)/x(5); % (13) Job Filling
                THss-x(5)/x(3); % (14) Slack
                x(14)-kap_c*Qss^(-omg); % (15) Kappa
                x(19)^(1-eps)-chi-(1-chi)*x(20)^(1-eps); % (16) Agg Price Level
                x(21)*(1-chi*x(19)^(eps))-(1-chi)*(x(19)/x(20))^(eps); % (17) Price Dispersion
                x(15)-x(9)*x(1)/x(2)-beta*chi*x(19)^(eps)*x(15); % (18) Numerator
                x(16)-x(1)/x(2)-beta*chi*x(19)^(eps-1)*x(16); % (19) Denominator
                x(20)-eps/(eps-1)*x(19)*x(15)/x(16); % (20) Reset
                x(18)-x(17)]; % (21) Real to Nash
            
x = fsolve(func,baseline);

Y  = x(1);
C  = x(2);
U = x(3);
N = x(4);
V = x(5);
K = x(6);
R = x(7);
r_k = x(8);
vp = x(9);
I = x(10);
M = x(11);
Tau = x(12);
Mu = x(13);
Kap = x(14);
x_aux_1 = x(15);
x_aux_2 = x(16);
W_real = x(17);
W_nash = x(18);
Pi = x(19);
Pi_star = x(20);
v_p = x(21);


func = @(x)[1-beta*x(7)/x(19); % (1) Consumption Euler equation
                1-beta*(x(8)+1-vth); % (2) Capital Euler equation
                (x(14)/x(13))-(x(14)/x(13))*beta*(1-delta)-(1-eta)*x(9)*Ass*(x(6)/x(4))^(eta) + x(17); % (3) Job creation
                x(8)-eta*x(9)*Ass*(x(6)^(eta-1))*x(4)^(1-eta); % (4) Rental rate
                x(18)-alpha*((1-eta)*x(9)*Ass*(x(6)/x(4))^(eta)+beta*(1-delta)*x(14)*THss)-(1-alpha)*psi*x(2); % (5) Nash wage
                x(1)-(Ass*x(6)^(eta)*x(4)^(1-eta))/x(21); % (6) Production
                x(10)-vth*x(6); % (7) Capital adjustment
                x(1)-x(2)-x(10)-x(14)*x(5); % (8) ARC
                x(11)-delta*x(4); % (9) Labor transition
                x(3)-1+(1-delta)*x(4); % (10) Unemployment
                x(11)-Zss*(x(3)^(gam))*x(5)^(1-gam); % (11) Matching
                x(12)-x(11)/x(3); % (12) Job Finding
                x(13)-x(11)/x(5); % (13) Job Filling
                THss-x(5)/x(3); % (14) Slack
                x(14)-kap_c*Qss^(-omg); % (15) Kappa
                x(19)^(1-eps)-chi-(1-chi)*x(20)^(1-eps); % (16) Agg Price Level
                x(21)*(1-chi*x(19)^(eps))-(1-chi)*(x(19)/x(20))^(eps); % (17) Price Dispersion
                x(15)-x(9)*x(1)/x(2)-beta*chi*x(19)^(eps)*x(15); % (18) Numerator
                x(16)-x(1)/x(2)-beta*chi*x(19)^(eps-1)*x(16); % (19) Denominator
                x(20)-eps/(eps-1)*x(19)*x(15)/x(16); % (20) Reset
                x(18)-x(17)]; % (21) Real to Nash
            
x2 = fsolve(func,x);

[x' x2' baseline']






