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

global  betac sigman chik epsilon_n deltac deltac_g alphac_b alphac_g omegab omegag epsilon_k epsilon_b epsilon_g ... 
        phi_1 phi_2 phi_3 d0 d1 d2 delta_x em_row delta_nw iota_w taubar ab_ss ag_ss
        %%gamma_w gamma_bank iota_w rho_r phi_pi phi_y rho_ab rho_ag r_ss pi_ss sigmac gamma_I gamma_k

number_of_parameters = M_.param_nbr;

for ii=1:number_of_parameters
    paramname=M_.param_names{ii};
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
check = 0;

%%
%%epsa_ab=exo(strmatch('epsa_ab',M_.exo_names,'exact'));
options=optimset('MaxFunEvals',10000,'TolX',1e-10,'TolFun',1e-10);

%%epsa_ab_bar=epsa_ab;
M=betac;
r=1/betac;
pi=1;
qb=1;
xb=1;
xb1=1;
ab=ab_ss; %%*exp(epsa_ab_bar/(1-rho_ab));  %%ab=ab_ss
epsa_r=0;


kb_g = 10; %%10;
nb_g=0.3; %%0.3;
em_g=0.5;
nw_g=9.40; %%2.40;

guess = [kb_g, nb_g, em_g, nw_g];
ss_sol = fsolve(@findsstiming, guess, options);


kb = ss_sol(1);
nb=ss_sol(2);
em=ss_sol(3);
nw=ss_sol(4);

ib=deltac*kb;

x=(em+em_row)/(1-delta_x);
dmg=d0+d1*x+d2*(x^2);
yb=ab*(1-dmg)*(kb^alphac_b)*(nb^(1-alphac_b));

%%capb=kb;
dep=qb*kb(1)-nw;
%%rkb=(epsilon_b/(epsilon_b-1))*r;
rkb=r;
%%rkb=alphac_b*(yb/kb)+1-deltac;
prf=(rkb-r)*kb+r*nw;

wb=(1-alphac_b)*(yb/nb);
mc=0.5; %%(1/(ab*(1-dmg)))*(((rkb/alphac_b)^alphac_b)*((wb/(1-alphac_b))^(1-alphac_b))); %%0.5
c=yb-ib;

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

number_of_endogenous_variables = M_.orig_endo_nbr;
for ii = 1:number_of_endogenous_variables
    varname = M_.endo_names{ii};
    eval(['ys(' int2str(ii) ') = ' varname ';']);
end