function [ys,params,check] = test1_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;
beta = 0.99;
betaE = 0.95;
alpha = 0.33;
eta = 0.66;
deltaA = 0.0295;
deltaB = 0.0100;
chi = 1.7;
gamma = 0.18;
theta = 0.75;
epsilon = 6;
N = 4;
tau = 0.10;

% Variables of the model in steady-state:
c_t = 0;
w_t = 0;
l_t = 0;
r_t = 0;
pi_t = 0;
d_t = 0;
ThatCPA_t = 0;
ThatCPB_t = 0;
ThatR_t = 0;
ThatB_t = 0;
phi_t = 0;
lambda2hatE_t = 0;
cE_t = 0;
nE_t = 0;
rL_t = 0;
lambda3hatE_t = 0;
rD_t = 0;
qA_t = 0;
yw_t = 0;
kA_t = 0;
x_t = 0;
mA_t = 0;
qB_t = 0;
kB_t = 0;
mB_t = 0;
bL_t = 0;
bD_t = 0;
a_t = 0;
iA_t = 0;
iB_t = 0;
zKA_t = 0;
Lambdahat_tt1 = 0;
zKB_t = 0;
F3hat_t = 0;
y_t = 0;
dG_t = 0;
muhat_t = 0;
PEDhat_t = 0;
Phihat_tt1 = 0;
MPKAhat_t = 0;
g_t = 0;

% Value of some variables in steady state:
Psiss = 1;
Piss = 1;
MAss = 0.8;
MBss = 0.5;
Ass = 1;
ZKAss = 1;
ZKBss = 1;
QAss = 1;
QBss = 1;
TCPAss = 0;
TCPBss = 0;
RDss = 1.005;
Xss = epsilon/(epsilon - 1);
Gss_over_Yss = 0.05;

% Variables that depends on other endogenous variables in steady-state:
Rss = Piss/beta;
F3ss = (1 - theta)/(1 - theta*(Piss)^(epsilon - gamma*epsilon));
Wss_times_Lss_over_Yss = (1 - alpha)/(Xss*(1 + tau));
lambda3Ess_times_CEss = 1 - betaE*RDss/Piss;
IAss_over_KAss = deltaA;
IBss_over_KBss = deltaB;
KBss_over_Ywss = (betaE*alpha*(1 - eta))/(Xss*(QBss - betaE*QBss*(1 - deltaB) - lambda3Ess_times_CEss*MBss*QBss*Piss*(1 - deltaB)/RDss));
Ywss_over_Yss = F3ss;
BDss_over_Yss = (MBss*QBss*Piss*(1 - deltaB)/RDss)*KBss_over_Ywss*Ywss_over_Yss;
TRss_over_Yss = ((1/F3ss) - (1/Xss))*Ywss_over_Yss;

