%--------------------------------------------------------------------------
% Banks, Nonbanks, and Business Cycles
% Yvan Becard and David Gauthier, February 2023
% Steady State File
%--------------------------------------------------------------------------


function [ys,params,check] = BG_steadystate(~,~,M_,~) 
          
%--------------------------------------------------------------------------
% Initialize
%--------------------------------------------------------------------------

% Global variables
global Fe Fi Gammae1 Gammai Ge Gi Gipr omegae Rk sigmae sigmai options grid_len

% Load parameters
check = 0;
eval(['[' (strjoin(cellstr(M_.param_names),' ')) '] = deal(' strjoin(cellstr(num2str(M_.params)),' , ') ');']);

% Options for numerical solver
options  = optimset('display','none','Diagnostics','off','TolX',1e-7,'TolFun',1e-6);
grid_len = 200;
params = M_.params; 
if isnan(Fe_p); nendo = M_.orig_endo_nbr; ys = NaN(nendo,1); return; end

% To 
Np = length(M_.param_names);
for ii = 1 : Np                                         % for each parameter
    paramname = deblank(M_.param_names{ii});            % get parameter name
    eval([paramname,' = M_.params(' int2str(ii),');']); % get parameter value
end


%--------------------------------------------------------------------------
% Shocks
%--------------------------------------------------------------------------

epsil    = 1;
gammae   = 1;
lambdap  = lambdap_p;
lambdaw  = lambdaw_p;
muup     = 1;
muz      = muz_p;
nu       = 1;
zetac    = 1;
zetah    = 1;
zetai    = 1;
zetaih   = 1;
nu_up    = 0;
etae     = 1;
etai     = 1;
sigmab  = 1;


%--------------------------------------------------------------------------
% Auxiliary Expressions
%--------------------------------------------------------------------------

pitarget = pi_p;
pi       = pitarget;  
piw      = pitarget * muz;
pitilde  = pi^iotap_p * pitarget^(1-iotap_p);
pitildew = pi^iotaw_p * pitarget^(1-iotaw_p);


%--------------------------------------------------------------------------
% Straightforward Equations
%--------------------------------------------------------------------------

% Entrepreneur leverage
Le = Le_p;

% Normalize hours
l = 1;

% Utilization
u = 1;

% Eqn 1, 2, and 3
mc = 1 / lambdap_p;

% Eqn 16: Return on housing
Rh = pi - 1;

% Eqn 28: Patient household FOC deposits
R = pi_p * muz_p / betap_p - 1;

% Eqn 25: Optimal capital investment 
qk = 1;

% Normalize housing price
qh = 1;

% Eqn 29: Patient household FOC securities
Rs = (1 + R) * nu - 1;


%--------------------------------------------------------------------------
% Block: Solve for rk
%--------------------------------------------------------------------------
 
try
    
% Solve block rk    
[rk,~,exit] = fzero(@(rk) solve_rk(delta_p,Fe_p,Le,muk_p,pi_p,qk,R,rk,tauk_p,upsil_p),[0.005 0.2],options);
if exit ~= 1; error('rk not found'); end


%--------------------------------------------------------------------------
% Production
%--------------------------------------------------------------------------

% Eqn 11: Rental rate of capital
k = upsil_p * muz * l * (alpha_p * mc / rk)^(1/(1-alpha_p));

% Eqn 12: Demand for patient household labor
wplp = mc * kappa_p * (1 - alpha_p) * epsil * upsil_p^-alpha_p * (k / muz)^alpha_p * l^(1-alpha_p);

% Eqn 13: Demand for impatient household labor
wili = mc * (1 - kappa_p) * (1 - alpha_p) * epsil * upsil_p^-alpha_p * (u * k / muz)^alpha_p * l^(1-alpha_p);

% Fixed cost in production
theta = (k / (muz_p * upsil_p))^alpha_p * l^(1-alpha_p) * (1 - 1 / lambdap_p); 

% Eqn 17: Aggregate output
yz = (k / (muz_p * upsil_p))^alpha_p * l^(1-alpha_p) - theta;

% Eqn 14: Capital accumulation 
i = k * (1 - (1 - delta_p) / (upsil_p * muz_p));


%--------------------------------------------------------------------------
% Entrepreneurs
%--------------------------------------------------------------------------

% Eqn 42: Entrepreneur leverage
ne = k / Le;

% Eqn 43: Entrepreneur debt
be = k - ne;

