Solving discrete time heterogeneous agent models with aggregate risk and many idiosyncratic states by perturbation

The problem is that once a heterogeneous-agent model has aggregate risk, the cross-sectional distribution must enter the aggregate state vector, and a rich individual planning problem makes that state vector too large to handle by brute force. Extending perturbation to these models “relies on writing the model in the form of a nonlinear difference equation that is function-valued instead of vector-valued (as in representative agent models),” linearized around the stationary equilibrium without aggregate risk. The two functionals that enter are the distribution of agents over idiosyncratic states and the value or policy function describing optimal individual behaviour, and the authors describe them as “replacements for the aggregate capital accumulation and consumption Euler equation in representative agent models.” The payoff of this framing is precise: “These replacements allow us to maintain all nonlinearity with respect to microeconomic shocks — yet obtaining a model that is linear in aggregate variables.” The binding constraint is then stated plainly: “the key practical issue is how to approximate the functionals involved because they need to be replaced by finite-dimensional objects for the actual computation of the model’s dynamics… The curse of dimensionality implies that it is hard to come up with a small enough finite-dimensional representation of the distribution function and the value/policy function without having any a priori knowledge of their shape.”

The paper’s own worked example: two assets plus idiosyncratic income at 50 × 50 × 9 grid points gives 22,500-entry vectors for both the distribution and the value function, Jacobian blocks larger than 45,000 × 45,000, and more than 7 GB per block stored as a full double-precision matrix. The authors add that even at that resolution “the precision is at the lower bound of what one would like to have.” The cost decomposes into two distinct problems: computing the Jacobian is “only quadratic in the number of grid points” and can be sped up by automatic differentiation, but “the matrix to be stored remains large,” and the qz-decomposition and generalized eigenvalue calculation “become very time-consuming (cubic in the number of grid points).”

The reduction happens after the stationary equilibrium is solved and before the Jacobian is computed — earlier than Reiter’s and Ahn et al.’s singular-value-decomposition reductions, later than Reiter’s and Winberry’s parametric approximations. The authors characterize the two existing camps and their costs. Reiter (2009) with splines and Winberry (2018) with parametric distribution families “rely on achieving dimensionality reduction ex ante, before solving for the stationary equilibrium, and hence impose a numerical constraint on this solution”; the downside is that “they might impose tight restrictions on the value function and distribution in the stationary equilibrium” and “no longer allow us to represent the Bellman equation and the distribution dynamics by conveniently linear systems.” Ahn et al. (2017) instead move to continuous time to make Jacobians sparse and then apply singular value decomposition, and — a feature the present paper adopts — perturb deviations of value and distribution functions from their stationary counterparts rather than the functions themselves, which “decouples the number of perturbed parameters from the number of parameters used in the approximation of the functions in the stationary equilibrium.” Against Ahn et al., the claimed advantages are three: avoiding the calculation of a very large Jacobian, applicability to discrete-time models “where the Jacobian would otherwise be too nonsparse to be efficiently stored in a PC’s memory,” and feasibility of second-order or higher perturbation. A further benefit is stated as a property of workflow rather than of accuracy: “the stationary equilibrium can be computed without taking into account that the goal is to solve for aggregate dynamics in the end.”

Because the magnitude of a DCT coefficient has a direct interpretation as that basis polynomial’s contribution to fit, which gives a principled ordering for what to keep. The DCT of a data array “yields the coefficients of the fitted (multidimensional) Chebyshev polynomial, where the polynomial is constructed such that the tensor grid for s is mapped to the Chebyshev knots,” and, citing Hu and Yu (1998), “the larger (in absolute value) a coefficient Θ(i) is, the more important is its corresponding Chebyshev polynomial for fitting ν.” The retained set I is chosen as the smallest set whose coefficients carry a target fraction of the total squared norm — the “energy” — and the inverse transform of the truncated coefficient vector “is the closest one to ν in a least squares sense among all potential inverse discrete cosine transforms of arrays of the same level of sparseness.” The crucial construction is that deviations are added to the full stationary coefficient array rather than replacing it, so at zero deviation the method “fully recovers the stationary equilibrium value function at the same precision as is used in the computation of the stationary equilibrium, that is, without creating any approximation error irrespective of the degree of sparseness.”

