function [ys,params,check] = NEW_HANK_JL_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = Gali_2010_steadystate(ys,exo,M_,options_)
% uses a numerical solver to get the steady state for the nonlinear version of Gali_2010_steadystate.mod 
% in order to set the required parameters matching the calibration
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% Copyright (C) 2025 Henri Graul and Matteo Samarani

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%% Caliobration of parameters

%% Step 1: Solve for tightness and implied stationary distribution

% 1a) Solve θ from supply-side aggregator (returns Φ, ζ_g, ζ_b)

[Phi_bar, theta_bar] = LM_eq(M_) ;
zeta_g_bar           = ( kappa*sig_g*((1-betta) + betta*sig_g/2) + eta*alpa + (1-eta)*Tu) / alpa;
zeta_b_bar           = zeta_b(M_, Phi_bar) ;


% Convenience aliases for later steps
theta   = theta_bar;          % labor-market tightness
Phi     = Phi_bar;            % θ^(1-γ)
Phi_g   = varrho * Phi_bar;       % success toward good jobs
Phi_b   = (1-varrho) * Phi_bar;   % success toward bad jobs
zeta_g_ss = zeta_g_bar;       % type-g SS cost term
zeta_b_ss = zeta_b_bar;       % type-b SS cost term

% 1b) Stationary distribution of (u,b,g) implied by θ
[u_bar, b_bar, g_bar, Pbar] = NEW_HANK_JL_stationary(M_);

% 1c) Immediate SS objects used in later calibrations
s_bar  = u_bar + sig_b*b_bar + sig_g*g_bar + zeta*(1 - sig_b)*b_bar;  % search pool
v_bar  = theta * s_bar;                                                      % vacancies
eta_g_bar = sig_g;                                                           % hires rate, good
eta_b_bar = sig_b + zeta*varrho*(1 - sig_b)*Phi;                       % hires rate, bad
rho_g_bar = 1 - sig_g;                                                       % retention, good
rho_b_bar = (1 - sig_b)*(1 - zeta*varrho*Phi);                         % retention, bad


%% Step 2: Back out all SS-calibrated parameters

% --- handy eps for safe division/clipping
tiny = 1e-12; clip01 = @(x) min(max(x,0),1);

% ---------- (a) MC aggregator weight ----------
if abs(varpi - 1) < 1e-8
    sbar_g = varrho;                         % Cobb–Douglas limit
else
    sbar_g = ( varrho * zeta_g_ss^(1 - varpi) ) ...
           / ( varrho * zeta_g_ss^(1 - varpi) + (1 - varrho) * zeta_b_ss^(1 - varpi) );
end
sbar_g   = clip01(sbar_g);
omegaA_g = sbar_g;
omegaA_b = 1 - omegaA_g;

% ---------- (b) Vacancy shares ----------
Hg = eta_g_bar * g_bar;                      % hires g
Hb = eta_b_bar * b_bar;                      % hires b
omegav_g = clip01( Hg / max(Hg + Hb, tiny) );  % common fill-rate baseline
omegav_b = 1 - omegav_g;

% ---------- (c) Hiring-cost weights (for dC^H) ----------
num_g = (eta_g_bar^2) * g_bar;
num_b = (eta_b_bar^2) * b_bar;
denCH = max(num_g + num_b, tiny);
omegaC_g = clip01( num_g / denCH );
omegaC_b = 1 - omegaC_g;

% ---------- (d) Search-input weights ----------
omegas_u = clip01( u_bar / max(s_bar, tiny) );
omegas_b = clip01( (sig_b + zeta*(1 - sig_b)) * b_bar / max(s_bar, tiny) );
omegas_g = clip01( sig_g * g_bar / max(s_bar, tiny) );

% ---------- (e) Stock law-of-motion weights ----------
omegap_eta_g = sig_g;                        % = \bar H_g / \bar g
omegap_eta_b = 1-(1-sig_b)*(1-zeta*varrho*Phi);                        % = \bar H_b / \bar b
omegap_rho_b = (1-sig_b)*(1-zeta*varrho*Phi);

% ---------- (f) Acceptance & retention sensitivities ----------
psi_b = clip01( ( zeta*(1 - sig_b) * (b_bar / max(s_bar,tiny)) ) ...
               / max( 1 - zeta*(1 - sig_b) * (b_bar / max(s_bar,tiny)), tiny ) );

xi_b  = clip01( ( zeta * varrho * Phi ) ...
               / max( 1 - zeta * varrho * Phi, tiny ) );

% ---------- (g) Type-level MC linearization coefficients ----------
% wages (levels) with alpha_bar = 1
w_g_bar = eta*phi_g + (1 - eta)*Tu;
w_b_bar = eta*phi_b + (1 - eta)*Tu;