% Eqn 40: Entrepreneur default cutoff
Re = omegae * (1 + Rk) * k / be - 1;

% Eqn 41: Entrepreneur net worth 
deltae_p  = (1 / (pi * muz) * (1 - Gammae1) * (1 + Rk) * k) / ne - 1;

% Eqn 44: Entrepreneur spread
Se = (1 + Re) / (1 + R);

% Resources for monitoring 
de = 0 * muk_p * Ge * (1 + Rk) * qk * k / (pi * muz);

% Derivative of entrepreneur probability of default
fe = normpdf((log(omegae) + (sigmae^2) / 2) / sigmae) / (sigmae * omegae);

% Eqn 47: Shadow bank FOC
qktilde = (1 - muk_p) * Ge / Fe * upsil_p / pi * qk * (1 + Rk) / (1 + R);


%--------------------------------------------------------------------------
% Impatient Households, Grid for omegai, First Try, Small Grid
%--------------------------------------------------------------------------

% Small grid
grid_ome = linspace(0.00001,0.99,grid_len);

% Solve Impatient FOCs
tbn = FOC_i(betai_p,Fi_p,muh_p,muz,grid_ome,pi,qh,R,Rh);

% Identify lowest omega that satisfies tbn < 0 (we want a change of sign to ensure a zero solution exists)
idx_ome = find(tbn<0, 1, 'last' );

% If no solution, display message, return NaN, and continue the search
if (isempty(idx_ome)) || (idx_ome==1)
nendo     = M_.orig_endo_nbr; 
ys        = NaN(nendo,1);
return
end


%--------------------------------------------------------------------------
% Impatient Households 
%--------------------------------------------------------------------------

% Solve for omegai
omegai = fzero(@(omegai) FOC_i(betai_p,Fi_p,muh_p,muz,omegai,pi,qh,R,Rh),[grid_ome(idx_ome-2) grid_ome(idx_ome+2)],options);

% Error message
if (isnan(omegai)); error('omegai is NaN'); end 

% Definition
Gipr = normpdf((log(omegai) + sigmai^2 / 2) / sigmai) / (sigmai);

% Eqn 47: Shadow bank FOC
qhtilde = (1 - muh_p) .* Gi ./ Fi .* qh / pi .* (1 + Rh) ./ (1 + R);

% Eqn 34: Bank participation constraint and Eqn 35: Impatient default cutoff
Ri = (1 + R) / (1 - Fi + Fi * pi * qhtilde / ((1 + Rh) * qh * omegai)) - 1;

% Eqn 30: Impatient household marginal utility of consumption 
cilambdazi = (muz - bci_p * betai_p) / ((muz - bci_p) * (1 + tauc_p));

% Eqn 31: Impatient household FOC physical housing 
rh = qh - omegai * (1 + Rh) * 1 / (1 + Ri) - betai_p / (pi * muz) * (1 - Gammai) * (1 + Rh) * qh;

% Impatient household FOC housing services
hilambdazi = zetah / rh;

% Eqn 35: Impatient household default cutoff
bilambdazi = omegai * (1 + Rh) * hilambdazi / (1 + Ri);

% Resources for monitoring impatient households
dilambdazi = 0 * muh_p * Gi * (1 + Rh) * hilambdazi / (pi * muz);

% Combine resource constraint, government budget constraint, and impatient household transfer from government to create an auxiliary variable
xi = (1 - kappa_p) * ((u * rk / upsil_p - 1 / pi * delta_p) * 1 / muz * k * tauk_p + (wili + wplp) * taul_p + (1 - etag_p) * (yz - 1 / (1 - etag_p) * i - de) * (tauc_p - etag_p / (1 - etag_p)) - i * etag_p / (1 - etag_p));

% Eqn 33: Impatient household budget constraint
lambdazi = ((1 + tauc_p) * cilambdazi +  hilambdazi * (qh - 1 / (pi * muz) * (1 - Gammai) * (1 + Rh) * qh) - bilambdazi + (1 - kappa_p) * (1 - etag_p) * (tauc_p - etag_p / (1 - etag_p)) * (dilambdazi)) / ((1 - taul_p) * wili + xi);

% Resource constraint and definitions of government spending and regulatory penalties on banks
clambdazi = (1 - etag_p) * (yz * lambdazi - de * lambdazi - dilambdazi) - i * lambdazi;

