Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
75,913 characters · 15 sections · 51 citation commands
Estimation of Latent Group Structures in Time-Varying Panel Data Models
\pagenumbering{Roman}
\paragraph{Abstract.} We consider panel data models where coefficients change smoothly over time and follow a latent group structure, being homogeneous within but heterogeneous across groups. To jointly estimate the group membership and group-specific coefficient trajectories, we propose FUSE-TIME, a pairwise adaptive group fused-Lasso estimator combined with polynomial spline sieves. We establish consistency, derive the asymptotic distributions of the penalized sieve estimator and its post-selection version, and show oracle efficiency. Monte Carlo experiments demonstrate strong finite-sample performance in terms of estimation accuracy and group identification. An application to the CO\textsubscript{2} intensity of GDP highlights the relevance of addressing both cross-sectional heterogeneity and time-variance in empirical exercises.
\paragraph{Keywords.} Panel data model, time-varying coefficients, clustering, basis splines.
\paragraph{JEL classification.} C22, C23, C38.
\pagenumbering{arabic} \doublespacing
In this paper, we study panel data models in which slope coefficients evolve smoothly over time and vary across an unobserved group structure, being homogeneous within groups but heterogeneous across groups. We propose a penalized sieve approach that leverages the pairwise adaptive group fused-Lasso (PAGFL; qian2016shrinkage,qian2016shrinkageET; mehrabani2023) to simultaneously estimate the time-varying coefficients and the latent grouping. We term our resulting estimator FUSE-TIME--Fused Unobserved group Spline Estimation of TIME varying coefficients--and derive its asymptotic properties when the number of groups, the cross-sectional dimension, and the time dimension jointly diverge: establishing estimation consistency, clustering consistency, and deriving the asymptotic distribution.
This paper builds on and extends two strands of the literature. The first strand studies panel data models with latent group structures, where coefficients are homogeneous within groups but vary across groups. Prior studies recover the unobserved groupings either by applying clustering algorithms such as $K$-means to preliminary estimates or by directly embedding the grouping into the estimation procedure through penalization schemes bonhomme2015grouped,ke2015homogeneity,su2016identifying,mehrabani2023,lumsdaine2023estimation,mehrabani2025. The second strand develops time-varying panel models, which describe temporal heterogeneity either through smoothly changing coefficients huang2004polynomial,cai2007trending,su2012sieve,robinson2012nonparametric or through discrete structural breaks bai2010common,qian2016shrinkage.
We contribute to these strands by proposing a single convex optimization problem that jointly identifies the latent group structure and estimates the smooth group-specific coefficient trajectories. Our approach does not require any prior specification of the number of groups or group compositions. Instead, the number of groups and their compositions are determined automatically in a data-driven way through a penalized sieve estimation procedure that adds a PAGFL penalty qian2016shrinkage,mehrabani2023 to the objective function. The penalty collects all pairwise differences of individual coefficient vectors, shrinking them to zero if two units belong to the same group. Following huang2004polynomial,su2019sieve, we approximate the time-varying coefficients using polynomial basis splines (B-splines), providing a flexible yet parsimonious representation of smooth temporal change. We propose a post-selection variant of our estimator, and introduce a consistent BIC-type information criterion (IC) to select the tuning parameter that regulates the number of groups and their compositions.
We provide primitive conditions for both the penalized sieve estimator and the post-selection variant of FUSE-TIME to achieve the oracle property, being asymptotically equivalent to an infeasible oracle estimator based on the true grouping. Our theoretical framework allows the number of groups, the cross-sectional dimension, and the time dimension to diverge jointly when deriving the asymptotic distributions of the estimators. FUSE-TIME thus provides a general framework for analyzing panels with both cross-sectional and temporal heterogeneity that is also robust to missing values and strikes a good balance between capturing heterogeneity and maintaining parsimony.
We also study the finite-sample performance of our estimator using both simulated and real data. Overall, our proposed estimator performs well across a variety of simulation settings in terms of recovering the latent group structure and estimating the time-varying coefficients accurately. FUSE-TIME outperforms the method of su2019sieve, which employs the Classifier Lasso (C-Lasso; su2016identifying), while also being computationally more efficient. Additionally, we illustrate the merits of the FUSE-TIME estimator by analyzing trends in the carbon dioxide (CO\textsubscript{2}) intensity of Gross Domestic Product (GDP), the amount of CO\textsubscript{2} emitted per unit of GDP produced. In a panel dataset of 92 countries spanning 1960 to 2023, the model identifies five groups, each with a unique trend of CO\textsubscript{2} intensity.
The plan for the remainder of this paper is as follows. Section (ref) introduces the model along with our estimation procedure. The asymptotic theory is developed in Section (ref). Section (ref) investigates the finite-sample performance via a simulation study. Section (ref) presents the empirical illustration. The \hyperref[sec:Conclusion]{final section} concludes. We provide a software implementation in our companion R-package PAGFL haimerl2025pagfl. Replication files of the simulation study and empirical application are available at \href{https://github.com/Paul-Haimerl/replication-tv-pagfl}{GitHub}.\footnote{R-Notebooks and a Dockerfile with replication material are published at \href{https://github.com/Paul-Haimerl/replication-tv-pagfl}{github.com/Paul-Haimerl/replication-fuse-time}.}
Notation. Throughout, we denote vectors by boldface small letters and matrices by boldface capital letters. For a real matrix $\boldsymbol{A}$, the Frobenius norm is written as $\|\boldsymbol{A}\|_F=\sqrt{\text{tr}(\boldsymbol{A}\boldsymbol{A}^\prime)}$. $\mu_{\max}(\boldsymbol{A})$ and $\mu_{\min}(\boldsymbol{A})$ denote the largest and smallest eigenvalues, respectively. $\| \boldsymbol{A} \|_{\text{sp}}= \sqrt{\mu_{\max} (\boldsymbol{A}\boldsymbol{A}^\prime)}$ gives the spectral norm. $\text{vec}(\boldsymbol{A})$ describes a column-wise vectorization of $\boldsymbol{A}$. $\|\boldsymbol{a}\|_2$ denotes the Euclidean norm for a real vector $\boldsymbol{a}$. The operators $\overset{P}{\to}$, $\overset{D}{\to}$, and plim signal convergence in probability, convergence in distribution, and the probability limit. $\otimes$ represents the Kronecker product. The superscript zero, or index when appropriate, marks a true quantity. $\boldsymbol{I}_a$ and $\boldsymbol{0}_a$ denote an $a \times a$ identity matrix and an $a \times 1$ vector of zeros. `With probability approaching one' is abbreviated as w.p.a.1. For two sequences of positive (random) numbers $a_n$ and $b_n$, $a_n \lesssim b_n$ indicates that $a_n / b_n$ is (stochastically) bounded and $a_n \asymp b_n$ signals that both $a_n \lesssim b_n$ and $b_n \lesssim a_n$ hold.
In this section, we introduce the time-varying panel data model (Section (ref)) and our penalized sieve estimation procedure (Section (ref)).
Consider the time-varying panel data model
where $y_{it}$ is the scalar response, $\boldsymbol{x}_{it}$ is a $p \times 1$ vector of explanatory variables, the $\gamma_i^0$'s denote the unobserved individual fixed effects, and $\epsilon_{it}$ represents a zero-mean idiosyncratic error. The superscript 0 is used to denote a true parameter value. The fixed effects $\gamma_i^0$ may correlate with some elements of $\boldsymbol{x}_{it}$ and are independently distributed across individuals. We assume that the $p$-dimensional vector of slope coefficients $\boldsymbol{\beta}_{it}^0$ varies smoothly over time:
Furthermore, we assume that the time-varying slope parameters adhere to the following unknown group structure
where $\boldsymbol{\alpha}_k^0 (t/T)$ is a $p \times 1$ vector of group-specific time-varying functional coefficients and $ \boldsymbol{1}\{\cdot\}$ denotes the indicator function. The latent group structure $\mathcal{G}_{K}^0=\{G_1^0, \dots, G_{K_0}^0\}$ partitions the cross-sectional units into disjoint sets; $\cup_{k=1}^{K_0} G_k^0 = \{1, \dots, N\}$ and $G_j^0 \cap G_k^0 = \emptyset$ for any $j \neq k$. We denote the cardinality of group $G_k^0$ as $N_k$.
Similarly to huang2004polynomial,su2019sieve, we approximate the $p \times 1$ functional coefficient vector $\boldsymbol{\beta}_{i}^0 (t/T)$ using an $M$-dimensional vector $\boldsymbol{b}(t/T)$ of polynomial spline basis functions:
where $\boldsymbol{\Pi}^0_i = (\boldsymbol{\pi}^0_{i1}, \dots, \boldsymbol{\pi}^0_{ip})$ denotes a $M \times p$ matrix of B-spline control points, and $\boldsymbol{\eta}_{it} = \boldsymbol{\beta}_i^0 (t/T) - \boldsymbol{\Pi}^{0 \prime}_i \boldsymbol{b} (t/T)$ is a $p \times 1$ sieve approximation error. B-splines offer two decisive advantages over typically employed kernel estimators. First, spline functions are computationally efficient and numerically stable while maintaining good approximation properties. Second, B-splines can be expressed as linear combinations of the basis functions $\boldsymbol{b}(t/T)$. The possibility of separating time-varying basis functions $\boldsymbol{b}(t/T)$, which do not have to be estimated, from individual time-constant control points $\boldsymbol{\Pi}^0_i$, which are estimated, makes B-splines particularly convenient to use in conjunction with the PAGFL penalty.
Note that $\boldsymbol{\Pi}_i^{0 \prime} \boldsymbol{b}(t/T)$ is a function in the sieve-space $\mathbb{B}_M$ spanned by the $M$ basis functions in $\boldsymbol{b}(t/T)$. By increasing the polynomial degree $d$ (polynomial order $d+1$) and the number of interior knots $M^*$ as $T$ diverges, we extend $M = M^* + d + 1$ along with $\mathbb{B}_M$ and obtain arbitrarily close approximations of $\boldsymbol{\beta}_{i}^{0} (t/T)$. Appendix (ref) contains full details on the spline basis functions.
Plugging (ref) into model (ref) yields
where $\boldsymbol{\pi}_i^0 = \text{vec} (\boldsymbol{\Pi}^{0}_i)$ is a $Mp \times 1$ coefficient vector, $\boldsymbol{z}_{it} = \boldsymbol{x}_{it} \otimes \boldsymbol{b}(t/T)$ is a $Mp \times 1$ vector of regressors, and $u_{it}$ collects the idiosyncratic error $\epsilon_{it}$ and the sieve approximation error $\boldsymbol{\eta}_{it}^\prime \boldsymbol{x}_{it}$. To obtain an estimate of the time-varying slope parameters $\boldsymbol{\beta}_{i}^0 (t/T)$ in model (ref), we estimate the time-constant vector $\boldsymbol{\pi}_i^{0}$ in model (ref) and take $\hat{\boldsymbol{\beta}}_{i} (t/T) = \hat{\boldsymbol{\Pi}}_i^\prime \boldsymbol{b}(t/T)$, where $\hat{\boldsymbol{\pi}}_i = \text{vec}(\hat{\boldsymbol{\Pi}}_i)$.
We jointly identify the latent group structure in (ref) and the time-varying coefficients $\boldsymbol{\beta}_{i}(t/T)$ through {a penalized sieve technique}, leveraging the PAGFL penalty qian2016shrinkage,mehrabani2023. In what follows, we introduce this estimation procedure in the context of the time-varying panel data model (ref). Beyond this baseline, our proposal can readily accommodate extensions, as discussed in Appendix (ref).
First, we concentrate out the individual fixed effects $\gamma_i^0$ in (ref) using $\tilde{y}_{it} = y_{it} - T^{-1} \sum_{t = 1}^{T} y_{it}$, with $\tilde{\boldsymbol{z}}_{it}$ defined analogously. Let $\boldsymbol{\pi} = (\boldsymbol{\pi}_1^\prime, \dots, \boldsymbol{\pi}_N^\prime)^\prime$ collect all individual parameters. We then take as objective function
where the first part is the usual sum of squared residuals and the second part is the PAGFL penalty term with tuning parameter $\lambda>0$ and adaptive penalty weights $\dot{\omega}_{ij}$. The penalty encourages sparsity in the differences of all $N(N-1)/2$ pairs of individual coefficient vectors. For large $\lambda$, some of these differences will be shrunken to exactly zero, implying that the corresponding cross-sectional units feature identical parameter estimates. The weights are set to $\dot{\omega}_{ij} = \| \dot{\boldsymbol{\pi}}_i - \dot{\boldsymbol{\pi}}_j \|_2^{-\kappa}$, where $\dot{\boldsymbol{\pi}}_i$ represents an initial consistent estimate $\dot{\boldsymbol{\pi}} = \arg \min_{\pi} T^{-1} \sum_{i = 1}^{N} \sum_{t=1}^{T} \left(\tilde{y}_{it} - \boldsymbol{\pi}_i^\prime \tilde{\boldsymbol{z}}_{it} \right)^2$ and $\kappa$ is specified by the user; for simplicity and in line with the standard choice in the adaptive Lasso literature, we maintain $\kappa=2$ throughout the paper.
Objective (ref) generalizes the objective function in mehrabani2023 to our FUSE-TIME approach. The penalized sieve estimator (PSE) $\hat{\boldsymbol{\pi}}$ of FUSE-TIME is then obtained by minimizing (ref), namely $$ \hat{\boldsymbol{\pi}}_\lambda = \arg \min_{\pi} \mathcal{F}_{NT}(\boldsymbol{\pi}, \lambda). $$ Cross-sectional units with identical slope estimates are assigned to the same cluster by collecting all $\hat{K}$ unique subvectors $\hat{\boldsymbol{\pi}}_i$ of $\hat{\boldsymbol{\pi}}_\lambda$ in the vector $\hat{\boldsymbol{\xi}} = (\hat{\boldsymbol{\xi}}_1^\prime, \dots, \hat{\boldsymbol{\xi}}_{\hat{K}}^\prime)^\prime$ and defining the set $\hat{G}_k = \{ i: \hat{\boldsymbol{\pi}}_i = \hat{\boldsymbol{\xi}}_k, \; 1 \leq i \leq N \}$ for each $k = 1, \dots, \hat{K}$. Subsequently, the group-specific PSE coefficient function equals $\hat{\boldsymbol{\alpha}}_{k} (t/T) = \hat{\boldsymbol{\Xi}}^\prime_k \boldsymbol{b}(t/T)$, where $\hat{\boldsymbol{\xi}}_k = \text{vec}(\hat{\boldsymbol{\Xi}}_k)$. As such, parameter estimation and identification of the latent group structure are performed simultaneously. The total number of clusters $\hat{K}$ is determined by the tuning parameter $\lambda$. We propose to select $\lambda$ using the consistent information criterion (IC) introduced in Section (ref). In the following, unless required, we suppress the dependence of $\hat{\boldsymbol{\pi}}$ on $\lambda$ to lighten the notation.
Objective function (ref) is convex in $\boldsymbol{\pi}$. To solve for $\boldsymbol{\pi}$, we employ a computationally efficient Alternating Direction Method of Multipliers (ADMM) algorithm, adapted from mehrabani2023 and detailed in Appendix (ref). An open-source implementation of this algorithm is provided in our companion R-package PAGFL haimerl2025pagfl.
Finally, to mitigate the bias coming from the fusion penalty term in (ref), we propose a post-selection fused-Lasso estimator of FUSE-TIME, labeled post-Lasso estimator for brevity, given the estimated group pattern $\hat{\mathcal{G}} = \{\hat{G_1}, \dots, \hat{G}_{\hat{K}} \}$:
Subsequently, the final post-Lasso time-varying coefficient estimates follow as $\hat{\boldsymbol{\alpha}}^p_{\hat{G}_k } (t/T) = \hat{\boldsymbol{\Xi}}^{p\prime}_{\hat{G}_k} \boldsymbol{b}(t/T)$, using $\hat{\boldsymbol{\xi}}^{p}_{\hat{G}_k} = \text{vec}(\hat{\boldsymbol{\Xi}}^{p}_{\hat{G}_k})$.
In this section, we study the consistency of the coefficient estimates and the clustering procedure (Sections (ref) and (ref)). Furthermore, we establish the limiting distribution of the PSE and the post-Lasso estimator (Section (ref)), and we propose a consistent BIC criterion to select the tuning parameter (Section (ref)). Formal proofs appear in Appendix (ref).
Let $\mathbb{C}^{\theta}[0,1]$ denote the space of functions that are $\theta$-times continuously differentiable on the unit interval. Moreover, let $\boldsymbol{x}_{it}^{(2)} = \boldsymbol{x}_{it}$ if no intercept is included in $\boldsymbol{x}_{it}$ and $\boldsymbol{x}_{it} = (1, \boldsymbol{x}_{it}^{(2)\prime})^\prime$ otherwise.
\phantomsection
(ref)(ref) is standard in the factor model literature and limits cross-sectional dependence bai2002determining. This condition is trivially satisfied if the idiosyncratic errors are independently distributed, an assumption frequently made for similar panel data models with latent group structures su2016identifying,qian2016shrinkage,su2019sieve.
(ref)(ref)-(ref) loosely follows su2019sieve. (ref)(ref) imposes a strong mixing process, which nests popular and extensively studied time series processes with geometrically decaying innovations, such as common ARMA, GARCH, Markov-switching, or threshold autoregressive models. Serial correlation as well as conditional heteroscedasticity in the error process are allowed. As su2019sieve point out, in dynamic panels Assumption \hyperref[line:A1]{1(ii)} implies the fixed effects $\gamma_i$ to be nonrandom. If the fixed effects are treated as random variables, the strong mixing condition in Assumption \hyperref[line:A1]{1(ii)} must hold conditional on each realized $\gamma_i$, with the distribution of $\gamma_i$ being bounded uniformly in $i$ hahn2011bias.
Assumptions \hyperref[line:A1]{1(i)-(ii)} can also be replaced by similar primitive bonhomme2015grouped or higher-level assumptions mehrabani2023 that limit the degree of dependence across time and the cross-section such that a central limit theorem can be applied.
Assumptions \hyperref[line:A1]{1(iii)-(iv)} place common moment conditions on the regressors and innovations. The conditions on the regressors in Assumptions \hyperref[line:A1]{1(iii)-(iv)} are redundant if $\boldsymbol{x}_{it}$ is deterministic. Assumption \hyperref[line:A1]{1(v)} allows group sizes to remain constant or to diverge at a speed slower than $N$. Assumption \hyperref[line:A1]{1(vi)} imposes smoothness on the true functional coefficients and ensures that they can be well approximated by splines of polynomial degree $d$.
Let $J_{\min}$ denote the minimum group separability in the sieve-space $\mathbb{B}_M$, $J_{\min} = \min_{i \in G_k^0, j \notin G_k^0} \| \boldsymbol{\pi}_i^0 - \boldsymbol{\pi}_j^0 \|_2$.
\phantomsection
Assumption \hyperref[line:A1]{2(i)} controls the size of the spline basis system. Assumption \hyperref[line:A2]{2(ii)} determines the rate at which the minimum group separability in $\mathbb{B}_{M}$ may shrink to zero. Assumption \hyperref[line:A2]{2(iii)} governs the speed at which the tuning parameter $\lambda$ must shrink to zero. Assumption \hyperref[line:A2]{2(iv)} places conditions on the relative rates of the number of coefficients to be estimated, coefficient convergence, and $\lambda$ so that group-specific trajectories can still be consistently estimated. In sum, $N$, $T$, $M$, and $K_0$ diverge to infinity, whereas $\lambda$ and $J_{\min}$ tend to zero in the limit.
Theorem \hyperref[sec:Theo_1]{3.1} establishes pointwise and mean-square convergence of $\hat{\boldsymbol{\pi}}_i$. The first term in the rates of Theorem \hyperref[sec:Theo_1]{3.1} reflects the stochastic error. The second term corresponds to the asymptotic bias of the sieve technique. Increasing the complexity of the spline system $M$ involves a bias-variance trade-off. On the one hand, the larger $M$, the slower the convergence due to the increased size of the coefficient vector.\footnote{Note that $O_p(MT^{-1/2} + M^{- \theta + 1/2}) = o_p(1)$ by Assumptions \hyperref[line:A2]{1(vi)} and \hyperref[line:A2]{2(i)}.} On the other, increasing $M$ reduces the asymptotic sieve bias. Moreover, the greater the order of continuous differentiability of the true coefficient functions $\theta$, the faster the sieve bias shrinks in $M$.
Corollary \hyperref[sec:Coro_1]{3.2} relates the results of Theorem \hyperref[sec:Theo_1]{3.1} to the actual functional coefficients. Note that the true coefficient functions are time-continuous. As a consequence, we report the supremum and integral over the unit interval. Interestingly, the pointwise rates of $\hat{\boldsymbol{\pi}}_i$ and $\hat{\boldsymbol{\beta}}_i(v)$ match, whereas the $L_2$ rate of $\hat{\boldsymbol{\beta}}_i(v)$ is faster in $M$. This result is due to the boundedness of the B-splines (see \hyperref[sec:Lemma_statements]{Lemma A.1(ii)}).
This subsection establishes clustering consistency.
Theorem \hyperref[sec:Theo_2]{3.3} establishes that, in the limit, all cross-sectional units belonging to the same group $G_k^0$ are jointly assigned to the same group; no other units are assigned to this group.
Based on Theorem \hyperref[sec:Theo_2]{3.3}, Corollary \hyperref[sec:Coro_2]{3.4} shows that the correct number of groups and group structure will be derived asymptotically. This result is intuitive since, given Theorem \hyperref[sec:Theo_2]{3.3}, $\hat{\boldsymbol{\pi}}$ can only hold $K_0$ distinct individual subvectors and all homogeneous individuals are assigned to the same group as $(N,T) \to \infty$. Since the latent grouping will be identified w.p.a.1, Corollary \hyperref[sec:Coro_2]{3.4} motivates the oracle property of the procedure. This implies that the estimation procedure is asymptotically equivalent to an infeasible oracle estimator based on the true grouping.
In the following, we derive the asymptotic distribution of the PSE and the post-Lasso.
\phantomsection
The penalty term in (ref) grows quadratically in $N$. Hence, Assumption \hyperref[line:A3]{3} strengthens Assumption \hyperref[line:A2]{2(iii)} and imposes conditions for the penalty term to vanish asymptotically, resulting in coinciding limiting distributions for the PSE and the post-Lasso. To this end, Assumption \hyperref[line:A3]{3} specifies a larger group separation or a faster rate at which $\lambda$ converges to zero. Assumption \hyperref[line:A3]{3} is only relevant for the PSE and need not hold for the asymptotic properties of the post-Lasso.
Let $\boldsymbol{\epsilon}_i = (\epsilon_{i1}, \dots, \epsilon_{iT})^\prime$. \phantomsection
Assumption \hyperref[line:A4]{4(i)} imposes a mild restriction on the error process, enabling the use of the Lindeberg condition to derive the limiting distribution in Theorem \hyperref[sec:Theo_3]{3.5}. Assumption \hyperref[line:A4]{4(ii)} provides an additional regularity condition that precludes the sieve approximation error from dominating the limiting distributions of the PSE and the post-Lasso. Let $\hat{\boldsymbol{\mathcal{Q}}}_{G_k^0, \tilde{z} \tilde{z}} = \sum_{i \in G_k^0} \hat{\boldsymbol{Q}}_{i, \tilde{z} \tilde{z}}$ with $\hat{\boldsymbol{Q}}_{i, \tilde{z} \tilde{z}} = T^{-1} \sum_{t = 1}^{T} \tilde{\boldsymbol{z}}_{it} \tilde{\boldsymbol{z}}_{it}^\prime$ and $\tilde{\boldsymbol{Z}}_i = (\tilde{\boldsymbol{z}}_{i1}, \dots, \tilde{\boldsymbol{z}}_{iT})$. Furthermore, define the variance-covariance matrix $\hat{\boldsymbol{\Omega}}_{G_k^0} = \hat{\boldsymbol{\nu}}_{G_k^0}^{\prime} \hat{\boldsymbol{\mathcal{E}}}_{G_k^0} \hat{\boldsymbol{\nu}}_{G_k^0}$ with $$ \hat{\boldsymbol{\nu}}_{G_k^0} = \left( M N_k^{-1} \hat{\boldsymbol{\mathcal{Q}}}_{G_k^0, \tilde{z} \tilde{z}} \right)^{-1} \left( \boldsymbol{I}_p \otimes \boldsymbol{b} (v) \right), $$ $$ \hat{\boldsymbol{\mathcal{E}}}_{G_k^0} = \frac{M}{N_k T} \sum_{i \in G_k^0} \tilde{\boldsymbol{Z}}_{i}^\prime E(\boldsymbol{\epsilon}_i \boldsymbol{\epsilon}_i^\prime) \tilde{\boldsymbol{Z}}_{i}, $$ and $\boldsymbol{q}_{G_k^0} = \sqrt{M / (N_k T)} \sum_{i \in G_k^0} \sum_{t = 1}^{T} E(\tilde{\boldsymbol{z}}_{it} \tilde{\epsilon}_{it})$.
$\boldsymbol{q}_{G_k^0}$ equals zero in the case of strictly exogenous regressors. However, commonly referred to as the Nickell bias, $\boldsymbol{q}_{G_k^0}$ is nonzero and of order $O \left( \sqrt{N_k / T} \right)$ for dynamic panel data models nickell1981biases,phillips2007bias. The bias emerges when $T$ remains fixed or grows slower than $N_k$. The within-transformation nets out the fixed effect $\gamma_i$ but simultaneously induces a contemporaneous correlation between the error term and the autoregressive regressor, thus biasing the coefficient estimate if $T$ does not grow fast enough relative to the number of individuals that feed into the group-specific autoregressive coefficient function nickell1981biases,kiviet1995bias,hahn2011bias,su2016identifying.
The PSE and the post-Lasso estimators are asymptotically equivalent to the infeasible oracle estimator with true grouping structure, hence both achieve oracle efficiency. Additionally, the B-splines yield arbitrarily close approximation of the true coefficient functions $\boldsymbol{\alpha}_k^0 (t/T)$. Given Assumption \hyperref[line:A3]{3}, both estimators feature the same asymptotic distribution. Nevertheless, despite their equivalence in the limit, we recommend using the post-Lasso in finite-sample applications. The penalty term of the PSE can lead to non-negligible bias, particularly in small samples. The simulation study in Section (ref) corroborates this finding.
Given a consistent estimator $\hat{\boldsymbol{\mathcal{E}}}_{\hat{G}_k}$ and exploiting the oracle property in Corollary \hyperref[sec:Coro_2]{3.4}, the variance of the limiting distribution can be estimated consistently. Potential techniques to derive $\hat{\boldsymbol{\mathcal{E}}}_{\hat{G}_k}$ are numerous. We follow the literature on heteroscedasticity and autocorrelation robust estimation of covariance matrices and take
where $H_T$ is the window size, $\hat{\boldsymbol{\mathcal{E}}}_{\hat{G}_k}^{(h)} = N_k^{-1} \sum_{i \in \hat{G}_k} T^{-1} \sum_{t = h + 1}^{T} \boldsymbol{\tilde{z}}_{it} \tilde{z}_{it-1}^\prime \hat{\epsilon}_{it} \hat{\epsilon}_{i t-1}$, and $\hat{\epsilon}_{it} = \tilde{y}_{it} - \hat{\boldsymbol{\alpha}}_{\hat{G}_k}^{p \prime} (t/T) \tilde{\boldsymbol{x}}_{it}$ newey1987simple,hahn2002asymptotically,pesaran2006estimation,muller2014hac. $w(h, H_T)$ denotes a weighting function subject to $\sup_{h} | w(h, H_T) | < \infty$ and $\lim_{T \to \infty} w(h, H_T) = 1$. It is common to specify $w(h,H_T) = 1 - h / (H_T + 1) \boldsymbol{1} \{ h \leq H_T \}$ newey1987simple. A rigorous derivation of the consistency of $\hat{\boldsymbol{\mathcal{E}}}_{\hat{G}_k}$ and required primitive conditions are beyond the scope of this paper and are left for future work su2012sieve.
In the remainder of this section, we make the dependence of $\hat{K}_\lambda$ and $\hat{G}_{k,\lambda}$ on $\lambda$ explicit. We choose the tuning parameter $\lambda$ that minimizes the following BIC-type criterion
where $\hat{\sigma}^2_{\hat{\mathcal{G}}_{\hat{K}, \lambda}} = (NT)^{-1} \sum_{k = 1}^{\hat{K}_{\lambda}} \sum_{i \in \hat{G}_{k, \lambda}} \sum_{t = 1}^{T} \left(\tilde{y}_{it} - \hat{\boldsymbol{\alpha}}_{\hat{G}_{k, \lambda}}^{p \prime} (t/T) \tilde{\boldsymbol{x}}_{it} \right)^2$ and $\rho_{NT}$ represents a tuning parameter qian2016shrinkage,su2019sieve,mehrabani2023. Define the set $\Lambda = [0, \lambda_{\max}]$ for some sufficiently large $\lambda_{\max} < \infty$ and partition $\Lambda$ into $\Lambda_{0,NT}$, $\Lambda_{-, NT}$, and $\Lambda_{+, NT}$, such that $\Lambda_{0,NT} = \{ \lambda \in \Lambda: \hat{K}_{\lambda} = K_0 \}$, $\Lambda_{-,NT} = \{ \lambda \in \Lambda: \hat{K}_{\lambda} < K_0 \}$, and $\Lambda_{+,NT} = \{ \lambda \in \Lambda: \hat{K}_{\lambda} > K_0 \}$.\footnote{We index the three subsets of $\Lambda$ with $NT$ to make it apparent that $\hat{K}$ is obtained from a random sample of dimension $(N,T)$.} In addition, assume that all $\lambda \in \Lambda_0$ comply with the regularity conditions stated in Assumptions \hyperref[line:A2]{2}-\hyperref[line:A3]{3}. Denote the set of all $K$-partitions of $\{1, \dots, N\}$ as $\mathbb{G}$.
\phantomsection
Assumption \hyperref[line:A5]{5} is standard in the literature and implies that the mean squared error (MSE) of an underfitted model is asymptotically larger than $\sigma_0^2$, the MSE of the true model.
\phantomsection
Assumption \hyperref[line:A6]{6} collects conditions such that the penalty term in (ref) either dominates the MSE or vanishes in the limit, depending on whether $\hat{K} > K_0$ or $\hat{K} \leq K_0$, respectively.
From Theorem \hyperref[sec:Theo_4]{3.6} follows that neither an underfitted nor an overfitted model maximizes the IC as $N$ and $T$ diverge to infinity. Consequently, the IC in (ref) uncovers the true model asymptotically.
When applying the IC, a practitioner must select $\rho_{NT}$. After some preliminary experiments and in line with previous literature, we recommend specifying $\rho_{NT} = c_\lambda \log(NT) (NT)^{-1/2}$ with $c_\lambda = 0.04$ mehrabani2023.
This section examines the finite-sample properties of FUSE-TIME using several Monte Carlo simulation studies. Section (ref) describes the DGPs, followed by implementation details and evaluation metrics in Section (ref). In Section (ref), we present the results.
We consider three DGPs for all $(N,T)$ combinations of $N=\{50, 100\}$ and $T=\{50, 100\}$. The cross-sectional individuals are sampled from $K_0=3$ groups with the proportions $0.3$, $0.3$, and $0.4$. Throughout the simulation study, we draw the individual fixed effects $\gamma_i^0$ and the idiosyncratic innovations $\epsilon_{it}$ from mutually independent standard normal distributions. Following su2019sieve, the DGPs are constructed as follows: \\ DGP 1: Trending panel data model \\ $$y_{it} = \gamma_i^0 + \beta_{i,0}^0 \left( \frac{t}{T} \right) + \epsilon_{it},$$ with the coefficient functions $$ \beta_{i,0}^0(v)=
$$ $F(v; a, b) = \left[1 + \exp(- (v - a) / b) \right]^{-1} $ denotes a cumulative logistic function. \\ \textbf{DGP 2}: Trending panel with an exogenous regressor\\ $$y_{it} = \gamma_i^0 + \beta_{i,1}^0 \left( \frac{t}{T} \right) + \beta_{i,2}^0 \left( \frac{t}{T} \right) x_{it} + \epsilon_{it},$$ that augments DGP 1 with a scalar exogenous explanatory variable $x_{it}$ which is generated from a standard normal distribution. Let $\beta_{i,1}^0 (t/T) = 0.5 \beta_{i,0}^0(t/T)$ (and hence $\alpha^0_{i,1}(t/T) = 0.5 \alpha_{i,0}^0(t/T)$), where $\beta_{i,0}^0$ and $\alpha_{i,0}^0$ are as defined in DGP 1, and $$ \beta_{i,2}^0(v)=
$$ \\ \textbf{DGP 3}: Dynamic panel data model $$y_{it} = \gamma_i^0 + \beta_{i,3}^0 \left( \frac{t}{T} \right) y_{it-1} + \epsilon_{it},$$ featuring a time-varying autoregressive functional relationship. In order to comply with the strong mixing condition in Assumption \hyperref[line:A1]{A.1(ii)}, $\sup_{v \in [0,1]} | \alpha_{k,3}^0(v) | < 1$ for all $k = 1, \dots, K_0$, the autoregressive coefficient is taken as $$ \beta_{i,3}^0(v)=
$$ Figure (ref) presents the sample paths of the simulated coefficient functions across the three DGPs.
The FUSE-TIME estimator is obtained as described in Section (ref). We select the number of interior knots $M^*$ of the spline system $\boldsymbol{b}(t/T)$ by taking $M^* = \max \left\{ \left\lfloor (NT)^{1/7} - \log(p)\right\rfloor, 1 \right\} $, where $\left\lfloor \cdot \right\rfloor$ rounds to the lower integer. When $T$ is small, large $M^*$ tend to overfit the individual time-varying coefficients and pollute the grouping. Conversely, a small $M^*$ caps the flexibility of the time-varying coefficients and may complicate the group differentiation in large panels with many latent groups. Furthermore, we set $d=3$ since cubic splines offer a good trade-off between flexibility and parsimony. This gives a $M^* + d + 1 = M$-dimensional vector of spline basis functions.
We compute the FUSE-TIME estimator for a grid of $\lambda$ values; the grids are specified in Table (ref) of Appendix (ref) for the three DGPs. The IC in equation (ref) selects the $\lambda$ tuning parameter. We hereby set $\rho_{NT} = c_\lambda \log(NT) (NT)^{-1/2}$ with $c_\lambda = 0.04$ as described in Section (ref). In large samples, the results do not vary substantially in $c_\lambda$.
To evaluate performance, we inspect grouping and estimation accuracy. The clustering performance is evaluated across $n_{\text{sim}}$ Monte Carlo experiments according to the following criteria:
The root mean square error (RMSE) quantifies the estimation accuracy of the time-varying coefficient functions and is computed as $$RMSE \left(\hat{\alpha}_l \right) = \frac{1}{N} \sum_{i = 1}^{N} \sqrt{ \frac{1}{T} \sum_{t = 1}^{T} \left[ \hat{\alpha}_{i,l} \left( \frac{t}{T} \right) - \alpha_{i,l}^0 \left( \frac{t}{T} \right) \right]^2},$$ where $\hat{\alpha}_{i,l} (t/T) = \hat{\alpha}_{k,l} (t/T)$ if $i \in \hat{G}_k$ and $\alpha_{i,l}^0 (t/T) = \alpha_{j,l}^0 (t/T)$ if $i \in G_j^0$, for $l = \{0, 1, 2, 3 \}$.
When reporting results on estimation accuracy, note that we explicitly distinguish between the proposed PSE and the post-Lasso estimator as defined in Section (ref); whereas throughout the rest of the paper we refer to the latter as the FUSE-TIME estimator.Their performance is compared with two benchmarks. The first is the infeasible oracle estimator that is based on the true latent grouping $\mathcal{G}_0$, averaged across all $n_{\text{sim}}$ experiments. The second is the time-varying C-Lasso by su2019sieve. Appendix (ref) lists details on the implementation of the time-varying C-Lasso and Appendix (ref) benchmarks the computational cost of the FUSE-TIME and \textit{C-Lasso} software implementations; the computational gains of the former are sizable.
We simulate all DGPs $n_{\text{sim}} = 300$ times and apply FUSE-TIME as well as the time-varying C-Lasso procedures. Table (ref) reports metrics on the grouping performance; note that there is no distinction between PSE and post-Lasso in terms of grouping performance, hence we report the results of our proposal under FUSE-TIME in Table (ref).
FUSE-TIME convincingly outperforms the C-Lasso benchmark, particularly in smaller samples and DGPs 2 and 3. FUSE-TIME's clustering accuracy improves quickly with increasing $T$, but not with $N$. This is intuitive since the clustering routine is driven by the estimated individual control points $\hat{\boldsymbol{\pi}}_i$, which are not consistent with respect to the cross-sectional dimension (see Theorem \hyperref[sec:Theo_1]{3.1}). Accordingly, when $T$ is small, the estimation of individual coefficient functions is highly noisy, complicating correct group assignment. This property is particularly evident when comparing DGP 1 and DGP 2. The latter involves estimating a second regression curve, which introduces additional uncertainty and subsequently does not lead to the near-perfect clustering observed in DGP 1. Nevertheless, FUSE-TIME still shows strong clustering performance for DGP 2. In contrast, the results for the dynamic panel data model of DGP 3 are markedly worse than for the other DGPs, especially when $T = 50$. Similar estimation devices have previously reported decreased performance for dynamic panel data models mehrabani2023. Nevertheless, as indicated by the $ARI$ measure, even when individual cross-sectional units are misclassified, the estimated grouping remains a close approximation of the true unobserved group structure, with FUSE-TIME substantially outperforming the C-Lasso: ARIs in the range of 0.83-0.98 for FUSE-TIME versus 0.49-0.54 for C-Lasso.
Table (ref) reports the RMSE for each time-varying coefficient, and split among PSE and post-Lasso estimators. Figure (ref) provides the estimated sample paths of the functional coefficients for each DGP and $N,T = 50$.
The RMSE results largely align with the clustering performance. This also holds for the comparison with the time-varying C-Lasso, which returns considerably poorer results regarding the RMSE as well. Unlike clustering accuracy, the post-Lasso RMSE also improves with $N$ due to the larger number of cross-sectional units available for pooling when estimating group-specific trajectories. Notably, the post-Lasso performs well relative to the oracle estimator, even in cases, such as DGP 3, where the clustering results may suggest a poor fit. This finding is likely because misclassified units tend to feature sample paths that are particularly similar to other groups, leading to a minor impact on the \textit{RMSE} compared to the oracle estimation. Precise identification of group-specific coefficients despite several misclassifications has been previously reported by bonhomme2015grouped. Table (ref) further highlights that the shrinkage penalty leads to a sizeable finite-sample bias in the \textit{PSE}, even though the \textit{PSE} and the \textit{post-Lasso} share the same asymptotic distribution (see Theorem \hyperref[sec:Theo_3]{3.5}). This motivates our usage of the \textit{post-Lasso} for empirical applications. Finally, Figure (ref) corroborates the previous findings: while minor deviations occur, partly due to misclassified units, the estimated trajectories closely follow the true underlying function.
Additionally, we also simulate the three DGPs with either idiosyncratic errors following an $AR(1)$ process with lag polynomial $a(L) = 1 - 0.3L$ or unbalanced panels with missing data, where we generate sample paths that randomly discard 30% of the observations. Detailed results are available in Appendix (ref). For unbalanced panels with missing data (see Appendix (ref)), the simulation results largely mirror the ones presented here. Discarding part of the sample has a similar effect to decreasing $T$ in a balanced panel. However, when the innovations follow an $AR(1)$ process (see Appendix (ref)), the performance decreases markedly, particularly when $T$ is small. Introducing autocorrelation to the errors increases the estimation uncertainty and thus contaminates the clustering mechanism. However, even when the grouping is correctly estimated, the coefficient trajectories offer a significantly poorer fit than the baseline case.
Many major economies pledge to reduce the emission of harmful greenhouse gases, a key driver of global warming. Achieving this objective without impeding economic growth requires emissions per unit of GDP to decrease. CO\textsubscript{2} accounts for approximately 66% of the anthropogenic contribution to global warming and is a ubiquitous by-product of economic activity and energy production bennedsen2023neural. Consequently, understanding the relationship between CO\textsubscript{2} pollution and economic growth is crucial for effective policy and climate action. We contribute by applying our FUSE-TIME procedure to identify trends in the CO\textsubscript{2} emission intensity of GDP, the CO\textsubscript{2} emitted per unit of GDP produced.
The CO\textsubscript{2} intensity of GDP is determined by the structural composition of an economy and the “greenness" of its energy production bella2014relationship. For instance, transitioning from a production-based to a service-oriented economy reduces energy consumption and direct emissions from high-pollution industries. Likewise, shifting from emission-intensive solid fuels, such as coal or biomass, to cleaner carriers like natural gas or renewables, lowers the elasticity between energy demand and associated pollution. The combined effects of economic restructuring and energy source optimization are often credited with driving the decoupling of GDP growth from production-based emissions in high-income countries jakob2012will,wang2017multi. This decoupling pattern has inspired an extensive literature on the inverted U-shaped environmental Kuznets curve, which posits a positive elasticity between income and emissions during the early stages of economic development, turning negative as technological progress and greater prosperity take hold grossman1995economic,wagner2015environmental,bennedsen2023neural. However, as (log) GDP, an integrated process, clearly does not satisfy the regularity conditions in Assumption \hyperref[line:A1]{1}, we instead study the relationship between emissions and income by estimating trends in CO\textsubscript{2} intensity. Trend functions capture the long-term behavior of CO\textsubscript{2} intensity by smoothing over short-run fluctuations such as business cycles, fuel price volatility, or temporary energy supply disruptions. Moreover, it is well-established that economies exhibit considerable heterogeneity in their developmental trajectories kose2003international,azomahou2006economic. Consequently, assuming a common CO\textsubscript{2} intensity trend across many countries is implausible. Nevertheless, much of the existing empirical literature imposes cross-sectional homogeneity or homogeneity conditional on an exogenous grouping churchill2018environmental,bennedsen2023neural. We relax this condition. Using our FUSE-TIME procedure, different economies only share a common trend if the consistent grouping mechanism identifies homogeneous coefficients, enabling us to exploit the cross-sectional dimension without imposing restrictive homogeneity assumptions. In addition, given the small number of group-specific trend functions relative to the cross-sectional dimension, it becomes feasible to interpret the long-term behavior of CO\textsubscript{2} intensity for a large number of economies.
To estimate trends in the CO\textsubscript{2} intensity of GDP, we employ production-based CO\textsubscript{2} emissions data from the Global Carbon Project Friedlingstein2024global.\footnote{Dataset National fossil carbon emissions v2024, available at \href{https://globalcarbonbudgetdata.org}{globalcarbonbudgetdata.org}. Accessed November 14}. GDP series are obtained from the World Bank Development Indicators database.\footnote{Series NY.GDP.MKTP.CD, available at \href{https://data.worldbank.org/indicator/NY.GDP.MKTP.CD}{data.worldbank.org/indicator/NY.GDP.MKTP.CD}. Accessed November 12, 2024.} Both datasets are in annual frequency and span the time period from 1960 to 2023 for 92 countries. Details on the countries included in the sample and descriptive statistics are provided in Appendix (ref). Using this data, we compute the CO\textsubscript{2} intensity and construct an unbalanced panel where the individual series range from 29 to 64 years in length.
Let $y_{it} = \text{CO}\textsubscript{2}_{it} / \text{GDP}_{it}$, with CO\textsubscript{2} measured in million tonnes and GDP in billion 2024 U.S. dollar. We formulate the time-varying panel data model subject to a latent grouping
where $\beta_i (t/T)$ denotes the trend function of interest (time-varying intercept).
Searching over a dense grid of $d$, $M^*$, and $\lambda$ values, $d = 2$, $M^* = 4$, and $\lambda = 0.72$ yields the lowest IC, implying five latent groups; the entire grid is provided in Appendix (ref). Figure (ref) presents the estimated group-specific trend functions $\hat{\alpha}^p_{\hat{G}_k} (t/T)$. Figure (ref) sketches the spatial pattern in the group structure. Table (ref) in Appendix (ref) provides a detailed record of which countries are assigned to which group.
Group 1 predominantly comprises middle- and high-income economies in Eastern Europe and Asia, which have exhibited a substantial and ongoing decline in CO\textsubscript{2} emission intensity since the 1990s. Structural changes and rapid technological advances following the fall of the USSR and the economic liberalization of China may drive this trend. Group 2 largely includes low-to-middle-income economies, with outliers such as Switzerland, Hong Kong, and Singapore. The trend of Group 2 features a much more gradual decrease compared to the first group due to more moderate growth in both income and emissions. Group 3 primarily consists of high-income countries that experienced a sharp decline in emission intensity pre-1980, followed by a continued but more dampened decline thereafter. This pattern reflects a stark reduction in emissions until the 1980s, accompanied by consistent GDP growth throughout the entire sample period, aligning with a transition to service-based economies and the offshoring of emission- and energy-intensive sectors to low-income countries. Group 4 contains low-income and developing economies, characterized by the absence of a trend. In these economies, technological progress has yet to reduce emission intensity significantly or, as in the case of Brazil, income growth is largely driven by the exploitation of natural resources. Group 5, consisting of low-income economies, exhibits a similar but slightly more attenuated trend compared to Group 2.
This article introduces a novel penalized sieve estimation technique--FUSE-TIME--to estimate panel data models with smoothly time-varying coefficients that are subject to an unobserved group structure. Our approach simultaneously identifies the functional coefficients, the number of latent groups, and group compositions. In addition, we propose a consistent BIC-type criterion to determine the penalty tuning parameter, introduce a post-Lasso estimator, and prove oracle-efficiency. Monte Carlo simulation studies confirm strong finite-sample performance in recovering the latent group structure and accurately estimating the time-varying coefficients across a wide range of scenarios. We apply our method to analyze trends in the CO\textsubscript{2} emission intensity of GDP.
An important yet largely unexplored extension regards inference methods for (time-varying) panel data models with a latent group pattern. The estimated grouping is an inherently noisy representation of the true unobserved structure. Subsequently, the clustering introduces uncertainty to the slope coefficients, analogous to the well-established problem of non-uniform convergence of post-selection estimation greenshtein2004persistence,leeb2008sparse,adamek2023lasso. To the best of our knowledge, dzemski2024confidence provide the only approach to creating confidence sets of the grouping by inverting poolability tests to date. Consequently, the development of inference methods for penalized panel data models that admit a latent grouping, including confidence sets for the group structure, as well as tests for coefficient poolability and constancy akin to friedrich2024sieve, remains an important direction for future research.
\cleardoublepage \setstretch{0.1}
\doublespacing