function [ys,params,check] = paper2_V8_steadystate(ys,~,M_,~)

%% ---------- 0. boiler-plate ----------
check  = 0;
params = M_.params;

% Parameters from Dynare
beta          = params(strcmp('beta',          M_.param_names));
h             = params(strcmp('h',             M_.param_names));
delta         = params(strcmp('delta',         M_.param_names));
gamma         = params(strcmp('gamma',         M_.param_names));
OmegaK_par    = params(strcmp('OmegaK_par',    M_.param_names));

alpha_L       = params(strcmp('alpha_L',       M_.param_names));
alpha_K       = params(strcmp('alpha_K',       M_.param_names));
omega         = params(strcmp('omega',         M_.param_names));
Omega_par     = params(strcmp('Omega_par',     M_.param_names));

alpha_B       = params(strcmp('alpha_B',       M_.param_names));
alpha_G       = params(strcmp('alpha_G',       M_.param_names));
eps_B         = params(strcmp('eps_B',         M_.param_names));
eps_G         = params(strcmp('eps_G',         M_.param_names));
Omega_B_par   = params(strcmp('Omega_B_par',   M_.param_names));
Omega_G_par   = params(strcmp('Omega_G_par',   M_.param_names));

gamma0        = params(strcmp('gamma0',        M_.param_names));
gamma1        = params(strcmp('gamma1',        M_.param_names));
gamma2        = params(strcmp('gamma2',        M_.param_names));
delta_M       = params(strcmp('delta_M',       M_.param_names));
phi           = params(strcmp('phi',           M_.param_names));
upsilon       = params(strcmp('upsilon',       M_.param_names));
phi1          = params(strcmp('phi1',          M_.param_names));
phi2          = params(strcmp('phi2',          M_.param_names));

phi_E         = params(strcmp('phi_E',         M_.param_names));
sigma         = params(strcmp('sigma',         M_.param_names));

g             = params(strcmp('g',             M_.param_names));
kappa         = params(strcmp('kappa',         M_.param_names));
varsigma      = params(strcmp('varsigma',      M_.param_names));
xi            = params(strcmp('xi',            M_.param_names));
delta_B       = params(strcmp('delta_B',       M_.param_names));
delta_G       = params(strcmp('delta_G',       M_.param_names));
delta_Xi      = params(strcmp('delta_Xi',      M_.param_names));
Omega_BK_par  = params(strcmp('Omega_BK_par',  M_.param_names));
Omega_GK_par  = params(strcmp('Omega_GK_par',  M_.param_names));

A_ss          = params(strcmp('A_ss',          M_.param_names));
IB_ss         = params(strcmp('IB_ss',         M_.param_names));
IG_ss         = params(strcmp('IG_ss',         M_.param_names));
tauZ_ss       = params(strcmp('tauZ_ss',       M_.param_names));
pi_ss         = params(strcmp('pi_ss',         M_.param_names));
piB_ss        = params(strcmp('piB_ss',        M_.param_names));
piG_ss        = params(strcmp('piG_ss',        M_.param_names));

%% ---------- 1. numerical core via fsolve ----------

% x = [Z , L , pB]
x0  = [0.02 ; 0.33 ; 1];
opt = optimoptions('fsolve','Display','off','TolFun',1e-14,'TolX',1e-14, ...
                   'MaxFunctionEvaluations',1e5,'MaxIterations',2e3);

[x_ss,~,exitflag] = fsolve(@core_residual,x0,opt);

if exitflag<=0
    check = 1;
    return
end

Z  = x_ss(1);
L  = x_ss(2);
pB = x_ss(3);

%% ---------- 2. rebuild the full steady state ----------

A    = A_ss;
IB   = IB_ss;
IG   = IG_ss;
tauZ = tauZ_ss;
pi   = pi_ss;
piB  = piB_ss;
piG  = piG_ss;

aI   = (varsigma - 1)/varsigma;
IE   = (kappa*IG^aI + (1-kappa)*IB^aI)^(1/aI);

OmegaK   = 0;
OmegaB   = 0;
OmegaG   = 0;
Omega    = 0;
OmegaBK  = 0;
OmegaGK  = 0;

MC   = (omega - 1)/omega;
MCB  = (eps_B - 1)/eps_B;
MCG  = (eps_G - 1)/eps_G;

IBK   = (1 - xi) * IB;
KB    = IBK/delta_B;
IBXi  = xi*IB + (tauZ*Z)/2;
Xi    = IBXi/delta_Xi;
phi_t = phi * exp(-upsilon*Xi);
U     = ((phi_t * tauZ)/(phi1 * phi2))^(1/(phi2 - 1));

M     = Z/delta_M;
D     = 1 - (gamma0 + gamma1*M + gamma2*M^2);

chiB    = A * D * KB^alpha_B;
Z_model = (1 - U) * phi_t * chiB;
CA      = phi1 * U^phi2 * chiB;
rB      = piB * MCB * alpha_B * chiB/KB;