% Resource constraint, government budget constraint, and definitions of government spending and regulatory penalties on banks
tlambdazi = (u * rk / upsil_p - 1 / pi * delta_p) * 1 / muz * tauk_p * k * lambdazi + (wili * lambdazi + wplp * lambdazi) * taul_p + clambdazi * (tauc_p - etag_p / (1 - etag_p)) - i * lambdazi * etag_p / (1 - etag_p);

% Deduce
ci       = cilambdazi / lambdazi;
di       = dilambdazi / lambdazi;
hi       = hilambdazi / lambdazi;
bi       = bilambdazi / lambdazi;
t        = tlambdazi  / lambdazi;
c        = clambdazi  / lambdazi;

% Eqn 36: Impatient household net worth
ni = 1 / (pi * muz) * (1 - Gammai) * (1 + Rh) * qh * hi;

% Eqn 37: Impatient household spread
Si = (1 + Ri) / (1 + R);

% Impatient household leverage 
Li = qh * hi / ni;

% Derivative of impatient household probability of default
fi = normpdf((log(omegai) + (sigmai^2) / 2) / sigmai) / (sigmai * omegai);


%--------------------------------------------------------------------------
% Production
%--------------------------------------------------------------------------

% Government spending
g = etag_p / (1 - etag_p) * (c + i);

% Eqn 22: GDP
y = g + c + i;

% Eqn 19: Aggregate consumption
cp = c - ci;

% Eqn 26: Patient household marginal utility of consumption 
lambdazp = 1 / cp * (muz - bcp_p * betap_p) / (muz - bcp_p) / (1 + tauc_p);


%--------------------------------------------------------------------------
% Prices and Wages
%--------------------------------------------------------------------------

% Find psil_p by combining price equations 4, 5, 6, 7, 8, 9 and aggregate hours
psil_p = ((l / (lambdaw_p / (lambdazi * (1 - taul_p) * wili))) / ((lambdazp / (1 / (1 - betap_p * xiw_p) / wplp * lambdaw_p * (1 - betap_p * xiw_p) / ((1 - taul_p))))^(1/(1-1/kappa_p))))^(1-kappa_p);

% Eqn 7, 8 and 9
li = (lambdaw_p * psil_p / (lambdazi * (1 - taul_p) * wili))^(-1/(1+sigmal_p));

% Eqn 20: Aggregate hours
lp = (l / li^(1-kappa_p))^(1/kappa_p);

% Deduce
wp = wplp / lp;
wi = wili / li;

% Eqn 1: FOC price 1
Fp = zetac * lambdazp * yz / (1 - betap_p * xip_p * (pitilde / pi)^(1/(1-lambdap_p)));

% Eqn 2: FOC price 2
Kp = zetac * lambdazp * lambdap_p * yz * mc / (1 - betap_p * xip_p * (pitilde / pi)^(lambdap_p/(1-lambdap_p)));

% Eqn 4: Patient household FOC wage 1
Fwp = zetac * lambdazp * lp * (1 - taul_p) / (lambdaw_p * (1 - betap_p * xiw_p * muz^(lambdaw_p/(1-lambdaw_p)) * (1 / piw)^(lambdaw_p/(1-lambdaw_p)) * pitildew^(1/(1-lambdaw_p)) / pi));

% Eqn 5: Patient household FOC wage 2
Kwp = lp^(1+sigmal_p) / (1 - betap_p * xiw_p);

% Eqn 7: Impatient household FOC wage 1
Fwi = zetac * lambdazi * li * (1 - taul_p) / (lambdaw_p * (1 - betai_p * xiw_p));

% Eqn 8: Impatient household FOC wage 2
Kwi = li^(1+sigmal_p) / (1 - betai_p * xiw_p);

% Aggregate wage
w = kappa_p * wp + (1 - kappa_p) * wi;


%--------------------------------------------------------------------------
% Housing and Financial Sector
%--------------------------------------------------------------------------

% Eqn 27: Patient household FOC housing
hp = muz / (muz * (1 - betap_p / muz) * lambdazp); 

% Eqn 21: Aggregate housing
h = hp + hi;

% Eqn 45: Shadow bank budget constraint
s = Fi / muz * qhtilde * hi + Fe / (upsil_p * muz) * qktilde * k;

% Eqn 49: Bank loans
b = be + bi;

% Eqn 48: Bank balance sheet
d = b;

% Eqn 50: Bank spread
Ss = (1 + Rs) / (1 + R);

