FIX: correct the Markov-jump Riccati recursion in LQMarkov (#982)

[mmcky]

This fixes the Markov-jump Riccati recursion in solve_discrete_riccati_system, which LQMarkov.stationary_values uses (#982).

The control $u_t$ is chosen before $s_{t+1}$ is realized, so the expectation over next-period states belongs inside the inverse. For each state $i$, the loop now forms $\bar P_i = \sum_j \Pi_{ij} P_j$ and applies the single-state update

$$ P_i = R_i + \beta A_i' \bar P_i A_i - (\beta A_i' \bar P_i B_i + N_i') (Q_i + \beta B_i' \bar P_i B_i)^{-1} (\beta B_i' \bar P_i A_i + N_i). $$

The F and d formulas in LQMarkov.stationary_values were already correct given Ps, so they are unchanged. The function's docstring now states the equation it solves.

One deliberate difference from the snippet in #982. The left factor $\beta A_i' \bar P_i B_i + N_i'$ is formed directly, not as (beta * Bs[i].T @ P_bar @ As[i] + Ns[i]).T. The two are equal only when $\bar P_i$ is exactly symmetric. In the transposed form, the rounding-level asymmetric part of Ps evolves as $E \mapsto \beta(A'\bar E A + (BF)'\bar E(BF))$, which has no closed-loop factor. It can then grow even though the optimal closed loop is stable, for example when $A$ is unstable. The iteration then diverges, or returns NaN, on well-posed problems that main solves, including single-regime ones. A comment in the code records this.

Behaviour change

LQMarkov.stationary_values returns different Ps, ds and Fs, and so compute_sequence returns different paths, whenever a row of $\Pi$ has more than one nonzero entry and $Q_i + \beta B_i' P_j B_i$ or $B_i' P_j A_i$ varies with $j$. On the test fixtures, the largest relative changes are 0.45% in P, 0.40% in d and 0.21% in F. On the example in #982, F changes by 1.9%. Problems with one regime, identical regimes, a permutation $\Pi$ or $B = 0$ are unaffected.

Tests

• Regenerated reference values. The values in test_stationary_scalar, test_stationary_mat and the two sequence tests came from the old recursion. The new values were reproduced independently by Howard policy iteration with exact policy evaluation. The closest uses 5.2% of its assertion's tolerance.
• New test_stationary_values_are_policy_costs. It checks that Ps equal the exact cost of following the returned Fs on the two existing fixtures and on the BUG: LQMarkov / solve_discrete_riccati_system averages over next-period states outside the inverse #982 example. On main it fails, with a relative deviation of 3.8e-3 against rtol=1e-6. Here it passes, with a relative error of about 1e-11.
• New test_identical_regimes_match_lq. With one regime, or two identical regimes, and an unstable $A$, LQMarkov must match LQ. It passes with about a 12,000× margin, and fails if the left factor is written as a transpose.

Verification

• I derived the Bellman equation from the model in the LQMarkov docstring and checked the code against it line by line, including F and d. On random models with m = 3 and nonzero N and C, I minimized the right-hand side directly: the minimizer is -F x and the minimum is x'Px + d.
• A Monte Carlo simulation of the documented model (400k paths) agrees with x'Px + d within sampling error. The old values are off by 11 to 244 standard errors.
• I ran a random sweep of 1,220 solvable problems with m from 1 to 4. The solver converges on all of them except two, which need more than max_iter=1000 iterations at β = 0.99, the same as any correct variant. On all 600 single-regime problems it agrees with LQ, as main does.
• Locally, the full test suite passes except for the three fetch_nb_dependencies tests, which need network access my environment doesn't have.

Once this is released, the notes about the discrepancy in markov_jump_lq (QuantEcon/lecture-python-advanced.myst#392) can be removed.

Closes #982

🤖 Generated with Claude Code

[Copilot]

🔵 Needs a closer look

It is a behaviour-changing numerical correctness fix to a core published library function with regenerated golden values whose matrix-case numbers could not be fully re-verified here, warranting expert human sign-off.

0 open findings

What changed in this PR

This PR fixes a correctness bug (#982) in solve_discrete_riccati_system, the Markov-jump Riccati solver used by LQMarkov.stationary_values. The previous recursion averaged the quadratic correction term over next-period Markov states outside the matrix inverse; since the control u_t is chosen before s_{t+1} is realized, the expectation P̄_i = Σ_j Π_ij P_j must enter inside the inverse. The fix makes LQMarkov return value matrices that are actually the cost of the returned policy and decision rules that are optimal.

Changes:

• Rewrote the per-state Riccati update in _matrix_eqn.py to form P̄_i once per sweep and apply the single-state update with the expectation inside the inverse; the left factor is formed directly (not as a transpose) to avoid amplifying rounding asymmetry.
• Updated the function docstring with the exact equation solved and corrected the Ps return type from ndim=2 to ndim=3.
• Regenerated the golden values in the scalar/matrix/sequence tests and added two regression tests (test_stationary_values_are_policy_costs, test_identical_regimes_match_lq) plus a _lq_markov_policy_cost helper.

File	Description
quantecon/​_matrix_eqn.py	Corrects the stacked Riccati recursion (expectation inside the inverse), documents the solved equation, and fixes the return-type docstring.
quantecon/​tests/​test_lqcontrol.py	Regenerates reference values to match the corrected recursion and adds policy-cost and identical-regime regression tests with a supporting helper.

🧠 Review effort: Balanced

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[coveralls]

coverage: 94.085% (-0.003%) from 94.088% — fix-lqmarkov-riccati-982 into main

[mmcky]

Leaving this open for now. #983 by @jadhavgaurav, opened the day before this one, makes the same core fix. Its code also builds the left factor directly, so the problem with the transposed form described above doesn't apply to it, and the two PRs give the same results up to rounding.

So one option is to push the extras from this PR onto #983:

• the docstring stating the solved equation, and the ndim=3 fix;
• the comment on why the left factor is formed directly;
• the _lq_markov_policy_cost helper and the two new tests.

They carry over with only the comment adapted to #983's inline expression. This PR's test file passes unchanged against #983's code, and test_identical_regimes_match_lq still catches the transposed form there. I've asked on #983 for the QEP-5 line first.

[mmcky]

Update: @jadhavgaurav offered to port this PR's extras (the docstring, the ndim=3 fix, the left-factor comment, the helper and the two tests) to #983, and we've taken up the offer (details). Since #983 came first with the same fix, this PR is pending close in favour of it.

[mmcky]

Closing in favour of #983, which has the same fix plus this PR's docstring and tests. It's released in v0.13.1.