function [ys,params,check] = design_cbdc_steadystate(ys,exo,M_,options_)
% =========================================================================
% External steady state for design_cbdc.mod
% - Computes the "macro core" steady state analytically
% - Solves the liquidity + CBDC/SC + foreign liquidity + issuer reserves
%   block numerically using FSOLVE
%
% Key robustness features:
%   - No eval() for parameter reading (avoids name conflicts like beta())
%   - Numerical block solved in logs (enforces positivity)
%   - ys fill is guarded: if a variable is not computed, keep incoming guess
%
% NOTE:
%   This file is written to match a .mod that contains:
%     - Household liquidity FOCs for M / CBDC / SC
%     - Foreign liquidity FOCs for Mf / DCf (with DCf = SCf)
%     - Issuer equations: lambdaSC definition + FOCs for MSCf and BSCf
%     - Clearing: SCs = SC + SCf, and Msf = Mf + MSCf
% =========================================================================

%% ---- Read parameters safely (avoid beta() / gamma() conflicts) ----------
p = cell2struct(num2cell(M_.params), cellstr(M_.param_names), 1);

% Explicit parameter unpacking (add any additional parameters you use below)
beta      = p.beta;
betaf     = p.betaf;

sigma     = p.sigma;
sigmaf    = p.sigmaf;

alpha     = p.alpha;
alphaf    = p.alphaf;

delta     = p.delta;
deltaf    = p.deltaf;

g         = p.g;
gf        = p.gf;

chi       = p.chi;
chif      = p.chif;

eta       = p.eta;
etaf      = p.etaf;

tau       = p.tau;
tauf      = p.tauf;

nu        = p.nu;
nuf       = p.nuf;

theta     = p.theta;
thetaf    = p.thetaf;

omega     = p.omega;
omegaf    = p.omegaf;

rho       = p.rho;
rhof      = p.rhof;

% Liquidity / CBDC / stablecoin parameters
mu        = p.mu;
muf       = p.muf;

chiL      = p.chiL;
chiLf     = p.chiLf;

rhoDC     = p.rhoDC;
rhoDCf    = p.rhoDCf;

chiDC     = p.chiDC;
chiDCf    = p.chiDCf;

Gamma     = p.Gamma;
Gammaf    = p.Gammaf;

omegaDC   = p.omegaDC;
varepsilon = p.varepsilon;

spread    = p.spread;

% Issuer parameters (foreign)
xi        = p.xi;
varrho    = p.varrho;
chiSCf    = p.chiSCf;
GammaSCf  = p.GammaSCf;

%% ---- Init indicator -----------------------------------------------------
check = 0;

%% ---- Normalizations / steady-state targets ------------------------------
P = 1;  Pf = 1;
PH = 1; PHf = 1;
PF = 1; PFf = 1;

PH_hat  = 1;  PHf_hat = 1;
PF_hat  = 1;  PFf_hat = 1;

NER = 1;
RER = 1;

pi  = 1;   pif  = 1;
piH = 1;   piHf = 1;
piF = 1;   piFf = 1;

A  = 1;    Af = 1;

% Monetary policy shocks in steady state
uR  = 0;
uRf = 0;

%% ---- Policy rates and returns -------------------------------------------
R  = 1/beta;
Rf = 1/betaf;

RD  = R;   RL  = RD;
RDf = Rf;  RLf = RDf;

R_CBDC = R - spread;

% Marginal costs (consistent with your Calvo block normalization)
MC  = (1+tau) *((nu -1)/nu);
MCf = (1+tauf)*((nuf-1)/nuf);

% Capital rental rates (gross real)
RK  = 1/beta  - (1-delta);
RKf = 1/betaf - (1-deltaf);

%% ---- Core steady state: K/N, W, N, C, K, I, ... -------------------------
k_n  = (MC *alpha /RK )^(1/(1-alpha));
k_nf = (MCf*alphaf/RKf)^(1/(1-alphaf));

W  = ((1-alpha)/alpha)  * k_n  * RK;
Wf = ((1-alphaf)/alphaf)* k_nf * RKf;

N  = (1/(((1-g )*k_n^(alpha -1)-delta )*(W /chi )^sigma *k_n ))^(1/(eta +sigma +1));
Nf = (1/(((1-gf)*k_nf^(alphaf-1)-deltaf)*(Wf/chif)^sigmaf*k_nf))^(1/(etaf+sigmaf+1));

C  = ((W /chi )*N ^(-eta ))^(-sigma );
Cf = ((Wf/chif)*Nf^(-etaf))^(-sigmaf);

K  = k_n  * N;
Kf = k_nf * Nf;

I  = delta  * K;
If = deltaf * Kf;

% Banking/entrepreneur identities in your model
L  = I;   D  = L;
Lf = If;  Df = Lf;

lambda  = C ^(-sigma );
lambdaf = Cf^(-sigmaf);

