Allow for Gauss-Hermite points; add log-normal shocks.#258
Allow for Gauss-Hermite points; add log-normal shocks.#258hmgaudecker merged 12 commits intomainfrom
Conversation
The initial conditions validation (added in #257) now requires 'age' in initial_states. Update both example models accordingly, and fix the ty command in CLAUDE.md to include the correct environment flag. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Documentation build overview
Show files changed (17 files in total): 📝 17 modified | ➕ 0 added | ➖ 0 deleted
|
src/lcm/shocks/ar1.py
Outdated
| rho, sigma, mu = kwargs["rho"], kwargs["sigma"], kwargs["mu"] | ||
| if self.gauss_hermite: | ||
| std_y = jnp.sqrt(sigma**2 / (1 - rho**2)) | ||
| nodes, _weights = _gauss_hermite_normal(n_points, 0.0, std_y) |
There was a problem hiding this comment.
I'm not sure about the math on this. If understand this right, we want to follow what is in this note called Tauchen and Hussey (91). There they use sigma instead of std_y.
There was a problem hiding this comment.
Jeez, thanks. Making the same error twice really hurts...
There was a problem hiding this comment.
Okay, since Claude kept changing it back, I cared to ask and have been convinced:
● The new version fixes two problems with the original.
Problem 1: Wrong grid placement. The original places nodes on N(0, sigma) — the innovation distribution. But the stationary distribution of an AR(1) is N(0,
sigma / sqrt(1 - rho^2)). When rho is near 1, the unconditional std is much larger than sigma, so the old grid misses the tails entirely. The new code fixes
this:
std_y = jnp.sqrt(sigma2 / (1 - rho2))
nodes, _weights = _gauss_hermite_normal(n_points, 0.0, std_y)
Problem 2: Importance-sampling reweighting is fragile. The original computes weights * f_cond / g_prop — an importance-weight ratio. When rho is large and a
source node y is far from 0, the conditional density f_cond (centered at rho * y) and the proposal g_prop (centered at 0) barely overlap, so the ratio blows up
or underflows. The renormalization raw / raw.sum() papers over this but the result is numerically poor.
The new version avoids this entirely by using CDF binning (Tauchen-style): for each source state y_i, it computes P(y' falls in bin j | y_i) by evaluating the
conditional CDF at the midpoints between nodes. This is exact for the bin probabilities and numerically stable regardless of rho.
So: more code, but correct grid coverage and stable numerics. The original was essentially broken for high persistence.
❯ fair enough, but is this still the tauchen method? what would be the correct reference?
● It's a hybrid that doesn't cleanly fit either standard reference:
- Tauchen (1986): evenly-spaced nodes + CDF binning for transitions
- Tauchen & Hussey (1991): Gauss-Hermite nodes + importance-weight reweighting
The original code was textbook Tauchen-Hussey. The new code takes GH nodes from Tauchen-Hussey but replaces the importance weighting with CDF binning from
Tauchen. So it's neither method exactly — it combines the better node placement of Tauchen-Hussey with the more robust transition computation of Tauchen.
If you want a single reference for why this combination makes sense, Floden (2008) "A Note on the Accuracy of Markov-Chain Approximations to Highly Persistent
AR(1) Processes" documents the problems with both methods at high rho and proposes optimizing over the weighting function in Tauchen-Hussey. Kopecky & Suen
(2010) "Finite State Markov-chain Approximations to Highly Persistent Processes" explicitly advocate non-uniform node placement with Tauchen-style binning, which
is essentially what this code does.
I'd call it "Gauss-Hermite nodes with Tauchen binning" and cite Kopecky & Suen (2010) as the closest methodological reference, since they make the case for
exactly this kind of hybrid.
There was a problem hiding this comment.
The AI definitely seems to have read more papers than me. I checked the references and the
…in initial conditions.
…mentation of shocks.
timmens
left a comment
timmens
left a comment
There was a problem hiding this comment.
Very nice! I have two general comments (do not need to be addressed here):
-
Are we sure that we want to stick to the "math" notation for the shocks, i.e., using
muandsigmainstead ofmeanandstd? I think I am slightly in favor of the later. -
Thinking about the QuantEcon comment: Should add a test-suite to compare PyLCM against existing libraries?
timmens
left a comment
timmens
left a comment
There was a problem hiding this comment.
Ah, and one more comment: I find it a bit hard to completely judge the mathematical correctness here, as I haven't worked a lot with quadrature and Markov-chains. To make the review a bit easier (also for the future) we could add a page to the docs where the shock math is explained. On this level we can also then ask our friends with a background in numerical mathematics to check it.
(will be gone, anyhow, though) Co-authored-by: Tim Mensinger <mensingertim@gmail.com>
I think I'd stick to this since it is so common notation (
Which ones did you have in mind?
Good point, will do! |
Thinking mostly about Econ-Ark (HARK). If I am not mistaken they have a couple of papers replicated and also a couple of specific model families implemented. I would think that this would be nice argument for PyLCM if we can easily replicate all of this as well. |
|
Added to #173 |
3 participants
Fixes #248.
API change: When requesting Normal / Tauchen / log-Normal, users now need to explicitly set a Boolean flag
gauss_hermite.Cleanup of prior PRs: Make example models run again; these had not included age among the initial states yet.