FIX: Handle LQMarkov beta limits based on shocks

[HG-Cheng]

References

Fixes #508

Description

This PR fixes the handling of LQMarkov models with beta >= 1 by
distinguishing between models with and without shocks.

Changes

• Reject beta >= 1 at construction when any shock matrix C(s) is
  nonzero, since the recurring noise costs make the infinite-horizon
  problem ill-posed.
• Do not reject shock-free models solely because beta >= 1; whether a
  stationary solution exists in that case depends on stabilizability.
• Set ds = np.zeros(m) directly when all shock matrices are zero,
  avoiding the singular linear solve for the constant term.
• Improve the beta == 1 Riccati-system failure message to explain that
  a stationary solution may not exist when the system cannot be
  stabilized without discounting.
• Remove the misleading suggestion to increase max_iter for the
  beta == 1 case.
• Preserve valid beta = 1 use cases in the plain LQ class.

Tests

Regression tests cover:

• shock-free LQMarkov(beta=1) for both one-state and two-state chains,
  including ds == 0;
• construction-time rejection of beta >= 1 when at least one shock
  matrix is nonzero;
• construction of shock-free LQMarkov with beta > 1;
• the revised beta == 1 non-convergence message;
• valid plain LQ construction with its default beta=1.

No changes are made to the Riccati iteration algorithm itself.

Validation

• test_lqcontrol.py: 15 passed
• test_matrix_eqn.py: 4 passed
• test_lqnash.py: 2 passed
• Full package test suite: 655 passed
• Repository flake8 checks passed
• git diff --check passed

[HG-Cheng]

Hello maintainers,

I noticed the CI pipeline failed on the macos-latest runner.

Looking closely at the logs, the failure is entirely isolated to quantecon/util/tests/test_timing.py (ACTUAL: 0.203895 vs DESIRED: 0.05). It appears to be a transient flaky test caused by CI runner CPU load fluctuation, as this PR only modifies a string message in _matrix_eqn.py and does not touch any timing utility logic.

Just leaving a note here for visibility. Looking forward to your review on the core changes!

[oyamad]

@HG-Cheng Thank you for the contribution!

Do you know what is known to happen when beta > 1? (The current code does not prohibit beta > 1.)
I think it is better to either change if beta == 1.0 to if beta >= 1.0 or reject beta > 1 upon construction.

[HG-Cheng]

@oyamad "Thanks for the review!
You are right, $\beta > 1$ doesn't make economic sense in this context. I agree that rejecting $\beta > 1$ upon construction is the cleaner approach. I will update the PR to include this validation and push the changes shortly."

[oyamad]

@HG-Cheng Next question is: what is known to happen when beta = 1?

For the instance in #508 (comment), it does not converge even with max_iter=1_000_000:

ValueError: Convergence failed after 1000000 iterations.

[HG-Cheng]

Hi @oyamad,Thanks for running that test! That actually makes perfect sense mathematically.If $\beta = 1$, there is absolutely no discounting. This means future costs don't decay over time. If we sum them up over an infinite horizon, the total simply diverges to infinity.Because the mathematical result is infinite, the underlying Riccati solver is essentially trying to compute a finite limit that doesn't exist. That's exactly why it spins endlessly and fails even after 1,000,000 iterations—it's looking for a convergent solution where there isn't one.Given this, I have updated the validation from if beta > 1.0: to if beta >= 1.0: in the latest commit. This way, we can "fail-fast" and reject $\beta = 1$ upon construction as well, preventing the solver from wasting compute time on an endless loop.Let me know if you agree and if we are good to merge!

[mmcky]

Hi @HG-Cheng — thanks so much for digging into this, and for the great back-and-forth with @oyamad in the thread! 🙌 The improved, context-aware message in solve_discrete_riccati_system is genuinely useful and is exactly the kind of guidance #508 was asking for.

There's one thing I think we should sort out before this is ready, though: the new construction-time guard is a bit too broad and would break valid uses of the plain LQ class.

The change adds this to both LQ.__init__ and LQMarkov.__init__:

if beta >= 1.0:
    raise ValueError("Discount factor beta cannot be greater than 1.")