Yes, and by a margin that becomes decisive as the coefficient count falls — below 35 retained coefficients the non-adaptive rule breaks the model outright. The comparison is against retaining coefficients corresponding to a complete polynomial of order N, which the authors concede “has a somewhat stronger theoretical underpinning (being a Taylor expansion).” At 101 coefficients the two are equivalent (max absolute difference of log capital 0.08 × 10⁻⁸ versus 0.10 × 10⁻⁸). At 41 coefficients the complete-polynomial rule is off by 37.37 × 10⁻⁸ max against 0.46 × 10⁻⁸ for DCT. And “for less than 35 retained coefficients, the selection based on forming a complete polynomial of given order yields such a bad approximation that we get a violation of the Blanchard-Kahn condition and the model fails to solve.” The stated reason is structural: “across different income states, the policy functions are relatively similar in the stationary equilibrium (think: one is an affine transformation of the other); the DCT method detects this, and this remains true even when prices change after a shock.”

As accurately as the unreduced Reiter method and more than 240 times faster than the original Krusell-Smith algorithm, while remaining the less accurate of the two on Den Haan’s test. The calibration is the JEDC comparison project’s: quarterly periods, β = 0.99, relative risk aversion ξ = 1, 2.5% quarterly depreciation, two-state Markov chains for idiosyncratic and aggregate productivity, 100 grid points for idiosyncratic capital, 1000-period simulations with productivity draws held fixed across methods. Twenty-five DCT coefficients conserve 99.99% of the energy. Against the Krusell-Smith algorithm, the mean absolute difference in log capital is 0.0324% for the reduced method and 0.0324% for the unreduced one; reduced against unreduced is 0.0003% mean, 0.0012% max. On the Den Haan metric the reduced method gives 0.0100% mean and 0.0191% max, the unreduced 0.0102% and 0.0193%, and the Krusell-Smith algorithm — “the most accurate algorithm in Den Haan, Judd, and Juillard (2010)” — 0.0051% and 0.0131%. The paper carries the scope condition rather than burying it: “as can be expected for a first-order perturbation, solution quality deteriorates with the variance of shocks.” Run times on a Dell laptop (Intel i7-7500U, 2.70 GHz, 4 cores; Matlab) are 0.38 seconds reduced, 1.19 unreduced, 91.61 for Krusell-Smith, with 7.05 seconds for the stationary equilibrium.

That the method solves a model the unreduced approach cannot solve at all — three dimensions of household heterogeneity, 120,000 states and 240,000 controls, reduced to 236 states and 1427 controls. The economy has a liquid nominal asset and illiquid capital, Rotemberg price adjustment costs and hence a Phillips curve, a Taylor rule with smoothing, a debt rule, entrepreneur households who receive pure rents, and capital adjustment costs; it follows Bayer, Luetticke, Pham-Dao and Tjaden (2019) and Luetticke (2018). The household problem uses 100 grid points for each asset and 12 for productivity. With those dimensions “it is infeasible to solve for the aggregate dynamics of the model on the full histogram.” Retaining DCT coefficients up to 99.9999% cumulative energy yields the 1427 controls; the fixed copula yields the 236 states; the solve takes 326 seconds plus 1311 seconds for the stationary equilibrium (roughly 5 minutes plus 22), in Julia, chosen because “with the richer model, some of the histogram entries contain very little mass and numerical derivatives become less precise” and a Julia automatic-differentiation package helps. Den Haan errors for the frictionless calibration are 0.033% mean and 0.092% max for capital, 0.081% and 0.617% for bonds; the authors note bonds “are only 10% of the capital stock in the steady state so that, relative to capital or output, the errors are comparable to the errors for capital.”

Business cycle moments, market clearing and the Sharpe ratio all move little across specifications, and in one case tightening the approximation makes things marginally worse. Across baseline, a version retaining more DCT coefficients (energy dropping only by 10⁻⁷ instead of 10⁻⁶), and a version that also perturbs 50 copula DCT coefficients, output volatility is 1.40, 1.42 and 1.31, consumption volatility 1.39, 1.42 and 1.28, investment volatility 4.49, 4.49 and 4.56. Mean absolute deviations from exact asset-market clearing are 0.03% on capital and 0.05% on bonds, with maxima of 0.29% and 0.65%. The authors report the one counterintuitive result rather than suppressing it: “Surprisingly, perturbing also the copula worsens the approximation quality marginally.” The Sharpe ratio is 0.70, 0.67 and 0.68 across the three; set against the Jordà, Knoll, Kuvshinov, Schularick and Taylor (2019) range of 0.6 for housing to 0.25 for equities, the paper’s claim is deliberately modest — the model “gets a long way in terms of being close to the observed Sharpe-ratios,” and “has somewhat to[o] stable asset returns if anything.”

