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.
105,832 characters
Empirical Bayes Estimation in Heterogeneous Coefficient Panel Models
\title{Empirical Bayes Estimation in Heterogeneous Coefficient Panel Models}
\author{\textsc{Myunghyun Song}\thanks{Department of Economics, Columbia University} \and \textsc{Sokbae Lee}\thanks{Department of Economics, Columbia University and Institute for Fiscal Studies} \and \textsc{Serena Ng}\thanks{Department of Economics, Columbia University and NBER\newline
We thank Jiaying Gu, Roger Koenker, Soonwoo Kwon, Bodhisattva Sen, Kaizheng Wang, and
seminar participants at Indiana University, the University of Pennsylvania, the University of Toronto, and Yale University for helpful and encouraging comments.
The third author would like to thank the National Science Foundation for financial support (SES: 2018369).}}
\date{February 6. 2026}
\maketitle
\thispagestyle{empty}
\setcounter{page}{0}
\begin{abstract}
We develop an empirical Bayes (EB) G-modeling framework for short-panel linear models with nonparametric prior for the random intercepts, slopes, dynamics, and non-spherical error variances. We establish identification and consistency of the nonparametric maximum likelihood estimator (NPMLE) under general conditions, and provide low-level sufficient conditions for several models of empirical interest. Conditions for regret consistency of the EB estimators are also established. The NPMLE is computed using a Wasserstein-Fisher-Rao gradient flow algorithm adapted to panel regressions. Using data from the Panel Study of Income Dynamics, we find that the slope coefficient for potential experience is substantially heterogeneous and negatively correlated with the random intercept, and that error variances and autoregressive coefficients vary significantly across individuals. The EB estimates reduce mean squared prediction errors relative to individual maximum likelihood estimates.
\\
\\
\noindent \textsc{Keywords}: G-modeling, nonparametric maximum likelihood, shrinkage estimation, Wasserstein-Fisher-Rao gradient flow, income dynamics
\end{abstract}
\newpage
\section{Introduction}
Understanding differences in individual behavior is the primary goal of many economic analyses. Why do workers with similar experience exhibit different earnings profiles? Why do economic agents with similar observed characteristics respond differently to economic shocks and policy interventions? A challenge for empirical researchers in addressing these questions is the presence of individual-level differences that affect behavior but are not captured by observed variables, commonly referred to as unobserved heterogeneity. While it is not hard to allow for unit-specific intercepts in linear regressions, entertaining heterogeneous slope parameters and innovation variances is more challenging, and in a frequentist setting, this usually requires parametric assumptions to gain tractability at the cost of robustness.
This paper considers a setting in which we observe $\{ (Y_{it}, X_{2,it}) : t = 1,\ldots,T \}$ for each unit $i = 1,\ldots,N$, assuming they are independent and identically distributed (i.i.d.) across $i$. We propose a framework for short-$T$ linear panel models with heterogeneous intercepts, slope coefficients, and error distributions. To fix ideas, consider the dynamic linear panel model with an exogenous covariate $X_{2,it}$ and AR(1) errors:
\begin{subequations}\label{def:model:hivxd}
\begin{align}
Y_{it} &= a_i + b_i X_{2,it} + u_{it}, \\
u_{it} &= \rho_i u_{it-1} + \sigma_i e_{it},
\quad e_{it} \overset{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0,1), \qquad t \ge 1,
\label{def:model:hivxd:2}
\end{align}
\end{subequations}
where
$a_i \in \mathbb{R}$ is a random intercept,
$b_i \in \mathbb{R}$ a random slope,
$\rho_i \in (-1,1)$ a random autoregressive coefficient,
$\sigma_i^2 > 0$ a random variance,
and $u_{i0}$ follows the stationary distribution $\mathcal{N}(0,\sigma_i^2/(1-\rho_i^2))$ conditional on $(\sigma_i^2,\rho_i)$. We refer to the model in \eqref{def:model:hivxd} as HIVDX since it admits heterogeneity in four dimensions: intercepts ($I$), error variances ($V$), dynamics ($D$) through the autoregressive coefficient, and the coefficient on the exogenous covariate $X_2$. The widely studied Gaussian location-scale model $Y_{it} = a_i + \sigma_i e_{it}$,
as well as models that impose $\rho_i = \rho$ and $b_i = b$, are special cases of HIVDX.
We analyze the HIVDX model and the special cases in Section~\ref{sec:model} using a general framework that partitions the random parameters into two sets: regression coefficients $\beta_i = (a_i, b_i)$ of dimension $d_\beta = 2$, and covariance parameters $\delta_i = (\sigma_i^2, \rho_i)$ of dimension $d_\delta = 2$. The parameters can be estimated separately for each $i$ by maximum likelihood. Instead, we assume that the random parameters $\theta_i \equiv (\beta_i, \delta_i)$ are drawn from a true but unknown common distribution $G_*$. Viewing $G_*$ as a prior, the optimal decision rule for each unit under compound squared-error loss is the posterior mean
${\theta}^{*}_i := \mathbb{E}_{G_*}[\theta_i \mid Y_i, X_{2,i}]$,
where $Y_i = (Y_{i1}, \ldots, Y_{iT})'$, $X_{2,i} = (X_{2,i1}, \ldots, X_{2,iT})'$, and $\mathbb{E}_{G_*}[\cdot]$ denotes expectation with respect to $\theta_i \sim G_*$.
Specifying that $\theta_i \stackrel{\text{i.i.d.}}{\sim} G_*$ is grounded in the principle of empirical Bayes (EB) introduced by \citet{robbins:56}. This approach is intrinsically frequentist because it estimates the prior distribution from the data. EB estimation improves on individual-level maximum likelihood estimation through shrinkage, borrowing strength across units to learn both the individual parameters $\theta_i$ and their population distribution $G_*$.
Following the terminology of \citet{efron-14}, the procedure is referred to as G-modeling when $G_*$ is estimated directly.
The resulting EB estimator is given by $\hat{\theta}^{\mathrm{EB}}_i := \mathbb{E}_{\hat{G}}[\theta_i \,|\, Y_i, X_{2,i}]$, where $\hat{G}$ denotes an estimator of $G_*$. An EB estimator $\hat \tau_i^{{\mathrm{EB}}}(\theta_i)=\mathbb E_{\hat G}[\tau_i \,|\, Y_i,X_{2,i}] $ for general unit-specific parameters $\tau_i$ defined as functions of $\theta_i$ can also be defined. While parametric G-modeling is convenient, nonparametric modeling guards against misspecification.
Recent reviews of EB and its applications in economics include \citet{KoenkerGu2024:JPEmicro}, \citet{KoenkerGu2025book}, and \citet{Walters:24}. To our knowledge, no nonparametric EB analysis has considered all four sources of heterogeneity present in the HIVDX model.
We study identification and estimation of the HIVDX model within an EB framework. The generality of our model with $d_\beta > 1$ and $d_\delta > 0$ introduces several challenges. First, while assuming the identification of $G_*$ in the Gaussian location model with i.i.d.\ errors without formal proof may be defensible, this assumption cannot be taken for granted in more complex settings.
Theorem~\ref{thm:identification} establishes that identification requires: (i) the design matrix to have full column rank almost surely, and (ii) the covariance structure to remain identifiable after removing the mean component by appropriately differencing the data.
These conditions imply that $d_\beta < T$, leaving $T - d_\beta$ degrees of freedom to identify the covariance structure; this imposes a constraint on the dimension of heterogeneous variance parameters ($d_\delta$).
We apply Theorem~\ref{thm:identification} to the location model ($d_\beta=1$) and illustrate how differencing facilitates identification. Specifically, we show that $T \ge 3$ suffices for AR(1) errors (Proposition~\ref{prop:dyp-identification}), whereas ARMA(1,1) errors require $T \ge 4$ (Proposition~\ref{prop:identification-ARMA}). The analysis is, however, more involved when $d_\beta > 1$, as the mean and variance components are convolved within the outcome distribution. Now whether and how $G_*$ can be identified depends critically on the specific time-series variation of the covariate.
We examine two canonical examples.
The first example includes models with a common time trend and heterogeneous initial conditions (e.g., potential experience in income dynamics) as a special case. For this example, Proposition \ref{prop:HIVDX-id-T=4} establishes that $G_*$ is nonparametrically identified with $T \ge 4$.
The second example includes a common mean shift as a special case. For this case, Proposition \ref{prop:HIVDX-id-T=5} shows that identification of $G_*$ is possible with $T\ge 5$. Identification must be worked out on a case by case basis.
The second challenge lies in establishing consistent estimation of the HIVDX model.
The classic proof of \citet{kiefer1956consistency} assumes compactness, but when $\beta_i$ is allowed to be unbounded, the space of mixing distributions is non-compact.
While the case of ${d_{\beta}} = 1$ can be resolved by embedding $\mathbb{R}$ into its compact completion $[-\infty, \infty]$, the multidimensional setting (${d_{\beta}} \ge 2$) presents unique difficulties because mass can diverge along infinitely many distinct directions.
To address this, Section~\ref{sec:consistency} adopts a metric induced by the vague topology to effectively compactify the space of prior distributions.
Theorem~\ref{thm:consistency} establishes almost sure consistency of the
nonparametric maximum likelihood estimator (\mbox{NP\hspace{0.05em}MLE}) with respect to this metric.
Building on this result, Theorem~\ref{thm:EB-consistency} in Section~\ref{subsec:EB-consistency} provides conditions for regret consistency (equivalently, asymptotic optimality) of the EB estimator for $\tau_i$, taking into account that $\tau_i$ may involve functions of the random coefficients and the data. For example, the optimal (but infeasible) one-step-ahead prediction under quadratic loss is ${\tau}^{*}_i = \mathbb{E}_{G_*}[a_i + b_i X_{2,i T+1} + \rho_i Y_{iT} \,|\, Y_{i}, X_{i},X_{2,i T+1}]
- \mathbb{E}_{G_*}[\rho_i a_i +\rho_i b_i X_{2,iT} \,|\, Y_{i}, X_{i}]$. The prediction error $\hat \tau^{{\mathrm{EB}}}_i-\tau_i^*$ thus involves interaction of the errors in estimating the random coefficients and $(Y_i,X_i)$. Proposition~\ref{prop:EB-HIVDX} provides low-level sufficient conditions for the EB prediction to be regret-consistent.
The third challenge arises from the computational complexity associated with multiple sources of heterogeneity.
In location models, the NPMLE of $G$ is typically computed by optimizing over a fixed grid of $m$ points.
The location of the grid points and the size of the grid are usually fixed during the iterations.
However, as $d_\theta$ increases, the grid required to maintain accuracy needs to grow exponentially.
It would be desirable to have an algorithm that allows the support points to move adaptively.
This motivates a shift from Euclidean optimization to optimization on the space of probability measures, where the updates are governed by gradient flows.
In this framework, the choice of geometry is critical.
For example, the Fisher-Rao gradient flow corresponds to reweighting fixed support points, whereas the Wasserstein gradient flow transports mass across the parameter space, effectively moving the support.
Recently, \citet{yan2024learning} demonstrate the theoretical and computational advantages of combining these geometries for the Gaussian location mixture model.
Motivated by this encouraging result, Section~\ref{sec:implementation} develops a Wasserstein-Fisher-Rao (WFR) gradient flow algorithm for panel regressions with multiple sources of heterogeneity.
We use the EB approach to analyze the income data in the Panel Study of Income Dynamics (PSID) in Section~\ref{sec:empirical}.
The model allows earnings $Y_{it}$ to have unit-specific mean $a_i$ and response $b_i$ to potential experience $X_{2,it}$, with a heteroskedastic error covariance structure that depends on internal dynamics $\rho_i$ and exposure to shocks through $\sigma_i^2$. We find that (i) the slope coefficient $b_i$ for potential experience exhibits substantial heterogeneity and is negatively correlated with the random intercept $a_i$, and (ii) there is pronounced cross-sectional heterogeneity in both the variance ($\sigma_i^2$) and the dynamics ($\rho_i$).
The EB estimates reduce mean squared prediction errors relative to individual maximum likelihood estimates.
Section~\ref{sec:Monte:Carlo} shows that features of the EB estimates and forecast error reduction are replicated in Monte Carlo experiments calibrated to the application.
Before turning to the main part of the paper, we briefly discuss recent work that, while not directly related to our approach, are conceptually connected to the broader EB literature.
\citet{KevinChen:2024} analyzes conditionally Gaussian settings with known variances \(\sigma_i^2\). He first studentized the data by estimating the conditional mean and variance of outcomes given \(\sigma_i^2\), and then estimate the prior to plug into EB decision rules.
The paper emphasizes the importance of nonparametric estimation in both stages.
In contrast, we allow for unobserved \(\sigma_i^2\) and heterogeneous dynamics across multiple time periods.
\citet{Kwon:2025} considers a multivariate Gaussian model with a hierarchical prior, while \citet{Cheng-et-al:2025} propose a parametric EB estimator for two-way effects with assortative matching.
Both methods rely on linear shrinkage rules selected to minimize unbiased risk estimates.
In contrast, we adopt a fully nonparametric EB approach via \mbox{NP\hspace{0.05em}MLE}.
\citet{Gaillac2024} develops posterior mean estimators for multiple random coefficient models for cross-sectional regressions.
In a panel extension, the paper constructs time averages for each unit and allows heterogeneity only in time-invariant regressors, while assuming homogeneous coefficients for time-varying regressors.
In contrast, our framework focuses on random coefficients associated with regressors that vary over time.
\citet{AGT:EBAdaptive} investigate EB methods in compound adaptive experiments, where the outcome of each arm follows a normal distribution with an unknown mean. They demonstrate that the risk guarantees for G-modeling, originally derived under i.i.d.\ sampling, continue to hold for adaptively collected data without knowledge of the sampling algorithm. In contrast, F-modeling is shown to yield biased estimates in this setting.
Recent work has extended nonparametric EB methods beyond the classical Gaussian location model to a variety of settings, including heteroskedastic models (e.g., \citeauthor{{Bodhi:JRSSB}}, \citeyear{{Bodhi:JRSSB}}),
high-dimensional or multivariate regression (e.g., \citeauthor{{Bodhi:EB-HD}}, \citeyear{{Bodhi:EB-HD}}; \citeauthor{{Wu:EB-HD}}, \citeyear{{Wu:EB-HD}}; \citeauthor{{Jiang:2025}}, \citeyear{{Jiang:2025}}), and variance estimation (e.g., \citeauthor{{IgnatiadisSen2025AoS}}, \citeyear{{IgnatiadisSen2025AoS}}). These works only consider a static model while our focus is the heterogeneous dynamic panel model presented above as HIVDX.
\section{The Econometric Setup and EB Modeling}\label{sec:model}
We consider a panel of observations on $N$ cross-sectional units over $T$ time periods, focusing on a short-panel framework in which $N$ grows while $T$ remains fixed.
Our main setup is a general heterogeneous coefficient (HC) model that nests the HIVDX model:
\begin{equation}
\label{eq:model}
Y_i = X_i \beta_i + P_i e_i ,
\end{equation}
where \(Y_i = (Y_{i1}, \ldots, Y_{iT})'\) is a \(T \times 1\) vector of outcomes, \(X_i = (X_{i1}, \ldots, X_{iT})'\) is a \(T \times d_{\beta}\) matrix of observed covariates, and $e_i = (e_{i1}, \ldots, e_{iT})'$ is a vector of independent Gaussian errors.
We assume that the first column of $X_i$ is $\iota_T$, the $T$-dimensional vector of ones.
For each $i$, $P_i$ is lower triangular, and the corresponding conditional covariance matrix $\Sigma_i \equiv P_i P_i'$ is positive definite. To capture its dependence on heterogeneous variance parameters $\delta_i \in \mathbb{R}^{{d_{\delta}}}$, we write $P_i = P(\delta_i)$, where $\delta \mapsto P(\delta)$ is a continuously differentiable function of $\delta$.
To see that the HIVDX model is a special case of the HC model for particular $(\beta_i, \delta_i)$, note that since $e_{it}$ are i.i.d. $\mathcal{N}(0,1)$, $u_i = (u_{i1}, \ldots, u_{iT})'$ is normally distributed with conditional covariance matrix $\Sigma_i$, whose $(t,s)$ entry is given by
\begin{equation}\label{def:Sigma:model}
\Sigma_{i,ts} = \sigma_i^2 \frac{\rho_i^{|t-s|}}{1-\rho_i^2}, \quad t, s = 1, \ldots, T.
\end{equation}
Hence, $P_i$ is defined via the Cholesky decomposition of $\Sigma_i = P_i P_i'$ with $\delta_i = (\sigma_i^2, \rho_i)$.
We develop our theoretical framework through the lens of the HC model since it provides a unified way to analyze a broader class of dynamics such as ARMA errors.
Let $I_T$ denote the $T \times T$ identity matrix. The following assumptions will be used throughout.
\begin{assumption}[DGP]\label{asm:DGP}
Let \((Y_i, X_i)_{i=1}^N\) be identically and independently distributed (i.i.d.) data, satisfying \eqref{eq:model}. In addition, the following conditions hold.
\begin{enumerate}
\item[(i)]
The pairs \((\beta_i, \delta_i)_{i=1}^N\) are i.i.d. draws from an unknown distribution \(G_*\) and are independent of \((X_i)_{i=1}^N\).
\item[(ii)]
The parameter space for $\delta_i$ is a compact subset $\mathcal{K}_{\delta} \subseteq \mathbb{R}^{{d_{\delta}}}$.
For each $\delta_i \in \mathcal{K}_{\delta}$, the covariance matrix $\Sigma_i \equiv P(\delta_i)P(\delta_i)'$ satisfies $\underaccent{\bar}{c}\, I_T \le \Sigma_i \le \bar{c}\, I_T$ for some $0 < \underaccent{\bar}{c} < \bar{c} < \infty$.
\item[(iii)]
$e_i \sim \mathcal{N}(0, I_T)$.
\end{enumerate}
\end{assumption}
Conditional on \((\beta_i, \delta_i)\) and \(X_i\), the outcome vector \(Y_i\) follows the distribution \(\mathcal{N}(X_i \beta_i, \Sigma_i)\). Individuals in this model are distinguished by types characterized by their heterogeneous features. Given the information structure, individual types cannot be directly uncovered from the data, but their distributions can be identified.
Let $\theta_i := (\beta_i, \delta_i)$ collect all heterogeneous parameters, lying in the space $\Theta := \mathbb{R}^{{d_{\beta}}}\times\mathcal{K}_{\delta} \subseteq \mathbb{R}^{{d_{\theta}}}$ with ${d_{\theta}} = {d_{\beta}} + {d_{\delta}}$.
If $\theta_i$ were observed, the conditional density of $Y_i$ given $X_i$ is
\begin{align*}
\ell(Y_{i} \,|\, X_i, \theta_i) := \frac{1}{({2\pi})^{T/2}|P(\delta_i)|} \exp \left( - \frac{1}{2}\|P(\delta_i)^{-1} (Y_i - X_i\beta_i )\|^2 \right).
\end{align*}
When $\theta_i$ is unobserved but its distribution is given as $G$, the marginal density of $Y_i$ given $X_i$, obtained by integrating out $\theta_i$ from the likelihood, is
\begin{equation*}
f_{G}(Y_i, X_i) := \int_{\Theta} \ell(Y_i \,|\, X_i, \theta) dG(\theta),
\end{equation*}
where the independence between $\theta_i$ and $X_i$ is used.\footnote{We use $G_*$ to denote the unknown true distribution and $G$ to denote a generic distribution in the class of possible distributions that includes $G_*$. The notation $f_{G}(Y_i, X_i)$ emphasizes that the marginal density of $Y_i$ given $X_i$ depends on $G$.}
For the estimation of $(\theta_i)_{i=1}^N$, two conceptual frameworks are available.
The first is the \emph{separate} decision framework that treats each $\theta_i$ as fixed and minimizes the individual risk $\mathbb{E}[\|\hat \theta_i - \theta_i\|^ 2 \,|\, \theta_i]$.
This yields the individual MLE,
\begin{equation*}
\hat {\theta}^{\operatorname{MLE}}_i := \operatornamewithlimits{arg\hspace{0.1em} max}_{\theta \in \Theta} \ell(Y_i \,|\, X_i, \theta),
\end{equation*}
which provides an asymptotically efficient estimator for $\theta_i$ as $T \to \infty$.
However, \citet{stein:54} and \citet{james-stein:61} demonstrate
for independent Gaussian observations
that the MLE is inadmissible under squared-error loss when the dimension of the parameter vector is at least three (here, $N \ge 3$). In this setting, the MLE is dominated by a shrinkage estimator that achieves lower total risk by pooling information across observations.
The second approach is the \emph{compound} decision framework. Unlike the separate framework, which minimizes risk for each unit individually, this approach seeks to minimize the aggregate risk over the ensemble of $N$ units.
While the strict compound decision formulation treats $(\theta_i)_{i=1}^N$ as fixed, \citet{robbins:56} demonstrates that this problem can be effectively solved by adopting an EB perspective.
By modeling the parameters as i.i.d.\ draws from a common latent distribution $G_*$, we can derive an optimal (oracle) decision rule, which takes the form of the posterior mean:
\begin{equation}
\label{eq:posterior-mean}
{\theta}^{*}_i := \mathbb{E}_{G_*}[\theta_i \,|\, Y_i, X_i] = \frac{\int_{\Theta} \theta \ell(Y_i \,|\, X_i, \theta) dG_*(\theta)}{\int_{\Theta} \ell(Y_i \,|\, X_i, \theta) dG_*(\theta)}.
\end{equation}
This estimator improves efficiency by pooling information across units, thereby reducing the variance relative to the individual maximum likelihood estimators $\hat {\theta}^{\operatorname{MLE}}_i$.
These efficiency gains come at the cost of introducing individual bias (i.e., $\mathbb{E}[{\theta}^{*}_i \,|\, \theta_i] \ne \theta_i$), but the reduction in the overall compound mean squared error dominates this cost.
In essence, ${\theta}^{*}_i$ shrinks $\hat {\theta}^{\operatorname{MLE}}_i$ toward the prior.
In practice, the compound decision rule is infeasible due to the unknown $G_*$.
In Bayesian analysis, the prior $G_*$ is typically specified as a parametric distribution,
either by fixing its hyperparameters or by assigning them their own distributions.
\citet{efron-morris:73} show that the James-Stein shrinkage estimator can be interpreted as a parametric EB estimator, with $G_*$ serving as a parametric prior whose hyperparameters are estimated from the data.
An Empirical Bayes analysis can approximate the infeasible oracle decision rule in one of two ways: F-modeling or G-modeling.
F-modeling directly computes the shrinkage term from the data without explicit knowledge of $G_*$.
For example, assuming $\Sigma_i \equiv \Sigma$ for some known homogeneous covariance matrix $\Sigma$ (that is, $\rho_i=\rho$ and $\sigma_i^2 = \sigma^2$ for some known $\rho$ and $\sigma^2$), the (infeasible) F-modeling EB estimator is
\begin{equation}\label{eq:Tweedie}
\hat {\beta}^{\operatorname{F-EB}}_i := \hat {\beta}^{\operatorname{MLE}}_i + (X_i'\Sigma^{-1} X_i)^{-1} \frac{\partial}{\partial \hat{\beta}}\log \hat f_{\hat{\beta}|X} (\hat {\beta}^{\operatorname{MLE}}_i, X_i),
\end{equation}
where $\hat {\beta}^{\operatorname{MLE}}_i = (X_i'\Sigma^{-1} X_i)^{-1}X_i'\Sigma^{-1} Y_i$
and $\hat f_{\hat{\beta}|X}(\hat{\beta}, x)$ is an estimator for $f_{\hat{\beta}|X}(\hat{\beta}, x)$ – the conditional density of $\hat {\beta}^{\operatorname{MLE}}_i$ given $X_i=x$ under $G_*$.
In practice, $\Sigma$ is typically unknown and must also be estimated.
Equation~\eqref{eq:Tweedie}, commonly known as \textit{Tweedie's formula} in the literature, shows that the posterior mean of $\beta_i$ adjusts the corresponding MLE by a shrinkage term derived from the cross-sectional distribution of the MLE.\footnote{In the simplest case when $Y \,|\, \theta \overset{\text{ind.}}{\sim} \mathcal{N}(\theta, 1)$ with $T=1$, the Tweedie formula is $\mathbb{E}_{G_*}[\theta| y] = y + \tfrac{d}{dy} \log f_{G_*}(y)$, where $f_{G_*}(\cdot)$ denotes the marginal density of $Y$ under $G_*$.}
The magnitude of the shrinkage term grows with the variance of $\hat {\beta}^{\operatorname{MLE}}_i$, as measured by $(X_i' \Sigma^{-1} X_i)^{-1}$, and with the informativeness of the prior, as reflected in $f_{\hat{\beta}|X}$.\footnote{ It can be shown that the oracle-decision errors and the shrinkage term are uncorrelated, that is,
$\operatorname{Cov}_{G_*} ({\beta}^{*}_i - \beta_i, \hat {\beta}^{\operatorname{MLE}}_i - {\beta}^{*}_i ) = 0$, where ${\beta}^{*}_i := \mathbb{E}_{G_*}[\beta_i \,|\, Y_i, X_i]$,
implying that the variance reduction is attributable to the shrinkage term in \eqref{eq:Tweedie}.}
F-modeling typically relies on the assumption that the likelihood belongs to the linear exponential family, offering computational convenience when this condition holds.
\citet{Walters:24} provides a survey of F-modeling applications in labor economics.
\citet{LMS:20} estimate the posterior mean of \(a_i\) by F-modeling using Tweedie's formula. However, this method does not readily extend to the HIVDX model.
We adopt the alternative of G-modeling and explicitly estimate the unknown prior $G_*$. \citet{shen2022empirical} demonstrate that G-modeling has superior statistical properties in Poisson mixture models than F-modeling.
While parametric G-modeling imposes a specific functional form on the prior, nonparametric G-modeling characterizes the true distribution $G_*$ as the maximizer of the population log-likelihood:
\begin{equation}
\label{eq:npmle-population-program}
G_* \in \operatornamewithlimits{arg\hspace{0.1em} max}_{G \in \mathcal{G}} F(G),
\end{equation}
where $\mathcal{G}$ is the class of all probability distributions on $\Theta$ and $F(G) := \mathbb{E}[\log f_{G}(Y_i, X_i)]$ denotes the expected marginal log-likelihood (i.e., the negative risk).
This formulation motivates estimating $G_*$ via its sample analog:
\begin{equation}
\label{eq:npmle-sample-program}
\hat G \in \operatornamewithlimits{arg\hspace{0.1em} max}_{G \in \mathcal{{G}}} F_N(G),
\end{equation}
where the empirical criterion is given by $F_N(G) := N^{-1} \sum_{i=1}^N \log f_{G}(Y_i, X_i)$.
We refer to $\hat G$ as the \mbox{NP\hspace{0.05em}MLE}\ of $G_*$.
Substituting $\hat G$ for $G_*$ in \eqref{eq:posterior-mean} yields the G-modeling EB estimator for $\theta_i$:
\begin{equation}\label{eq:posterior-mean-est}
\hat \theta_i^{{\mathrm{EB}}}:= \frac{\int_{\Theta} \theta \ell(Y_i \,|\, X_i, \theta) d\hat G(\theta)}{\int_{\Theta} \ell(Y_i \,|\, X_i, \theta) d\hat G(\theta)}.
\end{equation}
The \mbox{NP\hspace{0.05em}MLE}\ is built on the same principle as parametric MLE, namely that $G_*$ is the maximizer of the population log-likelihood defined in \eqref{eq:npmle-population-program}.
If $G_*$ is the unique maximizer of $F$ over $\mathcal{G}$, it is said to be (point-)identified.
However, it is known that the \mbox{NP\hspace{0.05em}MLE}\ in \eqref{eq:npmle-sample-program} may not be unique in multidimensional settings; in such cases, $\hat{G}$ refers to any maximizer.
Thus, before presenting an algorithm to compute $\hat G$ in Section \ref{sec:implementation}, we first establish conditions for the identification of $G_*$ in Section \ref{sec:identification}, and then show in Section \ref{sec:consistency} that consistent estimation is possible despite multiple sources of heterogeneity.
\section{Identification}\label{sec:identification}
This section establishes identification of the HC model defined in \eqref{eq:model}.
Let $\mathcal{G}$ denote the class of all distributions on $\Theta \equiv \mathbb{R}^{{d_{\beta}}}\times \mathcal{K}_{\delta}$, defined as
\begin{equation}\label{def:G-dist-class}
\mathcal{G} := \mathcal{P}(\mathbb{R}^{{d_{\beta}}}\times \mathcal{K}_{\delta}) = \{ G \in \mathcal{P}(\mathbb{R}^{d_{\theta}}) : \operatorname{supp}_G(\delta_i)\subseteq \mathcal{K}_{\delta}\},
\end{equation}
where $\mathcal{K}_{\delta}\subseteq\mathbb{R}^{{d_{\delta}}}$ is a compact subset in which $\delta_i$ is supported, as in Assumption~\ref{asm:DGP}.
Since any solution to the MLE in \eqref{eq:npmle-population-program} must lie in the set
\begin{equation}\label{eq:G-IDset-via-MLE}
\left\{G \in \mathcal{G} : f_G(Y_i,X_i) = f_{G_*}(Y_i, X_i)\ \ \text{a.s.}\right\},
\end{equation}
the identified set of distributions consists of distributions consistent with the observed marginal likelihood.
We now state the conditions on the covariates and the parametric form of $P(\cdot)$ required for the identification of $G_*$.
\begin{assumption}
\label{asm:id-conditions}
There exists a pair $({\pmb{X}},\pmb{M})$ of fixed matrices with ${\pmb{X}} \in \operatorname{supp}(X_i)$ and $\pmb{M} \in \mathbb{R}^{T\times (T-{d_{\beta}})}$ such that the following conditions hold:
\begin{enumerate}
\item[(i)] $\operatorname{rank}({\pmb{X}}) = {d_{\beta}} \le T-1$, $\operatorname{rank}(\pmb{M}) = T - {d_{\beta}} \ge 1$, and $\pmb{M}'{\pmb{X}} = 0$.
\item[(ii)] For any $\delta, \tilde{\delta} \in \mathcal{K}_{\delta}$,
\[
\pmb{M}'P(\delta)P(\delta)'\pmb{M} = \pmb{M}'P(\tilde{\delta})P(\tilde{\delta})'\pmb{M}
\quad \text{if and only if} \quad \delta = \tilde{\delta}.
\]
\end{enumerate}
\end{assumption}
Assumption~\ref{asm:id-conditions}(i) requires ${\pmb{X}}\in \operatorname{supp}(X_i)$ to have full column rank, paralleling the familiar rank condition in regression.
This ensures the existence of an annihilator matrix $\pmb{M}$ with rank $T - {d_{\beta}}$.\footnote{The choice of $\pmb{M}$ is unique (up to a nonsingular linear transformation) given ${\pmb{X}}$, since a different choice of $\pmb{M}$ can be written as $\tilde{\pmb{M}} = \pmb{M} A$ for some nonsingular $A$.}
We further require that ${d_{\beta}} \le T-1$ so that $\pmb{M}$ has non-null rank and can be used in part~(ii).
This restriction is needed because, in our setting, the covariances $\Sigma_i$ depend on unknown $\delta_i$.
Assumption~\ref{asm:id-conditions}(ii) imposes that $\delta\mapsto \pmb{M}'P(\delta)P(\delta)'\pmb{M}$ is a one-to-one function of $\delta \in \mathcal{K}_{\delta}$, meaning that the variance parameters remain identifiable after removing the effect of ${\pmb{X}} \beta_i$.
The intuition is as follows.
Since the concatenated matrix $({\pmb{X}}, \pmb{M})$ is of full rank, the information in $Y_i$, conditional on $X_i = {\pmb{X}}$, can be expressed as
\begin{equation*}
\left[ \begin{matrix}
\pmb{M}' Y_i \\
{\pmb{X}}' Y_i
\end{matrix} \right]
= \left[ \begin{matrix}
0 \\
{\pmb{X}}' {\pmb{X}} \beta_i
\end{matrix} \right]
+
\left[ \begin{matrix}
\pmb{M}' P_i e_i \\
{\pmb{X}}' P_i e_i
\end{matrix} \right],
\end{equation*}
where $\pmb{M} 'Y_i$ is a Gaussian mixture with mean zero and variance $\pmb{M}'P_iP_i'\pmb{M}=\pmb{M}'P(\delta_i)P(\delta_i)'\pmb{M}$. The distribution of $\pmb{M}'P(\delta_i)P(\delta_i)'\pmb{M}$ can therefore be identified from the distribution of $\pmb{M} 'Y_i$.
Then, the distribution of $\delta_i$ is identified in view of Assumption~\ref{asm:id-conditions}(ii).
To identify $G_*$, the remaining step is to recover the conditional distribution of $\beta_i$ given $\delta_i$.
Let $\mathcal{E} \in \mathbb{R}^{T\times(T-{d_{\beta}})}$.
The key idea is to express
\begin{align*}
\pmb{M}'Y_i + \mathcal{E}' {\pmb{X}}' Y_i & = (\pmb{M}' + \mathcal{E}' {\pmb{X}}')P(\delta_i) e_i + \mathcal{E}'({\pmb{X}}'{\pmb{X}}) \beta_i
\end{align*}
as a perturbation of $\pmb{M}'Y_i$ with noise $\mathcal{E}' {\pmb{X}}' Y_i$.
We measure the shift in the mixing distribution of $\pmb{M}'Y_i + \mathcal{E}' {\pmb{X}}' Y_i$ relative to the mixing distribution of $\pmb{M}'Y_i$ identified earlier.
This comparison reveals the conditional distribution of $\mathcal{E}'({\pmb{X}}'{\pmb{X}}) \beta_i$ given $\delta_i$ as $\mathcal{E}$ traverses $\mathbb{R}^{T\times(T-{d_{\beta}})}$, which then identifies $\beta_i$ given $\delta_i$ via the Cramer-Wold device.
The following theorem formalizes this result.
The proof of Theorem~\ref{thm:identification} as well as the proofs of all other theorems and propositions are given in Appendix~\ref{sec:proof}.
\begin{thm}
\label{thm:identification}
Under Assumptions~\ref{asm:DGP} and \ref{asm:id-conditions}, $G_*$ is identified in the class $\mathcal{G} = \mathcal{P}(\mathbb{R}^{{d_{\beta}}} \times \mathcal{K}_\delta)$.
\end{thm}
The primary challenge in identifying $G_*$ is to distinguish the mean and variance components that are convolved into the outcome distribution.
For instance, assume that $Y_i \sim \mathcal{N}(0,1)$ is a scalar observation.
The observed distribution of $Y_i$ can be rationalized by a continuum of models: $Y_i = a_i + \sigma_i e_i$, where $a_i \sim\mathcal{N}(0,v_a)$ and $\sigma_i\equiv \sqrt{1-v_a}$ for $v_a \in [0,1]$.
Thus, one cannot tell which model generates the observed distribution without imposing further restrictions on $a_i$ or $\sigma_i$.
\citet{bruni1985identifiability} establish general identification results for Gaussian mixtures under a compact support assumption for both the mean and variance components.\footnote{The results of \citet{bruni1985identifiability} are also used in identification of latent distributions in labor economic analyses (see, e.g., \cite{Pastorino:2024,bunting2024,dePaula2025}), but the context is different from identification of mixture models in econometrics such as discussed in \cite{Compiani:Kitamura:2016}.}
Using a type of identification-at-infinity argument, they show that both components can be separated by leveraging variation in the tails of $Y_i$.
However, this argument breaks down as soon as the mean component places a slight nonzero mass in the tails, producing observationally indistinguishable yet fundamentally different models, as illustrated above.
In contrast, our identification strategy rests on the rank and one-to-one restrictions that preclude this type of identification failure.
Theorem~\ref{thm:identification} extends to models that include homogeneous slope parameters:
\begin{equation}
\label{eq:extended-model}
Y_i = X_i \beta_i + W_i \gamma_* + P_i e_i,
\end{equation}
where \(W_i = (W_{i1}, \ldots, W_{iT})'\) is a \(T \times d_{\gamma}\) matrix of additional controls associated with common slopes $\gamma_* \in \mathbb{R}^{{d_{\gamma}}}$.
The common parameters $\gamma_*$ can often be identified without knowledge of $G_*$ by pooling cross-sectional information.
To see how $\gamma_*$ can be identified, rewrite \eqref{eq:extended-model} as a pooled-effects model:
\begin{equation*}
Y_{i} = X_{i} \beta_* + W_{i}\gamma_* + u_{i},
\end{equation*}
where $\beta_* := \mathbb{E}_{G_*}[\beta_i]$ and $u_i = X_i(\beta_i - \beta_*) + P_i e_i$ denotes the composite error.
The parameter $\gamma_*$ is identified as the familiar regression coefficient
\begin{equation*}
\gamma_* = \mathbb{E}[\tilde{W}_i' W_i]^{-1}\mathbb{E}[\tilde{W}_i' Y_i],
\end{equation*}
provided that $\mathbb{E}[\tilde{W}_i'\tilde{W}_i] > 0$, where $\tilde{W}_i := W_i - X_i \mathbb{E}[X_i'X_i]^{-1}\mathbb{E}[X_i' W_i]$.
Once $\gamma_*$ is identified, the model can be analyzed using the residuals $Y_i - W_i \gamma_*$.
Theorem~\ref{thm:identification} serves as a general tool for analyzing identification in HC models. Before applying Theorem~\ref{thm:identification} to the HIVDX model, we first consider a HIVD model without the covariate $X_{2i}$. We verify Assumption \ref{asm:id-conditions} by demonstrating how to choose a pair $({\pmb{X}},\pmb{M})$.
\subsection{Identification of the HIVD Model}
We consider two variants of the HIVD model: one with AR(1) errors and the other with ARMA(1,1) errors.
\subsubsection{AR(1) Errors}
\label{subsubsec:HIVD-AR1}
Consider the following AR(1) panel model without covariates:
\begin{equation}
\label{eq:simple-dynamic-panel}
\begin{aligned}
Y_{it} & = a_i + u_{it}, \\
u_{it} & = \rho_i u_{it-1} + \sigma_i e_{it}, \quad t = 1,\ldots, T,
\end{aligned}
\end{equation}
where the initial condition $u_{i0}$ is assumed to be drawn from the stationary distribution $\mathcal{N}(0, \sigma_i^2/(1-\rho_i^2))$.
The parameters are collected in $\theta_i = (a_i, \sigma_i^2, \rho_i) \overset{\mathrm{i.i.d.}}\sim G_*$.
\begin{prop}
\label{prop:dyp-identification}
Assume that $T \ge 3$.
Consider the panel AR(1) model given in \ref{eq:simple-dynamic-panel}.
$G_*$ is identified in $\mathcal{G} = \mathcal{P}(\Theta)$, where $\Theta = \mathbb{R} \times \mathcal{K}_{(\sigma^2,\rho)}$
and $\mathcal{K}_{(\sigma^2,\rho)}$ is a compact subset of $(0, \infty) \times (-1,1)$.
\end{prop}
The restriction on $\mathcal{K}_{(\sigma^2,\rho)}$ implies that $\rho_i$ is bounded away from $\{-1,1\}$ and that $\sigma_i^2$ is bounded away from $0$ and infinity.
The condition $T \ge 3$ is not only sufficient but also necessary for the identification of the three-dimensional parameter $\theta_i$.
The variance parameters are identified from the covariance structure after purging the fixed effect $a_i$ by first-differencing $Y_i$ via a suitably chosen $\pmb{M}$.
Finally, the compact support for $\delta_i$ ensures that $Y_i$ has non-degenerate and bounded covariance matrices across all $\delta_i$; in particular, this condition rules out unit root processes.
Our choice of $\pmb{M}$ can be illustrated with $T=3$.
Let ${\pmb{X}}$ be the column vector of ones whose rank is $1$.
Set $\pmb{M}$ equal to the first-difference matrix:
\begin{equation*}
\pmb{M} =
\begin{pmatrix}
-1 & 0 \\
1 & -1 \\
0 & 1
\end{pmatrix}.
\end{equation*}
This choice of $\pmb{M}$ yields
\begin{equation*}
\mathcal{V} := \pmb{M}' P(\delta)P(\delta)'\pmb{M}
= \frac{\sigma^2}{1+\rho}
\begin{pmatrix}
2 & -1+\rho \\
-1+\rho & 2
\end{pmatrix}.
\end{equation*}
Since $\rho = 1 + 2\mathcal{V}_{21}/\mathcal{V}_{11}$ and $\sigma^2 = \mathcal{V}_{11} + \mathcal{V}_{21}$, there exists a one-to-one relationship between $\mathcal{V}$ and $\delta$, which ensures that $\delta$ is identified from $\pmb{M}' P(\delta)P(\delta)'\pmb{M}$.
\subsubsection{ARMA(1,1) Errors}
\label{subsubsec:HIVD-ARMA}
Consider the ARMA(1,1) panel model:
\begin{equation*}
\begin{aligned}
Y_{it} &= a_i + u_{it}, \\
u_{it} &= \rho_i u_{it-1} + \sigma_i (e_{it} + \varphi_i e_{it-1}), \quad t = 1, \dots, T,
\end{aligned}
\end{equation*}
where $e_{it} \overset{\mathrm{i.i.d.}}\sim \mathcal{N}(0,1)$ and the initial condition $u_{i0}$ is drawn from the stationary distribution.
The additional random coefficient $\varphi_i \in [-1,1]$ captures heterogeneity in the MA component.
Under the given assumptions, $u_{i}$ is normally distributed with a conditional covariance matrix $\Sigma_i$, whose $(t,s)$ entry is given by
\begin{equation*}
\Sigma_{i,ts} = \frac{\sigma_i^2}{1-\rho_i^2} \times
\begin{cases}
1+\varphi_i^2 + 2 \varphi_i \rho_i & \text{if } t=s, \\
(\rho_i + \varphi_i) (1 + \rho_i \varphi_i) & \text{if } |t-s|=1, \\
\rho_i^{|t-s|-1} (\rho_i + \varphi_i) (1 + \rho_i \varphi_i) & \text{if } |t-s| \ge 2.
\end{cases}
\end{equation*}
\begin{prop}
\label{prop:identification-ARMA}
Let $\theta_i = (a_i, \sigma_i^2, \rho_i, \varphi_i)$ collect all heterogeneous coefficients in the panel model with ARMA(1,1) errors, with $\theta_i \overset{\mathrm{i.i.d.}}\sim G_*$.
Assume that $T \ge 4$ and let $G_*(\rho_i + \varphi_i = 0) = 0$ denote that the probability that $\rho_i + \varphi_i = 0$ equals zero under $G_*$.
Then $G_*$ is identified in $\mathcal{G} = \mathcal{P}(\Theta)$, where $\Theta = \mathbb{R} \times \mathcal{K}_{(\sigma^2,\rho,\varphi)}$
and $\mathcal{K}_{(\sigma^2,\rho,\varphi)}$ is a compact subset of $(0,\infty)\times (-1,1)\times [-1,1]$.
\end{prop}
The condition $T \ge 4$ is necessary to observe the process over a horizon sufficient to distinguish the two sources of persistence, namely $\rho_i$ and $\varphi_i$.
The restriction $G_*(\rho_i + \varphi_i = 0) = 0$ ensures that the two dynamic components do not perfectly offset one another.
We provide the specific choice of $({\pmb{X}},\pmb{M})$ for Assumption~\ref{asm:id-conditions} in Appendix~\ref{app:choice:x:m}.
\subsection{Identification of the HIVDX Model}
We now return to the HIVDX model given in \eqref{def:model:hivxd}.
In this model, the parameters are partitioned as $\beta_i = (a_i, b_i)$ and $\delta_i = (\sigma_i^2, \rho_i)$, yielding the full parameter vector $\theta_i = (a_i, b_i, \sigma_i^2, \rho_i)$ with dimension ${d_{\theta}} = 4$.
Define $\Delta x_t := x_t - x_{t-1}$ as the first difference operator applied to the sequence $(x_t)_{t=1}^T$.
Whether and how $G_*$ can be identified depends critically on the specific time-series variation of the covariates.
We illustrate this using two processes for $X_{2it}$. In the first case, the parameters of the HIVDX model are identified when $T\ge 4$, whereas identification requires $T\ge 5$ in the second case.
\subsubsection{Case I}
We start with the following proposition.
\begin{prop}
\label{prop:HIVDX-id-T=4}
Assume that $T\ge 4$ and either
$
\Delta X_{2,i2} = \pm \Delta X_{2,i3} \ne 0 \text{ or } \Delta X_{2,i2} \ne \Delta X_{2,i4}
$
holds with positive probability.
Then, $G_*$ is identified in $\mathcal{G} = \mathcal{P}(\Theta)$, where $\Theta = \mathbb{R} \times \mathcal{K}_{(\sigma^2,\rho)}$
and $\mathcal{K}_{(\sigma^2,\rho)}$ is a compact subset of $(0, \infty) \times (-1,1)$.
\end{prop}
To illustrate, suppose that $X_{2,it} = X_{2,i1} + (t - 1)$ is a linear trend with a random initial value. Then, we have $\Delta X_{2,it} \equiv 1$ for all $t$.
We need to find an appropriate $({\pmb{X}},\pmb{M})$.
Pick $X_{2,1}$ in $\operatorname{supp}(X_{2,i1})$ and define
\begin{equation*}
{\pmb{X}} := \left( \begin{matrix}
1 & X_{2,1} \\
1 & X_{2,1}+1 \\
1 & X_{2,1}+2 \\
1 & X_{2,1}+3
\end{matrix} \right) ,\quad
\pmb{M} := \left( \begin{matrix}
1 & 0 \\
-2 & 1 \\
1 & -2 \\
0 & 1
\end{matrix} \right),
\end{equation*}
where $\pmb{M}$ represents the second-difference operator in matrix form.
It is clear that ${\pmb{X}}$ has full column rank.
Furthermore, because $\pmb{M}$ annihilates both the intercept (column of ones) and $X_{2, it}$, the transformation $\pmb{M}' Y_i$ eliminates the fixed effects $\beta_i$ entirely.
Let
\begin{equation*}
\mathcal{V} := \pmb{M}' P(\delta)P(\delta)' \pmb{M} \equiv \mathbb{V}(\pmb{M}' Y_i \mid \delta_i = \delta, X_i = {\pmb{X}}),
\end{equation*}
where $\pmb{M}' Y_i \in \mathbb{R}^2$ collects the second differences of the errors, $(\Delta^2 {u}_{i3}, \Delta^2 {u}_{i4})$, conditional on $X_i={\pmb{X}}$.
We obtain the following equations:
\begin{align*}
\mathcal{V}_{11} &= \mathbb{V}(\Delta^2 {u}_{i3} \mid \delta_i = \delta) = \mathbb{V}(u_{i3}-2 u_{i2} + u_{i1} \mid \delta_i = \delta) = \frac{2(3-\rho)}{1+\rho}\sigma^2, \\
\mathcal{V}_{12} &= \operatorname{Cov}(\Delta^2 {u}_{i3}, \Delta^2 {u}_{i4} \mid \delta_i = \delta) = \operatorname{Cov}(u_{i3}-2 u_{i2} + u_{i1}, u_{i4}-2 u_{i3} + u_{i2} \mid \delta_i = \delta) \\
&= \frac{(3-\rho)(\rho-1)}{1+\rho}\sigma^2.
\end{align*}
Hence, the system can be inverted:
\begin{equation*}
(\sigma^2, \rho) = \left(\frac{\mathcal{V}_{11}}{2}\frac{1 + \mathcal{V}_{12}/\mathcal{V}_{11}}{1 - \mathcal{V}_{12}/\mathcal{V}_{11}} , \frac{2 \mathcal{V}_{12}}{\mathcal{V}_{11}} + 1\right).
\end{equation*}
This shows that the mapping $\delta \mapsto \mathcal{V}$ is one-to-one , thereby satisfying Assumption~\ref{asm:id-conditions}(ii).
However, as shown in Appendix \ref{app:choice:x:m},
Proposition~\ref{prop:HIVDX-id-T=4} will not identify the $X_2$ process in Case II that will now be discussed.
\subsubsection{Case II}
\begin{prop}
\label{prop:HIVDX-id-T=5}
Assume that $T \ge 5$ and that $\Delta X_{2,it}$ is not identically zero for some $t \in \{2,\ldots, T\}$.
Then, $G_*$ is identified in $\mathcal{G} = \mathcal{P}(\Theta)$, where $\Theta = \mathbb{R} \times \mathcal{K}_{(\sigma^2,\rho)}$
and $\mathcal{K}_{(\sigma^2,\rho)}$ is a compact subset of $(0, \infty) \times (-1,1)$.
\end{prop}
An example of $X_{2,it}$ that satisfies the condition in the proposition is a one-time universal level shift: $X_{2,it} = \mathbf{1}\{t \ge 3\}$,
as in a nationwide policy change.
Then,
the sequence of first differences is
$
(\Delta X_{2,i2}, \Delta X_{2,i3}, \Delta X_{2,i4}, \Delta X_{2,i5}) = (0, 1, 0, 0).
$
Then, $\pmb{M}$ can now be constructed as
\begin{equation*}
\pmb{M} = \begin{pmatrix}
-1 & 0 & 0 \\
1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 1 & -1 \\
0 & 0 & 1
\end{pmatrix},
\end{equation*}
which yields $\pmb{M}'Y_i = (\Delta u_{i2}, \Delta u_{i4}, \Delta u_{i5})'$. Since the $(1,1)$ and $(2,3)$ elements of the covariance matrix $\mathcal{V} = \mathbb{V}(\pmb{M}'Y_i \mid \delta_i = \delta)$ are given by
\begin{align*}
\mathcal{V}_{11} & = \mathbb{V}(\Delta {u}_{i2} \mid \delta_i = \delta) = \frac{2\sigma^2}{1+\rho},\\
\mathcal{V}_{23} & = \operatorname{Cov}(\Delta u_{i4}, \Delta u_{i5} \mid \delta_i = \delta) = \frac{(\rho-1)\sigma^2}{1+\rho},
\end{align*}
we can identify $\rho$ via the ratio $2\mathcal{V}_{23}/\mathcal{V}_{11} + 1$, and subsequently recover $\sigma^2$.
This establishes that the mapping $\delta \mapsto \mathcal{V}$ is one-to-one, as required by Assumption~\ref{asm:id-conditions}(ii).
The result in Proposition~\ref{prop:HIVDX-id-T=5} also extends to the analysis of heterogeneous treatment implemented at a known time $s \ge 2$, giving
\begin{align*}
Y_{it} = a_i + b_i\,\mathbf{1}\{t \ge s\} + u_{it}, \quad t = 1,\ldots,5.
\end{align*}
Then the term $\mathbf{1}\{t \ge s\}$ represents a common treatment or policy, and $b_i$ captures the individual-specific treatment effect.
As these examples illustrate, incorporating covariates necessitates a significantly more involved analysis.
This is because the transformation matrix $\pmb{M}$ that yields differenced data often obscures the autocovariance structure required to identify $\delta_i$.
It is worth emphasizing that while the sufficient conditions in Proposition~\ref{prop:HIVDX-id-T=5} are simple to verify, the underlying proof is combinatorially complex.
Establishing this general result requires an exhaustive analysis of all distinct patterns of $\Delta X_{2,it}$ to verify the one-to-one mapping between $\delta$ and $\mathcal{V}$ in every possible scenario, as shown in the Appendix.
\section{Consistency of the \mbox{NP\hspace{0.05em}MLE}}
\label{sec:consistency}
This section establishes that $G_*$ can be consistently estimated by the \mbox{NP\hspace{0.05em}MLE}\ under Assumption~\ref{asm:id-conditions} and the additional regularity assumptions stated below.
\begin{assumption}\label{asm:consistency}
\begin{enumerate}
\item [(i)] $\mathbb{E}[\|Y_i\|^2] < \infty$ and $\mathbb{E}[\|X_i\|^2]<\infty$.
\item [(ii)] $\mathbb{P}(\operatorname{rank}(X_i) = {d_{\beta}}) = 1$.
\end{enumerate}
\end{assumption}
Assumption~\ref{asm:consistency}(i) imposes standard moment conditions.
Assumption~\ref{asm:consistency}(ii) requires full column rank of each $X_i$. This ensures that each individual likelihood distinguishes between different values of $\beta_i$, so that the individual MLEs are well-defined.
Our proof of the consistency of the \mbox{NP\hspace{0.05em}MLE}\ builds upon \citet{kiefer1956consistency}, who consider the case where ${d_{\beta}} = 1$.
Since we are interested in cases where ${d_{\beta}} > 1$, we verify a multivariate version of Assumptions~1--5 in \citet{kiefer1956consistency}.
Distinct from the issue of identification, the multivariate setting presents a topological complexity not found in the scalar case: probability mass can diverge along multiple directions.
A primary technical challenge therefore lies in finding a metric that renders the space $\mathcal{G}$ totally bounded (thereby ensuring the compactness of its completion) while enabling the comparison of distributions with such diverging mass.
To this end, the vague topology is the natural choice. We adopt the canonical metrization of this topology, defined by
\begin{equation*}
d(G_0, G_1) = \sum_{r=1}^\infty \frac{1}{2^r} \left| \int_{\Theta} h_r(\theta)dG_0(\theta) - \int_{\Theta} h_r(\theta)dG_1(\theta) \right|,
\end{equation*}
where $\{h_r\}_{r=1}^\infty$ is a dense sequence in the unit ball of $C_c(\Theta)$ with respect to the supremum norm $\|h\|_{\infty}=\sup_{\theta \in \Theta} |h(\theta)|$.\footnote{Here, $C_c(\Theta)$ denotes the space of continuous real-valued functions on $\Theta$ with compact support.}
Since the test functions $h_r$ vanish at infinity, this metric treats any mass that drifts to infinity as if it were collapsed to a single ``point at infinity,'' irrespective of the direction of divergence.
This choice of metric is compatible with the model because the likelihood cannot distinguish between different parameter values at infinity (i.e., the likelihood vanishes as parameters diverge).
In effect, the metric compactifies the domain by collapsing all unbounded directions into a single point, mirroring the likelihood's asymptotic behavior.
Finally, we note that convergence in the vague topology coincides with standard weak convergence whenever the limiting distribution is itself a probability measure.
The next theorem presents the almost sure consistency of the \mbox{NP\hspace{0.05em}MLE}\ with respect to the weak metric $d$.
\begin{thm}
\label{thm:consistency}
Let Assumptions~\ref{asm:DGP}, \ref{asm:id-conditions}, and \ref{asm:consistency} hold.
Then, $d(\hat{G}, G_*) \overset{\operatorname{a.s.}}{\longrightarrow} 0$ as $N \to \infty$.
\end{thm}
Alternative metrics are discussed in the \mbox{NP\hspace{0.05em}MLE}\ literature, though they serve theoretical roles distinct from our metric $d$.
The first is the Hellinger distance between the marginal densities, $\frac{1}{2}\int (\sqrt{f_{\hat{G}}(y)} - \sqrt{f_{G_*}(y)})^2 dy$.
While Hellinger distance is frequently considered in mixture models (\citeauthor{{jiangGeneralMaximumLikelihood2009}}, \citeyear{{jiangGeneralMaximumLikelihood2009}}; \citeauthor{{Bodhi:JRSSB}}, \citeyear{{Bodhi:JRSSB}}), it measures proximity in the data space rather than the parameter space. In our setting, the convergence of the fitted marginal density $f_{\hat{G}}$ to the truth $f_{G_*}$ is of secondary importance.
Another alternative is the Wasserstein distance (e.g., $W_2$ as in \citeauthor{{Bodhi:JRSSB}}, \citeyear{{Bodhi:JRSSB}}).
However, $W_2$ requires the convergence of second moments and consequently becomes unbounded if probability mass escapes to infinity. This renders it unsuitable for our goal of establishing consistency through compactification.
The metric $d$ allows us to apply the Kiefer-Wolfowitz machinery. It strikes a necessary balance: it induces a topology sufficiently weak to ensure compactness (accommodating mass at infinity), yet strong enough to imply the consistency of downstream EB estimation.
The consistency of the \mbox{NP\hspace{0.05em}MLE}\ extends to the constrained \mbox{NP\hspace{0.05em}MLE}, defined as
\begin{equation*}
\tilde{G} \in \operatornamewithlimits{arg\hspace{0.1em} max}_{G \in \tilde{\mathcal{G}}} F_N(G),
\end{equation*}
where $\tilde{\mathcal{G}} \subseteq \mathcal{G}$ is a subset that contains $G_*$.
For instance, \citet{GK:2017} employ a profile likelihood approach to estimate a dynamic panel model similar to \eqref{eq:simple-dynamic-panel} under the restriction $\rho_i = \rho_*$. In our framework, imposing homogeneity ($\rho_i \equiv \rho$) corresponds to optimizing over a restricted subspace $\tilde{\mathcal{G}}$.
The estimator $\hat{G}$ allows for the recovery of various features of the data-generating process that can be expressed as continuous functionals.
In particular, consistent estimators of prior and posterior moments can be obtained from the \mbox{NP\hspace{0.05em}MLE}\ under suitable regularity conditions.
For example, the covariance matrix $\mathbb{V}_{G_*}(\theta_i)$ of the true prior can be estimated by $\mathbb{V}_{\hat{G}}(\theta_i)$, provided that $\|\theta\|^2$ is uniformly integrable with respect to the sequence of estimated priors.
See the discussion following Assumption~\ref{asm:EB} for sufficient conditions ensuring such uniform integrability.
\section{Regret Consistency and $\ell_p$ Convergence}
\label{subsec:EB-consistency}
Let $(\tau_i)_{i=1}^N = (\tau(\theta_i))_{i=1}^N$ denote the unit-level latent parameters of interest, each defined as a function of $\theta_i$.\footnote{We focus on the case of a scalar parameter here. However, the discussion and results in this section readily extend to the vector-valued case.}
The estimation problem consists of constructing a compound decision rule $(Y_i, X_i)_{i=1}^N \mapsto (\hat{\tau}_i)_{i=1}^N$. G-modeling yields the estimator:
\begin{equation*}
\hat{\tau}^{\mathrm{EB}}_i := \mathbb{E}_{\hat{G}}[\tau_i \mid Y_i, X_i] \equiv \frac{\int_{\Theta}\tau(\theta) \ell(Y_i \mid X_i,\theta) d\hat{G}(\theta)}{f_{\hat{G}}(Y_i,X_i)}.
\end{equation*}
This section establishes asymptotic guarantees for $\hat{\tau}^{\mathrm{EB}}_i$ in terms of regret
$R(\hat{\tau}^{\mathrm{EB}}, \tau) - R({\tau}^{*}, \tau)$,
where $R(\hat{\tau}, \tau) := \mathbb{E}_{G_*}[L(\hat{\tau},\tau)]$ is compound risk. Specifically, we consider quadratic loss $L(\hat{\tau}, \tau) :=\frac{1}{N}\sum_{i=1}^N (\hat\tau_i-\tau_i)^2$, and hence the \emph{oracle} posterior mean is defined by:
\begin{equation}
\label{eq:PM-tau}
{\tau}^{*}_i := \mathbb{E}_{G_*}[\tau_i \mid (Y_i,X_i)_{i=1}^N] = \mathbb{E}_{G_*}[\tau_i \mid Y_i, X_i], \quad i=1,\ldots, N,
\end{equation}
where the second equality follows from the independence across units. Write
\begin{equation*}
R(\hat{\tau}^{\mathrm{EB}}, \tau) - R({\tau}^{*}, \tau)
= \mathbb{E}_{G_*}\left[\frac{1}{N} \sum_{i=1}^N (\hat{\tau}^{\mathrm{EB}}_i-\tau_i)^2 - ({\tau}^{*}_i-\tau_i)^2\right] = \mathbb{E}\left[ \frac{1}{N} \sum_{i=1}^N (\hat{\tau}^{\mathrm{EB}}_i - {\tau}^{*}_i)^2\right] \ge 0,
\end{equation*}
since $\mathbb{E}_{G_*}[(\hat{\tau}^{\mathrm{EB}}_i - {\tau}^{*}_i)({\tau}^{*}_i - \tau_i)] = 0$ by the properties of conditional expectation.
The condition $\lim_{N \to \infty} [R(\hat{\tau}^{\mathrm{EB}}, \tau) - R({\tau}^{*}, \tau)] = 0$ is referred to as \emph{regret consistency} (see, e.g., \citet{AGT:EBAdaptive}) or \emph{asymptotic optimality} (see, e.g., \citet{Zhang:1997}).
Consider the case when $\tau_i$ is the prediction of $Y_{i T+1}$.
The oracle predictor ${\tau}^{*}_i$ in the HIVDX model is given by
\begin{align}\label{eq:forecast-formula}
{\tau}^{*}_i &= \mathbb{E}_{G_*}[a_i + b_i X_{2,i T+1} + \rho_i Y_{iT} \mid Y_{i}, X_{i}, X_{2,i T+1}] - \mathbb{E}_{G_*}[\rho_i a_i + \rho_i b_i X_{2,iT} \mid Y_{i}, X_{i}].
\end{align}
The EB one-period-ahead prediction $\hat{\tau}^{\mathrm{EB}}_i$ is then obtained by substituting the \mbox{NP\hspace{0.05em}MLE}\ for $G_*$ in \eqref{eq:forecast-formula}.
Evaluating the EB method in this setting requires controlling the compound error arising from the interaction between the estimation discrepancy $\hat{\theta}^{\mathrm{EB}}_i - {\theta}^{*}_i$ and these observed variables.
Consequently, the analysis necessitates controlling the $\ell_p$ distance (where $p \ge 2$, with $p=4$ being a particularly important case) between the EB and oracle estimates.
\begin{assumption}
\label{asm:EB}
Assume that $\mathbb{E}_{G_*}[|\tau_i|^p] < \infty$ for some $p \in [2, \infty)$. Moreover, the following conditions hold.
\begin{enumerate}
\item [(i)] The function $\tau(\theta)$ is continuous in $\theta \in \Theta$ and, for all $c > 0$, satisfies
\begin{equation*}
\lim_{\|\beta\|\to\infty} \exp(-c \|\beta\|^2) |\tau(\theta)| = 0.
\end{equation*}
\item [(ii)] For all $\epsilon > 0$, there exists $M < \infty$ such that
\begin{align*}
\limsup_{N\to\infty} \mathbb{E} \left[ \int_{\Theta} |\tau(\theta)|^p \mathbf{1}\{|\tau(\theta)| \ge M\} d\hat{G}_N(\theta) \right] &\le \epsilon,
\end{align*}
where $\hat{G}_N$ is the \mbox{NP\hspace{0.05em}MLE}\ defined in \eqref{eq:npmle-sample-program}, with $N$ denoting the sample size.
\end{enumerate}
\end{assumption}
Assumption~\ref{asm:EB}(i) requires that $\tau(\theta)$ depends continuously on $\theta$ and grows slower than any quadratic exponential function of $\|\beta\|$.
This ensures that $|\hat{\tau}^{\mathrm{EB}}_i - {\tau}^{*}_i| \to 0$ for each unit $i$ as $d(\hat{G}_N, G_*) \to 0$. The continuity requirement can be relaxed to allow $\tau(\theta)$ to be discontinuous on a set where $G_*$ places zero mass.
Assumption~\ref{asm:EB}(ii) imposes a high-level condition on the \mbox{NP\hspace{0.05em}MLE}\ sequence, requiring $|\tau_i|^p$ to be uniformly integrable with respect to the sequence of estimated priors.
This condition is automatically satisfied when $\tau_i$ has bounded support (e.g., $\tau_i = \delta_i$) but can be nontrivial for unbounded $\tau_i$ such as slope parameters.
For such cases, Assumption~\ref{asm:EB}(ii) requires specific moment conditions on $(Y_i,X_i)$;
for example, the $p$-th moment of $\beta_i$ in the HIVDX model is bounded by the $p$-th moment of $Y_i$, provided $X_i$ satisfies certain regularity conditions.
\begin{thm}
\label{thm:EB-consistency}
Let Assumptions~\ref{asm:DGP}, \ref{asm:id-conditions}, \ref{asm:consistency}, and \ref{asm:EB} hold for some $p \in [2, \infty)$.
Then $(\hat{\tau}^{\mathrm{EB}}_i)_{i=1}^N$ is regret-consistent, i.e.,
\begin{equation*}
R(\hat{\tau}^{\mathrm{EB}}, \tau) - R({\tau}^{*}, \tau) \to 0 \quad \text{as } N \to \infty.
\end{equation*}
Moreover, the EB estimators converge in $\ell_p$-norm:
\begin{equation*}
\mathbb{E}\left[\frac{1}{N}\sum_{i=1}^N |\hat{\tau}^{\mathrm{EB}}_i-{\tau}^{*}_i|^p\right] \to 0 \quad \text{as } N \to \infty.
\end{equation*}
\end{thm}
Theorem~\ref{thm:EB-consistency} provides the theoretical justification for using $(\hat{\tau}^{\mathrm{EB}}_i)_{i=1}^N$ in large samples.
The primary challenge in establishing this result is that the weak convergence of the priors, $d(\hat{G}, G_*) \overset{\operatorname{a.s.}}{\longrightarrow} 0$, implies the convergence of integrals only for \emph{bounded} continuous functions.
However, $\tau(\theta)$ and the squared-error loss are generally unbounded functions of the parameters.
Assumption~\ref{asm:EB}(ii) is critical because uniform integrability of the estimated prior moments ensures uniform integrability of the posterior moments across units.
This allows us to upgrade the weak convergence of $\hat{G}$ to the convergence of the posterior expectations (the EB estimators) in $\ell_p$-norm.
Notably, the second result of the theorem, that is, convergence in $\ell_p$-norm, is stronger than simple regret consistency ($p=2$) and is particularly useful when the latent parameters enter the decision problem nonlinearly or interact with unbounded covariates.
\subsection{Regret Consistency in the HIVDX Model}
\label{subsec:EB-estimation-HIVDX}
This subsection applies Theorem \ref{thm:EB-consistency} to the HIVDX model, focusing on several quantities of interest: the individual-level parameters $\theta_i$ and the one-period-ahead prediction. These examples provide the theoretical underpinnings for the empirical analysis presented in Section~\ref{sec:empirical}.
When $\tau_i = \theta_i$ and $\theta_i=(\beta_i',\delta_i')'$ in the HIVDX model, the EB estimator of $\theta_i$ can be partitioned as $\hat{\theta}^{\mathrm{EB}}_i = (\hat{\beta}^{\mathrm{EB}}_i, \hat{\delta}^{\mathrm{EB}}_i)$.
Regret consistency of $\hat{\delta}^{\mathrm{EB}}_i$ follows immediately from Theorem~\ref{thm:EB-consistency} and Assumption~\ref{asm:EB}, as $\delta_i$ lies in a compact set.
For $\hat{\beta}^{\mathrm{EB}}_i$, however, Assumption~\ref{asm:EB}(ii) requires $\|\beta_i\|^p$ to be uniformly integrable with respect to the sequence of estimated priors.
This condition is guaranteed by a moment condition on $(Y_i,X_i)$.
Once $(\theta_i)_{i=1}^N$ have been estimated, we turn to the one-period-ahead EB prediction $\hat{Y}^{\mathrm{EB}}_{i T+1}$ of $Y_{i T+1}$. Given the oracle estimator in \eqref{eq:forecast-formula},
the prediction error involves interaction of the errors in estimating the random coefficients and $(Y_i,X_i)$.
Let $\bar{X}_{2,i} = T^{-1}\sum_{t=1}^T X_{2,it}$ denote the time average of $X_{2,it}$ for unit $i$.
\begin{prop}
\label{prop:EB-HIVDX}
Let Assumption~\ref{asm:DGP} hold.
Assume that $|\bar{X}_{2,i}| \le M$ and $\sum_{t=1}^T (X_{2,it} - \bar{X}_{2,i})^2 \ge c$ for some constants $M, c > 0$ and for all $i \ge 1$.
\begin{enumerate}[leftmargin=0.3in]
\item [(i)] $(\hat{\sigma}^{2,\mathrm{EB}}_i)_{i=1}^N$ and $(\hat{\rho}^{\mathrm{EB}}_i)_{i=1}^N$ are regret-consistent.
\item [(ii)] In addition, assume that $\mathbb{E}[\|Y_i\|^{2+\varepsilon}] < \infty$ for some $\varepsilon > 0$. Then $(\hat{a}^{{\mathrm{EB}}}_i)_{i=1}^N$ and $(\hat{b}^{{\mathrm{EB}}}_i)_{i=1}^N$ are regret-consistent.
\item [(iii)] Assume further that $\mathbb{E}[\|Y_i\|^{4+\varepsilon}] < \infty$ for some $\varepsilon > 0$,
$\mathbb{E}[|X_{2,iT}|^{4}] < \infty$, and $\mathbb{E}[|X_{2,i T+1}|^{4}] < \infty$.
Then the predictions $(\hat{Y}^{{\mathrm{EB}}}_{i T+1})_{i=1}^N$ are regret-consistent.
\end{enumerate}
\end{prop}
Proposition~\ref{prop:EB-HIVDX}(iii) requires fourth (and slightly higher) moments of $(X_{2,iT}, X_{2,i T+1})$ and $Y_i$ to control the components of the prediction risk $\mathbb{E}[N^{-1} \sum_{i=1}^N (\hat{Y}^{{\mathrm{EB}}}_{i T+1} - \hat{Y}^{*}_{i T+1})^2]$.
For instance, applying the Cauchy-Schwarz inequality to the component involving the slope coefficient yields
\begin{equation*}
\mathbb{E}\left[ \frac{1}{N}\sum_{i=1}^N (\hat{b}^{{\mathrm{EB}}}_i-\hat{b}^{*}_i)^2 X_{2, i T+1}^2 \right] \le \sqrt{\mathbb{E}\left[ \frac{1}{N}\sum_{i=1}^N (\hat{b}^{{\mathrm{EB}}}_i-\hat{b}^{*}_i)^4\right]} \sqrt{\mathbb{E}\left[ \frac{1}{N}\sum_{i=1}^N X_{2, i T+1}^4\right]}.
\end{equation*}
The moment assumption on $Y_i$ ensures that the first term on the right-hand side vanishes asymptotically (by Theorem~\ref{thm:EB-consistency} with $p=4$), while the fourth moment of $X_{2, i T+1}$ ensures that the second term remains bounded.
\section{Computational Algorithms}
\label{sec:implementation}
A variety of computational methods have been proposed for implementing the \mbox{NP\hspace{0.05em}MLE}. Although the optimization problem is convex, computation is challenging because the optimization variable $G$ is an infinite-dimensional object, rendering direct implementation intractable. A common remedy is to approximate $G$ by a discrete distribution supported on $m$ atoms (with $m \in \mathbb{N}$), thereby reducing the problem to fitting an $m$-component mixture model. This approach is justified by the well-known result that a discrete \mbox{NP\hspace{0.05em}MLE}\ with at most $N_0$ support points exists (possibly among multiple maximizers), where $N_0$ is the number of distinct data points (\citeauthor{{lindsay1983geometry}}, \citeyear{{lindsay1983geometry}}).\footnote{\citet{shen2022empirical} show that the required number of components $m$ grows at the rate $O(\log N)$ in the one-dimensional case under subgaussian priors, and conjecture that in higher dimensions $m$ increases on the order of $O((\log N)^C)$, where $C$ depends on the dimension of $\theta_i$. For one-dimensional Gaussian location mixtures, \citet{Polyanskiy-Sellke:2025} establish computational guarantees, including an algorithm that computes the exact support size.}
Specifically, let $\pmb \theta = (\theta^j)_{j=1}^m$ denote a grid of support points for $G$, and let $\pmb \omega = (\omega^j)_{j=1}^m$ be the associated weights satisfying $\omega^j \ge 0$ for all $j$ and $\sum_{j=1}^m \omega^j = 1$.
Let $G^{\pmb \theta, \pmb \omega}$ denote the discrete distribution induced by $(\pmb \theta, \pmb \omega)$, i.e.,
\[
G^{\pmb \theta, \pmb \omega}(\Theta') = \sum_{j=1}^m \omega^j \, {\mathbf{1}\large[{\theta^j \in \Theta'}\large]}
\qquad \text{for all subsets } \Theta' \subseteq \Theta,
\]
where we use the superscript $j$ to denote the grid $(\theta^j)_{j=1}^m \in \Theta^m$, as distinct from the individual-specific parameters $(\theta_i)_{i=1}^N \in \Theta^N$. The \mbox{NP\hspace{0.05em}MLE}\ of the HC model is given by
\begin{equation*}
\operatornamewithlimits{arg\hspace{0.1em} max}_{\pmb \theta \in \Theta^m, \, \pmb \omega \ge 0}
F_N(G^{\pmb \theta, \pmb \omega})
= \operatornamewithlimits{arg\hspace{0.1em} max}_{\pmb \theta \in \Theta^m, \, \pmb \omega \ge 0}
\frac{1}{N} \sum_{i=1}^N \log f_{G^{\pmb \theta, \pmb \omega}}(Y_i, X_i)
\quad \text{s.t. } \sum_{j=1}^m \omega^j = 1,
\end{equation*}
where
\begin{equation*}
f_{G^{\pmb \theta, \pmb \omega}}(Y_i, X_i)
= \int_{\Theta} \ell(Y_i \,|\, X_i, \theta) \, dG^{\pmb \theta, \pmb \omega}(\theta)
= \sum_{j=1}^m \ell(Y_i \,|\, X_i, \theta^j) \, \omega^j.
\end{equation*}
By the Karush-Kuhn-Tucker theorem, this problem is equivalent to maximizing the Lagrangian
\begin{equation}
\label{eq:numerical-objective}
\mathcal{L}(\pmb \theta, \pmb \omega) =
\frac{1}{N} \sum_{i=1}^N \log \left( \sum_{j=1}^m \ell(Y_i \,|\, X_i, \theta^j) \omega^j \right)
- \lambda \left(\sum_{j=1}^m \omega^j - 1\right).
\end{equation}
The first-order conditions with respect to $\pmb \omega$ imply that the optimal Lagrange multiplier is fixed at $\lambda = 1$.
Most EB estimators are solutions to the discretized problem described above.
For example, \citet{jiangGeneralMaximumLikelihood2009} introduce a fixed-point EM algorithm that mirrors the standard EM routine and establish that it attains a global optimum up to an explicit error bound.
\citet{GK:2017} estimate a dynamic panel model with two-dimensional heterogeneity using the interior-point algorithm of \citet{koenkerConvexOptimizationShape2014} to fit weights on a fixed grid.
More recently, \citet{KCSA:2020}, \citet{ZCST:2024}, and \citet{WIM:2025} propose alternative methods based on general-purpose convex programming to improve scalability.
A feature shared by all these methods is that the location of the support points is fixed, so optimization is performed solely with respect to the weights $\pmb \omega$.
Unlike in simple mixture models, determining where to position the grid points is not obvious when the dimension of $\theta$ exceeds two.
\citet{ZCST:2024} suggest using observed data points as the grid when the dimension is three or higher.
The choice of grid size is also delicate: a grid dense enough to minimize approximation error becomes computationally intractable due to the curse of dimensionality, while a sparse grid risks significant discretization bias.
To overcome the limitations of fixed grids, we consider an algorithm that allows the support points to move adaptively.
This necessitates a shift from optimizing weights on a static grid to optimizing the probability distribution $G$ itself within the space of probability measures $\mathcal{P}(\Theta)$.
In this infinite-dimensional setting, the standard gradient descent used in Euclidean spaces is replaced by \emph{gradient flow}, where the specific trajectory depends critically on the geometric structure imposed on $\mathcal{P}(\Theta)$.
Because $\mathcal P(\Theta)$ can have a complex geometry, different metrics can generate markedly different gradient-flow trajectories.
Consider the gradient flow induced by two commonly used metrics: the Fisher-Rao (FR) and the Wasserstein (W).
The FR gradient flow updates $G_t^{\mathrm{FR}}$ according to
\[
G_{t + \Delta t}^{\mathrm{FR}} - G_t^{\mathrm{FR}} \approx G_t^{\mathrm{FR}}\big(\alpha_t - \mathbb{E}_{G_t^{\mathrm{FR}}}[\alpha_t]\big),
\]
where the reweighting strategy $\alpha_t: \Theta \to \mathbb{R}$ is chosen to maximize the likelihood gain $F_N(G_{t + \Delta t}^{\mathrm{FR}}) - F_N(G_t^{\mathrm{FR}})$ subject to the normalization $\mathbb{E}_{G_t^{\mathrm{FR}}}[(\alpha_t - \mathbb{E}_{G_t^{\mathrm{FR}}}[\alpha_t])^2] = 1$.
Intuitively, the FR flow increases probability mass at locations with high fitness but cannot create mass at new locations or move existing support points.
In contrast, the W gradient flow evolves according to the continuity (transport) equation
\[
G_{t + \Delta t}^{\mathrm{W}} - G_t^{\mathrm{W}} \approx - \mathrm{div}(G_t^{\mathrm{W}} V_t),
\]
where the velocity field $V_t: \Theta \to \mathbb{R}^{{d_{\theta}}}$ describes the movement of mass and is chosen to maximize the gain in $F_N$ under the normalization $\mathbb{E}_{G_t^{\mathrm{W}}}[\|V_t\|^2] = 1$.
This trajectory
effectively moves the support points across the parameter space $\Theta$.
\citet{yan2024learning} show for the Gaussian location mixture case that the FR flow is guaranteed to reach a global maximum when initialized with full support, while $W$ may be stuck at a local maximum. To also exploit the gains from transporting mass across space,
\citet{yan2024learning}
consider a discrete-time {Wasserstein-Fisher-Rao}\ (WFR) gradient flow algorithm that combines both reweighting of FR and transport dynamics of W when updating $G_t^{\mathrm{WFR}}$:
\[
G_{t+ \Delta t}^{\mathrm{WFR}} - G_{t}^{\mathrm{WFR}} \approx - \mathrm{div}(G_{t}^{\mathrm{WFR}} V_t) + G_{t}^{\mathrm{WFR}}\big(\alpha_t - \mathbb{E}_{G_t^{\mathrm{WFR}}}[\alpha_t]\big).
\]
The algorithm adaptively relocates support points so that $m$ remains moderate and need not exceed the sample size $N$. It thus inherits the global convergence guarantees of the FR flow, and adaptivity of W through mass transport. Again in a Gaussian location mixture setting, WFR is shown to be superior to the EM algorithm, and gradient-descent methods based solely on the FR or W flow. Motivated by these encouraging results,
we extend the WFR from a Gaussian location mixture to the HC model.
\subsection{WFR Algorithm for the HC Model}\label{subsec:alg:wfr}
The WFR gradient-descent algorithm updates $(\pmb \theta_n, \pmb \omega_n)_{n \ge 0}$ by alternately reweighting and moving the mass along the steepest descent path.
Let $\eta \in (0,1]$ denote the step size (learning rate), and define the posterior weights
\begin{equation*}
\pi_{ij}(\pmb \theta, \pmb \omega)
= \frac{\ell(Y_i \,|\, X_i, \theta^{j}) \, \omega^{j}}{\sum_{k=1}^m \ell(Y_i \,|\, X_i, \theta^{k}) \, \omega^{k}},
\qquad i=1,\ldots,N,\ j=1,\ldots,m,
\end{equation*}
which correspond to the E-step in the EM algorithm.
Holding $\pmb \theta_n$ fixed, the algorithm first updates the weights $\pmb \omega_n$ via
\begin{align}
\label{eq:FR-step}
\begin{aligned}
\omega_{n+1}^j
&= \omega_{n}^j \left(1 + \eta \left[ \frac{1}{N} \sum_{i=1}^N \frac{\ell(Y_i \,|\, X_i, \theta_n^j)}{\sum_{k=1}^m \ell(Y_i \,|\, X_i, \theta_n^{k}) \, \omega_n^{k}} - 1 \right]\right) \nonumber \\
&= (1-\eta) \omega_n^j + \frac{\eta}{N} \sum_{i=1}^N \pi_{ij}(\pmb \theta_n, \pmb \omega_n),
\qquad j=1,\ldots,m,
\end{aligned} \tag{FR-step}
\end{align}
which we refer to as the FR-step.
The reweighting function is chosen as the gradient of the objective function in \eqref{eq:numerical-objective} with respect to $\omega^j$, thereby ensuring descent.
Next, with $\pmb \omega_{n+1}$ fixed, the grid $\pmb \theta_n$ is updated according to
\begin{align}
\label{eq:W-step}
\begin{aligned}
\tilde{\theta}_{n+1}^j
&= \theta_n^j + \frac{\eta}{N} \sum_{i=1}^N \pi_{ij}(\pmb \theta_n, \pmb \omega_{n+1})
\frac{\partial}{\partial \theta_n^j} \log \ell(Y_i \,|\, X_i, \theta_n^j), \\
\theta_{n+1}^j
&= P_{\Theta}(\tilde{\theta}_{n+1}^j),
\qquad j=1,\ldots,m,
\end{aligned}
\tag{W-step}
\end{align}
where $P_{\Theta}(\cdot)$ denotes the metric projection onto $\Theta$, ensuring that $\theta_{n+1}^j \in \Theta$.\footnote{Formally, for any $\tilde{\theta} \in \mathbb{R}^{d_\theta}$, the metric projection is
$
P_\Theta(\tilde{\theta}) := \operatornamewithlimits{arg\hspace{0.1em} min}_{\theta \in \Theta} \|\tilde{\theta} - \theta\|.
$
It yields the closest point in $\Theta$ to $\tilde{\theta}$ under the Euclidean norm. For example, if $\Theta = [a,b]^{d_\theta}$ is a box constraint, then $P_\Theta(\tilde{\theta})$ simply clips each coordinate of $\tilde{\theta}$ to lie within $[a,b]$.
This ensures that each updated grid point $\theta_{n+1}^j$ remains feasible, even if the gradient step temporarily leaves the parameter space $\Theta$.}
This W-step moves each grid point in the direction of the posterior-weighted average score, corresponding to a discrete version of the W gradient flow.
Starting from \(n = 0\), the FR- and W-steps are applied alternately until \(n = \overline{n}\), the maximum number of iterations, or until the following convergence criterion is met: namely,
\[
\max_{j} \left\{ \frac{1}{N} \sum_{i=1}^N \pi_{ij}(\pmb{\theta}_n, \pmb{\omega}_n) - 1 \right\} \le \mathtt{tol}
\]
for some tolerance parameter \(\mathtt{tol}>0\).\footnote{A necessary and sufficient condition for \(\hat G\) to be \mbox{NP\hspace{0.05em}MLE}\ is
$
\frac{1}{N} \sum_{i=1}^N
\frac{\ell(Y_i \,|\, X_i, \theta)}{f_{\hat G}(Y_i, X_i)} \le 1
\text{ for all } \theta \in \Theta.
$
}
Let \(\hat n\) denote the stopping index (either the first \(n\) satisfying the convergence criterion or \(\overline n\) if it is not met). The resulting estimate is the \mbox{NP\hspace{0.05em}MLE}\ \(\hat G = G^{\pmb{\theta}_{\hat n}, \pmb{\omega}_{\hat n}}\).
Implementation of the HC model is summarized in Algorithm~\ref{alg:iteration}.
\begin{algorithm}[h]
\caption{{Wasserstein-Fisher-Rao}\ algorithm for the HC model}
\label{alg:iteration}
\KwIn{
Dataset $(Y_i, X_i)_{i=1}^N$; initial grid and weights $(\theta_0^j, \omega_0^j)_{j=1}^m$; step size $\eta \in (0,1]$; maximum iterations $\overline{n}$; tolerance parameter \(\mathtt{tol}\).
}
Set $n \leftarrow 0$.
\While{$n < \overline{n} \;\; \mathrm{and} \; \max_{j} \{ \frac{1}{N} \sum_{i=1}^N \pi_{ij}(\pmb{\theta}_n, \pmb{\omega}_n) - 1 \} > \mathtt{tol}$}{
\tcc{FR-step: reweighting}
\For{$j = 1,\ldots,m$}{
Compute posterior weights
$
\pi_{ij}(\pmb\theta_n,\pmb\omega_n)
= \frac{\ell(Y_i \mid X_i,\theta_n^j)\,\omega_n^j}
{\sum_{k=1}^m \ell(Y_i \mid X_i,\theta_n^k)\,\omega_n^k},\quad i=1,\ldots,N.
$
Update
$
\omega_{n+1}^j
= (1-\eta)\,\omega_n^j + \frac{\eta}{N}\sum_{i=1}^N \pi_{ij}(\pmb\theta_n,\pmb\omega_n).
$
}
\tcc{W-step: transport}
\For{$j = 1,\ldots,m$}{
Compute updated posteriors (with new weights)
$
\pi_{ij}(\pmb\theta_n,\pmb\omega_{n+1})
= \frac{\ell(Y_i \mid X_i,\theta_n^j)\,\omega_{n+1}^j}
{\sum_{k=1}^m \ell(Y_i \mid X_i,\theta_n^k)\,\omega_{n+1}^k},\quad i=1,\ldots,N.
$
Take a gradient step
$
\tilde{\theta}_{n+1}^j
= \theta_n^j + \frac{\eta}{N}\sum_{i=1}^N
\pi_{ij}(\pmb\theta_n,\pmb\omega_{n+1})\,
\frac{\partial}{\partial \theta_n^j}\log \ell(Y_i \,|\, X_i,\theta_n^j).
$
Project to the feasible set
$
\theta_{n+1}^j \leftarrow P_{\Theta}(\tilde{\theta}_{n+1}^j).
$
}
Set $n \leftarrow n+1$.
}
\KwOut{Final grid and weights $(\theta_{\hat{n}}^j, \omega_{\hat{n}}^j)_{j=1}^m$, where $\hat{n}$ denotes the stopping index (the first $n$ satisfying the convergence criterion, or $\overline{n}$ if it is not met).}
\end{algorithm}
\paragraph*{Starting values}
Choosing an initial grid that is well spread over $\Theta$ is essential for the WFR gradient-descent algorithm to converge to a global optimum (\citeauthor{{yan2024learning}}, \citeyear{{yan2024learning}}, Theorem~4).
If, instead, all atoms are initialized at the same location, i.e., $\theta_0^j = \tilde{\theta}$ for all $j=1,\ldots,m$, then the subsequent updates satisfy $\theta_n^1 = \cdots = \theta_n^m$ for every $n \ge 0$.
In this case, the distribution $G^{\pmb \theta_n, \pmb \omega_n}$ degenerates to a single point mass at $\theta_n^j = \tilde{\theta}$, failing to capture any heterogeneity in the $\theta_i$’s and hence unable to recover the \mbox{NP\hspace{0.05em}MLE}.
In practice, we generate a dispersed initialization by computing MLEs on random subsamples of the data.
The procedure is as follows.
Fix the subsample size $B$. If the time dimension $T$ is sufficient to ensure that individual MLEs are well-defined, one may choose $B=1$ to maximize the diversity of the grid points.
For $j=1,\ldots,m$:
\begin{itemize}
\item[(i)] Sample a set of indices $\mathcal{I}_j = \{i(1), \ldots, i(B)\}$ uniformly at random from $\{1,\ldots,N\}$.
\item[(ii)] Compute the subsample MLE using $\mathcal{I}_j$:
\[
\theta_0^j = \operatornamewithlimits{arg\hspace{0.1em} max}_{\theta \in \Theta} \sum_{b=1}^B \log \ell(Y_{i(b)} \,|\, X_{i(b)}, \theta)
\]
and assign a uniform weight $\omega_0^j = 1/m$.
\end{itemize}
The resulting collection $(\theta_0^j, \omega_0^j)_{j=1}^m$ serves as the initial grid and weights.
\subsection{Comparison with the EM Algorithm}
\label{subsec:alg:compare}
The EM algorithm is often used to find a (local) solution to mixture problems of the form \eqref{eq:numerical-objective}.
Let $D_{ij} = {\mathbf{1}\large[{\theta_i = \theta^j}\large]}$ denote an unobserved indicator of unit $i$’s membership in component $j$.
If the true $D_{ij}$ were known, the ``complete-data'' maximum likelihood problem would be
\begin{equation*}
\operatornamewithlimits{arg\hspace{0.1em} max}_{\pmb \theta \in \Theta^m}\frac{1}{N} \sum_{i=1}^N \sum_{j=1}^m D_{ij}\log \ell(Y_i \,|\, X_i, \theta^j)
= \operatornamewithlimits{arg\hspace{0.1em} max}_{\pmb \theta \in \Theta^m}\sum_{j=1}^m \sum_{i : D_{ij} = 1} \log \ell(Y_i \,|\, X_i, \theta^j),
\end{equation*}
which can be optimized directly with respect to $\pmb \theta$.
However, this is infeasible because the true $D_{ij}$ are unobserved.
The EM algorithm circumvents this difficulty by replacing $D_{ij}$ with its conditional expectation,
\begin{equation*}
\mathbb{E}_{\pmb \theta, \pmb \omega}[D_{ij} \,|\, Y_i, X_i]
= \frac{\ell(Y_i \,|\, X_i, \theta^{j}) \, \omega^{j}}{\sum_{k=1}^m \ell(Y_i \,|\, X_i, \theta^{k}) \, \omega^{k}} = \pi_{ij}(\pmb \theta, \pmb \omega),
\end{equation*}
evaluated at the current iterate $(\pmb \theta, \pmb \omega)$.
This step is referred to as the E-step.
Let $(\pmb \theta_n^{\mathrm{EM}}, \pmb \omega_n^{\mathrm{EM}})$ denote the $n$th iterate of the support points and their weights.
Given these values, the expected log-likelihood can be written as
\begin{equation*}
Q_N(\pmb \theta \,|\, \pmb \theta_n^{\mathrm{EM}}, \pmb \omega_n^{\mathrm{EM}})
:= \frac{1}{N} \sum_{j=1}^m \sum_{i=1}^N
\pi_{ij}(\pmb \theta_n^{\mathrm{EM}}, \pmb \omega_n^{\mathrm{EM}})
\log \ell(Y_i \,|\, X_i, \theta^j).
\end{equation*}
In the subsequent M-step, $(\pmb \theta_n^{\mathrm{EM}}, \pmb \omega_n^{\mathrm{EM}})$ are updated by solving
\begin{align*}
\theta_{n+1}^{j,\mathrm{EM}}
&= \operatornamewithlimits{arg\hspace{0.1em} max}_{\theta \in \Theta} \frac{1}{N} \sum_{i=1}^N
\pi_{ij}(\pmb \theta_n^{\mathrm{EM}}, \pmb \omega_n^{\mathrm{EM}})
\log \ell(Y_i \,|\, X_i, \theta), \\
\omega_{n+1}^{j,\mathrm{EM}}
&= \frac{1}{N} \sum_{i=1}^N \pi_{ij}(\pmb \theta_n^{\mathrm{EM}}, \pmb \omega_n^{\mathrm{EM}}),
\qquad j=1,\ldots,m.
\end{align*}
The E- and M-steps are iterated until convergence, $(\pmb \theta_n^{\mathrm{EM}}, \pmb \omega_n^{\mathrm{EM}}) \to (\hat{\pmb \theta}^{\mathrm{EM}}, \hat{\pmb \omega}^{\mathrm{EM}})$, yielding the EM estimator $\hat G^{\mathrm{EM}} = G^{\hat{\pmb \theta}^{\mathrm{EM}}, \hat{\pmb \omega}^{\mathrm{EM}}}$.
When $\eta$ in our algorithm is set to $1$, the FR-step coincides with the weight update in the EM algorithm.
For smaller values of $\eta$, the weights are adjusted more conservatively, leading to more stable updates than in EM.
Similarly, the W-step can be viewed as moving $\pmb \theta_n$ incrementally along the gradient of the $Q_N$-function.
Because it is formulated as a gradient-descent procedure, the WFR algorithm enjoys stronger convergence guarantees than EM.
By contrast, the EM algorithm generally requires more careful initialization to reach the global maximum (\citeauthor{{balakrishnanStatisticalGuaranteesEM2017}}, \citeyear{{balakrishnanStatisticalGuaranteesEM2017}}).
\section{Empirical Analysis of Income Dynamics}
\label{sec:empirical}
There is a large body of literature seeking to understand heterogeneity in earnings.
See \citet{Altoniji:2023} and \citet{browning-ejrnaes:13} for recent reviews of this literature.
We highlight several closely related methodological contributions.
\citet{Chamberlain:Hirano:1999} consider a parametric Bayesian approach that allow for heterogeneity in $\sigma_i^2$.
\citet{geweke-keane-joe:00} also adopt a Bayesian framework with a finite mixture model for the composite errors $\sigma_i e_{it}$, while \citet{hirano-ecma:02} extends this framework using an infinite mixture model.
\citet{GK:2017} propose an EB approach that was closest in spirit to ours.
These four papers did not allow for random $(b_i, \rho_i)$.
\citet{Giacomini-et-al:2025} develop an individual-weight shrinkage estimator that optimizes unit-level (rather than aggregate) accuracy by exploiting each individual's own past history, with feasible weights motivated by minimax regret.
\citet{moon-schorfheide-zhang} consider random components in all $(a_i, b_i, \rho_i, \sigma_i^2)$ within a parametric Bayesian framework.
Specifically, they adopt a spike-and-slab prior that accommodate either dense or sparse structures but assume independence across random components.
In contrast, we adopt a nonparametric EB approach for the prior distribution and allow for dependence among the random components.
In this section, we re-examine the application in \citet[GK hereafter]{GK:2017}, which analyzes log real earnings data studied in \citet{meghir2004income}. The extract consists of 938 individuals who have continuous earnings records from age 25 onward in the Panel Study of Income Dynamics (PSID) for the period 1968--1993.
The panel is unbalanced with varying numbers of time periods \(T_i\).\footnote{The only minor difference between the current model and the original HIVDX model in \eqref{def:model:hivxd} is that \(T_i\) may vary across \(i\). The likelihood function can be easily modified to accommodate unbalanced panels.}
As in \citet{GK:2017}, \(Y_{it}\) denotes residuals obtained from year-specific regressions of log real earnings on a vector of covariates.
Though the HIVDX model motivates GK's analysis, they assume $\rho_i$ is homogeneous and $b_i=0$ in the application. Furthermore, their likelihood, expressed in terms of sufficient statistics, is conditioned on the initial observation \(Y_{i1}\), which implicitly assumes that the random component is independent of \(Y_{i1}\). Thus, their marginal likelihood depends on $Y_{i1}$ only through $Y_{i2}-\rho_* Y_{i1}$.
Using GK's PSID sample extract and their definition of $Y_{it}$ (residualized incomes), we estimate the HIVDX model
with $\theta_i := (a_i, b_i, \sigma_i^2, \rho_i)$ and $X_{2,it} := \mathrm{Exp}_{it}/10$, so that $b_i$ is scaled up by a factor of $10$, and
$X_{2,it}$ is potential experience constructed as $\mathrm{Exp}_{it} \;:=\; \text{age} \;-\; \max\{\text{years of schooling},\,12\} \;-\; 6$.\footnote{An equivalent form used in the literature is $\mathrm{Exp}_{it} = \min\{ \text{age} - \text{years of schooling} - 6, \text{age} - 18 \}$.}
Furthermore, our analysis differs from GK in three other respects.
First, we model $Y_{i1}$ under a stationarity assumption without assuming independence between \(Y_{i1}\) and $\theta_i$.
Second, the likelihood is directly evaluated from the observed data since sufficient statistics are not available in the presence of $X_{2,it}$ and $\rho_i$. Third,
we fit the model and estimate the prior distribution of $\theta_i$ using the WFR algorithm described in Section~\ref{subsec:alg:wfr}. In empirical analyses, we used $\overline{n} = 100,000$ as the maximum number of iterations and did not implement the stopping criterion.
The initialization algorithm with $B = 1$ and $m=500$ and the step size of $\eta = 0.0005$ are used.
{
\begin{table}[h]
\normalfont
\centering
\caption{First and Second Moments of the Estimated Prior $\hat{G}$}
\label{tab:PSID:prior:moments}
\medskip
\begin{tabular}{C{0.125\textwidth}
S[table-format=1.4, table-column-width=0.125\textwidth]
S[table-format=1.4, table-column-width=0.125\textwidth]
S[table-format=1.4, table-column-width=0.125\textwidth]
S[table-format=1.4, table-column-width=0.125\textwidth]
}
\toprule
\multirow[c]{2.5}{*}{$\mathbb{V}_{\hat G}(\theta_i)$} & \multicolumn{4}{C{0.5\textwidth}}{$\theta_i$} \\
\cmidrule(lr){2-5}
& \multicolumn{1}{c}{$a_i$} & \multicolumn{1}{c}{$b_i$} & \multicolumn{1}{c}{$\sigma_i^2$} & \multicolumn{1}{c}{$\rho_i$} \\
\midrule
$a_i$ & 0.2052 & -0.1003 & -0.0105 & -0.0039 \\[0.8em]
$b_i$ & -0.1003 & 0.1007 & -0.0104 & 0.0074 \\[0.8em]
$\sigma_i^2$& -0.0105 & -0.0104 & 0.0229 & -0.0013 \\[0.8em]
$\rho_i$ & -0.0039 & 0.0074 & -0.0013 & 0.1052 \\[0.3em]
\cmidrule(lr){1-5}
$\mathbb{E}_{\hat G}[\theta_i]$ & 0.0134 & -0.0166 & 0.0787 & 0.4219 \\[0.25em]
\bottomrule
\end{tabular}
\medskip
{\begin{center}
\parbox{5.5in}{\footnotesize{Notes:
The first four rows form a $4\times 4$ covariance matrix.
The last row presents means.}}
\end{center}
}
\end{table}
}
Table~\ref{tab:PSID:prior:moments} shows the first and second moments from the estimated prior $\hat{G}$.
The marginal variances suggest that the coefficients are heterogeneous in all four dimensions.
The variance of \(a_i\) is substantial, especially relative to the variance of \(Y_{it}\) across all units and time ($0.22$). Regarding \(b_i\), the variance of \(b_i X_{2,it}\) when \(X_{2,it} = x\) is approximately \(0.1 x^2\). The value of $x$ ranges over $[0.2, 3.2]$ within the sample, indicating notable heterogeneity in experience profiles. The variance of \(\sigma_i^2\) is 0.0229 and is non-negligible.
Finally, the distribution of $\rho_i$ is also spread out, with approximately $73\%$ of probability mass lying in $[0.2,0.8]$.
Turning to the covariances of the random coefficients, the covariance between individual intercepts \((a_i)\) and variances \((\sigma_i^2)\) is negative, as in GK. However, allowing for additional sources of heterogeneity leads to some new findings. We find a negative correlation between $\sigma_i^2$ and $b_i$, and weak associations between $\rho_i$ with the other parameters ($a_i,b_i,\sigma_i^2$).
The most interesting result is the strong negative correlation between the individual intercepts
($a_i$) and slopes ($b_i$).
{
\begin{figure}[!htbp]
\caption{Shrinkage Effects: Comparison of $(\hat{\theta}^{\mathrm{MLE}}_i)_{i=1}^N$ and $(\hat{\theta}^{\mathrm{EB}}_i)_{i=1}^N$}
\label{fig:PSID:prior:dist}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_amle.jpg}
\subcaption{MLEs of $a_i$.}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_aeb.jpg}
\subcaption{EB estimates of $a_i$.}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_bmle.jpg}
\subcaption{MLEs of $b_i$.}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_beb.jpg}
\subcaption{EB estimates of $b_i$.}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_rhomle.jpg}
\subcaption{MLEs of $\rho_i$.}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_rhoeb.jpg}
\subcaption{EB estimates of $\rho_i$.}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_sigmle.jpg}
\subcaption{MLEs of $\sigma_i^2$.}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\centering
\includegraphics[width=0.85\linewidth,keepaspectratio]{hist_sigeb.jpg}
\subcaption{EB estimates of $\sigma_i^2$.}
\end{subfigure}
{\begin{center}
\parbox{0.95\textwidth}{\footnotesize Note: The left and right panels display the histograms of the individual MLE and EB estimates for each component of $(\theta_i)_{i=1}^N$.}
\end{center}
}
\end{figure}
}
A side-by-side comparison between the individual MLEs and EB estimates of $\theta_i$ is presented in Figure~\ref{fig:PSID:prior:dist}.
As we work with residual earnings, it is not surprising that the means of the individual MLEs and EB estimates of $a_i$ and $b_i$ are close to zero. A notable feature is that the individual MLEs of $\rho_i$ are highly dispersed, with a substantial fraction taking negative values, whereas the EB estimates are mostly positive and considerably more concentrated around their mean value of 0.422.
The cross-sectional variances of $\hat a_i^{\mathrm{EB}}$, $\hat b_i^{\mathrm{EB}}$, and $\hat \rho_i^{\mathrm{EB}}$ are reduced by approximately 63\%, 72\%, and 48\%, respectively, relative to the corresponding individual MLEs.
While no conspicuous shrinkage is observed for $\hat \sigma_i^{2,\mathrm{EB}}$, these reductions suggest substantial compound mean squared error gains from the EB approach.
Lastly, we consider one-period-ahead predictions of $(Y_{i,T_i})_{i=1}^N$ using a sample that holds out final-period observations.
Figure~\ref{fig:PSID:post:prec} compares the resulting predictions and prediction errors from the individual MLE and the EB methods.
The one-period-ahead prediction of $Y_{i,T_i}$ is computed as
\begin{align}\label{eq:forecast-formula:est}
\hat Y_{i,T_i} &= \hat a_i + X_{2,iT_i} \hat b_i + \hat \rho_i Y_{iT_i-1} - \widehat {\rho_i a_i} - X_{2,iT_i-1} \widehat {\rho_i b_i},
\end{align}
where the hats indicate either EB or individual-level MLE estimates.
As shown in Figure~\ref{fig:PSID:prior:dist}, the EB method exhibits substantial shrinkage, reducing the variance of the predicted values by approximately 25\% compared to the MLE. This translates into improved forecasting performance, leading to a variance reduction of approximately 20\% in the prediction errors.
\begin{figure}[h]
\caption{Pseudo-out-of-sample predictions and prediction errors}
\label{fig:PSID:post:prec}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\subcaption{Predictions: $\hat{Y}^{\mathrm{MLE}}_{i,T_i}$}
\includegraphics[width=0.875\linewidth,keepaspectratio]{fc_mle.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\subcaption{Predictions: $\hat{Y}^{\mathrm{EB}}_{i,T_i}$}
\includegraphics[width=0.875\linewidth,keepaspectratio]{fc_eb.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\subcaption{Errors: $\hat{Y}^{\mathrm{MLE}}_{i,T_i}-Y_{i,T_i}$}
\includegraphics[width=0.875\linewidth,keepaspectratio]{err_mle.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\subcaption{Errors: $\hat{Y}^{\mathrm{EB}}_{i,T_i} - Y_{i,T_i}$}
\includegraphics[width=0.875\linewidth,keepaspectratio]{err_eb.jpg}
\end{subfigure}
{\centering
\parbox{0.95\textwidth}{\footnotesize Notes: The reported standard deviations of the prediction errors are computed after excluding a single outlier observation whose prediction error falls below $-4.5$ under both methods; if this observation is included, the standard deviations for the MLE and EB errors are $0.456$ and $0.411$, respectively.}
}
\end{figure}
As shown in Table~\ref{tab:PSID:prior:moments}, a notable feature is the negative correlation between $a_i$ and $b_i$.
To further examine this relationship, the left panel of Figure~\ref{fig:exp-var-profile} plots the EB estimate pairs $(\hat a_i^{\mathrm{EB}}, \hat b_i^{\mathrm{EB}})$ for $N = 938$ individuals, revealing a clear negative association between the individual intercepts $\hat a_i^{\mathrm{EB}}$ and slopes $\hat b_i^{\mathrm{EB}}$.\footnote{ Figure~\ref{fig:PSID:post:dist} in Appendix~\ref{appendix:pairwise:EB}
presents all pairwise scatter plots of the EB estimates for \(\theta_i\).}
To explore the experience profile in the cross-section, we use the estimated moments in Table~\ref{tab:PSID:prior:moments} to compute the cross-sectional prior variance of \(a_i + b_i x\), which captures the deviation of an individual's income trajectory at \(X_{2,it} = x\) from the average.
The mapping \(x \mapsto \mathbb{V}_{\hat G}(a_i + b_i x)\), shown in the right panel of Figure~\ref{fig:exp-var-profile}, indicates that the variance declines up to about \(x = 1\) (ten years of experience) and increases thereafter.
This non-monotonic pattern arises from the negative covariance between \(a_i\) and \(b_i\).
Interestingly, a similar U-shaped relationship between the log variance of residual earnings and experience is documented in \citet[Chapter~6, p.104]{Mincer:1974}, who attributes it to ``a weak correlation between post-school investment ratios and earning capacity.''
{
\begin{figure}[h]
\caption{EB Estimates of $(a_i, b_i)$ and Cross-Sectional Variance of Earnings}
\label{fig:exp-var-profile}
\begin{center}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=0.99\linewidth,keepaspectratio]{ab_eb.jpg}
\subcaption{Pairs $(\hat a_i^{\mathrm{EB}}, \hat b_i^{\mathrm{EB}})$.}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=0.99\linewidth,keepaspectratio]{exp_prof.jpg}
\subcaption{Variance profile by experience.}
\end{subfigure}
\end{center}
\begin{center}
\parbox{0.95\textwidth}{\footnotesize
Note: The left panel plots the EB estimates $(\hat a_i^{\mathrm{EB}}, \hat b_i^{\mathrm{EB}})_{i=1}^N$.
The right panel displays the cross-sectional variance of earnings as a function of years of experience,
$\text{Exp} \mapsto \mathbb{V}_{\hat G}(a_i + b_i \text{Exp})$.
}
\end{center}
\end{figure}
}
For completeness, we also computed estimates using the EM algorithm discussed in Section~\ref{subsec:alg:compare}, initializing it with the same starting values used for WFR (as described in Section~\ref{subsec:alg:wfr}).
We find that the EM estimates are very similar to the WFR solution.
While the EM algorithm attains a slightly higher likelihood value by optimizing weights over the dense initial grid, the WFR algorithm yields an approximate solution with fewer mass points.
Furthermore, although the EM algorithm generally requires fewer iterations to converge, the computational cost per iteration is higher.
These results are encouraging, as they suggest WFR is well-suited for higher-dimensional problems where fixed-grid algorithms become computationally challenging.
\section{Monte Carlo Experiments}\label{sec:Monte:Carlo}
We conduct Monte Carlo experiments to evaluate the performance of the EB estimator in dynamic panel models relative to (i) the oracle decision rule and (ii) the individual MLE.
The following DGPs are considered.
\begin{enumerate}
\item The HIVD model
with $\theta_i = (a_i, \sigma_i^2, \rho_i)$.
The individual intercepts $(a_i)_{i=1}^N$ are drawn from the $\mathrm{Gamma}(1/2,\sqrt{2})$ distribution independently of $(\sigma_i^2, \rho_i)$.
The joint distribution of $(\sigma_i^2, \rho_i)$ assigns probability mass \(1/6\) to \((\sigma_i^2, \rho_i) = (0.1, 0.8)\) and \((0.3, 0.2)\), and probability mass \(1/3\) to \((\sigma_i^2, \rho_i) = (0.1, 0.2)\) and \((0.3, 0.8)\).
This design features a negative correlation of $-0.33$ between $\rho_i$ and $\sigma_i^2$.
The initial $u_{i0}$ is drawn from the stationary distribution $\mathcal{N}(0,\sigma_i^2/(1-\rho_i^2))$.
The sample size for estimation is $(N,T)=(1000,5)$. An additional observation $Y_{i T+1}$ is generated for each unit to compute prediction errors.
\item The Restricted HIVD model follows the HIVD model above but imposes $\rho_i \equiv 0.5$ and assumes $\sigma_i^2$ has a two-point distribution satisfying $G_*(\sigma_i^2 = 0.1) = G_*(\sigma_i^2 = 0.3) = 0.5$.
The initial distribution of $u_{i0}$ is $\mathcal{N}(0, \sigma_i^2/(1-0.5^2))$.
The distribution of $a_i$ and
the sample size is the same as in HIVD.
This DGP mimics a scenario in which a researcher does not know that $\rho_i$ is a constant and continue to estimate the HIVD model.
\item The HIDVX model generates $\theta_i = (a_i, b_i, \sigma_i^2, \rho_i)$ from the prior distribution estimated from the PSID sample. The covariate (interpreted as experience divided by 10) is assumed to evolve according to $X_{2,it} = X_{2,i,t-1} + 0.1$, with $(N,T) = (1000, 13)$, approximately as in the PSID sample.
The initial value $X_{2,i1}$ is drawn from a discrete distribution with probabilities $\mathbb{P}(X_{2,i1}=0.3) = 0.25$, $\mathbb{P}(X_{2,i1}=0.4) = 0.05$, $\mathbb{P}(X_{2,i1}=0.5) = \mathbb{P}(X_{2,i1}=0.6) = 0.1$, and $\mathbb{P}(X_{2,i1}=0.7) = 0.5$.
\end{enumerate}
For the first two models, the WFR parameters are $\eta = 0.1$, $\overline{n} = 2000$, and $\mathtt{tol} = 10^{-4}$.
For the more complex HIVDX, $\eta = 0.005$, $\overline{n} = 10000$ without applying the early stopping criterion to ensure stable and reliable convergence of the algorithm.
Let $\hat{\varepsilon}_i = \hat{\tau}_i - \tau_i$ denote the estimation error for unit $i$, where $\tau_i$ is one of the model parameters, or the one step ahead prediction.
For each method and replication, we compute
\begin{equation*}
\mathrm{Bias} = \frac{1}{N}\sum_{i=1}^N \hat{\varepsilon}_i,\quad
\mathrm{SD} = \left( \frac{1}{N}\sum_{i=1}^N \left(\hat{\varepsilon}_i - \frac{1}{N}\sum_{i=1}^N \hat{\varepsilon}_i\right)^2 \right)^{1/2},\quad
\mathrm{RMSE} = \left( \frac{1}{N}\sum_{i=1}^N \hat{\varepsilon}_i^{2} \right)^{1/2}.
\end{equation*}
Additionally, a measure of predictability, denoted by R$^2$, is computed as one minus the ratio of the mean squared prediction error to the sample variance of $Y_{i T+1}$.
Table~\ref{tab:MC-HIVD-HIVDR} reports results for the HIVD model and its restricted variant.
In both designs, the oracle estimator achieves the lowest estimation error, but the EB estimator improves upon the individual MLE with uniformly smaller RMSE across all parameters.
The gains are particularly pronounced for $\rho_i$.
The gains from EB are even larger in the restricted HIVD design, where the prior for $\rho_i$ is more informative. Most of the RMSEs of the EB estimator are within 10\% of those of the oracle.
Table~\ref{tab:MC-HIVDX} presents the results for the HIVDX model.
Similarly to Table~\ref{tab:MC-HIVD-HIVDR}, the EB estimator outperforms the individual MLE across all parameters.
For both $a_i$ and $b_i$, the EB method achieves MSEs that are approximately 55\% of those of the MLE.
The EB estimator substantially reduces the MSE for $\rho_i$ as well, primarily attributable to bias reduction.
It is interesting to note that the Monte Carlo results for the bias of the MLE of $\rho_i$ echo the pattern observed in the empirical application in Section~\ref{sec:empirical}. In terms of out-of-sample prediction accuracy, R$^2$ increases from 0.591 for the MLE to 0.662 for EB.
Finally, Table~\ref{tab:prior-moment} shows the prior moment estimation results in the HIVDX design. The first moments of $\theta_i$ are generally well estimated, with the exception of a modest downward bias in the estimate of $\mathbb{E}_{G_*}[\rho_i]$. The variances of $\theta_i$ tend to be somewhat overestimated, which is expected given that estimation noise generally inflates the estimated variance of $\theta_i$ relative to the true variance.
The negative covariance between $a_i$ and $b_i$ is estimated with a downward bias but the sign is reliably recovered.
\begin{table}[htbp!]
\caption{Errors from 500 Monte Carlo Replications}
\label{tab:MC-HIVD-HIVDR}
\centering
\begin{tabular}{
L{2.8cm}
C{1.6cm}
L{1.6cm}
S[table-format=1.3, table-column-width=2cm]
S[table-format=1.3, table-column-width=2cm]
S[table-format=1.3, table-column-width=2cm]}
\toprule
& & & \multicolumn{3}{c}{Estimator} \\
\cmidrule(lr){4-6}
Errors & DGP & Metric & \multicolumn{1}{c}{\text{Oracle}} & \multicolumn{1}{c}{\text{MLE}} & \multicolumn{1}{c}{\text{EB}} \\
\midrule
\multirow{6}{*}{$\hat a_i - a_i$}
& \multirow{3}{*}{HIVD} & Bias & -0.004 & 0.000 & -0.002 \\
& & SD & 0.378 & 0.449 & 0.414 \\
& & RMSE & 0.378 & 0.449 & 0.414 \\
\cmidrule(lr){2-6}
& \multirow{3}{*}{HIVDR} & Bias & 0.001 & 0.000 & -0.003 \\
& & SD & 0.303 & 0.346 & 0.339 \\
& & RMSE & 0.303 & 0.346 & 0.339 \\
\midrule
\multirow{6}{*}{$\hat\sigma_i^{2} - \sigma_i^{2}$}
& \multirow{3}{*}{HIVD} & Bias & 0.000 & -0.077 & -0.006 \\
& & SD & 0.078 & 0.117 & 0.085 \\
& & RMSE & 0.078 & 0.141 & 0.086 \\
\cmidrule(lr){2-6}
& \multirow{3}{*}{HIVDR} & Bias & 0.000 & -0.077 & -0.007 \\
& & SD & 0.079 & 0.118 & 0.085 \\
& & RMSE & 0.079 & 0.141 & 0.086 \\
\midrule
\multirow{6}{*}{$\hat\rho_i - \rho_i$}
& \multirow{3}{*}{HIVD} & Bias & -0.003 & -0.552 & -0.062 \\
& & SD & 0.273 & 0.482 & 0.323 \\
& & RMSE & 0.273 & 0.733 & 0.330 \\
\cmidrule(lr){2-6}
& \multirow{3}{*}{HIVDR} & Bias & 0.000 & -0.546 & -0.048 \\
& & SD & 0.000 & 0.447 & 0.174 \\
& & RMSE & 0.000 & 0.706 & 0.182 \\
\midrule
\multirow{8}{*}{\footnotesize $\hat Y_{i T+1} - Y_{i T+1}$}
& \multirow{4}{*}{HIVD} & Bias & -0.001 & 0.000 & 0.000 \\
& & SD & 0.485 & 0.560 & 0.490 \\
& & RMSE & 0.485 & 0.560 & 0.490 \\
& & $\mathrm{R}^2$ & 0.825 & 0.766 & 0.821 \\
\cmidrule(lr){2-6}
& \multirow{4}{*}{HIVDR} & Bias & -0.001 & 0.000 & -0.001 \\
& & SD & 0.471 & 0.562 & 0.477 \\
& & RMSE & 0.471 & 0.562 & 0.477 \\
& & $\mathrm{R}^2$ & 0.824 & 0.749 & 0.819 \\
\bottomrule
\end{tabular}
\begin{center}
\parbox{0.9\textwidth}{\footnotesize{Notes: Each row reports averages computed from $500$ replications.
Oracle refers to the true posterior mean, and MLE corresponds to the individual MLE.}}
\end{center}
\end{table}
\begin{table}[!htbp]
\caption{Errors from 500 Monte Carlo Replications (HIVDX)}
\label{tab:MC-HIVDX}
\centering
\begin{tabular}{L{2.5cm}L{2cm}
S[table-format=1.3, table-column-width=2.5cm]
S[table-format=1.3, table-column-width=2.5cm]
S[table-format=1.3, table-column-width=2.5cm]}
\toprule
& & \multicolumn{3}{c}{Estimator} \\
\cmidrule(lr){3-5}
Errors & Metric & \text{Oracle} & \text{MLE} & \text{EB} \\
\midrule
\multirow{3}{*}{$\hat a_i - a_i$}
& Bias & 0.001 & -0.001 & -0.001 \\
& SD & 0.244 & 0.426 & 0.318 \\
& RMSE & 0.244 & 0.426 & 0.318 \\
\midrule
\multirow{3}{*}{$\hat b_i - b_i$}
& Bias & -0.001 & 0.001 & 0.001 \\
& SD & 0.188 & 0.338 & 0.251 \\
& RMSE & 0.188 & 0.338 & 0.251 \\
\midrule
\multirow{3}{*}{$\hat\sigma_i^2 - \sigma_i^2$}
& Bias & 0.000 & -0.020 & -0.006 \\
& SD & 0.047 & 0.066 & 0.064 \\
& RMSE & 0.047 & 0.069 & 0.064 \\
\midrule
\multirow{3}{*}{$\hat\rho_i - \rho_i$}
& Bias & -0.002 & -0.326 & -0.082 \\
& SD & 0.212 & 0.308 & 0.275 \\
& RMSE & 0.212 & 0.448 & 0.288 \\
\midrule
\multirow{4}{*}{\footnotesize $\hat Y_{i T+1} - Y_{i T+1}$}
& Bias & 0.000 & 0.001 & 0.001 \\
& SD & 0.299 & 0.341 & 0.314 \\
& RMSE & 0.299 & 0.341 & 0.314 \\
& $\mathrm{R}^2$ & 0.692 & 0.591 & 0.662 \\
\bottomrule
\end{tabular}
\begin{center}
\parbox{0.9\textwidth}{\footnotesize{Notes: Each row reports averages computed from $500$ replications.
Oracle refers to the true posterior mean, and MLE corresponds to the individual MLE.}}
\end{center}
\end{table}
\begin{table}[!htbp]
\caption{Monte Carlo Results for Moment Estimates of $G_*$ (HIVDX)}
\label{tab:prior-moment}
\centering
\begin{tabular}{L{1.6cm} *{6}{S[table-format=-1.3]}}
\toprule\midrule
\multicolumn{7}{l}{\text{Panel A. Means}} \\
& & {$\mathbb{E}_{\hat{G}}[a_i]$} & {$\mathbb{E}_{\hat{G}}[b_i]$} & {$\mathbb{E}_{\hat{G}}[\sigma_i^2]$} & {$\mathbb{E}_{\hat{G}}[\rho_i]$} \\
\midrule
Truth & & 0.013 & -0.017 & 0.079 & 0.422 \\
Bias & &-0.001 & 0.001 & -0.006 & -0.082 \\
SD & & 0.022 & 0.017 & 0.005 & 0.016 \\
RMSE & & 0.023 & 0.018 & 0.011 & 0.080 \\
\midrule\midrule
\multicolumn{7}{l}{\text{Panel B. Variances}} \\
& &{$\mathbb{V}_{\hat{G}}(a_i)$} & {$\mathbb{V}_{\hat{G}}(b_i)$} & {$\mathbb{V}_{\hat{G}}(\sigma_i^2)$} & {$\mathbb{V}_{\hat{G}}(\rho_i)$} \\
\midrule
Truth & & 0.205 & 0.101 & 0.023 & 0.105 \\
Bias & & 0.054 & 0.030 & -0.007 & 0.017 \\
SD & & 0.037 & 0.025 & 0.006 & 0.009 \\
RMSE & & 0.065 & 0.039 & 0.009 & 0.020 \\
\midrule\midrule
\multicolumn{7}{l}{\text{Panel C. Pairwise Covariances}} \\
& {$(a_i,b_i)$} & {$(a_i,\sigma_i^2)$} & {$(a_i,\rho_i)$} & {$(b_i,\sigma_i^2)$} & {$(b_i,\rho_i)$} & {$(\rho_i,\sigma_i^2)$} \\
\midrule
Truth & -0.100 & -0.011 & -0.004 & -0.010 & 0.007 & -0.001 \\
Bias & -0.038 & -0.003 & 0.005 & 0.007 & -0.004 & 0.003 \\
SD & 0.028 & 0.011 & 0.009 & 0.010 & 0.006 & 0.002 \\
RMSE & 0.047 & 0.010 & 0.011 & 0.012 & 0.008 & 0.003 \\
\bottomrule
\end{tabular}
\begin{center}
\parbox{0.95\textwidth}{\footnotesize Notes: The row labeled `Truth' reports the true values of the corresponding moments.
Results are based on 500 replications.}
\end{center}
\end{table}
\section{Conclusions}\label{sec:conclusions}
We have proposed a nonparametric empirical Bayes framework for analyzing short-panel data with rich forms of unobserved heterogeneity. The analysis generalizes classical G-modeling to accommodate heterogeneous slopes and non-spherical error structures. Identification, consistency, and optimality results are established under general conditions, with primitive sufficient conditions derived for specific cases of interest.
However, we have not developed any accompanying inference methods.
Recent contributions such as \citet{AKPM:2022} and \citet{ignatiadis2022confidence} can be useful for future research.
\bibliographystyle{econometrica}
\bibliography{EB.bib}
\clearpage
\newpage