Add Monte Carlo vs Transition Matrix example notebooks and infrastructure#1750
Draft
Add Monte Carlo vs Transition Matrix example notebooks and infrastructure#1750
Conversation
- Add MC sampling statistics (mean, std, SE) to aggregate comparison - Add shaded confidence band (mean ± 2 std) to MC vs TM time-series plot - Rewrite "Precision vs Accuracy" section with bias-variance framing - Add quantitative MSE decomposition: Bias² + Variance for both methods - Add granular timing breakdown (grid setup / TM build / eigensolve) for both 2D and Harmenberg cases - Add tail-mass diagnostic for grid truncation checking - Add practical grid-tuning guidance section with rules of thumb - Fix 5 pre-existing malformed notebook outputs (missing stream 'name') Made-with: Cursor
Directory structure: - 9 pedagogical notebooks (01-09) building from basic Markov TM through Harmenberg, Krusell-Smith, and sequence-space Jacobians - mathematical-framework.ipynb: unified theory (perch notation, bias-variance tradeoff, Markov/Harmenberg/KS/SSJ extensions) - LESSONS-LEARNED.md: bugs, gotchas, and proposed HARK improvements - bibliography.md/.bib: annotated bibliography of simulation methods - README.md: reading order and file inventory - _archive/: scaffolding documents from development (plans, context, reference analyses) preserved for historical reference Made-with: Cursor
# Transition_Matrix_Example.ipynb — Change Log
Changes made to Will Du's `Transition_Matrix_Example.ipynb` notebook,
with diagnosis and rationale for each.
---
## 1. TM-initialized Monte Carlo via inverse-CDF sampling (cell 55)
**Problem:** The original notebook initialized the MIT-shock Monte Carlo
simulation using HARK's default `initialize_sim()`, which starts all agents
with near-zero assets and `pLvl=1`. This means the MC starts far from
steady state, requiring a long burn-in before the MC and TM paths can be
meaningfully compared.
**Fix:** After `initialize_sim()`, we overwrite agent states with draws from
the TM ergodic distribution. Specifically, we compute the CDF of the 1-D
Harmenberg ergodic distribution over normalized market resources, then use
piecewise-linear inverse-CDF sampling to draw continuous `mNrm` values for
each of the MC agents. Asset positions `aNrm` are interpolated from
the steady-state asset policy, and `pLvl` is set to 1 (Harmenberg neutral
measure). This ensures both MC and TM start from the identical steady-state
distribution, so any remaining discrepancy is purely MC sampling noise vs TM
discretization error.
**Cells affected:** 55
---
## 2. Newborn transitory-shock suppression (`t_age` fix) (cell 55)
**Problem:** After `initialize_sim()`, all agents have `t_age=0`, marking
them as "newborns." HARK's `get_shocks()` method
(`ConsIndShockModel.py`, lines 2191–2193) forces `TranShk=1.0` for all
newborn agents when the parameter `NewbornTransShk` is `False` (the
default). Since every agent is a newborn in the first simulated period,
*all* agents receive `TranShk=1.0` — completely suppressing transitory
income variation. This created a systematic +0.7% bias in first-period
aggregate consumption that decayed over ~3 periods.
**Diagnosis:** Inspecting `FinHorizonAgent.history['TranShk'][0,:]` revealed
all values were exactly 1.0, despite the `IncShkDstn` containing distinct
transitory shock values. Traced to:
```python
# HARK/ConsumptionSaving/ConsIndShockModel.py, lines 2191-2193
if not NewbornTransShk:
TranShkNow[newborn] = 1.0
```
**Fix:** After overriding `state_now` with the TM-sampled initial
conditions, set `FinHorizonAgent.t_age = np.ones(N_fin, dtype=int)` so
agents are not treated as newborns and draw proper transitory shocks from the
first period onward.
**Cells affected:** 55
---
## 3. Permanent income grid range (`max_p_fac`) (cell 13)
**Problem:** The default `max_p_fac=30` in `define_distribution_grid()`
creates a permanent income grid extending from ~10⁻¹⁰ to ~7.7 × 10⁹ (20
orders of magnitude). When computing aggregate assets in levels
(`grida = p × a`), the asset weights at extreme p-grid points reach
7.6 × 10¹³. The ergodic distribution has effectively zero mass at these
extremes (negative values at the last few points — numerical artifacts), but
machine-epsilon perturbations (~10⁻¹⁶) introduced by each transition-matrix
multiplication get amplified by these enormous weights into visible aggregate
drift, manifesting as a declining trend in the TM asset path figure.
**Diagnosis:** The drift scales with `max(|grida|)`:
| `max_p_fac` | p_max | max(\|grida\|) | 200-step drift |
|-------------|---------|----------------|----------------|
| 30 (default)| 7.7e+9 | 7.6e+13 | −3.1e−04 |
| 20 | 3.9e+6 | 3.8e+10 | −1.4e−06 |
| 10 | 2.0e+3 | 1.9e+07 | +6.4e−10 |
| 5 | 4.5e+1 | 4.4e+05 | −9.2e−11 |
**Fix:** Changed call to
`define_distribution_grid(num_pointsP=110, max_p_fac=10.0)`.
This keeps the grid within exp(7.6) ≈ 2000, still covering 99.99%+ of mass.
The 110 grid points (yielding 221 total p-grid points) are concentrated
where they matter rather than wasted on points in the empty tail.
**Cells affected:** 13
---
## 4. Consistent plot colors across all figures
**Problem:** Monte Carlo and Transition Matrix lines swapped colors between
figures — MC was sometimes blue (first plotted) and sometimes orange (second
plotted), making cross-figure comparison confusing.
**Fix:** Defined three color constants in the imports cell:
```python
COLOR_MC = "tab:blue"
COLOR_TM = "tab:orange"
COLOR_HARM = "tab:green"
```
Applied `color=COLOR_MC`, `color=COLOR_TM`, and `color=COLOR_HARM`
explicitly to every `plt.plot()` call across all plotting cells. Also
standardized legend labels (e.g., inconsistent capitalization like
"transition matrix", "transition Matrix", " Transition Matrices" → uniform
"Transition Matrix" or "Transition Matrix (Harmenberg)").
For the grid-comparison figure (`[aggregate_assets_grid_comparison]`), TM
lines at different grid sizes use an orange gradient (`plt.cm.Oranges`) to
visually associate them with the TM family while remaining distinguishable.
**Cells affected:** 3, 19, 23, 25, 29, 31, 39, 42, 61, 63
---
## 5. Increased MC agent count
**Problem:** The original `AgentCount` (in `Dict`) was lower, producing
noisier MC paths that made it harder to assess whether MC–TM discrepancies
were real or sampling noise.
**Fix:** Set `AgentCount = 200000` in the calibration dictionary.
**Cells affected:** 5, 47
---
## 6. Descriptive figure tags and axis labels
**Problem:** Figures lacked machine-readable tags and some had minimal axis
labels / titles.
**Fix:** Added bracketed tags (e.g., `# [mc_vs_tm_asset_paths]`) as the
first line of each plotting cell for cross-referencing. Added or improved
titles, axis labels, and MC summary statistics (mean, std, ±2σ band) to the
asset-path comparison figure.
**Cells affected:** 19 and other plotting cells
---
## 7. MSE decomposition diagnostic (cell 21)
**Problem:** No quantitative comparison of MC vs TM accuracy was provided.
**Fix:** Added a bias–variance–MSE decomposition for aggregate assets,
treating the MC long-run mean as the reference. This makes explicit the
tradeoff: MC has zero bias but positive variance; TM has zero variance but
nonzero discretization bias.
**Cells affected:** 21
---
## 8. Timing instrumentation (cell 13)
**Problem:** No visibility into computational cost of TM steps.
**Fix:** Added `time.time()` calls around `define_distribution_grid`,
`calc_transition_matrix`, and `calc_ergodic_dist`, printing elapsed times.
**Cells affected:** 13
---
## 9. Grid size and agent count in figure legends
**Problem:** Figures did not indicate the TM grid resolution or MC agent
count, making it impossible to judge whether discrepancies were due to
insufficient grid points, insufficient agents, or genuine bugs.
**Fix:** Legend labels now include the relevant parameters, e.g.:
- `"Monte Carlo (200,000 agents)"`
- `"Transition Matrix (100 m-pts, 2D)"`
- `"Transition Matrix — Harmenberg (1000 m-points)"`
- `"TM Harmenberg (500 m-points)"` (grid-comparison figure)
Values are computed from the live objects (`len(example1.dist_mGrid)`,
`n_agents`, etc.) so they stay accurate if parameters change.
**Cells affected:** 19, 23, 25, 29, 31, 39, 42, 61, 63
---
## 10. Terminology: "MIT shock" → "anticipated deviation" / "perfect-foresight transition path"
**Problem:** The notebook described the finite-horizon interest-rate
experiment as a "Perfect foresight MIT shock" and labeled the output plots
"MIT Shock IRF." Both terms are inaccurate for the exercise being performed.
An "MIT shock" in the heterogeneous-agent literature refers to an
*unanticipated* aggregate disturbance — a zero-probability event that
surprises agents at the moment it occurs (Boppart, Krusell, and Mitman, 2018,
JEDC). An "impulse response function" traces the economy's dynamics from
that moment of surprise onward.
In this notebook, the interest rate deviation at $t = t^*$ is fully
anticipated from $t=0$: the finite-horizon problem is solved by backward
induction over the full $T$-period horizon, so consumption and savings
policies at *every* date reflect knowledge of the upcoming rate change. This
corresponds to what the SSJ literature calls the response to "s-period-ahead
news" (Auclert, Bardóczy, Rognlie, and Straub, 2021, *Econometrica*,
Section 2.2) or, equivalently, a nonlinear perfect-foresight transition path
(McKay, 2023, lecture notes).
**Fix:** Replaced cell 32 with a thorough description of the experiment,
citing Auclert et al. (2021), Boppart et al. (2018), McKay (2023), and
Harmenberg (2021). Updated section headers and plot titles to use
"Perfect-Foresight Transition" and "Anticipated Rate Deviation" instead of
"MIT Shock IRF."
**Cells affected:** 32, 60, 61, 62, 63
---
## 11. Terminal condition bug: `cFunc_terminal_` vs `solution_terminal` (cell 49)
**Problem:** The notebook set
`FinHorizonAgent.cFunc_terminal_ = deepcopy(ss.solution[0].cFunc)` to give the
finite-horizon agent the steady-state consumption function as its terminal
condition. However, `cFunc_terminal_` is a *parameter attribute* that does **not**
propagate to `solution_terminal`, which is the actual object HARK's backward-
induction solver uses. `solution_terminal` retained its default "consume
everything" behavior (`c(m) = m`), meaning the last period's agent consumed
all market resources.
**Diagnosis:** Comparing `FinHorizonAgent.cFunc_terminal_(m)` vs
`FinHorizonAgent.solution_terminal.cFunc(m)` showed that the former matched
the SS consumption function while the latter was the identity function.
The solver reads `solution_terminal`, not `cFunc_terminal_`.
**Fix:** Replaced `cFunc_terminal_` assignment with:
```python
agent.solution_terminal = deepcopy(ss.solution[0])
```
This sets the full terminal solution (cFunc, vFunc, vPfunc, hNrm, MPCmin,
etc.) correctly. Now encapsulated in the `create_finite_horizon_agent()`
helper.
**Impact:** Small improvement near the terminal period, but did not resolve the
main asset-level discrepancy (see Fix 12).
**Cells affected:** 3 (helper function), 49
---
## 12. Critical bug: neutral-measure `IncShkDstn` used to solve finite-horizon agent
**Problem:** The finite-horizon agent's consumption problem was solved using
the **Harmenberg neutral-measure** income shock distribution instead of the
**original (true)** income shock distribution. The Harmenberg neutral measure
reweights shock probabilities to efficiently track the distribution via
transition matrices, but the agent's **optimization problem** (Euler equation)
must use the true income distribution.
The bug arose from ordering: `ss.neutral_measure = True` and
`ss.update_income_process()` were called *before* the finite-horizon agent was
created, overwriting `ss.IncShkDstn[0]` with the neutral-measure version.
Then `FinHorizonAgent.IncShkDstn = T_cycle * [ss.IncShkDstn[0]]` passed the
wrong distribution to the solver.
**Diagnosis:** The neutral-measure distribution has the same shock values but
different probability weights (pmv). This caused the Euler equation to
produce systematically wrong consumption policies, leading to aggregate
assets ~3% below the true steady state.
**Fix:** Save the original `IncShkDstn` *before* enabling the neutral measure:
```python
orig_IncShkDstn = deepcopy(ss.IncShkDstn[0])
ss.neutral_measure = True
ss.update_income_process()
```
Then use `orig_IncShkDstn` for the finite-horizon solve (now encapsulated
in `create_finite_horizon_agent()`). After solving, the neutral measure is
enabled on `FinHorizonAgent` for transition matrix computation (this is
correct — TMs should use the neutral measure).
**Cells affected:** 36, 49 (and the `create_finite_horizon_agent` helper in cell 3)
---
## 13. Recalibrated to cstwMPC baseline (Carroll, Slacalek, Tokuoka & White 2017)
**Problem:** The previous calibration appeared to be a mix of annual-frequency
and quarterly-frequency values. In particular, `DiscFac = 0.975` is a
typical *annual* discount factor (implying ~2.5 % annual discounting), but
when used at a quarterly frequency it implies ~10 % annual discounting—far
too impatient. Similarly, `UnempPrb = 0.0` omitted unemployment risk
entirely, and no single published calibration matched the parameter set.
**Fix:** Adopted the quarterly infinite-horizon "β-Point" calibration from the
cstwMPC paper (Carroll et al., Quantitative Economics 2017, 8(3), 977–1020;
DOI: 10.3982/QE694). Key parameter changes:
| Parameter | Old | New (cstwMPC) |
|-----------|-----|---------------|
| CRRA | 2 | 1.01 (near-log; 1.01 for numerical stability) |
| Rfree | 1.04^¼ ≈ 1.010 | 1.01 / LivPrb ≈ 1.016 (mortality-adjusted) |
| DiscFac | 0.975 | 0.9867 (estimated to match K/Y and Lorenz curve) |
| PermShkStd | 0.06 | (0.01×4/11)^½ ≈ 0.0603 (same, now with formula) |
| TranShkStd | 0.2 | (0.01×4)^½ = 0.2 (same, now with formula) |
| PermShkCount | 5 | 7 |
| TranShkCount | 5 | 7 |
| UnempPrb | 0.00 | 0.07 |
| IncUnemp | 0.0 | 0.15 (15 % of permanent income) |
Unchanged: `LivPrb = 0.99375`, `PermGroFac = 1.0`, `BoroCnstArt = 0.0`.
Income-process estimates originate from Sabelhaus & Song (2010) as used in
the cstwMPC paper.
**Also added:** A markdown calibration table immediately above the parameter
dictionary (cell 4) with a citation to the cstwMPC paper.
**Cells affected:** 4, 5
---
## 14. Documentation improvements
**Problem:** Several aspects of the notebook were unclear to new readers:
- Cell 32's experiment description stated `Δr = -0.05` and `t = 10`, but the
code uses `dx = -0.01` and `shock_t = 100` — a factual contradiction.
- ~15 HARK-specific attributes and methods were used without indicating which
source files define them.
- The class name `NewKeynesianConsumerType` would lead a reader to infer this
is a HANK/NK model, when in fact no NK features are present.
- The column-stochastic matrix convention (critical for correct forward
propagation: `T @ pi`, not `pi @ T`) was mentioned only once in passing.
- Two agent objects (`example1` and `ss`) were created from the same `Dict`
without explanation.
- The Bellman equation was referenced but never written down.
**Fixes (six changes):**
1. **Fixed stale parameter values in cell 32.** Replaced hardcoded values
with symbolic `Δr` and `t*`, with a parenthetical giving the actual code
values.
2. **Added HARK API source mapping (cell 1, introduction).** Maps key methods
and attributes to their defining source files.
3. **Clarified `NewKeynesianConsumerType` naming (cell 6).** Added a note
explaining the class is used solely for its TM infrastructure.
4. **Added column-stochastic convention note (cell 58).** Explicitly states
that HARK matrices are column-stochastic.
5. **Explained `example1` / `ss` distinction (cell 33).** Notes that `ss` is
a fresh instance because the distribution-grid methods mutate internal
state.
6. **Added the one-period Bellman equation (cell 6).** States the normalized
recursive problem with all variables defined and cites Carroll (2006) for
the EGM solution method.
**Cells affected:** 1, 6, 32, 33, 58
---
## 15. Lognormal newborn permanent income (`pLogInitStd = 0.4`)
**Problem:** With the default `pLogInitStd = 0.0`, all newborn MC agents enter
at exactly `pLvl = 1.0`. Combined with the 7-point discrete approximation to
permanent shocks (`PermShkCount = 7`), this creates a "multiplicative lattice"
of preferred pLvl values near 1.0. Each period ~0.625% of agents die
(Blanchard–Yaari) and are replaced at the degenerate point, replenishing the
lattice indefinitely. The TM does not show this artifact because it propagates
a continuous probability vector via grid interpolation.
**Root cause:** Degenerate newborn pLvl distribution + discrete permanent
shock approximation → lattice of products of shock values that never fully
smooths out in the MC histogram.
**Fix (two parts):**
1. **MC side (cell 5):** Set `pLogInitStd = 0.4` and `pLogInitMean = −0.08`
(Jensen correction so `E[pLvl] = 1.0`), with `pLvlInitCount = 25`.
The value 0.4 follows the cstwMPC life-cycle model
([Carroll, Slacalek, Tokuoka, and White (2017)](https://doi.org/10.3982/QE694)),
which calibrates it to match the total variance of income among young
households in the SCF 2004 data.
2. **TM side (cell 13):** HARK's `calc_transition_matrix` hardcodes
`NewBornDist` at `pLvl = 1.0`. We post-process the transition matrix
via the `correct_newborn_dist()` helper: subtract
`(1 − LivPrb) × old_NBD` and add `(1 − LivPrb) × new_NBD`
(lognormal-distributed across the p-grid). This preserves column sums
and ensures MC and TM solve the same economic model.
For the Harmenberg (1D) case, no correction is needed: the neutral measure
integrates out the permanent income dimension, and the Jensen correction
ensures `E[pLvl] = 1.0` for newborns.
**Cells affected:** 3 (helper), 5, 13
---
## 16. Distribution plots: probability density, not probability mass
**Problem:** All four cross-sectional distribution plots (normalized market
resources, permanent income, market resources in levels, liquid assets) plotted
**probability mass per bin** — `histogram_counts / total_counts` — rather than
**probability density**. With the double-exponential grids (`timestonest=2`),
bin widths vary by an order of magnitude. Narrower bins near the lower bound
(or near pLvl = 1.0 for the p-grid) naturally contain less mass, creating
spurious dips in the plotted distribution even when the true density is smooth.
This was the root cause of the "hole near pLvl = 1.0" that appeared in the
permanent income distribution.
**Fix:** For TM, divide probability mass at each grid point by the midpoint
bin width to obtain proper probability density. For MC, use
`np.histogram(..., density=True)` with 200 uniform bins. Change y-axis label
from "Probability" to "Probability Density."
**Additional issue in liquid assets plot:** `aPol_Grid` has ~22 entries at
exactly 0.0 (the borrowing constraint — for low-m agents, optimal savings
is a = 0). This creates zero-width bins, producing division-by-zero when
converting to density. Fix: skip to the first positive `aPol_Grid` entry
and report the mass point at the constraint separately as a printed
diagnostic.
**Cells affected:** 23 (mNrm), 25 (pLvl), 29 (mLvl), 31 (aLvl)
---
## 17. Code quality and structural improvements
**Problem:** The notebook accumulated various code quality issues during
iterative development: scattered imports, duplicated code blocks, Python
loops where vectorized numpy operations apply, variable shadowing bugs,
inconsistent plot formatting, and magic numbers.
**Fixes:**
1. **Consolidated imports.** All imports (`init_newkeynesian`,
`jump_to_grid_2D`, `LognormalDist`) moved to the single imports cell
(cell 3), eliminating scattered `from ... import` statements in cells 5
and 13.
2. **Extracted helper functions.** Two reusable functions defined in cell 3:
- `correct_newborn_dist(agent, param_dict)` — patches the TM newborn
distribution for lognormal pLvl (was ~20 inline lines in cell 13).
- `create_finite_horizon_agent(ss, Dict, T_cycle, shock_t, dx,
orig_IncShkDstn)` — creates, configures, and returns a finite-horizon
agent (was ~30 duplicated lines in cells 49 and 68).
- `jump_to_grid_fast(m_vals, probs, dist_mGrid)` — moved from cell 28 to
the utilities cell for earlier availability.
3. **Vectorized loops.**
- `gridc`/`grida` computation (cell 14): replaced `for j in range(len(p))`
loop with `np.outer(c_2d, p_grid_2d)`.
- `mLvl_vals`/`aLvl_vals` (cell 27): replaced double-nested append loops
with `np.outer(...).ravel()`.
4. **Fixed variable shadowing.**
- Cell 27's inner loop used `p` as the loop variable, overwriting
`p = example1.dist_pGrid` from cell 13. Eliminated by vectorization.
- Cell 37 assigned `c = ss.cPol_Grid` and `a = ss.aPol_Grid`, overwriting
`c` and `asset` from cell 13. Renamed to `c_harm` / `a_harm`; 2D-grid
variables renamed to `c_2d` / `asset_2d`.
5. **Named constants.** Introduced `BURNIN = 500`, `N_MC_BINS = 200`,
`N_P_DISC = 50`, `MAX_P_FAC = 10.0` in cell 3, replacing bare literals
throughout.
6. **Unified variable names.** Removed `num_pts` (duplicate of `n_m_grid`).
Defined `n_m_grid`, `n_p_grid`, and `n_agents` once in cell 13 after the
TM computation, rather than inside plotting cells. Renamed `N` to
`N_fin` in cell 55 for clarity.
7. **Standardized plot formatting.** Set `plt.rcParams` defaults for figure
size, label size, title size, legend size, and line width in cell 3.
Removed per-cell `linewidth=`, `fontsize=`, and `prop={"size": ...}`
overrides. Added `plt.tight_layout()` and explicit `figsize` to all
plotting cells.
8. **Standardized density plot ranges.** All four density plots now use
the same percentile-based range (0.5–99.5) for both MC and TM.
9. **Removed "Known HARK Issues and Workarounds" section.** The three
HARK workarounds (t_age suppression, cFunc_terminal_ vs
solution_terminal, neutral-measure IncShkDstn leakage) are documented
with brief inline comments at their respective workaround sites (cells
36, 49/helper, 55). A standalone section was unnecessary.
10. **Unhid cell 5.** Removed `"source_hidden": true` from the calibration
dictionary cell metadata so the parameters are visible by default.
**Cells affected:** 3, 5, 13, 14, 18, 19, 23, 25, 27, 28, 29, 31, 37, 39,
42, 49, 51, 55, 61, 63, 68, 71
Made-with: Cursor
New module ConsAggIndMarkovModel.py provides AggIndMrkvConsumerType, a MarkovConsumerType subclass that decomposes combined Markov states into aggregate (macro) and idiosyncratic (micro) components with a two-step draw: macro state from the economy, then micro states conditional on the macro state. KrusellSmithType now inherits from AggIndMrkvConsumerType instead of bare AgentType, using the hierarchical decomposition for its 2 macro (bad/good) x 2 micro (unemployed/employed) state structure. The internal shock variable is renamed from "Mrkv" to "MrkvAgg" to distinguish the aggregate state from the combined index. This class is also used by the HAFiscal project's AggFiscalType for aggregate-demand models with fiscal policy shocks. All existing KrusellSmith tests pass with identical numerical results. The implementation was guided by a design prompt and verified by Claude Opus 4.6 with extensive automated testing across both the Krusell-Smith and HAFiscal codebases. Made-with: Cursor
Set default_["model"] = "ConsIndShock.yaml" so initialize_sym() can build the YAML-driven AgentSimulator for NewKeynesianConsumerType (same dynamics as IndShockConsumerType). Needed for Transition_Matrix_Example new-API cells. Made-with: Cursor
- MarkovConsumerType / ConsAggShockModel: transition-matrix and related infrastructure (incl. hierarchical / KS-oriented paths as developed) - simulator.py: AgentSimulator improvements (e.g. explicit grid support, projection onto legacy-style grids where needed) - utilities.py: gen_tran_matrix_1D_markov for block-structured Markov TMs - distributions/base.py: small compatibility fix - tests: extend test_ConsMarkovModel.py - docs/reference/index.rst: index updates - Transition_Matrix_Example.ipynb: SSJ / TM vs MC example updates - Add debug notebook variants for TME investigation - sims-about: notebook and LESSONS-LEARNED / framework updates - .gitignore: ignore sims-about/** for local scratch; add .ragignore Made-with: Cursor
Drop 01-markov-tm-prototype and 05-tm-consolidation from the repo; copies retained outside the tree by author request. Made-with: Cursor
…nfra-and-examples
Consolidate seven sims-about/ notebooks into three well-organized examples in examples/MonteCarlovsTransitionMatrix/: - PE_MarkovConsumerType.ipynb: partial-equilibrium MC vs TM for 4-state unemployment, 5-state PermGroFac (2D grid), and Harmenberg neutral measure - GE_KrusellSmith.ipynb: general-equilibrium TM with aggregate shocks and TM forward propagation inside the KS loop - Validation_and_SSJ.ipynb: TM method validation and sequence-space Jacobians via Fake News Algorithm Also adds: - make_history_tm() on CobbDouglasMarkovEconomy for deterministic TM-based forward propagation through aggregate shock sequences - docs/guides/transition_matrix_methods.md guide - Registration in docs/example_notebooks/Include_list.txt Made-with: Cursor
Add missing "name": "stdout" to stream output and missing "metadata": {}
to display_data outputs. These caused the Sphinx docs build (Render CI
job) to fail with NotebookValidationError.
Made-with: Cursor
(Cherry-picked from fix/markov-newborn-permshk-missing-growth) Made-with: Cursor
Collaborator
Author
|
|
Contributor
There was a problem hiding this comment.
Pull request overview
This PR adds end-to-end infrastructure and documentation for comparing Monte Carlo (MC) simulation with transition-matrix (TM) methods in HARK, extending TM support to Markov and aggregate-shock (KS) settings and registering new example notebooks/docs.
Changes:
- Add TM construction/ergodic/steady-state/Jacobian methods to
MarkovConsumerType, plus a Numba TM builder utility. - Add documentation + example notebook registrations for TM methods and MC-vs-TM comparisons.
- Introduce hierarchical macro+micro Markov infrastructure for Krusell–Smith–style models (new
ConsAggIndMarkovModel, plus related updates).
Reviewed changes
Copilot reviewed 36 out of 44 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
HARK/ConsumptionSaving/ConsMarkovModel.py |
Adds TM grid/TM/ergodic/steady-state/Jacobian methods; fixes newborn PermShk growth logic. |
HARK/utilities.py |
Adds gen_tran_matrix_1D_markov (Numba) for Markov TM construction. |
tests/ConsumptionSaving/test_ConsMarkovModel.py |
Adds regression test for newborn PermShk growth + TM validation tests (column sums, ergodic fractions, etc.). |
HARK/ConsumptionSaving/ConsAggShockModel.py |
Refactors KS agent to hierarchical Markov base; adds make_history_tm; introduces HM reference implementations. |
HARK/ConsumptionSaving/ConsAggIndMarkovModel.py |
New hierarchical Markov consumer type + helper functions for macro/micro decomposition. |
examples/ConsAggShockModel/test_ks_hierarchical_markov.py |
Adds a comparison “test” script for original vs HM KS implementations. |
docs/guides/transition_matrix_methods.md |
New guide explaining TM methods and MC→TM conversion workflow. |
docs/guides/index.rst |
Registers the new TM methods guide. |
docs/example_notebooks/Include_list.txt |
Registers new MC-vs-TM example notebooks for docs builds. |
docs/reference/index.rst + docs/reference/ConsumptionSaving/ConsAggIndMarkovModel.rst |
Adds API reference page for the new module. |
docs/CHANGELOG.md |
Notes new module + KS Markov shock variable rename (Mrkv→MrkvAgg). |
HARK/simulator.py |
Adds support for explicit grids + multi-exponential grids, and improves grid casting for arbitrary grids. |
HARK/distributions/base.py |
Clarifies Markov transition-matrix convention (row-stochastic). |
HARK/ConsumptionSaving/ConsNewKeynesianModel.py |
Adds model metadata entry to defaults. |
.ragignore |
Adds a repo-level indexing ignore config. |
.gitignore |
Ignores sims-about/**. |
sims-about/** |
Adds development/provenance notebooks and reference docs (appears to be scratch/prototyping material). |
Comments suppressed due to low confidence (1)
.ragignore:213
.ragignorerepeats top-level YAML keys (ignore_patternsandcareful_processing) later in the file. In YAML, duplicate keys typically mean the earlier values are overwritten, so many ignore patterns / settings above will be silently dropped. Merge these sections into singleignore_patternsandcareful_processingmappings to avoid losing rules.
# Additional patterns based on HARK repository analysis
ignore_patterns:
# Test files (major noise reduction)
- "*/tests/"
- "test_*.py"
- "*_test.py"
# GitHub/CI files
- ".github/"
- ".pre-commit-config.yaml"
# Large calibration data
- "*/Calibration/*/life_tables/"
- "*/Calibration/*/SCF/"
- "*/Calibration/*/cpi/"
- "large_data/"
- "calibration_data/"
# Documentation build
- "Makefile"
- "make.bat"
- "conf.py"
# Enhanced notebook processing
careful_processing:
large_notebooks:
patterns: ["*.ipynb"]
max_size_mb: 2
strip_outputs: true
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+1057
to
+1060
| for j in range(self.MrkvArray[0].shape[0]): | ||
| these_nb = np.logical_and(newborn, j == MrkvNow) | ||
| if np.any(these_nb): | ||
| PermShkNow[these_nb] = self.PermGroFac[0][j] |
There was a problem hiding this comment.
Newborn PermShk is set using self.PermGroFac[0][j], which hard-codes period 0. In lifecycle / finite-horizon configurations (T_cycle>1), newborns created mid-simulation should use the current period’s growth factor (consistent with the t-1 indexing used above), e.g. indexing by self.t_cycle-1 (and the appropriate MrkvArray[t-1] shape) rather than [0].
Suggested change
| for j in range(self.MrkvArray[0].shape[0]): | |
| these_nb = np.logical_and(newborn, j == MrkvNow) | |
| if np.any(these_nb): | |
| PermShkNow[these_nb] = self.PermGroFac[0][j] | |
| # For newborns, apply only deterministic growth (no idiosyncratic ψ), | |
| # using the same t-1 indexing convention as above. | |
| for t in range(self.T_cycle): | |
| these_t_nb = np.logical_and(newborn, self.t_cycle == t) | |
| if not np.any(these_t_nb): | |
| continue | |
| for j in range(self.MrkvArray[t - 1].shape[0]): | |
| these_nb = np.logical_and(these_t_nb, j == MrkvNow) | |
| if np.any(these_nb): | |
| PermShkNow[these_nb] = self.PermGroFac[t - 1][j] |
Comment on lines
+1286
to
+1289
| shk_prbs = shk_k[0].pmv | ||
| perm_shks = shk_k[0].atoms[0] | ||
| tran_shks = shk_k[0].atoms[1] | ||
|
|
There was a problem hiding this comment.
calc_transition_matrix currently takes Markov-state shock distributions (shk_k) but then unconditionally uses shk_k[0] for shk_prbs/perm_shks/tran_shks. That means the TM is built assuming the same income shock distribution in every Markov state, which is incorrect when IncShkDstn differs by state (a common Markov setup). Either extend gen_tran_matrix_1D_markov/builder logic to accept per-state shock distributions, or add an explicit check/assertion that all shk_k[j] are identical and document the limitation.
Suggested change
| shk_prbs = shk_k[0].pmv | |
| perm_shks = shk_k[0].atoms[0] | |
| tran_shks = shk_k[0].atoms[1] | |
| # NOTE: gen_tran_matrix_1D_markov currently accepts a single income shock | |
| # distribution, so we must ensure that all Markov states share the same | |
| # shock distribution. If they differ, this builder does not yet support | |
| # truly state-varying income shocks and will raise an explicit error | |
| # instead of silently using shk_k[0] for all states. | |
| first_shk = shk_k[0] | |
| for _j_state, _shk in enumerate(shk_k[1:], start=1): | |
| if not ( | |
| np.array_equal(_shk.pmv, first_shk.pmv) | |
| and np.array_equal(_shk.atoms[0], first_shk.atoms[0]) | |
| and np.array_equal(_shk.atoms[1], first_shk.atoms[1]) | |
| ): | |
| raise NotImplementedError( | |
| "calc_transition_matrix/_build_one_tm currently assume an identical " | |
| "income shock distribution across all Markov states. State-varying " | |
| "IncShkDstn is not yet supported in this transition-matrix builder." | |
| ) | |
| shk_prbs = first_shk.pmv | |
| perm_shks = first_shk.atoms[0] | |
| tran_shks = first_shk.atoms[1] |
| if self.neutral_measure: | ||
| self.dist_pGrid = np.array([1]) | ||
| else: | ||
| self.dist_pGrid = np.array([1]) |
There was a problem hiding this comment.
When neutral_measure is False, define_distribution_grid still sets dist_pGrid = np.array([1]), silently behaving like the neutral-measure 1D case. This makes it easy for callers to think they’re getting a full (m,p,j) grid when they are not. Consider raising NotImplementedError (or implementing the full p-grid path) when neutral_measure is False.
Suggested change
| self.dist_pGrid = np.array([1]) | |
| raise NotImplementedError( | |
| "define_distribution_grid currently supports only the neutral " | |
| "measure (1D m-grid). Non-neutral measure with a full p-grid " | |
| "is not yet implemented." | |
| ) |
Comment on lines
+1621
to
+1623
| dstn_z = (z_tm if i == 0 else tranmat_ss) @ dstn_z | ||
| C_t_z[i] = np.dot(c_ss_flat, dstn_z) | ||
| A_t_z[i] = np.dot(a_ss_flat, dstn_z) |
There was a problem hiding this comment.
In the zeroth-column propagation loop, the code transitions the distribution before computing aggregates (dstn_z = TM @ dstn_z then A = a' @ dstn_z). This is an off-by-one convention relative to the usual “compute aggregates from current distribution, then transition” ordering and can systematically shift/scale the IRF/Jacobian validation. Consider swapping the order so A_t_z[i]/C_t_z[i] are computed from the pre-transition distribution at date i.
Suggested change
| dstn_z = (z_tm if i == 0 else tranmat_ss) @ dstn_z | |
| C_t_z[i] = np.dot(c_ss_flat, dstn_z) | |
| A_t_z[i] = np.dot(a_ss_flat, dstn_z) | |
| # Compute aggregates from the current (pre-transition) distribution | |
| C_t_z[i] = np.dot(c_ss_flat, dstn_z) | |
| A_t_z[i] = np.dot(a_ss_flat, dstn_z) | |
| # Then transition the distribution to the next period | |
| dstn_z = (z_tm if i == 0 else tranmat_ss) @ dstn_z |
Comment on lines
+2630
to
+2645
| mc_M = np.array(self.history["MaggNow"]) | ||
| mc_A = np.array(self.history["AaggNow"]) | ||
|
|
||
| TranMatrices = [] | ||
| aPols = [] | ||
| cPols = [] | ||
| for j in range(StateCount): | ||
| mask = MrkvHist == j | ||
| M_j = mc_M[mask].mean() if mask.any() else self.kSS | ||
| A_j = mc_A[mask].mean() if mask.any() else self.kSS | ||
| K_j = A_j | ||
| R_j = self.Rfunc(K_j) | ||
| W_j = self.wFunc(K_j) | ||
|
|
||
| c_j = consumer.solution[0].cFunc[j](dist_mGrid, M_j * np.ones(M)) | ||
| a_j = np.maximum(dist_mGrid - c_j, 0.0) |
There was a problem hiding this comment.
make_history_tm builds state-specific TMs using M_j/A_j computed as means from the existing MC history (self.history["MaggNow"], self.history["AaggNow"]). That makes the TM path depend on a prior MC run and prevents this method from being a standalone deterministic forward propagator given an aggregate shock sequence and the converged AFunc/policies (as described in the PR). Consider computing prices/policies each period from the current TM-implied distribution (or at least from the current MaggNow in the Markov history) and avoid reading MC history inside this method. Also, LivPrb is coerced to a scalar via LivPrb = LivPrb[0], which will be wrong if survival is state-dependent.
Comment on lines
+1
to
+15
| """ | ||
| Comparison test: original KrusellSmithType vs. AggIndMarkovConsumerType-based | ||
| KrusellSmithTypeHM. Both models use identical parameters and the same | ||
| aggregate Markov shock history. We verify that the converged aggregate | ||
| saving rules match within tolerance. | ||
| """ | ||
|
|
||
| import time | ||
|
|
||
| from HARK.ConsumptionSaving.ConsAggShockModel import ( | ||
| KrusellSmithType, | ||
| KrusellSmithEconomy, | ||
| KrusellSmithTypeHM, | ||
| KrusellSmithEconomyHM, | ||
| ) |
There was a problem hiding this comment.
This file is named test_*.py and contains top-level execution (prints + two full Krusell–Smith solves). With the repo’s current setup (no pytest collection restrictions), pytest will collect and run it, making the test suite extremely slow and non-deterministic. Please move this to a non-test location/name (e.g. examples/.../compare_ks_hierarchical_markov.py), guard it with if __name__ == '__main__':, or mark it as skipped for pytest (e.g. pytest.skip(..., allow_module_level=True)).
The transition-matrix infrastructure (define_distribution_grid, calc_transition_matrix, calc_ergodic_dist, make_history_tm, compute_pe_steady_state, calc_jacobian) added earlier on this branch is removed from MarkovConsumerType pending a cleaner integration via the AgentSimulator pipeline. Corresponding tests and unused imports (scipy.sparse, gen_tran_matrix_1D_markov, jump_to_grid_1D, make_grid_exp_mult) are also removed. Made-with: Cursor
For infinite-horizon models with finite T_age, all agents previously started at age 0 and hit T_age simultaneously, creating a periodic "cohort echo" that biases aggregates (e.g. ~1.25% in E[pLvl] for T_age=200, LivPrb=0.99375). New _initialize_ergodic_ages() draws initial ages from the truncated geometric distribution and scales pLvl by PermGroFac^age. Controlled by init_ages_ergodic attribute (default True); skipped for lifecycle models. Made-with: Cursor
When read_shocks is active, the newborn initialization loop warned about any state variable missing from newborn_init_history. Aggregate scalars like PlvlAgg are not per-agent arrays and are set elsewhere, so the warning was noise. Now only warns for idiosyncratic states (arrays of length AgentCount). Made-with: Cursor
Automatic reformat by nbformat: moves cell "id" field to canonical position per nbformat >= 5.7. No code or content changes. Made-with: Cursor
…eynesian Start from the default NewKeynesianConsumerType parameters and override only the values that define this KS calibration: preferences, income process (no unemployment), equilibrium prices, and wider grids. Made-with: Cursor
… names Automatic nbformat fixes: move cell "id" to canonical position, add required "name" field to stream outputs, clear stale execution counts. No code or content changes. Made-with: Cursor
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
examples/MonteCarlovsTransitionMatrix/that compare Monte Carlo (MC) and Transition Matrix (TM) simulation methods forMarkovConsumerType:PE_MarkovConsumerType.ipynb: Partial-equilibrium comparisons — 4-state unemployment (1D grid), 5-state PermGroFac (2D grid), and Harmenberg neutral measureGE_KrusellSmith.ipynb: General-equilibrium TM with aggregate shocks and TM forward propagation inside the Krusell-Smith loopValidation_and_SSJ.ipynb: Validates production TM methods and demonstrates sequence-space Jacobians via the Fake News Algorithmmake_history_tm()method added toCobbDouglasMarkovEconomyfor deterministic TM-based forward propagation through aggregate Markov shock sequences, enabling direct MC-vs-TM comparison of aggregate pathsdocs/guides/transition_matrix_methods.mdcovering the mathematical framework and learning pathAgentSimulator,initialize_sym()/symulate()/hystoryAPI, and grid-based distribution utilitiesBackground
These notebooks were developed iteratively in a temporary
sims-about/directory while building and debugging TM support. This PR consolidates the seven prototype notebooks into three well-organized examples suitable for permanent inclusion in the repository.Test plan
jupyter nbconvert --execute)make_history_tm()produces TM aggregate paths that track MC paths in KS economypytest tests/)Made with Cursor