Dual-measure simulation#

Dual-measure (P/Q) Monte Carlo simulation mixin for HARK.

The Harmenberg (2021) neutral measure Q reweights permanent income shock probabilities by psi/E[psi], eliminating permanent-income sampling noise from aggregate consumption estimates (3–17x variance reduction).

This module provides a mixin class DualMeasureMixin that, when composed with any IndShockConsumerType subclass, runs the standard P-measure simulation pipeline alongside a parallel Q-measure state update in a single pass. Markov transitions, mortality draws, and base uniform random numbers are shared between the two measures; only the shock magnitudes (and consequently pLvl, mNrm, cNrm, aNrm) diverge.

Usage:

from HARK.dual_measure import DualMeasureMixin
from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType

class DualIndShock(DualMeasureMixin, IndShockConsumerType):
    pass

agent = DualIndShock(**params)
agent.solve()
agent.setup_Q_measure()  # auto-computes IncShkDstn_Q
agent.initialize_sim()
agent.simulate()

P_cNrm = agent.history['cNrm']
Q_cNrm = agent.history_Q['cNrm']
Reference: Harmenberg (2021), “Aggregation with a permanent income shock”,

Journal of Economic Dynamics and Control.

HARK.dual_measure.make_Q_measure_dstn(dstn, warn=True)#

Reweight a DiscreteDistribution by psi/E[psi] (Harmenberg neutral measure).

Parameters:

  • dstn (DiscreteDistribution) – Joint (PermShk, TranShk) distribution under the physical measure. dstn.atoms[0] must be the permanent shock values.

  • warn (bool) – Whether to warn when no reweighting is possible. Callers reweighting a single distribution want the warning and get it by default. setup_Q_measure passes False because it maps this over every period of a lifecycle, where degenerate periods are the normal case rather than a symptom: retirement periods are built with n_approx_Perm = 1, so P equals Q there by construction. Warning once per period made a stock init_lifecycle emit 25 identical warnings from one call, which buries the aggregate warning that does mean something. It reports the count instead.

Returns:

New distribution with Q-measure probabilities and the same atoms. If the permanent shock has zero variance, or a non-positive mean, there is no neutral measure to construct and the original distribution is returned unchanged: the caller gets a Q measure identical to P, which is a no-op rather than the variance reduction the reweighting is asked for.

Return type:

DiscreteDistribution

class HARK.dual_measure.DualMeasureMixin#

Bases: object

Mixin that adds Harmenberg neutral-measure (Q) parallel tracking.

Compose with any IndShockConsumerType subclass via MRO:

class DualAgent(DualMeasureMixin, IndShockConsumerType):
    pass

When dual_measure=True (set by setup_Q_measure()), sim_one_period() runs the standard P-measure pipeline and then a parallel Q-measure state update that reuses the same mortality draws, Markov transitions, and base uniform random numbers.

Zero impact on base classes: AgentType, IndShockConsumerType, and MarkovConsumerType are not modified. The mixin overrides sim_one_period and simulate via MRO. Neither reimplements the P-pipeline: simulate delegates each period to super().simulate(1), and sim_one_period runs the base class’s own _sim_period_prologue() and _sim_period_epilogue() around the Q-step.

dual_measure = False#
setup_Q_measure()#

Auto-compute IncShkDstn_Q from IncShkDstn and enable dual mode.

For each period’s income shock distribution, the Q-measure reweights the probability mass by psi/E[psi]. The atoms (shock values) are unchanged; only their sampling probabilities differ.

Also turns on _cache_base_shock_draws so that get_shocks() records the uniforms the P-draw consumed, and registers IncShkDstn_Q with self.distributions so reset_rng rewinds it alongside the P-side. Both are needed for the shared base draws this class documents: without the flag the Q-side draws independently, and without the registration the Q generators keep advancing across initialize_sim() calls, so the coupling holds on the first run and decays afterwards. Set _cache_base_shock_draws = False after this call to get independent Q draws instead.

Warns when no distribution admits a neutral measure (every permanent shock degenerate), because dual mode then runs the whole Q pipeline to reproduce the P answer at twice the cost and none of the variance reduction. dual_measure is still set: the Q pipeline is well defined in that case, just not useful.

initialize_sim()#

Extend: allocate Q-state arrays after the P-pipeline initializes.

clear_history_Q()#

Allocate NaN-filled Q-history arrays for every tracked variable.

