function [ys,check] = m_w_cy_rp_steadystate(ys,exe)

% steady state file

% This version is based on deltau=tau_i-tau_c

global M_ beta gamma phi zeta nu varrhobar alpha etabar psi eps...
    delta tau_c rho_R r_infl nsig mc_f eps_w phi_w eRsig tau_i...
    sbar rhos ssig zbar rhoz zsig tmbar rhotm inflndx inflndx_w...
    r_ystar r_yprod r_dy r_y a deltau xibar rhoxi xsig tmsig inflbar...
    rhoarp rpsig rhomk sigpmk rhols gamma0 x m eD eE sD sE X dX RF...
    RFt R RB RBt RL l q MB MH rK dft theta lB dp dm d b k c i y t...
    h w v RS RW lH infl infl_w mc tm z s varrho mu eta xi yobs cobs...
    iobs qobs robs hobs lobs bobs rspr mobs inflobs arp mk ls rsprobs U

beta  = M_.params(1);
gamma  = M_.params(2);
phi  = M_.params(3);
zeta  = M_.params(4);
nu  = M_.params(5);
varrhobar  = M_.params(6);
alpha  = M_.params(7);
etabar  = M_.params(8);
psi  = M_.params(9);
eps = M_.params(10);
delta  = M_.params(11);
tau_c  = M_.params(12);
rho_R  = M_.params(13);
r_infl  = M_.params(14);
nsig  = M_.params(15);
mc_f  = M_.params(16);
eps_w  = M_.params(17);
phi_w  = M_.params(18);
eRsig  = M_.params(19);
tau_i = M_.params(20);
sbar  = M_.params(21);
rhos  = M_.params(22);
ssig  = M_.params(23);
zbar  = M_.params(24);
rhoz  = M_.params(25);
zsig  = M_.params(26);
tmbar  = M_.params(27);
rhotm  = M_.params(28);
inflndx  = M_.params(29);
inflndx_w = M_.params(30);
r_ystar  = M_.params(31);
r_yprod  = M_.params(32);
r_dy  = M_.params(33);
r_y  = M_.params(34);
a  = M_.params(35);
deltau  = M_.params(36);
xibar  = M_.params(37);
rhoxi  = M_.params(38);
xsig  = M_.params(39);
tmsig  = M_.params(40);
inflbar = M_.params(41);
rhoarp  = M_.params(42);
rpsig  = M_.params(43);
rhomk  = M_.params(44);
sigpmk  = M_.params(45);
rhols  = M_.params(46);
gamma0 = M_.params(47);

check = 0;

options = optimset('display','off',...
         'MaxIter',5000,'MaxFunEvals',5000,'TolX',1e-10,...
         'TolFun',1e-10,'algorithm','levenberg-marquardt');

m       = 0.12;
rB      = 0.015;
rK      = 0.07;
eD      = 0.5;	
tm	= tmbar;
varrho  = varrhobar;
eta     = etabar;
mu   = varrhobar/(1-varrhobar);
xi   = xibar;
x_ini   = [m;rB;rK;eD];

[x,fval,exitflag] = fsolve(@returns_deltau_ww,x_ini,options);

m       = x(1);
rB      = x(2);
rK      = x(3);
eD      = x(4);

tk = tm*(1-m)+tau_c*rK+tau_i*rB*(1-m)-tau_c*(delta+rB*(1-m));
RF = (1-tau_c)*rK+1-(1-tau_c)*delta;
RBt= 1+(1-tau_c)*rB;
eE = eD + (m+tm*(1-m))*(1+tk)/RF;
sD = 1/sbar*(log(eD)+0.5*sbar^2);
sE = 1/sbar*(log(eE)+0.5*sbar^2);
lB = 1+mu*normcdf(sE);

RFt= ((1+mu*normcdf(sE-sbar))/lB ...
   + (eD*normcdf(sD)-normcdf(sD-sbar))/((1-varrho)*lB))*RF;

X  = ((1-eta)*normcdf(sD-sbar)+deltau*eD*(1-normcdf(sD)))*RF;
dX = ((1-eta)*normpdf(sD-sbar)/(sbar*eD)...
   + deltau*(1-normcdf(sD)-normpdf(sD)/sbar))*RF;

rho_yk = rK/((1-alpha)*mc_f);
%rho_ck = rho_yk-delta;
dmk = -normcdf(sE)/(1-varrho)*(((normcdf(sE-sbar)-normcdf(sD-sbar))...
  / normcdf(sE)-eD)*RF+ tk -(m+tm*(1-m)));
rho_ck = rho_yk-delta-0*normcdf(sD-sbar)-0*varrho*dmk;
rho_hk = rho_yk.^(1/alpha);
c = ((1-a*beta)*(1-a)^-gamma*(alpha*mc_f/zeta)*rho_yk^((alpha-1-nu)/alpha)...
  *(rho_ck)^nu)^(1/(gamma+nu));
k = c/rho_ck;
y = rho_yk*k;
h = (rho_yk)^(1/alpha)*k;
q = 1;
l = q*k;
MB= beta;
MH= beta;
lB= 1+mu*normcdf(sE);
theta  = (1-tm)*lB./(1-(1-deltau)*beta*(1-normcdf(sD))); 
t = tk*k;
dp= (1-normcdf(sE))*(((1-normcdf(sE-sbar))/(1-normcdf(sE))-eD)*RF*k...
  + t - (m+tm*(1-m))*k);
dm= -normcdf(sE)/(1-varrho)*(((normcdf(sE-sbar)-normcdf(sD-sbar))...
  / normcdf(sE)-eD)*RF*k+ t -(m+tm*(1-m))*k);
d = dp-dm;
b = (1-m)*k;
i = delta*k;
w = alpha*mc_f*y/h;
v = d/(1-beta);
RS= (d+beta*v)/(beta*v)+varrho*dm/(beta*v);
RW= m*RS + (1+tm)*(1-m)*RBt;
RB= 1+rB;
R = 1/beta;
RL= R;
lH= (1-a*beta)*(1-a)^-gamma*c^-gamma;
infl=1;
infl_w=1;
z =1;
x = tm*(1-m);
s =sbar;
mc=mc_f;
W = (1/(1-gamma)*(c^(1-gamma)-1)-zeta/(1+psi)*h^(1+psi))/(1-beta);
dft   =(normcdf(sD));
yobs = 0;
cobs = 0;
iobs = 0;
qobs = 0;
robs  = 100*(1/beta*inflbar-1)*0;
hobs  = 0;
lobs  = 0;
mobs  = 0;
rspr  = 0;
bobs = 0;
inflobs = 100*(inflbar-1)*0;
arp = 1;
mk=0;
ls=0;
rsprobs=rspr;
U = (1/(1-gamma)*((c-a*c)^(1-gamma)-1)-zeta/(1+nu)*h^(1+nu));

ys = [x m eD eE sD sE X dX RF...
    RFt R RB RBt RL l q MB MH rK dft theta lB dp dm d b k c i y t ...
    h w v RS RW lH infl infl_w mc tm z s varrho mu eta xi yobs cobs ...
    iobs qobs robs hobs lobs bobs rspr mobs inflobs arp mk ls rsprobs U]';