function [ys,params,check] = Iacoviello_2010_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs:
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output:
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% Copyright © 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

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

TRENDY = TREND_AC + MUC/(1-MUC)*TREND_AK;  
TRENDH = (MUH+MUBB)*TREND_AC + (1-MUH-KAPPA-MUBB)*TREND_AH + MUC*(MUH+MUBB)/(1-MUC)*TREND_AK ;
TRENDK = TREND_AC + 1/(1-MUC)*TREND_AK ;
TRENDQ = (1-MUH-MUBB)*TREND_AC + MUC*(1-MUH-MUBB)/(1-MUC)*TREND_AK - (1-MUH-KAPPA-MUBB)*TREND_AH ;


% equations DEFINING CAPITAL TO PRODUCTION RATIOS (1) to (6) in "Steady state"
check=0;

% Trends
EXPTRENDY = exp ( TRENDY ) ;  
EXPTRENDK = exp ( TRENDK ) ;
EXPTRENDQ = exp ( TRENDQ ) ;
EXPTRENDH = exp ( TRENDH ) ;
gamma_k = exp ( TREND_AK );


% Interest rates
r = 1 / BETA ;
rkc = r*gamma_k - (1-DKC) ;
rkh = r - (1-DKH) ;
r1 =  r / EXPTRENDY - 1 ;

ZETA0 = BETA*EXPTRENDK*MUC/(gamma_k-BETA*(1-DKC))/X_SS ;
ZETA1 = BETA*EXPTRENDY*MUH/(1-BETA*(1-DKH));
ZETA2 = JEI/(1-BETA*EXPTRENDQ*(1-DH)) ;
ZETA3 = JEI/(1-BETA1*EXPTRENDQ*(1-DH)-EXPTRENDQ*(BETA-BETA1)*M) ;
ZETA4 = (r/EXPTRENDY-1)*M*EXPTRENDQ/r ;


% Usual...
R_SS = 1 / BETA ;
R1 = R_SS-1;

DH1 = 1 - (1-DH)/EXPTRENDH ;
DKC1 = 1 - (1-DKC)/EXPTRENDK ; 
DKH1 = 1 - (1-DKH)/EXPTRENDY ;

CHI1 = 1+DH1*ZETA2*(1-r1*ZETA1-KAPPA-ALPHA*(1-MUH-KAPPA-MUBB)) ;
CHI2 = (r1*ZETA1+KAPPA+ALPHA*(1-MUH-KAPPA-MUBB))*DH1*ZETA3+ZETA4*ZETA3 ;
CHI3 = (X_SS-1+r1*ZETA0*X_SS+ALPHA*(1-MUC))/X_SS ;
CHI4 = 1+DH1*ZETA3*(1-(1-ALPHA)*(1-MUH-KAPPA-MUBB))+ZETA4*ZETA3 ;
CHI5 = (1-ALPHA)*(1-MUH-KAPPA-MUBB)*DH1*ZETA2 ;
CHI6 = (1-ALPHA)*(1-MUC)/X_SS ;


CY = (CHI3*CHI4+CHI2*CHI6)/(CHI1*CHI4-CHI2*CHI5) ;
CYPRIME = (CHI1*CHI6+CHI3*CHI5)/(CHI1*CHI4-CHI2*CHI5) ;
QIY = DH1*ZETA2*CY + DH1*ZETA3*CYPRIME ;

RATION = (1-MUH-KAPPA-MUBB)/(1-MUC)*X_SS*QIY  ;
NHNC = RATION^(1/(1-NU)) ;
NHNC1 = RATION^(1/(1-NU1)) ;

nc = ( ((1-MUC)*ALPHA/CY/X_SS/XW_SS)/(1+RATION)^((ETA+NU)/(1-NU)) )^(1/(1+ETA)) ;
nh = NHNC*nc ;

nc1 = ( ((1-MUC)*(1-ALPHA)/CYPRIME/X_SS/XW_SS)/(1+RATION)^((ETA1+NU1)/(1-NU1)) )^(1/(1+ETA1)) ;
nh1 = NHNC1*nc1 ;

