function [ys,check] = FINAL_HDS_steadystate(junk,ys)


    %global M_


    j = get_param_by_name('j') ;
    tau = get_param_by_name('tau') ;
    alppha = get_param_by_name('alppha') ;
    mu = get_param_by_name('mu') ;
    nu = get_param_by_name('nu') ;
    sigma1 = get_param_by_name('sigma1') ;
    sigma2 = get_param_by_name('sigma2') ;
    del_k = get_param_by_name('del_k') ;
    
    rho = get_param_by_name('rho') ;
    pRpi = get_param_by_name('pRpi') ;
    pRy = get_param_by_name('pRy') ;
    pRq = get_param_by_name('pRq') ;
    %pRB = get_param_by_name('pRB') ;
    
    pMHy = get_param_by_name('pMHy') ;
    pMHq = get_param_by_name('pMHq') ;
    pMHB = get_param_by_name('pMHB');
    epsilon = get_param_by_name('epsilon') ;
    theta = get_param_by_name('theta') ;
    MK = get_param_by_name('MK') ;
    
    beta_p = get_param_by_name('beta_p') ;
    beta_i = get_param_by_name('beta_i') ;
    beta_r = get_param_by_name('beta_r') ;
    beta_e = get_param_by_name('beta_e') ;
    Xss = get_param_by_name('Xss') ;
    Yss = get_param_by_name('Yss') ;
    Gss = get_param_by_name('Gss') ;
    qss = get_param_by_name('qss') ;
    Bss = get_param_by_name('Bss') ; 
    MHss = get_param_by_name('MHss') ;
    Rss = get_param_by_name('Rss') ;
    tau_pss = get_param_by_name('tau_pss') ;
    
    rho_tau = get_param_by_name('rho_tau') ;
    del_y = get_param_by_name('del_y') ;
    
    tau_y = get_param_by_name('tau_y') ;
    I_r = get_param_by_name('I_r') ;
    I_p = get_param_by_name('I_p') ; 
    
    rho_az = get_param_by_name('rho_az') ;
    rho_aj = get_param_by_name('rho_aj') ;
    sigma_Z = get_param_by_name('sigma_Z') ;
    sigma_J = get_param_by_name('sigma_J') ;
    sigma_R = get_param_by_name('sigma_R') ;
    sigma_tau = get_param_by_name('sigma_tau') ;
    PIss = get_param_by_name('PIss') ;



    check = 0;


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                           STEADY STATE                               %                                                
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%% USEFUL CONSTANTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    R = 1/beta_p;
    rk = (1-beta_p)/beta_p + del_k;
    lambda1 = 1 - beta_i*R;
    lambda2 = 1 - beta_e*R;
    
    X = Xss;   % makeup equals to Xss which is derived from calibration or Iacoviello(2005)  
    PI = 1;    % P_t=P_{t-1} for all time  
    PI_x = 1;  % the ss value of PI_x is derived from equation 32 and PI=1  
    MH = MHss;
    tau_p = tau_pss;
    
    Q = beta_p*(1-tau_y)/(1-beta_p*(1-tau_p*(1-I_p*tau_y)) );
    Rk = (1-lambda2*MK)/beta_e - (1-del_k);
    Rv = (1-(lambda2*MH/R))/beta_e - (1-tau_p*(1-I_p*tau_y)) + (I_r*tau_y)/Q;

    WAGE = (1-alppha-nu)/X; 
     KEP = alppha*mu/(Rk*X);
     KPP = alppha*(1-mu)/(rk*X); 
    
    FQHR = j*Q;
    FCRY = (WAGE*(1-sigma1-sigma2))/(1+j);
    
    CIBI = (1-R)*MH/R;
    CIHI = I_r*tau_y/Q + (tau_p*(1-I_p*tau_y));
     QHI = 1 - lambda1*MH/R - beta_i*(1-tau_p*(1-I_p*tau_y)) + beta_i*I_r*tau_y/Q;
    FQHI = (beta_i*j)/QHI;
    FCIY = (sigma2*WAGE)/(1 - (CIBI - CIHI)*FQHI);

    CPHR = (1-tau_y)/Q - (tau_p*(1-I_p*tau_y));
     QHP = 1 - beta_p*(1-tau_p*(1-I_p*tau_y)) + beta_p*I_r*tau_y/Q; 
    FQHP = (beta_p*j)/QHP;

     QHE = 1 - lambda2*MH/R - beta_e*(1-tau_p*(1-I_p*tau_y)) + beta_e*I_r*tau_y/Q;
    FQHE = beta_e*nu/(X*QHE);

    CPYN = sigma1*WAGE + (X-1)/X - CIBI * FQHI * FCIY - CIBI * FQHE - (1-R)*MK*KEP + ((1-beta_p)/beta_p)* KPP + CPHR * FQHR * FCRY;
    FCPY = CPYN/(1 + CIHI*FQHP );

   % FCEY = (1/X - WAGE - rk*KPP + ((1-R)*MK - del_k)*KEP + (CIHI+CIBI)*FQHE);
    FCEY = 1/X  + (CIBI*FQHE + (1-R)*MK*KEP) - WAGE -del_k*KEP -rk*KPP -CIHI*FQHE;
    QHYP = FQHP * FCPY + FQHI * FCIY + FQHR * FCRY + FQHE;

    PNPN = sigma1*WAGE;
    PNPD = FCPY;
    PNIN = sigma2*WAGE;
    PNID = FCIY;
    PNRN = (1-sigma1-sigma2)*WAGE;
    PNRD = FCRY;

    %%%%%%%%%%%%%%%%%%%%%%%%% STEADY STATE VALUE %%%%%%%%%%%%%%%%%%%%%%%%%%
    % start with those variables which are consist of purely parameters %
    N_p = PNPN/(tau*PNPD + PNPN);
    N_i = PNIN/(tau*PNID + PNIN);
    N_r = PNRN/(tau*PNRD + PNRN);

    h_p = FQHP * FCPY / QHYP;
    h_i = FQHI * FCIY / QHYP;
    h_r = FQHR * FCRY / QHYP;
    h_e = FQHE / QHYP;
    
    Y = KPP^(alppha*(1-mu)/(1-alppha)) * KEP^(alppha*mu/(1-alppha)) * h_e^(nu/(1-alppha)) * N_p^(sigma1*(1-alppha-nu)/(1-alppha)) * N_i^(sigma2*(1-alppha-nu)/(1-alppha)) * N_r^((1-sigma1-sigma2)*(1-alppha-nu)/(1-alppha));
    k_p = KPP*Y;
    k_e = KEP*Y;
    b_i = MH*FQHI*FCIY*Y/R;
    b_e = MH*FQHE*Y/R + MK*KEP*Y;
    B = b_i + b_e; 
    w_p = sigma1*WAGE*Y/N_p;
    w_i = sigma2*WAGE*Y/N_i;
    w_r = (1-sigma1-sigma2)*WAGE*Y/N_r;

    q = QHYP*Y;  % The results under line 113 or line 114 are the same 
    % q = nu*Y/(Rv*h_e);
    rh = q/Q;

    c_p = FCPY*Y;
    c_i = FCIY*Y;
    c_r = FCRY*Y;
    c_e = FCEY*Y;

    F = ((X-1)/X)*Y;
    tax = (tau_p*(1-I_p*tau_y))*q*1 + tau_y*rh*(h_r + I_r*(h_p + h_i + h_e));
    G = tax;
     
    C = c_p + c_i + c_r + c_e;
    B = b_i + b_e;
    I = del_k*(k_p + k_e);
    N = N_p + N_i + N_r;
    
    Ygap = Y/Yss;
    qgap = q/qss;
    PIgap = PI/PIss;
    
    
    ucp = 1/c_p;
    uci = 1/c_i;
    ucr = 1/c_r;
    uce = 1/c_e;

    x1 = Y/(1-theta*beta_p);
    x2 = Y/(1-theta*beta_p);
    
    V_p = (log(c_p) + j*log(h_p) + tau*log(1-N_p))/(1-beta_p);
    V_i = (log(c_i) + j*log(h_i) + tau*log(1-N_i))/(1-beta_i);
    V_r = (log(c_r) + j*log(h_r) + tau*log(1-N_r))/(1-beta_r);
    V_e = (log(c_e))/(1-beta_e);
    
    welfare = (1-beta_p)*V_p + (1-beta_i)*V_i + (1-beta_r)*V_r + (1-beta_e)*V_e;
   
    A_z = 0; 
    A_j = 0;

    
    xxx = [ ...
    c_p   
    c_i   
    c_r   
    c_e   

    h_p   
    h_i   
    h_r   
    h_e  
    
    N_p    
    N_i    
    N_r   

    k_p    
    k_e   

    b_i    
    b_e  

    
    Y      
    F      
    G      
    tax    
    tau_p

    w_p    
    w_i    
    w_r    

    q     
    rh      
    rk     
    Rk   
    Rv

    R      
    X      
    PI     
    PI_x   
    MH     

    A_z  
    A_j


    x1    
    x2    

    ucp   
    uci   
    ucr   
    uce   
    
    V_p
    V_i
    V_r
    V_e
    welfare
    
    lambda1   
    lambda2
    
    C
    I
    N
    B
    
    Ygap
    qgap
    PIgap;] ;



ys = (xxx) ;