% M_s and N_s
M_g = kappa * eta_g_bar * ( rho_g_bar + 0.5*eta_g_bar );
M_b = kappa * eta_b_bar * ( rho_b_bar + 0.5*eta_b_bar );

N_g = kappa*eta_g_bar + w_g_bar - betta*M_g;
N_b = kappa*eta_b_bar + w_b_bar - betta*M_b;

% lambdas
lambda_eta_g = (kappa*eta_g_bar) / max(N_g,tiny);
lambda_w_g   = w_g_bar / max(N_g,tiny);
lambda_E_g   = (betta*M_g) / max(N_g,tiny);

lambda_eta_b = (kappa*eta_b_bar) / max(N_b,tiny);
lambda_w_b   = w_b_bar / max(N_b,tiny);
lambda_E_b   = (betta*M_b) / max(N_b,tiny);

% pipeline terms
chi_eta_g = (rho_g_bar + eta_g_bar) / max(rho_g_bar + 0.5*eta_g_bar, tiny);
chi_rho_g =  rho_g_bar              / max(rho_g_bar + 0.5*eta_g_bar, tiny);

chi_eta_b = (rho_b_bar + eta_b_bar) / max(rho_b_bar + 0.5*eta_b_bar, tiny);
chi_rho_b =  rho_b_bar              / max(rho_b_bar + 0.5*eta_b_bar, tiny);

% ---------- (h) Wage elasticities to TFP (optional, if used in model) ----------
omega_alpha_g = clip01( (eta*phi_g) / max(w_g_bar, tiny) );  % = ∂ln w_g / ∂ln α
omega_alpha_b = clip01( (eta*phi_b) / max(w_b_bar, tiny) );

% ---------- (i) Demand-side continuation weights (from transition matrix) ----------
Theta_gg = clip01( Pbar(3,3) );   % g -> g
Theta_gb = clip01( Pbar(3,2) );   % g -> b
Theta_gu = clip01( Pbar(3,1) );   % g -> u

Theta_bg = clip01( Pbar(2,3) );   % b -> g
Theta_bb = clip01( Pbar(2,2) );   % b -> b
Theta_bu = clip01( Pbar(2,1) );   % b -> u

Lambda_g = 1;                      % baseline under log utility
Lambda_b = 1;

% ---------- (j) Consumption weights and unemployment split ----------
omega_cg = clip01( g_bar );
omega_cb = clip01( b_bar );
omega_cu = clip01( u_bar );

omega_u_g = clip01( g_bar / max(u_bar, tiny) );
omega_u_b = clip01( b_bar / max(u_bar, tiny) );

% ---------- (k) Hiring-cost share in output ----------
CH_bar = 0.5*kappa*( eta_g_bar^2 * g_bar + eta_b_bar^2 * b_bar );
y_bar  = 1;                         % gaps normalization; change if you target a level
chiH   = CH_bar / max(y_bar, tiny);

% HH Laws of Motion

xi_bb   = b_bar*(sig_b*(1-varrho)*Phi + (1-sig_b)*(1-zeta*varrho*Phi)-...
          sig_b*(1-Phi));
xi_bg   = g_bar*((sig_b + (1-sig_b)*zeta)*varrho*Phi - sig_b*(1-Phi));
xi_bphi = Phi*((sig_b + (1-sig_b)*zeta)*varrho*g_bar + sig_b*(1-varrho)*b_bar...
          - sig_b*(1 - b_bar - g_bar) - (1-sig_b)*zeta*varrho*b_bar);
xi_gb   = b_bar*sig_g*((1-varrho)*Phi - (1-Phi));
xi_gg   = g_bar*(sig_g*varrho*Phi + (1-sig_g) - sig_g*(1-Phi));
xi_gphi = Phi*(sig_g*(1-varrho)*b_bar - sig_g*(1-g_bar - b_bar) + sig_g*varrho*g_bar);


% %% The conditions that need to be satisfied
% residuals=NaN(5,1);
% residuals(1)=N^(1-alfa)-(C+delta*N*Gamma*x^gammma); %(48)
% residuals(2)=(1-betta*(1-delta))*Gamma*x^gammma-xi*((1-alfa)*N^(-alfa)-chi*C*L^varphi);%(49)
% residuals(3)=(1-x)*xi*psi*chi*C*L^varphi-(1-xi)*Gamma*x^(1+gammma); %(50)
% residuals(4)=x*U-(1-x)*delta*N; %(51)
% residuals(5)=L-(N+psi*U); %(52)
% 
% if max(abs(residuals)>1e-6)
%     check=1; %set failure indicator
%     return; %return without updating steady states
% end


%% end own model equations

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') =  0 ;']);
end