"""Solvers for the WP19 model: seller values, stationary equilibrium with free entry, posted prices,
financing contingency extension, planner benchmark and transition dynamics.

Equations (model_specification.md):
  (S)  (r+delta) V(c) = -c + delta*w0 + Lam_s * beta * H(V(c)+U)
  (B)  r U = -kappa_b + Lam_b (1-beta) sum_c H(V(c)+U) N(c)
  free entry U = k; stocks n(c) = s f(c) / (eta(c)+delta), eta(c) = Lam_s (1 - G(V(c)+U)).
"""
from __future__ import annotations

from dataclasses import dataclass, field
from typing import Optional

import numpy as np
from scipy import optimize

from .model import Params


# ----------------------------------------------------------------------------------------------
# Seller block (partial equilibrium: Lam_s and U given)
# ----------------------------------------------------------------------------------------------
def seller_value(c: float, Lam_s: float, U: float, p: Params, xtol: float = 1e-12) -> float:
    """Unique root of Phi(V) = -c + delta*w0 + Lam_s*beta*H(V+U) - (r+delta) V (Lemma 2)."""
    G, r, d = p.G, p.r, p.delta

    def Phi(V):
        return -c + d * p.w0 + Lam_s * p.beta * float(G.H(V + U)) - (r + d) * V

    lo = (-c + d * p.w0) / (r + d) - 1e-9          # Phi(lo) >= 0 since H >= 0
    hi_support = G.hi if np.isfinite(G.hi) else G.mean() + 50.0
    hi = hi_support - U + 1e-9                     # H = 0 beyond the support => Phi(hi) < 0
    hi = max(hi, lo + 1.0)
    guard = 0
    while Phi(hi) >= 0 and guard < 60:             # unbounded (heavy-tailed) support safeguard
        hi *= 2.0; guard += 1
    return optimize.brentq(Phi, lo, hi, xtol=xtol, rtol=4 * np.finfo(float).eps, maxiter=500)


@dataclass
class SellerBlock:
    V: np.ndarray
    vstar: np.ndarray
    acc: np.ndarray        # conditional sale probability 1-G(v*)
    eta: np.ndarray        # sale hazard
    pbar: np.ndarray       # expected transaction price by type
    M: np.ndarray          # mean residual valuation at v*
    H: np.ndarray
    resid: np.ndarray      # residual of (S)


def solve_sellers(Lam_s: float, U: float, p: Params) -> SellerBlock:
    V = np.array([seller_value(float(ci), Lam_s, U, p, p.root_xtol) for ci in p.c])
    vstar = V + U
    acc = np.asarray(p.G.surv(vstar), dtype=float)
    H = np.asarray(p.G.H(vstar), dtype=float)
    M = np.where(acc > 0, H / np.maximum(acc, 1e-300), 0.0)
    eta = Lam_s * acc
    pbar = V + p.beta * M
    resid = -p.c + p.delta * p.w0 + Lam_s * p.beta * H - (p.r + p.delta) * V
    return SellerBlock(V, vstar, acc, eta, pbar, M, H, resid)


def stationary_stocks(eta: np.ndarray, p: Params):
    n = p.s * p.f / (eta + p.delta)
    S = n.sum()
    N = n / S
    Q = float((eta * n).sum())
    omega = eta * n / Q if Q > 0 else N
    return n, S, N, Q, omega


