function [ys,check] = main201808_steadystate(ys,junk,do_steady_print)

global M_

% read parameters values from params.m
[pi,rho, u0, h0, eta, alpha, delta_A, delta_phi, B, g, theta, beta, sigma, ...
Astar, phistar]=params;
% If some parameters have already been estimated in main201808.mod, then ...
% employ it instead of using values from params.m 
% see http://www.dynare.org/manual/index_18.html
if isnan(M_.params(:,:))==0
  pi=M_.params(1); 
  rho=M_.params(2); 
  u0=M_.params(3); 
  h0=M_.params(4); 
  eta=M_.params(5); 
  alpha=M_.params(6); 
  delta_A=M_.params(7); 
  delta_phi=M_.params(8); 
  B=M_.params(9); % share of liquid assets target pb*B/(pb*B+q*K)=0.12
  g=M_.params(10); 
  theta=M_.params(11); 
  beta=M_.params(12); 
  sigma=M_.params(13);
  Astar=M_.params(14); 
  phistar=M_.params(15);
  end  
  
  % values of endogenous variables at steady state
  % se i ce l c s b A phi q pb tau r w lamdae K
  Ass=Astar;
  %bss=B;
  lss=0.25/(1-pi); % total hours of work target=0.25
  % output
  %Fss=g/0.18; % govt spending/GDP target=0.18
  rss=alpha/(3.32*4); % ratio of capital to annual output target=3.32=K/(4*Y)
  Kss=(rss/alpha)^(1/(alpha-1))*lss*(1-pi);
  sss=Kss;
  Fss=Ass*(Kss)^(alpha)*((1-pi)*lss)^(1-alpha);
  % modify parameter g
  g=Fss*0.18; % govt spending/GDP target=0.18
  %rss=alpha*(Kss/(lss*(1-pi)))^(alpha-1);
  wss=(1-alpha)*(Kss/(lss*(1-pi)))^alpha;
  iss=(1-sigma)*Kss/pi;
  pbss=(1/1.02)^(0.25); % annualized net rate of returen on liquid assests to the target yields 1/(pb)^4-1=0.02
  % Set theta equal to phistar
  lamdaess=(pbss/beta-1)/pi; % bond pricing equation
  %qss=(lamdaess+1)/(1+theta*lamdaess); % entrepreneur's financing constraint is binding
  % Assume that theta=phistar and by equation (12) and (13) in Shi, 2015, JME
  qss=beta*(rss+pi*lamdaess*rss+sigma*pi*(1+lamdaess))/(1-beta*sigma+beta*sigma*pi);
  phiss=(lamdaess-qss+1)/(lamdaess*qss);
  % Modify parameter theta
  theta=phiss;
  % Modify parameter B
  B=0.12*qss*Kss/(0.88*pbss); %Target share of liquid assets pb*B/(pb*B+q*K)=0.12
  tauss=g-(pbss-1)*B; % govt budget constraint
  bss=B;
  cess=(rss+phiss*sigma*qss)*sss+bss-tauss-(1-theta*qss)*iss; % entrepreneur financing constraint equation (4) in the paper
  css=Fss-pi*iss-g;
  sess=(1-theta)*iss+(1-phiss)*sigma*sss;
  Yss=Fss;
  
  % modify several parameters values
  u0=((css-pi*cess)/((1-pi)*cess))^(-rho)*(1+lamdaess); 
  h0=1/eta*(lss)^(1-eta)*wss*((css-pi*cess)/(1-pi))^(-rho);
  phistar=phiss;
 
  % logarithmize steady state values
  Assl=log(Ass); bssl=log(bss); lssl=log(lss);
  Kssl=log(Kss); sssl=log(sss); rssl=log(rss); wssl=log(wss);
  issl=log(iss); pbssl=log(pbss); phissl=log(phiss); lamdaessl=log(lamdaess);
  qssl=log(qss); taussl=log(tauss); cessl=log(cess); cssl=log(css);
  sessl=log(sess); Yssl=log(Yss);
  
  % put steady state value into ys vector
  % se i ce l c s b A phi q pb tau r w lamdae K
  ys=NaN*ones(M_.orig_endo_nbr,1);
ys(                   1)=sessl                ;
ys(                   2)=issl               ;
ys(                   3)=cessl                 ;
ys(                   4)=lssl                ;
ys(                   5)=cssl             ;
ys(                   6)=sssl                ;
%ys(                   7)=bssl                ;
ys(                   7)=Assl               ;
ys(                   8)=phissl                ;
ys(                   9)=qssl          ;
ys(                  10)=pbssl                ;
ys(                  11)=taussl                 ;
ys(                  12)=rssl               ;
ys(                  13)=wssl             ;
ys(                  14)=lamdaessl                ;
ys(                  15)=Yssl                 ;
%ys(                  17)=Yssl;


% saving s.s. values that are used as params in static file
% pi,rho, u0, h0, eta, alpha, delta_A, delta_phi, B, g, theta, beta, sigma, ...
% Astar, phistar
M_.params(                   1)=pi;
M_.params(                   2)=rho;
M_.params(                   3)=u0            ;
M_.params(                   4)=h0            ;
M_.params(                   5)=eta              ;
M_.params(                   6)=alpha              ;
M_.params(                   7)=delta_A          ;
M_.params(                   8)=delta_phi          ;
M_.params(                   9)=B          ;
M_.params(                  10)=g          ;
M_.params(                  11)=theta           ;
M_.params(                  12)=beta             ;
M_.params(                  13)=sigma           ;
M_.params(                  14)=Astar          ;
M_.params(                  15)=phistar             ;  
  
  
check=0;

if iss<0,
    disp('Steady state i<0');
    check=1;
end

if lamdaess<0,
    disp('steady state lamdaess<0');
    check=1;
end
 

  
  
  
  
  
  
  
  
  
  