function [ys,params,check]=model_estimation_steadystate(ys,exo,M_,options_)
check = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set parameter values
  alp=M_.params(1);
  alpu=M_.params(2);
  bet=M_.params(3);
  deltak=M_.params(4);
  B0=M_.params(5);
  B1=M_.params(6);
  etac=M_.params(7);
  etal=M_.params(8);
  etals=M_.params(9);
  se=M_.params(10);  
  st=M_.params(11);  
  TAU_C=M_.params(12);
  TAU_H=M_.params(13);
  rs1=M_.params(14);
  rs2=M_.params(15);
  rs3=M_.params(16);
  pe=M_.params(17);
  pt=M_.params(18);

  chin=1;
  chii=1;
  tauh=TAU_H;
  tauc=TAU_C;
  A=1;
  B=1;

  R=(1/bet+deltak-1)/(1-tauh);
  Ru=R*(1-tauh);
  Hm_Km=(R/(alp*(1-tauc)))^(1/(1-alp));
  Hu_Ku=(Ru/(alpu*(1-tauc*st*pt)))^(1/(1-alpu));
  Km_Hm=1/Hm_Km;
  Ku_Hu=1/Hu_Ku;

  
  x=[0.2500    0.0555]';
  opt=optimset('Display','off');
  X=fsolve(@(x)obs_ss(x,alp, alpu, bet, pe, se, pt, st, deltak, etac, etal, etals, rs1, rs2, rs3, B1, B0, tauh, tauc)',x',opt);
  H=X(1);
  Hu=X(2);

  Hm=H-Hu;
  Km=Km_Hm*Hm;
  Ku=Ku_Hu*Hu;
  K=Ku+Km;
  Ym_Hm=A*(Km_Hm)^(alp);
  Yu_Hu=B*(Ku_Hu)^(alpu);
  Wm=(1-alp)*Ym_Hm/(1+rs1+rs2);
  Wu=((1-alpu)*(1-tauc*st*pt)*Yu_Hu)/(1+se*pe*(rs1+rs2+rs3)-rs3-tauc*(1+rs1+rs2)*st*pt);
  Iv=deltak*K;
  C=(1-tauh)*(Wm*Hm+R*Km)-Iv+Wu*Hu+Ru*Ku-rs3*(Wm*Hm+Wu*Hu);
  B0=((1-tauh-rs3)/H^etal)*C^(-etac)*Wm;
  B1=(-B0*(H)^etal+Wu*(1-rs3)*C^(-etac))/(Hu^etals);

  Ym=Ym_Hm*Hm;
  Yu=Yu_Hu*Hu;
  Y=Ym+Yu;


  gc=tauc*(Ym-Wm*Hm-rs1*Wm*Hm-rs2*Wm*Hm)+tauc*st*pt*(Yu-Wu*Hu-rs1*Wu*Hu-rs2*Wu*Hu);             
  gs=(rs1+rs2+rs3)*(Wm*Hm+se*pe*Wu*Hu);                                                         
  gh=tauh*(Wm*Hm+R*Km);                                                                           
  G=gc+gs+gh;
  Y=C+Iv+G; 
  evs=(1-pe)*Wu*Hu*(rs1+rs2+rs3)+(1-pt)*tauc*(Yu-Wu*Hu-(rs1+rs2)*Wu*Hu)+tauh*(Wu*Hu+Ru*Ku);   
  evs_per=evs/(evs+G)*100;
  ysy=Yu/Y*100;
  ysym=Yu/Ym*100;
  Ym_obs=0;
  C_obs =0;
  Hm_obs=0;
  Iv_obs=0;
  gh_obs=0;
  gc_obs=0;
  gs_obs=0;
 
ys=[Y         
  Ym        
  Yu
  ysy  
  K         
  Km        
  Ku        
  H         
  Hm        
  Hu 
  Wm        
  Wu        
  R         
  Ru        
  C
  Iv
  G
  gs        
  gh        
  gc
  pe
  pt        
  chii     
  chin
  log(A)
  log(B)
  tauc      
  tauh                
  evs      
  evs_per   
  C_obs       
  Ym_obs 
  Hm_obs
  Iv_obs      
  gs_obs       
  gh_obs       
  gc_obs
  ysym]; 

params=M_.params;