# ----------------------------------------------------------------------------------------------
# Stationary equilibrium with free entry (Nash bargaining, baseline screen phi)
# ----------------------------------------------------------------------------------------------
@dataclass
class Equilibrium:
    theta: float
    Lam_s: float
    Lam_b: float
    U: float
    sellers: SellerBlock
    n: np.ndarray
    S: float
    B: float
    Q: float
    omega: np.ndarray
    N: np.ndarray
    free_entry_resid: float
    converged: bool
    theta_roots: list = field(default_factory=list)
    protocol: str = "nash"
    extra: dict = field(default_factory=dict)

    # --- moments ----------------------------------------------------------------------------
    def moments(self, p: Params) -> dict:
        sb = self.sellers
        G = p.G
        pT = float((self.omega * sb.pbar).sum())
        # within-type variance of prices: beta^2 Var(w | w>=v*)
        cvar = np.array([G.cond_var(v) for v in sb.vstar])
        within = float((self.omega * (p.beta ** 2) * cvar).sum())
        between = float((self.omega * (sb.pbar - pT) ** 2).sum())
        mean_dur_listing = float((p.f * 1.0 / (sb.eta + p.delta)).sum())      # expected duration of a new listing
        sold_share = float((p.f * sb.eta / (sb.eta + p.delta)).sum())
        # mean duration of completed sales: E[tau | sale, c] = 1/(eta+delta) (exponential competing risks)
        mean_tom_sold = float((self.omega / (sb.eta + p.delta)).sum())
        pstar = frictionless_price(p)
        D = 1.0 - sb.pbar / sb.pbar[0]
        Dbar = float((self.omega * D).sum())
        welfare = welfare_flow(self, p)
        return {
            "theta": self.theta, "Lam_s": self.Lam_s, "Lam_b": self.Lam_b,
            "S": self.S, "B": self.B, "Q": self.Q, "sold_share": sold_share,
            "mean_listing_duration": mean_dur_listing, "mean_tom_sold": mean_tom_sold,
            "mean_tom_sold_months": 12.0 * mean_tom_sold,
            "buyer_mean_duration": 1.0 / (self.Lam_b * float((self.N * sb.acc).sum())),
            "price_mean": pT, "price_var_within": within, "price_var_between": between,
            "price_sd_total": float(np.sqrt(within + between)),
            "price_between_share": between / (within + between) if within + between > 0 else np.nan,
            "V_min": float(sb.V.min()), "V_max": float(sb.V.max()),
            "pbar_patient": float(sb.pbar[0]), "pbar_urgent": float(sb.pbar[-1]),
            "D_urgent": float(D[-1]), "D_mean": Dbar,
            "eta_patient": float(sb.eta[0]), "eta_urgent": float(sb.eta[-1]),
            "acc_patient": float(sb.acc[0]), "acc_urgent": float(sb.acc[-1]),
            "p_star": pstar, "Gamma": (pstar - pT) / pstar,
            "participation_ok": bool((sb.V >= p.w0 - 1e-9).all()),
            "welfare_total": welfare["total"], "welfare_sellers": welfare["sellers"],
            "welfare_buyers_net": welfare["buyers_net"], "welfare_holding_costs": welfare["holding_costs"],
            "free_entry_resid": self.free_entry_resid, "seller_resid_max": float(np.abs(sb.resid).max()),
            "converged": self.converged, "n_theta_roots": len(self.theta_roots),
        }


def frictionless_price(p: Params) -> float:
    """Walrasian free-entry benchmark: H(p*) = r k."""
    G = p.G
    target = p.r * p.k
    hi = G.hi if np.isfinite(G.hi) else G.mean() + 50.0
    return optimize.brentq(lambda v: float(G.H(v)) - target, G.lo - 10.0, hi, xtol=1e-13)


def free_entry_residual(theta: float, p: Params, U: Optional[float] = None) -> float:
    U = p.k if U is None else U
    Lam_s = p.phi * float(p.lam_s(theta))
    Lam_b = p.phi * float(p.lam_b(theta))
    sb = solve_sellers(Lam_s, U, p)
    n, S, N, Q, omega = stationary_stocks(sb.eta, p)
    return -p.kappa_b + Lam_b * (1.0 - p.beta) * float((N * sb.H).sum()) - p.r * U


def _build_equilibrium(theta: float, p: Params, roots, converged: bool) -> Equilibrium:
    Lam_s = p.phi * float(p.lam_s(theta))
    Lam_b = p.phi * float(p.lam_b(theta))
    sb = solve_sellers(Lam_s, p.k, p)
    n, S, N, Q, omega = stationary_stocks(sb.eta, p)
    resid = -p.kappa_b + Lam_b * (1.0 - p.beta) * float((N * sb.H).sum()) - p.r * p.k
    return Equilibrium(theta, Lam_s, Lam_b, p.k, sb, n, S, theta * S, Q, omega, N, resid, converged, list(roots))