When the shock hits the joint distribution directly rather than hitting everyone alike — idiosyncratic income-uncertainty shocks are the paper’s own counterexample. Under TFP shocks the authors reason that “as all households are similarly affected by the TFP shocks, there is no strong a priori reason for the copula to vary much over the cycle,” and measure the Jensen-Shannon distance between the reduced and unreduced joint asset-income distributions at 0.0005 — “an order of magnitude smaller than the distance between either solution and the stationary equilibrium distribution,” and “negligibly small,” with virtually no difference in the capital series. Under uncertainty shocks, which “affect the joint distribution of assets and income directly, so that the fixed-copula assumption has more potential to introduce approximation errors,” the gap “attains a significant order of magnitude,” and there is “some difference in the fluctuations of the capital stock that the model implies.” The fix is within the same framework: “perturbing the most important 41 coefficients (out of possible 2100) of the DCT of the copula virtually eliminates the already small difference to the full Reiter solution.”

The retained-coefficient set is a heuristic with no guarantee away from the stationary equilibrium, and the paper says so explicitly. The disadvantage “is that it is not guaranteed that the coefficients of the expansion around the stationary equilibrium value function that are shrunk to zero are unimportant for the shape of the value function outside the stationary equilibrium. They are only unimportant in the stationary equilibrium (and hence would have been left out in procedures that reduce the dimensionality entirely ex ante).” The remedy offered is empirical rather than theoretical: check with Den Haan’s simulation test. A footnote adds both the researcher’s obligation — “the researcher should check the robustness of her findings to the choice of the degree of sparseness” — and a tu quoque: the SVD-based reduction of Ahn et al. also requires choosing “the minimal singular value that is retained.” The solution remains a local, first-order (or second-order) approximation in aggregate variables throughout; nonlinearity is preserved only with respect to idiosyncratic shocks.

A proof of concept, explicitly labelled as such. Because the reduction keeps the derivative count low, second-order perturbation via Schmitt-Grohé and Uribe (2004) becomes feasible: for the Krusell-Smith model this “requires to calculate roughly 88 times the number of derivatives as for the first-order perturbation (in total 30,450),” parallelizable across cores, with the qz-decomposition itself not growing. The paper shows the impulse response of capital to a large TFP shock (ten standard deviations) and the ergodic capital distribution under both orders. Its own verdict: “We view this primarily as a proof-of-concept. For practical applications, one will need to further decrease the number of derivatives to be calculated by exploiting the economic structure of the problem, where, for example, the law of motion for the distributions is linear in the distribution.”

One, and it is presented as an illustration of what the two-asset model can do rather than as the paper’s finding: an idiosyncratic uncertainty shock is expansionary in the one-asset model and recessionary in the two-asset model. In response to a 54% increase in the standard deviation of idiosyncratic productivity, consumption falls in both models as households raise precautionary savings. “In the Krusell and Smith model, higher savings translate one-for-one into capital, which leads to an economic expansion. In the two-asset model, by contrast, households prefer to hold more liquid portfolios. They sell illiquid capital to save more in liquid assets. Higher uncertainty therefore causes a simultaneous fall in consumption, investment, and output. The recessionary effect is further amplified through sticky prices, which makes the economy demand-driven in the short run.” The authors tie the illiquidity premium that drives this to the “wealthy hand-to-mouth” concept of Kaplan and Violante (2014) and refer the reader to Bayer et al. (2019) for the portfolio-rebalancing channel itself.

It claims a faster and more scalable route to the same answers, not better answers. The conclusion restates the method as “an extension of Reiter’s method” whose reduction is achieved “by ’lossy compression’ of the value functions, which are control variables of the system, and by approximating the dynamics of the multidimensional distribution of individual characteristics by a distribution with an (almost) fixed copula and varying marginals.” On the Krusell-Smith benchmark it is “equally as precise as Reiter’s standard approach but faster” and “faster and slightly less precise than the Krusell and Smith algorithm in our example” — the “in our example” is the authors’ own qualifier. For the richer two-asset model no accuracy comparison against an unreduced discrete-time benchmark is possible, because none exists; what is shown instead is that the Den Haan test passes, that business cycle properties move little when more DCT coefficients are perturbed, and that the model produces “realistic asset return premia and a reasonably good approximation of asset market clearing.” Codes are provided with the paper, including a Python translation into the HARK toolkit.