function [ys,params,check] = MCBB_BIRF_EM_PHI_steadystate(ys,exo,M_,options_)

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

% Model equations
acm_ss = 1;
at_ss = 1;
an_ss = 1;
pcm_ss = PCMBAR;
ca_ss = 0;
dstar_ss = DBAR;
rstar_ss = RSTAR;
tb_ss = rstar_ss*dstar_ss;
r_ss = (1 - THETAB*(1+RSTAR) - NUB*PHIB)/(PHIB*THETAB) + RSTAR;
qcm_ss = 1+ETA*(r_ss/(1+r_ss));
qt_ss = 1+ETA*(r_ss/(1+r_ss));
qn_ss = 1+ETA*(r_ss/(1+r_ss));
ucm_ss = qcm_ss*(DELTA + rstar_ss);
ut_ss = qt_ss*(DELTA + rstar_ss);
un_ss = qn_ss*(DELTA + rstar_ss);
kohcm = (ucm_ss/(pcm_ss*acm_ss*ALFACM))^(1/(ALFACM-1));
wcm_ss = ((1-ALFACM)*pcm_ss*acm_ss*kohcm^ALFACM)/(1+ETA*(r_ss/(1+r_ss)));
THETAT = ((GAMMAT*ut_ss)/(ALFAT*pcm_ss*(1+ETA*(r_ss/(1+r_ss)))))^GAMMAT;
koht = (ut_ss/(at_ss*THETAT*ALFAT))^(1/(ALFAT+GAMMAT-1));
wt_ss = ((1-ALFAT-GAMMAT)*THETAT*at_ss*(koht^(ALFAT+GAMMAT)))/(1+ETA*(r_ss/(1+r_ss)));

%Calculating pn, hn and kohn
options = optimset('Display', 'off');
x0 = [1;10;1;1];
x = fsolve(@(x)NTSS_v8(x,ETA,ALFAN,OMEGAN,CHI,DELTA,MU,OMEGABAR,koht,kohcm,at_ss,acm_ss,wt_ss,wcm_ss,pcm_ss,ALFAT,GAMMAT,THETAT,OMEGAT,ALFACM,OMEGACM,r_ss,un_ss,ut_ss,tb_ss,an_ss,BETTA,b),x0,options); % Call solver
pn_ss = x(1);
kohn = x(2);
hn_ss = x(3);
reer_ss = x(4);

hcm_ss = (wcm_ss/(reer_ss*(OMEGABAR/(1-BETTA*b))))^(1/(OMEGACM -1));
ht_ss = (wt_ss/(reer_ss*(OMEGABAR/(1-BETTA*b))))^(1/(OMEGAT -1));
kcm_ss = kohcm * hcm_ss;
icm_ss = DELTA*kcm_ss;
ycm_ss = acm_ss*kohcm^ALFACM*hcm_ss;
kt_ss = koht * ht_ss;
it_ss = DELTA*kt_ss;
yt_ss = at_ss*THETAT*(koht^(ALFAT+GAMMAT))*ht_ss;
cmt_ss = ((GAMMAT*ut_ss)/(ALFAT*pcm_ss*(1+ETA*(r_ss/(1+r_ss)))))*kt_ss;
tbcm_ss = pcm_ss*(ycm_ss-cmt_ss);
wn_ss = reer_ss*(OMEGABAR/(1-BETTA*b))*(hn_ss^(OMEGAN-1));
kn_ss = kohn * hn_ss;
in_ss = DELTA*kn_ss;
yn_ss = an_ss* (kohn^ALFAN) * hn_ss;
ivt_ss = icm_ss + it_ss + in_ss;
h_ss = hcm_ss + ht_ss + hn_ss;
ct_ss = ((CHI*pn_ss/(1-CHI))^MU)*yn_ss;
tbt_ss = yt_ss - ct_ss - ivt_ss;
cn_ss = yn_ss;
y_ss =  tbcm_ss + yt_ss + pn_ss*yn_ss;
c_ss = (CHI * ct_ss^(1-1/MU ) + (1-CHI) * cn_ss^(1-1/MU))^(1/(1-1/MU));
Uprime_ss = ((c_ss-b*c_ss)-OMEGABAR*(((hcm_ss^OMEGACM)/OMEGACM)+((ht_ss^OMEGAT)/OMEGAT)+((hn_ss^OMEGAN)/OMEGAN)))^(-SIGG);
A1_ss = 1/reer_ss;
A2_ss = (1-CHI)*((c_ss/yn_ss)^(1/MU)); 
la_ss = (1-b*BETTA)*Uprime_ss*A1_ss;
dcm_ss = ETA*(wcm_ss*hcm_ss + icm_ss)/pn_ss;
dt_ss = ETA*(pcm_ss*cmt_ss + wt_ss*ht_ss + it_ss)/pn_ss;
dn_ss = ETA*(wn_ss*hn_ss + in_ss)/pn_ss;
d_ss = dcm_ss + dt_ss + dn_ss;
dstarb_ss = pn_ss*d_ss/(PHIB/(PHIB-1));
dstarh_ss = dstar_ss - dstarb_ss;
nw_ss = pn_ss*d_ss - dstarb_ss;
tby_ss = tb_ss/y_ss;
cay_ss = ca_ss/y_ss;
dstary_ss = dstar_ss/y_ss;
dy_ss = d_ss/y_ss;
gdp_def_ss = 1;
crt_ss = pn_ss*d_ss/reer_ss;
yr_ss = y_ss;


