Correlated Random Coefficient Distributions in Linear Panel Models

Correlated Random Coefficient Distributions in Linear Panel Models

\def\spacingset#1 \spacingset{1}

{\it Keywords:} irregular identification, correlated random coefficients, deconvolution.

\spacingset{1.9}

Introduction

Heterogeneity in individual behavior is pervasive in microeconometric applications, and panel data are a natural way to accommodate it. A leading framework is the correlated random coefficient (CRC) model (Chamberlain82), in which slope coefficients are heterogeneous across units and may covary with the regressors. In this setting, a natural target is the average partial effect (APE), given by the population mean of the random coefficient vector. Chamberlain1992 studied identification of the APE in the regular design, where the number of effective time periods $T$ exceeds the number of correlated coefficients $p$. GrahamPowell2012 showed that the APE remains identified in the irregular design, where $T=p$, by exploiting the distinct identifying content of movers and stayers. ArellanoBonhomme2012 focus on the regular design and show various moments of the correlated random coefficients, and even their density, can be identified under ARMA-type restrictions on the error process. In this paper, we derive sufficient conditions under which the density of the correlated random coefficients is identified in both regular and irregular designs, without restricting the time-series properties of the error terms.

We consider models with the following outcome equation. For $i=1,\dots,n$,

where the correlated random coefficients $\beta_i\in \mathbb R^p$ may depend arbitrarily on the observed regressors $(X_i,W_i)$, while the uncorrelated components $D_i\in\mathbb R^d$ are independent of $(\beta_i,X_i,W_i)$.\footnote{As explained in the next section, $D_i\in\mathbb R^d$ collects all remaining unobserved components, which includes both time-invariant random coefficients and time-varying idiosyncratic errors. The assumption can be weakened to conditional independence, e.g., \(\beta_i\perp D_i\mid (X_i,W_i)\) and \(D_i\perp X_i\mid W_i\). Then the results would be conditional on \(W_i\). We do not pursue this extension.} The coefficient $\beta_i$ represents the unit-level structural partial effect of $X_i$. Our parameter of interest is $f_\beta$, the cross-sectional density of $\beta_i$. Ours is a short-$T$ setting.

Identification proceeds in two steps. In the first step, we identify the characteristic function of $D_i$ by using a transformation, or annihilator matrix, that eliminates the $\beta_i$ component. The transformation depends on the relation between $T$ and $p$. In the regular design ($T > p$), there are multiple annihilation matrices; projection onto the orthogonal complement of the column space of $X_i$ achieves annihilation for all observations. In the irregular design ($T = p$), the annihilation matrix is the adjugate matrix of $X_i$; the $\beta_i$ component is then eliminated for stayers, i.e. units with $\det(X_i)=0$. The adjugate--determinant algebra used in the irregular design shares a common origin with GrahamPowell2012, but it is embedded here in an argument based on characteristic functions.

In the second step, we apply another transformation to the outcome equation that isolates $\beta_i$. We then recover $f_\beta$ by deconvolution, exploiting the independence restriction on $D_i$. Notably, we impose no restrictions on the time-series properties of the error terms.

We propose a sieve minimum-distance estimator for the density of the correlated coefficients and study its finite-sample properties via Monte Carlo simulations. The estimators we propose specialize to the scalar irregular design $(T,p)=(1,1)$ and to the low-dimensional regular design $(T,p)=(2,1)$, with $W_i$ being the identity matrix; these are the configurations that arise in our application. Our estimator is a two-step procedure. In the first step, which is design-specific, we construct a nonparametric estimator of the unconditional characteristic function of the correlated coefficients. This involves estimating a conditional characteristic-function ratio at each observation and then averaging over the sample to integrate out the covariates. In the irregular design, we partition the sample into stayers, used to estimate the characteristic function of $D_i$, and movers, used to estimate a trimmed ratio of characteristic functions, which is then averaged across the movers. In the regular design, no mover--stayer partition is needed: the characteristic function of $D_i$ is estimated by smoothing over normalized directions on the unit sphere, and the ratio is averaged over the full sample. In the second step, which is common across designs, we approximate the density by a finite-dimensional expansion in orthonormal Hermite functions and minimize a weighted distance between the Fourier transform of this expansion and the first-step characteristic function estimate, subject to a unit-mass normalization. The Hermite basis is a natural choice: Hermite functions are eigenfunctions of the Fourier transform, so the criterion is quadratic in the sieve coefficients and admits a closed-form constrained minimizer. We provide a heuristic rate decomposition that guides the choice of tuning parameters, a cross-validation procedure for selection of the bandwidth and sieve dimension, and Monte Carlo evidence on finite-sample performance.\footnote{Large-sample theory for the two-stage estimator, including data-driven tuning and bootstrap validity, is left for future work.}

We revisit the application in GrahamPowell2012, which uses panel data from rural Nicaragua collected as part of the evaluation of the conditional cash transfer program Red de Protección Social (RPS). In this application, the correlated random coefficient is the household-specific calorie--expenditure elasticity. We estimate its full cross-sectional density. The estimated densities reveal substantial heterogeneity: there is considerable dispersion, nontrivial mass near zero, and a non-negligible share of negative elasticities. These findings help reconcile mixed estimates in the calorie-demand literature, where average elasticities range from values close to zero BehrmanDeolalikar1987,BouisHaddad1992,Ravallion1990 to values around 0.3--0.5 SubramanianDeaton1996. They are also consistent with the view that households differ in how additional resources are allocated across quantity, quality, and nonfood margins Deaton1997,JensenMiller2008,Skoufias2011,DeatonDreze2009.

The experimental design of RPS gives the distributional estimates a causal interpretation. Since assignment was randomized at the village level and take-up was high, differences between the estimated elasticity distributions for recipient and nonrecipient households can be interpreted as causal contrasts between regime-specific densities, subject to the maintained CRC assumptions. The three survey waves also allow us to compare scalar irregular estimates separately for 2000--2001 and 2001--2002. Under the baseline model, which imposes a time-invariant household-specific elasticity, these estimates should recover the same density. The leftward shift we find among RPS recipients is therefore suggestive evidence that elasticities may vary with cumulative program exposure, while the regular stacked estimator should be viewed as a pooled benchmark under the common-\(\beta_i\) restriction.

Our paper contributes to the literature on CRC models in short panels. Existing work primarily focuses on low-dimensional functionals of the heterogeneity distribution. For example, GrahamPowell2012 study average partial effects in irregular designs, GrahamHahnPoirierPowell2018 study quantile effects under comonotonicity restrictions in both regular and irregular designs, while SasakiUra2026 study inference for average partial effects in irregular designs with slow movers. In regular designs, Verdier2020 considers identification of average treatment effects for stayers, recovering conditional means via linear extrapolation; Laage2024 studies average effects in the presence of time-varying endogeneity through control variables; while MurisWacker2026 analyze interaction effects, that is, how the correlated coefficient varies with observable regressors. We instead identify and estimate the marginal distribution of $\beta_i$ in both regular and irregular designs. This allows us to recover higher-order moments and shape features of the heterogeneity distribution, including dispersion, skewness, and multimodality, which are not accessible through mean- or quantile-based approaches.

The closest related result is ArellanoBonhomme2012, who also establish identification of the distribution of correlated random coefficients in regular designs. Their approach exploits restrictions on the time-series properties of the error terms. In contrast, we allow for unrestricted serial dependence in the error terms and instead exploit restrictions on the dependence between the errors and the observed regressors. This alternative source of variation delivers identification in both regular and irregular designs.

Building on ArellanoBonhomme2012, BotosaruLiu2025 uses empirical Bayes methods to recover posterior means of the correlated random coefficients. In contrast, we recover the marginal distribution of $\beta_i$ without imposing parametric restrictions, using a deconvolution-based estimator rather than empirical Bayes methods.

In a related but distinct framework, HoderleinWhite2012 study local average responses for stayers in nonseparable panel models under a time-invariant structural function. Chernozhukovetal2015 extend this to quantile effects for stayers. We maintain linearity of the outcome equation and identify the distribution of $\beta_i$ for the entire population.

Our estimation strategy builds on tools developed for random coefficient density estimation in cross-sectional settings. HoderleinKMammen propose a Radon-transform-based estimator, while Breunig2021 develops a sieve minimum-distance estimator based on Hermite functions. We also use a sieve-based weighted minimum-distance criterion, where the criterion function involves a ratio of characteristic functions evaluated at data-dependent frequencies.

More broadly, our approach is related to the deconvolution literature based on characteristic functions. In particular, KatoSasakiUra2021 reformulate Kotlarski's identity as a system of complex-valued moment restrictions and conduct inference via test inversion. This strategy yields inference that is robust to potential failure of the completeness condition. Our panel setting also implies moment restrictions derived from a characteristic-function factorization. Unlike their repeated-measurement model, which targets an unconditional characteristic function, our setting involves a conditional characteristic function indexed by the regressors. Moreover, the moments depend on a generated first-step estimator evaluated at data-dependent frequencies. These features, discussed in detail in Remark (ref), lead us to base estimation on an explicit regularized deconvolution step rather than moment inversion.

Finally, our results connect to applications with continuously measured treatments. Recent work studies average effects in such settings using difference-in-differences (DID), both when stayers are available deChaisemartinetal2024b and when all units change treatments between periods deChaisemartinetal2024a; see also CallawayGBSantAnna for a framework for continuous DID under parallel trends. Distributional effects have also been analyzed under unconfoundedness, e.g., GalvaoWang2015, CallawayHuang2020.

The remainder of the paper is organized as follows. Section (ref) introduces the model and identifying assumptions. Section (ref) presents the identification results for both regular and irregular designs, with worked examples. Section (ref) describes the sieve minimum-distance estimator. Section (ref) reports Monte Carlo evidence on the finite-sample performance of the estimator. Section (ref) applies the methods to estimate the distribution of calorie demand elasticities using the Nicaraguan household panel of Graham and Powell (2012).

Notation

For a matrix \(M\), \(\operatorname{adj}(M)\) denotes the adjugate matrix and \(\det(M)\) its determinant. For a vector $v$, $\|v\|$ denotes the Euclidean norm, for a matrix $M$, \(\|M\|\) denotes the Frobenius norm. The joint characteristic function of an arbitrary \(q\times 1\) random vector \(Z=(Z_1,\dots,Z_q)'\) is \[ \varphi_Z(s)\equiv \mathbb{E}\!\left[e^{\imath s'Z}\right],\qquad s\in\mathbb R^q, \qquad \imath=\sqrt{-1}. \] Note that we define characteristic functions on column vectors.

Whenever we condition on an event of probability zero, the conditional object is understood as a regular conditional distribution or conditional expectation whenever it exists. In particular, conditioning on the probability zero event \(\det(X_i)=0\) in the irregular design is interpreted through the limit as \(\det(X_i)\to 0\), under appropriate continuity assumptions stated in the paper.

For each $(s,u)\in\mathbb{S}^{d-1}\times\mathbb{R}$, let $\sigma_{s,u}$ denote the $(d-1)$-dimensional Hausdorff measure on the hyperplane $\{\delta\in\mathbb{R}^d:s'\delta=u\}$. Let $\mathcal{R}$ denote the Radon transform: \[ \mathcal{R}g(s,u) \equiv \int_{\{\delta:\,s'\delta=u\}} g(\delta)\,d\sigma_{s,u}(\delta), \] for any integrable $g:\mathbb{R}^d\to\mathbb{R}$.