def solve_equilibrium(p: Params, check_grid: bool = True) -> Equilibrium:
    """Find theta such that the free-entry residual is zero. Scans a log grid for sign changes (reports multiplicity)."""
    lo, hi = p.theta_bounds
    roots = []
    if check_grid:
        grid = np.geomspace(lo, hi, p.grid_check_points)
        vals = np.array([free_entry_residual(t, p) for t in grid])
        for i in range(len(grid) - 1):
            if vals[i] == 0.0:
                roots.append(float(grid[i]))
            elif vals[i] * vals[i + 1] < 0:
                roots.append(optimize.brentq(free_entry_residual, grid[i], grid[i + 1], args=(p,), xtol=1e-12))
        if not roots:
            # no sign change: report the end with smallest |residual|
            j = int(np.argmin(np.abs(vals)))
            return _build_equilibrium(float(grid[j]), p, [], converged=False)
        theta = roots[0] if len(roots) == 1 else roots[int(np.argmin(np.abs(np.log(np.array(roots)))))]
        return _build_equilibrium(theta, p, roots, converged=True)
    f_lo, f_hi = free_entry_residual(lo, p), free_entry_residual(hi, p)
    if f_lo * f_hi > 0:
        return _build_equilibrium(lo if abs(f_lo) < abs(f_hi) else hi, p, [], converged=False)
    theta = optimize.brentq(free_entry_residual, lo, hi, args=(p,), xtol=1e-12)
    return _build_equilibrium(theta, p, [theta], converged=True)


def partial_equilibrium(theta: float, p: Params) -> Equilibrium:
    """Hold tightness fixed (U = k), re-solve the seller block and stocks."""
    return _build_equilibrium(theta, p, [theta], converged=True)


# ----------------------------------------------------------------------------------------------
# Welfare
# ----------------------------------------------------------------------------------------------
def welfare_flow(eq: Equilibrium, p: Params) -> dict:
    """Steady-state flow surplus: sum_c eta n (E[w|w>=v*] - w0) - sum_c c n - kappa_b B - k Q."""
    sb = eq.sellers
    gains = float((sb.eta * eq.n * (sb.vstar + sb.M - p.w0)).sum())
    holding = float((p.c * eq.n).sum())
    search = p.kappa_b * eq.B
    entry = p.k * eq.Q
    total = gains - holding - search - entry
    sellers = p.s * float((p.f * (sb.V - p.w0)).sum())       # ex-ante value of a listing cohort net of own-use value
    buyers_net = 0.0                                          # free entry
    return {"total": total, "gains": gains, "holding_costs": holding, "search_costs": search,
            "entry_costs": entry, "sellers": sellers, "buyers_net": buyers_net}


def planner_optimum(p: Params, x0: Optional[np.ndarray] = None) -> dict:
    """Constrained-efficient steady state: choose theta and acceptance thresholds v*(c) to maximise flow surplus
    subject to the stock equations. Deterministic Nelder-Mead from the decentralised allocation."""
    G = p.G

    def objective(x):
        theta = np.exp(x[0]); vstar = x[1:]
        Lam_s = p.phi * float(p.lam_s(theta))
        acc = np.asarray(G.surv(vstar)); eta = Lam_s * acc
        n = p.s * p.f / (eta + p.delta); S = n.sum(); B = theta * S; Q = (eta * n).sum()
        gains = (eta * n * (np.asarray(G.cond_mean(vstar)) - p.w0)).sum()
        return -(gains - (p.c * n).sum() - p.kappa_b * B - p.k * Q)

    if x0 is None:
        eq = solve_equilibrium(p)
        x0 = np.concatenate([[np.log(eq.theta)], eq.sellers.vstar])
    res = optimize.minimize(objective, x0, method="Nelder-Mead",
                            options={"xatol": 1e-9, "fatol": 1e-12, "maxiter": 40000, "maxfev": 80000})
    res2 = optimize.minimize(objective, res.x, method="Powell", options={"xtol": 1e-10, "ftol": 1e-13, "maxiter": 20000})
    best = res2 if res2.fun <= res.fun else res
    return {"theta": float(np.exp(best.x[0])), "vstar": best.x[1:], "welfare": -float(best.fun),
            "success": bool(best.success), "nfev": int(best.nfev) + int(res.nfev)}