c = log(c_ss);
ct = log(ct_ss);
cn = log(cn_ss);
h = log(h_ss);
hcm = log(hcm_ss);
ht = log(ht_ss);
hn = log(hn_ss);
kcm = log(kcm_ss);
kt = log(kt_ss);
kn = log(kn_ss);
cmt = log(cmt_ss);
ycm = log(ycm_ss);
yt = log(yt_ss);
yn = log(yn_ss);
y = log(y_ss);
acm = log(acm_ss);
at = log(at_ss);
an = log(an_ss);
wcm = log(wcm_ss);
wt = log(wt_ss);
wn = log(wn_ss);
ucm = log(ucm_ss);
ut = log(ut_ss);
un = log(un_ss);
icm = log(icm_ss);
it = log(it_ss);
in = log(in_ss);
ivt = log(ivt_ss);
qcm = log(qcm_ss);
qt = log(qt_ss);
qn = log(qn_ss);
tbcm = tbcm_ss;
tbt = tbt_ss;
tb = tb_ss;
ca = ca_ss;
crt = log(crt_ss);
d = log(d_ss);
dcm = log(dcm_ss);
dt = log(dt_ss);
dn = log(dn_ss);
dstarb = dstarb_ss;
dstarh = dstarh_ss;
nw = nw_ss;
dstar = dstar_ss;
tby = tby_ss;
cay = cay_ss;
dstary = dstary_ss;
dy = dy_ss;
rstar = rstar_ss;
r = r_ss;
spread = r_ss - rstar_ss;
pn = log(pn_ss);
la = log(la_ss);
A1 = log(A1_ss);
A2 = log(A2_ss);
Uprime = log(Uprime_ss);
pcm = log(pcm_ss);
reer = log(reer_ss);
yr = log(yr_ss);
gdp_def = log(gdp_def_ss);

y_obs = 0;
i_obs = 0;
tby_obs = 0;
crt_obs = 0;
rstar_obs = 0;
reer_obs = 0;
pcm_obs = 0;

AUX_ENDO_LEAD_226 = 0;
AUX_ENDO_LEAD_271 = 0;
AUX_ENDO_LEAD_313 = 0;
AUX_ENDO_LAG_49_1 = 0;

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

end

function F = NTSS_v8(x,ETA,ALFAN,OMEGAN,CHI,DELTA,MU,OMEGABAR,koht,kohcm,at_ss,acm_ss,wt_ss,wcm_ss,pcm_ss,ALFAT,GAMMAT,THETAT,OMEGAT,ALFACM,OMEGACM,r_ss,un_ss,ut_ss,tb_ss,an_ss,BETTA,b)
% pn = x(1); kohn = x(2); hn = x(3); reer = x(4)
F = [un_ss - x(1)*ALFAN*an_ss*(x(2)^(ALFAN-1)); 
     (1+ETA*(r_ss/(1+r_ss)))*(x(4)*(OMEGABAR/(1-BETTA*b))*x(3)^(OMEGAN-1)) - x(1)*(1-ALFAN)*an_ss*(x(2)^ALFAN);
     x(1) - ((1-CHI)/CHI)*(((at_ss*THETAT*(koht^(ALFAT+GAMMAT))*(wt_ss/(x(4)*(OMEGABAR/(1-BETTA*b))))^(1/(OMEGAT -1)) - DELTA*koht*(wt_ss/(x(4)*(OMEGABAR/(1-BETTA*b))))^(1/(OMEGAT -1)) + pcm_ss*(acm_ss*kohcm^ALFACM*(wcm_ss/(x(4)*(OMEGABAR/(1-BETTA*b))))^(1/(OMEGACM -1))-((GAMMAT*ut_ss)/(ALFAT*pcm_ss*(1+ETA*(r_ss/(1+r_ss)))))*koht*(wt_ss/(x(4)*(OMEGABAR/(1-BETTA*b))))^(1/(OMEGAT -1))) - DELTA*kohcm*(wcm_ss/(x(4)*(OMEGABAR/(1-BETTA*b))))^(1/(OMEGACM -1)) - tb_ss-DELTA*x(2)*x(3))/(an_ss*(x(2)^ALFAN)*x(3)))^(1/MU));
     x(4) - (CHI^MU + ((1-CHI)^MU)*(x(1)^(1-MU)))^(1/(1-MU))];
end