function [ys,params,check] = test21_steadystate(ys,exo,M_,options)
% Inputs:
%   - ys        [vector] vector of initial values for the steady state of
%   the engogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Outputs:
%   - ys        [vector] vector of steady state values for the endogenous
%   variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] 0 if steady state computation worked and to
%                        1 if not (allows to impose restrictions on parameters)

%% Step 0: Initialize indicator and set options for numerical solver
check = 0;
options = optimset('Display', 'Final', 'TolX', 1e-10, 'TolFun', 1e-10);

%% Step 1: read out parameters to acess them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

%% Step 2: enter model equations here
% Calibration of some parameters:
phi = 1;
betta = 0.978;
betaE = 0.940;
alpha = 0.44;
deltaE = 0.025;
deltaB = 0.090;
hH = 0.74;
hE = 0.74;
mu = 0.20;
vB = 0.11;
kappaKB = 20;
chi = 3.42;
gamma = 0.33;
theta = 0.75;
epsilon = 6;
N = 5;
omegaB = 0.75;
phir = 0.79;
phipi = 2.43;
phiy = 0.16;
rhoGamma = 0.5;
rhoa = 0.9;
rhom = 0.75;
rhoKE = 0.5;
rhoKB = 0.9;
sigmaGamma = 0.0001;
sigmaa = 0.0025;
sigmam = 0.0001;
sigmaKE = 0.001;
sigmaKB = 0.01;
sigmaR = 0.0025;
sigma = 0.55;
phis = 0.10;

% Value of some variables in steady state:
Psi_t = 1;
Gamma_t = 1;
Pi_t = 1;
Pistar_t = 1;
M_t = 0.30;     % Internal calibration (and close to Ferreira(2013))
A_t = 1;
ZK_t = 1;
SKB_t = 1;
Q_t = 1;
TCP_t = 0;
X_t = epsilon/(epsilon - 1);

% Variables that depends on other endogenous variables in steady-state:
R_t = Pi_t/betta;
F3_t = (1 - theta)/(1 - theta*(Pi_t)^(epsilon - gamma*epsilon));
W_t_times_L_t_over_Yw_t = (1 - alpha)/X_t;
Lambda_tt1 = betta;
Phi_tt1 = betaE;
I_t_over_KE_t = deltaE;
Yw_t_over_Y_t = F3_t;
TR_t_over_Y_t = ((1/F3_t) - (1/X_t))*Yw_t_over_Y_t;

% Here, we need to compute KBss/Bss, RBss and Omegabarss simultaneously,
% because there's no way to isolate one variable as a function of another.
% That's why we need a function to compute these three values numerically:

% Initial guess to KBss/Bss, RBss and Omegabar_t:
initial_guess = [0.11, 1.0225, 0.50];

% Use fsolve to compute KBss/Bss and RBss:
options = optimoptions('fsolve', 'Display', 'iter');
solutions = fsolve(@complex_equations, initial_guess, options);

% The solutions of this non-linear system are:
KB_t_over_B_t_solution = solutions(1);
RB_t_solution = solutions(2);
Omegabar_t_solution = solutions(3);

% Define a function that returns the solutions:
function equations = complex_equations(x)
    KB_t_over_B_t = x(1);
    RB_t = x(2);
    Omegabar_t = x(3);
    equations(1) = KB_t_over_B_t - (1 - deltaB)*KB_t_over_B_t*(1/Pi_t) -...
        (omegaB/Pi_t)*(1 - normcdf((log(Omegabar_t)+sigma^2/2)/sigma))*RB_t + (omegaB/Pi_t)*R_t -...
        (omegaB/Pi_t)*KB_t_over_B_t*R_t + omegaB*KB_t_over_B_t +...
        (omegaB/Pi_t)*(kappaKB/2)*((KB_t_over_B_t - vB)^(2))*KB_t_over_B_t -...
        omegaB*normcdf(((log(Omegabar_t)+sigma^2/2)/sigma) - sigma)*mu*Q_t*(1 - deltaE)*(RB_t/(Omegabar_t*Pi_t*Q_t*(1 - deltaE)));
    equations(2) = RB_t*(1 - normcdf((log(Omegabar_t)+sigma^2/2)/sigma))*(1 - (1/N)*((1 + (1/(1 - alpha))*(M_t*(1 - normcdf(((log(Omegabar_t)+sigma^2/2)/sigma) - sigma))*Q_t*(1 - deltaE)*Pi_t/RB_t)/(Q_t - betaE*(1 - normcdf(((log(Omegabar_t)+sigma^2/2)/sigma) - sigma))*Q_t*(1 - deltaE) -...
        (M_t*(1 - normcdf(((log(Omegabar_t)+sigma^2/2)/sigma) - sigma))*Q_t*Pi_t*(1 - deltaE)/RB_t) + betaE*(1 - normcdf(((log(Omegabar_t)+sigma^2/2)/sigma) - sigma))*M_t*(1 - normcdf((log(Omegabar_t)+sigma^2/2)/sigma))*Q_t*(1 - deltaE)))^(-1))) - R_t +...
        kappaKB*(KB_t_over_B_t - vB)*((KB_t_over_B_t)^(2));
    equations(3) = Omegabar_t - M_t*(1 - normcdf(((log(Omegabar_t)+sigma^2/2)/sigma) - sigma));
end

KB_t_over_B_t = KB_t_over_B_t_solution;
RB_t = RB_t_solution;
Omegabar_t = Omegabar_t_solution;