post_state_hook()#

Extension point invoked between get_states() and get_controls() inside this mixin’s sim_one_period.

The default does nothing (beyond deferring to a base-class hook of the same name, if one ever exists), so composing this mixin changes no behavior. Cooperating mixins (e.g. a pLvl-normalization mixin) can override it to adjust states before controls are computed.

sim_one_period()#

Run the P-pipeline, then the Q-pipeline before time advancement.

Calling super().sim_one_period() and appending the Q-step does not work: the base class advances t_age and t_cycle at the end, and the Q-pipeline needs them at the pre-increment values the P-pipeline used. That is why the base class exposes the two halves separately, so the Q-step can sit between them. Everything before and after it is the base class’s own code, reached through the MRO, not a copy of it living here.

simulate(sim_periods=None)#

Extend: record Q-history alongside the base class’s P-history.

The P side is delegated to super().simulate() one period at a time rather than reimplemented here. An earlier version copied the base recording loop and dropped its final else branch, so any tracked variable that is a plain attribute instead of a key in state_now/shocks/controls (MPCnow is the common one) was silently left as NaN whenever dual mode was on. Delegating keeps that class of drift from recurring: turning dual mode on cannot change what the P pipeline records, by construction.

aggregate_Q(var='cNrm', burn=0, N=None, E_pLvl=None, pLvl_factor=None)#

Compute level-aggregate consumption from Q-measure history.

The Harmenberg identity gives:

C_P(t) = N * E_P[p] * F(t) * mean_Q(cNrm_Q(t))

where F(t) = pLvl_factor(t) tracks how E[p_t]/E[p_ss] evolves (equals 1 in a stationary economy).

Parameters:

  • var (str) – Variable name in history_Q to aggregate. Typically 'cNrm'. The variable must be in normalized (per-unit- permanent-income) space.

  • burn (int) – Burn-in periods to skip.

  • N (int or None) – Agent count for level scaling. Defaults to self.AgentCount.

  • E_pLvl (float or None) – Steady-state E[p]. If None, estimated empirically from the P-measure history.

  • pLvl_factor (np.ndarray or None) – Scaling E[p_t]/E[p_ss], one entry per simulated period. Pass the full, un-burned series of length self.T_sim: this method applies [burn:] itself, so a pre-burned array of length T_sim - burn would be trimmed a second time. If None, assumed to be 1 for all periods (stationary economy).

Returns:

Level-aggregate time series.

Return type:

np.ndarray of shape (T - burn,)

HARK.dual_measure.compute_mean_pLvl(agent, g=None)#

Analytical steady-state E[pLvl] for an infinite-horizon HARK agent.

Computes the ergodic cross-sectional mean of permanent income accounting for mortality-driven turnover and permanent income growth.

Parameters:

  • agent (AgentType) – Must have attributes LivPrb, PermGroFac, T_age (or defaults to 400), and pLogInitMean/pLogInitStd.

  • g (float or None) – Effective permanent income growth factor per period. Defaults to agent.PermGroFac[0] (scalar or [0] element). In models with unemployment, set g = (1-u)*G + u where u is the ergodic unemployment rate and G = PermGroFac.

Returns:

E[pLvl] in the stationary cross-section.

Return type:

float

HARK.dual_measure.compute_pLvl_factor(agent, unemployment_path, g_base=None)#

Compute pLvl_factor(t) = E[p_t] / E[p_ss] along a shock path.

In periods with elevated unemployment, average permanent income growth slows because unemployed agents get PermShk = 1 (no growth) while employed agents grow at PermGroFac. This AR(1) recurrence tracks the deviation from steady state:

F(t+1) = (1 - delta) * g_rec(t) * F(t)

  • [1 - (1 - delta) * g_base]

where delta = effective death rate, g_rec(t) = (1-u_t)*G + u_t, and g_base = (1-u_ss)*G + u_ss.

Parameters:

  • agent (AgentType) – Must have PermGroFac, LivPrb, T_age.

  • unemployment_path (array-like of shape (T,)) – Unemployment rate u_t at each period. For the baseline (no shock), pass a constant array at the ergodic rate.

  • g_base (float or None) – Steady-state growth factor. If None, uses the first entry of unemployment_path to compute it.

Returns:

pLvl_factor time series, starting at 1.0.

Return type:

np.ndarray of shape (T,)