Y  = (C + I )/(1-g );
Yf = (Cf+ If)/(1-gf);

G  = g  * Y;
Gf = gf * Yf;

% Q Euler in your .mod implies Q_ss = lambda_ss
Q  = lambda;
Qf = lambdaf;

% Home/foreign demand components
YH  = omega  *(PH /P )^(1/(rho -1))*(C + I + G );
XFf = (1-omega)*(PF /P )^(1/(rho -1))*(C + I + G );

YHf = omegaf *(PHf/Pf)^(1/(rhof-1))*(Cf+ If+ Gf);
XF  = (1-omegaf)*(PFf/Pf)^(1/(rhof-1))*(Cf+ If+ Gf);

% Calvo auxiliary recursions in steady state (given your MC choice, KH=FH etc.)
FH  = (lambda *YH )/(1-beta *theta );
KH  = FH;
FF  = (lambda *XF )/(1-beta *theta );
KF  = FF;

FHf = (lambdaf*YHf)/(1-betaf*thetaf);
KHf = FHf;
FFf = (lambdaf*XFf)/(1-betaf*thetaf);
KFf = FFf;

%% ---- Liquidity + CBDC/SC + issuer reserves: numerical block --------------
% Unknown vector in logs (positivity):
%   [M, CBDC, SC, Mf, SCf, MSCf, BSCf]
% Equations:
%   Domestic liquidity FOCs (8)(9)(10): 3 eqs
%   Foreign liquidity FOCs (18)(19):     2 eqs
%   Issuer reserve FOCs (22)(23):        2 eqs
%   (lambdaSC is pinned by its definition in your .mod: lambdaSC = 1 - betaf)

lambdaSC = 1 - betaf;  % steady state: lambdaf(+1)/lambdaf = 1

% helper: pull a positive initial guess from incoming ys if present
get0 = @(name, fallback) local_get_init_guess(ys, M_, name, fallback);

M0     = get0('M',     0.10*C);
CBDC0  = get0('CBDC',  0.05*C);
SC0    = get0('SC',    0.05*C);

Mf0    = get0('Mf',    0.10*Cf);
SCf0   = get0('SCf',   0.05*Cf);

MSCf0  = get0('MSCf',  0.05*(C+Cf));
BSCf0  = get0('BSCf',  0.05*(C+Cf));

u0 = log([M0; CBDC0; SC0; Mf0; SCf0; MSCf0; BSCf0]);

% fsolve options
if exist('optimoptions','file') == 2
    opt = optimoptions('fsolve', ...
        'Display','off', ...
        'FunctionTolerance',1e-12, ...
        'StepTolerance',1e-12, ...
        'MaxIterations',5000, ...
        'MaxFunctionEvaluations',100000);
else
    opt = optimset('Display','off','TolFun',1e-12,'TolX',1e-12,'MaxIter',5000,'MaxFunEvals',100000);
end

F = @(u) local_residuals_liq_issuer7(u, ...
    lambda,lambdaf,beta,betaf,pi,pif,R_CBDC,NER,Rf, ...
    mu,rhoDC,chiL,chiDC,Gamma,omegaDC,varepsilon, ...
    muf,rhoDCf,chiLf,chiDCf,Gammaf, ...
    xi,varrho,chiSCf,GammaSCf, ...
    lambdaSC);

[us,~,exitflag] = fsolve(F, u0, opt);

if exitflag <= 0 || any(~isfinite(us))
    check = 1;
    return;
end

x = exp(us);
M    = x(1);
CBDC = x(2);
SC   = x(3);
Mf   = x(4);
SCf  = x(5);
MSCf = x(6);
BSCf = x(7);

% Back out definitional objects
DC  = ( omegaDC*CBDC^((varepsilon-1)/varepsilon) + (1-omegaDC)*(NER*SC)^((varepsilon-1)/varepsilon) )^(varepsilon/(varepsilon-1));
Liq = ( mu*M^((rhoDC-1)/rhoDC) + (1-mu)*DC^((rhoDC-1)/rhoDC) )^(rhoDC/(rhoDC-1));

DCf  = SCf;
Liqf = ( muf*Mf^((rhoDCf-1)/rhoDCf) + (1-muf)*DCf^((rhoDCf-1)/rhoDCf) )^(rhoDCf/(rhoDCf-1));

SCs = SC + SCf;
Msf = Mf + MSCf;

%% ---- Return params (no parameter back-out in this file) ------------------
params = M_.params;

%% ---- Fill ys (guarded: keep incoming guess if not computed) -------------
NumberOfEndogenousVariables = M_.orig_endo_nbr;
missing_any = false;

for ii = 1:NumberOfEndogenousVariables
    varname = M_.endo_names{ii};
    if exist(varname,'var')
        ys(ii) = eval(varname); %#ok<EVLDIR>
    else
        % Keep incoming ys(ii) as fallback; flag for diagnostics
        missing_any = true;
    end
