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.
76,182 characters
Estimation of Random-Coefficient Dynamic Panel Data Models with a Fixed $T$
{\hypersetup{pdfborder={0 0 0}}\maketitle}
\vspace{1.5em}
\begin{abstract}
We study dynamic linear panel data models in which the lagged outcomes and strictly exogenous covariates carry individual-specific coefficients and the time-varying errors have a flexible covariance structure.
With a fixed number of time periods, we point-identify the joint distribution of the random coefficients and the structural errors under a distributional form of strict exogeneity, and propose a closed-form, multi-step estimator based on the inverse Radon transform.
We establish a uniform convergence rate for the estimator of the random coefficient density, as well as uniform consistency of the estimator for the conditional density of the time-varying structural errors.
Monte Carlo simulations demonstrate good finite-sample performance of the estimators.
\end{abstract}
\noindent\textbf{Keywords:} dynamic panel data; random coefficients; heterogeneous state dependence; fixed-$T$; inverse Radon transform
\newpage
\section{Introduction}\label{sec:intro}
Dynamic panel data models capture the persistence of economic outcomes, such as employment, earnings, or firm performance, by including lagged dependent variables among the covariates in structural equations.
They are instrumental for disentangling ``genuine'' state dependence, or the direct impact of past outcomes on current ones, from ``spurious'' dependence caused by persistent unobserved heterogeneity across individuals.
These distinct sources of dependence carry sharply different implications, because whether a transitory shock or a temporary policy intervention could leave a lasting imprint depends on the strength of true state dependence.
As such, dynamic panel models have become a workhorse in fields such as labor economics and industrial organization. Canonical applications include dynamic employment in \citet{ArellanoBond1991, BlundellBond1998}, production functions in \citet{blundell2000gmm}, corporate performance and competition in \citet{nickell1996competition}, firm innovation in \citet{aghion2005competition}, market shares and stock valuation in \citet{blundell1999market}, investments and financial factors in \citet{bond1994dynamic,bond2003financial}.
Earlier econometrics literature had studied dynamic panel models with constant autoregressive coefficients where the state dependence is homogeneous across units.
More recent works accommodated heterogeneous state dependence in workers, households, or firms by including random coefficients on the lagged dependent variables, and found it empirically salient in many settings, e.g., labor income dynamics in \citet{FernandezValGaoLiaoVella2022} in a large-$N$, large-$T$ setting.\footnote{
Relatedly, \citet{LiuMoonSchorfheide2020} and \citet{LuMiaoSu2024} applied dynamic panel models with random coefficients on strictly exogenous covariates in the contexts of forecasting and program evaluation.}
Bringing random coefficients to bear on dynamic panels raises two challenges that motivate this paper.
The first is to let the lagged dependent variable carry a random coefficient in a short panel.
Because the lagged outcome is predetermined rather than strictly exogenous, heterogeneous autoregressive dynamics are difficult to identify when the number of time periods $T$ is small.
Nonetheless, knowledge of the distribution of the autoregressive coefficient, rather than a single common value or its mean in the population, is needed for answering many questions.
Examples include measuring the mass of units with non- or near-unit-root persistence (which governs long-run responses), and quantifying the spread and skewness of persistence (which affects the shape of the distribution of forecasts or policy targets).
Whether large-$T$ or fixed-$T$ asymptotics is the more suitable framework depends on the application and the data at hand, and neither approach dominates the other; our aim is to complement this developing literature by offering a methodological alternative that is valid when $T$ is fixed.
The second challenge is to recover the joint distribution of the random coefficients \emph{together with} the time-varying errors (random intercepts in the \emph{structural} form), rather than selected moments or a marginal distribution. The distribution of random coefficients reveals how state dependence correlates with heterogeneous responses to the covariates, both of which are persistent; that of the random structural intercepts captures the distribution of time-varying shocks (permanent and transitory) underlying the process.
Knowledge of their joint distribution is needed not only to separate and quantify the stochastic contribution by permanent vs transitory components, but also to simulate counterfactual trajectories in outcomes such as earnings or productivity. Means or low-order moments from their marginal distributions alone cannot answer these questions.
We study the dynamic linear panel data model
\begin{equation*}
Y_{it} = \gamma_i Y_{it-1} + X_{it}'\beta_{it} + W_{it}'\delta_t + U_{it}, \qquad t = 1,\dots,T,
\end{equation*}
where the autoregressive coefficient $\gamma_i$ and the coefficients $\beta_{it}$ on the strictly exogenous $X_{it}$ are individual-specific, while the coefficients $\delta_t$ on the sequentially exogenous $W_{it}$ are common across units $i$.
The time-varying errors, or ``random structural intercepts'', $U_{it}$, absorb both fixed effects and transitory shocks (e.g., $U_{it}=\alpha_i+\varepsilon_{it}$).
Our object of interest is the joint distribution of $(\gamma_i,\beta_i,U_i)$ conditional on the initial $Y_{i0}$, along with $\delta_t$.
The model identification requires a \emph{distributional} form of strict exogeneity: $(\gamma_i,\beta_i,U_i)$ are jointly independent of the exogenous covariate history conditional on $Y_{i0}$ (Assumption~\ref{assn:indep}).
This condition permits flexible dependence among $\gamma_i$, $\beta_i$, and $U_i$, serial correlation in the errors, and dependence on $Y_{i0}$; it provides the source of variation to identify the full distribution of random coefficients and intercepts (rather than their means) in short panels with fixed $T$.
We impose no restriction on the covariance structure of the time-varying errors, which may be serially correlated and heteroskedastic in the initial condition $Y_{i0}$.
We show that both the joint distribution of $(\gamma_i,\beta_i,U_i)$ given $Y_{i0}$ and the coefficients for sequentially exogenous covariates $\delta_t$ are point identified.
Applying the analog principle to a constructive identification strategy, we propose a closed-form, multi-step estimator which does not involve any numerical optimization, simulation, or solution of an inverse problem beyond some deconvolution steps.
We establish uniform convergence of our estimator for the joint density of $(\gamma_i,\beta_i)$ and that of $U_i$ conditional on the random coefficients.
Our analysis departs from the literature on three fronts, which we preview here and develop in Section~\ref{sec:lit-review}. First, on the \emph{model}: we let the coefficient on the lagged dependent variable be individual-specific and identify its distribution under a \emph{fixed} $T$.
This is the regime in which \citet{Arellano2001} showed that even the mean of a heterogeneous autoregressive coefficient is under-identified when the covariates are only sequentially exogenous.
We circumvent this obstacle by bringing in strictly exogenous covariates under the distributional independence condition above, and we do so without any order condition relating $T$ to the number of covariates, and without the large-$T$ asymptotics that such panels often employ.
Second, on the \emph{target}: we recover the joint distribution of the random coefficients and time-varying errors, rather than a mean, an average partial effect, or a finite set of moments of the random coefficients alone, which the nearest dynamic random-coefficient panel literature targets, often via second-moment conditions, averaging, or partial identification.
In particular, we identify the distribution of the structural random intercepts, i.e., the time-varying errors in the \emph{structural} equation, and not merely that of the slope coefficients.
Triangular and simultaneous-equations models with random coefficients recover the coefficient distributions but leave that of the structural intercepts unidentified; we overcome this obstruction using a new argument that exploits the panel structure.
Third, on the \emph{method}: we adapt the inverse Radon transform used to estimate random-coefficient densities in static cross-sections \citep{hoderlein2010analyzing}, coupled with a deconvolution step, and apply them to an intermediate reduced-form density implied by our panel.
Recovering the structural objects from this reduced-form density is not immediate: it requires a Jacobian change of variables that extracts the joint density of $(\gamma_i,\beta_i)$, and an original argument, central to the proof of \Cref{theorem:two}, that identifies the distribution of $U_i$.
These steps constitute a large part of our methodological contribution.
\section{Related Literature} \label{sec:lit-review}
\noindent \textbullet \quad \textbf{Dynamic linear panel data models with random coefficients.}
Several seminal papers estimated dynamic linear panel data models with \emph{constant} coefficients under a fixed-$T$ setting.
These include \citet{AndersonHsiao1982}, \citet{HoltzEakin1988}, \citet{ArellanoBond1991}, \citet{ArellanoBover1995}, and \citet{BlundellBond1998}.
When the autoregressive coefficient is instead heterogeneous, identification becomes markedly harder: \citet{Arellano2001} showed that with only sequentially exogenous covariates and fixed $T$, even the \emph{mean} of that coefficient is under-identified. This motivates the strictly exogenous covariates and the distributional restriction we introduce in Section~\ref{subsec:iden_distr}.
Some papers studied dynamic linear panel models where strictly exogenous covariates have random coefficients but predetermined and sequentially exogenous ones have \textit{constant} coefficients.
\citet{ArellanoBonhomme2012} identified distributional features of the random coefficients for strictly exogenous covariates using second-moment conditions on time-varying errors.
In a forecasting framework, \citet{LiuMoonSchorfheide2020} constructed point predictors for the posterior mean of heterogeneous coefficients of strictly exogenous covariates and deterministic trends, under the assumption of Gaussian innovations and using Tweedie's formula.
We depart from these papers by letting the lagged dependent variable $Y_{it-1}$ carry a random coefficient, by dispensing with the Gaussian or second-moment conditions on the time-varying errors, and by recovering the joint distribution of random coefficients and intercepts.
\citet{Chamberlain2022} showed that in panel data models with sequentially exogenous covariates and a multi-dimensional vector of unobserved heterogeneity, the common finite-dimensional parameters (such as constant coefficients for other covariates) are not point identified.
\citet{Lee2026} studied the partial identification of a dynamic linear panel model with random coefficients for predetermined covariates and lagged dependent variables.
In contrast, we obtain point (rather than partial) identification of the joint distribution of random coefficients and intercepts by exploiting the full distribution (rather than a finite set of moments) of the outcomes.
Within a framework of potential outcomes and treatment effects, \citet{marx2025heterogeneous} obtain a causal interpretation of the classical IV/GMM estimands from the dynamic-panel literature \citep{AndersonHsiao1982, ArellanoBond1991} under a sequential exchangeability condition.\footnote{Sequential exchangeability is a joint restriction on the dynamics of potential outcomes and on treatment selection. \citet{marx2025heterogeneous} combine it with restrictions on treatment-effect heterogeneity to decompose the IV estimand into a convex-weighted aggregate of heterogeneous, history-dependent treatment effects.} They also relax that condition and identify certain average causal effects by imposing homogeneous autoregressive dynamics on the \emph{untreated} potential outcomes.
Our paper differs from \citet{marx2025heterogeneous} in two major ways. First, we study a distinct random-coefficient dynamic panel model and target a different parameter: the joint distribution of the random coefficients and time-varying errors in a model for \emph{observed} outcomes, rather than aggregates of treatment effects in a \emph{potential}-outcomes framework.
Second, we leverage a different source of variation. Our main identification results (Section~\ref{subsec:iden_distr}) hinge on a strong, distributional form of strict exogeneity of the covariates, whereas \citet{marx2025heterogeneous} include no strictly exogenous covariate by design, and instead use lagged instruments to handle covariates that are not strictly exogenous (e.g., treatments that depend on past observed outcomes).
\noindent \textbullet \quad \textbf{Multidimensional unobserved heterogeneity in linear panels.}
Other papers investigated linear panel models with vector, interactive, or correlated heterogeneity, typically focusing on the mean (instead of distribution) of such heterogeneity.
\citet{GrahamPowell2012} estimated average partial effects (or mean of random coefficients) in an ``irregular'' correlated random coefficient panel, using units whose covariates change little over time and requiring an order condition for just-identification (i.e., the number of periods $T$ equals the dimension of covariates).
Apart from the difference in the target parameter, their method does not apply in our case because the mean of composite random intercepts is not additive in a stationary function of the covariate history and a period-specific constant due to sequentially exogenous covariates.
Moreover, our method does not require any order condition on the number of time periods.
\citet{MoonWeidner2017} estimated a dynamic linear panel model that has a lagged dependent variable with a constant coefficient and interactive fixed effects.
\citet{LuSu2023, LuSu2025} studied linear panel models that allowed for random coefficients for sequentially exogenous covariates, but their theoretical frameworks do not include lagged dependent variables.
\citet{CaoJinLuSu2024} allowed for random coefficients on lagged dependent variables, but required a ``large-$N$, large-$T$'' setting where both the cross-section and time dimensions go to infinity.
In contrast with these papers, our method operates in a short panel ``fixed-$T$'' setting.
Furthermore, while these papers focus on estimating the mean (or a homogeneous/average partial effect) of random slope coefficients, we recover the joint distribution of the random coefficients on lagged dependent variables alongside time-varying errors.
Earlier and wider literature on random-coefficient panel models likewise targeted averages or low-order moments in static designs.
Examples include \citet{Swamy1970}, \citet{Chamberlain1992}, \citet{Wooldridge2005}, \citet{Murtazashvili2008}, and \citet{hsiao2008}.
\citet{Laage2024} estimated a correlated random-coefficient panel model with time-varying endogeneity, showing identification of the mean of random coefficients through control variables. Her model does not accommodate random coefficients for a lagged dependent variable.
\citet{li2026identification} estimated the average partial effect and the local average response in a correlated random coefficient panel data model, where regressors can be correlated with time-varying and individual-specific random coefficients.
\noindent \textbullet \quad \textbf{Triangular and simultaneous systems with random coefficients.}
Our dynamic linear panel model shares the structure of a triangular system, because $Y_{i1}$ feeds into $Y_{i2}$ as a predetermined regressor. Yet our paper differs from those on the triangular system, in terms of the target and the identification strategy. We review the core differences here, and relegate details of the comparison to Appendix~\ref{sec:Appendix_lit}.
\citet{HoderleinHolzmannMeister2017} identified the distribution of the random coefficients in a triangular model under two key restrictions: mutual independence between the coefficients and the instruments and exogenous covariates, and independence of the first-stage slopes from the outcome-equation coefficients (their Theorem 7). As we show in Appendix~\ref{sec:Appendix_lit}, these conditions are relaxed in our setting, and we additionally identify the distribution of random intercepts (time-varying errors) rather than that of the slope coefficients alone.
\citet{MastenTorgovitsky2016} identified an instrumental-variables correlated random coefficient model, relying on independence between the instruments and the coefficients together with a monotone relation between the endogenous covariate and a scalar latent control.
Our model does not fit their control-function framework: the first-stage equation is itself of random-coefficient form, which, as \citet{Imbens2007} and \citet{Kasy2011} showed, admits no such reduced-form, monotone representation with a scalar control.
\citet{masten2018random} treated triangular systems as a special case of simultaneous equations. He identified the joint distribution of the random coefficients for the endogenous regressor and the excluded instrument (his Proposition 4, Section 3.3), assuming the instrument is independent of all unobservables conditional on the shared covariates.
His analysis did not identify the distribution of the time-varying errors (random intercepts), nor that of the
random coefficients on the covariates shared in both equations.\footnote{
Indeed, \citet{masten2018random} showed the joint distribution of all structural unobservables is not point identified, because the reduced form is a seemingly-unrelated-regression system with common regressors (his Theorem 4).}
In comparison, our panel setting avoids this obstruction, because each period supplies its own exogenous covariates: a suitable linear combination of the outcomes behaves as a single-equation random-coefficient model in \emph{distinct} regressors. Building on this, and varying the combination weights, we obtain a new result---identification of the full joint distribution of all random coefficients \textit{and} the time-varying
errors. This requires an extended argument that exploits the conditional independence between the exogenous covariates and the random coefficients and intercepts given the initial condition.
\noindent \textbullet \quad \textbf{Estimation of linear RC models via inverse Radon transform.}
\citet{Beran1996} and \citet{hoderlein2010analyzing} estimated the joint density of the random coefficients in static cross-section regressions by inverting the Radon transform (which links the conditional density of the outcome to the coefficient density); the latter did so through a kernel estimator with a special Radon-transform kernel and sharp asymptotics. We use this method for the first step in our estimation.
\section{The Model and Identification} \label{sec:iden}
We consider a dynamic panel data model:
\begin{equation} \label{eq:seq-exo-w}
Y_{it} = \gamma_{i}Y_{it-1} + X_{it}'\beta_{it} + W_{it}'\delta_t + U_{it} \; \text{ for } t=1,2,\ldots,T,
\end{equation}
where \(\gamma_{i}\in\mathbb R\) and \(\beta_{it}\in\mathbb R^{J_x}\) are random coefficients, $ U_{it}\in\mathbb R$ are random intercepts, while $\delta_t\in\mathbb R^{J_w}$ are non-zero constant coefficients. There are no constant intercepts in $X_{it}, W_{it}$.
For a generic random array $\zeta_{it}$, let $ \zeta_i \equiv (\zeta_{i1}',\zeta_{i2}')'$ and $\zeta_i^t \equiv \{\zeta_{is}:s=1,\ldots,t\}$ denote the history up to time $t$, and $\Delta \zeta_{it} \equiv \zeta_{it} - \zeta_{it-1}$.
\medskip
\begin{assumption}[Strict and Sequential Exogeneity with Time Effects] \label{assn:seq-exo}
For $t=1,\ldots,T$,
$$E(U_{it} \mid X_i^T,W_i^t,Y_{i0},\gamma_i) = \mu_t(Y_{i0}) < \infty \; \text{ almost surely,}$$ where $\mu_t(\cdot)$ is an unknown, unrestricted function.
\end{assumption} \medskip
Under this condition, $X_{it}$ is strictly exogenous while $W_{it}$ is sequentially exogenous.
It permits period-specific fixed effects $\mu_t(Y_{i0})$.
\citet{ArellanoBonhomme2012}
investigated a similar specification where the coefficients for strictly exogenous covariates are random but those for predetermined or sequentially exogenous covariates are all \emph{constant}, and $\mu_t(Y_{i0})=0$ for all $t$.
In comparison, our specification differs by allowing a predetermined regressor, namely, the lagged dependent variable $Y_{it-1}$, to have a random coefficient.
While \Cref{assn:seq-exo} restricts the conditional mean of $U_{it}$ given $\gamma_i,Y_{i0}$, it does allow dependence between $U_i$ and $\gamma_i$ through higher moments (as in the correlated and scale designs in our simulation study).
The first part of our identification results (Section \ref{subsec:iden_moment}) does not require conditions on the second moments of time-varying errors used by \citet{ArellanoBonhomme2012}.
On the other hand, our method does rely on a stronger notion of strict exogeneity, i.e., the distributional independence in \Cref{assn:indep} instead of mean independence (Assumption 1) in \citet{ArellanoBonhomme2012}.
The parameters of interest are the joint distribution of random coefficients and intercepts $(\gamma_i, \beta_{i1}',\ldots,\beta_{iT}', U_{i1},\ldots,U_{iT})$ and the constant coefficients $(\delta_1,\ldots,\delta_T)$.
We prove point identification of these parameters via sequential steps.
\subsection{The conditional mean of random coefficients} \label{subsec:iden_moment}
First, we identify the constant coefficients $\delta_t$ for $t=1,\ldots,T$ and the conditional mean of random coefficients.
To fix ideas, focus on the case with $T=2$.
By recursive substitution and first-differencing,
\begin{align*}
\begin{split}
\Delta Y_{i2} & \equiv Y_{i2} - Y_{i1} \\
& = \gamma_i(\gamma_i-1)Y_{i0} + (\gamma_i-1)X'_{i1} \beta_{i1} +
X'_{i2} \beta_{i2} + (\gamma_i-1)W'_{i1} \delta_1 + W'_{i2}\delta_2 + \varpi_i,
\end{split}
\end{align*}
where \(\varpi_i \equiv \Delta U_{i2} + \gamma_iU_{i1}\).
Let $H_i$ be shorthand for the covariates that are exogenous with respect to $U_i\equiv (U_{i1},U_{i2})'$ conditional on the initial condition $Y_{i0}$. That is:
\[H_{i}\equiv (X_{i1}',X_{i2}',W_{i1}')'.\]
We also maintain a conditional mean independence condition on the random coefficients.
Let $\theta_i \equiv (\gamma_i,\beta_i')'$ with $\beta_i\equiv (\beta_{i1}',\beta_{i2}')'$.
\medskip
\begin{assumption} \label{assn:mean-indep}
$E(\theta_i\mid H_i,Y_{i0}) = E(\theta_i\mid Y_{i0})$ and $E(\theta_i\gamma_i\mid H_i, Y_{i0}) = E(\theta_i\gamma_i \mid Y_{i0})$.
\end{assumption} \medskip
Under Assumptions \ref{assn:seq-exo} and \ref{assn:mean-indep},
\[ E( \varpi_i \mid H_{i}, Y_{i0} ) = E\big[\gamma_i\mu_1(Y_{i0}) + \Delta \mu(Y_{i0}) \mid H_i, Y_{i0} \big] = m(Y_{i0}), \]
where $m(Y_{i0})\equiv E(\gamma_i\mid Y_{i0})\mu_1(Y_{i0}) + \Delta \mu(Y_{i0})$ with $ \Delta \mu \equiv \mu_2 - \mu_1 $, and
\begin{equation} \label{eq:levelreg}
E(Y_{i1}\mid H_i,Y_{i0}) = \pi_1(Y_{i0}) + X_{i1}'E(\beta_{i1}\mid Y_{i0}) + W_{i1}'\delta_1,
\end{equation}
with \(\pi_1(Y_{i0}) \equiv E(\gamma_i\mid Y_{i0})Y_{i0} + \mu_1(Y_{i0})\), and
\begin{equation} \label{eq:delta-y2}
E(\Delta Y_{i2} \mid H_i, Y_{i0}) = Z_i'\Psi(Y_{i0}),
\end{equation}
where $Z_i\equiv(\,1,H_i',E(W_{i2}'\mid H_i,Y_{i0})\,)'$ and
\begin{align} \label{defn:Psi}
\Psi(Y_{i0}) \equiv
\begin{pmatrix}
\psi_1(Y_{i0}) \\
\psi_2(Y_{i0}) \\
\psi_3(Y_{i0}) \\
\psi_4(Y_{i0}) \\
\psi_5(Y_{i0})
\end{pmatrix}
\equiv
\begin{pmatrix}
E(\gamma_i^2 - \gamma_i\mid Y_{i0}) Y_{i0} + m(Y_{i0}) \\
E[\beta_{i1}(\gamma_i - 1) \mid Y_{i0}] \\
E(\beta_{i2}\mid Y_{i0}) \\
\delta_1E(\gamma_i-1\mid Y_{i0}) \\
\delta_2
\end{pmatrix}.
\end{align}
That is, the means of $Y_{i1}$ and $\Delta Y_{i2}$ conditional on $H_i,Y_{i0}$ are both linear in $Z_i$. We maintain the following condition on the (conditional) support of $Z_i$. \medskip
\begin{assumption}[Non-singularity] \label{assn:rank}
$E(\,\| Z_i \|^2\mid Y_{i0}\,)< \infty$, and the support of $Z_i$ conditional on $Y_{i0}$ is not contained in any proper linear subspace of $\mathbb R^{2(J_x+J_w)+1} $ almost surely.
\end{assumption} \medskip
A sufficient condition for \Cref{assn:rank} is that
for any nonzero $c\in\mathbb R^{J_w}$, $c'E(W_{i2}\mid H_i,Y_{i0})$ is not almost surely affine in $H_i$ given $Y_{i0}$.
That is, nonlinearity in the conditional mean of $W_{i2}$ is essential.
Under \Cref{assn:rank}, both $E(Z_iZ_i'\mid Y_{i0})$ and $E(R_iR_i'\mid Y_{i0})$ are non-singular almost surely, where $R_i\equiv(1,X_{i1}',W_{i1}')'$ is a sub-vector of $Z_i$.
Therefore, we can identify $\pi_1(Y_{i0})$, $E(\beta_{i1}\mid Y_{i0})$, and $\delta_1$ from \eqref{eq:levelreg} using variation in $R_i$ given $ Y_{i0} $, and identify $\Psi(Y_{i0})$ from \eqref{eq:delta-y2} using variation in $Z_i$ given $Y_{i0}$.
It then follows from \eqref{defn:Psi} that
\begin{align}
& E(\beta_{i2}\mid Y_{i0}) = \psi_3(Y_{i0}),\quad \delta_2 = \psi_5(Y_{i0}), \quad E(\gamma_i\mid Y_{i0}) = \delta_1'\psi_4(Y_{i0})/\|\delta_1\|^2 + 1 , \\
& \mu_1(Y_{i0}) = \pi_1(Y_{i0}) - E(\gamma_i\mid Y_{i0})Y_{i0}, \quad
E(\gamma_i\beta_{i1}\mid Y_{i0}) = \psi_2(Y_{i0}) + E(\beta_{i1}\mid Y_{i0}).
\nonumber
\end{align}
The proposition below collects these identification results based on the mean of $Y_{it}$. \medskip
\begin{proposition} \label{pn:iden-rc-means}
Under Assumptions \ref{assn:seq-exo}, \ref{assn:mean-indep} and \ref{assn:rank}, $\delta_1$, $\delta_2$, $E(\beta_{i1}\mid Y_{i0})$, $E(\beta_{i2}\mid Y_{i0})$, $E(\gamma_i\mid Y_{i0})$, $E(\gamma_i\beta_{i1}\mid Y_{i0})$ and $\mu_1(Y_{i0})$ are identified from $E(Y_{it}\mid H_i,Y_{i0})$ for $t=1,2$.
\end{proposition} \medskip
\begin{remark}
The coefficients for sequentially exogenous covariates $\delta_1,\delta_2$ are all we need for the next step in \Cref{subsec:iden_distr}.
Nevertheless, we note the identification of the other parameters in \Cref{pn:iden-rc-means} is also useful, as they are obtained from conditional mean outcomes without invoking the distributional form of strict exogeneity in \Cref{assn:indep}.
\end{remark}
\begin{remark}
The second moment $E(\gamma_i^2\mid Y_{i0})$ and $\mu_2(Y_{i0})$ enter additively in $\psi_1(Y_{i0})$, and cannot be separately identified from the conditional means without further assumptions. We propose two approaches to identify them separately.
The first is to use a stationarity condition that $ \mu_t(Y_{i0}) $ is identical over $t=1,2$ almost surely. This allows us to recover $E(\gamma_i^2\mid Y_{i0})$ from $\psi_1(Y_{i0})$ using knowledge of parameters identified in \Cref{pn:iden-rc-means}.
The second approach is to strengthen the mean independence of random intercepts and coefficients in Assumptions \ref{assn:seq-exo} and \ref{assn:mean-indep} to a stronger form of distributional independence (\Cref{assn:indep}), which allows us to recover the distribution (and the second moment) of $\gamma_i$ given $Y_{i0}$ as in \Cref{theorem:two}. We can then recover $\mu_2(Y_{i0})$ from $\psi_1(Y_{i0})$, again using results from \Cref{pn:iden-rc-means}.
\end{remark}
\begin{remark}
If we strengthen \Cref{assn:mean-indep} with $E(\theta_i\mid Y_{i0})=E(\theta_i)$, $E(\gamma_i^2\mid Y_{i0})=E(\gamma_i^2)$, and $\mu_t(Y_{i0})=\mu_t$, then the conditional means of $Y_{i1}$ and $\Delta Y_{i2}$ are both \emph{unconditionally} linear in $H_i$ \emph{and} $Y_{i0}$. The unconditional expectations of $U_i,\gamma_i,\beta_i$ as well as $\delta_1,\delta_2$ are identified from simple regressions of $Y_{i1},\Delta Y_{i2}$ that pool over $Y_{i0}$ and include it as a regressor along with $H_i$.
\end{remark}
\begin{remark}
Our method does not need $T$ to be at least as large as the dimension of strictly exogenous covariates; this differs from \citet{ArellanoBonhomme2012} (with $T$ strictly larger than the latter) and \citet{GrahamPowell2012} (with $T$ equal to the latter).
Furthermore, our model differs from these two papers by allowing the distribution of the random coefficients for strictly exogenous covariates to vary over time, and not to be additive in a stationary function of covariate history. (See Appendix \ref{sec:Appendix_lit} for details.)
\end{remark}
\Cref{pn:iden-rc-means} uses conditional mean outcomes for identification, and invites comparison with Section 3 of \citet{marx2025heterogeneous}.
Beyond the difference in target parameters, their identifying assumptions and ours are non-nested. To compare the two, map their treatment $D_{it}$ into our sequentially exogenous covariate $W_{it}$, and their composite error $\theta_t + \alpha_i + \varepsilon_{it}$ into our random intercept $U_{it}$.
First, \citet{marx2025heterogeneous} focus on the IV/GMM estimand and include \emph{no} strictly exogenous covariate, whereas we require at least one.
Second, sequential exchangeability in \citet{marx2025heterogeneous} permits past outcomes to feed into the current treatment, carrying information about the fixed effects beyond $Y_{i0}$ and resulting in $E(U_{i2} \mid D_{i1}, D_{i2}, Y_{i0}) \neq \mu_2(Y_{i0})$ in general, thus violating our \Cref{assn:seq-exo}.
Moreover, the two approaches draw on different variation: we exploit the (conditional) linearity of $E(Y_{it} \mid H_i, Y_{i0})$ in exogenous covariates, whereas \citet{marx2025heterogeneous} use lagged outcomes as instruments in a first-differenced equation.
\subsection{Distribution of random coefficients and intercept} \label{subsec:iden_distr}
Next, we identify the joint distribution of random coefficients and intercepts, using a stronger form of conditional independence.
Throughout this section, we treat $\delta_1,\delta_2$ as known, having been identified in \Cref{pn:iden-rc-means}.
\medskip
\begin{assumption} \label{assn:indep}
\((\gamma_i,\beta_i,U_i) \perp H_i \mid Y_{i0}\),
and the conditional distribution of $(\gamma_i,\beta_i,U_i)$ given $ Y_{i0}$ admits a Lebesgue density almost surely.
\end{assumption} \medskip
\Cref{assn:indep} allows flexible dependence among $\gamma_i$, $\beta_i$ and $U_i$, and is compatible with Assumptions \ref{assn:seq-exo} and \ref{assn:mean-indep}.
It also allows the distribution of random coefficients and intercepts to be heterogeneous in the initial condition $Y_{i0}$.
In the special case with additive fixed effects $U_{it} = \alpha_i + \varepsilon_{it}$, \Cref{assn:indep} accommodates general dependence between $Y_{i0},\gamma_i,\beta_i,\alpha_i$ and $\varepsilon_i\equiv(\varepsilon_{i1},\varepsilon_{i2})$, as well as serial correlation of $\varepsilon_{it}$.
\medskip
Define \(\widetilde Y_{it} \equiv Y_{it} - W'_{it}\delta_t \) for \( t = 1, 2\).
For any $d\equiv(d_1,d_2)\in\mathbb R^2$,
\begin{align}
\label{eq:ty}
d_1\widetilde Y_{i1} + d_2\widetilde Y_{i2}
= C_i(d) + X_{i1}'\underset{S_{i1}(d)}{\underbrace{(\beta_{i1}d_1 + \beta_{i1}\gamma_id_2)}} +
X_{i2}'\underset{S_{i2}(d)}{\underbrace{\beta_{i2}d_2}} + W_{i1}'\underset{S_{i3}(d)}{\underbrace{\delta_1\gamma_id_2}},
\end{align}
where \[C_i(d) \equiv (d_1\gamma_i + d_2\gamma_i^2)Y_{i0} + (d_1+d_2\gamma_i)U_{i1} + d_2U_{i2}.\]
\noindent Let $S_i(d) \equiv (S_{ik}(d):k=1,2,3)$.
Under \Cref{assn:indep},
\begin{equation} \label{eq:cond-indep}
(\;C_i(d),S_{i}(d)\;)\perp H_i \mid Y_{i0} \text{ for any }d\in\mathbb R^2.
\end{equation}
The {\it first} step in identifying the distribution of random coefficients is to recover the joint distribution of $(\,C_i(d),S_i(d)\,)$ conditional on $Y_{i0}$. \medskip
\begin{assumption}\label{assn:richSupp}\leavevmode
\begin{enumerate}[(i)]
\item \label{assn:richSupp:i} The support of \(H_i\mid Y_{i0}\) contains an open ball in $\mathbb R^{2J_x+J_w}$ almost surely.
\item \label{assn:richSupp:ii} For every $d\in\mathbb R^2$, the conditional distribution of $(C_i(d), S_i(d))$ given $Y_{i0}$ is uniquely determined by its moments, and has finite absolute moments of all orders almost surely.
\end{enumerate}
\end{assumption} \medskip
For linear regressions with random coefficients and intercept in reduced form, \citet{masten2018random} (Lemma 2) showed that the moment-determinacy and finite-moment conditions such as (\ref{assn:richSupp:ii}) are sufficient for point-identification, provided these are independent of the regressors whose joint support satisfies an open-ball condition as in (\ref{assn:richSupp:i}).\footnote{
The generic pair $(A,B)$ and regressors $Z$ in \citet{masten2018random} correspond to $ (C_i(d),S_i(d)) $ and $H_i$ in our case, respectively. \citet{masten2018random} also showed that if $\mathrm{supp}(H_i\mid Y_{i0})$ is bounded, then the conditions in (\ref{assn:richSupp:ii}) are necessary for identification.}
A primitive sufficient condition for \eqref{assn:richSupp:ii} is that, conditional on $Y_{i0}$, the components of $(\gamma_i,\beta_i,U_i)$ have sub-Gaussian tails, which permits them to have non-compact support.
Every linear combination of $(C_i(d),S_i(d))$, whose entries include the quadratic term $\gamma_i^2 Y_{i0}$ as well as the products $\gamma_i U_{i1}$ and $\gamma_i \beta_{i1}$, is then sub-exponential
and hence satisfies Carleman's condition, ensuring moment determinacy of the joint conditional distribution for every $d \in \mathbb{R}^2$ \citep{Petersen1982}.
The following lemma follows immediately from the conditional independence we established in (\ref{eq:cond-indep}), \Cref{assn:richSupp}, and Lemma~2 of \citet{masten2018random}.
\medskip
\begin{lemma}\label{lm:one}
Suppose Assumptions \ref{assn:indep} and \ref{assn:richSupp} hold. For all $d\in\mathbb R^2$, the distribution of $(C_i(d),S_i(d))$ given $Y_{i0}$ is identified from that of $(d_1\widetilde Y_{i1} + d_2\widetilde Y_{i2},H_i')$ given $Y_{i0}$.
\end{lemma}
\medskip
The {\it second} step is to recover the joint density of $(C_i(d),\gamma_i,\beta_i)$ given $Y_{i0}$ from that of $(C_i(d),S_i(d))$ using Jacobian transformation. Specifically, choose $j\in\{1,\ldots,J_w\}$ so that the $j$-th component of $\delta_1$, denoted by $\delta_{1,j}$, is nonzero, and define a Jacobian matrix:
\begin{equation} \label{defn:Jacobian}
J_i(d) \equiv
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & d_2\beta_{i1} & (d_1+d_2\gamma_i)I & 0 \\
0 & 0 & 0 & d_2I \\
0 & d_2\delta_{1,j} & 0 & 0
\end{pmatrix}, \ \text{ where } I \text{ is a $J_x$-by-$J_x$ identity matrix.}
\end{equation}
Suppose for $d$ with $d_2\neq 0 $, the determinant of $J_i(d)$ is non-zero almost surely.
(Because the conditional distribution of $\gamma_i$ given $Y_{i0}$ is atomless under \Cref{assn:indep}, $\Pr\{d_1+d_2\gamma_i=0 \mid Y_{i0}\} = 0 $ for any $d_2\neq 0$ almost surely.)
Recover the joint density of $C_i(d),\gamma_i,\beta_i$ conditional on $Y_{i0}=y_0$ as:
\begin{align}\label{eq:iden_CRC}
\begin{split}
&f_{C_i(d),\gamma_i,\beta_{i1},\beta_{i2} \vert Y_{i0}=y_0}(c,\gamma,b_1,b_2) \\
= & f_{C_i(d),S_{i1}(d),S_{i2}(d),S_{i3,j}(d)\vert Y_{i0}=y_0}(c,b_1(d_1+d_2\gamma),b_2d_2,d_2\gamma\delta_{1,j}) \times
\vert \delta_{1,j}d^{J_x+1}_2(d_1 + d_2\gamma)^{J_x}\vert.
\end{split}
\end{align}
This identifies the joint distribution of $(\gamma_i,\beta_i)$ given $Y_{i0}$.
In the {\it third} step, for any $a \equiv (a_1,a_2) \in \mathbb R^2$ and realization of random coefficients $\gamma$, define the map \(\bar d(a,\gamma) \equiv (a_1-a_2\gamma,\,a_2)\), with components \(\bar d_1 \equiv a_1-a_2\gamma\) and \(\bar d_2 \equiv a_2\). Then for any $b\equiv(b_1',b_2')'$ and initial condition $y_0$,
\begin{align}\label{eq:iden_tU}
\begin{split}
& C_i(\bar d) \mid \gamma_i = \gamma, \beta_i=b, Y_{i0} = y_0 \\
\sim & \ (\bar d_1 + \bar d_2\gamma)U_{i1} + \bar d_2 U_{i2} + \psi^* \mid \gamma_i = \gamma, \beta_i=b, Y_{i0} = y_0 \\
\sim & \ a_1U_{i1} + a_2 U_{i2} + \psi^* \mid \gamma_i = \gamma, \beta_i=b, Y_{i0} = y_0,
\end{split}
\end{align}
where \(\psi^* \equiv (\bar d_1 \gamma + \bar d_2\gamma^2) y_0 = \gamma a_1 y_0\) is a known constant.
{\it Lastly}, combine the three steps above and note for any non-zero $a_1,a_2\in\mathbb R$, the distribution of \(a_1U_{i1} + a_2U_{i2}\) conditional on \((\gamma_i,\beta_i,Y_{i0})\) is identified from that of \((\widetilde Y_{i1}, \widetilde Y_{i2}, H_i')\) given $Y_{i0}$.
Therefore, the characteristic function (and the distribution) of \(U_i=(U_{i1},U_{i2})'\) given $(\gamma_i,\beta_i,Y_{i0})$ is identified.\footnote{
We show identification of the characteristic function over nonzero $(a_1,a_2)$. Recovery of the characteristic function at $a$ s.t. $a_1a_2=0$ follows from uniform continuity of characteristic functions.}
In particular, the intercept means $\mu_t(Y_{i0})=E(U_{it}\mid Y_{i0})$ and $E(\gamma_i^2\mid Y_{i0})$---left non-separable by the conditional-mean step---are recovered here as features of this joint law.
The next theorem formalizes the identification result.
\medskip
\begin{theorem}\label{theorem:two}
Suppose Assumptions \ref{assn:seq-exo}, \ref{assn:mean-indep}, \ref{assn:rank}, \ref{assn:indep} and \ref{assn:richSupp} hold. Then $\delta_1,\delta_2$ and the joint distribution of \((U_i,\gamma_i,\beta_i)\) given \(Y_{i0}\) are identified.
\end{theorem}
\medskip
\begin{remark}
\Cref{theorem:two} uses \Cref{assn:seq-exo}--\ref{assn:rank} only to identify $\delta_1,\delta_2$ via \Cref{pn:iden-rc-means}.
Given $\delta_1,\delta_2$, the distributional argument rests on \Cref{assn:indep}--\ref{assn:richSupp} alone.
\Cref{assn:mean-indep} is implied by \Cref{assn:indep} (as long as the conditional means exist).
\end{remark}
The identification strategy above is constructive; we apply the analog principle to define closed-form estimators in the next section.
In Appendix~\ref{sec:T3}, we extend this identification strategy to any fixed $T\ge3$ even without variation in sequentially exogenous $W_{it}$.
\section{Closed-Form Estimators}\label{sec:est}
We propose multi-step estimators for $\delta_1,\delta_2$ and the joint density of $ (\gamma_i, \beta_i, U_i)$.
The identification results are conditional on the initial condition $Y_{i0}$; we fix it at a constant $y_0\neq 0$ throughout this section, and suppress it from the notation. As in Section \ref{sec:iden}, we focus on the case with $T=2$.
\subsection{Moments of random coefficients and \texorpdfstring{$\delta_t$}{deltat}}
\label{subsec:est-delta}
Let \(T=2\) and use lower-case letters to denote the realizations in the sample, e.g., $h_i \equiv (x_{i1}',x_{i2}',w_{i1}')'$.
Let $\widehat w_{i2}$ denote a Nadaraya-Watson estimator for the mean of $W_{i2}$ conditional on $h_i$;
let $\widehat z_i \equiv (1,h_i',\widehat w_{i2}')'$, and
\[\widehat \Psi \equiv \left(\sum\nolimits_i \widehat z_i\widehat z_i'\right)^{-1}
\left(\sum\nolimits_i\widehat z_i\Delta y_{i2}\right),\]
where the $k$-th component $\widehat \psi_k$ is an estimator of $\psi_k$ in (\ref{defn:Psi}), with the initial condition $y_0$ fixed and suppressed in notation.
Let $r_i \equiv (1,x_{i1}',w_{i1}')'$ and estimate the period-1 level regression \eqref{eq:levelreg} by
\begin{equation} \label{eq:est_condY1}
\widehat\Phi \equiv \left(\sum\nolimits_i r_ir_i'\right)^{-1}\left(\sum\nolimits_i r_i\, y_{i1}\right) = (\widehat \pi_1,\ \widehat E(\beta_{i1})',\ \widehat\delta_1')'.
\end{equation}
Together with $\widehat\Psi$ from the regression on $\Delta Y_{i2}$, we also estimate
\begin{align} \label{est_condMean}
\widehat E(\beta_{i2}) \equiv & \ \widehat\psi_3, \quad \widehat\delta_2 \equiv \widehat\psi_5, \quad
\widehat E(\gamma_i) \equiv 1 + \widehat\delta_1'\widehat\psi_4/\|\widehat\delta_1\|^2.
\end{align}
We can also estimate the following parameters (even though they are not needed for subsequent estimation of the joint density of $(\gamma_i,\beta_i',U_i')'$):
\begin{align*}
& \widehat\mu_1 \equiv \widehat \pi_1 - \widehat E(\gamma_i)\,y_0, \quad
\widehat E(\gamma_i^2) \equiv \int \gamma^2\,\widehat f_{\gamma_i}(\gamma)\,d\gamma, \quad \widehat E(\gamma_i\beta_{i1})\equiv \widehat\psi_2+\widehat E(\beta_{i1}), \nonumber \\
& \widehat\mu_2 \equiv \widehat\psi_1 - \big(\widehat E(\gamma_i^2)-\widehat E(\gamma_i)\big)y_0 - \widehat E(\gamma_i)\widehat\mu_1 + \widehat\mu_1, \nonumber
\end{align*}
where $\widehat f_{\gamma_i}(\gamma) \equiv \int \widehat f_{\gamma_i,\beta_i}(\gamma,b)\,db$ uses the density estimator of Section \ref{subsec:est_RC}.
With $\delta_1\neq 0$, the probability limit of the denominator $\|\widehat\delta_1\|^2$ in $\widehat E(\gamma_i)$ is strictly bounded away from zero.
Consequently, a standard argument using Slutsky's Theorem ensures that division by $\|\widehat\delta_1\|^2$ does not affect the $\sqrt{n}$-consistency and asymptotic normality of the estimators in \eqref{eq:est_condY1} and \eqref{est_condMean}.
\subsection{The joint density of random coefficients} \label{subsec:est_RC}
\noindent
For the rest of this section and Section~\ref{sec:asymp}, assume $W_{i1}\in\mathbb R$ for simplicity.
Fix $d=(0,1)$, so that (\ref{eq:ty}) reduces to
\begin{equation*}
\widetilde Y_{i2} = \underbrace{\gamma_i^2 y_0 + \gamma_iU_{i1}+U_{i2}}_{C_i^*}+ X'_{i1}\underbrace{\beta_{i1}\gamma_i}_{S_{i1}^*}+X'_{i2}\underbrace{\beta_{i2}}_{S_{i2}^*}+W'_{i1}\underbrace{\delta_1\gamma_i}_{S_{i3}^*}.
\end{equation*}
Following \cite{hoderlein2010analyzing}, we define an (inverse) Radon transform estimator (RTE) for the density of $(C_i^*,S_i^*)$.
Recall $ H_i \equiv (X_{i1}',X_{i2}',W_{i1}')' $. For $i=1,\dots,n$, define
\begin{equation*}
Q_i \equiv \| (1, H_i')\|^{-1}(1, H_i')' \in \mathbb{S}_+^{l-1} \ \text{ and } \
V_i \equiv \|(1, H_i')\|^{-1} \left(Y_{i2} - W_{i2}'\widehat\delta_2\right) \in \mathbb{R},
\end{equation*}
where $\|\cdot\|$ denotes the Euclidean norm, and $\mathbb{S}_+^{l-1} = \{z \in \mathbb{R}^l: z_1>0, \|z\|=1 \}$ is the upper hemisphere of the unit sphere in $\mathbb{R}^l$ with $l\equiv\dim[(1, H_i')]=2J_x+2$.
Define \(
\mathbb{S}(\underline q_n)\equiv\{q\in\mathbb{S}_+^{l-1}: q_1\ge \underline q_n\}\), where $\underline q_n\to 0$ as $n\to\infty$.
Estimate the density of $(C_i^*,S_i^*)$ by
\begin{equation} \label{eq:est_cstar}
\widehat f_{C_i^*,S_i^*}(c,s) \equiv \frac{2}{n} \sum\nolimits_i \frac{\mathbf 1\{q_{i1}\ge \underline q_n\}}{\widehat f_{Q_i}(q_i)}K_\nu(q_i'(c,s')'-v_i),
\end{equation}
where $\mathbf 1(\cdot)$ is the indicator function, and $K_\nu$ is the kernel defined by Equations (7)--(8) of \cite{bissantz2014confidence}.\footnote{\label{note:defn_K}
The kernel is defined through its Fourier transform $\mathcal F K_\nu(t)=\frac{1}{2}(2\pi)^{-l+1}|t|^{l-1}\mathcal L(\nu|t|)$, which by symmetry of $\mathcal L$ admits an explicit form:
\(K_\nu(w)\equiv (2\pi)^{-l}\int^\infty_0\cos(tw)\,t^{l-1}\mathcal L(\nu t)\,dt\).
The symmetric filter is \(\mathcal L(t)\equiv(1-|t|^r)\mathbf{1}_{[-1,1]}(t)\), with an order parameter $0<r<\infty$.
The kernel $K_\nu$ can be written, via a change of variable $\tilde t=\nu t$, as \(K_\nu(w)=\nu^{-l}(2\pi)^{-l}\!\int_0^\infty\!\cos\!\big(\tilde t\,\tfrac{w}{\nu}\big)\tilde t^{\,l-1}\mathcal L(\tilde t)\,d\tilde t=\nu^{-l}K(\tilde w)\big|_{\tilde w=w/\nu}\), where \( K(\tilde w)\equiv (2\pi)^{-l}\!\int_0^\infty\!\cos(t\tilde w)\,t^{l-1}\mathcal L(t)\,dt \). }
The denominator is a kernel estimator of the density of $Q_i$ on its spherical support:
\begin{equation*}
\widehat f_{Q_{i}}(\tilde q) \equiv \frac{1}{n}\sum\nolimits_i \varkappa(\tau)\widetilde K\left(\tau^{-2}(1-\tilde q'q_i)\right),
\end{equation*}
where $\widetilde K(\cdot)$ is a kernel function, $\tau > 0$ is a smoothing parameter, and $\varkappa(\tau)$ is a normalization constant defined in \cite{hoderlein2010analyzing}.
Estimate the density of $(C_i^*,\gamma_i,\beta_i)$ by
\[\widehat f_{C_i^*,\gamma_i,\beta_i}(c,\gamma,b) \equiv \widehat f_{C_i^*,S_i^*}(c, b_1\gamma,b_2,\widehat\delta_1\gamma) \vert \widehat\delta_1\gamma^{J_x}\vert,\]
and estimate the density of $(\gamma_i,\beta_i)$ by
\[\widehat f_{\gamma_i,\beta_i}(\gamma,b) \equiv \int_{\mathcal C_n} \widehat f_{C_i^*,\gamma_i,\beta_i}(c,\gamma,b)\, d c ,\]
where $\mathcal C_n\equiv[-\bar c_n,\bar c_n] \subset \mathbb{R}$, with $\bar c_n\to\infty$ as $n\to\infty$.
\subsection{The joint density of \texorpdfstring{$U_i$}{Ui} conditional on random coefficients} \label{subsec:est_U}
For any $a = (a_1,a_2) \in \mathbb{R}^2$ and any given $\gamma$ on the support of $\gamma_i$, define $\bar d \equiv (a_1-a_2\gamma, a_2)$.
With slight abuse of notation, write $C_i(\bar d) = C_i(a,\gamma)$ and $S_i(\bar d) = S_i(a,\gamma)$. Define:
\[V_i(a,\gamma) \equiv \|(1, H_i')\|^{-1} [\ (a_1-\gamma a_2)\widehat Y_{i1} + a_2 \widehat Y_{i2}\ ], \]
where $\widehat Y_{it} \equiv Y_{it} - W_{it}'\widehat\delta_t$ for $t=1,2$.
Estimate the density of $(C_i(a,\gamma), S_i(a,\gamma))$ by replacing $V_i$ with $V_i(a,\gamma)$ in (\ref{eq:est_cstar}), and denote that estimator by $\widehat f_{C_i(a,\gamma), S_i(a,\gamma)}$.
For simplicity, we suppress the arguments $(a,\gamma)$ for the rest of this subsection.
Estimate the density of $(C_i, \gamma_i, \beta_i)$ at $(c,\gamma,b)$ by
\begin{align*}
\widehat f_{C_i, \gamma_i, \beta_{i}} (c,\gamma,b) \equiv \widehat f_{C_i, S_{i1}, S_{i2}, S_{i3}}(c,b_1a_1,b_2a_2, \widehat\delta_1a_2\gamma)|\widehat\delta_1a_1^{J_x}a_2^{J_x+1}|,
\end{align*}
and estimate the conditional density of $C_i$ given $\gamma_i = \gamma, \beta_i=b$ as
\begin{align*}
\widehat f_{C_i\mid \gamma_i = \gamma, \beta_i=b}(c) = \frac{\widehat f_{C_i, \gamma_i, \beta_{i}} (c,\gamma,b)}{\widehat f_{\gamma_i,\beta_i}(\gamma,b)}.
\end{align*}
For any $a\in\mathbb R^2$, let $\widetilde U_i(a) \equiv C_i - \psi^* = a_1U_{i1}+a_2U_{i2}$, where $\psi^*=a_1\gamma y_0$ as defined in (\ref{eq:iden_tU}).
Estimate the conditional density of $\widetilde U_i(a)$ by
\begin{align*}
\widehat f_{\widetilde U_i(a) \mid \gamma_i = \gamma, \beta_i=b}(\tilde u)
= \widehat f_{C_i \mid \gamma_i = \gamma, \beta_i=b}(\tilde u+a_1\gamma y_0).
\end{align*}
Estimate the joint density of $U_i=(U_{i1},U_{i2})'$ conditional on $\gamma_i = \gamma, \beta_i=b$ following steps similar to those in \cite{masten2018random}.
First, estimate the conditional characteristic function of $U_i$ by
\begin{equation*}
\widehat \phi_{U_{i1},U_{i2}\mid \gamma_i = \gamma, \beta_i=b}(a_1,a_2)=\int_{\mathcal U_n} \exp(i\tilde u)\, \widehat f_{\widetilde U_i(a) \mid \gamma_i = \gamma, \beta_i=b}(\tilde u)\, d\tilde u,
\end{equation*}
where $\mathcal U_n$ is the interval implied by $\bar c_n$, which expands as $n\to\infty$; see the proof of \Cref{theorem:u} in Appendix~\ref{sec:AppendixA}.
Then, use the inverse Fourier transform to estimate the conditional density of $U_i$ as
\begin{align*}
&\widehat f_{U_{i1},U_{i2}\mid \gamma_i= \gamma, \beta_i=b}(u_1,u_2) \\
= &\mathrm{Re}\!\left(\frac{1}{(2\pi)^2} \int_{\mathcal{A}_n}\!\exp[-i(a_1u_1+a_2u_2)]\widehat \phi_{U_{i1},U_{i2}\mid \gamma_i = \gamma, \beta_i=b}(a_1,a_2)da_1da_2\right),
\end{align*}
where $\mathrm{Re}(\cdot)$ is the real part of a complex number and
\[\mathcal{A}_n\equiv\{a\in\mathbb R^2:\ \underline{\vartheta}_n\le|a_1|\le \bar{\vartheta}_n,\ \underline{\vartheta}_n\le|a_2|\le \bar{\vartheta}_n\},\]
with tuning sequences $\bar{\vartheta}_n\to\infty$ and $\underline{\vartheta}_n\to0$ as $n\to\infty$.
Truncating the Fourier inversion to a bounded region is standard in deconvolution \citep{fan1991optimal,meister2009deconvolution}; \citet{masten2018random} applied the same technique in estimation.
\section{Asymptotic Properties} \label{sec:asymp}
Define a Sobolev space \[\mathcal{W}^m(\mathbb{R}^l) \equiv \{f \in L^2(\mathbb{R}^l): (1+\|\cdot \|^2)^{m/2}\mathcal{F}f(\cdot) \in L^2(\mathbb{R}^l)\}\] with a semi-norm
\[\|f\|^2_m \equiv \int_{\mathbb R^l}(1+\|\xi\|^2)^m | \mathcal F f(\xi)|^2d\xi,\]
where $\mathcal{F}f(\xi)\equiv \int_{\mathbb R^l}f(x) \exp(-ix'\xi)dx$ is the Fourier transform of $f(\cdot)$, and $m>0$ is the {\it order of smoothness} of the Sobolev space.
For any constant $B>0$, define
\[\mathcal{W}^m(\mathbb{R}^l; B) \equiv \{ f \in \mathcal{W}^m(\mathbb{R}^l): \| f\|_m \leq B\}.\]
We treat $\delta_t$ as known.
The first-step estimators $\widehat \delta_1,\widehat \delta_2$ in Section \ref{subsec:est-delta} are $\sqrt n$-consistent (with $Y_{i0}$ fixed at a constant $y_0$ in the sample).
With $m > l/2$, the uniform convergence rate of the intermediate density estimator is strictly slower than $n^{-1/2}$ (\Cref{lemma:rte}). The plug-in error from replacing $\delta_t$ by $\widehat\delta_t$ in $\widetilde Y_{it}$ is
asymptotically negligible relative to the inverse Radon transform estimation.
For any $d = (d_1,d_2)\in\mathbb R^2$, define a Radon transform estimator \(\widehat f_{C_i(d),S_i(d)}\) for the joint density of $(C_i(d),S_i(d))$ by replacing $v_i$ in (\ref{eq:est_cstar}) with
$v_i(d) \equiv \| (1,h_i')\|^{-1} (d_1 \tilde y_{i1} + d_2 \tilde y_{i2})$, where $\tilde y_{it} \equiv y_{it} - w_{it}'\delta_t$.
We establish uniform convergence of our estimators for the joint density of
$(\gamma_i,\beta_i)$ and the conditional density of $U_i$ given $(\gamma_i,\beta_i)$.
Recall these are all conditional on $Y_{i0}$, which is fixed at a constant value $y_0$ in the sample, and suppressed in notation.
We maintain the conditions for identification in \Cref{theorem:two} throughout, and do not reiterate them below for brevity.
In addition, we maintain the following assumptions.
\medskip
\begin{assumption}
\label{assn:rc}
Let $\mathcal T$ be a compact subset of the joint support of $(\gamma_i,\beta_i)$ such that
\begin{enumerate}[(i)]
\item \label{assn:rc:jac} the Jacobian determinant $|J_i(d)|$ in \eqref{defn:Jacobian} is bounded away from zero for any fixed $d \in \mathbb R^2$ with $d_2\neq0$.
\item \label{assn:rc:bound} $\sup_{(\gamma,b)\in\mathcal T}f_{\gamma_i,\beta_i}(\gamma,b)<\infty$.
\item \label{assn:rc:positive}
there exists a constant $M_b>0$ such that $\inf_{(\gamma,b)\in\mathcal T} f_{\gamma_i,\beta_i}(\gamma,b)\ge M_b$.
\end{enumerate}
\end{assumption}
\medskip
For a fixed $d$, the Jacobian $|J_i(d)|=|\delta_{1}d_2^{J_x+1}(d_1+d_2\gamma)^{J_x}|$ vanishes at $\gamma=-d_1/d_2$.
\Cref{assn:rc}\eqref{assn:rc:jac} excludes this value from the compact set $\mathcal T$.
\Cref{assn:rc}\eqref{assn:rc:positive} bounds $f_{\gamma_i,\beta_i}$ away from zero on $\mathcal T$, so that $\widehat f_{\gamma_i,\beta_i}$ stays away from zero asymptotically.
\medskip
\begin{assumption}
\label{assn:f_cs}
For any $d \in \mathbb R^2$ with $d_2\neq0$, the joint density $f_{C_i(d),S_i(d)}$ satisfies
\begin{enumerate}[(i)]
\item \label{assn:f_cs:smooth} $f_{C_i(d),S_i(d)} \in \mathcal W^m(\mathbb R^l; B)$ for some $B > 0$ and smoothness order $m > l/2$, and
\item \label{assn:f_cs:tail} there exist constants $M_c>0$ and $\kappa>l-1$ such that $f_{C_i(d),S_i(d)}(z)\le M_c\,(1+\|z\|)^{-\kappa}$ for all $z\in\mathbb R^l$.
\end{enumerate}
\end{assumption}
\medskip
The Sobolev smoothness in \Cref{assn:f_cs}\eqref{assn:f_cs:smooth} controls the bias of the Radon transform
estimator, and $m>l/2$ ensures that $f_{C_i(d),S_i(d)}$ is bounded and continuous.
This coincides with the condition of \cite{bissantz2014confidence} for the uniform bound on the bias.
The decay exponent $\kappa>l-1$
controls the variance of the Radon transform estimator.
\citet{holzmann2020} impose a similar polynomial tail bound that determines the optimal rate of their Fourier-based estimator.
\citet{goldenshluger2021deconvolution} impose a comparable tail bound on the target density.
\medskip
\begin{assumption}
\label{assn:f_Q}\leavevmode
\begin{enumerate}[(i)]
\item \label{assn:f_Q:bound} The density $f_{Q_i}$ is uniformly bounded.
\item \label{assn:f_Q:smooth} For some even integer $\zeta>(l-1)^2/(2m-l+1-2\varsigma)$ with some constant $\varsigma\in(0,m-l/2)$, the partial derivatives of $f_{Q_i}$ up to order $\zeta$ exist and are bounded.
\item \label{assn:f_Q:positive} There exist constants $M_q>0$ and $k >1$ such that $f_{Q_i}(q)\ge M_q\,q_1^{\,k}$ for any $q\in\mathbb{S}_+^{l-1}$.
\item \label{assn:f_Q:bw} The smoothing parameter of $\widehat f_{Q_i}$ satisfies $\tau\asymp(\ln n/n)^{1/(2\zeta+l-1)}$.
\item \label{assn:f_Q:sphK} The spherical kernel $\widetilde K$ is bounded, H\"older continuous, and of order at least $\zeta$.
\end{enumerate}
\end{assumption}
\medskip
\Cref{assn:f_Q}\eqref{assn:f_Q:bound} corresponds to the upper boundedness condition in Assumption~3 of \citet{hoderlein2010analyzing}; \eqref{assn:f_Q:positive} adapts the lower bound in their Assumption~4 to the trimmed setting, which yields \(M_{q,n}\equiv\inf_{q\in\mathbb S(\underline q_n)}f_{Q_i}(q)\ge M_q\underline q_n^k\).
The smoothness condition~\eqref{assn:f_Q:smooth}, together with the spherical kernel condition~\eqref{assn:f_Q:sphK}, gives the uniform rate of the spherical kernel estimator $\widehat f_{Q_i}$
under \(\tau\asymp(\ln n/n)^{1/(2\zeta+l-1)}\); see \citet[Appendix~B]{hoderlein2010analyzing}.
The next assumption is on the kernel function $K_\nu(\cdot)=\nu^{-l}K(\cdot/\nu)$ in Footnote~\ref{note:defn_K}.
\medskip
\begin{assumption}
\label{assn:K_nu}\leavevmode
\begin{enumerate}[(i)]
\item \label{assn:K_nu:bound} The kernel $K$ is uniformly bounded.
\item \label{assn:K_nu:lip} There exists a constant $M_k > 0$ such that, for all $z,z' \in \mathbb{R}$,
$|K(z)- K(z')| \leq M_k\|z- z'\|$.
\item \label{assn:K_nu:bw} The smoothing parameter $\nu$ is of the order $\nu\asymp(\underline q_n^{1-k}\ln n/n)^{1/(2m+l-1-2\varsigma)}$.
\item \label{assn:K_nu:order} The order parameter $r$ of the kernel $K_\nu$ satisfies $r\ge m-l/2$.
\end{enumerate}
\end{assumption}
\medskip
\Cref{assn:K_nu}\eqref{assn:K_nu:bound} and \eqref{assn:K_nu:lip} require the kernel $K$ to be bounded and Lipschitz continuous. They are standard for uniform convergence of kernel density estimators, as in \cite{li2007nonparametric} and \cite{masten2018random}.
\Cref{assn:K_nu}\eqref{assn:K_nu:bw} and \eqref{assn:K_nu:order} specify the rate of the bandwidth \(\nu\) and the order of the kernel.
Following Lemma~1 of \cite{bissantz2014confidence},
the order $r\ge m-l/2$ ensures the kernel does not limit the bias rate of the intermediate Radon transform estimator $f_{C_i(d),S_i(d)}$, which is then determined by the smoothness $m$.
\medskip
\begin{assumption}
\label{assn:f_u}
For any $(\gamma,b)$ on the joint support of $(\gamma_i,\beta_i)$,
\begin{enumerate}[(i)]
\item \label{assn:f_u:tail} there exists $p>1$ such that $\sup_{(\gamma,b)\in\mathcal T} E(|U_{ij}|^{p}\mid\gamma_i=\gamma,\beta_i=b)<\infty$ for $j=1,2$.
\item \label{assn:f_u:bound} $\max \{\int\sup_{u_2 \in \mathbb R}f_{U_i\mid\gamma_i=\gamma,\beta_i=b}(u_1,u_2)\,du_1, \int\sup_{u_1\in \mathbb R}f_{U_{i}\mid\gamma_i=\gamma,\beta_i=b}(u_1,u_2)\,du_2\}<\infty$.
\item \label{assn:f_u:smooth} The conditional density $f_{U_{i}\mid\gamma_i=\gamma,\beta_i=b}$ is continuously differentiable, and both
\(\int\sup_{u_1} (1+|u_2|)\,|\partial_{u_1} f_{U_{i}\mid\gamma,b}(u_1,u_2)|\,du_2\) and \(\int\sup_{u_2}(1+|u_1|)\,|\partial_{u_2} f_{U_{i}\mid\gamma,b}(u_1,u_2)|\,du_1\) are finite.
\item \label{assn:f_u:char} There exist constants $M_{\phi}>0$ and $\rho>2$ such that $|\phi_{U_{i}\mid\gamma_i=\gamma,\beta_i=b}(a)|\le M_{\phi}(1+\|a\|)^{-\rho}$ for all $a\in\mathbb R^2$, where $\phi_{U_i\vert \gamma_i,\beta_i}(\cdot)$ denotes the characteristic function.
\end{enumerate}
\end{assumption}
\medskip
\Cref{assn:f_u}\eqref{assn:f_u:tail} and~\eqref{assn:f_u:char} bound the conditional moments of the random intercepts and the tail of their characteristic function.
Conditions (\ref{assn:f_u:tail}), (\ref{assn:f_u:bound}), and (\ref{assn:f_u:char}) are instrumental for dealing with the trimming devices $\mathcal C_n $, $\mathcal U_n$, and $\mathcal A_n $ in asymptotics.
The smoothness condition~(\ref{assn:f_u:smooth}) adapts Assumption~E1.3 of \citet{masten2018random} to non-compact support, yielding a bound on $\|\nabla_a f_{\widetilde U_i(a)\mid\gamma,b}\|$ that is analogous to his Lemma~S3 under compact support.
\Cref{assn:f_u}\eqref{assn:f_u:char} is the ordinary-smooth condition of the deconvolution literature \citep{fan1991optimal,meister2009deconvolution}. Since $(1+\|a\|)^{-\rho}$ is integrable on $\mathbb R^2$ for $\rho>2$, the conditional characteristic function is absolutely integrable, which is Assumption~E1.4 of \citet{masten2018random} for the Fourier inversion. The decay helps to account for the bias from truncating the inversion on $\mathcal A_n$.
Define $r_{1n}\equiv(\underline q_n^{1-k}\ln n/n)^{\frac{m-l/2-\varsigma}{2m+l-1-2\varsigma}}$, $r_{2n}\equiv(\ln n/n)^{\frac{\zeta}{2\zeta+l-1}}$, and let $$r_n\equiv r_{1n}+\underline q_n^{1/2} +\underline q_n^{1-k}r_{2n}r_{1n}^{-\frac{2(l-1)}{2m-l-2\varsigma}}.$$
\begin{assumption}\label{assn:tuning}
$\underline q_n\to0$,
$\bar c_n \asymp r_n^{-1/(p+1)}$, and
$$ \qquad \underline q_n^{-1}(\ln n/n)^{\min\left\{\frac{2m+l-2-2\varsigma}{k(4m+2l-3-4\varsigma)+1},\ \frac{\zeta}{k(2\zeta+l-1)}, \ \frac{\zeta(2m-l+1-2\varsigma)-(l-1)^2}{2(k-1)(m-\varsigma)(2\zeta+l-1)}\right\}}\to 0.$$
\end{assumption}
The sequence $r_n$ characterizes the uniform convergence rate of an intermediate density estimator $\widehat f_{C_i(d),S_i(d)}$ for a fixed index $d$ (\Cref{lemma:rte}).
\Cref{assn:tuning} restricts the trimming sequence $\underline q_n$ and the truncation sequence $\bar c_n$, and ensures $r_n \to 0$.
The next proposition establishes uniform convergence of our estimator for the density of $(\gamma_i,\beta_i)$.
Recall that $l = 2J_x+2$ denotes the dimension of $(1,H_i')'$ (with $W_{i1}\in\mathbb R$).
\medskip
\begin{proposition}\label{prop:margin}
Fix $d\in\mathbb R^2$ with $d_2\neq 0$, and suppose Assumptions \ref{assn:rc}\eqref{assn:rc:jac}--\eqref{assn:rc:bound}, \ref{assn:f_cs}, \ref{assn:f_Q}, \ref{assn:K_nu}, \ref{assn:f_u}\eqref{assn:f_u:tail} and~\ref{assn:tuning} hold.
Then, for any moment order $p$ in \Cref{assn:f_u}\eqref{assn:f_u:tail},
$$\sup_{(\gamma,b) \in \mathcal T} \left|\widehat f_{\gamma_i,\beta_i}(\gamma,b) - f_{\gamma_i,\beta_i}(\gamma,b)\right| = O_p\!\left(r_n^{p/(p+1)}\right).$$
\end{proposition}
\begin{remark}
The adjustment factor $p/(p+1)$ is due to the accommodation of non-compact support of random intercepts $U_i$. The convergence rate is $r_n$ under compact support, or arbitrarily close to $r_n$ if the conditional moments are finite for all orders $p> 1$ (e.g., a Gaussian, a Gaussian mixture, or a logistic distribution).
\end{remark}
\begin{remark}
\citet{hoderlein2010analyzing} estimate the joint density of random coefficients in a cross-sectional linear model using a Radon transform estimator, including a trimmed version for models with an intercept.
In our setting, the Radon transform estimator is applied to the intermediate reduced-form density \(f_{C_i(d),S_i(d)}\).
\Cref{prop:margin} provides a sup-norm rate required for the conditional density of $U_i$ in \Cref{theorem:u}.
\end{remark}
We next establish uniform convergence of our estimator for the
joint density of $U_i \equiv (U_{i1},U_{i2})'$ conditional on
$(\gamma_i, \beta_i)$, which we recover by Fourier inversion of
the characteristic function estimator
$\widehat \phi_{U_{i1},U_{i2} \mid \gamma_i,\beta_i}$.
The result is obtained under Assumptions S.1 and S.2 in Online Appendix S2.
Assumption S.1 imposes tail and smoothness bounds on the density of $(C_i(\bar d),S_i(\bar d))$ uniformly over $a \in \mathcal A_n$ almost everywhere.
Assumptions S.2 imposes restrictions on the tuning sequences $\bar \vartheta_n$ and $\underline \vartheta_n$.
\medskip
\begin{theorem}\label{theorem:u}
Suppose Assumptions~\ref{assn:rc}, \ref{assn:f_cs}, \ref{assn:f_Q}, \ref{assn:K_nu}, \ref{assn:f_u}, \ref{assn:tuning}, S.1, S.2 hold,
and that $E(|\widetilde Y_{i1}|)<\infty$ and $E(|\widetilde Y_{i2}|)<\infty$. Then, for any $(\gamma,b)$ on the support of $(\gamma_i,\beta_i)$,
\begin{align*}
\sup_{(u_1,u_2)\in \mathbb{R}^2}\left|\widehat f_{U_{i1},U_{i2}|\gamma_i = \gamma, \beta_i=b}(u_1,u_2) - f_{U_{i1},U_{i2}|\gamma_i = \gamma, \beta_i=b}(u_1,u_2)\right| = o_p(1).
\end{align*}
\end{theorem}
\begin{remark}
\citet[Section~C.4]{masten2018random} proves sup-norm consistency of this Fourier inversion under compact support;
\Cref{theorem:u} permits non-compact support.
The proof bounds the error of the Fourier inversion of $\widehat\phi$ over $\mathbb R^2$ by decomposing it into $R_1+R_2+R_3$, where
\begin{equation*}
R_1=\frac{1}{(2\pi)^2}\int_{\mathcal A_n}|\widehat\phi(a)-\phi(a)|\,da,\qquad
R_2+R_3=\frac{1}{(2\pi)^2}\int_{\mathbb R^2\setminus\mathcal A_n}|\phi(a)|\,da,
\end{equation*}
where $R_2$ and $R_3$ decompose further over $[-\bar{\vartheta}_n,\bar{\vartheta}_n]^2\setminus\mathcal A_n$ and $\mathbb R^2\setminus[-\bar{\vartheta}_n,\bar{\vartheta}_n]^2$, respectively.
We derive the convergence rates of these components.
The estimation error $R_1$ is due to three sources: the inverse Radon transform estimation, the trimming over $\mathcal A_n$, and the truncation of the non-compact intercept support.
The truncation errors satisfy $R_2 + R_3 = O(\underline{\vartheta}_n+\bar{\vartheta}_n^{-(\rho-2)})$, which follows from the characteristic function decay in \Cref{assn:f_u}\eqref{assn:f_u:char}.
\Cref{theorem:u} then follows from Assumption~S.2.
\end{remark}
\section{Monte Carlo Simulations}\label{sec:simulations}
In this section, we evaluate the finite-sample performance of the multi-step estimators for the joint density of $(\gamma_i,\beta_i)$ and the conditional density of $U_i \equiv (U_{i1},U_{i2})'$ given $(\gamma_i,\beta_i)$. We focus on the case with $T=2$ and fix the initial condition at $y_0=1$ throughout this simulation study.
The data-generating process (DGP) follows \eqref{eq:seq-exo-w}, with $(\delta_1,\delta_2)=(1,1)$. The regressors are generated independently of $(\gamma_i,\beta_i,U_i)$ as follows:
\[
X_i \equiv (X_{i1},X_{i2})'\sim \mathcal{N}(0,\Sigma_x),
\qquad
\Sigma_x=
\begin{pmatrix}
2 & 1 \\
1 & 2
\end{pmatrix},
\]
and
\(
W_{i2}=0.5\,W_{i1}+0.3\,W^2_{i1}+\varepsilon_{w,i},
\)
where $W_{i1}\sim\mathcal{N}(0,2.5)$ and $\varepsilon_{w,i}\sim\mathcal{N}(0,2.5)$.
The autoregressive coefficient is generated as $\gamma_i \sim \text{Beta}(6,3).$
The marginal distributions of $\beta_i$ and $U_i$ are given by
\begin{equation}\label{eq:marginals}
\beta_i \sim \mathcal{N}(\mu_\beta,\Sigma_\beta),
\qquad
U_i \sim \mathcal{N}(0,\Sigma_u),
\end{equation}
where
\[
\mu_\beta=(1,1)', \qquad
\Sigma_\beta=
\begin{pmatrix}
0.5 & 0.25 \\
0.25 & 0.5
\end{pmatrix},
\qquad
\Sigma_u=
\begin{pmatrix}
0.3 & 0.15 \\
0.15 & 0.3
\end{pmatrix}.
\]
\noindent We consider the following four designs.
\paragraph{Baseline.} The random coefficients and intercepts $(\gamma_i, \beta_i, U_i)$ are mutually independent, with marginals as specified above.
\paragraph{Correlated.} This design introduces dependence among $(\gamma_i, \beta_i, U_i)$ while preserving the conditional mean restriction in Assumption~\ref{assn:seq-exo}.
While $\gamma_i$ has the same Beta marginal as in the baseline design, $\beta_i$ and $U_i$ are constructed via
$$
\beta_i = \mu_\beta + \exp(\lambda_\beta \gamma_i)\, L_\beta\, \eta_i^\beta, \qquad U_i = \exp(\lambda_u \gamma_i)\, L_u\, \eta_i^u,
$$
where $L_\beta$ and $L_u$ are the lower-triangular matrices from the Cholesky decompositions of $\Sigma_\beta$ and $\Sigma_u$, respectively.
The shocks $\eta_i = (\eta_i^\beta, \eta_i^u)'$ are drawn
independently of $\gamma_i$ from $\mathcal{N}(0, \Sigma_{\mathrm{corr}})$, and the $4$--by--$4$ correlation matrix $\Sigma_{\mathrm{corr}}$ has identity diagonal blocks and cross-block matrix
\[
\Sigma^{\beta,u} = \begin{pmatrix} 0.25 & 0.10 \\ 0.10 & 0.25 \end{pmatrix}.\]
We set $\lambda_\beta=\lambda_u=0.3$.
Under this design, conditional on $\gamma_i$, the vectors $\beta_i$ and $U_i$ have the same means as in the baseline design but covariance matrices that scale with $\gamma_i$:
$$\beta_i \mid \gamma_i \sim \mathcal{N}\left(\mu_\beta,\,
\exp(2\lambda_\beta\gamma_i)\, \Sigma_\beta\right), \qquad
U_i \mid \gamma_i \sim \mathcal{N}\left(0,\,
\exp(2\lambda_u\gamma_i)\, \Sigma_u\right).$$
\paragraph{Scale-dependence.}
This design is simpler than the preceding correlated design in that there is no direct correlation between $\beta_i$ and $U_i$; their dependence is solely driven by $\gamma_i$ through the scale:
\[
\beta_i \mid \gamma_i \sim \mathcal{N}\left(\mu_\beta,\, v(\gamma_i)\,\Sigma_\beta\right),
\qquad
U_i \mid \gamma_i \sim \mathcal{N}\left(0,\, w(\gamma_i)\,\Sigma_u\right),
\]
where
\(
v(\gamma)=\exp\left[2\lambda_\beta(\gamma-\mathbb{E}[\gamma])\right]/Z(\lambda_\beta)\) and
\(w(\gamma)=\exp\left[2\lambda_u(\gamma-\mathbb{E}[\gamma])\right]/Z(\lambda_u),
\)
with $Z(\cdot)$ being a normalization constant such that $\mathbb{E}[v(\gamma_i)]=\mathbb{E}[w(\gamma_i)]=1$.
\paragraph{Bimodal.} The random coefficients $\gamma_i$ and $\beta_i$ are independent with the same marginal distributions as in the baseline design, but $U_i$ follows a symmetric Gaussian mixture:
\[
U_i \sim \tfrac12\,\mathcal{N}(\mu_u,\Sigma_u)
+ \tfrac12\,\mathcal{N}(-\mu_u,\Sigma_u),
\qquad \mu_u=(0.7,-0.7)'.
\]
\citet{hoderlein2010analyzing} and \citet{masten2018random} use similar Gaussian-mixture designs to investigate estimator performance for multi-modal densities.
For each design, we vary the sample sizes $n \in \{500, 2000, 8000\}$ and generate $N_{\mathrm{rep}} = 500$ Monte Carlo replications.
The estimators for the joint density $\widehat f_{\gamma_i,\beta_i}$ and for the conditional density $ \widehat f_{U_{i1},U_{i2} \mid \gamma_i=\gamma,\beta_i=b}$ are defined as in Sections~\ref{subsec:est_RC} and~\ref{subsec:est_U}.
Following \citet{hoderlein2010analyzing}, we select the bandwidths $(\tau, \nu)$ by minimizing the mean density-weighted integrated squared error.
Following \citet{masten2018random}, we obtain an optimal $\nu$ for a single $ a \in [-\bar{\vartheta}_n,\bar{\vartheta}_n]^2$ and then rescale it to the others. In the baseline design at $n=8000$, the average $\nu$ used over all grid points $a$ is $0.65$.
As the asymptotic theory prescribes, the cutoff $\bar{\vartheta}_n$ increases with the sample size, and we scale the bandwidths for the other sample sizes accordingly.
We first present the performance of $\widehat f_{\gamma_i, \beta_i}$ through one-dimensional conditional densities. Table~\ref{tab:jgb_mise} reports the mean integrated squared error (MISE) of $\widehat f_{\gamma_i \mid \beta_i}$, $\widehat f_{\beta_{i1} \mid \gamma_i, \beta_{i2}}$, and $\widehat f_{\beta_{i2} \mid \gamma_i, \beta_{i1}}$.
For every design and conditional density reported, the MISE declines with $n$, confirming the consistency of our estimator.
\begin{table}[H]
\begin{center}
\caption{Conditional densities of the estimated joint density of $(\gamma_i, \beta_i)$}
\label{tab:jgb_mise}
\scalebox{1}{
\begin{tabular}{ll ccc}
\toprule
& & \multicolumn{3}{c}{MISE} \\
\cmidrule(lr){3-5}
& $n$ & $\widehat f_{\gamma_i \mid \beta_i}$ & $\widehat f_{\beta_{i1} \mid \gamma_i, \beta_{i2}}$ & $\widehat f_{\beta_{i2} \mid \gamma_i, \beta_{i1}}$ \\
\midrule
\multirow{3}{*}{Baseline}
& 500 & 0.2359 & 0.1054 & 0.0654 \\
& 2000 & 0.2001 & 0.0790 & 0.0476 \\
& 8000 & 0.1749 & 0.0655 & 0.0369 \\
\midrule
\multirow{3}{*}{Correlated}
& 500 & 0.2699 & 0.0664 & 0.0414 \\
& 2000 & 0.2372 & 0.0468 & 0.0306 \\
& 8000 & 0.2052 & 0.0374 & 0.0227 \\
\midrule
\multirow{3}{*}{Scale}
& 500 & 0.2417 & 0.1031 & 0.0607 \\
& 2000 & 0.2046 & 0.0792 & 0.0459 \\
& 8000 & 0.1747 & 0.0644 & 0.0355 \\
\midrule
\multirow{3}{*}{Bimodal}
& 500 & 0.2434 & 0.1066 & 0.0639 \\
& 2000 & 0.2096 & 0.0816 & 0.0483 \\
& 8000 & 0.1860 & 0.0695 & 0.0395 \\
\bottomrule
\end{tabular}}\end{center}
\vspace{12pt}
{\noindent\footnotesize \textit{Notes:} The table reports, across $N_{\mathrm{rep}} = 500$ Monte Carlo replications, the mean integrated squared error (MISE) of the estimated joint density $\widehat f_{\gamma_i, \beta_i}$, summarized through its three one-dimensional conditional densities; each varies one coordinate of $(\gamma_i, \beta_{i1}, \beta_{i2})$ with the other two held at their medians.}
\end{table}
Figures~\ref{fig:marginal_Q50} and~\ref{fig:heatmap_Q50} show the estimated conditional density of $U_i$ along each axis and as contour plots, and Table~\ref{tab:cond_u} reports its MISE and the bias, standard deviation (SD), and mean squared error (MSE) of the estimated interquartile range (IQR) of its two marginals.
For the three unimodal designs (Baseline, Correlated, and Scale), the MISE falls with $n$ and the Monte Carlo mean density converges toward the true density, with the bias small in large samples.
The conditional location is recovered accurately: along each axis the estimated density stays centered on the true density even in small samples.
\citet[Figure~2 and Table~2]{masten2018random} likewise finds the location well recovered even when the full shape is hard to estimate.
The IQR in Table~\ref{tab:cond_u} is estimated less precisely in small samples but approaches its true value as $n$ grows. For the three unimodal designs it tends to be over-estimated, possibly due to over-smoothing at small $n$ which blurs the estimate: it moves probability mass away from the peak into the sides, so the estimated density is too flat and too wide, and its IQR too large (Figure~\ref{fig:marginal_Q50}).
\citet[Table~2]{masten2018random} reports the same upward IQR bias for his unimodal designs, where the estimated density is more spread out than the true density.
The Bimodal design is the most demanding. The MISE declines slowly ($0.054$ to $0.035$). Along each axis the estimate for the marginal density does not register a clear bimodal pattern even at $n=8000$ (Figure~\ref{fig:marginal_Q50}), although the contour separates its two modes by $n=2000$ (Figure~\ref{fig:heatmap_Q50}). Relatedly, the spread and hence the IQR are \emph{under}-estimated. Again, this pattern is consistent with \citet{masten2018random}'s observation that the shape is harder to estimate for multi-modal densities.
\begin{table}[ht]
\begin{center}
\caption{Estimator for the conditional density of $U_i$}
\label{tab:cond_u}
\scalebox{1}{
\begin{tabular}{ll c ccc ccc}
\toprule
& & & \multicolumn{3}{c}{$\widehat{\mathrm{IQR}}(U_{i1} \mid \gamma_i, \beta_i)$} & \multicolumn{3}{c}{$\widehat{\mathrm{IQR}}(U_{i2} \mid \gamma_i, \beta_i)$} \\
\cmidrule(lr){4-6} \cmidrule(lr){7-9}
& $n$ & MISE & Bias & SD & MSE & Bias & SD & MSE \\
\midrule
\multirow{3}{*}{Baseline}
& 500 & 0.0841 & 0.4571 & 0.0306 & 0.2099 & 0.4353 & 0.0294 & 0.1903 \\
& 2000 & 0.0277 & 0.2313 & 0.0171 & 0.0538 & 0.2157 & 0.0172 & 0.0468 \\
& 8000 & 0.0053 & 0.1039 & 0.0099 & 0.0109 & 0.0903 & 0.0099 & 0.0082 \\
\midrule
\multirow{3}{*}{Correlated}
& 500 & 0.0516 & 0.4449 & 0.0293 & 0.1988 & 0.4256 & 0.0297 & 0.1820 \\
& 2000 & 0.0163 & 0.2117 & 0.0171 & 0.0451 & 0.1920 & 0.0180 & 0.0372 \\
& 8000 & 0.0044 & 0.0602 & 0.0125 & 0.0038 & 0.0435 & 0.0137 & 0.0021 \\
\midrule
\multirow{3}{*}{Scale}
& 500 & 0.0833 & 0.4552 & 0.0311 & 0.2082 & 0.4331 & 0.0292 & 0.1884 \\
& 2000 & 0.0271 & 0.2297 & 0.0182 & 0.0531 & 0.2130 & 0.0179 & 0.0457 \\
& 8000 & 0.0055 & 0.1024 & 0.0097 & 0.0106 & 0.0871 & 0.0107 & 0.0077 \\
\midrule
\multirow{3}{*}{Bimodal}
& 500 & 0.0538 & $-$0.2158 & 0.0738 & 0.0520 & $-$0.2236 & 0.0775 & 0.0560 \\
& 2000 & 0.0421 & $-$0.2136 & 0.0626 & 0.0496 & $-$0.2129 & 0.0734 & 0.0507 \\
& 8000 & 0.0346 & $-$0.1057 & 0.0489 & 0.0136 & $-$0.1246 & 0.0594 & 0.0191 \\
\bottomrule
\end{tabular}}\end{center}
\vspace{12pt}
{\noindent\footnotesize \textit{Notes:} The table reports, across $N_{\mathrm{rep}} = 500$ Monte Carlo replications, the mean integrated squared error (MISE) of the estimated conditional density $\widehat f_{U_{i1},U_{i2}\mid \gamma_i, \beta_i}$; and the bias, standard deviation (SD), and mean squared error (MSE) of the estimated interquartile range $\widehat{\mathrm{IQR}}(U_{ij} \mid \gamma_i, \beta_i)$, $j=1,2$, of the two marginal distributions. All quantities are evaluated at the median value of $(\gamma_i, \beta_i)$.}
\end{table}
\begin{figure}[p]
\begin{center}
\scalebox{0.9}{
\includegraphics[width=\textwidth,height=0.92\textheight,keepaspectratio]{figures/marginal_slices_v2_runs_routeT4_N8000_T4.25_Q50.png}}
\caption{Estimated conditional density of $U_i$ along each axis}
\label{fig:marginal_Q50}
\end{center}
\vspace{12pt}
{\noindent\footnotesize \textit{Notes:} $\widehat f_{U_{i1}, U_{i2} \mid \gamma_i, \beta_i}$ along $u_1$ (at $u_2=0$) and along $u_2$ (at $u_1=0$), at the median value of $(\gamma_i,\beta_i)$. Each row shows these two curves for one design; columns are the sample sizes $n \in \{500, 2000, 8000\}$. Black dashed: the true density; navy solid: the Monte Carlo mean across $N_{\mathrm{rep}} = 500$ replications; shaded band: the 95\% pointwise Monte Carlo confidence band.}
\end{figure}
\begin{figure}[p]
\begin{center}
\scalebox{0.9}{
\includegraphics[width=\textwidth,height=0.92\textheight,keepaspectratio]{figures/contour_overlay_runs_routeT4_N8000_T4.25_Q50.png}}
\caption{Contour plots of the estimated conditional density of $U_i$}
\label{fig:heatmap_Q50}
\end{center}
\vspace{12pt}
{\noindent\footnotesize \textit{Notes:} $\widehat f_{U_{i1}, U_{i2} \mid \gamma_i, \beta_i}$ for the four designs (rows) at sample sizes $n \in \{500, 2000, 8000\}$ (columns), conditioning on the median value of $(\gamma_i, \beta_i)$. Gray fill with black dashed contours: the true density; navy solid contours: the Monte Carlo mean across $N_{\mathrm{rep}} = 500$ replications. The contours of the true density enclose 25\%, 50\%, and 75\% of its probability mass, and the Monte Carlo contours use the same levels.}
\end{figure}
\FloatBarrier
\section{Conclusion}
We recover the entire joint distribution of the random coefficients and the time-varying errors. This is new relative to the literature on random-coefficient triangular and simultaneous equation models which left the distribution of time-varying errors unidentified. This is possible because a linear combination of the outcome history behaves as a single-equation random-coefficient regression in distinct regressors: inverting a Radon transform recovers a reduced-form density, and varying the combination and using deconvolution recover the distribution of the time-varying errors. We apply the analog principle to this constructive identification strategy, and propose a closed-form, multi-step estimator.
The identifying condition, a distributional form of strict exogeneity under which the covariates are independent of the coefficients and errors given the initial condition, is both the source of this reach and the main limitation of the approach: it permits arbitrary dependence among the coefficients and errors, but lets the covariates and the unobserved heterogeneity be related only through the initial condition.
Several questions remain for future research, including inference, data-driven bandwidth choices, and a continuously distributed initial condition.
\vspace{1cm}