function [ys, params, check] = modelo_final_steadystate(ys, exo, M_, options_)
% Steady-state file for modelo_final.mod — NK model with search and vacancies (consistent and robust)
% Author: César M. Martínez 
% Robust version: scalar theta + fixed-point on mc + fixed-point on Rk

check = 0;

%% === 1) Load parameters ===
p = struct();
for ii = 1:M_.param_nbr
    pname = strtrim(M_.param_names{ii});
    p.(pname) = M_.params(ii);
end

% Aliases (for readability)
beta    = p.beta;     delta   = p.delta;    alpha   = p.alpha;
eps     = p.eps;      mu_p    = p.mu_p;     rho     = p.rho;
phi     = p.phi;      mu      = p.mu;       omega   = p.omega;
b       = p.b;        kappa_h = p.kappa_h;  Abar    = p.Abar;
Xbar    = p.Xbar;     Sbar    = p.Sbar;     sigma   = p.sigma; %#ok<NASGU>
eta     = p.eta;      varphi  = p.varphi;   kappa   = p.kappa;
phi_pi  = p.phi_pi;

%% === 2) Exogenous variables at steady state ===
A = Abar;   % productivity
X = Xbar;   % exogenous cost
S = Sbar;   % subsidies

%% === 3) Numerical helper functions ===
clippos = @(x) max(real(x),1e-12);

% Y/n given (mc, Rk)
yell_of_mc = @(mc,Rk) ( A * ((mc*alpha)/Rk).^alpha ).^(1/(1-alpha));

% Nash wage w(mc,theta,Rk)
w_of_mc_theta = @(mc,th,Rk) (1-omega)*b + omega*((1-alpha).*mc.*yell_of_mc(mc,Rk) - kappa_h.*th);

% Z(mc,Rk) = (1-α)mc*yl
Z_of_mc = @(mc,Rk) (1-alpha).*mc.*yell_of_mc(mc,Rk);

% Fixed point on mc (for a given theta and Rk), solving: mc = [ w(mc,θ)^(1-α) * Rk^α ] / A
function [mc, ok] = mc_of_theta(theta, Rk)
    mc0   = max(1/mu_p, 1e-3);
    damp  = 0.7;
    tol   = 1e-12;
    itmax = 500;
    mc    = mc0;
    ok    = false;
    for it = 1:itmax
        w  = w_of_mc_theta(mc, theta, Rk);
        if ~isfinite(w) || w<=0
            mc = NaN; return;
        end
        mc_new = (w.^(1-alpha) .* Rk.^alpha) / A;
        if ~isfinite(mc_new) || mc_new<=0
            mc = NaN; return;
        end
        if abs(mc_new - mc) < tol
            mc = clippos(mc_new);
            ok = true; return;
        end
        mc = clippos(damp*mc_new + (1-damp)*mc);
    end
    % If it did not converge:
    ok = false; mc = NaN;
end

% Free-entry residual function for root in theta (with mc(theta))
% F_entry(theta;Rk) = (Z(mc)-w(mc)) / [1-(1-ρ)β] - κ_h θ^μ / φ
function [F, mc, w, ok] = entry_residual(theta, Rk)
    [mc, ok] = mc_of_theta(theta, Rk);
    if ~ok || ~isfinite(mc) || mc<=0
        F = NaN; w = NaN; return;
    end
    w   = w_of_mc_theta(mc, theta, Rk);
    Z   = Z_of_mc(mc, Rk);
    lhs = (Z - w) / (1 - (1 - rho)*beta);
    rhs = kappa_h * theta^mu / phi;
    F   = lhs - rhs;
end

% Solve for theta satisfying F_entry(theta)=0 (robust: bracketing + fallback)
function [theta, mc, w, ok] = solve_theta(Rk)
    ok = false;
    % 1) Broad grid to detect sign change
    th_grid = logspace(-6, 2, 200);  % from 1e-6 to 1e2
    Fvals = NaN(size(th_grid));
    mcvals = NaN(size(th_grid));
    wvals  = NaN(size(th_grid));
    for i=1:numel(th_grid)
        [Fvals(i), mcvals(i), wvals(i), ok_i] = entry_residual(th_grid(i), Rk);
        if ~ok_i || ~isfinite(Fvals(i)), Fvals(i)=NaN; end
    end
    % 2) Try bracketing
    idx = find(isfinite(Fvals));
    if numel(idx)>=2
        sgn = sign(Fvals(idx));
        chg = find(diff(sgn)~=0);
        if ~isempty(chg)
            % Found a bracket [a,b]
            a = th_grid(idx(chg(1)));
            b = th_grid(idx(chg(1)+1));
            % Use fzero with bracket
            try
                theta = fzero(@(th) entry_residual(th,Rk), [a,b]);
                [~, mc, w, ok] = entry_residual(theta, Rk);
                if ok && isfinite(theta) && theta>0
                    return;
                end
            catch
                % Go to fallback
            end
        end
    end
    % 3) Fallback: minimize |F| using fminbnd in a reasonable range
    a = 1e-6; b = 1e2;
    absF = @(th) abs(entry_residual(th,Rk));
    try
        theta = fminbnd(absF, a, b);
        [Fbest, mc, w, ok] = entry_residual(theta, Rk);
        if ok && isfinite(Fbest) && abs(Fbest) < 1e-8 && theta>0
            ok = true; return;
        end
    catch
    end
    % 4) Last resort: best grid point
    if ~ok && ~isempty(idx)
        [~, j] = min(abs(Fvals(idx)));
        theta = th_grid(idx(j));
        mc    = mcvals(idx(j));
        w     = wvals(idx(j));
        ok    = isfinite(theta) && isfinite(mc) && isfinite(w);
    end
    % DEBUG: inspect sign of F_entry over grid
    th_grid = logspace(-6, 2, 200);
    Fvals = NaN(size(th_grid));
    for i=1:numel(th_grid)
        Fvals(i) = entry_residual(th_grid(i), Rk);
    end
    %fprintf('[DEBUG] F_entry(min)=%.4e, F_entry(max)=%.4e\n', min(Fvals), max(Fvals));