# ----------------------------------------------------------------------------------------------
# Posted-price protocol (Proposition 8)
# ----------------------------------------------------------------------------------------------
def posted_threshold(a: float, p: Params) -> float:
    """y solving y - a = (1-G(y))/g(y) on (a, hi) (unique under log-concave g). Falls back to direct maximisation."""
    G = p.G
    hi = G.hi if np.isfinite(G.hi) else G.mean() + 50.0

    def foc(y):
        s = float(G.surv(y)); g = float(G.pdf(y))
        return (y - a) * g - s

    lo_y = max(a, G.lo) + 1e-10
    try:
        if foc(lo_y) < 0 and foc(hi - 1e-10) > 0:
            return optimize.brentq(foc, lo_y, hi - 1e-10, xtol=1e-12)
    except ValueError:
        pass
    res = optimize.minimize_scalar(lambda y: -float(G.surv(y)) * (y - a), bounds=(lo_y, hi), method="bounded",
                                   options={"xatol": 1e-12})
    return float(res.x)


def seller_value_posted(c: float, Lam_s: float, U: float, p: Params) -> tuple[float, float]:
    G, r, d = p.G, p.r, p.delta

    def J(V):
        y = posted_threshold(U + V, p)
        return float(G.surv(y)) * (y - U - V), y

    def Phi(V):
        return -c + d * p.w0 + Lam_s * J(V)[0] - (r + d) * V

    lo = (-c + d * p.w0) / (r + d) - 1e-9
    hi_support = G.hi if np.isfinite(G.hi) else G.mean() + 50.0
    hi = hi_support - U - 1e-9
    V = optimize.brentq(Phi, lo, hi, xtol=p.root_xtol, maxiter=500)
    return V, J(V)[1]


def solve_sellers_posted(Lam_s: float, U: float, p: Params) -> SellerBlock:
    out = [seller_value_posted(float(ci), Lam_s, U, p) for ci in p.c]
    V = np.array([o[0] for o in out]); y = np.array([o[1] for o in out])
    acc = np.asarray(p.G.surv(y), dtype=float); H = np.asarray(p.G.H(y), dtype=float)
    M = np.where(acc > 0, H / np.maximum(acc, 1e-300), 0.0)
    eta = Lam_s * acc
    pbar = y - U                                   # deterministic posted price
    resid = -p.c + p.delta * p.w0 + Lam_s * acc * (y - U - V) - (p.r + p.delta) * V
    return SellerBlock(V, y, acc, eta, pbar, M, H, resid)


def free_entry_residual_posted(theta: float, p: Params) -> float:
    Lam_s = p.phi * float(p.lam_s(theta)); Lam_b = p.phi * float(p.lam_b(theta))
    sb = solve_sellers_posted(Lam_s, p.k, p)
    n, S, N, Q, omega = stationary_stocks(sb.eta, p)
    # buyer surplus per meeting with a type-c seller: E[(w - y(c))^+] = H(y)
    return -p.kappa_b + Lam_b * float((N * sb.H).sum()) - p.r * p.k


def solve_equilibrium_posted(p: Params) -> Equilibrium:
    lo, hi = p.theta_bounds
    grid = np.geomspace(lo, hi, p.grid_check_points)
    vals = np.array([free_entry_residual_posted(t, p) for t in grid])
    roots = [optimize.brentq(free_entry_residual_posted, grid[i], grid[i + 1], args=(p,), xtol=1e-12)
             for i in range(len(grid) - 1) if vals[i] * vals[i + 1] < 0]
    converged = bool(roots)
    theta = roots[0] if roots else float(grid[int(np.argmin(np.abs(vals)))])
    Lam_s = p.phi * float(p.lam_s(theta)); Lam_b = p.phi * float(p.lam_b(theta))
    sb = solve_sellers_posted(Lam_s, p.k, p)
    n, S, N, Q, omega = stationary_stocks(sb.eta, p)
    resid = -p.kappa_b + Lam_b * float((N * sb.H).sum()) - p.r * p.k
    eq = Equilibrium(theta, Lam_s, Lam_b, p.k, sb, n, S, theta * S, Q, omega, N, resid, converged, roots, "posted")
    return eq


