HabitPortfolioConsumerType: habit formation with portfolio choice#1748
Draft
HabitPortfolioConsumerType: habit formation with portfolio choice#1748
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
CAVEAT EMPTOR: I am skeptical of the model's results late in life, and particularly skeptical of the response to the "tragic" wealth shock. The model shows a -19% consumption drop and +11% jump in risky share at the point of tragedy, and the gap from the no-shock path *widens* through retirement rather than recovering. The narrative that agents are "gambling for resurrection" — taking on more risk to rebuild wealth relative to their habit stock — is probably plausible-sounding BS that masks a real bug in the solver or the simulation dynamics. This needs careful review before anyone should trust the late-life results. What was built: A new AgentType combining the 2D habit formation model (ConsHabitModel) with risky portfolio choice (ConsPortfolioModel). The agent has CRRA preferences over consumption relative to a habit stock u(c/h^α), and allocates savings between risk-free and risky assets. The state space is (m_t, h_t) with two controls (c_t, s_t). Files: - ConsHabitPortfolioModel.py: solver + AgentType class (603 lines) - ConsHabitPortfolio.yaml: YAML simulation dynamics - HabitPortfolioConsumerType.ipynb: math derivations + demo notebook - HabitPortfolioConsumerType.md: AI-readable mathematical description - test_HabitPortfolio_convergence.ipynb: validates convergence to standard PortfolioConsumerType as habit weight α → 0 - test_ConsHabitPortfolioModel.py: basic regression tests - __init__.py: registers the new type in ConsumptionSaving exports Development process: This model was built interactively with Claude Opus 4.6 in Cursor (conversation 6d75051f-53ac-44e9-9d21-9d7943a7ab19). The process was: 1. Mathematical derivations: started from the habit model's Bellman equation and FOCs, added portfolio choice following the two-stage optimization approach from ConsPortfolioModel. Key insight: habits make risk tolerance state-dependent. 2. Implementation: 3-stage solver (share search over risky returns, habit EGM with FOC-inverter, marginal values via income integration). Followed HARK conventions (dict-based solutions, Curvilinear2DInterp, BilinearInterp for marginal values, YAML-based simulation). 3. Iterative debugging: fixed ShareLimit list/scalar handling, terminal period marginal values, -inf singularities in direct evaluation of next-period marginal values (replaced with dvdH_cont interpolant), and empty simulation histories (T_sim initialization). 4. Validation: convergence test shows clean α → 0 convergence for moderate α values, with expected plateau at very small α due to 2D vs 1D interpolation differences. Known issues: - Late-life tragic shock results are suspect (see caveat above) - Convergence test shows plateau at α < 0.01 (numerical, not economic) - Default CRRA=2 gives 100% risky share for most of life; notebook uses CRRA=5 for more interesting portfolio dynamics Made-with: Cursor
Contributor
|
I'm not vouching for the validity of the code sight unseen, but the late-life behavior that you describe doesn't strike me as facially wrong. I'm not sure about the AI's interpretation of it, but here's how I would think about it: After a big negative wealth shock, the risky asset share will jump up whether or not there is habit formation. As long as the consumers were operating on the downward sloping (w.r.t. The "wiggly" risky share function at lower |
…tation - Calibrate DiscFac (0.849) for portfolio model to match habit-only retirement wealth, enabling apples-to-apples comparison - Add calibrate_DiscFac.ipynb with bisection search - Use SSA life tables for mortality and CGM (2005) income process with averaged HS/College replacement rate at retirement - Fix ShareFunc solver: ClippedBilinearInterp prevents extrapolation errors; monotonize Share_opt to eliminate FOC noise; use BilinearInterp on regular (m,h) grid instead of Curvilinear2DInterp - Use stop_dead=False and np.mean consistently across all simulations - Add named figures: [lifecycle-asset-accumulation], [lifecycle-consumption-and-portfolio-allocation], [habit-only-vs-habit-portfolio-consumption], [tragic-wealth-shock-at-age-67], [consumption-dynamics-after-tragic-shock] - Add explanatory text on DiscFac calibration, equity premium effects on consumption levels, and human wealth effect on portfolio rebalancing Made-with: Cursor
llorracc
changed the title
…gic shock Add markdown on Carroll et al. habit weight and new code cells comparing wealth and consumption after the tragic shock for HabitWgt=0.7 and 0.9 vs baseline. Refresh outputs for earlier cells from re-execution. Made-with: Cursor
(Cherry-picked from fix/markov-newborn-permshk-missing-growth) Made-with: Cursor
Collaborator
Author
|
|
These should not be in the PR, and they don't all pertain to this model anyway.
These seem to be from other PRs.
Claude's code is not wildly wrong, but it did some funny things. This PR cleans up the representation of the risky share function a bit, but there's more work to be done. The bit of code that enforces that risky share is non-increasing in wNrm produces nice results, but I worry that it's just brushing problems under the carpet when they should be addressed more directly.
Contributor
|
I'm working on reviewing and testing what Claude produced, and it's pretty good. The next bit of improvement would be easier / cleaner if #1753 were merged into main (and then into this branch). |
Contributor
|
I've done a bunch more work to clean up the AI's solver. There were a few issues, some of which were due to it copying small mistakes from the base habit-formation solver (now fixed). But other problems were due to badly chosen grids, or an incorrect assumption about what the solution "should" look like. I'm still not 100% sure it's correct, but it seems like the risky share function might not be monotonic after all. As early (going backward in time) as period T-2 it develops a smooth wiggle between roughly m=1 and m=2. That wiggle expands and distorts as you go further back in time. I'll look at the code some more to make sure there isn't some other problem, but I think we need to really work through the theory on this and see if a non-monotone risky share should be possible. I'll also do a version of the code fully on my own and see if I get the same solution. The other bit that makes me concerned that there's still a serious bug in the code or method is that things go wildly wrong if you set |
Fixed small bugs, implemented re-interpolation, adjusted grids, used pseudo-inverse trick, removed forced monotonicity.
Contributor
|
Also: For the demo notebook, I feel like the more relevant comparison is with and without habits. I.e. compare outcomes from |
Also adjusted test target so it will pass. Change is due to adjustments to solution method and parameters.
This notebook has a few parameter experiments with the habit-portfolio model. There's no commentary yet, and I only play with HabitWgt, not HabitRte.
Mostly just added comments and changed some coding style and variable names.
ConsHabitModel.py now has an alternative version of HabitPortfolioConsumerType, which is fully untested. This one uses a modular solver with standalone functions for the portfolio allocation and consumption-saving problems. The latter is just the solver for the basic habit formation model. If it works, then it can be used for other two-state models with portfolio choice.
It looks like the solver is breaking with even moderately faster habit rate. This might be due to an extrapolation issue when consumption is large (near T), which will be more of a problem when lambda is higher because the habit stock updates faster.
It worked straight away, other than these two dumb typos.
This actually speeds it up slightly by avoiding wasted work.
Contributor
There was a problem hiding this comment.
Pull request overview
Adds a new lifecycle Habit+Portfolio consumer type to HARK, including YAML simulator support, solver implementation, tests, and accompanying calibration/validation materials.
Changes:
- Introduces
HabitPortfolioConsumerType(habit formation + risky/risk-free portfolio choice) with a new solver and YAML model. - Updates habit-model grid defaults (and regression test value) and expands
ConsumptionSavingexports to include the new type. - Adds notebooks and a mathematical write-up for validation/convergence and discount-factor calibration.
Reviewed changes
Copilot reviewed 10 out of 13 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
HARK/ConsumptionSaving/ConsHabitPortfolioModel.py |
New solver + AgentType for habit formation with portfolio choice, plus defaults and helper utilities. |
HARK/models/ConsHabitPortfolio.yaml |
YAML simulator dynamics for the new model with ShareFunc(m,h). |
HARK/models/ConsHabitPortfolioAlt.yaml |
Alternative YAML specification where ShareFunc is defined over (w,H). |
HARK/ConsumptionSaving/__init__.py |
Exports HabitPortfolioConsumerType at package level. |
HARK/ConsumptionSaving/ConsHabitModel.py |
Adds portfolio-related helper functions/types and tweaks habit-grid defaults (plus additional solver/class). |
tests/ConsumptionSaving/test_ConsHabitPortfolioModel.py |
Adds regression + smoke tests for solving and simulating the new habit-portfolio type. |
tests/ConsumptionSaving/test_ConsHabitModel.py |
Updates regression target for HabitConsumerType after grid/default changes. |
examples/ConsHabitModel/test_HabitPortfolio_convergence.ipynb |
Validation notebook comparing convergence to standard portfolio model as habit weight → 0. |
examples/ConsHabitModel/calibrate_DiscFac.ipynb |
Standalone bisection calibration of DiscFac to match retirement wealth across models. |
examples/ConsHabitModel/HabitPortfolioConsumerType.md |
Mathematical description of the model and algorithm. |
.gitignore |
Ignores sims-about/**. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| mNrmMin = wNrmMin | ||
| PermShkVals = IncShkDstn.atoms[0, :] | ||
| TranShkVals = IncShkDstn.atoms[1, :] | ||
| kNrmMin_cand = (mNrmMin - TranShkVals) / Rfree * (PermShkVals * PermGroFac) |
There was a problem hiding this comment.
Natural borrowing constraint computation divides by Rfree, but this model’s state transition uses m = k/G + θ (no Rfree factor; portfolio return is already embedded in k via the twist). This makes kNrmMin inconsistent with the simulated dynamics and can yield an incorrect lower bound and interpolation domain. Compute kNrmMin_cand without /Rfree (k ≥ (mMin-θ)*G).
Suggested change
| kNrmMin_cand = (mNrmMin - TranShkVals) / Rfree * (PermShkVals * PermGroFac) | |
| kNrmMin_cand = (mNrmMin - TranShkVals) * (PermShkVals * PermGroFac) |
Comment on lines
+157
to
+165
| 4. Build policy functions on the endogenous $(\hat{m}, \hat{h})$ grid: | ||
| - `cFunc(m, h)` as a `Curvilinear2DInterp` of $\hat{c}$ on $(\hat{m}, \hat{h})$. | ||
| - `ShareFunc(m, h)` as a `Curvilinear2DInterp` of $s^*$ on the same grid. | ||
|
|
||
| 5. Add borrowing constraint: `cFunc = LowerEnvelope2D(cFuncUnc, cFuncCnst)`. | ||
|
|
||
| **HARK tools:** | ||
| - `HabitFormationInverter` (from `ConsHabitModel.py`): Pre-computed lookup table that inverts the nonlinear FOC mapping $(H, \chi) \to (c, h)$. | ||
| - `Curvilinear2DInterp`: Interpolation on a warped (endogenous) 2D grid. |
There was a problem hiding this comment.
Stage 2 description says both cFunc and ShareFunc are built as Curvilinear2DInterp on the endogenous (m̂,ĥ) grid, but the implementation re-interpolates cFunc onto a regular (m,h) grid via BilinearInterp and builds ShareFunc via LinearInterpOnInterp1D + clipping. Update the algorithm description to match the implemented interpolation strategy (and note the reason for using a regular grid for ShareFunc).
Suggested change
| 4. Build policy functions on the endogenous $(\hat{m}, \hat{h})$ grid: | |
| - `cFunc(m, h)` as a `Curvilinear2DInterp` of $\hat{c}$ on $(\hat{m}, \hat{h})$. | |
| - `ShareFunc(m, h)` as a `Curvilinear2DInterp` of $s^*$ on the same grid. | |
| 5. Add borrowing constraint: `cFunc = LowerEnvelope2D(cFuncUnc, cFuncCnst)`. | |
| **HARK tools:** | |
| - `HabitFormationInverter` (from `ConsHabitModel.py`): Pre-computed lookup table that inverts the nonlinear FOC mapping $(H, \chi) \to (c, h)$. | |
| - `Curvilinear2DInterp`: Interpolation on a warped (endogenous) 2D grid. | |
| 4. Build policy functions using the endogenous grid and then re-interpolating onto a regular $(m, h)$ grid: | |
| - First construct $\hat{c}$ on the endogenous $(\hat{m}, \hat{h})$ grid, then re-interpolate it onto a regular $(m, h)$ grid via `BilinearInterp` to obtain the final `cFunc(m, h)`. | |
| - Construct `ShareFunc(m, h)` on the same regular $(m, h)$ grid using `LinearInterpOnInterp1D` over portfolio shares and then clip the result to the admissible share bounds. Using a regular grid for $(m, h)$ ensures that the share search and interpolation are numerically stable and compatible with the portfolio-optimization routine. | |
| 5. Add borrowing constraint: `cFunc = LowerEnvelope2D(cFuncUnc, cFuncCnst)`. | |
| **HARK tools:** | |
| - `HabitFormationInverter` (from `ConsHabitModel.py`): Pre-computed lookup table that inverts the nonlinear FOC mapping $(H, \chi) \to (c, h)$. | |
| - `BilinearInterp`: Interpolation of $\hat{c}$ from the endogenous $(\hat{m}, \hat{h})$ grid onto a regular $(m, h)$ grid. | |
| - `LinearInterpOnInterp1D`: Builds `ShareFunc` on the regular $(m, h)$ grid by interpolating over portfolio shares, with clipping to the feasible range. |
Comment on lines
+263
to
+267
| "cFunc": Curvilinear2DInterp, # c(m, h) | ||
| "ShareFunc": Curvilinear2DInterp, # s(m, h) | ||
| "dvdkFunc": MargValueFuncCRRA, # ∂v/∂k(k, hPre) | ||
| "dvdhFunc": ValueFuncCRRA, # ∂v/∂hPre(k, hPre) | ||
| "kNrmMin": float, # minimum capital |
There was a problem hiding this comment.
The documented solution object types list cFunc/ShareFunc as Curvilinear2DInterp, but the solver currently returns cFunc as a BilinearInterp (after re-interpolation) and ShareFunc as a Clipped2DInterp wrapping a LinearInterpOnInterp1D. Please update these type annotations so the docs reflect what users actually receive in solution[t].
Suggested change
| "cFunc": Curvilinear2DInterp, # c(m, h) | |
| "ShareFunc": Curvilinear2DInterp, # s(m, h) | |
| "dvdkFunc": MargValueFuncCRRA, # ∂v/∂k(k, hPre) | |
| "dvdhFunc": ValueFuncCRRA, # ∂v/∂hPre(k, hPre) | |
| "kNrmMin": float, # minimum capital | |
| "cFunc": BilinearInterp, # c(m, h) | |
| "ShareFunc": Clipped2DInterp, # s(m, h), wrapping a LinearInterpOnInterp1D | |
| "dvdkFunc": MargValueFuncCRRA, # ∂v/∂k(k, hPre) | |
| "dvdhFunc": ValueFuncCRRA, # ∂v/∂hPre(k, hPre) | |
| "kNrmMin": float, # minimum capital |
Comment on lines
+166
to
+169
| "# Bisection search\n", | ||
| "lo, hi = 0.80, 0.96\n", | ||
| "tol = 0.005 # match within 0.5% of target wealth\n", | ||
| "max_iter = 20\n", |
There was a problem hiding this comment.
The bisection tolerance is defined as tol = 0.005 (0.5%), but the convergence check later hard-codes abs(gap_pct) < 0.5 and never uses tol. Use tol in the stopping rule (or rename it to a pct tolerance) so the parameter actually controls convergence and the comment stays accurate.
The previously skeletal "dynamic experiments" notebook has now been improved. The simulations use grid methods rather than Monte Carlo, and beta & rho have been calibrated to match long run averages across types. There are also now section headings and a description of what's going on. Discussion of results is still missing, but I will add that tonight or tomorrow morning. There is also a 2 line change to the grid simulation method that I forgot in the main PR: using a sparse matrix for grid sim. I also attached the "modular" version of HabitPortfolioConsumerType in ConsumptionSaving.init so that it is imported by default.
Turned off permanent income growth, turned down beta and rho, added analysis/description for 3 of 4 experiments. Fourth is pending.
One target changed because default gridpoints changed. I turned up CRRA on HabitPortfolio test so that the Share test wouldn't be at 1. Also commented out a line in a calibration notebook so that it doesn't run forever when tested.
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
ClippedBilinearInterpwrapper prevents ShareFunc extrapolation errors;Share_optmonotonized to eliminate FOC noise at low wealth;BilinearInterpon regular (m,h) grid replacesCurvilinear2DInterpto remove oscillations at the constrained/interior kinkDiscFac=0.849calibrated via bisection to match habit-only model's pre-retirement wealth[lifecycle-asset-accumulation],[tragic-wealth-shock-at-age-67], etc.), consistentstop_dead=False/np.meanacross all simulations, explanatory text on equity premium effects and human wealth rebalancingcalibrate_DiscFac.ipynb— standalone bisection search for the DiscFac that matches retirement wealth across modelsTest plan
examples/ConsHabitModel/HabitPortfolioConsumerType.ipynbend-to-endexamples/ConsHabitModel/calibrate_DiscFac.ipynbto verify DiscFac convergencetests/ConsumptionSaving/test_ConsHabitPortfolioModel.py