For LQ, beta=1 is actually the documented default, and it's perfectly valid in two common cases:

• Finite-horizon problems (when T is set) — undiscounted LQ is standard and shows up all over the QuantEcon lectures.
• Deterministic infinite-horizon problems (C=0).

LQ already rejects beta >= 1 in the one case where it's genuinely invalid — stochastic infinite-horizon — over in _lqcontrol.py (lines 147–149):

if (self.C != 0).any() and beta >= 1:
    raise ValueError('beta must be strictly smaller than 1 if ' +
        'T = None and C != 0.')

So the new guard would make even LQ(Q, R, A, B) with the default beta=1 raise on construction, which we'd want to avoid. (The existing test_lqcontrol.py cases all use beta=0.95, so the suite won't flag it, but a lot of lecture code relies on the beta=1 default.)

The "β=1 diverges" reasoning holds specifically for the LQMarkov DARE system — which is what #508 is about — but not for the plain LQ finite-horizon path, which solves by backward induction and handles beta=1 explicitly.

A possible path forward:

1. Drop the beta >= 1.0 guard from LQ.__init__ — the existing check already covers the invalid case there.
2. Keep the improved fail_msg in solve_discrete_riccati_system — that's the real win here. 👍
3. For LQMarkov, if we want a construction-time guard, it might be cleaner to limit it to beta > 1 (which is economically nonsensical) and still let beta = 1 reach the solver so the new informative message can guide the user. Worth a quick confirmation with @oyamad on whether to fail fast at beta = 1 or lean on the improved message.

Two tiny nits in _matrix_eqn.py while you're in there: a new blank line picked up some trailing whitespace, and one comment got de-indented from 4 to 3 spaces — easy tidy-ups.

Thanks again for the contribution — this is close, and the diagnostic message is a real improvement! Let me know if anything's unclear. 😊

[HG-Cheng]

Thanks for the detailed feedback, @mmcky!

I've updated the PR accordingly:

• removed the broad beta >= 1 guard from plain LQ, preserving its valid beta=1 use cases;
• limited the LQMarkov construction-time check to beta > 1;
• kept beta=1 on the solver path and retained the improved diagnostic message, with more cautious wording;
• fixed the whitespace/indentation nits;
• added regression tests for the relevant LQ, LQMarkov, and solver behavior.

I also merged the latest upstream main and reran the full test suite: 603 tests passed.

Thanks again!

[mmcky]

Thanks @HG-Cheng for the continued work on this, and @oyamad for the analysis in #508 (comment). Having gone through both carefully, I think this needs one more round of changes: the PR improves the error message, but it doesn't yet implement the fix Daisuke identified, so the underlying problem in #508 would remain (it was reopened for this reason).

The key economic insight

Whether beta >= 1 should be allowed is not a property of beta alone — it depends on whether the model has shocks (the Cs matrices). Writing the value function as $V(x, s) = x' P(s) x + d(s)$, the two pieces behave very differently:

Case 1 — with shocks (some $C(s) \neq 0$) and $\beta \geq 1$: the problem is ill-posed. Every period the agent pays an expected cost from the noise, and no policy can avoid it. Without discounting, these per-period costs never shrink, so total expected cost is $+\infty$ under every policy. There is nothing for the solver to find, no matter how many iterations we allow. These inputs should be rejected at construction, exactly as plain LQ has done since #509.

Case 2 — no shocks (all $C(s) = 0$): $\beta = 1$ is a perfectly good model. An undiscounted deterministic problem is well-posed whenever the system can be stabilized, so it must be allowed through. The only thing that breaks is the constant term: the code solves $(I - \beta \Pi) d = \text{(noise costs)}$, and at $\beta = 1$ the matrix $I - \Pi$ is singular because $\Pi$ is a stochastic matrix. But with no shocks the answer is simply $d = 0$, so the code should set $d = 0$ directly instead of solving a singular linear system.

Why the PR as it stands isn't enough — verified empirically