# ----------------------------------------------------------------------------------------------
# Financing-contingency extension (cash vs mortgage buyers, Proposition 7)
# ----------------------------------------------------------------------------------------------
@dataclass
class FinancingEquilibrium:
    theta: float
    lam_s: float
    lam_b: float
    U_c: float
    U_m: float
    pi_tilde: float
    V: np.ndarray
    vstar_c: np.ndarray
    vstar_m: np.ndarray
    eta_c: np.ndarray
    eta_m: np.ndarray
    n: np.ndarray
    m_pending: np.ndarray
    S: float
    B_c: float
    B_m: float
    Q: float
    Q_c: np.ndarray
    Q_m: np.ndarray
    pbar_c: np.ndarray
    pbar_m: np.ndarray
    Delta_b: np.ndarray
    sel_term: np.ndarray
    outside_term: float
    resid: dict
    converged: bool
    iterations: int


def _fin_state(theta: float, U_c: float, U_m: float, pit: float, p: Params) -> dict:
    """Seller block, stocks and buyer-side objects for given (theta, U_c, U_m, pi_tilde)."""
    G = p.G
    lam_s = float(p.lam_s(theta)); lam_b = float(p.lam_b(theta))
    nu, phi, r, d = p.nu, p.phi_m, p.r, p.delta
    kappa = nu * phi / (r + nu)

    def vstar_m_of(Vc, c):
        return Vc + U_m + (c + p.kappa_b + r * (Vc + U_m)) / (nu * phi)

    V = np.empty_like(p.c)
    for i, c in enumerate(p.c):
        def Phi(Vc):
            return (-c + d * p.w0 + lam_s * p.beta * (pit * float(G.H(Vc + U_c)) + (1 - pit) * kappa * float(G.H(vstar_m_of(Vc, c))))
                    - (r + d) * Vc)
        lo = (-c + d * p.w0) / (r + d) - 1e-9
        hi = max((G.hi if np.isfinite(G.hi) else G.mean() + 50.0) - min(U_c, U_m) + 1e-9, lo + 1.0)
        guard = 0
        while Phi(hi) >= 0 and guard < 50:
            hi += 1.0; guard += 1
        V[i] = optimize.brentq(Phi, lo, hi, xtol=p.root_xtol, maxiter=500)
    vs_c = V + U_c
    vs_m = np.array([vstar_m_of(V[i], c) for i, c in enumerate(p.c)])
    acc_c = np.asarray(G.surv(vs_c)); acc_m = np.asarray(G.surv(vs_m))
    eta_c = lam_s * pit * acc_c; eta_m = lam_s * (1 - pit) * acc_m
    n = p.s * p.f / (eta_c + phi * eta_m + d); S = n.sum(); N = n / S
    R_c = r * U_c - (-p.kappa_b + lam_b * (1 - p.beta) * float((N * np.asarray(G.H(vs_c))).sum()))
    R_m = r * U_m - (-p.kappa_b + lam_b * (1 - p.beta) * kappa * float((N * np.asarray(G.H(vs_m))).sum()))
    Q = float(((eta_c + phi * eta_m) * n).sum()); b = Q
    xi_c = lam_b * float((N * acc_c).sum()); xi_m = lam_b * float((N * acc_m).sum())
    B_c = p.pi * b / max(xi_c, 1e-300); B_m = (1 - p.pi) * b / max(phi * xi_m, 1e-300)
    pit_cons = B_c / (B_c + B_m) if (B_c + B_m) > 0 else p.pi
    R_theta = theta - (B_c + B_m) / S
    return dict(lam_s=lam_s, lam_b=lam_b, V=V, vs_c=vs_c, vs_m=vs_m, eta_c=eta_c, eta_m=eta_m, n=n, S=S, N=N, Q=Q, B_c=B_c, B_m=B_m,
                acc_c=acc_c, acc_m=acc_m, xi_c=xi_c, xi_m=xi_m, R_c=R_c, R_m=R_m, R_pit=pit_cons - pit, R_theta=R_theta, pit_cons=pit_cons)