KG    = (IG + (tauZ*Z)/2)/delta_G;
chiG  = A * D * KG^alpha_G;
rG    = piG * MCG * alpha_G * chiG/KG;

E     = (phi_E * chiG^((sigma - 1)/sigma) + (1 - phi_E) * chiB^((sigma - 1)/sigma))^(sigma/(sigma - 1));
q     = (phi_E/(1-phi_E)) * (chiG/chiB)^(-1/sigma);
Fq    = (phi_E^sigma * q^(1 - sigma) + (1 - phi_E)^sigma)^(1/(1 - sigma));
piE   = 1;

r     = pi * (1/beta - (1 - delta));
K_Y   = (MC * alpha_K)/(r/pi);
Y     = (D * A * L^alpha_L * K_Y^alpha_K * E^(1 - alpha_L - alpha_K))^(1/(1 - alpha_K));
K     = K_Y * Y;
I     = delta * K;
W     = MC * alpha_L * (Y/L);
PE    = MC * (1 - alpha_L - alpha_K) * (Y/E);

Psi   = Y - W*L - r*K/pi - PE*E - Omega*Y;

pG    = q * pB;

PsiB  = pB*chiB - pB*(rB*KB/piB) - CA - tauZ*Z - pB*OmegaB*chiB;
PsiG  = pG*chiG - pG*(rG*KG/piG) - pG*OmegaG*chiG;

G     = g * Y;
T     = G + IE - pB*(rB*KB/piB) - pG*(rG*KG/piG) - PsiB - PsiG;
C     = r*K/pi + W*L + Psi - T - I;

lambda  = (1 - beta*h)/((1 - h) * C);
lambdaK = lambda;

% store back to Dynare's parameter vector
params(strcmp('r_ss',    M_.param_names)) = r;
params(strcmp('Y_ss',    M_.param_names)) = Y;
params(strcmp('piE_ss',  M_.param_names)) = piE;

%% ---------- 3. load results into ys ----------
ys(strcmp('C',        M_.endo_names)) = C;
ys(strcmp('I',        M_.endo_names)) = I;
ys(strcmp('r',        M_.endo_names)) = r;
ys(strcmp('K',        M_.endo_names)) = K;
ys(strcmp('pi',       M_.endo_names)) = pi;
ys(strcmp('W',        M_.endo_names)) = W;
ys(strcmp('L',        M_.endo_names)) = L;
ys(strcmp('Psi',      M_.endo_names)) = Psi;
ys(strcmp('T',        M_.endo_names)) = T;
ys(strcmp('OmegaK',   M_.endo_names)) = OmegaK;
ys(strcmp('lambda',   M_.endo_names)) = lambda;
ys(strcmp('lambdaK',  M_.endo_names)) = lambdaK;

ys(strcmp('chiB',     M_.endo_names)) = chiB;
ys(strcmp('A',        M_.endo_names)) = A;
ys(strcmp('M',        M_.endo_names)) = M;
ys(strcmp('Z',        M_.endo_names)) = Z;
ys(strcmp('U',        M_.endo_names)) = U;
ys(strcmp('phi_t',    M_.endo_names)) = phi_t;
ys(strcmp('CA',       M_.endo_names)) = CA;
ys(strcmp('PsiB',     M_.endo_names)) = PsiB;
ys(strcmp('rB',       M_.endo_names)) = rB;
ys(strcmp('piB',      M_.endo_names)) = piB;
ys(strcmp('tauZ',     M_.endo_names)) = tauZ;
ys(strcmp('MCB',      M_.endo_names)) = MCB;
ys(strcmp('OmegaB',   M_.endo_names)) = OmegaB;

ys(strcmp('chiG',     M_.endo_names)) = chiG;
ys(strcmp('KG',       M_.endo_names)) = KG;
ys(strcmp('PsiG',     M_.endo_names)) = PsiG;
ys(strcmp('rG',       M_.endo_names)) = rG;
ys(strcmp('piG',      M_.endo_names)) = piG;
ys(strcmp('MCG',      M_.endo_names)) = MCG;
ys(strcmp('OmegaG',   M_.endo_names)) = OmegaG;

ys(strcmp('E',        M_.endo_names)) = E;
ys(strcmp('q',        M_.endo_names)) = q;
ys(strcmp('piE',      M_.endo_names)) = piE;

ys(strcmp('Y',        M_.endo_names)) = Y;
ys(strcmp('PE',       M_.endo_names)) = PE;
ys(strcmp('MC',       M_.endo_names)) = MC;
ys(strcmp('Omega',    M_.endo_names)) = Omega;

ys(strcmp('G',        M_.endo_names)) = G;
ys(strcmp('IE',       M_.endo_names)) = IE;
ys(strcmp('IG',       M_.endo_names)) = IG;
ys(strcmp('IB',       M_.endo_names)) = IB;
ys(strcmp('IBK',      M_.endo_names)) = IBK;
ys(strcmp('IBXi',     M_.endo_names)) = IBXi;
ys(strcmp('KB',       M_.endo_names)) = KB;
ys(strcmp('Xi',       M_.endo_names)) = Xi;
ys(strcmp('OmegaBK',  M_.endo_names)) = OmegaBK;
ys(strcmp('OmegaGK',  M_.endo_names)) = OmegaGK;

