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.
104,709 characters · 23 sections · 39 citation commands
Double Descent and Benign Overfitting in Macroeconomic Forecasting
\onehalfspacing
A striking finding in modern machine learning theory is that models that fit the training data perfectly can nonetheless forecast well out of sample. As the number of parameters grows, the forecast error traces a characteristic double-descent curve, first spiking at the interpolation threshold $N=T$ and then falling again, sometimes below the level attained by the best classical estimator. This phenomenon, known as benign overfitting, defies the classical bias-variance tradeoff intuition: more parameters than observations does not necessarily imply worse out-of-sample performance belkin2019reconciling. The conditions under which this occurs are now well understood in the linear regression case, as established by bartlett2020benign (henceforth, BLLT). The key requirement is a specific factor structure on the predictors: a small number of strong principal components must capture most of the signal, while the remaining components each explain a similar, negligible share of the variance. When this residual spectrum is sufficiently flat, the noise absorbed by the interpolating estimator is spread thinly and evenly across many components, and its contribution to forecast error becomes negligible.
Macroeconomic datasets provide a natural setting to test these conditions. A large body of work in empirical macroeconomics has established that panels of macroeconomic and financial indicators are well approximated by factor models, in which a small number of latent common factors drive the co-movement of many observable series stock2002forecasting, stock2002macroeconomic, bai2002determining. The implied covariance structure, a few spiked eigenvalues from the factors plus a residual idiosyncratic component, bears a strong resemblance to the spectral conditions that the benign overfitting literature requires. Whether these conditions hold in practice depends on the dimensionality of the panel: the noise dilution argument above requires the number of small principal components, roughly $N$ (the number of predictors in the panel) minus the number of factors, to be large relative to the sample size $T$. However, existing macroeconomic panels are only moderately high-dimensional, with typical cross-sectional dimensions of $N \sim 100$--$200$ and time series lengths of $T \sim 200$--$700$. Whether these conditions are met in practice, and whether the moderate dimensions of macroeconomic panels allow them to be exploited for forecast gains, are the central questions this paper addresses.
Our paper makes the following contributions. First, we provide the first illustration of the double-descent risk curve in macroeconomic forecasting, focusing on four representative macroeconomic series: CPI inflation, real personal income, industrial production, and the unemployment rate. In particular, we show how, in each case, the mean squared error (MSE) changes as the effective number of parameters in the forecasting model increases, starting from a simple AR(1) model and adding an increasing number of principal components until reaching the “overparameterized” regime, pushing the effective parameter count beyond the saturation/interpolation point where the number of parameters equals the sample size. A characteristic double-descent pattern emerges robustly across these four series: the MSE rises sharply near the interpolation point, then falls and stabilizes in the overparameterized regime.
Second, having documented that the double-descent pattern arises in macroeconomic data, we move to study when the factor structure of macroeconomic panels formally satisfies the conditions for benign overfitting. The theoretical results of bartlett2020benign were derived under i.i.d. data, and regardless of the dimension of the panel, it is not obvious that they carry over to the serially dependent, factor-structured data typical of macroeconomic panels. We show that they do, under both the exact and approximate factor model assumptions commonly maintained in the macroeconomic literature. Under the exact factor model, where idiosyncratic errors are cross-sectionally uncorrelated, the first two BLLT conditions hold automatically as consequences of the factor model assumptions, and the third condition holds whenever the idiosyncratic variances are sufficiently similar across series. Under the approximate factor model, which permits residual cross-sectional correlations, the first two conditions still hold, but the third becomes harder to satisfy: correlations among the idiosyncratic components can cause their variances to become more dispersed across series, making benign overfitting less likely.
Third, our first two contributions raise a practical question: if the BLLT conditions are satisfied in macroeconomic panels, can we actually exploit them given that these panels have only moderate dimensions ($N \sim 100$--$200$, $T \sim 200$--$700$), so that the BLLT regime $N/T \to \infty$ does not hold directly? We address this by proposing to augment the predictors with $B$ synthetic copies of the data generated using a factor model trained on the original time series, artificially increasing the effective $N/T$ ratio. The key insight is that each synthetic copy preserves the estimated factor structure while drawing fresh idiosyncratic noise, so that two time periods with similar factor values produce similar synthetic predictors; in the limit, this implicit similarity measure is precisely a kernel. We prove that, as the number of synthetic replications $B \to \infty$, this procedure converges to kernel ridge regression with a factor-structured kernel that measures similarity between time periods through their estimated factor realizations, and a ridge parameter controlled by the total idiosyncratic variance. Thus, the synthetic expansion works not by adding information, but by constructing an implicit kernel whose structure is inherited from the factor model.
These contributions together support a central conclusion: benign overfitting, when it works, succeeds because overparameterization implicitly constructs a well-behaved kernel, not because overparameterization is intrinsically desirable.
We then turn to verify this claim empirically, evaluating how a semiparametric factor kernel model (obtained by combining a linear AR specification with a nonparametric kernel component) fares against the standard factor model benchmark stock2002forecasting, stock2002macroeconomic. We perform this comparison using all available series in the FRED-MD dataset of mccracken2016fred, forecasting $h = 1, 3, 6, 12$ months ahead using a rolling-window setup where the AR(1) dynamics is estimated via generalized least squares (GLS). We find that the factor kernel consistently outperforms the factor model benchmark across all series and forecast horizons. The gains are pervasive --- the factor kernel beats the benchmark for 92--99% of monthly targets and 96--99% of quarterly targets --- and become increasingly large and statistically significant as the forecast horizon lengthens. At short horizons, improvements are modest but highly consistent; at longer horizons, more than half of the wins are statistically significant by the Diebold--Mariano test. We interpret these findings as empirical evidence that the semiparametric factor kernel provides regularization naturally adapted to the data's factor structure, connecting the modern interpolation literature to the classical tradition of factor-based forecasting stock2002forecasting, stock2002macroeconomic and to kernel methods in economics.
This paper contributes to multiple strands of literature. First, we add to the recent empirical literature that has documented how ML methods can outperform standard benchmarks for macroeconomic forecasting medeiros2021forecasting, coulombe2022machine, chi2025macroeconomic. We complement the existing work by identifying the spectral mechanism that determines when and why such gains can manifest. Our kernel-based approach is related to the nonlinear and kernel methods applied to factor-augmented forecasting by exterkate2016nonlinear and boot2019macroeconomic, though our kernel arises endogenously from the factor structure rather than being imposed exogenously. We also add to the literature on the theory behind benign overfitting. bartlett2020benign establish sufficient and necessary conditions under i.i.d.\ data; nakakita2022benign extend the analysis to dependent processes; hastie2022surprises and bunea2021interpolating provide further characterizations in related settings. We contribute to this theory by mapping the BLLT conditions onto the factor structure of macroeconomic panels, characterizing when they hold and when they fail, and deriving an implementable kernel estimator that achieves benign overfitting when the native $N/T$ ratio is insufficient. The factor model forecasting tradition of stock2002forecasting, stock2002macroeconomic provides both our empirical benchmark and the structural motivation for the factor kernel. Our work builds on the static factor model framework; the dynamic factor model literature forni2000generalized, forni2005generalized and the targeted predictor approach of bai2008forecasting offer complementary perspectives on dimension reduction in large panels. Finally, our data augmentation strategy is related in spirit to the Bayesian shrinkage approach of demol2008forecasting, which also exploits the factor structure for regularization, though through a different mechanism.
The remainder of the paper is organized as follows. Section (ref) introduces the ridgeless estimation framework and characterizes the double-descent phenomenon. Section (ref) reviews the BLLT conditions for benign overfitting and their extension to dependent data. Section (ref) establishes when factor-structured data satisfy the BLLT conditions, showing that the first two conditions hold automatically while the third depends on the dispersion of idiosyncratic variances. Section (ref) develops the synthetic factor augmentation strategy and proves its equivalence to kernel ridge regression with a factor-structured kernel. Section (ref) presents the out-of-sample forecasting evaluation comparing the factor kernel to the Stock-Watson benchmark. Section (ref) concludes. The appendices contain proofs (Appendix (ref)), the double-descent analysis using quarterly FRED-QD data (Appendix (ref)), an illustration of bias and variance behavior under alternative spectral shapes (Appendix (ref)), a derivation of the relationship between tail index and concentration ratio (Appendix (ref)), and the GLS formulation of the semiparametric kernel regression (Appendix (ref)).
Macroeconomic forecasters routinely face many more candidate predictors than observations. Standard datasets such as FRED-MD contain over 120 monthly series against sample sizes of $T \sim 400$--$800$, and including all series directly in a forecasting regression inflates estimation variance and degrades out-of-sample performance. The literature's response has been to reduce the effective parameter count: extract a small number of common factors stock2002forecasting, stock2002macroeconomic, shrink coefficients via ridge regression or the LASSO medeiros2021forecasting, coulombe2022machine, or select variables via information criteria or Bayesian model averaging. All of these approaches keep the number of estimated parameters well below $T$, appealing to the familiar bias-variance tradeoff to justify parsimony.
Inspired by recent developments in the machine learning literature, this paper investigates a different approach: rather than shrinking toward more parsimonious models, we ask what happens when the number of regressors $N$ is pushed past the interpolation threshold $N = T$ (the point at which the model can fit every observation exactly) and into the overparameterized regime. We formalize this setup below and characterize the conditions under which this approach can deliver good out-of-sample forecasts.
Let $\mathbf{x}_t \in \mathbb{R}^N$ be a vector of macroeconomic predictors at date $t$ and let $y_t$ be a scalar target series. We assume that the forecasting model is linear:
where $\boldsymbol{\beta}^* \in \mathbb{R}^N$ is the vector of population projection coefficients and $\varepsilon_t$ is the forecast error, with $\mathbb{E}[\varepsilon_t \mid \mathbf{x}_t] = 0$ and $\mathbb{E}[\varepsilon_t^2] = \sigma^2$. In matrix notation,
where $\mathbf{y} \in \mathbb{R}^{T}$ and $\mathbf{X} \in \mathbb{R}^{T \times N}$ stacks the $T$ observations row-wise.
The behavior of the least-squares estimator depends critically on the relationship between $N$ and $T$. When $N < T$ (under-parameterized regime), $\mathbf{X}$ has full column rank, $\mathbf{X}'\mathbf{X}$ is invertible, and OLS delivers the unique unbiased estimator $\hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}$. The model is identified in the usual sense and the standard bias-variance tradeoff applies. When $N = T$ (interpolation point), $\mathbf{X}$ is square and, if full rank, the system $\mathbf{X}\boldsymbol{\beta} = \mathbf{y}$ has a unique solution that fits the training data with zero residuals, $\hat{\boldsymbol{\beta}} = \mathbf{X}^{-1} \mathbf{y}$. In contrast, when $N > T$ (overparameterized regime), the system $\mathbf{X}\boldsymbol{\beta} = \mathbf{y}$ is underdetermined: in other words, there are infinitely many coefficient vectors $\hat{\boldsymbol{\beta}}$ that produce a perfect in-sample fit, and since $\mathbf{X}'\mathbf{X}$ is rank-deficient standard OLS cannot be applied directly.
Among the infinitely many solutions that perfectly fit the training data when $N > T$, a natural choice is the one that solves the normal equations with the smallest coefficient norm, i.e. the minimum norm estimator:
which has closed form solution
where $(\mathbf{X}'\mathbf{X})^{\dagger}$ is the pseudoinverse of $(\mathbf{X}'\mathbf{X})$. This estimator turns out to correspond, in the limit, to the ridge estimator $\hat{\boldsymbol{\beta}}_\lambda = (\mathbf{X}'\mathbf{X} + \lambda \mathbf{I}_N)^{-1}\mathbf{X}'\mathbf{y}$ as the penalty parameter $\lambda$ approaches zero. For this reason, it is also referred to as the ridgeless estimator.
Depending on whether $\mathbf{X}$ is full row or full column rank, we end up with the following three cases:
The underparameterized ($N<T$) and overparameterized ($N>T$) regimes imply a qualitatively different behavior for the mean squared forecast error (MSFE), and the transition from one regime to the next can give rise to the double-descent pattern we mentioned in the Introduction. We will characterize this formally in Section (ref). Before turning to that, we first examine what happens as we transition across the interpolation threshold ($N = T$) in terms of the identifiability of $\boldsymbol{\beta}^*$.
As it is well known, when $N < T$ the model is identified in the usual sense: $\mathbf{X}$ has full column rank $N$, $\mathbf{X}'\mathbf{X}$ is invertible and OLS uniquely determines $\hat{\boldsymbol{\beta}}$, with $\mathbb{E}[\hat{\boldsymbol{\beta}} \mid \mathbf{X}] = \boldsymbol{\beta}^*$. The estimator is unbiased and $\boldsymbol{\beta}^*$ is point-identified. The interpolation threshold marks a fundamental change in identification. When $N>T$, the data matrix $\boldsymbol{X}$ has full row rank $T$, it admits only a right inverse and can no longer uniquely identify $\boldsymbol{\beta}^*$. The minimum norm estimator will be biased:
where $\bm{P}_r = \boldsymbol{X}'(\boldsymbol{X}\boldsymbol{X}')^{-1}\boldsymbol{X}$ is the orthogonal projector onto the row space of $\boldsymbol{X} \subset \mathbb{R}^N$. Intuitively, with $T$ observations and $N > T$ predictors, the training data can pin down at most $T$ linearly independent combinations of the $N$ parameters. The remaining $N - T$ combinations fall in the null space of $\mathbf{X}$ and are invisible to any estimator. $\boldsymbol{\beta}^*$ is therefore set-identified and the identified set is the affine subspace
\[ \mathcal{I}(\boldsymbol{\beta}^*) = \{\tilde{\boldsymbol{\beta}} \in \mathbb{R}^N : \mathbf{X}\tilde{\boldsymbol{\beta}} = \mathbf{X}\boldsymbol{\beta}^*\} = \boldsymbol{\beta}^* + \mathrm{null}(\mathbf{X}), \] The data can only help identify the component of $\boldsymbol{\beta}^*$ in the row space of $\mathbf{X}$, which spans the $T$ directions that can be inferred from the estimation sample. The remaining $N - T$ directions, lying in the null space of $\mathbf{X}$, are unrecoverable regardless of the available sample size. We can then decompose $\boldsymbol{\beta}^*$ as\footnote{The null-space membership $\boldsymbol{\beta}^*_\perp\in\mathrm{null}(\mathbf{X})$ follows from $\mathbf{X}(\mathbf{I}-\mathbf{P}_r) = \mathbf{X} - \mathbf{X}\mathbf{X}'(\mathbf{X}\mathbf{X}')^{-1}\mathbf{X} = \mathbf{X} - \mathbf{X} = \mathbf{0}$.} \[ \boldsymbol{\beta}^* = \mathbf{P}_r\boldsymbol{\beta}^* + \boldsymbol{\beta}^*_{\perp}, \qquad \boldsymbol{\beta}^*_{\perp} := (\mathbf{I} - \mathbf{P}_r)\boldsymbol{\beta}^* \in \mathrm{null}(\mathbf{X}), \] The minimum-norm estimator recovers $\mathbf{P}_r\boldsymbol{\beta}^*$ and sets the null-space component to zero,\footnote{Substitute $\mathbf{y}=\mathbf{X}\boldsymbol{\beta}^*+\boldsymbol{\varepsilon}$ into $\hat{\boldsymbol{\beta}}=\mathbf{X}'(\mathbf{X}\mathbf{X}')^{-1}\mathbf{y}$ to get $\hat{\boldsymbol{\beta}} = \mathbf{P}_r\boldsymbol{\beta}^* + \mathbf{X}'(\mathbf{X}\mathbf{X}')^{-1}\boldsymbol{\varepsilon}$; taking $\mathbb{E}[\,\cdot\mid\mathbf{X}]$ yields $\mathbb{E}[\hat{\boldsymbol{\beta}}\mid\mathbf{X}]=\mathbf{P}_r\boldsymbol{\beta}^*$.} so its bias is precisely the unidentified component: \[ \mathbb{E}[\hat{\boldsymbol{\beta}} \mid \mathbf{X}] - \boldsymbol{\beta}^* = -\boldsymbol{\beta}^*_{\perp}. \]
Since every element of $\mathcal{I}(\boldsymbol{\beta}^*)$ produces the same in-sample fit and the same forecasts for test points within the training span, any selection rule that picks one element over another is essentially an identifying restriction on the unidentified component. The ridgeless estimator sets all unidentifiable coefficient directions to zero: \[ (\mathbf{I} - \mathbf{P}_r)\hat{\boldsymbol{\beta}} = \mathbf{0}. \] The choice of the selection rule has no consequence for in-sample fit, since the fitted values are an unbiased estimator of the conditional mean of $Y$: $\mathbb{E}[\boldsymbol{X}\boldsymbol{\hat{\beta}}\mid \boldsymbol{X}]=\boldsymbol{X}\bm{P}_r\boldsymbol{\beta}=\boldsymbol{X}\boldsymbol{\beta}$; It matters only through its effect on the bias component $\boldsymbol{\beta}^*_\perp$ when forecasting future observations whose predictor values fall outside the span of the training data.
Suppose we forecast a new observation $(\mathbf{x}_0, y_0)$ with $\mathbf{x}_0 \sim (0, \boldsymbol{\Sigma})$ and $y_0 = \mathbf{x}_0'\boldsymbol{\beta}^* + \varepsilon_0$. The MSFE expands as
and decomposes into bias and variance components:\footnote{Write $\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}^* = (\hat{\boldsymbol{\beta}} - \mathbb{E}[\hat{\boldsymbol{\beta}} \mid \mathbf{X}]) + (\mathbb{E}[\hat{\boldsymbol{\beta}} \mid \mathbf{X}] - \boldsymbol{\beta}^*)$ and expand the quadratic form in (ref). Taking $\mathbb{E}[\,\cdot \mid \mathbf{X}]$, the cross-term vanishes because $\mathbb{E}[\hat{\boldsymbol{\beta}} - \mathbb{E}[\hat{\boldsymbol{\beta}} \mid \mathbf{X}] \mid \mathbf{X}] = \mathbf{0}$. The first remaining term is $\mathbb{E}_\mathbf{X}[\operatorname{tr}(\boldsymbol{\Sigma}\,\text{Var}(\hat{\boldsymbol{\beta}} \mid \mathbf{X}))]$ (the variance component) and the second is $\mathbb{E}_\mathbf{X}[\|\boldsymbol{\Sigma}^{1/2}(\mathbb{E}[\hat{\boldsymbol{\beta}} \mid \mathbf{X}] - \boldsymbol{\beta}^*)\|^2]$ (the squared bias).}
Importantly, the $\boldsymbol{\Sigma}$-weighting reflects the fact that forecast accuracy depends not on parameter error per se but on how that error maps into prediction error through the predictor covariance. This is the key to understanding how double descent works. The two components take qualitatively different forms on either side of the interpolation threshold. In the underparameterized regime ($N < T$), OLS is unbiased (Bias$^2 = 0$) and \[ \mathrm{Variance} = \mathbb{E}[\sigma^2\,\mathrm{tr}\!\big(\boldsymbol{\Sigma}(\mathbf{X}'\mathbf{X})^{-1}\big)], \] In the overparameterized regime ($N > T$), using $\mathbb{E}[\hat{\boldsymbol{\beta}} \mid \mathbf{X}] = \mathbf{P}_r\boldsymbol{\beta}^*$,
where $\boldsymbol{\beta}^*_\perp = (\mathbf{I} - \mathbf{P}_r)\boldsymbol{\beta}^*$ is the unrecovered component defined in the preceding subsection. The derivation of these variance expressions is given in Appendix (ref).
The spike in MSFE at the interpolation threshold is a guaranteed consequence of the variance divergence as $N \to T$. Seeing this is easy under the (standard) assumption that the rows of $\bm X$ follow a Gaussian distribution with zero mean and variance $\bm \Sigma$, in which case we can write the variance term in the $N < T$ case as: \[ \mathbb{E}[\sigma^2\,\mathrm{tr}\!\big(\boldsymbol{\Sigma}(\mathbf{X}'\mathbf{X})^{-1}\big)] = \sigma^2 \frac{\mathrm{tr}(\bm \Sigma \bm \Sigma^{-1})}{T-N-1} = \sigma^2 \frac{N}{T-N-1} \] which diverges as $N \to T-1$. A similar pattern occurs in the $N \to T^+$ case. To see this, again assuming Gaussianity for simplicity, write $\bm X = \bm Z \bm\Sigma^{1/2}$ where $\bm Z$ is $T \times N$ with i.i.d.\ $\mathcal{N}(0,1)$ entries. Since $\bm X\bm X^\top = \bm Z\bm\Sigma \bm Z^\top \preceq \lambda_{\max}(\bm\Sigma)\,\bm Z\bm Z^\top$, we have $(\bm X\bm X^\top)^{-1} \succeq \lambda_{\max}(\bm\Sigma)^{-1}(\bm Z\bm Z^\top)^{-1}$, where $\lambda_{s} (\bm \Sigma)$ (for $s \in \{\text{min}, \text{max} \}$) denotes the minimum/maximum eigenvalue of $\bm \Sigma$.\footnote{We stick to this notation throughout the remainder of the paper} Since $\bm Z\bm Z^\top \sim W_T(N, \bm I_T)$, taking traces and expectations gives the lower bound: \[ \mathbb{E}\!\big[\mathrm{tr}\{(\bm X\bm X^\top)^{-1}\}\big] \;\geq\; \frac{1}{\lambda_{\max}(\bm\Sigma)}\,\frac{T}{N-T-1} \] for $N > T+1$, which diverges as $N \to T+1$. The bound discards all spectral structure of $\bm\Sigma$ except its largest eigenvalue, but suffices to establish that the variance diverges at the interpolation threshold for any $\bm\Sigma$ with bounded $\lambda_{\max}$. At $N = T$, both the underparameterized formula $\sigma^2 \frac{N}{T-N-1}$ and the overparameterized bound $\frac{1}{\lambda_{\max}(\bm\Sigma)}\,\frac{T}{N-T-1}$ are undefined, since both denominators are negative. This is the interpolation threshold at which the double descent peak occurs.
If bias and variance are both bounded as $N/T$ grows beyond the interpolation point, the MSFE will exhibit a double descent, that is a pattern in which it reduces twice: a first descent in the underparameterized regime, a spike at the interpolation threshold, and a second descent as the model moves deeper into the overparameterized regime.
Whether the second descent would eventually produce a lower MSFE than the underparameterized regime, i.e. benign overfitting, is not guaranteed. To achieve benign overfitting it is not sufficient for the bias and variance terms to be bounded, they have to approach zero, which means the MSFE will have to achieve the noise floor $\sigma^{2}$. As we shall see, this can happen when $\boldsymbol{\Sigma}$ has a spectral shape characterized by a few large eigenvalues and many small, but non-vanishingly small, eigenvalues. In such a situation as $N$ grows, the row space of $\mathbf{X}$ progressively aligns with the high-variance directions, pushing the unrecovered component $\boldsymbol{\beta}^*_{\perp}$ into low-variance tail directions where it contributes little to forecast error. Simultaneously, the training noise gets absorbed thinly across many regressors, reducing the variance term.
Section (ref) established that the MSFE spike at the interpolation threshold is guaranteed, and that a second descent follows whenever the bias term (ref) and variance term (ref) remain bounded as $N$ grows beyond $T$. We now ask whether this pattern is visible in macroeconomic data. Using the FRED-MD dataset, which contains $N = 122$ monthly series spanning 1959:01 to 2025:09 ($T = 798$ observations), we trace the full MSFE curve from the underparameterized region, through the interpolation threshold, and into the overparameterized regime. We split the sample in half, with $T_{\mathrm{train}} = 399$ and $T_{\mathrm{test}} = 399$.\footnote{In this motivating exercise the single-split design is meant to make the double-descent MSFE curve visible. Whether it also delivers systematic forecasting gains is addressed in Section (ref), where we use a full rolling-window out-of-sample design across all targets and horizons. A parallel analysis using the quarterly FRED-QD dataset is presented in Appendix (ref).} All models include an AR(1) lag as a baseline predictor.
In the notation of equation (ref), $\bm{y}$ is the $T_{\mathrm{train}} \times 1$ vector of the target variable and $\mathbf{X}$ is the $T_{\mathrm{train}} \times N$ matrix stacking the remaining $N = 122$ FRED-MD series as predictors, with each variable standardized using training-window moments and transformed using the FRED-MD provided transformation codes. We consider four target series: CPI inflation (CPIAUCSL), real personal income (RPI), industrial production (INDPRO), and the unemployment rate (UNRATE). As with the predictors, all four target series were transformed using the FRED-MD transformation codes.
We begin by taking an AR(1) model (one parameter), and gradually increasing model size toward $N = 122$ by adding $r$ principal components (PCs) of $\mathbf{X}$, ranked by eigenvalue, one at a time, where $r$ runs from $1$ to $N$. This yields an ordered path from the AR(1)+1 PC model to the AR(1)+all $N$ PCs model (which is algebraically equivalent to using all $N$ original predictors). The first set of PCs captures the dominant common factors, while later PCs increasingly capture idiosyncratic dynamics.
At $r = N$ the model uses all available variables and we have reached the boundary of the original data. To push into the overparameterized regime (where the number of effective parameters exceeds the number of training observations), we augment the design matrix with synthetic copies of $\bm X$. Each copy preserves the estimated common component (extracted using $\hat{k}$ Bai--Ng factors) but replaces the idiosyncratic errors with fresh Gaussian draws. Appending $B$ such copies creates a design matrix with $N(B+1)$ columns, so the effective number of parameters grows in discrete jumps of $N$. The construction and its theoretical properties are developed in Section (ref); here we simply use it as a device to extend the MSFE curve beyond the interpolation threshold.
Throughout, $N_{\mathrm{eff}}$ denotes the total number of slope coefficients being estimated (cross-sectional predictors plus the AR(1) lag): $N_{\mathrm{eff}}=r+1$ along the PC path and $N_{\mathrm{eff}}=(B+1)N+1$ with augmentation. The interpolation threshold is reached when $N_{\mathrm{eff}}=T_{\mathrm{train}}$. Since we only have $N=122$ PCs, we start augmenting the dataset with artificial covariates after we have included all PCs. After crossing the interpolation point this augmentation path then extends the analysis into the deeply overparameterized regime.
We apply this procedure to each of the four targets and plot the resulting MSFE against the effective number of parameters. Figure (ref) shows a clear and consistent double-descent pattern across all four series: MSFE rises sharply near the interpolation threshold and then declines in the overparameterized region, confirming that the second descent is a robust feature of macroeconomic data.
The natural next question is whether this decline is deep enough to constitute benign overfitting: in a true out of sample setting, does the MSFE on the right branch eventually approach the noise floor $\sigma^2$, or does it stabilize at a higher level? A rigorous assessment is carried out in Section (ref). Before turning to that, we characterize the theoretical conditions for benign overfitting and show how the factor structure of macro panels naturally satisfies them.
This section characterizes the conditions under which benign overfitting occurs. The material draws on the papers cited and is included for completeness, providing necessary background for Section (ref). Readers already familiar with the benign overfitting literature may skip directly to Section (ref).
We start with analyzing separately the bias and the variance term in expression ((ref)), and we discuss heuristically how and in which circumstances they may converge to zero. Then we formally introduce the bartlett2020benign (BLLT) conditions for i.i.d. data. Finally, we discuss how the BLLT conditions apply also to dependent data, drawing on the results of nakakita2022benign.
Before proceeding, we introduce terminology that will be used throughout. The spectrum of the $N \times N$ predictor covariance matrix $\boldsymbol{\Sigma}$ refers to its ordered eigenvalues $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_N > 0$. We partition the spectrum into two regions: the head consists of the first $k$ eigenvalues, which are large and capture the dominant directions of variation (in factor models, these correspond to the common factors); the tail consists of the remaining $N - k$ eigenvalues, which are small and represent residual variation (in factor models, these correspond to the idiosyncratic components). The split point $k$ is chosen to separate eigenvalues of order $O(N)$ from those of order $O(1)$. The key to benign overfitting is that the head must be low-dimensional (small $k$) while the tail must be both flat (eigenvalues of similar magnitude) and large-dimensional (large $N - k$), so that noise can be dispersed across many weak directions.
As we discussed above, double descent merely requires that the MSFE eventually decreases past the interpolation point. It says nothing about the level to which it decreases. Benign overfitting is a strictly stronger statement: it requires that the MSFE converge to the noise floor $\sigma^2$ as $N/T \to \infty$, meaning that both the bias and the variance components must converge to zero.
\paragraph{Bias term.} As $N/T \to \infty$ the unrecovered component $\boldsymbol{\beta}^*_{\perp}$ does not vanish in Euclidean norm, which means the minimum norm estimator remains biased. However, the crucial observation is that the bias term in the MSFE rescales the bias by $\boldsymbol{\Sigma}$, i.e., each direction of $\boldsymbol{\beta}^*_{\perp}$ is weighted by the corresponding predictor variance. This means that even though the bias does not vanish in Euclidean norm, the $\boldsymbol{\Sigma}$-weighted bias might, depending on the eigenvalue structure of $\boldsymbol{\Sigma}$ and the alignment of $\boldsymbol{\beta}^*$ with the factor space.
The key mechanism is as follows. Suppose $\boldsymbol{\Sigma}$ features a strong head of large eigenvalues and a flat tail of small, non-vanishing eigenvalues, and suppose the signal $\boldsymbol{\beta}^*$ is concentrated in the head directions. As $N$ increases, the projection operator $\bm{P}_r$ evolves because $\mathbf{X}$ gains new columns. The row space becomes progressively better aligned with the head space where $\boldsymbol{\beta}^*$ has support, and the residual $(I - \bm{P}_r)\boldsymbol{\beta}^*$ shrinks geometrically. Moreover, since the remaining residual lies in low-variance tail directions with tiny eigenvalues $\lambda_j$, the $\boldsymbol{\Sigma}$-weighted bias $\|\boldsymbol{\Sigma}^{1/2}(\bm I - \bm{P}_r)\boldsymbol{\beta}^*\|^2$ decreases even faster than the unweighted residual.
This mechanism, however, requires that the tail eigenvalues do not vanish too quickly. If they decay exponentially fast, adding new predictors contributes negligible variance to the row space, and $\bm{P}_r$ stops evolving before fully aligning with the head space. In this case, a fixed portion of $\boldsymbol{\beta}^*$ remains unrecovered regardless of $N$, and the bias does not vanish. Conversely, if the tail is too heavy (eigenvalues decay too slowly), the unrecovered component $\boldsymbol{\beta}^*_{\perp}$ may lie in directions with non-negligible variance, preventing the weighted bias from vanishing.
Thus, benign overfitting requires the tail to be flat: disperse enough (many directions) to allow the row space to expand and cover the signal support, yet light enough (small eigenvalues) to ensure the weighted bias vanishes. Appendix (ref) illustrates these requirements through three concrete spectral configurations. In a factor-structured panel, this structure arises naturally. The high-variance head directions correspond to the common factors that drive cross-sectional covariation, while the flat tail corresponds to the idiosyncratic components with similar small variances. If the variable we are forecasting loads predominantly on the factors, then $\boldsymbol{\beta}^*$ is concentrated in the head, $\boldsymbol{\beta}^*_{\perp}$ ends up in the idiosyncratic tail directions, and its contribution to MSFE is negligible even if $\|\boldsymbol{\beta}^*_{\perp}\|$ is not.\footnote{This requires the predictive signal $\boldsymbol{\beta}^*$ to have significant support in fewer than $T$ directions. If it is spread across more than $T$ directions, the $T$-dimensional row space of $\mathbf{X}$ cannot cover them all, and a fixed portion of the signal remains unrecovered regardless of $N$.}
\paragraph{Variance term.} The variance component measures how much in-sample noise is amplified into out-of-sample forecast error by the minimum-norm interpolator. When $N > T$, the interpolator fits the training noise exactly but distributes it across all $N$ coefficient directions. Here, again, the eigenvalues of $\boldsymbol{\Sigma}$ play a key role.
If $\boldsymbol{\Sigma}$ has a flat tail of small but non-vanishing eigenvalues, then as $N$ grows, the minimum-norm solution distributes the fit over more columns, keeping the $\ell_2$ norm of the coefficient vector small. Since each tail direction contributes little variance, the noise gets absorbed thinly across many regressors, and the out-of-sample variance decreases. Conversely, if the tail is concentrated on a few large eigenvalues, noise accumulates in high-variance directions and the variance remains large. In a factor-structured panel, the idiosyncratic components provide precisely this flat spectral tail: many series with similar small variances across which the interpolated noise is spread evenly.
The BLLT conditions formalize exactly what is the spectral shape of $\bm{\Sigma}$ for which benign overfitting may happen: conditions (i) and (ii) ensure the spectrum has a well-defined factor/idiosyncratic head-tail structure with the signal confined to the head, while condition (iii) ensures the idiosyncratic tail is flat enough to disperse noise. We now state these conditions formally, following bartlett2020benign. Let $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_N > 0$ be the eigenvalues of $\bm{\Sigma}$. Define two measures of effective rank for the tail beyond coordinate $k$:
The first, $r_k$, measures the total tail variance relative to the largest tail eigenvalue; the second, $R_k$, measures how uniform the tail is. When all tail eigenvalues are equal, $r_k = R_k = N - k$. When the tail is concentrated on a few directions, $R_k \ll r_k$. Note that both $r_k$ and $R_k$ are scale-invariant: they depend only on the ratios among the tail eigenvalues. Consequently, the BLLT ratios below are unaffected by an overall rescaling of $\bm{\Sigma}$.
Define the split index \[ k^* := \min\big\{k \geq 0 : r_k(\bm{\Sigma}) \geq b\,T\big\}, \] for a universal constant $b > 0$. The split index identifies the point in the ordered spectrum where the tail becomes sufficiently spread out that noise can be dispersed across roughly $T$ weak directions. It separates the “head” (strong directions) from the “tail” (weak, noise-dispersing directions).
The three ratios capture:
These quantities and conditions follow bartlett2020benign (their Definition 3, Definition 4, and Theorem 1). Their determinant normalization is without loss of generality: since $r_k$, $R_k$, and $k^*$ depend only on eigenvalue ratios, rescaling $\bm{\Sigma}$ does not affect whether (ref)--(ref) hold.
The BLLT theory assumes that the rows of $\bm{X}$ are drawn independently from a sub-Gaussian distribution. In macroeconomic applications the rows of $\bm{X}$ are consecutive time-series observations and therefore serially dependent. nakakita2022benign extend the theory to dependent (sub-Gaussian) data and show that the BLLT spectral conditions still apply, with an adjustment for temporal dependence.
The key idea is a decorrelation argument: a multivariate whitening transformation removes the autocorrelation structure and restores the i.i.d.\ setting to which BLLT applies directly. The results then translate back to the original scale, and the BLLT bounds apply directly to the untransformed data, with an inflated variance term that accounts for the additional uncertainty from temporal dependence. The only additional requirement is a coherence condition (Assumption (ref)) ensuring that the temporal dependence structure is such that an appropriate decorrelation transformation exists.
Let $\{\mathbf{x}_t\}_{t=1}^T$ be a stationary $N$-dimensional process with autocovariance function \[ \bm{\Gamma}_X(h) := \mathbb{E}[\mathbf{x}_t \mathbf{x}_{t-h}'] \in \mathbb{R}^{N \times N}, \qquad h = 0, 1, 2, \ldots \] Let $\bm{\Sigma} = \bm{U}\bm{\Lambda}_\Sigma \bm{U}'$ be the spectral decomposition of $\bm{\Sigma} = \bm{\Gamma}_X(0)$, with eigenvectors $\bm{u}_1, \ldots, \bm{u}_N$ and eigenvalues $\lambda_1 \geq \cdots \geq \lambda_N > 0$. Each eigenvector $\bm{u}_i$ defines a linear combination of the $N$ original series; the index $i = 1, \ldots, N$ thus runs over these $N$ directions, which in a factor model correspond to the $k$ common factors (large eigenvalues) and the $N-k$ idiosyncratic components (small eigenvalues). The projected series $\tilde{\mathbf{x}}_{it} := \bm{u}_i' \mathbf{x}_t$ is the scalar time series along the $i$-th such direction, with autocovariance
For each $i$, let $\bm{\Xi}_{i,T}$ denote the $T \times T$ autocovariance matrix of the projected series $(\tilde{\mathbf{x}}_{i1}, \ldots, \tilde{\mathbf{x}}_{iT})'$. When $\gamma_i(h)$ does not depend on $i$, we can write $\bm{\Xi}_{i,T} = \bm{\Xi}_T$ for all $i$ and the covariance of $\mathrm{vec}(\bm{X})$ has the separable Kronecker form $\bm{\Sigma} \otimes \bm{\Xi}_T$. In general, $\gamma_i(h)$ varies across directions, and the autocovariance takes the non-separable form \[ \mathrm{Cov}\big(\mathrm{vec}(\bm{X})\big) = \sum_{i=1}^N \Xi_{i,T} \otimes (u_i u_i'). \]
\paragraph{Decorrelation.} nakakita2022benign model the design matrix through the factorization \[ \bm{X} = \big[\bm{\Xi}_{1,T}^{1/2}\,z_1,\;\ldots,\;\bm{\Xi}_{N,T}^{1/2}\,z_N\big]\,\bm{\Lambda}_\Sigma^{1/2}\,\bm{U}^\top, \] where $\bm{z}_i \sim \mathcal{N}(0, \bm{I}_T)$ are i.i.d.\ In the separable case $\bm{\Xi}_{i,T} = \bm{\Xi}_T$, this reduces to $\bm{X} = \bm{\Xi}_T^{1/2}\,\bm{Z}\,\bm{\Lambda}_\Sigma^{1/2}\,\bm{U}^\top$, and pre-multiplying by $\bm{\Xi}_T^{-1/2}$ yields $\tilde{\bm{X}} = Z\Lambda_\Sigma^{1/2}U^\top$, which has i.i.d.\ rows with covariance $\bm{\Sigma}$, exactly the BLLT setting. The minimum-norm interpolant of the original problem coincides with that of the whitened problem, so the BLLT spectral conditions apply directly after decorrelation. In the non-separable case, each projected series requires its own whitening transformation $\bm{\Xi}_{i,T}^{-1/2}$, and no single transformation decorrelates all directions simultaneously. The following assumption ensures these direction-specific structures are sufficiently similar that a single reference transformation approximately decorrelates all of them.
Coherence bounds how much $\gamma_i(h)$ (defined in (ref)) can vary across eigendirections: it rules out situations where some directions are near-integrated while others are white noise. This condition is satisfied by any stationary VAR with companion matrix $\bm A$ and spectral radius $\rho(\bm A)<1-\delta$ for some $\delta>0$ and eigenvalues bounded away from zero, conditions met by standard macroeconomic VARs after applying the McCracken--Ng stationarity transformations. More broadly, the condition holds for any stationary VARMA process with bounded spectral density.
Crucially, nakakita2022benign show that when coherence is satisfied, there is no need to actually perform any transformation: the usual BLLT bounds hold directly on the untransformed data, although with an inflated variance term on account of the temporal dependence.
This Section specializes the general setup of Section (ref) and Section (ref) to factor-structured predictors. We derive the key spectral implications of both exact and approximate factor structures and verify when the BLLT conditions hold.
As detailed in Section (ref), benign overfitting requires three conditions: (i) the overall effective rank must be small relative to the sample size; (ii) the signal must be concentrated in a low-dimensional subspace; (iii) the residual spectrum must be sufficiently flat to disperse noise. This Section shows that with factor-structured predictors, the first two conditions hold automatically as consequences of the factor structure, while the third condition, the tail condition, depends on the dispersion of the idiosyncratic variances. The tail flatness condition is easier to satisfy under exact factor models than under approximate factor models. In the exact case, the tail eigenvalues coincide with the idiosyncratic variances, so flatness requires only that no series has idiosyncratic variance orders of magnitude larger than others. In the approximate case, residual correlations among idiosyncratic components can concentrate variance into a few large tail eigenvalues, creating spikes that reduce the effective dimension of the tail and violate the flatness requirement. Readers not interested in the technical details can skip to Section (ref).
Our analysis is related to but distinct from bunea2021interpolating, who study benign overfitting in a factor regression model, where the response variable depends directly on the latent factors. By contrast, we study linear regression where the predictors have a factor structure but the response may depend on all $N$ predictors. This distinction matters for the bias condition: in a factor regression the signal is confined to the factor space by assumption, whereas in our setting the signal alignment with the factor space is an empirical question. Moreover, bunea2021interpolating is based on an exact factor model within an i.i.d. setting, while we allow for an approximate factor structure and temporal dependence, both of which are typical features of macroeconomic data.
Recall the linear regression model shown in equation (ref): \[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta}^* + \boldsymbol{\varepsilon}. \] Our goal is to estimate this model directly, without first extracting factors or reducing dimensionality. That is, the factor structure of the regressors is not used in the estimation step; rather, it provides the spectral configuration that enables benign overfitting. Specifically, we assume each row of $\mathbf{X}$ has a factor structure: $\mathbf{x}_t = \bm{\Lambda} \mathbf{f}_t + \mathbf{e}_t, t = 1, \ldots, T,$ where $\mathbf{f}_t \in \mathbb{R}^k$ is a vector of $k$ latent common factors, $\bm{\Lambda} \in \mathbb{R}^{N \times k}$ is the matrix of factor loadings, and $\mathbf{e}_t \in \mathbb{R}^N$ is the idiosyncratic error. The population covariance is \[ \bm{\Sigma} = \bm{\Lambda}\bm{\Lambda}^\prime + \bm{\Psi}, \] where $\bm{\Psi} = \mathbb{E}[\mathbf{e}_t \mathbf{e}_t^\prime]$ is the idiosyncratic covariance matrix. The distinction between exact and approximate factor models revolves around the structure of $\bm{\Psi}$: the exact model assumes idiosyncratic errors are cross-sectionally uncorrelated ($\bm{\Psi}$ diagonal), while the approximate model permits cross-sectional correlations ($\bm{\Psi}$ non-diagonal but with bounded eigenvalues).
The diagonal structure of $\bm{\Psi}$ in part (ii) defines the exact factor model: idiosyncratic errors are cross-sectionally uncorrelated. This can be relaxed:
The approximate factor structure, introduced by ChamberlainRothschild1983, relaxes the exact model by allowing weak cross-sectional correlation in the idiosyncratic errors. Under either assumption, the eigenvalues of $\bm{\Sigma}$ consist of $k$ spike eigenvalues driven by the factors and $N - k$ tail eigenvalues driven by the idiosyncratic component. Weyl's inequality gives
where $\lambda_j(\bm{\Lambda}\bm{\Lambda}') = 0$ for $j > k$, and $\psi_{\min}^2 := \lambda_{\min}(\bm{\Psi})$, $\psi_{\max}^2 := \lambda_{\max}(\bm{\Psi})$ denote the smallest and largest eigenvalues of the idiosyncratic covariance matrix $\bm{\Psi}$ (these satisfy $\underline{\psi}^2 \leq \psi_{\min}^2 \leq \psi_{\max}^2 \leq \bar{\psi}^2$, where $\underline{\psi}^2$ and $\bar{\psi}^2$ are the assumption-level bounds in Assumptions (ref)--(ref)). The key spectral implication of a factor structure is a sharp separation between factor-driven and idiosyncratic eigenvalues, formalized in the following lemma.
The proof is given in Appendix (ref). To verify the BLLT conditions, we define the tail moments \[ \bar{\lambda}_{\text{tail}} := \frac{1}{N-k}\sum_{j=k+1}^{N} \lambda_j, \qquad \overline{\lambda}^2_{\text{tail}} := \frac{1}{N-k}\sum_{j=k+1}^{N} \lambda_j^2, \] and the tail concentration ratio
with equality if and only if all tail eigenvalues are identical. The concentration ratio $\mathcal{C}$ measures the dispersion of the tail eigenvalues: $\mathcal{C} = 1$ indicates perfect tail flatness (all tail eigenvalues equal), while $\mathcal{C} \to \infty$ indicates tail concentration on a few dominant directions.
The effective ranks satisfy\footnote{From the definitions (ref): $r_k = \sum_{i>k}\lambda_i / \lambda_{k+1} = (N-k)\bar{\lambda}_{\text{tail}} / \lambda_{k+1}$. For $R_k$: $R_k = (\sum_{i>k}\lambda_i)^2 / \sum_{i>k}\lambda_i^2 = ((N-k)\bar{\lambda}_{\text{tail}})^2 / ((N-k)\overline{\lambda}^2_{\text{tail}}) = (N-k)(\bar{\lambda}_{\text{tail}})^2 / \overline{\lambda}^2_{\text{tail}} = (N-k)/\mathcal{C}$.}
These objects are the key inputs for verifying the BLLT ratios under exact and approximate factor structures.
We now verify the BLLT conditions under the exact and approximate factor structures.
The proof is given in Appendix (ref).
Conditions (i) and (ii) are satisfied by construction in a factor model, by virtue of Lemma (ref). The key to condition (ii) is that $k^* = k$: the factor structure induces a sharp spectral jump from $O(N)$ eigenvalues (the factor spikes) to $O(1)$ eigenvalues (the idiosyncratic tail) at position $k$, so the natural split index is exactly the number of factors. In finite samples this gap can be weak but asymptotically the separation diverges under both factor model assumptions.
For condition (iii), the critical requirement is that the idiosyncratic variances are not too dispersed across series. Under the exact factor model, where $\bm{\Psi}$ is diagonal, the concentration ratio $\mathcal{C}$ is determined directly by the dispersion of the diagonal entries $\psi_j^2$. If $\bar\psi^2/\underline\psi^2 = O(1)$, then $\mathcal{C} = O(1)$ and the flatness condition holds automatically as $N/T \to \infty$. Homoscedasticity ($\bm{\Psi} = \psi^2 I_N$) is the cleanest special case, giving $\mathcal{C} = 1$ and $R_k = N - k$.
Under the approximate factor model, $\bm{\Psi}$ can have off-diagonal entries reflecting residual cross-sectional correlations, and $\mathcal{C}$ is determined by the eigenvalues $\mu_1, \ldots, \mu_N$ of $\bm{\Psi}$ rather than its diagonal entries (so that $\psi_{\min}^2 = \min_j \mu_j$ and $\psi_{\max}^2 = \max_j \mu_j$). The condition $\bar{\psi}^2/\underline{\psi}^2 = O(1)$ now bounds the ratio of the largest to smallest eigenvalues of $\bm{\Psi}$. Residual correlations can inflate $\mathcal{C}$ relative to the exact factor case: if idiosyncratic components are strongly positively correlated, their joint variance concentrates into fewer directions, making the tail less flat and benign overfitting less likely. In the extreme case where all idiosyncratic components are perfectly correlated, $\bm{\Psi}$ has rank one and $\mathcal{C} \to \infty$, violating the flatness condition entirely. In practice, if the factor structure captures most of the cross-sectional dependence, the residuals will be approximately uncorrelated and $\mathcal{C}$ will be close to its exact-factor-structure value.
In Theorem (ref), the first two BLLT conditions from Condition (ref) — overall effective rank (i) and low-dimensional head (ii) — are automatically satisfied whenever the data-generating process follows a factor structure, an assumption that is reasonable in macroeconomic panels. The only empirically critical condition is the third BLLT condition (iii), the tail-flatness requirement $T/R_{k^*} \to 0$, which depends on the concentration ratio (ref).
The concentration ratio $\mathcal{C}$ governs how effectively the spectral tail disperses interpolation noise. When tail eigenvalues are uniform ($\mathcal{C}=1$), noise is spread across many comparable directions and averages out at rate $T/(N-k)$. When the tail is heterogeneous ($\mathcal{C}>1$), noise concentrates in a few dominant directions and averaging is weaker. Equivalently, $R_{k^*}=(N-k)/\mathcal{C}$ is the number of effective tail directions, and the flatness condition requires $T\mathcal{C}/(N-k)\to 0$.
Which of cases (a)--(c) in condition (iii) of Theorem (ref) applies depends on how the spectrum of $\bm{\Psi}$ behaves as $N \to \infty$. As shown in the proof of Theorem (ref), the interlacing bounds imply that the tail eigenvalues of $\bm{\Sigma}$ are asymptotically governed by the ordered eigenvalues $\{\mu_j\}_{j=1}^N$ of $\bm{\Psi}$, so the concentration ratio is asymptotically determined by the cross-sectional moments of the $\mu_j$: \[ \mathcal{C} \;\to\; \frac{\frac{1}{N-k}\sum_{j>k} \mu_j^2}{\Big(\frac{1}{N-k}\sum_{j>k} \mu_j\Big)^2} \qquad \text{as } N \to \infty. \] In case (a) we have $\mathcal{C} = 1$ and the tail is perfectly flat. We are in case (b) if and only if this limit remains bounded, and in case (c) if it diverges but slowly enough that $\mathcal{C} = o(N/T)$. If none of these cases apply---that is, if $\mathcal{C}$ grows faster than $N/T$---then condition (iii) fails and benign overfitting does not obtain.
Appendix (ref) shows that a sufficient condition for case (b) is that the eigenvalues $\{\mu_j\}$ have a Pareto tail with index $\alpha > 2$. This ensures $\mathcal{C} = O(1)$. If $1 < \alpha \leq 2$, $\mathcal{C}$ diverges but slowly enough to remain $o(N/T)$, placing us in case (c). The case $\alpha \leq 1$ violates the bounded eigenvalue assumptions of the factor model and is ruled out.
We apply these theoretical diagnostics to the FRED-MD and FRED-QD panels in Appendix (ref), estimating the tail index $\alpha$ using Hill hill1975simple and log-rank gabaix2011rank methods at various split indices $k^*$. The finite-sample evidence is suggestive but inconclusive: at moderate splits ($k^* = 15$--$30$), tail-index estimates enter the range $1 < \alpha \leq 2$ consistent with case (c), but the flatness ratio $T/R_{k^*}$ remains far from the asymptotic regime required for benign overfitting. In other words, with the available sample sizes ($N \sim 100$--$240$, $T \sim 260$--$800$), there is not enough data to reach the benign overfitting regime directly. This motivates the feature-expansion analysis in Section (ref), which achieves $N \gg T$ by construction and provides a direct implementation of the factor kernel without relying on asymptotic approximations.
Section (ref) showed that, when the data-generating process follows a factor structure, benign-overfitting conditions are satisfied in high dimensions, but in standard macro panels $N$ and $T$ are of similar order, so the condition $N\gg T$ must be engineered rather than assumed. This section addresses the issue in two steps. First, Section (ref) introduces a synthetic data augmentation procedure that artificially expands the feature dimension to achieve $N \gg T$, thereby enabling benign overfitting in finite samples. Second, Section (ref) shows that as the number of synthetic replicates grows, this augmentation converges to a closed-form kernel estimator in a Reproducing Kernel Hilbert Space (RKHS). This limiting representation is valuable because it delivers the benefits of benign overfitting without requiring practitioners to generate or store synthetic data. More broadly, the kernel duality reveals that benign overfitting, when it works, succeeds because overparameterization implicitly constructs a well-behaved kernel, not because overparameterization is intrinsically desirable.
The procedure operates in two steps. First, we estimate the factor model parameters from the observed data. Second, we use the estimated model to generate $B$ synthetic datasets by drawing fresh idiosyncratic noise while holding the estimated factor path fixed, thereby expanding the cross-sectional dimension without altering the time-series structure.
Recall from Section (ref) the factor model
where $\mathbf{f}_t \in \mathbb{R}^k$ contains the $k$ common factors, $\boldsymbol{\Lambda} \in \mathbb{R}^{N \times k}$ is the loadings matrix, and $\mathbf{e}_t \in \mathbb{R}^N$ is the idiosyncratic error with $\mathbb{E}[\mathbf{e}_t \mathbf{e}_t'] = \boldsymbol{\Psi}$. We estimate $\hat{\boldsymbol{\Lambda}}$ and the factor path $\hat{\mathbf{f}}_1,\ldots,\hat{\mathbf{f}}_T$ by principal components, with the number of factors $k$ selected via the information criterion of bai2002determining. Given the estimated factors, the idiosyncratic covariance is estimated from the residuals $\hat{\mathbf{e}}_t = \mathbf{x}_t - \hat{\boldsymbol{\Lambda}}\hat{\mathbf{f}}_t$:
giving $\hat{\boldsymbol{\Psi}} = \mathrm{diag}(\hat{\psi}_1^2,\ldots,\hat{\psi}_N^2)$.
Given $\hat{\boldsymbol{\Lambda}}$, $\hat{\boldsymbol{\Psi}}$, and the estimated factor path $\hat{\mathbf{f}}_1,\ldots,\hat{\mathbf{f}}_T$, we generate $B$ synthetic datasets by keeping the common component fixed and resampling only the idiosyncratic noise: \[ \mathbf{x}_t^{*(b)} = \hat{\boldsymbol{\Lambda}}\hat{\mathbf{f}}_t + \mathbf{e}_t^{*(b)}, \qquad \mathbf{e}_t^{*(b)} \sim \mathcal{N}(\mathbf{0}, \hat{\boldsymbol{\Psi}}), \] for $b = 1, \ldots, B$ and $t = 1, \ldots, T$. The common component $\hat{\boldsymbol{\Lambda}}\hat{\mathbf{f}}_t$ is identical across replications and carries all temporal dependence through the estimated factor path. Only the idiosyncratic noise is redrawn, independently across $t$ and $b$. We have also experimented with a bootstrap approach to resampling the errors, which allows for any correlation pattern across the residuals, and found the results to be very similar.
Each synthetic dataset $\mathbf{X}^{*(b)} = [\mathbf{x}_1^{*(b)}, \ldots, \mathbf{x}_T^{*(b)}]' \in \mathbb{R}^{T \times N}$ has the same dimensions as the original data matrix $\mathbf{X} = [\mathbf{x}_1, \ldots, \mathbf{x}_T]' \in \mathbb{R}^{T \times N}$. The augmented predictor matrix stacks the original and all $B$ synthetic datasets horizontally:
This expands the feature dimension from $N$ to $(B+1)N$ while keeping the sample size fixed at $T$. For sufficiently large $B$, the ratio $(B+1)N/T$ can be made arbitrarily large, placing the regression problem firmly in the overparameterized regime where $N \gg T$.
A Reproducing Kernel Hilbert Space (RKHS) is a Hilbert space of functions in which evaluation at any point is a continuous linear functional. Concretely, an RKHS is defined by a positive definite kernel $k(\cdot, \cdot)$: for any set of observations $\mathbf{z}_1, \ldots, \mathbf{z}_T$, the kernel matrix $\bm{K}$ with entries $K_{ts} = k(\mathbf{z}_t, \mathbf{z}_s)$ encodes pairwise similarities, and the minimum-norm interpolator in the RKHS predicts via $\hat{y}_{\mathrm{new}} = \mathbf{k}'(\mathbf{z}_{\mathrm{new}}) \bm{K}^{-1} \mathbf{y}$, where $\mathbf{k}(\mathbf{z}_{\mathrm{new}})$ collects the kernel evaluations between the new point and each training point. The key insight of this subsection is that the minimum-norm interpolator applied to the augmented system of Section (ref) is equivalent to kernel regression with a specific factor-structured kernel, and that this equivalence becomes exact as $B \to \infty$.
The most direct approach to exploiting the augmented matrix $\mathbf{Z}$ in (ref) is to estimate the regression model \[ y_t = \mathbf{z}_t'\boldsymbol{\beta} + \varepsilon_t, \qquad t = 1, \ldots, T, \] where $\mathbf{z}_t \in \mathbb{R}^{(B+1)N}$ is the $t$-th row of $\mathbf{Z}$, stacking $\mathbf{x}_t$ with the $B$ synthetic replicates $\mathbf{x}_t^{*(1)}, \ldots, \mathbf{x}_t^{*(B)}$, using the minimum-norm interpolator:
which is the ridgeless estimator (ref) applied to the $(B+1)N$-dimensional augmented system. However, this direct approach becomes computationally demanding as $B$ grows: $\hat{\boldsymbol{\beta}}$ lives in a space that expands with every new replicate. A key duality result shows that predictions from (ref) can be computed equivalently via kernel regression operating only in the $T$-dimensional observation space, without ever forming $\hat{\boldsymbol{\beta}}$ explicitly. This section establishes that equivalence and derives the limiting kernel as $B \to \infty$.
Lemma (ref) establishes the kernel duality: predictions can be computed without ever forming the $(B+1)N$-dimensional coefficient vector $\hat{\boldsymbol{\beta}}$, depending only on the $T \times T$ Gram matrix $\mathbf{Z}\mathbf{Z}'$ and kernel evaluations between the new observation and the augmented data from the whole estimation sample. For the synthetic augmentation in (ref), the Gram matrix entries are \[ (\mathbf{Z}\mathbf{Z}')_{ts} = \mathbf{x}_t^\prime\mathbf{x}_s + \sum_{b=1}^B \mathbf{x}_t^{*(b)\prime} \mathbf{x}_s^{*(b)}. \] As $B \to \infty$, the sum of synthetic inner products converges to an expectation over the idiosyncratic noise distribution, yielding a closed-form kernel that depends only on the estimated factor structure. Proposition (ref) characterizes this limiting kernel.
Proposition (ref) reveals that benign overfitting in factor models works through the implicit construction of a factor-structured kernel. The factor structure ensures the kernel preserves the low-dimensional signal (the common component) while the tail flatness condition ensures noise is dispersed across many effective directions. The split index $k^*$ identifies the boundary between the signal subspace (which the kernel must preserve) and the noise-dispersing tail (which the kernel must spread evenly), while the tail flatness measure $R_{k^*}$ quantifies how effectively the kernel disperses noise. This explains why overparameterization succeeds when it does: not because adding parameters is intrinsically beneficial, but because it constructs a kernel that separates signal from noise in an appropriate way.
The kernel (ref) measures temporal proximity through the estimated factor structure rather than calendar time: periods $t$ and $s$ receive high weight when their factor realizations $\hat{\mathbf{f}}_t$ and $\hat{\mathbf{f}}_s$ are close in the metric induced by $\hat{\boldsymbol{\Lambda}}'\hat{\boldsymbol{\Lambda}}$. The ridge parameter $\lambda = \textup{tr}(\hat{\boldsymbol{\Psi}})$ provides automatic regularization scaled by idiosyncratic variance. This is different from the kernel estimators typically used in the econometrics literature, as those are local methods which borrow from adjacent dates, while this kernel estimator is a global method which borrows from economically similar dates, regardless of calendar distance. This property is valuable in regime-dependent environments where comparable macroeconomic conditions reappear non-consecutively.
Finally, since the augmentation strategy converges to kernel ridge regression with a closed-form kernel, practitioners can implement the kernel estimator directly without generating synthetic data. This avoids the computational cost and Monte Carlo variability of finite-$B$ augmentation. A detailed RKHS formulation and the connection to the forecasting implementation are provided in Appendix (ref).
We now evaluate the factor kernel of Proposition (ref) in a comprehensive rolling-window out-of-sample forecasting exercise covering all available series in the FRED-MD dataset. Unlike the motivating example of Section (ref), which used a single holdout period to trace out the full double-descent curve, this will constitute a proper pseudo real-time forecasting evaluation: every forecast will be computed using only information available at the time it is made, the window will roll forward one observation at a time, and performance will be measured by the mean squared forecast error (MSFE) accumulated over the full evaluation sample.
We use the February 2026 vintage of the FRED-MD monthly dataset mccracken2016fred, which contains $N_{\mathrm{raw}} = 127$ monthly macroeconomic and financial series. All series are transformed using the mccracken2016fred transformation codes (log differences, second differences, etc., as indicated by tcode). We end the sample on October 2025 to avoid real-time data revision issues that are most severe at the end of the sample.\footnote{Series with more than 20% missing observations are excluded; remaining missing predictor values are imputed with column medians computed in the training window. After cleaning, a typical target uses about $N \approx 120$ predictors.}
For each of the $N_{\mathrm{raw}}$ series and each forecast horizon $h \in \{1, 3, 6, 12\}$ months, we compute direct $h$-step-ahead forecasts. The general forecasting model is
where $y_t$ is the target variable available at forecast origin $t$, and $\mathbf{x}_t$ is the predictor vector observed at the same origin, obtained by pulling together all the remaining series in the panel. The AR(1) term $\phi\, y_t$ is common to all models; they differ only in the specification of $f(\cdot)$. Because the forecast is direct (we are forecasting $y_{t+h}$), the model is estimated separately for each horizon $h$. We use a rolling window estimation with $W = 120$ months. At forecast origin $t$, the model is estimated on observations $s=t-W-h+1,\ldots,t-h$. With data spanning January 1959 to October 2025, the first forecast origin is January 1969 for $h=1$ (forecasting February 1969), providing approximately 56 years of out-of-sample evaluation. The number of out-of-sample forecasts is $P=T-h-W$.
At each forecast origin, variables are standardized using moments from the current training window only: predictor means and standard deviations for $\mathbf{x}_s$, and separate mean and standard deviation for the dependent variable and its lag. Standardization is recomputed at each origin. Squared forecast errors are accumulated in standardized units, so MSFE is dimensionless and comparable across targets.
For each model $M$ and each target--horizon pair $(y, h)$, the out-of-sample MSFE is
computed in standardized units. The relative performance of a model $M$ against a benchmark $B$ is summarized by the RMSE ratio
where values below one indicate that model $M$ improves over the benchmark $B$. In our evaluation, $M = \mathrm{FK}$ and $B = \mathrm{FM}$.
We compare two models. The first one, which will serve as our benchmark, is a factor model estimated via OLS (FM). This is a special case of (ref) with $f(\mathbf{x}_t) = \boldsymbol{\gamma}_h'\hat{\mathbf{f}}_t$, a linear function of the factors extracted from the panel of predictors. This is the diffusion-index forecast of stock2002forecasting, augmented with an AR(1) term. For each origin $t$ and horizon $h$, we estimate by OLS \[ y_{s+h} = \phi_h y_s + \boldsymbol{\gamma}_h'\hat{\mathbf{f}}_s + u_{s+h}, \qquad s=t-W-h+1,\ldots,t-h, \] where $\hat{\mathbf{f}}_s$ are the factors extracted from the rolling window using principal components. The forecast is \[ \hat{y}^{\mathrm{FM}}_{t+h\mid t} = \hat\phi_h y_t + \hat{\boldsymbol{\gamma}}_h'\hat{\mathbf{f}}_t. \] At each forecast origin, the number of factors is selected with the Bai--Ng $\mathrm{IC}_{p2}$ criterion bai2002determining from the standardized predictor panel in the rolling window, with $k_{\max}=20$. This factor-selection step determines the number of factors used in both models.
Our second model is a semiparametric factor-based kernel regression estimated via GLS (FK). This is also a special case of (ref) with $f(\mathbf{x}_t)$ replaced by a nonparametric kernel component, the closed-form limit of the synthetic augmentation strategy from Proposition (ref). As derived below, this yields $f(\mathbf{x}_t) = \hat{\mathbf{k}}_{t,h}'\hat{\boldsymbol{\alpha}}_h$, where $\hat{\mathbf{k}}_{t,h}$ collects the kernel evaluations between origin $t$ and the training-window dates; the RKHS derivation is in Appendix (ref).
For a fixed target, horizon $h$, and forecast origin $t$, define \[ \mathbf{y}_h = (y_{t-W+1+h},\ldots,y_{t+h})', \qquad \mathbf{y}_{-1} = (y_{t-W+1},\ldots,y_t)'. \] Let $\hat{\mathbf{K}}_h$ be the $W\times W$ estimated factor kernel matrix over window dates, with entries $\hat{\mathbf{K}}_h(s,r)=\hat{\mathbf{f}}_s'\hat{\boldsymbol{\Lambda}}'\hat{\boldsymbol{\Lambda}}\hat{\mathbf{f}}_r$, and let \[ \hat{\mathbf{K}}_h^{(\lambda)} = \hat{\mathbf{K}}_h + \lambda \mathbf{I}_W, \] where $\lambda = \mathrm{tr}(\hat{\boldsymbol{\Psi}})$ is the ridge parameter determined by the trace of the estimated idiosyncratic covariance matrix, as established in Proposition (ref). The semiparametric factor kernel model with GLS AR(1) is \[ \mathbf{y}_h = \phi_h \mathbf{y}_{-1} + \hat{\mathbf{K}}_h\boldsymbol{\alpha}_h + \mathbf{u}_h. \] Using the Moore--Penrose inverse of $\hat{\mathbf{K}}_h^{(\lambda)}$, the GLS estimates are
For forecast origin $t$, define the kernel vector to the window dates, \[ \hat{\mathbf{k}}_{t,h} = \big(\hat{\mathbf{K}}_h(t,t-W+1),\ldots,\hat{\mathbf{K}}_h(t,t)\big)'. \] Then the direct forecast is
Larger idiosyncratic dispersion (larger $\lambda$) implies stronger regularization automatically. The key difference from the benchmark is that the factor kernel measures temporal proximity through the estimated factor structure: periods $t$ and $s$ receive high weight when their factor realizations $\hat{\mathbf{f}}_t$ and $\hat{\mathbf{f}}_s$ are close, regardless of calendar distance. This state-dependent pooling is valuable in regime-dependent environments where comparable macroeconomic conditions reappear non-consecutively.
Figure (ref) shows the cross-sectional distribution of $\mathrm{Imp}_{y,h}$ across all FRED-MD series, separately for each horizon. A distribution skewed to the left of one indicates that the factor kernel systematically outperforms the FM benchmark in the cross-section. The vertical dashed line marks the break-even point.
\paragraph{Does the factor kernel consistently beat the FM benchmark?} The FM model is a strong baseline: it uses the same factors and the same AR(1) term as FK, but combines them through linear OLS rather than kernel-weighted GLS. Any FK advantage therefore comes from the kernel metric and from state-dependent pooling (dates with similar factor realizations receive more weight), as discussed in Section (ref) and formalized in equations (ref)--(ref).
Figure (ref) shows that, across all FRED-MD variables, the distribution of $\mathrm{Imp}_{y,h}$ is centered to the left of one at every horizon considered: the factor kernel delivers a lower MSFE than the FM benchmark for 92% of targets at $h=1$ and 99% at $h=6$ and $h=12$. The mean RMSE ratio is 0.99 at $h=1$, indicating a 1% mean improvement. While the magnitude of improvement is modest, it is highly consistent across targets and becomes increasingly pervasive at longer horizons. Statistical significance testing using the Diebold--Mariano test (two-sided, 5% level) reveals that 34% of FK wins at $h=1$ are statistically significant, rising to 54% at $h=6$ and 50% at $h=12$. The increase in both win rate and statistical significance at longer horizons indicates that FK's advantage becomes both more reliable and more pervasive as the forecast horizon lengthens.
\paragraph{Which variables benefit most from the factor kernel?} To interpret and understand the heterogeneity in the performance of the factor kernel, we organize the FRED-MD variables in terms of their persistence (lag-1 autocorrelation of the standardized variable in the second half of the sample) and map this against relative out-of-sample performance. Figure (ref) shows a clear negative relationship between persistence and RMSE ratio, with the correlation strengthening from $-0.34$ at $h=1$ to $-0.54$ at $h=6$. Highly persistent variables benefit substantially more from the factor kernel's state-dependent pooling: the kernel identifies economically similar states across time, which is particularly valuable when the target exhibits strong temporal dependence. FK's advantage becomes increasingly pervasive (winning for nearly all targets) and statistically robust at longer horizons, even as the magnitude of improvement becomes more modest and uniform across targets. Persistence maintains its role as the primary predictor of FK advantage even at $h=12$.
In this section, we repeat the forecasting evaluation exercise above using the FRED-QD quarterly dataset mccracken2016fred. The quarterly data provide a complementary perspective on the factor kernel's performance in a setting where the cross-sectional dimension systematically exceeds the time-series dimension, placing the estimation problem squarely in the interpolation regime throughout the evaluation period.
The February 2026 vintage of FRED-QD contains $N_{\mathrm{raw}} = 248$ quarterly macroeconomic and financial series. After excluding series with more than 20% missing observations and imputing remaining missing values with training-window medians, we are left with approximately $N \approx 240$ predictors. The sample spans 1960:Q1 to 2025:Q3 and the rolling window length is set to $W = 40$ quarters (10 years). Forecast horizons are $h \in \{1, 2, 4, 8\}$ quarters. The first forecast origin is 1969:Q1, providing approximately 56 years of out-of-sample evaluation. Because $N \approx 240 \gg W = 40$, every rolling window lies in the overparameterized regime $N > T$, so the ridgeless estimator fits the training data exactly at each forecast origin. This contrasts with the monthly data, where $N \approx 120 < W = 120$ in most windows, placing the estimation in the underparameterized regime.
The models, estimation procedures, and performance metrics are identical to those described in Sections (ref) and (ref). We rely on the Bai--Ng $\mathrm{IC}_{p2}$ criterion to select the number of factors at each origin with $k_{\max} = 20$. The factor kernel is constructed from the estimated factors and loadings as in Proposition (ref), and the ridge parameter $\lambda = \mathrm{tr}(\hat{\boldsymbol{\Psi}})$ is determined by the idiosyncratic variance. The relative performance measure $\mathrm{Imp}_{y,h} = \mathrm{RMSE}_{\mathrm{FK}} / \mathrm{RMSE}_{\mathrm{FM}}$ is computed for each target and horizon.
Figure (ref) displays the cross-sectional distribution of $\mathrm{Imp}_{y,h}$ across all FRED-QD series, separately for each horizon. As in the monthly case, a distribution concentrated to the left of one indicates systematic outperformance by the factor kernel. The results show that the factor kernel delivers lower MSFE than the FM benchmark for 96% of targets at $h=1$, $h=2$, and $h=8$, and 99% at $h=4$, with the mean RMSE ratio of 0.95 at $h=1$ indicating a 5% mean improvement. The mean improvement is most pronounced at short horizons and attenuates to approximately 3% at longer horizons, consistent with the well-documented decline in factor-model predictability as the forecast horizon lengthens. Statistical significance testing reveals that 9% of FK wins at $h=1$ are statistically significant by the Diebold--Mariano test, rising to 30% at $h=8$, indicating that the improvements become more statistically robust at longer horizons.
Figure (ref) shows the relationship between persistence and relative out-of-sample performance. Target persistence emerges as the primary predictor of FK advantage, with the correlation between persistence and RMSE ratio strengthening from $-0.17$ at $h=1$ to $-0.32$ at $h=8$. Highly persistent targets benefit more from the factor kernel's state-dependent pooling, consistent with the pattern observed in the monthly data. The relationship is strongest at longer horizons, where the kernel's ability to identify economically similar states becomes most valuable.
To assess whether the factor kernel's advantage is stable across different economic regimes, we partition the 56-year evaluation period into three 20-year windows: 1969--1989 (pre-Great Moderation), 1989--2009 (Great Moderation and Financial Crisis), and 2009--2025 (post-Crisis era). For each target, horizon, and time period, we recompute the MSFE using only forecasts whose target dates fall within that period, and we test for statistical significance using the Diebold--Mariano test applied to the subsample of forecast errors.
Figure (ref) shows the results for FRED-MD. Each panel corresponds to one forecast horizon, and each point represents one target variable. The horizontal axis reports target persistence (lag-1 autocorrelation), and the vertical axis reports the RMSE ratio $\mathrm{RMSE}_{\mathrm{FK}} / \mathrm{RMSE}_{\mathrm{FM}}$ computed over the specified time period. Points are color-coded by period: orange for 1969--1989, blue for 1989--2009, and green for 2009--2025. Filled circles indicate that the difference is not statistically significant at the 5% level; crosses indicate statistical significance.
The results reveal that the factor kernel's advantage is broadly stable across time periods, with the negative relationship between persistence and RMSE ratio evident in all three subsamples. At $h=1$, the FK win rates are 91%, 94%, and 89% for the three periods respectively, indicating consistent performance over time. The mean improvements are 0.3%, 0.5%, and 1.1%, suggesting that the magnitude of FK's advantage has grown in the most recent period. Statistical significance is lower in recent periods, with 12% of wins significant in 2009--2025 compared to 32% in 1969--1989, reflecting the shorter sample period and lower volatility in the post-Crisis era.
Figure (ref) presents the corresponding results for FRED-QD. The quarterly data, which place the estimation in the interpolation regime throughout ($N \gg W$), show similar patterns. The FK win rates are 87%, 85%, and 95% for the three periods at $h=1$, with mean improvements of 0.7%, 1.2%, and 6.7%. The persistence-based predictability of FK advantage remains strong across all periods, with correlations ranging from $-0.17$ to $-0.32$. Notably, the mean improvement is substantially larger in the most recent period (2009--2025), suggesting that the factor kernel's state-dependent pooling has become increasingly valuable in the post-Crisis era.
Overall, the subsample analysis confirms that the factor kernel's systematic outperformance is not driven by a particular economic regime or historical episode. The advantage persists across the high-volatility pre-Great Moderation period, the relatively stable Great Moderation, the Financial Crisis, and the post-Crisis era. This robustness across diverse macroeconomic conditions strengthens the case that the factor kernel's state-dependent pooling mechanism provides a genuine improvement over linear factor models, rather than exploiting regime-specific features of the data.
This paper makes three contributions. First, it documents clear double-descent patterns in macroeconomic forecasting data. In the in-sample motivating exercise, the MSFE rises near interpolation and then falls in the overparameterized region, with meaningful differences across variables depending on how strongly each series aligns with the common-factor structure. This provides direct empirical evidence that the double descent phenomenon is relevant in macro panels and not only in canonical machine-learning datasets.
Second, it develops a theoretical link between benign overfitting and factor models. Under exact factor structure, the spectral conditions for benign overfitting are naturally satisfied. Under approximate factor structure, the same mechanism can still operate, but its strength depends on the concentration of the idiosyncratic spectrum. This clarifies when interpolation can be benign in macroeconomic environments and why the answer is target dependent.
Third, it shows that synthetic data augmentation has a precise kernel interpretation and therefore yields an implementable forecasting method. The key insight is that benign overfitting in factor models works by implicitly defining a factor-structured kernel: the augmentation strategy converges to kernel ridge regression with a closed-form kernel that measures temporal proximity through the estimated factor structure rather than calendar time. This kernel can be implemented directly without generating synthetic data, avoiding computational cost and Monte Carlo variability.
In the out-of-sample rolling-window evaluation, the factor-kernel model systematically outperforms the standard Stock-Watson factor model across both monthly and quarterly datasets at all horizons, with win rates exceeding 90%. The improvements are modest but highly consistent: mean gains are approximately 1% in the monthly data (underparameterized regime) and 5% in the quarterly data (interpolation regime). Statistical significance increases substantially with forecast horizon, with approximately half of wins becoming statistically significant at longer horizons in the monthly data and 30% in the quarterly data. The evaluation period spans approximately 56 years (1969-2025), providing robust evidence of the factor kernel's performance across diverse macroeconomic conditions. The quarterly results provide direct evidence that the factor kernel performs well in the interpolation regime, where the minimum-norm estimator fits training data exactly yet generalizes effectively when the predictor covariance has the spectral structure implied by factor models.
Across both datasets, we find that highly persistent target series benefit substantially more from the factor kernel's state-dependent pooling, which identifies economically similar states across time regardless of calendar distance. This suggests that the kernel's ability to pool information across similar economic conditions is particularly valuable when the target exhibits strong temporal dependence.
An interesting direction for future research concerns the relationship between signal alignment and forecast performance. While the BLLT theory predicts that targets with signal concentrated in the leading eigenspace should benefit most from the factor kernel, the empirical relationship in our data is weak. Understanding in which situations the factor kernel outperforms linear methods, and reconciling the role of signal alignment versus state-dependent pooling in driving forecasting performance, remains an open question.