function [ys,params,check] = test2_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.66;
betta = 0.995;
betaE = 0.970;
alpha = 0.33;
eta = 0.2;
deltaA = 0.0295;
deltaB = 0.0100;
chi = 10;
gamma = 0.18;
theta = 0.75;
epsilon = 6;
N = 3;
tau = 0.10;

% Value of some variables in steady state:
Psi_t = 1;
Pi_t = 1;
Pistar_t = 1;
Gamma_t = 1;
MA_t = 0.8;
MB_t = 0.5;
A_t = 1;
ZKA_t = 1;
ZKB_t = 1;
QA_t = 1;
QB_t = 1;
TCPA_t = 0;
TCPB_t = 0;
RD_t = 1.010;
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_Y_t = (1 - alpha)/(X_t*(1 + tau));
lambda3E_t_times_CE_t = 1 - betaE*RD_t/Pi_t;
Lambda_tt1 = betta;
Phi_tt1 = betaE;
IA_t_over_KA_t = deltaA;
IB_t_over_KB_t = deltaB;
KB_t_over_Yw_t = (betaE*alpha*(1 - eta))/(X_t*(QB_t - betaE*QB_t*(1 - deltaB) - lambda3E_t_times_CE_t*MB_t*QB_t*Pi_t*(1 - deltaB)/RD_t));
Yw_t_over_Y_t = F3_t;
BD_t_over_Y_t = (MB_t*QB_t*Pi_t*(1 - deltaB)/RD_t)*KB_t_over_Yw_t*Yw_t_over_Y_t;
TR_t_over_Y_t = ((1/F3_t) - (1/X_t))*Yw_t_over_Y_t;

% Using the expressions for Ywss/KAss and MPKAss in PEDss, we obtain PEDss as
% function of RLss; substituting muss by an expression of FEDss in RLss =
% muss*Rss, we can write a (very large) expression of RLss using only
% parameters.
% We use fzero function to compute the value for RLss:
[RL_t,fval,exitflag]=fsolve(@(RL_t) RL_t*N + RL_t*N*((1 + alpha*eta - alpha)/(1 - alpha))*((MA_t*QA_t*Pi_t*(1 - deltaA))/...
    (RL_t*QA_t - betaE*RL_t*QA_t*(1 - deltaA) - MA_t*QA_t*Pi_t*(1 - deltaA) + betaE*RL_t*MA_t*QA_t*(1 - deltaA)))...
    - RL_t - N*R_t - N*R_t*((1 + alpha*eta - alpha)/(1 - alpha))*(MA_t*QA_t*Pi_t*(1 - deltaA)/(RL_t*QA_t - RL_t*betaE*QA_t*(1 - deltaA) ...
    - MA_t*QA_t*Pi_t*(1 - deltaA) + RL_t*betaE*MA_t*QA_t*(1 - deltaA))), 1.010, options);
if exitflag <1
    %indicate the SS computation was not sucessful; this would also be detected by Dynare
    %setting the indicator here shows how to use this functionality to
    %filter out parameter draws
    check=1; %set failure indicator
    return; %return without updating steady states
end


% Continue with variables in steady-state that depends on RLss:
S_t = RL_t - R_t;
PED_t = 1 + ((1 + alpha*eta - alpha)/(1 - alpha))*(MA_t*QA_t*Pi_t*(1 - deltaA)/(RL_t*QA_t - RL_t*betaE*QA_t*(1 - deltaA) ...
    - MA_t*QA_t*Pi_t*(1 - deltaA) + RL_t*betaE*MA_t*QA_t*(1 - deltaA)));
mu_t = 1/(1 - 1/(N*PED_t));
lambda2E_t_times_CE_t = 1 - betaE*(RL_t/Pi_t);
KA_t_over_Yw_t = (betaE*alpha*eta)/(X_t*(QA_t - betaE*QA_t*(1 - deltaA) - lambda2E_t_times_CE_t*MA_t*QA_t*Pi_t*(1 - deltaA)/RL_t));
MPKA_t = alpha*eta*(1/X_t)*((KA_t_over_Yw_t)^(-1));
BL_t_over_Yw_t = (MA_t*QA_t*Pi_t*(1 - deltaA)/RL_t)*KA_t_over_Yw_t;
TB_t_over_Y_t = (1/Pi_t)*(RL_t - R_t)*(BL_t_over_Yw_t)*(Yw_t_over_Y_t);
CE_t_over_Yw_t = ((1 - betaE)/betaE)*(QA_t*(KA_t_over_Yw_t) + QB_t*(KB_t_over_Yw_t) - BL_t_over_Yw_t - BD_t_over_Y_t);
NE_t_over_Yw_t = (1/(1 - betaE))*CE_t_over_Yw_t;

% Re-writing some steady-state ratios:
IA_t_over_Y_t = (IA_t_over_KA_t)*(KA_t_over_Yw_t)*(Yw_t_over_Y_t);
IB_t_over_Y_t = (IB_t_over_KB_t)*(KB_t_over_Yw_t)*(Yw_t_over_Y_t);
CE_t_over_Y_t = (CE_t_over_Yw_t)*(Yw_t_over_Y_t);
NE_t_over_Y_t = (NE_t_over_Yw_t)*(Yw_t_over_Y_t);
BL_t_over_Y_t = (BL_t_over_Yw_t)*(Yw_t_over_Y_t);
KA_t_over_Y_t = (KA_t_over_Yw_t)*(Yw_t_over_Y_t);
KB_t_over_Y_t = (KB_t_over_Yw_t)*(Yw_t_over_Y_t);

% Aggregate equilibrium
C_t_over_Y_t = 1 - CE_t_over_Y_t - IA_t_over_Y_t - IB_t_over_Y_t;
% Household's budget constraint
D_t_over_Y_t = ((W_t_times_L_t_over_Y_t) + (TR_t_over_Y_t) + (TB_t_over_Y_t) - (C_t_over_Y_t))/...
    (1 - (R_t/Pi_t));
% BNDES's budget constraint
DG_t_over_Y_t = (tau*(W_t_times_L_t_over_Y_t) + (RD_t/Pi_t)*(BD_t_over_Y_t)...
    - (BD_t_over_Y_t))/((R_t/Pi_t) - 1);
% 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 + phi));
% In production function of entrepreneur, we can compute Lss/Yss:
L_t_over_Yw_t = (A_t*((KA_t_over_Yw_t)^(alpha*eta))*((KB_t_over_Yw_t)^(alpha*(1 - eta))))^(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)/CE_t;
lambda3E_t = (lambda3E_t_times_CE_t)/CE_t;
KA_t = (KA_t_over_Y_t)*Y_t;
KB_t = (KB_t_over_Y_t)*Y_t;
Yw_t = (Yw_t_over_Y_t)*Y_t;
BL_t = (BL_t_over_Y_t)*Y_t;
BD_t = (BD_t_over_Y_t)*Y_t;
IA_t = (IA_t_over_Y_t)*Y_t;
IB_t = (IB_t_over_Y_t)*Y_t;
TB_t = (TB_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;
DG_t = (DG_t_over_Y_t)*Y_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