% Household leverage
Lh = (qh * h + b) / (qh * h + b - bi);


%--------------------------------------------------------------------------
% Observation Equations
%--------------------------------------------------------------------------

gdp_obs              = 1;
consumption_obs      = 1;
investment_obs       = 1;
hours_obs            = 1;
inflation_obs        = 1;
nominalrate_obs      = 1;
householdcredit_obs  = 1;
firmcredit_obs       = 1;
householdspread_obs  = 1;
firmspread_obs       = 1;
houseprice_obs       = 1;
shadowbank_obs       = 1;


%--------------------------------------------------------------------------
% Percent Deviation from Steady State
%--------------------------------------------------------------------------

b_       = log(b)                            * 100;
be_      = log(be)                           * 100;
bi_      = log(bi)                           * 100;
c_       = log(c)                            * 100;
ci_      = log(ci)                           * 100;
cp_      = log(cp)                           * 100;
Fi_      = log(Fi)                           * 100;
Fe_      = log(Fe)                           * 100;
h_       = log(h)                            * 100;
hi_      = log(hi)                           * 100;
hp_      = log(hp)                           * 100;
i_       = log(i)                            * 100;
qh_      = log(qh)                           * 100;
qhtilde_ = log(qhtilde)                      * 100;
qk_      = log(qk)                           * 100;
qktilde_ = log(qktilde)                      * 100;
k_       = log(k)                            * 100;
l_       = log(l)                            * 100;
li_      = log(li)                           * 100;
lp_      = log(lp)                           * 100;
Le_      = log(Le)                           * 100;
Lh_      = log(Lh)                           * 100;
Li_      = log(Li)                           * 100;
s_       = log(s)                            * 100;
ne_      = log(ne)                           * 100;
ni_      = log(ni)                           * 100;
w_       = log(w)                            * 100;
wi_      = log(wi)                           * 100;
wp_      = log(wp)                           * 100;
y_       = log(y)                            * 100;
pi_      = 1  * 100;
R_       = 1  * 100;
Re_      = 1  * 100;
Ri_      = 1  * 100;
Rs_      = 1  * 100;
Se_      = Se * 100;
Si_      = Si * 100;


%--------------------------------------------------------------------------
% Wrap Up
%--------------------------------------------------------------------------
params = (cellfun(@eval,cellstr(M_.param_names)));
ys     = (cellfun(@eval,cellstr(M_.endo_names)));

for jj = 1 : length(M_.endo_names)% for each variable
    varname = deblank(M_.endo_names{jj});       % get variable name
    eval(['ys(' int2str(jj) ') = ' varname ';']); 
    varnameval(jj) = eval(varname);               % get variable value
    assignin('base',varname,(varnameval(jj)));    % assign variable in workplace
%     if ~ismember(varname,{'Kp','Fp','Kw','Fw'}),disp([varname,' = ',num2str(varnameval((jj)),15),';']); end
end

%--------------------------------------------------------------------------
% If Error During "Try"
%--------------------------------------------------------------------------

% Catch error
catch lasst
display(lasst.message)

% Fill steady state vector with NaN
nendo     = M_.orig_endo_nbr;
ys        = NaN(nendo,1);
params    = NaN(M_.param_nbr,1); 
end


function [tbn] = FOC_i(betai_p,Fi_p,muh_p,muz,omegai,pi,qh,R,Rh)
global Fi fi Gammai Gammaipr Gi Gipr sigmai 

% Find sigmai for a given omegai
sigmai = findsigma(Fi_p,omegai);
 
% Definitions
Fi       = Fi_p;
fi       = normpdf((log(omegai) + (sigmai.^2) ./ 2) ./ sigmai) ./ (sigmai .* omegai);
Gi       = normcdf((log(omegai) + sigmai.^2 ./ 2) ./ sigmai - sigmai);
Gipr     = normpdf((log(omegai) + sigmai.^2 ./ 2) ./ sigmai) ./ (sigmai);
Gammai   = omegai .* (1 - Fi) + Gi;
Gammaipr = 1 - Fi;

% Eqn 47: Shadow bank FOC
qhtilde = (1 - muh_p) .* Gi ./ Fi .* qh ./ pi .* (1 + Rh) ./ (1 + R);

% Eqn 34: Bank participation constraint and Eqn 35: Impatient default cutoff
Ri = (1 + R) ./ (1 - Fi + Fi .* pi .* qhtilde ./ ((1 + Rh) .* qh .* omegai)) - 1;