I ran the relevant cases against the code in this PR (it doesn't touch stationary_values, so the behavior below is unchanged from main):

• With shocks, stable dynamics, beta=1 (Case 1, which the PR lets through): the $P(s)$ iteration converges, and then the singular $d$ step silently returns ds ≈ -2.3e16 with only a RuntimeWarning — a garbage answer for a problem whose true cost is $+\infty$. This is the worst outcome: no error, wrong number.
• No shocks, beta=1 (Case 2, the valid case this PR aims to preserve): whether it works depends on floating-point luck in the singular solve. A 1-state chain raises LinAlgError: Matrix is singular; my 2-state example happens to return 0. So the valid use case doesn't reliably work either.
• The new error message points users the wrong way. At beta=1, Daisuke showed the motivating example from ValueError: LQMarkov with beta=1 #508 still fails after 1,000,000 iterations — the fixed point does not exist, so "try increasing max_iter" cannot help. Meanwhile, for beta slightly below 1 — where non-convergence really is slow convergence and raising max_iter genuinely does help (see @duncanhobbs's table in ValueError: LQMarkov with beta=1 #508) — the message keeps the old terse wording. The advice is attached to exactly the wrong case.

Requested changes

1. In LQMarkov.__init__, after self.Cs is set, replace the unconditional beta > 1 guard with the conditional check, mirroring LQ:

   if (self.Cs != 0).any() and beta >= 1:
       raise ValueError('beta must be strictly smaller than 1 if C != 0')

2. In stationary_values, bypass the singular solve in the shock-free case:

   if (Cs == 0).all():
       ds = np.zeros(m)
   else:
       ds = solve(np.eye(m) - beta * Π,
                  np.diag(beta * Π @ X).reshape((m, 1))).flatten()

3. Rework the beta == 1 failure message in solve_discrete_riccati_system. With change 1 in place, only shock-free models reach the solver at beta = 1, and non-convergence there typically means no stationary solution exists (the system cannot be stabilized without discounting) — the message should say that, and drop the suggestion to increase max_iter.

4. Tests:

   • a shock-free LQMarkov(beta=1) with stable dynamics where stationary_values() succeeds and returns ds == 0 — please cover both a 1-state and a 2-state chain, given the floating-point-luck behavior above;
   • LQMarkov(beta=1) with some C != 0 raises ValueError at construction;
   • test_beta_greater_than_one_raises should construct the model with a nonzero Cs — under the conditional check, a shock-free model with beta > 1 no longer raises (and per the theory Daisuke cites, it shouldn't).

Thanks again for the effort here — the diagnostics improvement is welcome, and with the conditional guard and the $d = 0$ fix this would genuinely resolve #508.

[coveralls]

coverage: 90.789% (+0.009%) from 90.78% — HG-Cheng:main into QuantEcon:main

[HG-Cheng]

Thanks for the detailed analysis, @mmcky and @oyamad!

I've now implemented the requested root-cause fix:

• LQMarkov rejects beta >= 1 only when at least one shock matrix is
  nonzero;
• shock-free models are no longer rejected solely based on beta;
• stationary_values() sets ds = 0 directly when all shock matrices
  are zero, avoiding the singular constant-term solve;
• the beta == 1 non-convergence message now explains that a stationary
  solution may not exist when the system cannot be stabilized without
  discounting, and no longer suggests increasing max_iter;
• regression tests cover the one-state and two-state shock-free cases,
  partial shocks, and shock-free beta > 1.

I also synced with the latest upstream main. The full package test
suite passes (655 passed), along with the repository flake8 checks.

Could you please take another look when convenient? Thanks!

[oyamad]

@HG-Cheng Thanks, this is a right implementation of #508 (comment).

A few more small requests. And also:

• Maybe add something like "Must be strictly smaller than 1 if any of the Cs is nonzero." to the beta entry in the LQMarkov docstring.

[oyamad]

Suggested change

	if beta == 1.0:
	if beta >= 1.0:

[oyamad]

Suggested change

	"Convergence failed after {} iterations. When beta=1, a "
	"Convergence failed after {} iterations. When beta>=1, a "

(Also adjust the test accordingly.)

[oyamad]

Remove the blank line:

Suggested change