ys(strcmp('pG',       M_.endo_names)) = pG;
ys(strcmp('pB',       M_.endo_names)) = pB;

%% ========== NESTED RESIDUAL FUNCTION =================================
    function F = core_residual(x)

        Z  = x(1);
        L  = x(2);
        pB = x(3);

        if ~isreal(Z) || ~isreal(L) || ~isreal(pB) || Z<=0 || L<=0 || L>=1 || pB<=0
            F = 1e8 * ones(3,1);
            return
        end

        if phi_E<=0 || phi_E>=1 || sigma<=0 || varsigma<=0 || delta_B<=0 || delta_G<=0 || delta_Xi<=0 || delta_M<=0
            F = 1e8 * ones(3,1);
            return
        end

        A    = A_ss;
        IB   = IB_ss;
        IG   = IG_ss;
        tauZ = tauZ_ss;
        pi   = pi_ss;
        piB  = piB_ss;
        piG  = piG_ss;

        aI   = (varsigma - 1)/varsigma;
        IE   = (kappa*IG^aI + (1-kappa)*IB^aI)^(1/aI);

        MC   = (omega - 1)/omega;
        MCB  = (eps_B - 1)/eps_B;
        MCG  = (eps_G - 1)/eps_G;

        IBK  = (1 - xi) * IB;
        if IBK<=0
            F = 1e8 * ones(3,1);
            return
        end

        KB    = IBK/delta_B;
        IBXi  = xi*IB + (tauZ*Z)/2;
        Xi    = IBXi/delta_Xi;
        phi_t = phi * exp(-upsilon*Xi);
        if phi_t<=0
            F = 1e8 * ones(3,1);
            return
        end

        inside_U = (phi_t * tauZ)/(phi1 * phi2);
        if inside_U<=0
            F = 1e8 * ones(3,1);
            return
        end

        U = inside_U^(1/(phi2 - 1));

        M = Z/delta_M;
        D = 1 - (gamma0 + gamma1*M + gamma2*M^2);
        if ~isreal(D) || D<=0
            F = 1e8 * ones(3,1);
            return
        end

        chiB = A * D * KB^alpha_B;
        if ~isreal(chiB) || chiB<=0
            F = 1e8 * ones(3,1);
            return
        end

        Z_model = (1 - U) * phi_t * chiB;
        CA      = phi1 * U^phi2 * chiB;
        rB      = piB * MCB * alpha_B * chiB/KB;

        KG = (IG + (tauZ*Z)/2)/delta_G;
        if ~isreal(KG) || KG<=0
            F = 1e8 * ones(3,1);
            return
        end

        chiG = A * D * KG^alpha_G;
        if ~isreal(chiG) || chiG<=0
            F = 1e8 * ones(3,1);
            return
        end

        rG = piG * MCG * alpha_G * chiG/KG;

        E = (phi_E * chiG^((sigma - 1)/sigma) + (1 - phi_E) * chiB^((sigma - 1)/sigma))^(sigma/(sigma - 1));
        if ~isreal(E) || E<=0
            F = 1e8 * ones(3,1);
            return
        end

        q  = (phi_E/(1-phi_E)) * (chiG/chiB)^(-1/sigma);
        Fq = (phi_E^sigma * q^(1 - sigma) + (1 - phi_E)^sigma)^(1/(1 - sigma));

        r   = pi * (1/beta - (1 - delta));
        K_Y = (MC * alpha_K)/(r/pi);
        if ~isreal(K_Y) || K_Y<=0
            F = 1e8 * ones(3,1);
            return
        end

        Y = (D * A * L^alpha_L * K_Y^alpha_K * E^(1 - alpha_L - alpha_K))^(1/(1 - alpha_K));
        if ~isreal(Y) || Y<=0
            F = 1e8 * ones(3,1);
            return
        end

        K  = K_Y * Y;
        I  = delta * K;
        W  = MC * alpha_L * (Y/L);
        PE = MC * (1 - alpha_L - alpha_K) * (Y/E);

        Psi  = Y - W*L - r*K/pi - PE*E;

        pG   = q * pB;
        PsiB = pB*chiB - pB*(rB*KB/piB) - CA - tauZ*Z;
        PsiG = pG*chiG - pG*(rG*KG/piG);

        G = g * Y;
        T = G + IE - pB*(rB*KB/piB) - pG*(rG*KG/piG) - PsiB - PsiG;
        C = r*K/pi + W*L + Psi - T - I;
        if ~isreal(C) || C<=0
            F = 1e8 * ones(3,1);
            return
        end

        lambda = (1 - beta*h)/((1 - h) * C);

        F    = zeros(3,1);
        F(1) = Z - Z_model;
        F(2) = L - (1 - gamma/(lambda*W));
        F(3) = pB - PE/Fq;

    end
% ========== end of nested function ====================================
end