function [ys,check] = trigari_estimation_steadystate(ys,exe)
global M_ lgy_
global alpha beta delta sigma rho h epsilonp eta epsilong ry rpi etak b lambda lambdap gamma rhos gammaz kappa psinu rhob rhoi rhoz rhoeta bb rhop rhor rhog csi sigmab sigmai sigmaz sigmaeta sigmap sigmar sigmag y k n c i g x u s e Gamma pi rk qk pw a w chi mu H yc yi yg ynu yx htilda etanu betatilda kappaa kappaw kappalambda tau tau1 tau2 taup phia phis phix phib phichi gammab gammao gammaf fi smallc smallcp G iotab iotao iotaf fip psi phieta ;



if isfield(M_,'param_nbr') == 1
NumberOfParameters = M_.param_nbr;
for i = 1:NumberOfParameters
  paramname = deblank(M_.param_names(i,:));
  eval([ paramname ' = M_.params(' int2str(i) ');']);
end
check = 0;
end



rk = gammaz/beta-1+delta;

x = 1-rho;

n = (s/x)/(1+(s/x));

u = 1-n;

qk = 1;

k = ((rk/alpha)^(-1/(1-alpha)))*n;

a = (1-alpha)*(k/n)^alpha;

y = n*(k/n)^alpha;

i = y*(1-(1-delta)/gammaz)*gammaz*(k/n)^(1-alpha);

g = yg*y;

mu = 1/(1-lambda*beta);

e = 1/(1-rho*lambda*beta);

chi = eta/(eta+(1-eta)*mu/e);

pw = 1;    % to check this...

% fsolve
options = optimset('Display','off','TolFun',1e-50);
x0 = [5.4; y*(1-0.33); 0.87];   %initial values for fsolve
xd = fsolve(@kappa_w_trigari_test, x0, options, x, a, bb, beta, rho, chi, s, gammaz, pw);

kappa = xd(1);

w = xd(2);

b = xd(3);

c = (1-(g/y+i/y+kappa/2*x^2*n/y))*y;

Gamma = 1;

yc = c/y;

yi = i/y;

ynu = rk*k/y;

yx = (kappa/2)*(x^2*n/y);

htilda = h/gammaz;

etanu = (1-psinu)/psinu; 

betatilda = beta/gammaz;

kappaa = (kappa*x)^(-1)*pw*a;

kappaw = (kappa*x)^(-1)*w;

kappalambda = beta*(1+rho)/2;

phia = chi*pw*a*w^(-1);

psi = chi*beta*lambda*mu+(1-chi)*rho*beta*lambda*e;

tau = psi/(1+psi);

% i take the initiative for finding H_upperbar, from equations 
% 41 and B3 of the paper

H = chi/(1-chi)*kappa*x;

phis = (1-chi)*s*beta*H*w^(-1);

phix = chi*beta*kappa*x^2*w^(-1);

phib = (1-chi)*b*w^(-1);

phichi = chi*(1-chi)^(-1)*kappa*x*w^(-1);

phieta = phichi*(1-chi)*(1-eta)*y^(-1);

taup = 1+(epsilonp-1)*csi;

smallcp = (1-lambdap)*(1-lambdap*beta)*(lambdap)^(-1);

% what is steady state inflation?
pi = 1;

% should be correct, but no upperbar in appendix
gammap = pi;

fip = 1+beta*gammap;

iotaf = beta*(fip)^(-1);

iotao = (smallcp/taup)*(fip)^(-1);

iotab = gammap*(fip)^(-1);

G = (1-eta*x*beta*lambda*mu)*eta^(-1)*mu*kappaw;

tau1 = (kappaw*mu*phix+phichi*(1-chi)*(x*beta*lambda)*(kappaw*mu)*mu*(rho*beta)+phis*G)*(1-tau);

tau2 = -(kappaw*mu)*phichi*(1-chi)*(x*beta*lambda)*mu*(1-tau);

smallc = (1-lambda)*(1-tau)*lambda^(-1);

fi = (1+tau2)+smallc+(tau*lambda^(-1)-tau1);

gammab = (1+tau2)*fi^(-1);

gammao = smallc*fi^(-1);

gammaf = (tau*lambda^(-1)-tau1)*fi^(-1);



yt=0; 
kt=0; 
nt=0; 
ct=0; 
it=0; 
gt=0; 
nut=0; 
xt=0; 
mt=0; 
ut=0; 
vt=0; 
qt =0;
st =0;
kpt =0;
et=0;
Gammat=0; 
rt =0;
pit =0;
rkt =0;
qkt =0;
pwt =0;
at =0;
wt =0;
chit =0;
mut =0;
wot=0;
bt =0;
fyt =0;
thetat =0;
e_wt =0;
e_bt =0;
e_it =0;
e_zt =0;
e_etat =0;
e_pt =0;
e_rt=0;
e_gt =0;
dy=100*(gammaz-1);
dinve=100*(gammaz-1);
dc =100*(gammaz-1);
dw =100*(gammaz-1);
labobs =0;
pinfobs =0;
robs=0;
AUX_ENDO_LAG_3_1=0;





for iter = 1:length(M_.params)
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

if isfield(M_,'param_nbr') == 1

if isfield(M_,'orig_endo_nbr') == 1
NumberOfEndogenousVariables = M_.orig_endo_nbr;
else
NumberOfEndogenousVariables = M_.endo_nbr;
end
ys = zeros(NumberOfEndogenousVariables,1);
for i = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(i,:));
  eval(['ys(' int2str(i) ') = ' varname ';']);
end
else
ys=zeros(length(lgy_),1);
for i = 1:length(lgy_)
    ys(i) = eval(lgy_(i,:));
end
check = 0;
end
end