The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
106,815 characters
Beyond Aggregate VARs: A Bayesian Benchmark for HANK Models
\maketitle
\vspace*{-6em}
\begin{center}
\end{center}
\begin{center}
\begin{minipage}{.32\textwidth}
\centering\normalsize Florian \MakeUppercase{Huber}\\[0.25em]
\small \textit{University of Salzburg}
\end{minipage}
\begin{minipage}{.32\textwidth}
\centering\normalsize Gary \MakeUppercase{Koop}\\[0.25em]
\small \textit{University of Strathclyde}
\end{minipage}
\begin{minipage}{.32\textwidth}
\centering\normalsize Christian \MakeUppercase{Matthes}\\[0.25em]
\small \textit{University of Notre Dame}
\end{minipage}
\end{center}
\begin{center}
\small This version: \today
\end{center}
\begin{abstract}
\noindent
Heterogeneous-agent New Keynesian (HANK) models characterize how entire cross-sectional distributions respond to structural shocks. Traditional representative-agent models are routinely disciplined by impulse responses from aggregate vector autoregressions (VARs). HANK models have no comparable established empirical benchmark because they make predictions not only about aggregates, but also about distributions of micro-level data. We propose a Bayesian benchmark that jointly models macroeconomic aggregates and several marginal distributions from repeated cross sections, including distributions observed in different surveys. Our approach can use both standard structural VAR identification approaches on macroeconomic aggregates and identification restrictions imposed on micro-level data. The model delivers a joint posterior of the distributional effects of shocks, without the need for household panel data or a separate first-stage density estimate.
\end{abstract}
\vspace*{1em}
\begin{center}
\begin{minipage}{0.85\textwidth}
\noindent\small\textbf{\sffamily JEL}: C11, C32, D31, E52 \\
\textbf{\sffamily KEYWORDS}: Heterogeneous-agent models; functional VAR; Gaussian mixtures; Bayesian inference; distributional impulse responses; monetary and fiscal policy
\end{minipage}
\end{center}
\vfill\noindent{\footnotesize\textit{Acknowledgements}: We would like to thank Ludwig Straub for very useful comments.}
\thispagestyle{empty}\onehalfspacing\normalsize\newpage
\section{Introduction}\label{sec:introduction}
Heterogeneous-agent New Keynesian (HANK) models predict that aggregate variables and household distributions respond together to structural shocks. Building on the incomplete-markets economies of \citet{krusell1998income}, \citet{kaplan2018monetary} show that household heterogeneity changes how monetary policy works. Most of the consumption response operates through general-equilibrium effects on labor income rather than through intertemporal substitution. Heterogeneity also weakens forward guidance \citep{mckay2016power}, turns redistribution into a transmission channel \citep{auclert2019monetary}, and makes fiscal multipliers depend on the distribution of marginal propensities to consume \citep{auclert2024intertemporal}. \citet{kaplan2018microeconomic} survey this research program.
Vector autoregressions (VARs) are the standard benchmark for aggregate predictions of dynamic stochastic general equilibrium (DSGE) models, but they do not recover the distributional responses that distinguish HANK models from their representative-agent counterparts. This paper delivers such a benchmark, while using both the Bayesian toolkit and VAR methods that are familiar to macroeconomists \citep{canova2007methods}.
We keep the VAR's transparent identification and add the distributions the theory is designed to explain. The Joint Aggregate--Micro Mixture VAR (JAMM-VAR) augments a Bayesian structural VAR (SVAR) with several marginal distributions constructed from repeated cross sections. Because the JAMM-VAR models marginal rather than joint distributions, the cross sections may come from separate surveys: earnings can come from the Current Population Survey (CPS) and consumption from the Consumer Expenditure Survey (CEX) even though the surveys neither follow nor contain the same households. We represent each marginal distribution as a Gaussian mixture, allowing a small number of normal components to approximate fat tails, skewness, and multimodality \citep{marron1992exact,FS_book}. Estimation requires neither household panels nor a preliminary density estimate.
The JAMM-VAR accommodates standard SVAR identification schemes for shocks in the aggregate block, including recursive orderings, sign and narrative restrictions, and instruments included in the system \citep{plagborg2021local}. It also identifies an aggregate shock represented by an innovation to a common latent factor, using restrictions on the factor's distributional and aggregate effects. We call it a common micro shock because cross-sectional information helps identify it, not because it is an idiosyncratic household shock. The joint posterior therefore delivers both the distributional responses to shocks identified within the aggregate SVAR and the aggregate responses to the common micro shock.
The JAMM-VAR links the aggregate and distributional blocks in both directions. Contemporaneous and lagged aggregate variables, together with common latent factors, determine the mixture weights and thus the shape of each distribution. The latent factors capture movements shared across distributions that observable aggregates do not explain. Conversely, lagged cross-sectional quantiles enter the macro block, allowing distributional conditions to shape subsequent aggregate dynamics. This feedback allows nonlinear aggregate dynamics; setting the quantile-feedback coefficients to zero yields a linear aggregate SVAR.
The posterior bands for impulse responses account for uncertainty from estimating the cross-sectional distributions. We draw the mixture components, weights, factors, and VAR coefficients from a joint posterior rather than estimating densities first and treating them as data. Two-step procedures can omit this source of uncertainty.
Our approach uses familiar Bayesian and VAR tools. Shrinkage priors regularize the high-dimensional weight and coefficient blocks, truncated priors impose the identifying restrictions, and information criteria compare specifications. The Gibbs sampler is modular: missing survey waves and stochastic volatility require additional blocks but leave the core algorithm unchanged. \autoref{sec:extensions} develops both extensions. Because the aggregate block remains an SVAR, the model also supports variance decompositions, historical decompositions, and conditional forecasts.
We evaluate the method in two simulation exercises. The first uses the model itself as the data-generating process. We verify that the sampler recovers the weights, densities, and impulse response functions (IRFs) when the model nests the data-generating process. The second exercise is more challenging, but also much more relevant for macroeconomists. We estimate our model on data simulated from a HANK economy that it does not strictly nest. Despite this deliberate misspecification, and conditional on matching the policy-variable impact, the model recovers the broad propagation of a monetary policy shock across aggregates, earnings quantiles, and consumption quantiles in a long simulated sample.
Finally, we apply our approach to U.S. data. We combine a standard policy VAR with CPS earnings and CEX consumption to estimate posterior response targets for identified fiscal and monetary policy shocks and for a common micro shock. We use the monetary-policy targets to evaluate the estimated HANK model of \citet{bayer2024shocks}. The model reproduces the cross-quantile pattern and the magnitude of the earnings responses on impact. Its aggregate responses are too large, and both its aggregate and earnings responses fade too quickly.
Our paper continues a long tradition of using VAR evidence to discipline equilibrium models. \citet{gali1999technology} uses the estimated responses to a long-run-identified technology shock to discriminate between competing business-cycle mechanisms. \citet{christiano2005nominal} choose the structural parameters of a New Keynesian model to match the VAR responses to an identified monetary policy shock, and \citet{altig2011firm} extend this matching strategy to technology and investment-specific shocks. \citet{delnegro2004priors} turn the mapping around and use an equilibrium model as a prior for a VAR. \citet{delnegro2007fit} use the resulting hybrid to measure how far New Keynesian models fall short of VAR fit. \citet{canova2011business} evaluate what sign-restricted VARs can recover from equilibrium models, and \citet{loria2022economic} use a VAR as a common empirical framework for confronting several structural theories with the data. Throughout this tradition the evidence is aggregate. We extend it to distributions.
The closest work is the functional-VAR literature, which jointly models aggregates and cross-sectional distributions. \citet{chang2024heterogeneity} estimate a period-specific sieve representation of each distribution and embed the resulting coefficients in a Bayesian state-space model. They also use functional-VAR responses to discipline heterogeneous-agent models. \citet{chang2024monetary} extend this framework to several marginal distributions drawn from separate data sources. The JAMM-VAR differs in two respects. First, a joint likelihood combines the micro cross sections with the aggregate dynamics, so no estimated functional coefficients are passed to a second stage and the posterior retains the uncertainty from fitting the distributions. Second, the same posterior contains both distributional responses to aggregate shocks identified with standard SVAR restrictions and aggregate responses to a common shock identified through restrictions on micro data.
In related work, \citet{nagasaka2026identifying} uses a two-step procedure. The first step approximates the time series of cross-sectional densities with a finite set of functional principal components, following \citet{chang2024functional}. The second places the resulting functional principal-component loadings and aggregate variables in a mixed autoregression and uses direct effects estimated from microeconometric research designs as prior restrictions to identify aggregate shocks. \citet{nagasaka2026identifying} builds on \citet{matthes2025missing}, who identify aggregate shocks through units' heterogeneous exposure but require a long panel. The density-based approach of \citet{nagasaka2026identifying} works with the repeated cross sections available at the household level. Unlike panel and pseudo-panel approaches, including \citet{baumeister2026havar} and \citet{koop2026pseudo}, we neither model individual transitions nor require stable household groups. Naturally, if panel data were available, one could still use our approach by neglecting the panel dimension of the data. \citet{ettmeier2024functional} discuss this distinction between distributional and household-level questions.
Other work estimates heterogeneous-agent models directly. \citet{liu2023full} develop full-information Bayesian inference that combines aggregate time series with repeated micro cross sections, while \citet{parraalvarez2023estimation} construct the cross-sectional likelihood from the model's Fokker--Planck equation. \citet{auclert2021using} use sequence-space Jacobians to make large heterogeneous-agent models tractable to estimate, and \citet{bayer2024shocks} estimate a HANK model with state-space methods. These approaches estimate structural parameters under the full restrictions of a particular equilibrium model. The JAMM-VAR is complementary: it summarizes joint aggregate--distributional evidence without imposing any one model's structure. The resulting posterior response objects can therefore be used to evaluate several structural models under common measurement and identification choices.
The paper proceeds as follows. \autoref{sec:hank_targets} describes how the method fits into the workflow of a researcher using heterogeneous-agent models. It presents the posterior targets it delivers and a step-by-step protocol for using them. \autoref{sec:econometrics} develops the econometric framework. It covers the mixture representation of the cross sections, the macro block, identification, priors, and the posterior sampler. \autoref{sec:validation} assesses the method on simulated data, first with the model itself as the data-generating process and then with a HANK economy that the model does not nest. \autoref{sec:empirical} presents the U.S. application, namely the distributional effects of identified fiscal and monetary policy shocks and the aggregate effects of a common micro shock. \autoref{sec:hankbench} uses the empirical estimates to benchmark an estimated HANK model. \autoref{sec:extensions} develops two extensions, cross sections observed at only some dates and stochastic volatility. \autoref{sec:conclusions} concludes.
\section{Our Method Delivers Targets for HANK Models}\label{sec:hank_targets}
This section explains how researchers can use the posterior targets of the JAMM-VAR to evaluate heterogeneous-agent models. We focus on HANK models, but the protocol applies more broadly whenever a model generates predictions for distributions observed in repeated cross sections. It applies, for example, to models of firm heterogeneity that predict how the distributions of productivity, investment, or employment respond to aggregate shocks \citep[e.g.,][]{liu2023full,winberry2018method,ottonello2020financial,marcellino2025firm}.
The JAMM-VAR supplies aggregate and distributional response targets for evaluating a HANK model, much as an aggregate SVAR supplies targets for a representative-agent DSGE model. The additional evidence consists of responses of the marginal distributions HANK models are designed to explain. Unlike direct structural estimation, the JAMM-VAR neither imposes the structure of a particular equilibrium model nor estimates household policy rules or deep parameters, although structural models can inform its priors. The targets below are not exhaustive; researchers can add others suited to their application.
\subsection{Posterior Targets}
The posterior delivers at least three sets of targets:
\begin{enumerate}[label=(\roman*), nosep]
\item aggregate responses to shocks identified with the specified SVAR restrictions;
\item responses of marginal distributions, including their quantiles and inequality measures, to the same shocks;
\item aggregate and distributional responses to any shock identified using restrictions on micro data.
\end{enumerate}
Density and quantile IRFs summarize the same underlying distributional response and therefore contain overlapping information. We report both because density IRFs show how probability mass shifts across the support, whereas responses of selected quantiles make magnitudes and cross-quantile patterns easy to compare across horizons and with structural models.
A candidate HANK model should be evaluated against the joint posterior of these response objects, not against point estimates assembled from separate procedures. The joint posterior preserves estimation uncertainty and dependence across the targets. Separate procedures may also use different samples, information sets, or identifying assumptions, producing targets that need not be mutually coherent.
To give some examples where our approach can be useful, in \citet{kaplan2018monetary}, intertemporal substitution accounts for a small share of the consumption response to a monetary policy shock and general-equilibrium labor-income effects account for the rest. The two transmission regimes have different implications for our targets, because direct transmission front-loads the consumption-quantile responses, which jump on impact and decay, while indirect transmission makes them inherit the hump shape and peak timing of the earnings-quantile responses. One diagnostic could then compare the peak horizon and peak size of each consumption-quantile response with those of the corresponding earnings quantiles, with credible bands from one posterior. Because the surveys are separate repeated cross sections, corresponding earnings and consumption quantiles need not contain the same households, so such comparisons are diagnostics to be read through a structural model rather than mechanisms identified by the quantile responses alone. The earnings-heterogeneity channel of \citet{auclert2019monetary} operates through unequal incidence of aggregate labor-income movements, and one implication is the cross-quantile ordering of the earnings responses. After a contractionary shock, the bottom quantiles fall by more than the top, so the P90--P10 spread (the 90th minus the 10th percentile) rises. In \citet{auclert2024intertemporal}, intertemporal marginal propensities to consume (MPCs) determine fiscal multipliers, and high-MPC households spend when income arrives rather than when news arrives. Because the fiscal shock in our application is a news shock, the announcement-versus-realization timing of the consumption-quantile responses provides a timing diagnostic related to the model's intertemporal-MPC profile. Finally, \citet{bayer2024shocks} estimate impulse responses of income and consumption inequality to structural shocks inside a HANK model. Such model responses can be compared with our posterior inequality bands horizon by horizon.
The benchmark does not identify household-level transitions. Repeated cross sections reveal how a marginal distribution changes, not which households move across its regions. Similarly, combining CPS earnings with CEX consumption does not identify the household-level joint distribution of earnings and consumption. The JAMM-VAR instead models the dynamics and comovement of these marginal distributions. The framework could be extended to model joint distributions, but datasets containing the required joint outcomes are relatively rare.
\subsection{A Protocol for HANK Researchers}
We now describe a prototypical workflow for disciplining a heterogeneous-agent model.
\begin{steps}
\item \textbf{Choose the empirical objects.} Select the aggregate variables and the marginal distributions of earnings, consumption, wealth, expectations, or other outcomes that the HANK model seeks to explain.
\item \textbf{Construct repeated cross sections.} Choose a common sample period for the aggregate and cross-sectional data. The cross sections may come from separate surveys and need not contain the same households. \autoref{sub:ext_missing} relaxes the common-sample requirement when cross sections are observed in only some aggregate periods.
\item \textbf{Identify the shocks.} If the comparison uses impulse responses, apply a standard SVAR identification scheme for the aggregate shocks. If desired, state separately the restrictions that identify a shock using the micro data.
\item \textbf{Estimate the joint posterior.} Use posterior draws to construct the selected targets, including impulse responses for aggregates, densities, quantiles, and inequality measures when shocks are identified.
\item \textbf{Make the HANK objects comparable.} Simulate the candidate model using the same variable definitions, transformations, survey observation rules, shock size, and response horizons as in the empirical benchmark.
\item \textbf{Diagnose the model.} Compare the model-implied targets with their posterior distributions and identify the aggregate or distributional margins on which they differ. The comparison need not be limited to responses to identified shocks; it can also use second and higher moments of the aggregate time series or cross-sectional distributions. Existing methods for comparing calibrated-model predictions with statistical estimates can guide this step \citep{smith1993estimating,canova1994statistical,dejong1996bayesian}.
\end{steps}
\section{Econometric Framework}\label{sec:econometrics}
We jointly model $M$ macroeconomic aggregates $\bm{Q}_t$ and $S$ marginal cross-sectional distributions. At date $t$, cross section $s$ contains
\[
\bm{y}_{t,s} = (y_{1t,s}, \dots, y_{n_{t,s}t,s})',
\]
where $n_{t,s}$ is the sample size and $p_s(y_{it,s} \mid \bm{\vartheta}_s)$ is the density for observation $i$. Sample sizes may vary across dates and distributions. The data need not track the same units over time, and different distributions may come from separate surveys. Conditional on the common variables and distribution-specific parameters, we treat observations as independent within each cross section. Cross section $(t,s)$ therefore contributes $\prod_{i=1}^{n_{t,s}}p_s(y_{it,s} \mid \bm{\vartheta}_s)$ to the joint likelihood.
We now provide a stylized description of the most common approach in the literature. We add this description here to make clear how our approach differs. One could approximate each cross-sectional distribution by a univariate Gaussian, $p_s \approx \mathcal{N}(\mu_{t,s}, \sigma^2_{t,s})$, estimated separately at each date. Then stack the macro aggregates $\bm{Q}_t$ with the cross-section-specific means and log variances into
\[
\bm{z}_t = (\bm{Q}'_t, \mu_{t,1}, \log\sigma^2_{t,1}, \dots, \mu_{t,S}, \log\sigma^2_{t,S})'.
\]
A natural statistical model would then be a VAR for $\bm{z}_t$ that links $\bm{Q}_t$ to all cross-sectional distributions:
\begin{equation}
\bm{z}_t = \bm{A}_1 \bm{z}_{t-1} + \dots + \bm{A}_P \bm{z}_{t-P}
+ \bm{\varepsilon}_t, \label{eq:fVAR_simple}
\end{equation}
where $\bm{A}_j$ is a $(M + 2S) \times (M + 2S)$ coefficient matrix and $\bm{\varepsilon}_t \sim \mathcal{N}(\bm{0}, \bm{\Sigma}_z)$. Responses of $\bm{z}_t$ to shocks driving $\bm{Q}_t$ trace the effects on the first two moments of each distribution. A Gaussian approximation, however, cannot capture the pronounced asymmetry and long right tails of income and wealth distributions. Adding a small number of higher moments relaxes the Gaussian approximation, but it still compresses each distribution into a few chosen statistics: different shifts in probability mass can produce the same responses of those statistics. Adding many moments or quantiles increases the dimension of the VAR, and separately modeled quantiles need not remain ordered.
\citet{chang2024heterogeneity} address these limitations with a flexible spline representation. In their comparisons, the functional approach yields substantially tighter response bands than a VAR containing selected percentiles. A VAR containing inequality measures also produces implausible long-horizon responses. They approximate each time-varying cross-sectional density with spline basis functions and include the resulting coefficients in a VAR similar to \autoref{eq:fVAR_simple}. The JAMM-VAR instead estimates a parametric mixture representation jointly with the aggregate dynamics using the underlying micro observations. The mixture yields densities and quantiles directly and accommodates several marginal distributions without adding every basis coefficient to the VAR state.
\subsection{Marginal Distributions from Repeated Cross Sections}
We represent each cross-sectional distribution with a finite Gaussian mixture \citep{FS_book}. The fully parametric mixture requires no preliminary kernel-density estimate and is estimated jointly with the aggregate dynamics. The functional form is common across distributions $s \in \{1, \dots, S\}$, but the component parameters and number of components may differ. We therefore present the specification for a generic distribution $s$.
Rather than approximating $p_s$ by a single Gaussian, we use a mixture of $G_s$ Gaussians:
\begin{equation}
p_s(y_{it,s} \mid \bm{\vartheta}_s) \approx \sum_{g=1}^{G_s}
w_{tg,s} \mathcal{N}(y_{it,s} \mid \mu_{g,s}, \sigma^2_{g,s}),
\quad i = 1, \dots, n_{t,s}, \label{eq: mixture_1}
\end{equation}
where $w_{tg,s}$, $\mu_{g,s}$, and $\sigma^2_{g,s}$ are the component weights, means, and variances for cross section $s$. We allow the number of components $G_s$ to differ across cross sections. Only the weights vary over time. Component locations and scales are fixed. The weights follow a multinomial logistic specification:
\begin{equation}
w_{tg,s} = \frac{e^{\eta_{tg,s}}}{\sum_{j=1}^{G_s - 1}
e^{\eta_{tj,s}} + 1}, \quad g = 1, \dots, G_s - 1, \label{eq:softmax}
\end{equation}
with component $G_s$ as the reference, meaning $\eta_{tG_s,s}=0$, so $w_{tG_s,s} = \big(\sum_{j=1}^{G_s-1} e^{\eta_{tj,s}} + 1\big)^{-1}$. The unnormalized log weight $\eta_{tg,s}$ is a linear function of contemporaneous and lagged macro aggregates and common latent factors:
\begin{equation}
\eta_{tg,s} = \beta_{0g,s} + \bm{\beta}'_{g,s} \bm{x}_t
+ \bm{\lambda}'_{g,s} \bm{f}_t.
\label{eq:log_weights}
\end{equation}
Here $\bm{x}_t = (\bm{Q}'_t, \dots, \bm{Q}'_{t-P})'$, $\beta_{0g,s}$ is a component- and cross-section-specific intercept, and $\bm{\beta}_{g,s}$ is a $(P+1)M$-dimensional vector of macro loadings for component $g$ in cross section $s$. These loadings vary freely across $g$ and $s$, allowing each distribution to respond differently to the same macroeconomic conditions.
The $R$ static factors $\bm{f}_t \sim \mathcal{N}(\bm{0}, \bm{I}_R)$ are common to all $S$ distributions, with component- and distribution-specific loadings $\bm{\lambda}_{g,s}$. They induce comovement by shifting the mixture weights of each distribution on which they load. Absent additional restrictions, the factors are statistical objects without a structural interpretation. \autoref{sub:identification} describes the restrictions on their loadings and realizations that identify selected factor innovations as structural shocks.
The mixture likelihood is invariant to permutations of the component labels within each distribution, creating the standard label-switching problem. We resolve it by imposing the ordering $\mu_{1,s} < \dots < \mu_{G_s,s}$ for each $s$ throughout estimation.
\paragraph{A two-component example.}
To make the mixture mechanism concrete, consider one distribution with two components ($G=2$) and one macro aggregate $Q_t$. We suppress the distribution subscript $s$. The density is
\begin{equation*}
p(y_{it}\mid\bm \vartheta)
= w_{t1}\mathcal{N}(y_{it}\mid\mu_1,\sigma_1^2)
+ (1-w_{t1})\mathcal{N}(y_{it}\mid\mu_2,\sigma_2^2).
\end{equation*}
The ordering restriction $\mu_1 < \mu_2$ means that component $1$ is centered at lower outcomes and component $2$ at higher outcomes. Because the component means and variances are fixed, changes in $w_{t1}$ generate all time variation in the distribution.
To see how $Q_t$ changes the distribution through $w_{t1}$, consider a logit specification without latent factors:
\begin{equation*}
w_{t1} = \frac{\exp(\eta_{t1})}{1+\exp(\eta_{t1})},
\qquad
\eta_{t1} = \alpha_1 + \beta_1 Q_t.
\end{equation*}
Suppose $\beta_1<0$. Then
\begin{equation*}
\frac{\partial w_{t1}}{\partial Q_t}
= \beta_1 w_{t1}(1-w_{t1}) < 0.
\end{equation*}
An increase in $Q_t$ therefore lowers $w_{t1}$ and shifts probability mass toward the higher-outcome component. The shift is largest when the two components have equal weight ($w_{t1}=0.5$) and approaches zero as either component's weight approaches one.
The implied cross-sectional mean is
\begin{equation*}
\mathbb{E}_t(y_{it})
= w_{t1}\mu_1 + (1-w_{t1})\mu_2.
\end{equation*}
Differentiating with respect to $Q_t$ gives
\begin{equation*}
\frac{\partial \mathbb{E}_t(y_{it})}{\partial Q_t}
= (\mu_1-\mu_2)\beta_1 w_{t1}(1-w_{t1}) > 0.
\end{equation*}
The derivative is positive because $\mu_1-\mu_2$ and $\beta_1$ are both negative. The shift toward the higher-outcome component therefore raises the cross-sectional mean. The cross-sectional variance also depends on the mixture weight:
\begin{equation*}
\mathrm{Var}_t(y_{it})
= w_{t1}\sigma_1^2 + (1-w_{t1})\sigma_2^2
+ w_{t1}(1-w_{t1})(\mu_1-\mu_2)^2.
\end{equation*}
Thus changes in $w_{t1}$ can alter both the location and dispersion of the cross-sectional distribution even though the component means and variances remain fixed.
\subsection{The Aggregate Time Series Model}
The macro block is a structural VAR($P$) for $\bm Q_t$ augmented by two links to the cross-sectional distributions: lagged quantiles and the common latent factors. The lagged quantiles allow the distributions to affect subsequent aggregate dynamics, while the factors capture contemporaneous comovement between the macro and mixture blocks. Let $\bm A_0$ denote the $M \times M$ contemporaneous coefficient matrix. The structural form is
\begin{equation}
\bm A_0 \bm Q_t = \bm c + \bm A_1 \bm Q_{t-1} + \dots + \bm A_P \bm Q_{t-P} + \sum_{s=1}^{S} \sum_{r \in \mathcal{R}} \bm \alpha_{r,s} \mathcal{Q}_r(\bm y_{t-1,s}) + \bm \Lambda_q \bm f_t + \bm u_t, \quad \bm u_t \sim \mathcal{N}(\bm 0, \bm D). \label{eq: VAR}
\end{equation}
We normalize the diagonal of $\bm A_0$ to one and parameterize it as $\bm A_0 = \bm I_M - \bm W$, where $\bm W$ has a zero diagonal and contains the free contemporaneous coefficients. The structural-shock covariance matrix is $\bm D = \text{diag}(d_1,\dots,d_M)$. Because $\bm D$ is diagonal, the elements of $\bm u_t$ are mutually uncorrelated. Premultiplying the structural system by $\bm A_0^{-1}$ yields the reduced form
\begin{align}
\bm Q_t
&=\bm A_0^{-1}\bm c
+\sum_{j=1}^{P}\bm A_0^{-1}\bm A_j\bm Q_{t-j}
+\bm A_0^{-1}\sum_{s=1}^{S}\sum_{r\in\mathcal{R}}
\bm\alpha_{r,s}\mathcal{Q}_r(\bm y_{t-1,s}) \notag\\
&\quad+\bm A_0^{-1}\bm\Lambda_q\bm f_t
+\bm\varepsilon_t,
\qquad
\bm\varepsilon_t=\bm A_0^{-1}\bm u_t,
\qquad
\operatorname{Var}(\bm\varepsilon_t)=\bm A_0^{-1}\bm D\bm A_0^{-\prime}.
\label{eq:VAR_reduced}
\end{align}
The reduced-form innovations $\bm\varepsilon_t$ are correlated across equations. The data pin down only their covariance $\bm\Sigma=\bm A_0^{-1}\bm D\bm A_0^{-\prime}$. The restrictions in \autoref{sub:identification} resolve the rotation from $\bm\varepsilon_t$ to the orthogonal structural shocks $\bm u_t$, exactly as in a standard SVAR. The vector $\bm c$ collects the intercepts, the matrices $\bm A_1,\dots,\bm A_P$ are the $M \times M$ structural lag coefficients, $\mathcal{R}$ is a finite set of $n_q$ quantile levels in $(0,1)$, and $\mathcal{Q}_r(\bm y_{t-1,s})$ is the $r$-quantile of lagged cross section $s$.
The static factors enter both blocks, with macro loadings $\bm \Lambda_q$. Absent additional restrictions, a factor has no structural interpretation. It links the macro series to the micro distributions. In the empirical application, restrictions on $\bm A_0$ (equivalently $\bm W$) identify the aggregate shocks, and restrictions on the loadings $\bm \Lambda_q$ and $\bm \lambda_{g,s}$ label the common micro shock. In our implementation, we think of the aggregate structural shocks as being elements of $\bm u_t$.
The lagged quantiles make the joint system nonlinear. The log-weight indices in \autoref{eq:log_weights} are linear in the macro variables and factors, but the multinomial-logit map converts these indices into mixture weights nonlinearly. The quantiles implied by those weights then feed back into $\bm Q_t$ through \autoref{eq: VAR}. Impulse responses can therefore depend on a shock's sign and size and on the initial state. Setting $\bm\alpha_{r,s}=\bm 0$ for every $r$ and $s$ and $\bm \Lambda_q = \bm 0$ removes this feedback and nests a linear structural VAR for the aggregate block. If we fix $\bm\alpha_{r,s}=\bm 0$ while allowing the quantiles to enter the VAR, the aggregate block is a linear factor-augmented SVAR. The priors in \autoref{sub:priors} center the quantile loadings at zero, so the posterior departs from the linear submodel only when the likelihood supports distributional feedback.
Unlike functional VARs that include basis coefficients or functional principal-component loadings among the dependent variables, our aggregate block keeps selected quantiles and factors on the right-hand side and therefore retains $M$ equations. Excluding contemporaneous coefficients, each aggregate equation contains $MP$ lag coefficients and an intercept. The quantiles and factors add $Sn_q+R$ regressors per equation, or $M(Sn_q+R)$ coefficients to the system, while $\bm A_0$ and $\bm D$ are unchanged.
In the empirical application, $M=4$, $P=4$, $S=2$, $n_q=5$, and $R=1$. The aggregate-lag block alone has $17$ regressors per equation, including the intercept. The ten quantiles and one factor add $11$ regressors per equation, or $44$ coefficients across the system. Excluding intercepts, the aggregate block therefore contains $64$ lag coefficients, $40$ quantile loadings, and $4$ factor loadings, for a total of $108$.
A functional VAR that stacks $K$ basis coefficients per distribution among the dependent variables has $M+SK$ equations and $(M+SK)^2P$ lag coefficients. At $K=10$, it becomes a $24$-variable VAR with $2{,}304$ lag coefficients. The $108$ count pertains only to our aggregate block; the JAMM-VAR also estimates the log-weight equations and the mixture-component parameters. The comparison therefore concerns the size of the VAR block rather than the total number of model parameters. Shrinkage priors regularize the quantile, factor, and log-weight coefficients (\autoref{sub:priors}).
\subsection{Structural Identification}\label{sub:identification}
We identify aggregate shocks and the common micro shock separately. The aggregate block supports the same identification schemes as a standard SVAR, including recursive, sign, magnitude, narrative, and instrument-based restrictions. These restrictions identify selected elements of $\bm u_t$ as aggregate shocks. Because the mixture weights depend on contemporaneous and lagged aggregates, the mixture block maps each identified aggregate shock into responses of the marginal distributions.
We identify the common micro shock as an innovation to a selected factor in $\bm f_t$. The loadings $\bm\lambda_{g,s}$ map this innovation into changes in mixture weights across marginal distributions, while $\bm\Lambda_q$ maps it into the aggregate equations. The normalization $\bm f_t\sim\mathcal{N}(\bm 0,\bm I_R)$ fixes the shock's scale. Sign and zero restrictions on both sets of loadings, together with narrative restrictions on selected factor realizations, give the shock an economic interpretation. Because these restrictions impose some impact responses, the empirical analysis distinguishes the imposed impact signs from the subsequent responses learned from the posterior.
For both types of shock, we report posterior responses of the aggregate variables and of each marginal distribution.
\subsection{Prior Distributions}\label{sub:priors}
The model contains many coefficients relative to the short macroeconomic sample, so the priors serve two roles: they regularize estimation and encode the identifying restrictions. We use an ordered Normal--inverse-Gamma prior for the mixture-component means and variances, horseshoe shrinkage for the log-weight and VAR coefficients, inverse-Gamma priors for the structural variances, and Gaussian priors for the free contemporaneous coefficients. Sign and magnitude restrictions truncate the relevant prior support, zero restrictions fix selected coefficients, and narrative restrictions truncate selected factor realizations. Each restriction therefore holds at every posterior draw \citep{baumeister2015sign}.
\subsubsection{Priors on the Mixture Components and Weights}
For each component $g$ of cross section $s$, we use a conjugate prior that links its location and scale parameters:
\begin{equation*}
p(\mu_{g, s}, \sigma_{g, s}^2) \propto p(\mu_{g, s}| \sigma_{g, s}^2) ~ p(\sigma_{g,s}^2)
\end{equation*}
where the marginal prior for $\sigma^2_{g,s}$ and the conditional prior for $\mu_{g,s}$ are
\begin{equation}
\sigma^2_{g,s} \sim \mathcal{IG}(a_0,b_0),\qquad
\mu_{g,s} | \sigma^2_{g,s},\kappa_{0,s}\sim
\mathcal{N}(m_{0,s},\tfrac{\sigma^2_{g,s}}{\kappa_{0,s}}) ~ \mathbb{I}(\mu_{1,s}<\dots<\mu_{G_s,s}),
\label{eq:prior_mu}
\end{equation}
where $\mathcal{IG}$ denotes the inverse-Gamma distribution. The ordering constraint resolves the label-switching problem discussed in \autoref{sec:econometrics}. We set $a_0=b_0=0.01$ and center the component means for cross section $s$ on its pooled sample mean,
\begin{equation*}
m_{0,s}
=
\left(\sum_{t=1}^{T}n_{t,s}\right)^{-1}
\sum_{t=1}^{T}\sum_{i=1}^{n_{t,s}}y_{it,s}.
\end{equation*}
The precision $\kappa_{0,s}$ controls how tightly the component means cluster around $m_{0,s}$: larger values pull them toward the pooled mean, while smaller values permit greater separation across component means. We place a heavy-tailed hyperprior on $\kappa_{0,s}$ and update it within the sampler. Posterior draws can place component means close together or assign negligible weight to some components. Thus $G_s$ is an upper bound on the effective number of distinct components.
For each nonreference component $g=1,\dots,G_s-1$ in cross section $s$, collect the intercept, macro loadings, and factor loadings in
$\bm b_{g,s}=(\beta_{0g,s},\bm\beta_{g,s}',\bm\lambda_{g,s}')'$. Each cross section contains one such log-weight equation for every nonreference component, and the coefficients vary freely across components and cross sections. We regularize this high-dimensional block with a horseshoe prior \citep{carvalho2010}:
\begin{equation}
b_{l,g,s} | \varphi_{l,g,s},\tau_{g,s} \sim
\mathcal{N}(0,\varphi_{l,g,s}^2\tau_{g,s}^2),\qquad
\varphi_{l,g,s}\sim \mathcal{C}^+(0,1),\qquad
\tau_{g,s}\sim \mathcal{C}^+(0,1),
\label{eq:prior_B}
\end{equation}
where $l$ indexes the elements of $\bm b_{g,s}$, $\tau_{g,s}$ is a component-specific global scale, $\varphi_{l,g,s}$ is a coefficient-specific local scale, and $\mathcal{C}^+$ denotes the half-Cauchy distribution. We sample the scale parameters using the inverse-Gamma augmentation of \citet{makalic2016}. The horseshoe strongly shrinks small coefficients toward zero while leaving large coefficients weakly penalized.
Some elements of $\bm b_{g,s}$ carry identifying restrictions. Restrictions on the contemporaneous macro loadings in $\bm\beta_{g,s}$ govern the distributional impact of shocks identified in the aggregate block, while restrictions on the factor loadings in $\bm\lambda_{g,s}$ help identify the common micro shock. We impose these sign and magnitude restrictions by multiplying the joint horseshoe prior for $\bm b_{g,s}$ and its scales by an indicator for the admissible region and normalizing the joint density once. Conditional on an admissible coefficient draw, the indicator and joint normalizer are constant with respect to the scales, so their inverse-Gamma full conditionals are unchanged. In the augmented sampler, the coefficient full conditional is Gaussian truncated to the admissible region. Every posterior draw therefore satisfies the restrictions without a scale-dependent normalizing-constant correction. We apply the same joint-truncation convention to inequality-restricted coefficients in the macro block, treating $d_i$ as part of the joint prior; zero restrictions instead fix coefficients at zero. The empirical application specifies the exact restrictions.
\subsubsection{Priors on the Macro Block}
The structural macro block in \autoref{eq: VAR} contains three parameter groups: the regression coefficients in $\bm\Phi$, the structural variances in $\bm D=\text{diag}(d_1,\dots,d_M)$, and the free contemporaneous coefficients in $\bm W$, where $\bm A_0=\bm I_M-\bm W$. Let $\bm m_t$ collect the $P$ lags of $\bm Q_t$, the lagged cross-sectional quantiles $\mathcal{Q}_r(\bm y_{t-1,s})$, the common factors $\bm f_t$, and a constant. Stacking the corresponding coefficients in $\bm\Phi$ gives
$\bm A_0\bm Q_t=\bm\Phi'\bm m_t+\bm u_t$.
The $i$th column $\bm\phi_i$ is the coefficient vector for equation $i$; it contains the coefficients on aggregate lags, lagged quantiles, factors, and the intercept. Conditional on its structural variance $d_i$, we assign the Gaussian prior
\begin{equation}
\bm\phi_i\mid d_i\sim\mathcal{N}\big(\underline{\bm\phi}_i,
d_i\underline{\bm V}_i\big),
\label{eq:prior_phi}
\end{equation}
where the prior mean $\underline{\bm\phi}_i$ follows the Minnesota convention \citep{litterman1986}. The own first lag is centered at $\delta$; cross-lags, higher-order own lags, quantile and factor loadings, and the intercept are centered at zero. We set $\delta=0.8$ in the empirical application. The prior precision is
\[
\underline{\bm V}_i^{-1}
=
\text{diag}\left(\frac{1}{\tau_i^2\varphi_{l,i}^2}\right),
\]
where $\tau_i$ is an equation-specific global scale and $\varphi_{l,i}$ is a coefficient-specific local scale, both following the horseshoe hierarchy in \autoref{eq:prior_B}. This prior centers each equation on a persistent univariate process and shrinks the distributional-feedback and factor loadings toward zero.
Identifying restrictions apply to the factor loadings in $\bm\Phi$. Following the joint-truncation convention above, sign restrictions truncate the corresponding multivariate Gaussian conditional. A zero restriction instead fixes the selected factor loading at zero and draws the remaining coefficients from their Gaussian conditional given that restriction.
For the structural variances, we use the conditionally conjugate prior
\begin{equation}
d_i\sim\mathcal{IG}\Big(\tfrac{M+2}{2},\tfrac{\hat\sigma^2_i}{2}\Big),
\label{eq:prior_D}
\end{equation}
where $\hat\sigma^2_i$ denotes the residual variance obtained by estimating the reduced-form of the model on an equation-by-equation basis using a unit ridge penalty. The prior is weakly informative, with mean
\[
\mathbb{E}[d_i]=\frac{\hat\sigma^2_i}{M},
\]
and, for $M>2$, variance
\[
\operatorname{Var}(d_i)
=
\frac{2(\hat\sigma^2_i)^2}{M^2(M-2)}.
\]
The free off-diagonal elements of $\bm W$ follow independent Gaussian priors
\begin{equation*}
W_{ij}\sim\mathcal{N}(0,\underline{l}_W),
\end{equation*}
restricted to intervals $[\underline{w}_{ij},\overline{w}_{ij}]$ that encode the identifying restrictions on $\bm A_0$. A sign restriction sets one endpoint to zero; a magnitude restriction sets a nonzero bound calibrated using the corresponding prior in \citet{baumeister2018inference}; and an exclusion restriction sets both endpoints to zero. We set $\underline{l}_W=2$. We normalize the shared factors as $\bm f_t\sim\mathcal{N}(\bm 0,\bm I_R)$, as stated in \autoref{sec:econometrics}. Narrative restrictions impose sign constraints on selected factor realizations.
\subsection{Posterior Simulation}\label{sub:posterior}
In the unaugmented posterior, three features rule out direct conjugate updates: the multinomial-logit likelihood for the log-weight coefficients, the determinant term in the conditional for $\bm A_0$, and the shared factors' entry in every component index. P\'olya--Gamma augmentation \citep{polson2013} yields conditionally Gaussian Gibbs updates for the log-weight coefficients. We update $\bm A_0$ and the shared factors with Metropolis--Hastings steps that use Gaussian proposals and target their exact conditional distributions; all remaining blocks are sampled by Gibbs. The steps below describe one MCMC sweep. \autoref{app:technical} derives the full conditionals and Metropolis--Hastings acceptance ratios.
\begin{steps}
\item \textbf{Component allocation.} For each micro observation $y_{it,s}$, draw a label $z_{it,s}\in\{1,\dots,G_s\}$ from its multinomial full conditional, with probabilities proportional to
\[
w_{tg,s}\mathcal{N}(y_{it,s}\mid\mu_{g,s},\sigma^2_{g,s}),
\]
using the current weights from \autoref{eq:softmax}.
\item \textbf{Component parameters and mean-shrinkage precision.} Conditional on the allocations, the observations assigned to each component form a Gaussian subsample. For each component, draw $\sigma^2_{g,s}$ from its inverse-Gamma conditional given the current mean and then $\mu_{g,s}$ from its Gaussian conditional truncated to the interval between the neighboring means, which preserves $\mu_{1,s}<\dots<\mu_{G_s,s}$. The two draws form an exact Gibbs update of the ordered Normal--inverse-Gamma conditional implied by \autoref{eq:prior_mu}. Then draw the cross-section-specific precision $\kappa_{0,s}$ from its generalized-inverse-Gaussian ($\mathcal{GIG}$) full conditional.
\item \textbf{Log-weight coefficients.} For each cross section $s$ and nonreference component $g=1,\dots,G_s-1$, update $\bm b_{g,s}$ using the one-vs-rest representation of the multinomial-logit likelihood. Introducing $\omega_{tg,s}\sim\mathrm{PG}(n_{t,s},\psi_{tg,s})$ makes the full conditional for $\bm b_{g,s}$ Gaussian \citep{polson2013}, where $n_{t,s}$ is the sample size of cross section $s$ at date $t$ and $\psi_{tg,s}$ is the current one-vs-rest log odds. Draw each $\omega_{tg,s}$ using the saddlepoint approximation of \citet{windle2014}. Then draw $\bm b_{g,s}$ from its Gaussian full conditional, truncated to the admissible region when the coefficients carry sign or magnitude restrictions. Finally, update the local and global horseshoe scales using the inverse-Gamma augmentation of \citet{makalic2016}. Under the joint-truncation convention in \autoref{sub:priors}, these restrictions leave the scale full conditionals unchanged. \autoref{app:technical} gives the full conditional formulas.
\item \textbf{Macro coefficients.} Conditional on $\bm A_0$, $\bm D$, the shared factors, and the lagged quantiles, the structural system separates by equation. For each equation $i$, combining the likelihood with the Gaussian prior in \autoref{eq:prior_phi} yields a Gaussian full conditional for $\bm\phi_i$. Draw unrestricted coefficient vectors directly from this conditional and coefficient vectors with sign-restricted factor loadings from the corresponding truncated multivariate Gaussian. If a factor loading is fixed at zero, set it to zero and draw the remaining coefficients from their Gaussian conditional given that restriction. Finally, update the local and global horseshoe scales using the same inverse-Gamma augmentation as in the log-weight block.
\item \textbf{Structural variances.} Conditional on $\bm\Phi$, $\bm A_0$, the shared factors, and the lagged quantiles, draw each $d_i$ from its inverse-Gamma full conditional. The update combines the prior in \autoref{eq:prior_D}, the structural residual sum of squares for equation $i$, and the quadratic term from the $d_i$-scaled coefficient prior in \autoref{eq:prior_phi}. \autoref{app:technical} gives the exact parameters.
\item \textbf{Contemporaneous matrix.} The full conditional for each row of $\bm A_0=\bm I_M-\bm W$ contains the Jacobian $|\det \bm A_0|^{T}$ and is therefore non-Gaussian. Because $\det \bm A_0$ is linear in a single row, a second-order Laplace expansion adds a rank-one term to the Gaussian precision, which the Sherman--Morrison identity evaluates in $O(M^2)$. Use the resulting approximate Gaussian, truncated to the sign, magnitude, and zero restrictions on $\bm A_0$, as the proposal in a Metropolis--Hastings step targeting the exact conditional. The sign restrictions on the impact responses in $\bm A_0^{-1}$ involve all rows of $\bm A_0$ at once. They enter the target of every row update as an indicator, so a proposed row is rejected whenever the implied impact responses violate them.
\item \textbf{Common factors.} Draw the factor history on a $t$ by $t$ basis. The P\'olya--Gamma representation is valid for the log-weight coefficients because it holds the competing component indices fixed while $\bm b_{g,s}$ is updated. The factors are shared. They enter \emph{every} component index, so moving $\bm f_t$ also moves the log-sum-exp offset of each one-vs-rest representation, and the resulting kernel is no longer Gaussian in $\bm f_t$. We therefore build a Gaussian proposal that pools information from the macro shocks, through $\bm\Lambda_q$ and $\bm D$, with information from the log-weight equations, through $\bm\lambda_{g,s}$ and {the conditional means of the P\'olya--Gamma weights evaluated at the value the proposal is centered on}. We accept the proposal with a Metropolis--Hastings probability computed from the exact multinomial logistic likelihood{, with the proposal density evaluated in both directions}. This leaves the target posterior invariant. Acceptance rates are between $0.92$ and {$0.95$} across chains in the empirical application. At dates with a narrative restriction, the proposal is drawn from the corresponding truncated Gaussian, so the restricted sign of (elements of) $\bm f_t$ is imposed exactly. \autoref{app:technical} gives the exact conditional and the acceptance ratio.
\end{steps}
We discard burn-in draws and retain every $k$th remaining draw. Thinning reduces storage and the cost of the simulation-based impulse responses. We run several independent chains in parallel and pool the thinned draws to compute impulse responses for each identified shock. The reported MCMC budget for each exercise includes burn-in, retained draws, and thinning. For each retained draw, we propagate the estimated system forward from the sample mean of the state to construct the structural impulse responses. The individual mixture parameters mix more slowly than the fitted distributions and impulse responses. We discuss the convergence diagnostics in \autoref{app:mcmc_diagnostics}.
The sampler is modular. The two extensions developed in \autoref{sec:extensions}---allowing cross sections to be observed at only some dates and introducing stochastic volatility---each add or replace one sampling block; the remaining updates are unchanged or modified only through known weights. Because each retained draw contains the full structural system, the posterior can also be used to construct variance decompositions, historical decompositions, and conditional forecasts in addition to the reported impulse responses.
Before turning to the empirical applications, a brief word on computation times is in order. The estimation times are based on running a single chain with an Apple M4 Max processor and include the simulation of the impulse responses. For the HANK exercise, with $T=500$ dates, twelve aggregate series, and two cross sections of $1{,}000$ {simulated} households per date summarized by $500$ quantile points, a chain of $10{,}000$ draws ($5{,}000$ burn-in and $5{,}000$ retained) takes about $35$ minutes. For the controlled DGP of \autoref{sub:rf_dgp}, a chain of $2{,}500$ draws takes about five minutes. For the empirical application of \autoref{sec:empirical}, with $T=71$ dates, a chain of $15{,}000$ draws ($5{,}000$ burn-in and $10{,}000$ retained) takes about $15$ minutes. Additional chains run in parallel on separate CPU cores at little extra cost in elapsed time.
\section{Assessing Our Approach Using Simulated Data}\label{sec:validation}
We assess how well our approach does in two exercises. The first uses a data-generating process (DGP) that the JAMM-VAR nests. We think of this as a low bar to pass, but nonetheless a bar we want to check. The second exercise uses a HANK DGP outside the model class and asks whether the statistical model still recovers dynamics, {cross sections}, and impulse responses. The HANK economy's known shocks and responses provide a controlled laboratory for assessing the JAMM-VAR's identification and specification choices.
\subsection{Our Model as the Data-Generating Process}\label{sub:rf_dgp}
The controlled DGP combines a two-variable VAR(1), two cross sections, and one common standard-normal factor. The VAR includes the lagged median of each cross section and the factor as regressors. At each of $T=250$ dates, each cross section contains $500$ observations drawn from a three-component Gaussian mixture whose weights depend on the contemporaneous and lagged macro state and on the factor. With only two cross sections, the factor cannot be separated from the observed macro state without additional structure. We identify it using sign restrictions of the same form as those in the empirical application. The design contains two structural channels: an innovation to the macro block and an innovation to the common factor. \autoref{app:validation} reports the full parameterization.
Posterior median weights track the true weights almost exactly, and the fitted densities are nearly indistinguishable from the true densities (\autoref{fig:rf_weights,fig:rf_densities} in the appendix). For both structural shocks, the posterior median aggregate, density, and quantile responses track the truth closely. The 90\% credible bands contain the true quantile responses at every percentile and horizon, the true density responses at every point of the support and horizon shown except a small part of the support at impact for the factor shock, and the true aggregate responses at every horizon except the impact response of $Q_{2t}$ to the factor shock (\autoref{fig:rf_dist_quant_macro,fig:rf_dist_quant_micro} in the appendix). We next turn to a HANK model that is not nested in the JAMM-VAR model.
\subsection{HANK Data-Generating Process}
We use the one-asset HANK model of \citet{auclert2021using} as the DGP. Households face incomplete markets, trade a single liquid asset subject to a zero borrowing limit, and supply labor with idiosyncratic productivity that follows an AR(1) in logs. Firms produce with a linear technology, and Rotemberg price adjustment yields a New Keynesian Phillips curve. Monetary policy follows a Taylor rule. The government keeps the supply of real bonds fixed and levies lump-sum taxes to finance interest payments. \autoref{app:hank_calibration} reports the full calibration and discretization.
We solve the model with the sequence-space Jacobian method of \citet{auclert2021using}, using a modified version of their public toolkit.\footnote{\url{https://github.com/shade-econ/sequence-jacobian}.} Two aggregate AR(1) disturbances drive the economy. A shock to the Taylor-rule intercept has persistence $0.61$ and an innovation standard deviation of $25$ basis points; a TFP shock has persistence $0.8$ and an innovation standard deviation of $1$ percent. We discard the first $2{,}000$ quarters of the simulation and retain the next $T=500$ quarters.
The statistical model observes twelve aggregate series: the exogenous policy-rate shock (MP) and the model's output, consumption, labor demand, real interest rate, real wage, inflation, aggregate assets, hours, effective labor, dividends, and taxes. Because MP is the observed Taylor-rule disturbance, we identify the monetary policy shock as the innovation to its equation. TFP is not among the observables, so one of the two aggregate shocks is hidden from the statistical model. This is a deliberate source of misspecification beyond the mixture approximation itself. Aggregate assets equal the fixed bond supply and are therefore constant.
We estimate the model with $G=(6,6)$ components and $R=2$ latent factors. With two latent factors and unrestricted loadings, the factor configuration is identified only up to rotation. We therefore anchor each factor to one equation in the spirit of \citet{geweke1996}. The first factor is excluded from the equation of the policy instrument, and its loading on output is restricted to be negative. The second factor is excluded from the output equation, and its loading on the policy instrument is restricted to be positive. Both factors are excluded from the asset equation, which is constant in this simulation because bonds are in fixed supply. These restrictions are labeling conventions, analogous to the ordering of the mixture means, and match the sparse loading pattern an unrestricted run recovers.
At each date, we draw separate cross sections of consumption and gross labor earnings, each containing $1{,}000$ households. As in the empirical application, each simulated cross section is summarized by $500$ equally weighted quantile points before estimation. Gross labor earnings equal the product of the real wage, individual productivity, and individual hours. We obtain both cross sections by inverse-CDF sampling from the model-implied distribution over productivity and assets. We use one simulated history and fixed random seeds throughout. We compare the estimates with the HANK model's true impulse responses to a one-time $25$-basis-point monetary policy shock.
Because the HANK DGP is not a finite mixture, it has no true mixture components or weights to recover. The estimated components and weights are approximation devices. We therefore evaluate the cross-sectional distributions implied jointly by the components and weights. The fitted six-component mixtures reproduce the broad support, modes, skewness, and tails of consumption and earnings throughout the sample, although the fit is not exact in every bin (\autoref{fig:hank_weights,fig:hank_densities} in the appendix). The impulse-response analysis therefore focuses on changes in the fitted distributions rather than in individual components or weights.
\subsubsection{Responses to the Monetary Policy Shock}
\autoref{fig:hank_irfs} compares the estimated responses to a monetary policy shock with the true HANK responses. Panel~(a) reports responses for the twelve aggregate variables, while panels~(b) and~(c) report responses for the 10th, 25th, 50th, 75th, and 90th percentiles of consumption and earnings.
\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.80\textwidth}
\centering
\caption{Macro IRFs}
\label{fig:hank_macro}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{HANK_multi_results_MP/macro_irfs.pdf}
\end{subfigure}\\[0.4em]
\begin{subfigure}[b]{0.80\textwidth}
\centering
\caption{Quantile IRFs: consumption}
\label{fig:hank_quant_cons}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{HANK_multi_results_MP/consumption_irfs.pdf}
\end{subfigure}\\[0.4em]
\begin{subfigure}[b]{0.80\textwidth}
\centering
\caption{Quantile IRFs: earnings}
\label{fig:hank_quant_earn}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{HANK_multi_results_MP/earnings_irfs.pdf}
\end{subfigure}
\caption*{\footnotesize \textbf{Notes}: Panel (a) reports responses of the twelve macro variables to a monetary policy shock on a $3\times4$ grid. Aggregate assets ($A$) equal the fixed bond supply and appear as a flat line on a {narrow} symmetric scale{ ($\pm 10^{-3}$)}. Panels (b) and (c) report responses of the 10th, 25th, 50th, 75th, and 90th percentiles of consumption and earnings. Navy lines are posterior medians, shaded areas are 68\% and 90\% credible bands, and dashed red lines are the true HANK responses. IRFs are normalized so that the estimated impact response of the policy variable equals the truth at horizon~1.}
\caption{HANK simulation: macro and quantile IRFs}
\label{fig:hank_irfs}
\end{figure}
After normalizing the estimated impact response of the policy variable to match the truth, the aggregate responses in panel~(a) reproduce the signs and decay of the HANK responses. Posterior medians closely track the truth for output, consumption, the real rate, inflation, wages, dividends, and taxes. The labor-market responses are slightly attenuated on impact, but the discrepancies are concentrated at short horizons and the posterior bands contain the true paths.
For consumption (see panel~(b)), the model recovers the negative impact response across the distribution, its decay, and the larger responses in the tails than at the median. Posterior medians nevertheless attenuate the impact contraction at every percentile except the median and return to zero too quickly. Earnings responses, shown in panel~(c), are recovered more closely. The medians reproduce the negative impact response and its increasing absolute magnitude toward the top of the distribution. The main localized error is a small positive overshoot in the true 25th-percentile response at intermediate horizons. The 90\% credible bands, which are wide at short horizons, generally contain the true consumption and earnings paths.
Overall, our approach recovers the broad aggregate and distributional propagation of the HANK economy. The discrepancies discussed above concern the point estimates, summarized by the posterior medians. Once we account for posterior uncertainty, the credible bands generally contain the true aggregate and distributional response paths.
\section{Structural Evidence from U.S. Repeated Cross Sections}\label{sec:empirical}
We combine a four-variable aggregate VAR with CPS earnings and CEX consumption distributions. Earnings and consumption enter as separate marginal distributions linked only through the aggregate state and the latent factor. The CPS and CEX samples therefore need not contain the same households, and we do not link individual records across surveys. We report distributional responses to identified fiscal and monetary policy shocks and aggregate responses to the common micro shock. The monetary policy responses provide the targets for the HANK model comparison in \autoref{sec:hankbench}; the responses to the fiscal and common micro shocks are additional empirical results.
\subsection{Data and Identification}\label{sub:emp_data_spec}
We follow the small-scale macroeconomic VAR in \citet{baumeister2018inference}, adding only the Ben Zeev--Pappa fiscal-news series \citep{benzeev2017}, denoted $\text{BZP}_t$. The aggregate vector contains $M=4$ endogenous variables:
\begin{equation*}
\bm Q_t = (\text{BZP}_t, OG_t, \pi_t, R_t)'.
\end{equation*}
Here $OG_t = 100 \cdot \log(\text{GDPC1}_t / \text{GDPPOT}_t)$ is the output gap based on Congressional Budget Office potential output, $\pi_t$ is annual log-difference inflation in personal consumption expenditures (\texttt{PCECTPI}), and $R_t$ is the federal funds rate (FFR). All series are quarterly, and the sample runs from $1990$Q2 through $2007$Q4. Availability of the BZP series determines the end date, which also keeps the zero lower bound and pandemic periods outside the sample.
The BZP series measures news about future defense spending. \citet{benzeev2017} identify the corresponding shock by requiring it to be orthogonal to current defense spending and to best explain subsequent movements in defense spending. We obtain the series from the harmonized structural-shock dataset maintained by Jonathan Adams\footnote{\url{https://github.com/jonathanjadams/structuralshocks}.} and documented in \citet{adams2025empirical} and \citet{adams2026ricardian}.
We construct the earnings and consumption cross sections as in \citet{chang2024heterogeneity} and \citet{chang2024monetary}. For earnings, we combine the monthly CPS weekly-earnings variable, annualized to a yearly rate, with the CPS employment indicator. Following \citet{chang2024heterogeneity}, we define relative earnings as
\begin{equation*}
z_{i,t} = \frac{\text{annual earnings}_{i,t}}{\tfrac{2}{3}\cdot\text{nominal per-capita GDP}_t},
\end{equation*}
where $2/3$ approximates labor's share of gross domestic product (GDP). The mixture model uses the inverse hyperbolic sine (IHS) transformation,
\begin{equation*}
x_{i,t} = \operatorname{asinh}(z_{i,t}) = \log\left(z_{i,t}+\sqrt{z_{i,t}^{2}+1}\right),
\end{equation*}
which is approximately linear near zero and behaves like $\log z_{i,t}+\log 2$ for large positive values of $z_{i,t}$. The transformation compresses the right tail of the CPS earnings distribution. It is the same transformation used in \citet[eq.~(30)]{chang2024heterogeneity} and \citet[eq.~(26)]{chang2024monetary}. The earnings cross sections consist of employed individuals with reported weekly earnings. They contain no mass at zero, so the earnings distribution describes the intensive margin among the employed. The quantiles of our processed earnings data closely match those of \citet{chang2024heterogeneity}.
For consumption, we use quarterly household expenditures from the CEX. Following \citet{chang2024monetary}, we normalize household consumption by nominal National Income and Product Accounts (NIPA) consumption per person aged 16 or older:
\[
z^{c}_{i,t} = \frac{\text{consumption}_{i,t}}{C^{\text{NIPA}}_t / \text{pop}^{16+}_t}.
\]
The mixture model uses $z^{c}_{i,t}$ directly. Consumption has a less extreme right tail than earnings, so we do not apply the IHS transformation. The quarterly cross sections contain about $12{,}000$ to $16{,}000$ individual observations for earnings and $4{,}500$ to $8{,}200$ households for consumption. We represent each quarterly cross section using $500$ equally weighted quantile points, which closely preserve the shape of the empirical distribution. The common grid prevents differences in survey sample sizes across dates and outcomes from mechanically changing their relative influence on the estimates. The empirical results are conditional on this fixed-grid representation, and the micro likelihood treats the grid points as the cross-sectional observations. Such an approximation helps to lower the influence of extreme outliers in the micro data. Alternatively, one could clean the data from outliers, i.e. trim or winsorize the raw microdata, as in related applications using earnings and consumption data \citep[e.g.,][]{heathcote2010macro,moffitt2012trends,coibion2021consumption}. In our Monte Carlo exercises, we instead directly use simulated micro data.
Both micro variables normalize household outcomes by aggregate quantities: earnings by labor-share-adjusted nominal GDP per capita and consumption by nominal NIPA consumption per person aged 16 or older. We report density, quantile, and inequality responses in these normalized units. An earnings quantile response of $+0.05$ means that the corresponding percentile of $x_{i,t}$ rises by $0.05$. Since $\operatorname{asinh}'(1)=1/\sqrt{2}$, a $0.05$ change in $x_{i,t}$ corresponds locally to a change of about $0.07$ in $z_{i,t}$ when $z_{i,t}=1$, or roughly a 7\% increase in earnings relative to two-thirds of nominal GDP per capita. Consumption is not transformed, so a response of $0.05$ raises the ratio to per-capita NIPA consumption by five percentage points. We do not convert these responses into raw earnings or consumption levels because doing so would also require the response of each aggregate normalizer and, for earnings, inversion of the IHS transformation.
The multinomial logistic weights depend on a constant, contemporaneous $\bm Q_t$, and its first four lags, matching the VAR lag order $P=4$. A single latent factor ($R=1$), denoted $f_t$, captures distributional movements not explained by the macro aggregates and enters the aggregate VAR through the loading vector $\bm\Lambda_q$. We set its loading in the BZP equation to zero, so $f_t$ does not enter that equation contemporaneously, and restrict its loadings in the output-gap and inflation equations to be positive. For each cross section, we order the component means to resolve label switching.
We choose the number of components $G_s$ using the Widely Applicable Information Criterion (WAIC) \citep{watanabe2010,gelman2014waic}, computed from posterior draws for $G_s\in\{4,5,6,8,10\}$. \autoref{app:model_fit_selection} reports the comparison. The fitted six-component mixtures closely track the observed consumption distributions and reproduce the location, spread, and skewness of the earnings distributions at the dates shown in \autoref{fig:emp_densities_ot}. The estimated mixture weights vary substantially over time (\autoref{fig:emp_weights}). As in the HANK exercise, the components are approximation devices, so we interpret the density and quantile responses implied jointly by all components rather than any individual weight path.
We combine seven sets of restrictions. They identify four shocks in the VAR block and the common micro shock. The four VAR shocks are supply, demand, monetary policy, and fiscal policy shocks. We jointly identify them using the sign pattern of \citet{baumeister2018inference}, the exogeneity restriction on BZP, and impact sign restrictions. The common micro shock is the latent-factor innovation. Restrictions on its aggregate effects and on its values at selected policy events identify its direction. \autoref{tab:restrictions} summarizes these restrictions, and \autoref{app:technical} provides details.
\begin{table}[htbp]
\centering
\caption{Identifying restrictions in the U.S. application}
\label{tab:restrictions}
\small
\begin{tabular}{@{}p{0.30\textwidth}p{0.46\textwidth}p{0.16\textwidth}@{}}
\toprule
Object restricted & Restriction & Type \\
\midrule
Contemporaneous $(OG,\pi,R)$ block of $\bm A_0$ & Sign and magnitude bounds calibrated using the priors in \citet{baumeister2018inference} (\autoref{tab:Wbounds}) & Sign, magnitude \\
First row of $\bm W$ (BZP equation) & BZP does not respond contemporaneously to the other VAR variables (exogeneity) & Zero \\
Impact responses to fiscal and monetary shocks & The fiscal shock raises the output gap; the contractionary monetary shock raises the federal funds rate and lowers the output gap and inflation & Sign \\
Latent factor $f_t$ at six policy events & The sign of $f_t$ is restricted at the 1993 Omnibus Budget Reconciliation Act, the minimum-wage increases to \$4.75, \$5.15, and \$5.85, the Bush rebates, and the 2003 Jobs and Growth Tax Relief Reconciliation Act & Narrative \\
Aggregate loadings of $f_t$ in $\bm\Lambda_q$ & The BZP loading is set to zero; the output-gap and inflation loadings are restricted to be positive & Zero, sign \\
Macro loadings in the log-weight equations & A higher output gap or inflation shifts mass toward higher-mean components. A higher policy rate shifts mass toward lower-mean components. The coefficient floors are $0.02$ and $0.08$. & Sign, magnitude \\
Factor loadings in the log-weight equations & A positive factor innovation moves mass from both tails toward the middle components. The magnitude floors are $0.05$ for the tails and $0.08$ for the middle. & Sign, magnitude \\
\bottomrule
\end{tabular}
\end{table}
\subsection{Distributional Effects of Monetary and Fiscal Policy Shocks}\label{sub:emp_macro_shocks}
We focus on monetary and fiscal policy shocks. The joint identification also includes supply and demand shocks, but we do not report their responses. {We scale each reported shock to three standard deviations.} For each shock, one figure reports aggregate responses and another reports density and quantile responses; a separate figure collects the inequality responses for both shocks.
\subsubsection{Responses to a Fiscal Policy Shock}
The fiscal-shock identification imposes an increase in the output gap on impact. The response paths at later horizons are learned from the posterior. The median output-gap response is hump-shaped between roughly quarters two and five, while the federal funds rate and inflation also rise (\autoref{fig:emp_fiscal_macro}). The inflation response is less precisely estimated.
\begin{figure}[!htb]
\centering
\includegraphics[width=\textwidth, trim=0 3bp 0 10bp, clip]{results_BZP_zeroF/irf_macro_FiscalBenZeevPappa_size3.pdf}
\caption*{\footnotesize \textbf{Notes}: Aggregate responses (BZP, the output gap, inflation, and the federal funds rate) to a three-standard-deviation fiscal shock. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Macro IRFs to the fiscal shock}
\label{fig:emp_fiscal_macro}
\end{figure}
The posterior median output gap peaks at about $0.3$ percent around quarter three, while inflation rises by about $0.23$ percentage points over the same interval. The federal funds rate response peaks near $0.7$ percentage points around quarter five and then declines slowly. The joint increases in activity, inflation, and the policy rate are consistent with the central bank leaning against a demand expansion.
\begin{figure}[!t]
\centering
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Density IRF: earnings}
\label{fig:emp_fiscal_dens_earn}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_densities_cs1_FiscalBenZeevPappa_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Density IRF: consumption}
\label{fig:emp_fiscal_dens_cons}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_densities_cs2_FiscalBenZeevPappa_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Quantile IRF: earnings}
\label{fig:emp_fiscal_quant_earn}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_quantiles_cs1_FiscalBenZeevPappa_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Quantile IRF: consumption}
\label{fig:emp_fiscal_quant_cons}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_quantiles_cs2_FiscalBenZeevPappa_size3.pdf}
\end{subfigure}
\caption*{\footnotesize \textbf{Notes}: Panels (a) and (b) report the density change $\Delta f(x;h)=f_{\mathrm{shocked}}(x;h)-f_{\mathrm{baseline}}(x;h)$ at selected horizons in response to a three-standard-deviation fiscal shock, and panels (c) and (d) report the 10th, 25th, 50th, 75th, and 90th percentiles of earnings and consumption. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Density and quantile IRFs to the fiscal shock}
\label{fig:emp_fiscal_dist}
\end{figure}
The posterior median density responses in panels (a) and (b) of \autoref{fig:emp_fiscal_dist} trace how the fiscal shock redistributes probability mass. For earnings, mass initially shifts toward below-average values. At longer horizons, it moves from the lower to the upper part of the distribution, consistent with positive responses across the earnings quantiles. For consumption, mass accumulates at low values around quarter four, but this shift reverses at later horizons.
The posterior median earnings-quantile responses in panel (c) of \autoref{fig:emp_fiscal_dist} oscillate at short horizons. They dip around quarter four, coinciding with the peak in the federal funds rate, and are positive at every reported percentile from quarter eight onward. At their peak, they reach about $0.005$ in the IHS-transformed earnings ratio, which near $z_{i,t}=1$ corresponds to roughly a $0.7$ percent increase in earnings relative to two-thirds of nominal GDP per capita. The posterior median consumption-quantile responses in panel (d) fall around quarter five by about two percentage points in the normalized consumption ratio at the median and higher percentiles, then exhibit a small positive hump. The 90\% posterior intervals contain zero at most horizons. Taken together, the posterior medians indicate an eventual rightward shift in the earnings distribution and a temporary consumption decline concentrated in its upper half.
\subsubsection{Responses to a Monetary Policy Shock}
The sign restrictions define a contractionary monetary policy shock as one that raises the federal funds rate and lowers the output gap and inflation on impact. The response paths beyond impact are learned from the posterior. The median policy-rate response rises by about $0.20$ percentage points on impact and {declines substantially within roughly five quarters}. The median output-gap and inflation responses fall by about $0.08$ percent and $0.11$ percentage points, respectively, and {recover most of the decline} within roughly five to seven quarters. All median aggregate responses are close to zero {by around quarters ten to twelve} (\autoref{fig:emp_mp_macro}).
\begin{figure}[!htb]
\centering
\includegraphics[width=\textwidth, trim=0 3bp 0 10bp, clip]{results_BZP_zeroF/irf_macro_MonetaryPolicy_size3.pdf}
\caption*{\footnotesize \textbf{Notes}: Aggregate responses (BZP, the output gap, inflation, and the federal funds rate) to a three-standard-deviation monetary policy shock. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Macro IRFs to the monetary policy shock}
\label{fig:emp_mp_macro}
\end{figure}
The posterior median density responses in panels (a) and (b) of \autoref{fig:emp_mp_dist} trace a short-lived redistribution of probability mass in both distributions. On impact, earnings mass moves from the upper part of the distribution toward its center. By quarter four, the direction reverses, with mass shifting from below-average toward above-average earnings. Consumption mass moves on impact from the middle and upper parts of the distribution toward low values. This shift begins to unwind around quarter four and has largely disappeared by quarter twelve.
\begin{figure}[!t]
\centering
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Density IRF: earnings}
\label{fig:emp_mp_dens_earn}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_densities_cs1_MonetaryPolicy_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Density IRF: consumption}
\label{fig:emp_mp_dens_cons}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_densities_cs2_MonetaryPolicy_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Quantile IRF: earnings}
\label{fig:emp_mp_quant_earn}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_quantiles_cs1_MonetaryPolicy_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Quantile IRF: consumption}
\label{fig:emp_mp_quant_cons}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_quantiles_cs2_MonetaryPolicy_size3.pdf}
\end{subfigure}
\caption*{\footnotesize \textbf{Notes}: Panels (a) and (b) report the density change $\Delta f(x;h)=f_{\mathrm{shocked}}(x;h)-f_{\mathrm{baseline}}(x;h)$ at selected horizons in response to a three-standard-deviation monetary policy shock, and panels (c) and (d) report the 10th, 25th, 50th, 75th, and 90th percentiles of earnings and consumption. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Density and quantile IRFs to the monetary policy shock}
\label{fig:emp_mp_dist}
\end{figure}
At the posterior medians, the monetary tightening compresses the earnings distribution on impact: the tenth percentile rises slightly, while the ninetieth and other upper quantiles fall (panel (c) of \autoref{fig:emp_mp_dist}). All reported consumption quantiles also fall on impact (panel (d)). Across most percentiles, both distributions then exhibit a small positive hump between quarters three and six. The cross-quantile location of the consumption decline provides a diagnostic for the mechanism in \citet{kaplan2018monetary}. If transmission were concentrated among hand-to-mouth households and these households sat mainly in the lower part of the consumption distribution, the decline would concentrate in the lower and middle consumption quantiles, whereas our posterior medians place it at and above the median. The credible sets admit both patterns, so neither transmission mechanism is ruled out. Consumption ranks do not identify liquidity status, so this comparison is informative only through the mapping a given model implies.
\subsubsection{Inequality IRFs}
To summarize the distributional responses, \autoref{fig:emp_inequality_main} reports impulse responses of four inequality measures: the Gini coefficient, the standard deviation, the interquartile range, and the P90--P10 spread, each computed from the fitted mixture on the scale used in estimation (IHS-transformed earnings and normalized consumption). Panels (a) and (b) report earnings and consumption responses to the fiscal shock; panels (c) and (d) report the corresponding responses to the monetary policy shock.
\begin{figure}[!t]
\centering
\begin{subfigure}[b]{\textwidth}\centering
\caption{Fiscal shock: earnings}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_inequality_cs1_FiscalBenZeevPappa_size3.pdf}\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}\centering
\caption{Fiscal shock: consumption}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_inequality_cs2_FiscalBenZeevPappa_size3.pdf}\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}\centering
\caption{Monetary-policy shock: earnings}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_inequality_cs1_MonetaryPolicy_size3.pdf}\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}\centering
\caption{Monetary-policy shock: consumption}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_inequality_cs2_MonetaryPolicy_size3.pdf}\end{subfigure}
\caption*{\footnotesize \textbf{Notes}: Panels (a) and (b) report responses to a three-standard-deviation fiscal shock, panels (c) and (d) responses to a three-standard-deviation monetary policy shock. Each panel reports the responses of the Gini coefficient, standard deviation, interquartile range, and P90--P10 spread (90th minus 10th percentile) of earnings or consumption. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Inequality IRFs to the fiscal and monetary policy shocks}
\label{fig:emp_inequality_main}
\end{figure}
The fiscal shock produces little precisely estimated movement in either earnings or consumption inequality (panels (a) and (b)). The posterior medians of all four earnings measures decline in quarter two, rise in quarters four and five, and fade after about six quarters, mirroring the oscillation in the earnings quantiles. For consumption, the posterior median responses rise around quarter three, briefly reverse, and remain slightly positive through about quarter ten. The measures differ in timing and magnitude. Only the rise in consumption dispersion around quarter three is estimated with some precision; the 90\% credible bands include zero at almost all horizons.
The monetary tightening compresses both distributions on impact (panels (c) and (d)). For earnings, the posterior medians of the standard deviation, interquartile range, and P90--P10 spread fall on impact. The 68\% credible bands of the standard deviation and the P90--P10 spread lie below zero. These responses return to zero within about five quarters. The compression in consumption is larger and more precisely estimated. The P90--P10 spread falls by about $0.8$ percentage points of the consumption ratio, and the 68\% credible bands for all four consumption measures lie below zero on impact. The bands of the Gini coefficient and the standard deviation remain below zero during the first three quarters. These responses return to zero around quarter five.
The inequality measures summarize both impact responses as a compression. The quantile responses in \autoref{fig:emp_mp_dist} show that the incidence differs: the tenth earnings percentile rises while the upper earnings percentiles fall, whereas consumption falls across all reported percentiles, with the largest declines in the upper tail. The summary measures alone do not reveal these cross-quantile patterns, which matter for evaluating the transmission mechanisms of heterogeneous-agent models.
Conditional on our model and priors, the earnings responses place an empirical bound on the quantitative strength of the earnings-heterogeneity channel in \citet{auclert2019monetary}. In that channel, the unequal incidence of labor-income changes causes a contractionary shock to widen the earnings distribution. Our posterior median instead shows a transitory compression in the P90--P10 spread. At each horizon, the upper edge of the 90\% credible band provides a posterior bound on widening.
\subsection{Aggregate Effects of a Common Micro Shock}\label{sub:emp_micro_shock}
The common micro shock is an aggregate innovation to the latent factor $f_t$, which captures movements shared by the earnings and consumption distributions beyond the macro block. One structural interpretation comes from heterogeneous-agent theory. In \citet{bayer2019risk}, a rise in idiosyncratic income risk raises precautionary saving and the demand for liquid assets, causing output and the policy rate to fall. We orient our shock in the reverse direction for dispersion and output: a positive innovation compresses both distributions on impact, lowering realized dispersion, while raising output and inflation.
For models with time-varying idiosyncratic risk, these responses provide two empirical diagnostics:
\begin{enumerate}
\item A dynamic multiplier from dispersion to activity, defined as the peak output gap response divided by the impact response of either the P90--P10 spread or the standard deviation.
\item The policy-rate path. Output and inflation rise on impact by construction, but the sharp impact decline in the federal funds rate is not imposed because its factor loading is unrestricted. The rate path provides a diagnostic for how a model's monetary rule responds to distributional shocks.
\end{enumerate}
\autoref{fig:emp_micro_macro} and \autoref{app:micro_inequality} contain the responses needed to construct these diagnostics; \autoref{fig:emp_micro_dist} reports the underlying density and quantile responses.
\subsubsection{Macro, Density and Quantile Responses}
Beyond the impact restrictions described above, the macro, density, and quantile response paths are learned from the posterior. The factor is excluded contemporaneously from the BZP equation, so the common micro shock does not affect identification of the fiscal shock.
\begin{figure}[!htb]
\centering
\includegraphics[width=\textwidth, trim=0 3bp 0 10bp, clip]{results_BZP_zeroF/irf_macro_Micro_size3.pdf}
\caption*{\footnotesize \textbf{Notes}: Aggregate responses (BZP, the output gap, inflation, and the federal funds rate) to a three-standard-deviation shock to $f_t$. BZP is zero on impact by construction. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Macro IRFs to the micro shock}
\label{fig:emp_micro_macro}
\end{figure}
The posterior median shows a sharp impact decline in the federal funds rate (\autoref{fig:emp_micro_macro}). The aggregate responses then return toward zero and are generally close to zero after eight to ten quarters.
\begin{figure}[!t]
\centering
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Density IRF: earnings}
\label{fig:emp_micro_dens_earn}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_densities_cs1_Micro_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Density IRF: consumption}
\label{fig:emp_micro_dens_cons}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_densities_cs2_Micro_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Quantile IRF: earnings}
\label{fig:emp_micro_quant_earn}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_quantiles_cs1_Micro_size3.pdf}
\end{subfigure}\\[0.2em]
\begin{subfigure}[b]{\textwidth}
\centering
\caption{Quantile IRF: consumption}
\label{fig:emp_micro_quant_cons}
\includegraphics[width=\textwidth, trim=0 0 0 10bp, clip]{results_BZP_zeroF/irfs_quantiles_cs2_Micro_size3.pdf}
\end{subfigure}
\caption*{\footnotesize \textbf{Notes}: Panels (a) and (b) report density changes $\Delta f(x;h)$ at selected horizons in response to a three-standard-deviation shock to $f_t$, and panels (c) and (d) report the 10th, 25th, 50th, 75th, and 90th percentiles of earnings and consumption. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Density and quantile IRFs to the micro shock}
\label{fig:emp_micro_dist}
\end{figure}
The density responses in panels (a) and (b) of \autoref{fig:emp_micro_dist} show the impact compression in both distributions as probability mass moves from the tails toward the center. At the posterior medians, the tenth earnings percentile rises on impact, while the median and upper percentiles fall (panel (c)). All reported consumption percentiles fall on impact, with the largest declines in the upper tail (panel (d)). Declines at the median and upper percentiles are relatively precise at short horizons, but posterior uncertainty widens rapidly thereafter.
The inequality responses in \autoref{fig:emp_micro_inequality} summarize the same impact compression. For both earnings and consumption, all four measures fall sharply on impact, as implied by the identifying restrictions on the factor loadings, and then rebound before fading. Posterior uncertainty widens rapidly after impact. Because the initial compression is imposed, the strength and persistence of the rebound are the informative features.
\section{Benchmarking an Estimated HANK Model}\label{sec:hankbench}
We benchmark the JAMM-VAR posterior responses to a monetary policy shock against the estimated HANK model of \citet{bayer2024shocks}. We take the model directly from the authors' public repository and do not alter its economic structure.\footnote{Code and parameter values are those distributed at \url{https://github.com/BASEforHANK/HANK_BusinessCycleAndInequality}. We re-estimate nothing and choose no parameter ourselves. The repository ships a solved steady state and the linear system that governs the response to each shock, so we load the model rather than re-solve it.} Two features make the model well suited to this exercise: monetary policy affects its earnings distribution, and its shock processes were estimated using data that include measures of inequality. Following \autoref{sec:hank_targets}, we compare the model-implied responses with our posterior bands horizon by horizon.
The empirical analysis in \autoref{sec:empirical} implements the first four steps of the protocol. To begin the fifth step, we scale the monetary policy shock in the HANK model so that the policy rate rises on impact by the JAMM-VAR posterior median of about $0.20$ percentage points.
\subsection{Making the HANK Objects Comparable}\label{sec:bblstep5}
Households in the model face idiosyncratic productivity risk, hold liquid and illiquid assets, and pay a cost to move funds between them. Prices and wages are sticky, and a small group of entrepreneurs receives profit income rather than labor income.
Having matched the shock size, we align the remaining model and empirical objects in their variable definitions, transformations, measurement conventions, and response horizons. We construct model earnings $z_{i,t}$ for each productivity state from the model's income function: a productivity-dependent share of the wage bill plus an equal allocation of union rents. We then apply the same inverse hyperbolic sine transformation used for the survey data, $x_{i,t}=\operatorname{asinh}(z_{i,t})$. We exclude the entrepreneurial state because its income consists of profits rather than labor earnings. Finally, we convert the model's annualized quarterly inflation rate into the four-quarter measure used in the empirical analysis by averaging it over the current and previous three quarters.
The nonlinear inverse hyperbolic sine transformation makes the level of $z_{i,t}$ relevant for comparing responses, so we make one level adjustment. The empirical mean of $z_{i,t}$ is $1.328$, whereas the model mean is $0.934$. We therefore multiply model earnings by their ratio, $1.328/0.934=1.42$. This normalization aligns the steady-state means rather than fitting the model's responses.\footnote{Because the inverse hyperbolic sine transformation is nonlinear, the level at which it is applied affects the scale of the transformed responses.}
Before comparing responses, we check whether the model's steady-state earnings distribution resembles the empirical cross section. In units of $x_{i,t}$, the tenth and ninetieth percentiles are $0.35$ and $1.68$ in the survey and $0.54$ and $1.65$ in the model. The model closely matches the ninetieth percentile but places the tenth percentile too high. Its P90--P10 range is therefore $83$ percent of the empirical range. The two distributions are broadly comparable, with the main discrepancy concentrated in the lower tail.
Two limitations restrict the comparison. First, the model has no unemployment, so it cannot reproduce movements into and out of employment. Because the CPS earnings cross sections cover employed individuals, the comparison concerns the intensive margin among the employed on both sides. Empirical earnings responses may nevertheless reflect changes in the composition of the employed, which the model does not capture. Second, the public code does not expose the joint distribution of income and wealth needed to recover the model's consumption distribution. We therefore restrict the comparison to earnings.
\subsection{Diagnosing the Model}\label{sec:bblstep6}
The sixth step evaluates the HANK model against the JAMM-VAR posterior of the response objects rather than only their posterior medians. We use the credible bands to assess whether differences in point estimates are large relative to posterior uncertainty. \autoref{fig:bblmacro} compares the aggregate and inequality responses. \autoref{fig:bblfan} compares the earnings-quantile responses on impact and the P90--P10 response across horizons, while \autoref{fig:bbldyn} traces each reported earnings-quantile response over twelve quarters.
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{figs/bbl_bh_macro.pdf}
\caption*{\footnotesize \textbf{Notes}: Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands. The HANK model is in dark green, scaled so that the policy rate rises $0.197$ percentage points on impact. The last two panels report the two inequality measures the model itself reports, computed on our side from the fitted mixture. Both the empirical and the model measures are computed on the IHS-transformed earnings scale, $x_{i,t}=\operatorname{asinh}(z_{i,t})$. In the last panel the green line is the top 10\% earnings share, which is the object the survey measures, and the amber line is the share of total income including income from capital.}
\caption{Aggregate and inequality IRFs to the monetary policy shock}
\label{fig:bblmacro}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{figs/bbl_bh_fan.pdf}
\caption*{\footnotesize \textbf{Notes}: Earnings quantile responses on impact, and the P90--P10 spread over quarters. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Earnings quantile IRFs on impact and the P90--P10 spread}
\label{fig:bblfan}
\end{figure}
\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth]{figs/bbl_bh_quantdyn.pdf}
\caption*{\footnotesize \textbf{Notes}: One panel per reported percentile, over twelve quarters. Navy lines are posterior medians, and shaded areas are 68\% and 90\% credible bands.}
\caption{Earnings quantile IRFs by horizon}
\label{fig:bbldyn}
\end{figure}
On impact, the HANK model reproduces the cross-quantile pattern in the JAMM-VAR earnings responses (\autoref{fig:bblfan}). It raises earnings at lower quantiles and lowers them at upper quantiles, matching the ordering of the posterior medians. The model-implied responses also match the magnitudes. At every reported percentile, the model lies within $0.72$ posterior standard deviations of the posterior median. The model-implied P90--P10 impact response equals $1.08$ times the posterior median. The impact responses are therefore not a dimension along which the model and the evidence disagree.
The model generates this cross-quantile pattern through wage-setting rents. Because wages adjust slowly, a contraction widens the gap between what firms pay and what workers would accept. The resulting rents rise $4.3$ percent on impact and reach $7.0$ percent in the second quarter. They are distributed equally across households rather than in proportion to earnings, so this flat transfer lifts lower earnings quantiles relative to upper quantiles. Thus, the unequal incidence of monetary policy across workers stems from wage setting rather than labor supply. The comparison links the empirical distributional pattern to a specific model mechanism that aggregate responses alone cannot reveal.
The HANK model also reproduces the impact compression in two summary measures of earnings inequality (last two panels of \autoref{fig:bblmacro}). Neither response was directly targeted in the model's estimation or calibration. The model's standard deviation of IHS-transformed earnings falls $0.57$ percent on impact, close to the JAMM-VAR posterior median of $0.48$ percent and within the 68\% credible band of $0.25$ to $0.95$ percent. The top 10\% earnings share falls $0.26$ percent in the model, compared with a posterior median of $0.59$ percent. Both measures therefore imply the same transitory compression reported in \autoref{sub:emp_macro_shocks}, though the model overstates the fall in dispersion and understates the fall in the top share. By contrast, the model's top 10\% share of total income, which includes capital income, rises $0.12$ percent. The model's labor-earnings distribution is therefore broadly consistent with the evidence, but its capital-income channel reverses this response when total income is considered.
The main disagreement concerns strength and persistence: the HANK model's contraction is too large on impact and too short-lived (\autoref{fig:bblmacro}). Output falls about $0.20$ percent on impact, compared with a JAMM-VAR posterior median of $0.08$ percent, while inflation falls $0.25$ percentage points, compared with a posterior median of $0.11$ percentage points. The model-implied aggregate responses leave the credible bands at several horizons, most clearly for output. The policy-rate response reverses within two quarters, so the model's contraction is largely over by the fifth quarter, when the JAMM-VAR posterior medians remain near their troughs. The same lack of persistence appears in the earnings distribution (\autoref{fig:bbldyn}). The JAMM-VAR posterior places the largest earnings responses between the fourth and sixth quarters, whereas the model's quantiles return to their starting values by the fourth quarter because its wage-setting rents dissipate once the policy-rate response reverses.
Since the model does not reproduce the persistence of the aggregate responses, the comparison does not show whether the earnings response would remain too short-lived in a version of the model with more persistent aggregate dynamics. \citet{auclert2020microjumps} show that sticky household expectations can generate hump-shaped aggregate responses while preserving the front-loaded response of consumption to transitory household income shocks.
\subsection{Implications of the Comparison}\label{sec:bblneeds}
The comparison points to one requirement for the model's monetary transmission mechanism and one choice about how to map empirical earnings observables to model-based objects.
\begin{itemize}
\item The model needs weaker but more persistent monetary transmission. One possibility is a policy rule that generates a more persistent policy-rate response. This would extend both the aggregate contraction and the cross-sectional earnings response, since both currently end when the policy rate reverses. Matching their magnitudes is harder. The degree to which wages lag prices governs both the strength of aggregate transmission and the size of the earnings-distribution response. Weakening this channel to improve the aggregate fit would therefore also attenuate the cross-sectional response. The model needs an additional channel that affects workers unequally without amplifying aggregate transmission.
\item Because the model abstracts from unemployment, the comparison concerns the intensive margin among employed workers. The CPS earnings cross sections cover employed individuals, so the empirical statistics are already conditional on employment. Compositional changes in the pool of the employed can still move the empirical distribution without a model counterpart, and we interpret discrepancies in light of this difference.
\end{itemize}
\section{Extensions}\label{sec:extensions}
We now highlight two extensions that can be added in a straightforward manner to our model, namely cross sections observed at only some dates and stochastic volatility. Each adds or replaces one block in the sampler of \autoref{sub:posterior} and requires only limited adjustments to the remaining steps. This section states what each extension adds. \autoref{app:extensions} provides the modified samplers.
\subsection{Cross Sections Observed in Some Periods Only}\label{sub:ext_missing}
Cross-sectional surveys often start later than aggregate series, contain gaps, or are fielded less frequently. This mismatch binds in our application: the consumption cross sections begin in 1990, although the aggregate series extend several decades further back. The extension lets the macro block use those earlier aggregate observations while treating the corresponding cross sections as missing. It also accommodates lower-frequency surveys, such as annual income data paired with quarterly national accounts or triennial Survey of Consumer Finances (SCF) data paired with monthly aggregate variables.
Conditional on the macro state, common factors, and parameters, the mixture model defines a population density for each cross section at every date, whether or not the cross section is observed. Missing micro data matter only when their sample quantiles enter an in-sample macro equation through the feedback term in \autoref{eq: VAR}. We therefore augment the posterior with the missing cross sections needed to construct those quantiles and update them in a Metropolis--Hastings imputation block \citep{tanner1987}. Conditional on the completed cross sections, the existing blocks retain their form. Missing cross sections whose quantiles do not enter an in-sample macro equation are integrated out. If all quantile-feedback coefficients are set to zero, no imputation is needed, and the micro likelihood runs only over observed dates.
The posterior therefore yields model-implied cross-sectional densities, with credible bands, at the frequency of the aggregate block, including dates between survey waves. At dates without micro data, the aggregate likelihood informs the common factors through their loadings in the macro equations; the macro state and factors then determine the mixture weights. For example, pairing the triennial SCF with monthly aggregate variables produces a monthly path for the wealth distribution between SCF waves. These paths are model-based interpolations rather than direct survey measurements.
One modeling choice matters. When a survey was fielded but its micro data are unavailable, the latent cross section inherits the survey's design sample size. When no survey was fielded, the latent sample size is a modeling choice. It controls the sampling noise in the latent quantiles and therefore the model itself.
\subsection{Stochastic Volatility}\label{sub:ext_sv}
Structural shock variances change over long macroeconomic samples. The Great Moderation is the leading example. Such variation may matter for the normalization of shocks and the interpretation of magnitude restrictions. Because the structural shocks are orthogonal and $\bm D$ is diagonal, volatility enters equation by equation. Replace each constant variance $d_i$ with $d_{i,t}=\exp(h_{i,t})$ and let the log volatility follow an autoregression. One per-equation block, built on the mixture sampler of \citet{kim1998stochastic}, draws the volatility path and its parameters. This block replaces the structural-variance step, exactly as in VARs with time-varying volatility \citep{primiceri2005time}. Every likelihood term that previously depended on $d_i$ now depends on $d_{i,t}$, so the coefficient conditionals take generalized-least-squares form. The coefficient prior must also be modified to accommodate time-varying volatility. \autoref{app:ext_sv} provides the details.
\section{Conclusions}\label{sec:conclusions}
Aggregate SVARs provide empirical benchmarks for structural macroeconomic models. Our JAMM-VAR extends this role to HANK models by adding several marginal distributions to a Bayesian SVAR. The model retains standard aggregate-shock identification, separately identifies a common micro shock, and estimates aggregate and distributional responses within one posterior. The household observations are repeated cross sections, possibly from separate surveys. Household histories and first-stage density estimates are unnecessary. Our approach jointly estimates all relevant objects, taking into account the estimation uncertainty in the aggregate and distributional blocks jointly, while using the Bayesian and VAR toolkits macroeconomists are familiar with.
The main payoff of our approach is a direct, posterior-based comparison between estimated aggregate and distributional responses and their counterparts in heterogeneous-agent equilibrium models. The HANK simulation shows that this comparison remains informative even when the JAMM-VAR does not nest the equilibrium model. Conditional on matching the policy variable's impact response, the JAMM-VAR recovers the broad propagation of a monetary policy shock across aggregates, earnings quantiles, and consumption quantiles, and the posterior bands generally contain the true paths. In the U.S. application, aggregate data, CPS earnings, and CEX consumption yield a joint posterior of response objects that we compare with HANK-model responses horizon by horizon. This comparison distinguishes features the HANK model reproduces, such as the impact cross-quantile pattern, from dimensions along which its propagation differs, especially strength and persistence. More generally, our approach lets researchers use repeated cross sections to assess and refine heterogeneous-agent models against a coherent joint posterior of response objects rather than isolated point estimates.
\newpage
{\setstretch{1.15}\bibliographystyle{custom}\bibliography{references}}
\clearpage