% 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:
[RLss,fval,exitflag]=fzero(@(RLss) RLss*N + (N/(1 - alpha*eta))*((MAss*QAss*Piss*(1 - deltaA))/...
    (QAss - betaE*QAss*(1 - deltaA) - MAss*QAss*Piss*(1 - deltaA)/RLss + betaE*MAss*QAss*(1 - deltaA)))*...
    (1 + (alpha*eta*(alpha - alpha*eta))/((1 - alpha)*(1 - alpha + alpha*eta))) - RLss - N*Rss ...
    - (N*Rss/(1 - alpha*eta))*(MAss*QAss*Piss*(1 - deltaA)/(RLss*QAss - RLss*betaE*QAss*(1 - deltaA) ...
    - MAss*QAss*Piss*(1 - deltaA) + RLss*betaE*MAss*QAss*(1 - deltaA)))*(1 + (alpha*eta*(alpha - alpha*eta)/...
    ((1 - alpha*eta)*(1 - alpha + alpha*eta)))), 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
RLss

% Continue with variables in steady-state that depends on RLss:
PEDss = 1 + (1/(1 - alpha*eta))*(MAss*QAss*Piss*(1 - deltaA)/(RLss*QAss - RLss*betaE*QAss*(1 - deltaA) ...
    - MAss*QAss*Piss*(1 - deltaA) + RLss*betaE*MAss*QAss*(1 - deltaA)))*(1 + (alpha*eta*(alpha - alpha*eta)/...
    ((1 - alpha*eta)*(1 - alpha + alpha*eta))));
muss = 1/(1 - 1/(N*PEDss));
lambda2Ess_times_CEss = 1 - betaE*(RLss/Piss);
KAss_over_Ywss = (betaE*alpha*eta)/(Xss*(QAss - betaE*QAss*(1 - deltaA) - lambda2Ess_times_CEss*MAss*QAss*Piss*(1 - deltaA)/RLss));
MPKAss = alpha*eta*(1/Xss)*((KAss_over_Ywss)^(-1));
BLss_over_Ywss = (MAss*QAss*Piss*(1 - deltaA)/RLss)*KAss_over_Ywss;
TBss_over_Yss = (1/Piss)*(RLss - Rss)*(BLss_over_Ywss)*(Ywss_over_Yss);
CEss_over_Ywss = ((1 - betaE)/betaE)*(QAss*(KAss_over_Ywss) + QBss*(KBss_over_Ywss) - BLss_over_Ywss - BDss_over_Yss);
NEss_over_Ywss = (1/(1 - betaE))*CEss_over_Ywss;

% Re-writing some steady-state ratios:
IAss_over_Yss = (IAss_over_KAss)*(KAss_over_Ywss)*(Ywss_over_Yss);
IBss_over_Yss = (IBss_over_KBss)*(KBss_over_Ywss)*(Ywss_over_Yss);
CEss_over_Yss = (CEss_over_Ywss)*(Ywss_over_Yss);
NEss_over_Yss = (NEss_over_Ywss)*(Ywss_over_Yss);
BLss_over_Yss = (BLss_over_Ywss)*(Ywss_over_Yss);
KAss_over_Yss = (KAss_over_Ywss)*(Ywss_over_Yss);
KBss_over_Yss = (KBss_over_Ywss)*(Ywss_over_Yss);

% Aggregate equilibrium
Css_over_Yss = 1 - CEss_over_Yss - IAss_over_Yss - IBss_over_Yss - Gss_over_Yss;
% Household's budget constraint
Dss_over_Yss = ((Wss_times_Lss_over_Yss) + (TRss_over_Yss) + (TBss_over_Yss) - (Css_over_Yss))/...
    (1 - (Rss/Piss));
% Government's budget constraint
DGss_over_Yss = (tau*(Wss_times_Lss_over_Yss) + (RDss/Piss)*(BDss_over_Yss) - (Gss_over_Yss) ...
    - (BDss_over_Yss))/((Rss/Piss) - 1);
% In the FOC of household, we can compute labor:
Yss_over_Css = (Css_over_Yss)^(-1);
Lss = ((Wss_times_Lss_over_Yss)*(Yss_over_Css))^(1 + phi);
% In production function of entrepreneur, we can compute Lss/Yss:
Lss_over_Ywss = (Ass*((KAss_over_Ywss)^(alpha*eta))*((KBss_over_Ywss)^(alpha*(1 - eta))))^(1/(alpha - 1));
Lss_over_Yss = (Lss_over_Ywss)*(Ywss_over_Yss);
% We can compute Yss:
Yss = Lss/(Lss_over_Yss);
% With Yss, we can compute all variables in steady state from the ratios:
Wss = (Wss_times_Lss_over_Yss)*(Yss/Lss);
CEss = (CEss_over_Yss)*Yss;
NEss = (NEss_over_Yss)*Yss;
lambda2Ess = (lambda2Ess_times_CEss)/CEss;
lambda3Ess = (lambda3Ess_times_CEss)/CEss;
KAss = (KAss_over_Yss)*Yss;
KBss = (KBss_over_Yss)*Yss;
Ywss = (Ywss_over_Yss)*Yss;
BLss = (BLss_over_Yss)*Yss;
BDss = (BDss_over_Yss)*Yss;
IAss = (IAss_over_Yss)*Yss;
IBss = (IBss_over_Yss)*Yss;
TBss = (TBss_over_Yss)*Yss;
TRss = (TRss_over_Yss)*Yss;
Css = (Css_over_Yss)*Yss;
Dss = (Dss_over_Yss)*Yss;
DGss = (DGss_over_Yss)*Yss;
Gss = (Gss_over_Yss)*Yss;

% Print the steady state variables
Rss
Xss
F3ss
Lss
Css
lambda2Ess
lambda3Ess
KAss
KBss
Wss
BLss
BDss
Ywss
IAss
IBss
Yss
TRss
Dss
TBss
RLss
muss
PEDss
MPKAss
DGss
Gss
CEss
NEss
            

% check if these equations were used:
%Lss = (Wss/Css)^(1/phi);
%Css = (Rss/Piss)*Dss + (Wss*Lss) + TCPAss + TCPBss + TRss + TBss - Dss;
%lambda2Ess = (1 - betaE*(RLss/Piss))/CEss;
%lambda3Ess = (1 - betaE*(RDss/Piss))/CEss;
%KAss = betaE*alpha*eta*Ywss/(Xss*(QAss - betaE*QAss*(1 - deltaA) - (CEss*lambda2Ess*MAss*QAss*Piss*(1 - deltaA)/RLss)));
%KBss = betaE*alpha*(1 - eta)*Ywss/(Xss*(QBss - betaE*QBss*(1 - deltaB) - (CEss*lambda3Ess*MBss*QBss*Piss*(1 - deltaB)/RDss)));
%Wss = ((1 + tau)*Xss*Lss)/((1 - alpha)*Ywss);
%NEss = alpha*(Ywss/Xss) + QAss*KAss*(1 - deltaA) + QBss*KBss*(1 - deltaB) - (RLss*BLss/Piss) - (RDss*BDss/Piss);
%CEss = (1 - betaE)*NEss;
%betaE*NEss = QAss*KAss + QBss*KBss - BLss - BDss;
%BLss = MAss*QAss*Piss*(1 - deltaA)*KAss/RLss;
%BDss = MBss*QBss*Piss*(1 - deltaB)*KBss/RDss;
%Ywss = Ass*((KAss)^(alpha*eta))*((KBss)^(alpha*(1 - eta)))*(Lss)^(1 - alpha);
%IAss = deltaA*KAss;
%IBss = deltaB*KBss;
%Yss = Ywss/F3ss;
%TRss = ((1/F3ss) - (1/Xss))*Ywss;
%Dss = BLss + DGss; % redundant equation; not used
%TBss = (1/Piss)*(RLss - Rss)*BLss;
%RLss = muss*Rss;
%muss = (1 - (1/N)*(PEDss^(-1)))^(-1);
%PEDss = 1 + (1/(1 - alpha*eta))*((MAss*QAss*Piss*(1 - deltaA))/(betaE*RLss*MPKAss))*(1 + (alpha*eta*(alpha - alpha*eta))/((1 - alpha*eta)*(1 - alpha + alpha*eta)));
%MPKAss = alpha*eta*(Ywss/KAss)*(1/Xss);
%DGss = (tau*(Wss*Lss) + (RDss/Piss)*BDss - Gss - BDss)/((Rss/Piss) - 1);
%Gss = Yss - Css - CEss - IAss - IBss;


%% 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

