function [ys,params,check] =rentalhousing_steadystate(ys,exe,M_,options_)
% function [ys,params,check] = BP2020_CES_steadystate(ys,exe,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)
 

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


rh = 1/BETAH; 
rm = rh - (1-DELTA);

lambdab = (1-BETAB*rh)/(1-BETAB*RHOD) ;

re = (1/BETAB)*(1-((1-BETAB)*RHOD+(1-RHOD)*GAMMAE)*lambdab) ; 

rs = (1/BETAB)*(1-((1-BETAB)*RHOD+(1-RHOD)*GAMMAS)*lambdab) ; 


lambdae = (1-BETAE*re)/(1-BETAE*RHOE);
lambdas = (1-BETAS*rs)/(1-BETAS*RHOS);

rk = 1/BETAE*(1-lambdae*(1-RHOE)*MK) - (1-DELTA) ;
rv = 1/BETAE*(1-lambdae*MH/re*(1-RHOE)) - 1   ;

oo1 = JEI/(1-BETAH) ;
oo2 = JEI/(1-BETAS-lambdas*(1-RHOS)*MS/rs) ;  
CAPITALGAMMASTAR = ( (1-THETAS)/(OMEGAR*THETAS) ) * ( OMEGAR/(1-BETAH)*(1-BETAS-lambdas*(1-RHOS)*MS/rs)  )^(KAPPAS)  ;
T1 = (THETAS)^(1/KAPPAS)/( (THETAS)^(1/KAPPAS) + (1-THETAS)^(1/KAPPAS)*(OMEGAR*CAPITALGAMMASTAR)^(1-1/KAPPAS) )  ;
T2 = oo2*T1 ;
T3 = 1/( 1 - T2*( (1-rs)*MS/rs - CAPITALGAMMASTAR*(1-BETAH) ) ) ; 
T4 = T2*T3*( GAMMAS*MS/rs*(rh-1) + CAPITALGAMMASTAR*(1-BETAH) ) ;
oo4 = GAMMAE*(rh-1)*( MH*NU/(re*X_SS*rv) + MK*ALPHA*MIU/(X_SS*rk) - MN*(1-ALPHA-NU)/( X_SS*(1+MN*(1-RHOE)*lambdae)))  + (X_SS-1)/X_SS ;
oo6 = rm - DELTA ;
oo7 = (1-ALPHA-NU)*(1-SIGMA)/(X_SS*(1+MN*(1-RHOE)*lambdae)) ;
oo8 = (1-MIU)*ALPHA/(X_SS*rm) ;

ns = 1/(1+T3*TAUS) ;

zeta1 = oo6*oo8/oo7 + 1 + oo4/oo7 + T4*SIGMA/(1-SIGMA)  ; %has changed

nh = 1/(1+TAUH*zeta1) ;

ctoys = T3*(1-ALPHA-NU)/(X_SS*(1+MN*(1-RHOE)*lambdae))*SIGMA ;
ctoyh = zeta1*oo7 ;

%iktoye = DELTA*ALPHA*MIU/(rk) ;
%iktoyh = DELTA*ALPHA*(1-MIU)/(rm) ;

denh = X_SS*oo1*ctoyh + X_SS*T2*ctoys*(1+CAPITALGAMMASTAR) + NU/(rv) ; %has changed
hh = X_SS*oo1*ctoyh/denh ;
hs = X_SS*T2*ctoys/denh ;
he = (NU/(rv))/denh ;
hr = hs*CAPITALGAMMASTAR ;
Z = OMEGAR*hr ;

hstilda = ( (THETAS)^(1/KAPPAS)*(hs)^(1-1/KAPPAS) + (1-THETAS)^(1/KAPPAS)*(Z)^(1-1/KAPPAS) )^(KAPPAS/(KAPPAS-1)) ;

y = (ALPHA*(1-MIU)/(rm*X_SS))^(ALPHA*(1-MIU)/(1-ALPHA))*(ALPHA*MIU/(rk*X_SS))^(ALPHA*MIU/(1-ALPHA))*he^(NU/(1-ALPHA))*(nh^(1-SIGMA)*ns^SIGMA)^((1-ALPHA-NU)/(1-ALPHA)) ;

kh = ALPHA*(1-MIU)*y/(X_SS*rm) ;
ke = ALPHA*MIU*y/(X_SS*rk) ;