% Eqn 32: Impatient household FOC leverage 
tbn = betai_p ./ (pi .* muz) .* (1 + R + fi .* (1 + Ri) .* (omegai - pi .* qhtilde ./ ((1 + Rh) ./ qh))) - 1;


function [tbn] = FOC_e(Fe_p,Le,muk_p,omegae,pi,qk,R,Rk,upsil_p)

% Find sigmae for a given omegae
sigmae = findsigma(Fe_p,omegae);
 
% Definitions
Fe        = Fe_p;
fe        = normpdf((log(omegae) + (sigmae.^2) / 2) ./ sigmae) / (sigmae .* omegae);
Ge        = normcdf((log(omegae) + sigmae.^2 / 2) ./ sigmae - sigmae);
Gammae1   = omegae .* (1 - Fe) + Ge;
Gammae1pr = 1 - Fe;
Gammae2   = omegae .* (1 - Fe) + Fe;
Gammae2pr = 1 - Fe + fe .* (1 - omegae);

% Eqn 47: Shadow bank FOC
qktilde = (1 - muk_p) .* Ge ./ Fe .* qk .* upsil_p ./ pi .* (1 + Rk) ./ (1 + R);

% Eqn 39: Bank participation constraint for entrepreneurs
Ra = (1 + Rk) .* (Gammae2 + (pi .* qktilde ./ (upsil_p .* (1 + Rk) .* qk) - 1) .* Fe) .* Le ./ (Le - 1) - 1;

% Eqn 38: Entrepreneur FOC
tbn = (1 - Gammae1) .* (1 + Rk) ./ (1 + Ra) - Gammae1pr ./ (Gammae2pr + (pi .* qktilde / (upsil_p .* (1 + Rk) * qk) - 1) .* fe) * 1 ./ Le;


function [sigma] = findsigma(F_p,omega)

% Express sigma as a function of F and omega
det = norminv(F_p).^2 - 2 .* log(omega); % Discriminant Delta = b^2 - 4ac

% Grid for sigma, two roots, choose highest
sigma = max(norminv(F_p) + sqrt(det),norminv(F_p) - sqrt(det));


function [tbn_rk] = solve_rk(delta_p,Fe_p,Le,muk_p,pi_p,qk,R,rk,tauk_p,upsil_p)
global Fe Gammae1 Gammae1pr Gammae2 Gammae2pr Ge omegae Ra Rk sigmae options grid_len

% Eqn 15: Return on capital 
Rk  = ((1 - tauk_p) * rk + 1 - delta_p) * pi_p / upsil_p + tauk_p * delta_p - 1;

% Grid for omegae
grid_ome =  linspace(0.001,0.9999,grid_len);

% FOC Entrepreneurs
tbn = FOC_e(Fe_p,Le,muk_p,grid_ome,pi,qk,R,Rk,upsil_p);

% Identify lowest omegae that satisfies tbn < 0
idx_ome = min(find(tbn<0)); if (isempty(idx_ome)) || (idx_ome==1); return;end

% Solve for omegae
[omegae,~,~] = fzero(@(omegae)FOC_e(Fe_p,Le,muk_p,omegae,pi,qk,R,Rk,upsil_p),[grid_ome(idx_ome-2) grid_ome(idx_ome+2)],options);

% Find sigmae given omegae
sigmae =  findsigma(Fe_p,omegae);

% Definitions
Fe        = Fe_p;
fe        = normpdf((log(omegae) + (sigmae^2) / 2) / sigmae) / (sigmae * omegae);
Ge        = normcdf((log(omegae) + sigmae.^2 / 2) ./ sigmae - sigmae);
Gammae1   = omegae * (1 - Fe) + Ge;
Gammae1pr = 1 - Fe;
Gammae2   = omegae * (1 - Fe) + Fe;
Gammae2pr = 1 - Fe + fe * (1 - omegae);

% Eqn 47: Shadow bank FOC
qktilde = (1 - muk_p) * Ge / Fe * upsil_p / pi * qk * (1 + Rk) / (1 + R);

% Eqn 39: Bank participation constraint for entrepreneurs
Ra = (1 + Rk) * (Gammae2 + (pi * qktilde / (upsil_p * (1 + Rk) * qk) - 1) * Fe) * Le / (Le - 1) - 1;

tbn_rk = Ra - R;