Y = (nc^ALPHA)*(nc1^(1-ALPHA)) *  ZETA0^(MUC/(1-MUC)) / EXPTRENDK^(MUC/(1-MUC)) ;


I = (nh^(ALPHA*(1-MUH-KAPPA-MUBB))) *  ...
    (nh1^((1-ALPHA)*(1-MUH-KAPPA-MUBB))) * ZETA1^MUH * (Y*QIY)^MUH / EXPTRENDY^(MUH) * (MUBB*Y*QIY)^MUBB ;

q = QIY*Y / I ;

QI = QIY*Y ;

kc = ZETA0*Y ;

kh = ZETA1*QI ;

BRI = MUBB*QI ;

c = CY*Y ;
c1 = CYPRIME*Y ;

h = ZETA2*c/q ;
h1 = ZETA3*c1/q ;


b = M*q*EXPTRENDQ*h1/r ;

lm = (1-BETA1/BETA)/c1 ;

uc = 1/c;
uc1 = 1/c1;
uh = JEI/h;
uh1 = JEI/h1;

CC = c + c1 ;
IH = I  ;
IK = DKC1 * kc + DKH1* kh  ;

YY_SS = CC + IK + BRI ;
QI_SS = q*I ;

WEALTH = (kc+kh) + q*(h+h1) ;

BB_SS = log(b) ;
CC_SS = log(CC) ;
IH_SS = log(IH) ;
IK_SS = log(IK) ;
QQ_SS = log(q) ;
RR_SS = log(R_SS) ;

NC_SS = ALPHA*log(nc) + (1-ALPHA)*log(nc1) ;
NH_SS = ALPHA*log(nh) + (1-ALPHA)*log(nh1) ;


dp = 1 ;

wc = ALPHA*(1-MUC)*Y/X_SS/nc ;
wh = ALPHA*(1-MUH-KAPPA-MUBB)*QI/nh ;
wc1 = (1-ALPHA)*(1-MUC)*Y/X_SS/nc1 ;
wh1 = (1-ALPHA)*(1-MUH-KAPPA-MUBB)*QI/nh1 ;

WC_SS = log(wc) ;
WC1_SS = log(wc1) ;
WH_SS = log(wh) ;
WH1_SS = log(wh1) ;


rr = 1 ;
c_ss = c ;

H_SS = log(h) ;
H1_SS = log(h1) ;

v = uh/uc*h/(h+h1) + uh1/uc1*h1/(h+h1) ;
V_SS = log(v) ;



xwc = XW_SS;
xwc1 = XW_SS;
xwh = XW_SS;
xwh1 = XW_SS;

WP_SS = log(wc + wh + wc1 + wh1) ;



% X = X_SS;
% data_CC = EXPTRENDY;
% data_DP = 1        ;
% data_GDP = 1       ;
% data_IH = EXPTRENDH;
% data_IK = EXPTRENDK;
% data_NC = 1        ;
% data_NH = 1        ;
% data_QQ = EXPTRENDQ;
% data_RR = 1        ;
% data_WC = 1        ;
% data_WH = 1        ;
X=X_SS;
a_c=1;
a_h=1;
a_j=1;
a_k=1;
a_s=1;
a_t=1;
a_z=1;

data_CC=EXPTRENDY^1; % C
data_DP=1;
data_IH=EXPTRENDH^1;
data_IK=EXPTRENDK^1;
data_NC=1;
data_NH=1;
data_QQ=EXPTRENDQ^1;
data_RR=1;
data_WC=1;
data_WH=1;
zata_GDP=1;
zkc=1;
zkh=1;

xxx = [ a_c     
a_h
a_j
a_k
a_s
a_t
a_z
b  
c  
c1 
data_CC
data_DP
data_IH
data_IK
data_NC
data_NH
data_QQ
data_RR
data_WC
data_WH
dp
h  
h1 
I  
kc 
kh 
lm 
nc 
nc1
nh 
nh1
q  
r  
rkc
rkh
uc 
uc1
wc 
wc1
wh
wh1
X_SS
xwc
xwc1
xwh
xwh1
Y   
zata_GDP
zkc
zkh];


if (any(xxx<=0))
  check = 1;
end


%% end own model equations

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) ') = log(' varname ');']);
end