% Continue with variables in steady-state that depends on RBss:
S_t = RB_t - R_t;
F_t = normcdf((log(Omegabar_t)+sigma^2/2)/sigma);
G_t = normcdf(((log(Omegabar_t)+sigma^2/2)/sigma) - sigma);
lambda2E_t_times_CE_t_habit = 1 - betaE*(1 - F_t)*(RB_t/Pi_t);
KE_t_over_Yw_t = (betaE*alpha)/(X_t*(Q_t - betaE*(1 - G_t)*Q_t*(1 - deltaE) - lambda2E_t_times_CE_t_habit*M_t*(1 - G_t)*Q_t*Pi_t*(1 - deltaE)/RB_t));
MPKE_t = alpha*(1/X_t)*((KE_t_over_Yw_t)^(-1));
PED_t = 1 + (1/(1 - alpha))*(M_t*(1 - G_t)*Q_t*(1 - deltaE)*Pi_t/RB_t)/(betaE*MPKE_t);
B_t_over_Yw_t = (M_t*(1 - G_t)*Q_t*Pi_t*(1 - deltaE)/RB_t)*(KE_t_over_Yw_t);
KB_t_over_Yw_t = ((omegaB/Pi_t)*(B_t_over_Yw_t)*((1 - F_t)*RB_t - R_t) + omegaB*G_t*mu*Q_t*(1 - deltaE)*KE_t_over_Yw_t)/...
    (1 - (1 - deltaB)*(1/Pi_t) - (omegaB/Pi_t)*R_t + omegaB + (omegaB/Pi_t)*(kappaKB/2)*((KB_t_over_B_t - vB)^2));
TB_t_over_Yw_t = (1 - F_t)*(1/Pi_t)*RB_t*(B_t_over_Yw_t) - (1/Pi_t)*R_t*(B_t_over_Yw_t) + (R_t/Pi_t)*(KB_t_over_Yw_t) -...
    KB_t_over_Yw_t - (kappaKB/(2*Pi_t))*((KB_t_over_B_t - vB)^2)*(KB_t_over_Yw_t);

% Re-writing some steady-state ratios:
I_t_over_Y_t = (I_t_over_KE_t)*(KE_t_over_Yw_t)*(Yw_t_over_Y_t);
B_t_over_Y_t = (B_t_over_Yw_t)*(Yw_t_over_Y_t);
KE_t_over_Y_t = (KE_t_over_Yw_t)*(Yw_t_over_Y_t);
W_t_times_L_t_over_Y_t = (W_t_times_L_t_over_Yw_t)*(Yw_t_over_Y_t);
KB_t_over_Y_t = (KB_t_over_B_t)*(B_t_over_Y_t);
TB_t_over_Y_t = (TB_t_over_Yw_t)*(Yw_t_over_Y_t);

% Balance sheet identity:
D_t_over_Y_t = B_t_over_Y_t - KB_t_over_Y_t;
% Household budget constraint:
C_t_over_Y_t = ((R_t/Pi_t) - 1)*D_t_over_Y_t + W_t_times_L_t_over_Y_t + TR_t_over_Y_t + (1 - omegaB)*(TB_t_over_Y_t + mu*G_t*Q_t*(1 - deltaE)*KE_t_over_Y_t) +...
    (1 - mu)*G_t*Q_t*(1 - deltaE)*KE_t_over_Y_t; 
% Aggregate equilibrium:
CE_t_over_Y_t = 1 - C_t_over_Y_t - I_t_over_Y_t - (1/Pi_t)*(kappaKB/2)*((KB_t_over_B_t - vB)^(2))*KB_t_over_Y_t;
% Entrepreneurs' borrowing:
NE_t_over_Y_t = CE_t_over_Y_t + Q_t*KE_t_over_Y_t - B_t_over_Y_t;
% In the FOC of household, we can compute labor:
Y_t_over_C_t = (C_t_over_Y_t)^(-1);
L_t = ((W_t_times_L_t_over_Y_t)*(Y_t_over_C_t)*(1/(1 - hH)))^(1/(1 + phi));
% In production function of entrepreneur, we can compute Lss/Yss:
L_t_over_Yw_t = (A_t*((KE_t_over_Yw_t)^(alpha)))^(1/(alpha - 1));
L_t_over_Y_t = (L_t_over_Yw_t)*(Yw_t_over_Y_t);
% We can compute Yss:
Y_t = L_t/(L_t_over_Y_t);
% With Yss, we can compute all variables in steady state from the ratios:
F1_t = Y_t/(X_t*(1 - theta*betta));
F2_t = Y_t/(1 - theta*betta);
W_t = (W_t_times_L_t_over_Y_t)*(Y_t/L_t);
CE_t = (CE_t_over_Y_t)*Y_t;
NE_t = (NE_t_over_Y_t)*Y_t;
lambda2E_t = (lambda2E_t_times_CE_t_habit)/(CE_t - hE*CE_t);
KE_t = (KE_t_over_Y_t)*Y_t;
KB_t = (KB_t_over_Y_t)*Y_t;
Yw_t = (Yw_t_over_Y_t)*Y_t;
B_t = (B_t_over_Y_t)*Y_t;
I_t = (I_t_over_Y_t)*Y_t;
TR_t = (TR_t_over_Y_t)*Y_t;
C_t = (C_t_over_Y_t)*Y_t;
D_t = (D_t_over_Y_t)*Y_t;
TB_t = (TB_t_over_Y_t)*Y_t;
Lev_t = B_t/NE_t;
Rr_t = R_t/Pi_t;
log_Y_t = log(Y_t);
log_Pi_t = log(Pi_t);
% Welfare
WH_t = (1/(1 - betta))*(log(C_t - hH*C_t) - (L_t^(1 + phi))/(1 + phi));
WE_t = (1/(1 - betaE))*(log(CE_t - hE*CE_t));
WS_t = (1 - betta)*WH_t + (1 - betaE)*WE_t;

%% Step 3: Update parameters and variables
params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
end