Model

The outcome equation (ref) is our point of departure. We think of it as having been generated by, for example,

where $x_{it}\in\mathbb{R}^p$ and $w_{it}\in\mathbb{R}^q$ are observed covariates; $\beta_i\in\mathbb{R}^p$ and $\delta_i\in\mathbb{R}^q$ are vectors of unit-specific random coefficients; and $u_{it}$ is an idiosyncratic error. The random coefficients $(\alpha_i,\beta_i')$ may be arbitrarily correlated with the covariates $\{x_{it},w_{it}\}_{t=0}^{T}$. The random coefficients $\delta_i$ are independent of the covariates. Equation (ref) is obtained by first-differencing (ref) across adjacent periods to eliminate $\alpha_i$, stacking over the effective time periods $t=1,\dots,T$, and letting $Y_{it}\equiv y_{it}-y_{i,t-1}$, $X_{it}\equiv x_{it}-x_{i,t-1}$, $W_{it}\equiv w_{it}-w_{i,t-1}$, $U_{it}\equiv u_{it}-u_{i,t-1}$, and $$W_i \equiv \bigl(\,\Delta{w}_i \;\; I_T\bigr), \; \Delta{w}_i \equiv (W_{i1},\ldots,W_{iT})' \in \mathbb{R}^{T\times q},\; D_i \equiv (\delta_i',\, U_{i1},\ldots,U_{iT})'\in\mathbb R^d, \; d = q + T.$$

Assumption (ref)(i) ensures that the correlated random coefficients \(\beta_i\) admit a density function. This can be relaxed if interest is in the distribution function or the moments of $\beta_i$ instead.\footnote{Our identification strategy recovers the characteristic function of \(\beta_i\), which can be inverted via Fourier inversion to obtain either the density function or via L\'evy inversion to obtain the cumulative distribution function; see Lukacs1970.} Assumption (ref)(ii) allows arbitrary correlation between $\beta_i$ and $(X_i,W_i)$, while restricting the joint distribution of $(\beta_i,D_i,X_i,W_i)$. In particular, the assumption has the following implications. First, $D_i$ is statistically independent of $(X_i,W_i)$.\footnote{Heteroskedasticity \(\operatorname{Var}(D_i\mid W_i)\) can be allowed by weakening the statistical independence assumption. Identification results would then be stated conditional on $W_i$, which we do not pursue.} Second, $\beta_i$ and $D_i$ are conditionally independent given $(X_i,W_i)$. These implications allow us to, first, identify the characteristic function of $D_i$, and then that of $\beta_i$ via a deconvolution argument.

Assumption (ref)(ii) has observable implications when $q\neq0$. Let $x\in\mathbb{R}^{T\times p}$ and $w\in\mathbb{R}^{T\times d}$ denote matrix realizations of $X_i,\;W_i$, respectively. For finite second moments, for almost every \((x,w)\) in the support of \((X_i,W_i)\), it follows that \[ \operatorname{Var}(Y_i\mid X_i=x,W_i=w) = x\,\operatorname{Var}(\beta_i\mid x,w)\,x' + w\,\operatorname{Var}(D_i)\,w', \] where we used $\operatorname{Var}(D_i\mid X_i,W_i)=\operatorname{Var}(D_i),$ and $\operatorname{Cov}(\beta_i,D_i\mid X_i,W_i)=0.$ Therefore, the difference \[ \operatorname{Var}(Y_i\mid x,w) - w\,\operatorname{Var}(D_i)\,w'\in\mathbb R^{T\times T} \] must be positive semidefinite for almost every \((x,w)\). A negative eigenvalue of the associated matrix provides evidence against Assumption (ref)(ii).

Identification

The parameter of interest is $f_\beta$, the density of the correlated random coefficients. We sketch the logic behind our identification strategy below; the main result is in Section (ref).

Identification of $f_\beta$ proceeds in two steps. In the first step, we apply a transformation to (ref) that allows us to remove the $\beta_i$ component, which then yields the characteristic function of $D_i$. In the second step, we apply a different transformation to (ref) that isolates $\beta_i$. Given the identification of the characteristic function of $D_i$, this then allows us to recover $f_\beta$ by deconvolution. The choice of transformations depends on the relation between $T$ and the number of correlated random coefficients $p$. Concrete examples are given in Section (ref).

Let $\tau_1(X_i)\in\mathbb R^{T\times T}$ denote a measurable transformation. In the regular design $T>p$, it satisfies

so that, premultiplying (ref) by $\tau_1(X_i)$, yields:

We use (ref) to identify the characteristic function of $D_i$ in the regular design.

In the irregular design $T=p$, we take $\tau_1(X_i)=\operatorname{adj}(X_i)$. Using the identity

we premultiply (ref) by $\tau_1(X_i)$ to obtain:

The $\beta_i$ component vanishes on the stayer set $\{\det(X_i)=0\}$. We use (ref) to recover the characteristic function of $D_i$ in the irregular design, for observations with $\det(X_i)=0$.

Let $\tau_2(X_i)\in\mathbb R^{p \times T}$ be any measurable mapping satisfying

Premultiplying (ref) by $\tau_2(X_i)$ gives

Given the characteristic function of $D_i$ from the first step, we use (ref) to identify the density of $\beta_i$ via a deconvolution argument.\footnote{ The role of the transformations is not algebraic recovery of the random coefficients. Generally, neither (ref) nor (ref) identifies $D_i$ or $\beta_i$. For example, $\widetilde W_i=\tau_1(X_i)W_i$ contains at most $\operatorname{rank}(\tau_1(X_i))$ linearly independent equations for $d$ unknown components. Since $\operatorname{rank}(\widetilde W_i)\le \operatorname{rank}(\tau_1(X_i))<d$ in the designs of interest, the system is underdetermined for each $i$. Similarly, $\widetilde{\widetilde Y}_i=\beta_i+\widetilde{\widetilde W}_iD_i$ provides $p$ equations in the $p+d$ latent components $(\beta_i,D_i)$ and is likewise underdetermined at the unit level. }

Identification of the density function \texorpdfstring{$f_\beta$}{f_beta}

Throughout the analysis, we maintain the following assumptions, where $\tOne Y_i\in\mathbb R^T,\;\tOne W_i\in\mathbb R^{T\times d}$ are as defined in (ref) and $\tTwo Y_i\in\mathbb R^p,\;\tTwo W_i\in\mathbb R^{p\times d}$ are as defined in (ref). Note that $d\geq 1$ in the irregular design, while $d\geq 2$ in the regular design.\footnote{Consider the regular design with $d=1$. Then $T+q=1$, so it must be that $q=0$ and $T=1$, which implies that $p=0$.}

This is a standard assumption in deconvolution, ensuring that division in the Fourier domain is well defined. It rules out, e.g., uniform and truncated normal distributions.

Assumption (ref) is the support condition that makes the first-step characteristic-function argument informative in all relevant directions. Whenever the annihilated equation satisfies $\widetilde Y_i=\widetilde W_iD_i,$ we have, for any measurable projection vector \(\mu_i\), \[ \mu_i'\widetilde Y_i = \mu_i'\widetilde W_iD_i = \lambda_i'D_i . \] Since \(\lambda_i\) is measurable with respect to \((X_i,W_i)\) and \(D_i\) is independent of \((X_i,W_i)\), it follows that, for every \(v\in\mathbb R\), \[ \mathbb E\!\left[ \exp\{ \iota v\mu_i'\widetilde Y_i\} \,\middle|\, \lambda_i=\lambda \right] = \varphi_D(v\lambda), \] with the conditioning interpreted conditionally on the stayer event in the irregular design. Thus the first step identifies \(\varphi_D\) on the following set \[ \mathcal C_1 := \{v\lambda:\ v\in\mathbb R,\ \lambda\in\operatorname{supp}(\lambda_i)\}. \] Equivalently, it identifies \(\varphi_D\) along the rays generated by the normalized directions \(S_i\). Assumption (ref) requires these directions to be rich enough that \[ \overline{\operatorname{supp}(S_i)} = \mathbb S^{d-1} \] in the regular design, and \[ \overline{\operatorname{supp}(S_i\mid \det(X_i)=0)} = \mathbb S^{d-1} \] in the irregular design. Since characteristic functions are continuous, knowledge of \(\varphi_D\) on a dense collection of rays identifies \(\varphi_D\) on all of \(\mathbb R^d\). This is the sense in which Assumption (ref) ensures compatibility between the frequencies at which \(\varphi_D\) is recovered in the first step and the frequencies at which it is needed in the deconvolution step.

In the irregular design, \(\tau_1(X_i)=\operatorname{adj}(X_i)\). On the stayer set \(\{\det(X_i)=0\}\), and under the generic rank condition \[ \operatorname{rank}(X_i)=p-1, \] the adjugate matrix has rank one. Hence $\widetilde W_i$ has row space of dimension at most one. Each stayer therefore contributes at most one direction, up to scale. If a fixed row of \(\widetilde W_i\) is nonzero on the relevant support, one may use that row. For instance, if the first row is nonzero on the stayer set, one may take \(\mu_i=e_1\), giving \[ \lambda_i = \widetilde W_i'e_1 = W_i'\operatorname{adj}(X_i)'e_1 . \] A sufficient support condition is then \[ \overline{\operatorname{supp}\left( \frac{\lambda_i} {\|\lambda_i\|} \,\middle|\, \det(X_i)=0 \right)} = \mathbb S^{d-1}. \] If no fixed row is nonzero on the relevant support, one can instead define a measurable row-selection rule. Thus, in the irregular design, directional variation must come from the conditional distribution of \((X_i,W_i)\) on, or near, the singular set.

In the regular design, \(\tau_1(X_i)X_i=0\), and the annihilated equation $\widetilde Y_i=\widetilde W_iD_i$ holds for every observation. If \(\operatorname{rank}(X_i)=p\), the standard annihilator has rank \(T-p\), so \[ \operatorname{rank}(\widetilde W_i)\leq T-p. \] When the row space of \(\widetilde W_i\) has dimension greater than one, a single observation generates a family of possible directions, and the projection vector \(\mu_i\) selects one direction from that family. When the row space is one-dimensional, each observation contributes only one direction up to scale, and the required richness must come from cross-sectional variation in \((X_i,W_i)\).

A fixed projection vector need not exploit all available directional variation. For example, if \(\mu_i\equiv e_1\), then $ \lambda_i=\widetilde W_i'e_1$ uses only one fixed linear combination of the rows of \(\widetilde W_i\). The resulting directions may lie in a strict subset of \(\mathbb S^{d-1}\), even if the collection of row spaces of \(\widetilde W_i\) is itself rich. Allowing \(\mu_i=\mu(X_i,W_i)\) to vary with the data enlarges the attainable set of directions by selecting different linear combinations of the rows of \(\widetilde W_i\). Assumption (ref) holds when the normalized selected directions are dense in the unit sphere.

Assumption (ref) is an integrability requirement ensuring that the conditional density \(f_{\beta\mid X,W}(\cdot\mid x,w)\) exists and that the characteristic function of $\beta_i$ is integrable so that $f_\beta$ can be obtained via Fourier inversion.

Assumption (ref)(i) implies $\Pr(\det(X_i)=0)=0$, so stayers occur with probability zero. Identification in the irregular case therefore exploits the fact that $\det(X_i)$ can be arbitrarily close to zero with positive probability. The second part of (i) ensures that the joint conditional distribution used in the proof of Theorem (ref) (specifically, the distribution of $\mu_i'\widetilde{Y}_i$ given both $\widetilde{W}_i'\mu_i=\lambda$ and $\det(X_i)=d_x$) has a well-defined limit as $d_x\to 0$.\footnote{The conditioning notation $\mathbb{E}[\cdot\mid\widetilde{W}_i'\mu_i=\lambda,\,\det(X_i)=0]$ in Theorem (ref) is understood as the limit $\lim_{d_x\to 0} \mathbb{E}[\cdot\mid\widetilde{W}_i'\mu_i=\lambda,\,\det(X_i)=d_x]$.} This rules out pathological configurations in which the conditional relationship between $\widetilde{Y}_i$ and $\widetilde{W}_i'\mu_i$ changes discontinuously as the regressor matrix approaches singularity. A sufficient primitive condition is that the transformation \[ (X_i,W_i,\beta_i) \mapsto \bigl(\mu_i'\beta_i,\widetilde W_i'\mu_i,\det(X_i)\bigr) \] induce a joint density that is continuous in a neighborhood of \(\{\det(X_i)=0\}\), and that the marginal density of \((\widetilde W_i'\mu_i,\det(X_i))\) be continuous and bounded away from zero on the relevant conditioning set.

Assumption (ref)(ii) is the standard full-rank condition ensuring that $X_i$ has a well-defined $p$-dimensional column space for all observations in the regular design. This ensures that both $\tau_1(X_i)$ and $\tau_2(X_i)$ exist for all units.

Define \[ Z_i^* := \frac{\mu_i' \tOne Y_i}{\|\lambda_i\|}, \qquad \lambda_i := \tOne W_i'\mu_i \in \mathbb R^d, \qquad S_i := \frac{\lambda_i}{\|\lambda_i\|}\in\mathbb S^{d-1}, \] on the event $\{\lambda_i\neq 0\}$. On the event where the annihilated equation satisfies $\tOne Y_i=\tOne W_iD_i$, we have \[ Z_i^* = \frac{\mu_i'\tOne W_iD_i}{\|\lambda_i\|} = \frac{\lambda_i'D_i}{\|\lambda_i\|} = S_i'D_i . \] In the regular design this equality holds for every observation, whereas in the irregular design it holds on the stayer set $\{\det(X_i)=0\}$.

Since $S_i$ is a measurable function of $(X_i,W_i)$ and $D_i$ is independent of $(X_i,W_i)$, it follows that $D_i$ and $S_i$ are independent. Hence, conditional on $S_i=s$, the distribution of $Z_i^*$ coincides with the distribution of the scalar projection $s'D_i$. If $D_i$ admits a density $f_D$, then the density of this scalar projection is \[ \mathcal R f_D(s,u) := \int_{\{d\in\mathbb R^d:s'd=u\}} f_D(d)\,d\sigma(d), \] where $d\sigma$ denotes the induced Lebesgue measure on the hyperplane $\{d\in\mathbb R^d:s'd=u\}$. Thus, for observed directions $s$, \[ f_{Z^*\mid S}(u\mid s)=\mathcal R f_D(s,u), \] with the corresponding conditional-on-stayers version in the irregular design. Under the stronger directional support conditions stated in Corollary (ref), the Radon transform is known for enough directions to identify $f_D$; see, for example, HoderleinKMammen.

Examples

The examples in this section give concrete choices of the first- and second-step transformations in irregular and regular designs. The examples also show that, in low-dimensional settings, a fixed $\mu_i$ (rather than a data-dependent one) is often convenient.

Irregular design

The logic of the scalar irregular design extends to higher-dimensional irregular designs, as illustrated by the next example.

Regular design.

When $T>p$, the regressor matrix $X_i$ has full column rank with probability one under Assumption (ref)(ii). In this case the annihilator $\tau_1(X_i)$ and the left inverse $\tau_2(X_i)$ can be constructed for every observation.

A convenient choice for $\tau_1(X_i)$ is the orthogonal projector onto the null space of $X_i'$,

which satisfies $\tau_1(X_i)X_i=0$ identically. Thus the transformed equation $\widetilde{Y}_i=\widetilde{W}_iD_i$ holds for all units, so the first step can exploit the full sample. Any measurable $\tau_2(X_i)$ satisfying $\tau_2(X_i)X_i=I_p$ suffices for the second-step decomposition in (ref), since $X_i$ has full column rank and a left inverse exists for every observation.

In contrast to the irregular design, where first-step identification relies on a restricted subset of observations, in the regular design all observations contribute to the first step. Accordingly, the relevant support requirements are unconditional rather than conditional on a stayer-type event.

In higher-dimensional regular designs, the same logic applies, but the support condition involves the joint variation of $(X_i$, $W_i)$, and the measurable choice of $\mu_i$. When $q > 0$, variation in $X_i$ alone generally does not suffice. Example (ref) illustrates this for $(T,p,q)=(2,1,1)$.

Estimation of \texorpdfstring{$f_\beta$}{f_beta}

The parameter of interest is $f_\beta$, with characteristic function $\varphi_\beta(u)$, $u\in\mathbb R$. The results in Section (ref) imply that there exists an identified function $m_0:\mathbb R\to\mathbb{C}$ such that

Given (ref), define the weighted population criterion

where $d\nu(u)=\nu_0(u)\,du$ is a finite weighting measure with $\nu_0$ having full support on $\mathbb R$, so that $Q(\varphi)=0$ if and only if $\varphi(u)=\varphi_\beta(u)$ $\nu$-almost everywhere. \footnote{To see this from (ref): the ratio of the two conditional characteristic functions appearing in that expression identifies the conditional characteristic function $\varphi_{\beta|X,W}$, and $m_0(\cdot)$ is obtained by integrating this ratio with respect to the distribution of $(X,W)$.}

Estimation proceeds in two stages. The first stage, which is design-specific, constructs an estimator $\widehat{m}_N(u)$ of $m_0(u)$. The second stage, which is common across designs, recovers $f_\beta$ by minimizing the sample analog of (ref) over a finite-dimensional sieve. We describe the common second stage first, then specialize the first stage to the scalar irregular design $(T,p,q)=(1,1,0)$ and the regular design $(T,p,q)=(2,1,0)$. These are the configurations that arise in our application.

The identification results in Section (ref) apply to both regular and irregular designs in arbitrary dimensions. The estimators developed in this section specialize to the scalar irregular design and a low-dimensional regular design. The first-step estimation of $\varphi_D$ and second-step sieve minimum distance extend to higher-dimensional designs by replacing the scalar kernel regressions with their multivariate counterparts and the directional smoothing on $\mathbb{S}^1$ with smoothing on $\mathbb{S}^{d-1}$, but the additional tuning complexity and the curse of dimensionality in the first stage make a general-purpose implementation less practical. We therefore focus on the cases with empirical relevance.

Sieve minimum distance estimator

Let $\{q_s:s=0,\ldots,S-1\}\subset L^1(\mathbb R)$ be a vector of basis functions, $q^S(b)=(q_0(b),\ldots,q_{S-1}(b))'$, and for $\pi\in\mathbb R^S$ set $f_S(b;\pi)=q^S(b)'\pi$. The characteristic function of $f_S(\cdot\,;\pi)$ is then

Enforcing the unit-mass restriction $\varphi_S(0;\pi)=1$ defines the parameter space

The sieve minimum distance estimator minimizes the sample analog of (ref) over $\Pi_S$:

and the estimator of the density function is given by:

\paragraph{Hermite sieve.} For implementation we take $\{q_s\}$ to be the orthonormal Hermite functions on $L^2(\mathbb R)$: \[ q_s(v)=c_s^{-1/2} e^{-v^2/2}H_s(v), \qquad c_s=2^s s!\sqrt{\pi}, \] where $H_s$ denotes the physicists' Hermite polynomial of degree $s$. These functions are eigenfunctions of the Fourier transform $(\mathcal{F}q_s)(u)=\sqrt{2\pi}\,\imath^s q_s(u)$, so that

Since $f_\beta$ is real-valued and each $q_s$ is real, we restrict $\pi\in\Pi_S$. Under this restriction, the criterion $\widehat{Q}_N(\pi)$ is quadratic in $\pi$ and admits a closed-form constrained minimizer given by (ref) below.

\paragraph{Closed-form solution.} Approximating the integral in (ref) by a quadrature rule $\int g(u)\,d\nu(u)\approx\sum_{\ell=1}^L w_\ell g(u_\ell)$, define

Expanding the squared modulus in (ref) gives \[ \widehat{Q}_N(\pi) = \pi'\widehat{\Omega}\,\pi - 2\,\widehat{V}'\pi + \sum_{\ell=1}^L w_\ell\,\bigl|\widehat{m}_N(u_\ell)\bigr|^2, \] with unconstrained first-order condition $\widehat{\Omega}\pi=\widehat{V}$. Setting $A=z^S(0)\in\mathbb R^S$, and assuming $\widehat{\Omega}$ is invertible, \footnote{A sufficient condition is that the quadrature weights $w_\ell$ are positive and the $L\times S$ matrix with rows $z^S(u_\ell)'$ has full column rank $S$. By construction $\widehat{\Omega}$ is real symmetric and positive semidefinite; under the stated rank condition it is positive definite, and hence invertible.} the constrained minimizer subject to $A'\pi=1$ is

which satisfies $A'\widehat{\pi}=1$ by construction. The sieve estimator of $f_\beta$ is (ref) with $\widehat{\pi}$ given by (ref).

The remainder of this section constructs the design-specific first-stage estimator $\widehat{m}_N(u)$.

Scalar irregular design: \texorpdfstring{$(T,p,q)=(1,1,0)$}{(T,p,q)=(1,1,0)}

In the scalar irregular design, \[ Y_i=X_i\beta_i+D_i, \] where $D_i$ is scalar. Let $\tau_x>0$ be a mover--stayer threshold and define \[ \mathcal{S}_N(\tau_x)=\{i:\ |X_i|<\tau_x\}, \qquad \mathcal{M}_N(\tau_x)=\{i:\ |X_i|\ge \tau_x\}. \] Observations in $\mathcal{S}_N(\tau_x)$ are stayers, used to estimate the disturbance characteristic function $\varphi_D$, while observations in $\mathcal{M}_N(\tau_x)$ are movers, used in the deconvolution step.

For movers, the second-step transformation is \[ \tTwo Y_i=\frac{Y_i}{X_i} = \beta_i+\frac{D_i}{X_i}, \qquad i\in\mathcal{M}_N(\tau_x). \]

We estimate $\varphi_{\tTwo Y\mid X}(u\mid X_i)=\mathbb{E}[e^{\imath u\tTwo Y_i}\mid X_i]$ by kernel regression on the mover sample:

where \[ \omega_{k,i}^{(x)} = \frac{K_{h_x}(X_k-X_i)} {\sum_{s\in\mathcal{M}_N(\tau_x)}K_{h_x}(X_s-X_i)}. \] We estimate $\varphi_D$ by local smoothing at $X_i=0$, where $Y_i=D_i$: \footnote{The sample average in (ref) is taken over movers only, since $\tTwo Y_i=Y_i/X_i$ is defined only for $i\in\mathcal{M}_N(\tau_x)$; the exclusion of stayers is asymptotically negligible as $\tau_x\to 0$.}

For each mover $i$, the trimmed conditional ratio is

where $\tau_{\rm den}>0$ is a denominator trimming threshold. Averaging over movers gives the first-stage estimator

The tuning parameters $h_x$, $h_0$, $\tau_x$, and $\tau_{\rm den}$ are discussed in Section (ref).

Regular design: \texorpdfstring{$(T,p,q)=(2,1,0)$}{(T,p,q)=(2,1,0)}

In the regular design, \[ Y_i=X_i\beta_i+D_i, \qquad X_i=

, \qquad D_i=

. \] We use the same transformations as in Example (ref). Define \[ \tau_1(X_i) =

, \qquad \tau_2(X_i)=\frac{X_i'}{X_i'X_i},\qquad X_i'X_i>0. \]

Setting $\tOne Y_i=\tau_1(X_i)Y_i$ and $Y_i^*:=e_1'\tOne Y_i$, the condition $\tau_1(X_i)X_i=0$ gives \[ Y_i^*=X_{i2}Y_{i1}-X_{i1}Y_{i2}=\lambda_i'D_i, \qquad \lambda_i=

. \] Each observation therefore identifies $\varphi_D$ along the line generated by $\lambda_i$: for any $v\in\mathbb R$ and any $\lambda$ in the support of $\lambda_i$, \[ \varphi_D(v\lambda) = \mathbb{E}\!\left[ e^{\imath vY_i^*} \,\middle|\, \lambda_i=\lambda \right]. \]

The second-step transformation isolates $\beta_i$: \[ \tTwo Y_i = \tau_2(X_i)Y_i = \frac{X_i'Y_i}{X_i'X_i} = \beta_i+\frac{X_i'D_i}{X_i'X_i}. \]

We estimate $\varphi_{\tTwo Y\mid X}(u\mid X_i)=\mathbb{E}[e^{\imath u\tTwo Y_i}\mid X_i]$ by kernel regression:

where \[ \omega_{k,i}^{(X)} = \frac{K_{h_X}(X_k-X_i)} {\sum_{s=1}^N K_{h_X}(X_s-X_i)}. \]

To estimate $\varphi_D$, normalize the first-step direction: \[ S_i=\frac{\lambda_i}{\|\lambda_i\|}\in\mathbb{S}^1, \qquad \|\lambda_i\|=\sqrt{X_{i1}^2+X_{i2}^2}. \] Since $Y_i^*/\|\lambda_i\|=S_i'D_i$ depends only on the unit direction $S_i$, for any $\xi\in\mathbb R^2\setminus\{0\}$ with polar decomposition $\xi=\|\xi\|\,s(\xi)$, $s(\xi)=\xi/\|\xi\|$, \[ \mathbb{E}\!\left[ e^{\imath \|\xi\|\,Y_i^*/\|\lambda_i\|} \,\middle|\, S_i=s(\xi) \right] = \varphi_D(\xi). \] We therefore estimate $\varphi_D$ by smoothing over directions on $\mathbb{S}^1$: \[ \widehat\varphi_D(\xi) =

\] where $\omega_j^{(S)}(s)=K(\|S_j-s\|/h_S)/\sum_{k=1}^N K(\|S_k-s\|/h_S)$, and $\widehat\varphi_D(0)=1$ by the property $\varphi_D(0)=1$.

For each observation $i$, we form the trimmed conditional ratio:

Since $X_i'X_i>0$ for all $i$, the average is taken over the full sample: \footnote{In the irregular design the average is restricted to movers $\mathcal{M}_N(\tau_x)$ because $\tTwo Y_i=Y_i/X_i$ requires $|X_i|\ge\tau_x$. Here $\tTwo Y_i=X_i'Y_i/(X_i'X_i)$ is well-defined for all $i$ under the maintained assumption $X_i'X_i>0$, so no such restriction is needed.}

Implementation

Algorithms (ref) and (ref) in Section (ref) summarize the first-stage implementation for the scalar irregular and regular designs, respectively. In both cases the output is a vector of estimated characteristic function values $\{\widehat{m}_N(u_\ell)\}_{\ell=1}^L$ on a prescribed frequency grid $\{u_\ell\}_{\ell=1}^L \subset [-U_N,U_N]$. Algorithm (ref) in Section (ref) describes the common second stage, which takes these values as input and returns $\widehat{f}_\beta$.

The estimator depends on a number of tuning parameters. Three are common to all algorithms: the denominator trimming threshold $\tau_{\mathrm{den}}$, the frequency truncation $U_N$, and the sieve dimension $S$. The remaining parameters are design-specific. In the scalar irregular design, the first stage uses the bandwidth $h_0$ for estimating $\varphi_D$ from stayers, the bandwidth $h_x$ for estimation $\varphi_{\tTwo Y | X}$ from movers, and the mover--stayer threshold $\tau_x$. The role of $\tau_x$ is distinct from that of the bandwidths: it defines the stayer sample and thereby localizes the denominator estimator in regressor space. Crucially, $\widehat\varphi_D$ is constructed once from the stayer subsample and subsequently evaluated at the mover-specific arguments $u/X_i$. In the regular design, the first stage uses the bandwidth $h_X$ for kernel regression of the second-step transformed outcome on $X_i$, and the bandwidth $h_S$ for smoothing over normalized directions $S_i\in\mathbb{S}^1$; no mover--stayer partition is needed.

The integrals entering both the first-stage frequency evaluations and the second-stage criterion are approximated on the grid $\{u_\ell\}_{\ell=1}^L$ with positive weights $\{w_\ell\}_{\ell=1}^L$, which determine the quadrature objects $\widehat{\Omega}$ and $\widehat{V}$ in (ref). The choice of all tuning parameters is discussed next.

Heuristic rate discussion and tuning parameter selection in the scalar irregular case

This subsection provides heuristic guidance for selecting $(\tau_x, h_0, U_N, \tau_{\mathrm{den}})$. Conditional on these choices, the remaining tuning parameters $(h_x, S)$ are selected by cross-validation. The tuning parameters play different roles in the estimator. The first group of tuning parameters primarily determines how the stayer conditioning event is approximated and how stable the deconvolution step is, while the second group primarily determines the degree of smoothing and the complexity of the sieve approximation.

Let the numerator and denominator estimation errors be defined as \[ \rho_N := \sup_{\substack{|u|\le U_N \\ |x|\ge \tau_x}} \left| \widehat{\varphi}_{\widetilde{\widetilde{Y}}\mid X}(u\mid x) - \varphi_{\widetilde{\widetilde{Y}}\mid X}(u\mid x) \right|, \qquad b_N := \sup_{|v|\le U_N/\tau_x} \left| \widehat{\varphi}_D(v) - \varphi_D(v) \right|, \] and define the stability factor \[ \Delta_N := \inf_{\substack{|u|\le U_N \\ |x|\ge \tau_x}} |\varphi_D(u/x)|. \]

To separate truncation from estimation error, let \(f_{\beta,U_N}\) denote the population density obtained by restricting Fourier inversion to the frequency region \(|u|\le U_N\), and let \(f_{\beta,S}\) denote its sieve approximation. A heuristic decomposition then yields

The final term captures the bias induced by restricting inversion to a bounded frequency region and is decreasing in \(U_N\). The middle terms reflect the ratio structure of the estimator: numerator errors enter linearly through multiplication by \(1/\varphi_D\), while denominator errors enter through the derivative of \(1/\varphi_D\), leading to amplification factors of order \(\Delta_N^{-1}\) and \(\Delta_N^{-2}\), respectively.

The decomposition in (ref) highlights how the tuning parameters enter the different components of the estimation error. The threshold $\tau_x$ and bandwidth $h_0$ determine the accuracy of the denominator estimator through $b_N$, while the trimming parameter $\tau_{\mathrm{den}}$ and frequency cutoff $U_N$ jointly determine the stability factor $\Delta_N$. The cutoff $U_N$ also governs the truncation bias induced by restricting inversion to a bounded frequency region. The numerator bandwidth $h_x$ controls the estimation error $\rho_N$, and the sieve dimension $S$ determines the sieve approximation error. The discussion below explains how these parameters are chosen.

In the scalar irregular design, $\varphi_D(v)=\mathbb{E}[e^{\imath v Y_i}\mid X_i=0]$. If $P(X_i=0)>0$, this quantity could be estimated by an average over $\{i:X_i=0\}$, so that no bandwidth would be required. In the continuous case, however, $P(X_i=0)=0$, and conditioning on $X_i=0$ must be approximated using observations with $|X_i|$ close to zero. We therefore estimate $\varphi_D(v)$ using the stayer sample $\{|X_i|<\tau_x\}$ together with kernel weights $h_0$. To ensure that the resulting estimator mimics the infeasible conditional average at $X_i=0$, we impose

so that $X_i/h_0 \to 0$ uniformly over $|X_i|<\tau_x$. Under this regime, the kernel weights are approximately constant over the stayer window, and the estimator behaves like an averaging estimator over $\{|X_i|<\tau_x\}$. Consequently, the relevant localization scale is $\tau_x$, while $h_0$ serves only to ensure regularity of the estimator.

The numerator bandwidth $h_x$ is selected by cross-validation. To ensure that numerator smoothing is not affected by the thresholding operation, our heuristics presumes that $\tau_x = o(h_x)$, in which case, smoothing in the numerator step is governed by $h_x$.

The parameters $U_N$ and $\tau_{\mathrm{den}}$ must be chosen jointly, as both govern the stability of the deconvolution ratio. Since $|u|\le U_N$ and $|x|\ge \tau_x$ imply $|u/x|\le U_N/\tau_x$, the denominator is evaluated over the frequency region \[ |v| \le \frac{U_N}{\tau_x}, \] so that \[ \Delta_N \asymp \inf_{|v|\le U_N/\tau_x} |\varphi_D(v)|. \] Increasing $U_N$ expands this region and reduces truncation bias, but also decreases $\Delta_N$ as $\varphi_D(v)$ decays in $|v|$, thereby amplifying denominator noise. Stability therefore requires that $\Delta_N$ remains at least of the same order as $\tau_{\mathrm{den}}$. This implies that $U_N$ cannot grow faster than allowed by the decay of $\varphi_D$. In particular, if $|\varphi_D(v)| \asymp (1+|v|)^{-\alpha}$ (ordinary smooth), then \[ \Delta_N \asymp (U_N/\tau_x)^{-\alpha}, \quad\text{so that}\quad U_N \lesssim \tau_x\,\tau_{\mathrm{den}}^{-1/\alpha}, \] whereas if $|\varphi_D(v)| \asymp \exp(-c|v|^\gamma)$ (supersmooth), \[ \Delta_N \asymp \exp\!\left(-c(U_N/\tau_x)^\gamma\right), \quad\text{so that}\quad U_N \lesssim \tau_x \bigl[\log(1/\tau_{\mathrm{den}})\bigr]^{1/\gamma}. \] Thus $U_N$ and $\tau_{\mathrm{den}}$ are linked through the requirement that the smallest denominator value over the working frequency region remains of the same order as the trimming threshold.

These considerations lead to the following rules of thumb. Let $s_0$ index the smoothness of $x \mapsto \mathbb{E}[e^{\imath vY_i}\mid X_i=x]$ at $x=0$. Under $$\tau_x \asymp N^{-\kappa},$$ the averaging logic with $\tau_x=o(h_0)$ suggests the choice $$h_0 \asymp (N\tau_x)^{-1/(2s_0+1)}\asymp N^{-(1-\kappa)/(2s_0+1)},\qquad \kappa>\frac{1}{2(s_0+1)}.$$

The frequency trimming tuning parameter $U_N$ is pinned down by the stability constraint and the tail behavior of $\varphi_D$. Under \[ \tau_{\mathrm{den}} \asymp N^{-\eta},\quad \eta >0, \] in the ordinary smooth case with, e.g., $|\varphi_D(v)| \asymp (1+|v|)^{-\alpha}$ for some $\alpha>0$, the stability condition yields polynomial growth for $U_N$: \[ U_N \asymp N^{-\kappa + \eta/\alpha},\qquad \frac{\eta}{\alpha} > \kappa. \]

In the supersmooth case with, e.g., $|\varphi_D(v)| \asymp \exp(-c|v|^\gamma)$ for some $c,\gamma>0$, to obtain a diverging $U_N$, the trimming threshold must decay faster than polynomially. For example, if \[ \tau_{\mathrm{den}} \asymp \exp(-N^{r}), \quad r>0, \] then \[ U_N \asymp \tau_x N^{r/\gamma} \asymp N^{-\kappa + r/\gamma}, \qquad \frac{r}{\gamma} > \kappa. \]

Under these choices, each term in (ref) is controlled: the denominator error $b_N$ is governed by $\tau_x$, the stability factor is maintained at the order of $\tau_{\mathrm{den}}$, and the truncation bias is reduced by allowing $U_N$ to grow at the maximal rate compatible with stability.

Cross-Validation

A practical way to select the remaining tuning parameters is via $K$-fold cross-validation. Conditional on $(\tau_x, h_0, U_N, \tau_{\mathrm{den}})$ choices, cross-validation is used to select the numerator bandwidth $h_x$ and the sieve dimension $S$.

We develop the cross-validation procedure below for the irregular design. In the regular design, the directional bandwidth $h_S$, which governs the estimation of $\varphi_D$ by smoothing over $S_i\in\mathbb{S}^{d-1}$, plays the role of $h_0$ and is fixed; the numerator bandwidth $h_X$ plays the role of $h_x$ and is selected jointly with $S$ by cross-validation.

In the irregular design we fix $\tau_x$ using a rule motivated by

where $\kappa$ governs the size of the stayer window. In the simulations and application, we use $\kappa = 1/3$.

In the regime $\tau_x = o(h_0)$, the bandwidth $h_0$ is set

Setting $s_0=1$, yields $h_0 \asymp N^{-2/9}$ when $\kappa=1/3$.

To construct the candidate grid for $h_x$, we compute a reference bandwidth

where $X_i^{\mathcal M}$ denotes the mover subsample and $N_m = |\mathcal M_N(\tau_x)|$. The candidate grid is $\Theta_{h_x} = \{c \cdot h_x^{\mathrm{ref}} : c \in \mathcal C\}$ for a finite grid of positive multipliers $\mathcal C$.

Finally, we set $\tau_{\mathrm{den}}$ at a small value (e.g., $10^{-4}$), which determines the minimal magnitude of the denominator that is deemed numerically stable. Conditional on this choice, the cutoff $U_N$ is selected using the stability considerations in Section (ref). In particular, $U_N$ is increased conservatively until the empirical denominator $|\widehat{\varphi}_D(v)|$ becomes small relative to $\tau_{\mathrm{den}}$, so that the working frequency region remains restricted to values for which the ratio is well behaved.

For each sample size $N$, fix a frequency grid $\{u_j\}_{j=1}^J \subset [-U_N,U_N]$ with associated nonnegative weights $\{\varpi_j\}_{j=1}^J$, where the cutoff $U_N$ is determined as above and is held fixed across candidate values of $\theta=(h_x,S)$ in the cross-validation search.

Stabilized cross-validation

To describe the cross-validation procedure, let $\mathcal I_1, \dots, \mathcal I_K$ be a partition of the sample indices $\{1,\dots,N\}$, and let $\mathcal I_{-k} = \{1,\dots,N\}\setminus \mathcal I_k$ denote the training sample associated with fold $k$.

For each fold $k$ and candidate $\theta$, define the training and validation stayer and mover sets \[ \mathcal S_{-k}(\tau_x) = \{i\in \mathcal I_{-k}: |X_i|<\tau_x\}, \qquad \mathcal M_{-k}(\tau_x) = \{i\in \mathcal I_{-k}: |X_i|\ge \tau_x\}, \] \[ \mathcal S_k(\tau_x) = \{i\in \mathcal I_k: |X_i|<\tau_x\}, \qquad \mathcal M_k(\tau_x) = \{i\in \mathcal I_k: |X_i|\ge \tau_x\}, \] with cardinalities $N_{m,-k}=|\mathcal M_{-k}(\tau_x)|$ and $N_{m,k}=|\mathcal M_k(\tau_x)|$. Since $\tau_x$ is fixed, this partition does not vary across candidates $\theta$.

On the training sample, estimate

and on the validation sample, estimate

Since $h_0$ and $\tau_x$ are held fixed, $\widehat\varphi_D^{(-k)}$ and $\widehat\varphi_D^{(k)}$ are computed once per fold and do not vary across candidate values of $\theta=(h_x,S)$.

For each candidate bandwidth $h_x$ and each mover $i\in \mathcal M_{-k}(\tau_x)$, estimate

Form the trimmed training ratio

and average over the training movers:

Apply the second-stage sieve minimum distance estimator with sieve dimension $S$ to $\{\widehat m_N^{(-k)}(u_j)\}_{j=1}^J$ to obtain the training-fold sieve coefficient vector $\widehat\pi^{(-k)}$, and define the implied training-fold characteristic function estimator

Standard cross-validation performs poorly in deconvolution problems because the validation target inherits the instability of the estimator. When the candidate bandwidth $h_x$ is small, the validation ratio $\widehat{R}_N^{(k)}(u \mid X_i)$ is dominated by high-frequency noise amplified through division by $\widehat{\varphi}_D^{(k)}(u/X_i)$. If the training estimator uses the same small $h_x$, both sides of the cross-validation comparison are noisy at high frequencies, and the criterion effectively compares noise to noise.\footnote{This issue is well documented in the bandwidth selection literature; see ScottTerrell1987 for the distinction between biased and unbiased cross-validation, and HallMarronPark1992 for a discussion of how presmoothing affects the bias--variance tradeoff. Kent2024 develops a related stabilized cross-validation approach for smoothness-penalized deconvolution, showing that standard CV leads to severe undersmoothing in ill-posed problems.}

We adopt a stabilized cross-validation approach in which the validation target is constructed using a fixed pilot bandwidth $g_{\mathrm{pilot}}$ that does not vary with the candidate $h_x$. This produces an oversmoothed, low-variance validation target, so that the CV criterion compares candidate estimators against a stable reference rather than against an equally noisy holdout estimate. The pilot bandwidth is set to

where $c_{\mathrm{pilot}} \in [1.5, 3]$ and $r > 5$. Since $r > 5$ implies $N^{-1/r} > N^{-1/5}$, the reference bandwidth $h^{\mathrm{ref}}$ decays more slowly than the standard Silverman rate, so that $g_{\mathrm{pilot}}$ is larger than a typical bandwidth choice. This produces an oversmoothed validation target, which provides a stable reference against which candidate estimators are compared.

\paragraph{Heuristic justification.} Let $\theta=(h_x,S)$ denote the smoothing and sieve parameters. The infeasible finite-grid risk is \[ R_k(\theta) = \sum_{j=1}^J \varpi_j \left| \widehat\varphi_{\beta,S}^{(-k)}(u_j;\theta) - \varphi_\beta(u_j) \right|^2 . \] Cross-validation replaces the unknown target $\varphi_\beta$ by a validation estimate constructed on the hold-out fold. In the present deconvolution problem, the natural validation object is itself a ratio estimator, \[ \widehat m_N^{(k)}(u;h) = \frac{1}{N_{m,k}} \sum_{i\in\mathcal M_k} \frac{ \widehat\varphi_{\tilde Y\mid X}^{(k)}(u\mid X_i;h) }{ \widehat\varphi_D^{(k)}(u/X_i) } \mathbf{1}\bigl\{ |\widehat\varphi_D^{(k)}(u/X_i)|>\tau_{\rm den} \bigr\}. \] For a generic bandwidth $h$, write heuristically \[ \widehat m_N^{(k)}(u;h) \approx \varphi_\beta(u)+b_h(u)+\xi_h^{(k)}(u), \] where $b_h(u)$ is smoothing bias, $\xi_h^{(k)}(u)$ is stochastic error, and the approximation reflects a first-order expansion of the ratio estimator. The difficulty is that $\xi_h^{(k)}(u)$ is amplified by division by $\widehat\varphi_D^{(k)}(u/X_i)$, especially at frequencies for which the denominator is small. Thus a validation target constructed with a small bandwidth can be dominated by high-frequency noise.

We therefore construct the validation target using a fixed oversmoothed pilot bandwidth, \[ \widehat m_N^{(k)}(u; g_{\mathrm{pilot}}) \approx \varphi_\beta(u)+b_g(u)+\xi_g^{(k)}(u). \] The purpose of $g_{\mathrm{pilot}}$ is to reduce the numerator noise before that noise is amplified by the inverse characteristic function. Since the denominator $\widehat{\varphi}_D^{(k)}$ does not depend on $g_{\mathrm{pilot}}$, a larger pilot bandwidth reduces the variance of the numerator and hence the variance of the ratio.

Substitution into the cross-validation criterion gives

The second term on the right hand side is common across candidate values of $\theta$ and therefore does not affect the ranking of tuning parameters. The last term is the relevant distortion. Sample splitting mitigates its stochastic component in the sense that there is no systematic bias, since $\widehat\varphi_{\beta,S}^{(-k)}$ and $\xi_g^{(k)}$ are computed on disjoint folds. The remaining pilot-bias component is the price paid for stabilization: the pilot bandwidth is chosen to trade a small, smooth bias against a substantial reduction in the variance of the validation target.

Thus the stabilized criterion is intended to approximate the infeasible risk up to a candidate-invariant term and a reduced cross term. This does not establish optimality of the selected tuning parameters, but it explains why stabilizing the validation target is preferable to comparing the estimator with another noisy inverse estimate. A related approach with formal consistency guarantees is developed by Kent2024; that method estimates risk for a smaller hypothetical sample size and requires knowledge of asymptotic rates to rescale the selected parameter.

For each validation mover $i\in \mathcal M_k(\tau_x)$, estimate the numerator using the pilot bandwidth:

Form the validation ratio

and average over the validation movers:

This construction uses the pilot bandwidth $g_{\mathrm{pilot}}$ for the numerator smoothing while the denominator $\widehat\varphi_D^{(k)}$ continues to be estimated using the fixed bandwidth $h_0$. Since the validation target does not depend on the candidate $h_x$, it can be computed once per fold outside the loop over candidate bandwidths, which provides a computational advantage.

Repeated cross-validation and selection

To reduce sensitivity to the random fold assignment, we employ repeated cross-validation. Let $n_{\mathrm{rep}}$ denote the number of repetitions. For each repetition $r = 1, \ldots, n_{\mathrm{rep}}$, we draw an independent random partition of the sample into $K$ folds and compute the fold-specific holdout discrepancy

where $\theta = (h_x, S)$. The repetition-specific cross-validation score is $\mathrm{CV}^{(r)}(\theta) = K^{-1}\sum_{k=1}^{K}\mathrm{CV}_k^{(r)}(\theta)$.

Aggregating across repetitions, define

and the standard error $\mathrm{SE}(\theta) = \mathrm{SD}(\theta) / \sqrt{n_{\mathrm{rep}}}$.

One-standard-error rule

In ill-posed inverse problems, favoring simpler models can improve stability. We adopt a one-standard-error rule: let $\theta^* = \arg\min_\theta \overline{\mathrm{CV}}(\theta)$ denote the minimizer of the mean CV score, and let $\mathrm{SE}^* = \mathrm{SE}(\theta^*)$ be the standard error at the minimum. Define the set of near-optimal candidates

Among candidates in $\Theta_{\mathrm{1SE}}$, we select the pair with the largest $h_x$ (most smoothing) and, among those, the smallest $S$ (simplest sieve). This lexicographic ordering favors stability: larger bandwidths reduce variance in the numerator estimation, and smaller sieve dimensions reduce the complexity of the second-stage projection.

The selected tuning parameters are

Algorithm (ref) in Section (ref) summarizes the procedure.

Bootstrap Inference

We assess sampling uncertainty using a nonparametric pairs bootstrap.\footnote{The method here is provisional.} For each bootstrap replication $b = 1, \ldots, B$:

Inference is therefore for the post-processed density estimator. In the application, we use $B = 499$ replications.

Confidence intervals

We construct pointwise confidence intervals using the basic (reverse-percentile) bootstrap method. For each evaluation point $b$ in the density grid:

This is equivalent to the formula $\mathrm{CI}_{1-\alpha}(b) = [2\widehat{f}_\beta(b) - \widehat{q}_{1-\alpha/2}^*(b),\; 2\widehat{f}_\beta(b) - \widehat{q}_{\alpha/2}^*(b)]$, where $\widehat{q}_\tau^*(b)$ denotes the $\tau$-quantile of the bootstrap estimates $\{\widehat{f}_\beta^{*1}(b), \ldots, \widehat{f}_\beta^{*B}(b)\}$. After construction, confidence bands are truncated below at zero to respect the non-negativity constraint on densities.

Standard errors and moments

Bootstrap standard errors are computed as \[ \widehat{\mathrm{SE}}(b) = \mathrm{SD}\bigl(\widehat{f}_\beta^{*1}(b), \ldots, \widehat{f}_\beta^{*B}(b)\bigr). \] Since the standard deviation is translation-invariant, this quantity does not depend on whether the bootstrap distribution is centered.

Moments of the random coefficient distribution are computed by numerical integration over the post-processed density: \[ \widehat{\mathbb{E}}[\beta_i] = \int b\, \widehat{f}_\beta(b)\, db, \qquad \widehat{\mathrm{Var}}(\beta_i) = \int b^2\, \widehat{f}_\beta(b)\, db - \bigl(\widehat{\mathbb{E}}[\beta_i]\bigr)^2. \] Bootstrap standard errors and confidence intervals for moments are constructed analogously, applying the same integration to each bootstrap density $\widehat{f}_\beta^{*b}$.

Monte Carlo: Finite sample performance

Irregular design $(T,p,q)=(1,1,0)$

We consider the following data generating process. For each $i=1,\dots,N$,

where the covariates, random coefficients, and error terms are generated as follows:

with $\varepsilon_{\beta,i}\sim F_{\varepsilon_\beta}$, $\varepsilon_{1i}\sim F_{\varepsilon_{1}}$, and $\varepsilon_{2i}\sim F_{\varepsilon_{2}}$, where $\mathrm{Var}(\varepsilon_{1i})=\sigma^2_{\varepsilon_1}=1$ and $\mathrm{Var}(\varepsilon_{2i})=\sigma^2_{\varepsilon_2}=2$. The innovations $\varepsilon_{1i}$ and $\varepsilon_{2i}$ are independent and mean-zero; their distributions $F_{\varepsilon_{1}}$, $F_{\varepsilon_{2}}$, and $F_{\varepsilon_\beta}$ vary across specifications as described below. After first-differencing, the model reduces to \[ Y_i = X_i\beta_i + D_i, \] where $X_i = z_i\sim\mathcal N(0,4)$ and $D_i = \varepsilon_{2i}+(\theta-1)\varepsilon_{1i}$. The random coefficient $\beta_i$ is a scale mixture: since $\beta_i = \zeta_i\varepsilon_{\beta,i}$ with $\zeta_i = 1+\delta X_i$, $\beta_i$ and $X_i$ are statistically dependent. The error terms $u_{it}$ are serially correlated through the parameter $\theta$.

We consider four specifications, varying the distribution of the innovation $\varepsilon_{\beta,i}$, hence the distribution of the CRC, and the distribution of the error terms:

In specifications (a)--(c), the disturbance $D_i = \varepsilon_{2i}+(\theta-1)\varepsilon_{1i}$ is a linear combination of Gaussians, so its characteristic function $\varphi_D$ decays exponentially (supersmooth). In specification (d), the Laplace innovations yield $\varphi_D(t) = O(|t|^{-4})$ (ordinary smooth), while the distribution of $\beta_i$ is unchanged. Comparing (c) and (d) isolates the effect of the smoothness class of $D_i$ on the estimator's ability to recover multimodal features.

In each case, the unconditional density of $\beta_i$ is the scale mixture \[ f_\beta(b) = \int |s|^{-1}\,f_{\varepsilon_\beta}(b/s)\,f_\zeta(s)\,ds, \] where $f_\zeta$ is the density of $\zeta_i\sim\mathcal N(1,4\delta^2)$. When $\delta$ is small, the scaling factor $\zeta_i$ remains close to one with high probability, and the shape of $f_{\varepsilon_\beta}$ is largely preserved in the unconditional density $f_\beta$.

\paragraph{Implementation.} We generate $N=2000$ observations and run $100$ Monte Carlo replications. The stayer threshold is set to $\tau_x = c_\tau\cdot\tau_x^{\mathrm{ref}}$ with $\tau_x^{\mathrm{ref}}$ defined in (ref). In specifications (a) and (b), $c_\tau=4$, yielding approximately $25\%$ of observations classified as stayers across replications; in the bimodal specifications (c) and (d), $c_\tau=5$, yielding approximately $31\%$ stayers. The bandwidth $h_0$ for the stayer denominator estimator is set proportional to the threshold, $h_0 = c_0\cdot\tau_x$ with $c_0 = 1$.\footnote{This deviates from the discussion in Section (ref). The discussion there motivates the roles of \(\tau_x\) and \(h_0\), while the implementation localizes the denominator estimator on the same near-stayer region used to define the stayer sample.} The denominator trimming threshold is $\tau_{\mathrm{den}}=10^{-4}$, and the frequency grid consists of $L=101$ equally spaced points on $[-U_N,U_N]$ with $U_N=4$.

The tuning parameters $(h_x,S)$ are selected by $5$-fold cross-validation as described in Section (ref) and summarized in Algorithm (ref). In specifications (a) and (b), the bandwidth $h_x$ is global (fixed across movers), and the candidate grid consists of the multiples $$\Theta_{h_x}=\{0.5,\,0.75,\,1,\,1.5,\,2\}\times h_x^{\mathrm{ref}},$$ where $h_x^{\mathrm{ref}}$ is Silverman's rule on the mover subsample. The candidate grid for $S$ is $\{3,5,\dots,15\}$ in specification (a) and $\{3,5,\dots,19\}$ in specification (b). The weighting function $\nu$ is the standard normal density.

In the bimodal specifications (c) and (d), the global bandwidth is replaced by a $k$-nearest-neighbor ($k$-NN) adaptive rule (see, e.g., GaoOhViswanath2017). For each mover $i$, the local bandwidth $h_x(X_i)$ is set equal to the distance from $X_i$ to its $k$-th nearest neighbor in the mover subsample. Cross-validation selects $(k,S)$ jointly over $k\in\{5,10,15,20,30\}$ and $S\in\{7,9,\dots,19\}$. The weighting function is the Student-$t$ density with $3$ degrees of freedom, which places more mass on the higher frequencies needed to resolve bimodality.

In all four specifications, the feasibility thresholds are $\Gamma_{\max}=100$ and $\rho_{\max}=0.5$ (see Remark (ref)). After estimation, the sieve density is truncated below at zero and renormalized to integrate to one.

\paragraph{Results.} Figure (ref) displays the estimation results for all four specifications. In each panel, the black solid line is the true density $f_\beta$, the dashed line is the pointwise average of $\hat f_\beta$ across replications, the solid gray line is the pointwise median, and the shaded region is the pointwise interquartile range.

In the symmetric case (panel (a)), the estimator performs well. The pointwise median tracks the true density closely, and the interquartile range is narrow throughout the support. The cross-validation procedure selects $S=3$ in the large majority of replications ($79$ out of $100$), with occasional selections of $S=5$ ($13$ replications) and $S=7$ ($7$ replications). The selected bandwidth $h_x$ has a median of approximately $0.71$, with an interquartile range of $[0.36,\,0.95]$.

In the skewed case (panel (b)), the estimator captures the asymmetry and the location of the mode accurately. The pointwise median follows the true density closely, while the pointwise average shows a slight negative bias near the peak and some positive mass in the far left tail where the true density is near zero. The interquartile range is wider than in the symmetric case, reflecting the greater difficulty of the estimation problem. The cross-validation selects a median sieve dimension of $S=9$ (interquartile range $[7,11]$) across replications, with the full distribution spread over $S\in\{5,\dots,19\}$. The selected bandwidth $h_x$ has a median of $0.93$ with an interquartile range of $[0.48,\,0.96]$. This variability reflects the sensitivity of the numerator characteristic-function estimator to the local density of movers, as discussed in Section (ref).

The bimodal case is considerably more demanding. Panel (c) of Figure (ref) displays the results under specification (c) (Gaussian errors), where the $k$-NN adaptive bandwidth and the sieve dimension are selected jointly by cross-validation. Across $100$ replications, the procedure selects $k=30$ in $88$ replications and $S=7$ in $93$ replications. While the estimator identifies the general location and spread of the density, it does not fully resolve the trough between the two modes: the pointwise median and average estimates are flatter than the true density in the bimodal region. The $k$-NN adaptive bandwidth and the Student-$t$ weighting help, but the attenuation of the bimodal feature persists.

The difficulty can be traced to the supersmooth character of the Gaussian error distribution. Bimodality in $f_\beta$ manifests as oscillations in the characteristic function $\varphi_\beta$ at moderate to high frequencies. Recovering these oscillations requires accurate estimation of the ratio $\hat\varphi_{\tTwo Y\mid X}(u\mid X_i)/\hat\varphi_D(u/X_i)$ at frequencies where the denominator $\varphi_D$ is small. Under Gaussian errors, $\varphi_D$ decays exponentially: for the baseline DGP, $|\varphi_D(v)| = \exp(-\mathrm{Var}(D_i)\,v^2/2)$, which at the argument $v = u/X_i = 3$ equals approximately $8\times 10^{-5}$---close to the trimming threshold $\tau_{\mathrm{den}} = 10^{-4}$. The deconvolution ratio is therefore dominated by noise at precisely the frequencies that encode the bimodal structure.

The need for the adaptive bandwidth can also be understood from this perspective. With a global bandwidth, the kernel regression averages over movers whose deconvolution problems have very different conditioning: movers with large $|X_i|$ produce well-conditioned ratios (since $|u/X_i|$ is small, keeping $\hat\varphi_D$ bounded away from zero), while movers with small $|X_i|$ produce ratios dominated by noise. A single bandwidth that stabilizes the latter necessarily oversmooths the former, attenuating the high-frequency content of $\hat m_N(u)$.

The $k$-NN rule addresses this trade-off indirectly. Because movers with small $|X_i|$ concentrate near the threshold $\tau_x$, where the mover density is highest, they receive tight bandwidths under the $k$-NN rule; movers with large $|X_i|$, located in the sparser tails, receive wider bandwidths. The tight bandwidths in the dense region preserve the oscillatory structure of the conditional characteristic function that encodes bimodality, which a global bandwidth would smooth away. At the same time, the wider bandwidths in the tails provide additional smoothing for movers that, while well-conditioned for deconvolution, are few in number and would otherwise contribute noisy numerator estimates. Thus, the $k$-NN rule primarily improves resolution where data are plentiful, rather than stabilizing the noisiest deconvolution ratios directly. Stabilization of the small-$|X_i|$ ratios is instead provided by the denominator trimming threshold $\tau_{\mathrm{den}}$ and the averaging over movers in (ref).

Specification (d) isolates the role of the smoothness class of $D_i$. Since the Laplace characteristic function decays polynomially, $\varphi_D(t) = O(|t|^{-4})$ remains well above the trimming threshold at the moderate-to-high frequencies where the bimodal information resides: at the frequency argument $v = u/X_i = 3$, we have $|\varphi_D(v)| \approx 0.07$ under Laplace errors, compared to $8\times 10^{-5}$ under Gaussian errors---a ratio of approximately $865$. All tuning parameters are as in specification (c), and cross-validation again selects the tuning parameters jointly. Across $100$ replications, the procedure selects $S=7$ in $91$ replications (with occasional values of $9$ or $11$) and $k=30$ in $84$ replications, with approximately $30.7\%$ of observations classified as stayers.

Panel (d) of Figure (ref) displays the results under specification (d). The bimodality is recovered substantially better: the pointwise median tracks both modes and the trough between them, and the interquartile range is markedly tighter than under specification (c). The comparison between panels (c) and (d) is particularly clean because the cross-validation procedure selects essentially the same tuning parameters in both specifications---$k=30$ and $S=7$ in the large majority of replications---yet produces strikingly different outcomes. Since $f_\beta$ is identical in both specifications, this confirms that the attenuation observed under specification (c) is a consequence of the exponential decay of $\varphi_D$---a well-known limitation of deconvolution with supersmooth noise Fan1991---rather than a deficiency of the sieve estimator or the adaptive bandwidth procedure. Under Gaussian errors, the deconvolution ratio $\hat\varphi_{\tTwo Y\mid X}/\hat\varphi_D$ is dominated by noise at moderate-to-high frequencies, so the CV-selected $S=7$ reflects the highest sieve dimension that avoids fitting this noise. Under Laplace errors, the polynomial decay of $\varphi_D$ preserves the high-frequency content of $\widehat m_N(u)$ that encodes the bimodal structure, so the same sieve dimension $S=7$ now captures genuine features of $f_\beta$.

Application

We revisit the Nicaraguan calorie-demand application in GrahamPowell2012. The data consist of a balanced panel of $N=1{,}358$ poor rural households observed in 2000, 2001, and 2002 in communities covered by the conditional cash transfer program Red de Protecci\'on Social (RPS). Data construction is described in GrahamPowell2012 and the references therein.

Model and object of interest

Let \(r \in \{0,1\}\) denote the RPS regime, where \(r=1\) corresponds to assignment to receipt of RPS transfers and \(r=0\) corresponds to the no-RPS regime. For each household \(i\), define the potential structural calorie equation

where \(x\) denotes a counterfactual value of log real per-capita total household expenditure. The coefficient \(\beta_i(r)\) is the household-specific elasticity of calorie intake with respect to expenditure under regime \(r\). It is a causal derivative of the potential outcome schedule with respect to expenditure, holding the RPS regime fixed.

Within each RPS stratum, we suppress the regime argument, so that \(\beta_i\) should be read as \(\beta_i(r)\) for households observed under regime \(r\). Consequently, when the estimator is applied separately to RPS recipients and non-recipients, the corresponding densities estimate \(f_{\beta(1)}\) and \(f_{\beta(0)}\), respectively.

The raw panel has three time periods, \(s=2000,2001,2002\). After fixed-effect elimination, for each adjacent pair of periods and within a fixed RPS regime \(r\),

where \(t=1\) corresponds to 2000--2001 and \(t=2\) to 2001--2002, with

Within each RPS stratum, we again suppress the regime argument and write \[ Y_{it}=X_{it}\beta_i+U_{it}. \] The regressor is the change in log expenditure between adjacent periods, so that \[ X_{it} \approx \frac{\text{Exp}_{is}-\text{Exp}_{i,s-1}} {\text{Exp}_{i,s-1}} \] represents the approximate percentage change in household expenditure. Each year pair, taken separately, is a scalar irregular design with \((T,p,q)=(1,1,0)\). Stacking the two differenced periods, \(Y_i=(Y_{i1},Y_{i2})'\), yields the regular design with \((T,p,q)=(2,1,0)\).

Our objects of interest are the regime-specific cross-sectional density functions \(f_{\beta(r)}\), \(r\in\{0,1\}\). These densities summarize heterogeneity in household-specific structural calorie-expenditure elasticities under each RPS regime. In the full sample, the corresponding object is the distribution of the observed-regime elasticity \[ \beta_i^{\mathrm{obs}} = \beta_i(R_i), \] where \(R_i\) denotes the observed RPS status of household \(i\). Thus the full-sample density is an observed-regime mixture, whereas the subsample densities estimate \(f_{\beta(1)}\) and \(f_{\beta(0)}\) separately.

We estimate the density of the observed-regime elasticity in the full sample and estimate the regime-specific densities separately for RPS recipients, \(\mathrm{RPS}=1\), and non-recipients, \(\mathrm{RPS}=0\). For the full sample, we report both the scalar irregular estimator, applied separately to the two adjacent year pairs, and the regular estimator, which pools the two differenced periods. For the subsample analysis, we apply the scalar irregular estimator separately within each RPS stratum.

Because RPS assignment was randomized at the community level and take-up was high, comparisons of the estimated densities across RPS strata can be interpreted as causal contrasts between the regime-specific distributions, subject to the maintained CRC assumptions. Formally, the experimental design identifies an assignment-induced contrast in structural elasticity distributions; high take-up supports interpreting this contrast as receipt-regime contrast.

Allowing \(\beta_i(r)\) to be correlated with the path of expenditure growth \(X_i(r)=(X_{i1}(r),X_{i2}(r))'\) is natural in this application. Transfer amounts varied with household composition, so households receiving larger transfers experienced larger expenditure changes MaluccioFlores2005. The same households plausibly faced greater unmet caloric needs and therefore had stronger incentives to allocate additional resources toward calorie-dense staple consumption. In addition, households closer to caloric subsistence likely had both larger expenditure shocks in proportional terms and larger marginal calorie responses to additional resources. These considerations motivate a correlated random coefficient specification rather than an exogenous random coefficient model.

Within each RPS regime, Assumption 1(ii) requires that the transitory component \(U_{it}(r)\) is independent of \((\beta_i(r),X_i(r))\), where \(X_i(r)=(X_{i1}(r),X_{i2}(r))'\) denotes the vector of expenditure changes under regime \(r\). In the realized stratum-specific equation, where the regime argument is suppressed, this restriction becomes the requirement that \(U_{it}\) is independent of \((\beta_i,X_i)\).

This restriction should be interpreted as a decomposition of the outcome into systematic and idiosyncratic components. The model allows for arbitrary dependence between the structural elasticity \(\beta_i(r)\) and the expenditure path \(X_i(r)\), so that systematic co-movement between expenditure dynamics and calorie demand may operate through heterogeneous behavioral responses. The disturbance \(U_{it}(r)\) is then the residual component that remains after accounting for this heterogeneous structure, and is interpreted as an idiosyncratic noise term.

Under this interpretation, Assumption 1(ii) rules out residual dependence between the unexplained component and the regressors within each RPS regime, but does not preclude economically meaningful dependence between expenditure and calorie demand more generally. Instead, it requires that such dependence be captured by the distribution of \(\beta_i(r)\), rather than by the disturbance \(U_{it}(r)\).

Data

Figure (ref) reports histograms of changes in log expenditure and log calorie intake over the two differenced periods. A central feature of the data for the scalar irregular design is that the distribution of $X_i=\Delta\log(\text{Exp}_i)$ is centered near zero, with substantial mass in a neighborhood of the origin. This is important because first-step identification in the scalar irregular design relies on stayers, that is, households with $|X_i|$ close to zero. Beyond this concentration near zero, the expenditure-change distributions are heavy-tailed. The distributions of calorie changes are unimodal and fat-tailed.

Estimation

We implement both the scalar irregular estimator, which treats the two adjacent year pairs separately as $(T, p, q) = (1, 1, 0)$ designs, and the regular estimator, which pools them as a $(T, p, q) = (2, 1, 0)$ design.

\paragraph{Scalar irregular estimator.} The scalar irregular estimator exploits the presence of stayers, or households whose expenditure changed by less than a threshold $\tau_x$ in absolute value, to identify the characteristic function of the transitory component $D_i$ along a ray. This characteristic function is then used to deconvolve the distribution of $\beta_i$ from the conditional characteristic function of the transformed outcome for movers ($|X_i| \geq \tau_x$).

In this application, the stayer condition $|X_{it}| < \tau_x$ has a direct economic interpretation: households are classified as stayers if their expenditure changed by less than approximately $100 \times \tau_x$ percent between adjacent survey waves. The threshold $\tau_x$ is set according to a data-driven rule that scales with sample size and the dispersion of $X_i$; the resulting stayer proportions range from 7 to 10 percent across specifications, as reported in Table (ref).

\paragraph{Regular estimator.} When both differenced periods are stacked, the regressor matrix $X_i = (X_{i1}, X_{i2})' \in \mathbb{R}^2$ has $T = 2 > p = 1$, placing the problem in the regular design. The first-step transformation annihilates $\beta_i$ for every observation, yielding a scalar equation from which the characteristic function $\varphi_D$ is estimated by smoothing over normalized directions on $\mathbb{S}^1$. The second step recovers $f_\beta$ by deconvolution via the sieve minimum distance procedure.

\paragraph{Tuning parameter selection.} In all cases, tuning parameters are selected by the cross-validation procedure of Algorithm (ref). The selected values are summarized in Table (ref). Details on the implementation, including the formula for $\tau_x$ and the cross-validation grid, can be found in Section (ref).

\paragraph{Inference.} Sampling uncertainty is assessed using a nonparametric pairs bootstrap with $B = 499$ replications. In the scalar irregular design, each bootstrap draw resamples the differenced observations corresponding to the relevant year pair and recomputes the estimator. The stayer threshold $\tau_x$ and the first-stage bandwidth $h_0$ are held fixed at their original-sample values, while the second-stage tuning parameters $(h_x, S)$ are also fixed at their original cross-validation choices. In the regular design, each bootstrap draw resamples the household-level observations $(Y_{i1}, Y_{i2}, X_{i1}, X_{i2})$ and recomputes the estimator with the first-stage tuning parameter $h_S$ fixed at its original value and $(h_X, S)$ fixed at their original cross-validation selections.

Confidence intervals for the density are constructed using the basic (reverse-percentile) bootstrap.\footnote{We emphasize that this procedure is a practical approximation. Theoretical validity of the bootstrap for deconvolution density estimators has not been established in the literature, and extends a fortiori to the two-stage sieve minimum distance estimator with post-processing employed here. The reported confidence bands should therefore be interpreted as indicative measures of sampling variability rather than as intervals with guaranteed asymptotic coverage.} For each evaluation point $b$, let $\hat{f}_\beta(b)$ denote the estimator of $f_\beta(b)$ and let $\{\hat{f}_\beta^{*(r)}(b)\}_{r=1}^B$ denote the bootstrap draws. The $(1 - \alpha)$ confidence interval for $f_\beta(b)$ is given by \[ \left[\, \hat{f}_\beta(b) - q_{1-\alpha/2}(b), \; \hat{f}_\beta(b) - q_{\alpha/2}(b) \,\right], \] where $q_p(b)$ denotes the $p$-th quantile of the centered bootstrap distribution $\hat{f}_\beta^*(b) - \hat{f}_\beta(b)$. The reported confidence bands are obtained by applying this procedure pointwise over $b$ and therefore do not provide uniform coverage. Inference is conducted conditional on the selected tuning parameters and does not account for additional uncertainty arising from their estimation.

Results

\paragraph{Moments.} Table (ref) reports estimates of the mean and variance of the cross-sectional distribution of $\beta_i$, together with bootstrap standard errors, for the full sample and by RPS status.

Across specifications, the mean elasticity is positive and statistically significant, consistent with GrahamPowell2012. In the full sample, the regular estimator yields $\widehat{\mathbb{E}}[\beta_i] = 0.677$ (s.e.\ $0.035$), which lies between the irregular estimates for 2000--2001 ($0.684$) and 2001--2002 ($0.467$) but is much closer to the former. The regular estimator is roughly twice as precise as either irregular estimate, reflecting both pooling across periods and the imposition of time-invariance: rather than identifying $f_\beta$ separately from each pair of differenced periods using only stayer households for the denominator characteristic function, the regular estimator uses all $N$ households and the full collinearity structure of $X_i'X_i$, eliminating the leading source of variance in the irregular design.

All specifications exhibit substantial dispersion. In the full sample, the regular estimator gives $\widehat{\mathrm{Var}}[\beta_i] = 0.591$ (s.e.\ $0.017$), while the irregular estimates yield $0.765$ and $0.535$, albeit with markedly larger sampling variability (s.e.\ $0.188$ and $0.224$, respectively). The implied standard deviations range from $0.73$ to $0.88$ across designs, indicating considerable heterogeneity in marginal calorie responses across households, including within each RPS stratum.

The irregular estimates suggest an intertemporal shift concentrated among recipients. In the full sample, the mean declines from $0.684$ to $0.467$, but this difference is not statistically significant at conventional levels. Among non-recipients, the mean declines modestly from $0.726$ (s.e.\ $0.161$) to $0.606$ (s.e.\ $0.082$), a change that is also not statistically significant. In contrast, for RPS recipients, the mean falls sharply from $0.787$ (s.e.\ $0.086$) to $0.518$ (s.e.\ $0.099$), a difference that is statistically significant.

Across strata, the regular estimator places the recipient mean above the non-recipient mean ($\widehat{\mathbb{E}}[\beta_i \mid \mathrm{RPS}=1] = 0.698$ vs.\ $0.625$ for $\mathrm{RPS}=0$), with a difference of $0.073$ that is not statistically significant (s.e.\ of the difference $\approx 0.065$ under the standard independent-subsamples approximation, since the two pairs bootstraps draw from disjoint household sets). The regular variance estimates are $\widehat{\mathrm{Var}}[\beta_i \mid \mathrm{RPS}=1] = 0.635$ (s.e.\ $0.027$) and $\widehat{\mathrm{Var}}[\beta_i \mid \mathrm{RPS}=0] = 0.586$ (s.e.\ $0.028$); the cross-stratum difference of $0.049$ is likewise within one standard error of zero. Thus, while the point estimates suggest that recipients have a slightly higher mean and a slightly more dispersed distribution of marginal calorie responses than non-recipients, the regular estimator does not detect a statistically distinguishable shift in either moment.

\paragraph{Densities.} Figures (ref) and (ref) report the corresponding density estimates.

Irregular estimates. In the full sample (panel c), both period-specific densities are unimodal with modes near $0.5$, and exhibit substantial overlap, consistent with the failure to reject time-invariance at conventional significance levels. The 2001--2002 density is modestly shifted leftward relative to the 2000--2001 density. This shift reflects a change in location of the distribution, with the mode and central mass shifting leftward, rather than changes confined to the tails.

Among recipients (panel b), the shift is more pronounced: the 2000--2001 density has a mode near $0.6$, while the 2001--2002 density shifts leftward with a mode near $0.4$--$0.5$. The bootstrap bands for the two periods show less overlap than in the full sample, consistent with the statistically significant decline in the mean. Among non-recipients (panel a), the densities are similar across periods, with substantially overlapping bootstrap bands, consistent with time-invariant elasticities in this group.

Regular estimates. The regular estimator yields a unimodal density centered near $0.5$ with substantially narrower bootstrap bands than the irregular estimates, reflecting the efficiency gains from pooling and from imposing time-invariance. Comparing strata (Figure (ref)), the two density estimates are remarkably close: both peak in the interval $[0.4, 0.6]$ at heights near $0.50$, and the $90\%$ pointwise bootstrap bands overlap across essentially the entire support. The recipient density has slightly more mass in the right tail and slightly less mass at the mode, consistent with its larger mean and variance reported in Table (ref), but neither feature is large relative to sampling uncertainty.

The regular design provides little evidence of a large regime-level shift in the pooled elasticity distribution. Its main implication is instead that calorie--expenditure elasticities are heterogeneous. In light of the period-specific irregular estimates, however, the regular estimator should be interpreted as a common-\(\beta_i\) benchmark: if elasticities change with cumulative RPS exposure, the regular estimator necessarily averages over, and may mask, those exposure-specific changes.

\paragraph{Reconciling the literature.} The estimated density of the calorie--expenditure elasticity $\beta_i$ reveals substantial heterogeneity in household responses to changes in expenditure. While the average elasticity is positive and of the same order as those commonly reported in the literature, the distribution exhibits significant mass near zero as well as non-negligible probability on negative values.

This pattern helps reconcile seemingly conflicting findings in the development literature. Early studies document small calorie--expenditure elasticities, often close to zero BehrmanDeolalikar1987,BouisHaddad1992,Ravallion1990, whereas SubramanianDeaton1996 report elasticities in the range of 0.3--0.5 and argue against the view that calorie responses are negligible. Our results suggest that these differences may reflect aggregation: the mean elasticity masks substantial dispersion in individual responses.

A natural interpretation of the heterogeneity in $\beta_i$ is that households adjust along both quantity and quality margins. For some households, particularly those facing caloric constraints, higher expenditure translates into increased calorie intake, corresponding to positive values of $\beta_i$. For others, however, expenditure increases are primarily allocated toward higher-quality, more diverse foods rather than additional calories. This compositional shift, which is well documented in the literature Deaton1997,JensenMiller2008,Skoufias2011, weakens the relationship between income and calorie intake and can generate elasticities close to zero or even negative.

The presence of negative support in the estimated density indicates that, for a subset of households, increases in expenditure are associated with reductions in calorie intake, consistent with substitution away from inexpensive, calorie-dense staples toward foods with lower caloric density per unit of expenditure. Such behavior is also consistent with broader evidence on dietary transitions and changing consumption patterns as living standards improve DeatonDreze2009.

Overall, the distributional evidence underscores the limitations of focusing exclusively on average elasticities. By recovering the full distribution of $\beta_i$, the analysis reveals that the relationship between expenditure and nutrition is fundamentally heterogeneous, possibly reflecting differences in constraints, preferences, and consumption margins across households.

\FloatBarrier

Conclusion

This paper develops identification and estimation results for the distribution of correlated random coefficients in short panel data models, without imposing restrictions on serial dependence in the error terms. We show that the structure of the panel model leads to two distinct identification strategies, depending on whether the design is regular or irregular, and we provide corresponding estimators based on characteristic-function methods. In irregular designs, identification relies on limit arguments near singular regressor realizations, while in regular designs it follows from orthogonal projections. In both cases, estimation proceeds via a regularized deconvolution step combined with a Hermite sieve minimum-distance approximation.

Our analysis highlights that the statistical difficulty of the problem is driven by the first-stage deconvolution, which involves estimating a characteristic function at random, regressor-dependent frequencies and controlling instability arising from near-zero denominators. The resulting estimator exhibits a nonstandard bias--variance tradeoff, governed by the interaction between the stayer approximation, Fourier-domain regularization, and sieve approximation.

Our empirical application to the elasticity of calorie expenditure contributes to an ongoing debate in the development economics literature on whether households are calorie constrained or instead engage in quality substitution as income rises. By recovering the full distribution of heterogeneous elasticities, rather than focusing on average effects alone, our framework allows for a richer characterization of household behavior. In particular, it permits distinguishing between populations for which calorie intake responds strongly to expenditure changes and those for which additional resources are allocated toward higher-quality, non-caloric food consumption. This distributional perspective complements existing approaches and provides new evidence on the extent and nature of heterogeneity underlying the aggregate elasticity estimates commonly reported in the literature.

An important direction for future work is the development of inference procedures that accommodate both generated moment conditions and random frequency dispersion. In our setting, moment conditions naturally depend on first-step estimates of the disturbance characteristic function, introducing additional sampling variation and bias. At the same time, the relevant characteristic functions are evaluated at effective frequencies of the form $u/X_i$, so that information is dispersed according to the distribution of $1/X_i$. This dispersion governs the stability and informativeness of the empirical moments and complicates inference relative to standard deconvolution settings. Designing inference procedures that remain robust to these features is a promising avenue for further research.