end

%% === 4) External fixed-point on Rk (since pi depends on mc via NKPC) ===
Rk = 1/beta - 1 + delta;      % initial value with pi=1
tol_Rk  = 1e-10;
maxitPF = 100;
damp_Rk = 0.7;
theta = 0.5; mc = 1/mu_p; w = (1-omega)*b;  % soft initial guesses

for it = 1:maxitPF
    % Given Rk, solve for theta and mc
    [theta_new, mc_new, w_new, ok_theta] = solve_theta(Rk);
    if ~ok_theta || ~isfinite(mc_new) || mc_new<=0 || ~isfinite(theta_new) || theta_new<=0
        error('Steady state could not be closed: try different initial values or parameters.');
    end
    theta = theta_new; mc = mc_new; w = w_new;

    % NKPC (in levels): pi given mc
    gap_mc = log(mc) - log(1/mu_p);
    pi = exp( (kappa/(1 - beta)) * gap_mc );
    if ~isfinite(pi) || pi<=0
        error('NKPC produced non-finite or non-positive pi. Check calibration.');
    end

    % Capital Euler equation → Rk consistent with pi
    Rk_new = (pi/beta) - 1 + delta;

    % Convergence in Rk
    if abs(Rk_new - Rk) < tol_Rk
        Rk = Rk_new;
        break
    end
    % Damping
    Rk = damp_Rk*Rk_new + (1-damp_Rk)*Rk;
end
if it==maxitPF && abs(Rk_new - Rk) >= tol_Rk
    error('Steady state could not be closed: fixed-point iteration in Rk did not converge.');
end

%% === 5) Derived objects with final (mc,theta,Rk,pi) ===
yl    = yell_of_mc(mc,Rk);     % Y/n
Z     = Z_of_mc(mc,Rk);

% Labor market: ρ n = φ (1-n) θ^{1-μ}
n   = (phi * theta.^(1-mu)) / (rho + phi * theta.^(1-mu));
n   = min(max(n, 1e-9), 1-1e-9);
u   = 1 - n;
m   = rho * n;
f   = m / u;
v   = clippos(theta * u);
q   = phi * theta^(-mu);
J   = kappa_h / q;

% Nominal and NK variables
R     = pi/beta;
pibar = 1;
Rbar_new = R / (pi/pibar)^phi_pi;

% Update Rbar in M_.params
idx_Rbar = find(strcmp(strtrim(M_.param_names), 'Rbar'));
if ~isempty(idx_Rbar)
    M_.params(idx_Rbar) = Rbar_new;
end
params_vec = M_.params;
for ii = 1:M_.param_nbr
    if strcmp(strtrim(M_.param_names{ii}), 'Rbar')
        params_vec(ii) = Rbar_new;
        break
    end
end

%% === 6) Production, stocks, and allocation ===
KY  = (mc*alpha)/Rk;   % K/Y from capital FOC
Y   = yl * n;
K   = KY * Y;
I   = delta * K;
C   = Y - I - Sbar;

% Profits
Pi  = Y - mc*Y - X + S;

%% === 7) Pack endogenous variables ===
ys = zeros(M_.endo_nbr,1);
for ii = 1:M_.endo_nbr
    vname = strtrim(M_.endo_names{ii});
    try
        ys(ii) = eval(vname);
    catch
        ys(ii) = 0;
    end
end

% Diagnostic log
fprintf('[SS] mc=%.6f, theta=%.6f, w=%.6f, n=%.6f, pi=%.6f, R=%.6f, Rk=%.6f\n', ...
        mc, theta, w, n, pi, R, Rk);

%% === 8) Outputs ===
params = real(params_vec);
ys     = real(ys);

if all(isfinite(ys)) && all(isfinite(params)) && all(abs(imag(ys))==0)
    check = 0;
else
    disp('⚠️  Steady state has non-finite or complex values, setting check=1');
    check = 1;
end

ys     = double(ys(:));
params = double(params(:));
return
end
