% computes the steady state given parameters
% called by PH_steadystate, which uses this repeatedly to calibrate
% parameters

function [ys,params,check] = PH_steadystate_sublevel(ys,exe, M_)

% Here we load the values of the deep parameters in a loop.

for ii = 1:M_.param_nbr                                 
  eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);         %    Get the value of parameter i.
end                                                        % End of the loop.  
check = 0;
% END OF THE FIRST MODEL INDEPENDENT BLOCK.

%% THIS BLOCK IS MODEL SPECIFIC.

% temporary parameter
gv = ((1-gam)/vsig);

% Here we define the steady state.
A = Abar;
zeta = 1;
eta = etabar;
pie = pibar;
iSR = (1+pie) / bet - 1;
iT = iSR;
iLR = iSR;
iF = iSR;
kappa = (1+s)*(1-eta)/eta;
SDF = bet; 
ytil = 0;
Ktil = 0;

z = 0;
zI = 0;
zp = 0;
zw = 0;
alp = 0;
alpU = 0;

%% PH variables all zero
dPW = 0;
dPB = 0;
dPB_1 = 0;
PsiH = 0;

% depreciation vars
u = 1;
del = del0; 

% Taylor approximations
ups = -UT1;
delpie = DT1;

Ups = 0;

% solve for lamk given SDF
lamK = 1/SDF;

% hence we have mu
mu = lamK - (1-del)*SDF*lamK;

% given alp = 0 in steady state, we have mc
mc = (1+s)*(eta-1)/eta;

% given mc and mu, solve for w
%w = A*(1-the) * (((mu/the)^the)*mc)^(1/(1-the));

w = A * (1-the) * (mc * (the/mu)^the) ^ (1/(1-the));

% given w we have wH
wH = ((etaW-1)/etaW)*w;

% guess n to solve K
tsolve = @(guess) solve_steady_core(guess, mu, the, w, wH, del, A, h, psy);
options = optimoptions('fsolve', 'Display', 'off', 'StepTolerance', 1e-12, 'OptimalityTolerance', 1e-8);
n = fsolve(tsolve, 0.5, options);

K = n * (the/(1-the)) * (w/mu);

% given K we have I
I = K * del;

% this gives us output and consumption
y = K^the * (A*n)^(1-the);
c = y - I;

% c gives us cbar
cbar = (1-h)*c;

% now we solve profits

vpiM = y*(1+s-mc);
vpiI = mu*K-I;
vpiU = n*(w-wH);

VM = vpiM * bet / (1-bet);
VI = vpiI * bet / (1-bet);

WW = VM + VI + bbar*y; % stocks + government debt

% bond prices

QSR = 1/(1+iSR);
QLR = 1/(iLR+1-k);

% government debt
BLRG = (vthe*WW)/QLR;
BSRG = (bbar * y - vthe*WW) / QSR;

% bond clearing
BLR = 0;
BSR = 0;

BLRH = BLRG-BLR;
BSRH = BSRG-BSR;
BLR0 = 0;
BLR0_e = 0;

% government
CBS = 0;
t = s*y - CBS - QSR*BSRG - QLR*BLRG + (BSRG+(1+k*QLR)*BLRG) /(1+pie);

% utility normalised
U = 1;

UU = (cbar ^ psy) * (1-n) ^ (1-psy);

% derive norc implied by U
norc = U / ((1-bet)*UU^gv+bet*U^gv)^(1/gv); 

% given U, we have EtV and lam
EtV= U^(1-gam);
%lam = (1-bet)*zeta*(cbar/U)^(vsig-1);
lam = bet;

wal = y - c - I - zp; % walras law, must be zero

% observation variables
I_obs = 0;
y_obs = 0;
n_obs = 0;
QE_obs = QLR*BLR/y;
i_obs = iSR;
pie_obs = pie;





% END OF THE MODEL SPECIFIC BLOCK.

%% DO NOT CHANGE THIS PART.

% Here we define the steady state values of the endogenous variables of
% the model.
%
params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

for ii = 1:M_.orig_endo_nbr %auxiliary variables are set automatically
  eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
end
check = 0;
%%
%% END OF THE SECOND MODEL INDEPENDENT BLOCK.