% function [u_bar, b_bar, g_bar, P] = NEW_HANK_JL_stationary(theta, ys, exo, M_, options_)
function [u_bar, b_bar, g_bar, P] = NEW_HANK_JL_stationary(M_)

% Compute stationary distribution (ū, b̄, ḡ) given theta
% Uses the 3x3 transition matrix Π_t in your notes, with:
%   Φ  = θ^(1-γ),  Φ_g = ϱ Φ,  Φ_b = (1-ϱ) Φ,
% and OJS intensity ζ on surviving bad.

    % --- read parameters by name ---
    for ii = 1:M_.param_nbr
        eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');' ]);
    end

    % Get Theta
    [~, theta] = LM_eq(M_) ;

    % aliases (match your .mod names)
    rho   = varrho;       % share of good-job offers
    sigb  = sig_b;        % separation: bad
    sigg  = sig_g;        % separation: good
    gamma_= gama;       % matching elasticity

    % --- success rates ---
    Phi   = theta^(1 - gamma_);
    Phi_g = rho * Phi;
    Phi_b = (1 - rho) * Phi;

    % --- transition matrix Π (rows: current state u,b,g; cols: next state u,b,g) ---
    % Row u:     [1-Φ,      Φ_b,                              Φ_g]
    % Row b:     [σ_b(1-Φ), σ_b Φ_b + (1-σ_b)(1 - ζ Φ_g),     σ_b Φ_g + (1-σ_b) ζ Φ_g]
    % Row g:     [σ_g(1-Φ), σ_g Φ_b,                          σ_g Φ_g + (1-σ_g)]
    P = [ ...
        1 - Phi,                      Phi_b,                                   Phi_g; 
        sigb*(1 - Phi),               sigb*Phi_b + (1 - sigb)*(1 - zeta*Phi_g), sigb*Phi_g + (1 - sigb)*(zeta*Phi_g);
        sigg*(1 - Phi),               sigg*Phi_b,                              sigg*Phi_g + (1 - sigg) ...
    ];

    % (optional) small numerical clean-up
    P = max(P, 0);                        % clip tiny negatives
    P = P ./ sum(P,2);                    % row-normalize to sum exactly to 1

    % --- stationary distribution π solves π = π P, sum(π)=1 ---
    A = [P.' - eye(3); 1 1 1];            % 4x3 system (last row is normalization)
    b = [0;0;0;1];
    pi = A \ b;                           % least-squares solve

    % ensure real/normalized
    pi = real(pi); 
    pi = pi / sum(pi);

    % return in order (u,b,g)
    u_bar = pi(1);
    b_bar = pi(2);
    g_bar = pi(3);
end