def solve_equilibrium_financing(p: Params, x0: Optional[np.ndarray] = None) -> FinancingEquilibrium:
    """Stationary equilibrium of the financing extension as a 3-equation root problem in (log theta, U_c, logit pi_tilde),
    with U_m = (k - pi U_c)/(1 - pi) imposed by ex-ante free entry. Equations: buyer HJBs for cash and mortgage buyers and
    stock-consistency of the cash share of meetings. (The theta = B/S identity is implied by the flow accounting and is
    reported as a diagnostic.)"""
    pi = p.pi

    def unpack(x):
        theta = float(np.exp(x[0])); U_c = float(x[1]); pit = float(1 / (1 + np.exp(-x[2])))
        U_m = (p.k - pi * U_c) / (1 - pi)
        return theta, U_c, U_m, pit

    def F(x):
        theta, U_c, U_m, pit = unpack(x)
        st = _fin_state(theta, U_c, U_m, pit, p)
        return np.array([st["R_c"] / p.r, st["R_m"] / p.r, st["R_pit"]])

    if x0 is None:
        eq0 = solve_equilibrium(p.copy(phi=pi + (1 - pi) * p.phi_m))
        x0 = np.array([np.log(eq0.theta), p.k, np.log(pi / (1 - pi))])
    best = None
    for method in ("hybr", "lm"):
        try:
            sol = optimize.root(F, x0, method=method, tol=1e-12, options={"maxfev": 4000} if method == "hybr" else {"maxiter": 4000})
        except (ValueError, FloatingPointError, ZeroDivisionError, OverflowError):
            continue
        if best is None or np.abs(sol.fun).max() < np.abs(best.fun).max():
            best = sol
        if np.abs(sol.fun).max() < 1e-9:
            break
    sol = best
    theta, U_c, U_m, pit = unpack(sol.x)
    st = _fin_state(theta, U_c, U_m, pit, p)
    G = p.G; nu, phi, r = p.nu, p.phi_m, p.r
    V = st["V"]
    Delta_b = ((1 - p.beta) * (p.c + r * V) - p.beta * (p.kappa_b + r * U_m)) / (nu * phi)
    outside = p.beta * (U_c - U_m)
    Mc = np.asarray(G.M(st["vs_c"])); Mm = np.asarray(G.M(st["vs_m"]))
    pbar_c = V + p.beta * Mc
    pbar_m = V + p.beta * (st["vs_m"] + Mm - V - U_m) + Delta_b
    sel = p.beta * ((st["vs_m"] + Mm) - (st["vs_c"] + Mc))
    m_pending = st["eta_m"] * st["n"] / nu
    Q_c = st["eta_c"] * st["n"]; Q_m = phi * st["eta_m"] * st["n"]
    resid = {"R_c": st["R_c"], "R_m": st["R_m"], "R_pit": st["R_pit"], "free_entry": pi * U_c + (1 - pi) * U_m - p.k,
             "theta_identity": st["R_theta"], "max_abs": float(np.abs(sol.fun).max())}
    converged = bool(np.abs(sol.fun).max() < 1e-8)
    return FinancingEquilibrium(theta, st["lam_s"], st["lam_b"], U_c, U_m, pit, V, st["vs_c"], st["vs_m"], st["eta_c"], st["eta_m"],
                                st["n"], m_pending, st["S"], st["B_c"], st["B_m"], st["Q"], Q_c, Q_m, pbar_c, pbar_m, Delta_b, sel,
                                outside, resid, converged, int(sol.nfev))


# ----------------------------------------------------------------------------------------------
# Transition dynamics (Experiment 9): unanticipated shock to A or phi at t = 0
# ----------------------------------------------------------------------------------------------
@dataclass
class Transition:
    t: np.ndarray
    theta: np.ndarray
    V: np.ndarray          # (T, n_c)
    eta: np.ndarray
    n: np.ndarray
    pbar: np.ndarray
    Q: np.ndarray
    S: np.ndarray
    P_obs: np.ndarray      # transaction-weighted mean price
    P_fix: np.ndarray      # fixed (initial) composition index
    omega: np.ndarray
    share_urgent_sales: np.ndarray
    mean_tom_sold: np.ndarray
    converged: bool
    iterations: int
    max_update: float


