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

global  betac sigman chik epsilon_n deltac deltac_g alphac_b alphac_g omegab omegag epsilon_k epsilon_y epsilon_b epsilon_g ... 
        phi_1 phi_2 phi_3 d0 d1 d2 delta_x em_row delta_nw iota_w taubar gamma_w taub taug gamma_wg iota_wg mu_dmg
        %%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;

%%

options=optimset('MaxFunEvals',10000,'TolX',1e-10,'TolFun',1e-10);

M=betac;
r=1/betac;
pi=1;
qb=1;
qg=1;
xb=1;
xb1=1;
xg=1;
xg1=1;
ab=1;
ag=1;
rkb=(1/betac)+deltac-1;
rkg=(1/betac)+deltac_g-1;
epsa_ab=0;
epsa_ag=0;
epsa_r=0;


kb_g = 8.10; 
kg_g = 5.5;  
nb_g= 0.8;   
ng_g= 0.6;    
em_g= 0.4;    

guess = [kb_g, kg_g, nb_g, ng_g, em_g];
ss_sol = fsolve(@findssmc_tax_ramsey, guess, options);


kb = ss_sol(1);
kg = ss_sol(2);
nb=ss_sol(3);
ng=ss_sol(4);
em=ss_sol(5);


n=(nb^(1+epsilon_n)+ng^(1+epsilon_n))^(1/(1+epsilon_n)); 
ib=deltac*kb;
ig=deltac_g*kg;
x=(em+em_row)/(1-delta_x);
dmg=1-exp(-mu_dmg*x);
yb=ab*(1-dmg)*(kb^alphac_b)*(nb^(1-alphac_b));
abt=1-(em/(phi_1*yb));
z=phi_2*(abt^phi_3)*yb;
tau=(phi_2*phi_3/phi_1)*(abt^(phi_3-1));
yg=ag*(1-dmg)*(kg^alphac_g)*(ng^(1-alphac_g));
y=((omegab^(1/epsilon_k))*yb^((epsilon_k-1)/epsilon_k)+(omegag^(1/epsilon_k))*yg^((epsilon_k-1)/epsilon_k))^(epsilon_k/(epsilon_k-1));
pb=(omegab*(y/yb))^(1/epsilon_k);
pg=(omegag*(y/yg))^(1/epsilon_k);


mc=(epsilon_y-1)/epsilon_y;
wb=(1-alphac_b)*(yb/nb)*(pb-phi_1*tau*(1-abt)-phi_2*(abt^phi_3));
wg=(1-alphac_g)*(yg/ng)*pg;
i=ib+ig;
c=y-i;


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