q = NU*y/(X_SS*rv*he) ;

qr = q*(1-BETAH)/OMEGAR ; 

loe = MH*q/re*he + MK*ke - (MN*(1-ALPHA-NU)/(1+MN*(1-RHOE)*lambdae))*(y/X_SS) ;

wh = ((1-SIGMA)*(1-ALPHA-NU)/(1+MN*(1-RHOE)*lambdae))*((y/X_SS)/nh) ;
ws = (SIGMA*(1-ALPHA-NU)/(1+MN*(1-RHOE)*lambdae))*((y/X_SS)/ns);

los = MS/rs*T2*T3*ws*ns ;

d = GAMMAE*loe + GAMMAS*los ;

ch = (rm-DELTA)*kh + (rh-1)*d + wh*nh + qr*OMEGAR*hr + y*(X_SS-1)/X_SS ;     

cb = (re-1)*loe + (rs-1)*los - (rh-1)*d ;

ce = y/X_SS - DELTA*ke - (re-1)*loe - wh*nh - ws*ns - rm*kh ;

cs = T3*ws*ns ;

dp = 1;

ctot = ch + cs + ce + cb ;
itot = DELTA*(kh+ke) ;

kb = loe + los - d;

ktot = ke+kh;
ntot =  nh^(1-SIGMA)*ns^SIGMA  ;

GDP = itot+ctot;

inflation = 1;

inflation_SS = inflation; 


hr_ratio = (hr/(hs + hh)); 

price_rent_ratio = q/qr;


acdb=0;
acdh=0;
ache=0;
achh=0;
achs=0;
acke=0;
ackh=0;
acloeb=0;
aclosb=0;
acloee=0;
acloss=0;


Y_SS = (y);
ITOT_SS = (itot);
CTOT_SS = (ctot);
Q_SS = (q);
LTOT_SS = loe+los;
NTOT_SS = ntot;
KE_SS = ke;
KH_SS = kh;
RK_SS = rk;
D_SS = d;
LOE_SS = loe;
LOS_SS = los;
HH_SS = hh;
HS_SS = hs ;
HR_SS = hr ;
HE_SS = he ;
HMORT_SS = hs + he ;
HHO_SS = hh + hs ;
QR_SS = qr ;
HR_HTOT_SS = (hr_ratio) ;
QQ_QR_SS = price_rent_ratio;
HOWNER_SS = (hh + hs + he);


intshare = 1;

loetoy=loe/4/y;
lostoy=los/4/y;
lotoy=(loe+los)/4/y;

cbtoy = cb/y ;
cetoy = ce/y ;
chtoy = ch/y ;
cstoy = cs/y ;

inflation = 1;

inflation_SS = inflation; 

X = X_SS ;

div = y*(X_SS-1)/X_SS ;

kab=1-(d/(loe+los));

%data_LLE=0;
%data_LLH=0;
data_CC=0;
data_CH=0;
data_CS=0;
data_KAB=kab;
data_II=0;
data_LL=0;
data_LO=0;
data_LOE=0;
data_LOH=0;
data_QQ=0;
data_YY=0;
data_YL=0;
data_TFP=0;
data_NN=0;
data_HS=0;
data_HR=0;
data_HR_stock = 0;
data_HE=0;
data_QR=0;
data_HMORT=0;
data_HHO = 0; 
data_HH =0;
data_RS =0;
data_QQ_QR = 0;
data_HOWNER = 0;
data_HR_growth= 0;
data_QR_growth= 0;
data_QQ_growth= 0;
data_RH = 0 ;
data_infl = 0;

zke=1;
zkh=1;
dke=DELTA;
dkh=DELTA;

CE_SS = ce;
CH_SS = ch;
RE_SS = re;
LE_SS = lambdae;
LB_SS = lambdab;
CB_SS = cb;
CS_SS = cs;

ucb = cb^-SIGMAB;
uce = ce^-SIGMAE;
uch = ch^-SIGMAH;
ucs = cs^-SIGMAH;

%a_be=0;
%a_bh=0;
a_j=0;
a_k=0;
a_me=0;
a_mh=0;
a_p=0;
a_z=0;
a_s=0;

%e_a_be=0;
%e_a_bh=0;

%% end own model equations
for iter = 1:length(M_.params) %update parameters set in the file
    eval([ 'params(' num2str(iter) ',1) = ' M_.param_names{iter,:} ';' ])
end

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