def transition(p_old: Params, p_path, T: Optional[float] = None, dt: Optional[float] = None) -> Transition:
    """p_path(t) -> Params in force at time t (piecewise constant shocks). Terminal condition = steady state of p_path(T)."""
    T = p_old.T if T is None else T
    dt = p_old.dt if dt is None else dt
    nt = int(round(T / dt)) + 1
    t = np.linspace(0.0, T, nt)
    eq0 = solve_equilibrium(p_old)
    pT = p_path(T)
    eqT = solve_equilibrium(pT)
    n_c = len(p_old.c)
    theta = np.full(nt, eqT.theta)
    # parameter paths
    A_t = np.array([p_path(ti).A for ti in t]); phi_t = np.array([p_path(ti).phi for ti in t])
    pp = pT  # structural constants other than A, phi are assumed unchanged along the path
    G = pp.G; r, d, beta, k = pp.r, pp.delta, pp.beta, pp.k
    V = np.tile(eqT.sellers.V, (nt, 1)); n = np.tile(eq0.n, (nt, 1)); eta = np.zeros((nt, n_c))
    converged = False; it = 0; upd = np.inf
    for it in range(1, pp.dyn_max_iter + 1):
        # backward HJB (semi-implicit): V_j = [V_{j+1} + dt(-c + d w0 + Lam_s beta H(V_{j+1}+k))]/(1+(r+d)dt)
        V[-1] = eqT.sellers.V
        for j in range(nt - 2, -1, -1):
            Lam_s = phi_t[j] * A_t[j] * theta[j] ** pp.alpha
            V[j] = (V[j + 1] + dt * (-pp.c + d * pp.w0 + Lam_s * beta * np.asarray(G.H(V[j + 1] + k)))) / (1 + (r + d) * dt)
        # forward KFE with hazards from V_j
        n[0] = eq0.n
        for j in range(nt):
            Lam_s = phi_t[j] * A_t[j] * theta[j] ** pp.alpha
            eta[j] = Lam_s * np.asarray(G.surv(V[j] + k))
            if j < nt - 1:
                n[j + 1] = (n[j] + dt * pp.s * pp.f) / (1 + dt * (eta[j] + d))   # implicit Euler
        # free entry pointwise: r k + kappa_b = phi lam_b(theta) (1-beta) sum H(V_j + k) N_j
        theta_new = np.empty(nt)
        for j in range(nt):
            N = n[j] / n[j].sum(); Hj = float((N * np.asarray(G.H(V[j] + k))).sum())
            rhs = (r * k + pp.kappa_b) / (phi_t[j] * A_t[j] * (1 - beta) * max(Hj, 1e-300))
            # lam_b = theta^(alpha-1) => theta = rhs^(1/(alpha-1))
            theta_new[j] = rhs ** (1.0 / (pp.alpha - 1.0))
        upd = float(np.abs(theta_new - theta).max())
        theta = (1 - pp.damping) * theta + pp.damping * theta_new
        if upd < pp.dyn_tol:
            converged = True
            break
    # outputs
    vstar = V + k
    M = np.asarray(G.M(vstar)); pbar = V + beta * M
    Q = (eta * n).sum(axis=1); S = n.sum(axis=1)
    omega = eta * n / Q[:, None]
    P_obs = (omega * pbar).sum(axis=1)
    P_fix = (eq0.omega[None, :] * pbar).sum(axis=1)            # pre-shock (predetermined) transaction composition
    urgent = p_old.c >= np.median(p_old.c)
    share_urgent = omega[:, urgent].sum(axis=1)
    mean_tom = (omega / (eta + d)).sum(axis=1)
    return Transition(t, theta, V, eta, n, pbar, Q, S, P_obs, P_fix, omega, share_urgent, mean_tom, converged, it, upd)