end

if missing_any
    % If anything is missing, force Dynare to stop (better than silently wrong SS)
    check = 1;
end

end

% =========================================================================
% Local helpers
% =========================================================================
function val = local_get_init_guess(ys, M_, name, fallback)
idx = strmatch(name, M_.endo_names, 'exact'); %#ok<MATCH2>
if isempty(idx)
    val = fallback;
    return;
end
val = ys(idx(1));
if ~isfinite(val) || val <= 0
    val = fallback;
end
val = max(val, 1e-10);
end

function F = local_residuals_liq_issuer7(u, ...
    lambda,lambdaf,beta,betaf,pi,pif,R_CBDC,NER,Rf, ...
    mu,rhoDC,chiL,chiDC,Gamma,omegaDC,varepsilon, ...
    muf,rhoDCf,chiLf,chiDCf,Gammaf, ...
    xi,varrho,chiSCf,GammaSCf, ...
    lambdaSC)

x = exp(u);
M    = x(1);
CBDC = x(2);
SC   = x(3);
Mf   = x(4);
SCf  = x(5);
MSCf = x(6);
BSCf = x(7);

% Definitions
DC  = ( omegaDC*CBDC^((varepsilon-1)/varepsilon) + (1-omegaDC)*(NER*SC)^((varepsilon-1)/varepsilon) )^(varepsilon/(varepsilon-1));
Liq = ( mu*M^((rhoDC-1)/rhoDC) + (1-mu)*DC^((rhoDC-1)/rhoDC) )^(rhoDC/(rhoDC-1));

DCf  = SCf;
Liqf = ( muf*Mf^((rhoDCf-1)/rhoDCf) + (1-muf)*DCf^((rhoDCf-1)/rhoDCf) )^(rhoDCf/(rhoDCf-1));

SCs = SC + SCf;

% RHS terms in steady state
rhs_M   = lambda  - beta *(lambda /pi);
rhs_C   = lambda  - beta *(lambda /pi)*R_CBDC;
rhs_SC  = lambda  - beta *(lambda /pi)*(NER/NER); % = lambda*(1-beta) in SS

rhs_Mf  = lambdaf - betaf*(lambdaf/pif);
rhs_DCf = lambdaf - betaf*(lambdaf/pif);

% Issuer RHS
rhs_iss_MSCf = 1 - betaf;            % since lambdaf(+1)/lambdaf = 1 in SS
rhs_iss_BSCf = 1 - betaf*Rf/pif;     % = 0 if Rf=1/betaf and pif=1

% Domestic penalty helper
term_pen_dom = (Gamma - M/(M+DC));
penM  = 2*chiDC*DC/((M+DC)^2)*term_pen_dom;
penDC = 2*chiDC*M /((M+DC)^2)*term_pen_dom;

% Domestic FOCs (8)(9)(10)
eq8  = (chiL/Liq)*mu*(Liq/M)^(1/rhoDC) + penM - rhs_M;

chainCBDC = omegaDC*(DC/CBDC)^(1/varepsilon);
eq9  = (chiL/Liq)*(1-mu)*(Liq/DC)^(1/rhoDC)*chainCBDC - penDC*chainCBDC - rhs_C;

% Keep the same precedence as your .mod if you wrote DC/NER*SC (NOT DC/(NER*SC))
chainSC = (1-omegaDC)*(DC/NER*SC)^(1/varepsilon);
eq10 = (chiL/Liq)*(1-mu)*(Liq/DC)^(1/rhoDC)*chainSC - penDC*chainSC - rhs_SC;

% Foreign FOCs (18)(19)
term_pen_for = (Gammaf - Mf/(Mf+DCf));
eq18 = (chiLf/Liqf)*muf*(Liqf/Mf)^(1/rhoDCf) + 2*chiDCf*DCf/((Mf+DCf)^2)*term_pen_for - rhs_Mf;
eq19 = (chiLf/Liqf)*(1-muf)*(Liqf/DCf)^(1/rhoDCf) - 2*chiDCf*Mf/((Mf+DCf)^2)*term_pen_for - rhs_DCf;

% Issuer FOCs (22)(23) with fixed lambdaSC
eq22 = lambdaSC*xi^(1/varrho)*(SCs/MSCf)^(1/varrho) ...
       + 2*chiSCf*BSCf/((MSCf+BSCf)^2)*(GammaSCf - MSCf/(MSCf+BSCf)) - rhs_iss_MSCf;

eq23 = lambdaSC*(1-xi)^(1/varrho)*(SCs/BSCf)^(1/varrho) ...
       - 2*chiSCf*MSCf/((MSCf+BSCf)^2)*(GammaSCf - MSCf/(MSCf+BSCf)) - rhs_iss_BSCf;

F = [eq8; eq9; eq10; eq18; eq19; eq22; eq23];

end
