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.
70,735 characters · 16 sections · 40 citation commands
Probabilistic Targeted Factor Analysis
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
Keywords: Dynamic models, Expectation-Maximization algorithm, Factor extraction, High-dimensional data, Missing data.
JEL Classification: C38, C53, C55, G12, G17
\spacingset{1.7}
Empirical economic analysis increasingly relies on large panels of predictors to measure latent economic quantities such as financial conditions, aggregate risk, or macroeconomic slack (e.g. JuradoLudvigsonNg2015,Petrella2017, Stock2016). While modern datasets provide rich information, the relevant economic object is often not the dominant source of variation in the data, creating a fundamental measurement challenge. Dimensionality-reduction methods such as Principal Component Analysis (PCA) summarize overall variance, but are silent about which components are economically relevant for a given outcome of interest.
Partial Least Squares (PLS), originally developed by wold_soft_1975, addresses this limitation by constructing latent factors that maximize covariance between predictors and a set of target variables. This targeted nature has led to widespread adoption of PLS in economics and finance, particularly in high-dimensional forecasting environments where multicollinearity is pervasive welch_comprehensive_2007, kelly_market_2013, kelly_three-pass_2015, giglio_systemic_2016, groen_revisiting_2016, goyal_comprehensive_2024. However, standard PLS remains an algorithmic, deterministic procedure, offering no probabilistic interpretation of the latent factors or a principled treatment of uncertainty, missing data, or time-series structure.
In this paper, we develop Probabilistic Targeted Factor Analysis (PTFA), a likelihood-based framework for supervised factor extraction that provides a probabilistic foundation for PLS. PTFA learns latent factors jointly from predictors ($\mathbf{X}$) and targets ($\mathbf{Y}$) by maximizing the observed-data likelihood via an expectation–maximization (EM) algorithm. Unlike the two-step procedure underlying classical PLS regression,\footnote{The Nonlinear Iterative Partial Least Squares (NIPALS; wold_soft_1975) algorithm extracts latent components and subsequently regresses $\mathbf{Y}$ on these components; see butler_peculiar_2002.} PTFA explicitly models the stochastic data-generating process of both latent factors and observed variables. This probabilistic formulation yields latent representations that can be interpreted as economically relevant measurement objects, while remaining computationally efficient and scalable in high-dimensional settings.
Much like the probabilistic reformulation of PCA proposed by tipping_probabilistic_1999, PTFA embeds a widely used dimensionality-reduction technique into a coherent statistical framework. This allows uncertainty to be formally quantified and naturally accommodates features commonly encountered in economic data, including missing observations, mixed-frequency sampling, stochastic volatility, and time-series persistence. Estimation proceeds via a fast EM algorithm, and we provide an open-source implementation to facilitate empirical adoption \footnote{currently available in Python: \href{https://pypi.org/project/ptfa/}{pypi.org/project/ptfa}.}. Probabilistic variants of PLS have previously been proposed in other literatures, particularly in Chemometrics and Statistics gustafsson_probabilistic_2001, li_probabilistic_2011, zheng_probabilistic_2016, el_bouhaddani_probabilistic_2018. These approaches primarily focus on covariance maximization or algorithmic reconstruction in cross-sectional settings. In contrast, PTFA is designed as a targeted measurement framework for economic time series, directly optimizing the predictive likelihood of economically meaningful targets. This distinction allows PTFA to integrate uncertainty, time dependence, and incomplete data in a unified probabilistic structure, while retaining interpretability and computational tractability.
Our contribution is threefold. First, we provide a probabilistic characterization of targeted factor extraction, establishing formal parallels with probabilistic PCA and identifying solution structures unique to the supervised setting. Second, the probabilistic framework enables extensions that are central to applied economic analysis, including missing data, mixed-frequency designs, stochastic volatility, and dynamic factor behavior. Third, we demonstrate the empirical relevance of PTFA through simulation evidence and applications in macroeconomics and finance, showcasing its in-sample and out-of-sample outperformance.
Over the past decades, factor analysis has become a cornerstone of financial and macroeconometric research Stock2016. Building on the EM-based estimation of approximate dynamic factor models doz_two-step_2011, doz_quasimaximum_2012, barigozzi2024, our framework is most closely related to chan_bayesian_2019, from which it departs by explicitly targeting a subset of variables of economic interest and jointly modeling predictors and targets. We adopt a block-banded representation of the factor state equation, accommodating local cross-sectional dependence and enhancing computational performance in large systems, while remaining fully compatible with EM-based inference jungbacker_koopman_2015, coroneo2016.
In the empirical section, we explore novel insights in three applications to Economics and Finance allowed by the targeted factor setup we introduce. The first application uses the dynamic PTFA framework to construct targeted Financial Conditions Indices, directly addressing identification and interpretability critiques of conventional FCIs. We provide evidence that this novel approach to measuring financial conditions addresses key critiques in the literature related to identification of FCIs. In addition, we show indices are easy (and fast) to construct, and can target the dynamics of any key variable (or combination of variables) of interest to policymakers. In a second application, we employ PTFA to forecast three key macroeconomic variables: industrial production, consumer price index (CPI) inflation and unemployment. We conduct a forecasting exercise that harvests the information contained in 126 Federal Reserve Economic Data monthly time series (FRED-MD) mimicking the setup in mccracken_fred-md_2016. Lastly, we use our model to predict the equity premium using 26 signals made available by goyal_comprehensive_2024. In all applications PTFA generally outperforms PLS and PCA at out-of-sample forecasting of key targets, at similar computational costs.
The remainder of the paper is organized as follows. Section (ref) provides a refresher on PLS and outlines our probabilistic foundation resulting in the PTFA method. Section (ref) extends the baseline framework and EM algorithm to deal with real-world data complexities. Section (ref) presents the setup for our simulation exercises and discusses the findings. In Section (ref), we provide three applications of PTFA using popular datasets in Economics and Finance. Section (ref) provides our concluding remarks. The appendices provide additional theoretical and implementation details on our method: Appendix (ref) provides a theoretical analysis of the properties of the maximum likelihood estimator resulting from the observed-data likelihood of our probabilistic model; Appendix (ref) derives an EM algorithm based on iterative maximization of this likelihood; Appendix (ref) provides pseudo-code for all algorithms implemented in the \verb|PTFA| package; and additional results are placed in Appendix (ref).
Throughout this section, we let $\bm{x}$ be a $p$-dimensional vector of features and $\bm{y}$ be a $q$-dimensional vector of prediction targets. Our main assumption is that there are $k$ common components or factors collected in vector $\bm{f}$, where one typically expects $p \gg k$. This is generally the case when one aims to extract signals from high-dimensional data with a large feature space.
To understand the motivation behind the development of the PTFA, it is helpful to compare it with PCA. Both PTFA and PCA aim to reduce the dimensionality of a large set of variables, via a dense low-rank representation of the data-generating process (DGP). However, they do so with different objectives in mind. PCA transforms a set of possibly correlated variables into a set of linearly uncorrelated latent vectors called principal components. These components are constructed in an unsupervised way, such that the first principal component captures the maximum variance in the data, the second captures the maximum variance after projecting out the first component, and so on. PTFA, on the other hand, constructs latent vectors that maximize the covariance between the predictor ($\mathbf{X}$) and response ($\mathbf{Y}$) variables, targeting the former. Unlike PCA, which focuses solely on the predictors, PTFA considers both predictors and responses, aiming to find the directions in the features that best predict the responses. This makes PTFA particularly attractive in contexts where interpretability of the latent vectors is important.
Our setup is akin to the DGP assumed in the theoretical results supporting the popular forecasting technique Three-Pass Regression Filter kelly_three-pass_2015, of which PLS is a special case. We additionally assume throughout the paper that both $x$ and $y$ have been standardized and purged of the influence of any common observable effects.\footnote{Let all common observable effects between targets and features be captured by $\bm{z}$, an $r$-dimensional vector of controls (including a constant for de-meaning). We can then assume that $\bm{x}$ and $\bm{y}$ are the residuals from a linear projection of the original features and targets onto $z$.} Before proceeding with PTFA, we also provide a short summary of the intuition underlying PLS, for which PTFA provides a probabilistic foundation.
Partial Least Squares regression aims to find $k$ factors from $\bm{x}$ that best predict $\bm{y}$. While techniques like PCA also perform dimensionality reduction on the set of potential predictors $\bm{x}$, they are silent about the relevance of the principal components to predict $\bm{y}$. On the contrary, PLS directly recovers scores from $\bm{x}$ with predictability of $\bm{y}$ in mind. To be specific, PLS regression searches for two sets of scores $\bm{f}_x$ and $\bm{f}_y$ that perform a simultaneous decomposition of $\bm{x}$ and $\bm{y}$ such as to maximize the correlation between $\bm{y}$ and the recovered $\bm{f}_x$. By focusing on recovering only $k$ scores, we project the feature space to the directions that maximize predictability of $\bm{y}$ in the mean-squared error sense. Therefore, the independent and dependent variables are decomposed as linear transformations of the scores with loadings $\mathbf{P}$ and $\mathbf{Q}$ such as
Loadings $\mathbf{P}$ and $\mathbf{Q}$ are chosen so as to maximize $\operatorname{Cov}(\bm{x}, \bm{y}) = \operatorname{Cov}(\mathbf{P} \bm{f}_x, \mathbf{Q} \bm{f}_y)$. The Nonlinear Iterative Partial Least Squares (NIPALS) algorithm wold_soft_1975 --- commonly used in the literature to estimate and motivate PLS --- can then be used to efficiently recover the loadings, and these can be used to forecast $\bm{y}$ for any value of $\bm{x}$. Note that, in contrast to PCA, the loadings recovered from PLS are not necessarily orthogonal.
Importantly, note that this representation does not have a probabilistic foundation in mind as there is no randomness embedded into the factor or loading recovery process. Additionally, the factor extraction process for these techniques is usually thought of in an algorithmic or geometric manner, rather than a statistical one. This also means a standard PLS approach does not acknowledge additional sources of variation in the data such as noise or incomplete data patterns.
We now provide a simple statistical formulation for performing factor extraction that embodies the objective of PLS. This framework, which we denote as Probabilistic Targeted Factor Analysis (PTFA), provides a simple unifying setting to probabilistic extraction of factors from features $\bm{x}$ that optimally predict a pre-specified target $\bm{y}$. We assume the following model for $\bm{x}$ and $\bm{y}$ as generated from some common latent components $\bm{f}$ as
where $\mathbf{P}$ and $\mathbf{Q}$ again represent loadings, and $\bm{e}_x \sim \mathcal{N}_p(\mathbf{0}_p, \sigma^2_x \mathbf{I}_p)$ and $\bm{e}_y \sim \mathcal{N}_q(\mathbf{0}_q, \sigma^2_y \mathbf{I}_q)$ are isotropic Gaussian noise terms.\footnote{Zero-mean errors without a constant in the model are without loss of generality as both $\bm{x}$ and $\bm{y}$ are centered prior to any processing.}
The structure provided by Eqs. (ref)--(ref) and the normality of the errors is not to be taken as a literal distributional assumption, but rather a probabilistic framework to embed PLS and similar targeted factor extraction techniques. The assumption of isotropic Gaussian errors is made for convenience as it results in a particularly simple likelihood structure; it can be easily extended to full covariance matrices. If we assume the factor decomposition to be correctly specified in expectation ---i.e., that there exists $\mathbf{P} \in \mathbb{R}^{p \times k}$ and $\mathbf{Q} \in \mathbb{R}^{q \times k}$ such that $\mathbb{E}[\mathbf{X} \mid \mathbf{F}] = \mathbf{F}\mathbf{P}^{\top}$ and $\mathbb{E}[\mathbf{Y} \mid \mathbf{F}] = \mathbf{F}\mathbf{Q}^{\top}$ --- then, under regularity conditions, our estimator will remain consistent even if the data is subject to other kind of more complex non-Gaussian error processes gong_pseudo_1981, gourieroux_pseudo_1984. Furthermore, results by doz_two-step_2011, doz_quasimaximum_2012, barigozzi2024 guarantee consistency of maximum likelihood estimation (MLE) of factors even in the case of approximate factor models in which the decomposition does not hold in expectation, provided the dimensionality of predictors ($p$ in our notation) is large or growing with sample size.
Similar to the 3PRF, PTFA considers predictors $\bm{x}$ and targets $\bm{y}$ to be correlated solely through the latent common factors they share, such that $\operatorname{Cov}(\bm{e}_x, \bm{e}_y) = \mathbf{0}_{p \times q}$. This choice of correlation structure results in a particularly simple yet consistent (quasi-)maximum likelihood estimator (QMLE) of the loadings $\mathbf{P}$ and $\mathbf{Q}$ (see Appendix (ref) for details). Other models, such as that considered by groen_revisiting_2016, explicitly allow for additional predictors to be directly correlated with the target variable. While we can easily allow for additional correlation between features and targets conditional on the factors, for simplicity of exposition this paper focuses on the case where no additional correlation is assumed.
The latent scores are assumed to be normally distributed $\bm{f} \sim \mathcal{N}_k(\bm{0}_k, \mathbf{V}_F)$ with positive-definite prior variance $\mathbf{V}_F$. We simply set $\mathbf{V}_F = \mathbf{I}_k$ in the absence of any prior information on the latent scores or if they are meant to represent structurally uncorrelated components. Letting $d \coloneqq p + q$, and using the conditional independence between $\bm{x}$ and $\bm{y}$, the conditional likelihood is given by
a multivariate Gaussian with $d$-dimensional mean vector $\bm{\mu} \coloneqq [\bm{f}^{\top} \mathbf{P}^{\top}, \bm{f}^{\top} \mathbf{Q}^{\top}]^{\top}$ and $d \times d$ variance-covariance matrix $\bm{\Sigma} \coloneqq \operatorname{diag}(\sigma^2_x \mathbf{I}_p, \sigma^2_y \mathbf{I}_q)$. To derive the posterior distribution $p(\bm{f}|\bm{x},\bm{y})$ we use Bayes' rule to find
Completing the squares, we can derive the factor posterior as $\bm{f} \mid \bm{x}, \bm{y} \sim \mathcal{N}_k(\bm{m}, \bm{\Omega})$ with posterior mean and covariance matrix given as
Details on implementation can be found in Appendix (ref).
In this section, we preview a host of possible extensions that become simple to implement once one has a probabilistic framework for targeted factor analysis. Specifically, we provide extensions to incomplete data (both under at-random and mixed-frequency designs), stochastic volatility in the features and targets, and time-series persistence in the latent factors driving co-movements. All of these are of particular interest in areas such as economics and finance, where data is likely to exhibit such patterns.
While we focus on providing simple and computationally efficient extensions for PTFA through the EM approach, PTFA could also be augmented to variational or Bayesian inference by specifying priors over $\theta$, providing access to all the benefits of these frameworks. This means, for example, sparsity in the loadings via shrinkage priors, structural consideration of stochastic volatility, variable selection and model uncertainty, among many others. Some have been partially considered in the literature klami2013bayesian, vidaurre_bayesian_2013, li_process_2018, zheng_semisupervised_2018, xie_fault_2019, yang_robust_2021, el_bouhaddani_statistical_2022, and we will continue exploring such extensions in future research.
PTFA offers a natural approach to the estimation of the principal axes in cases where some of the data in $\mathbf{X}$ and $\mathbf{Y}$ are missing at random. We follow standard methodology for maximizing the likelihood of a Gaussian model in the presence of missing values little_statistical_2019. We explain now the changes that we make to the standard algorithm to account for incomplete data.
For any row $t \in \{1, \ldots, T\}$ and feature $j \in \{1, \ldots, p\}$, we let $\tau_{tj}^{(X)}$ be an indicator for whether that particular observation is missing in the feature matrix $\mathbf{X}$. That is, $\tau_{tj}^{(X)} = 1$ only when observation $tj$ is missing, and is 0 otherwise. We can analogously define a missing indicator for the target matrix $\mathbf{Y}$ and denote it as $\tau_{tj}^{(Y)}$, where $j \in \{1, \ldots, q\}$.
Note that, after standardization, the unconditional mean for all columns of $\mathbf{X}$ and $\mathbf{Y}$ is zero. Therefore, a natural initial imputation strategy is to replace all missing observations (those with $\tau_{tj}^{(X)} = 1$ or $\tau_{tj}^{(Y)} = 1$) by 0. Let the matrices with imputed values be denoted as $\widetilde{\mathbf{X}}$ and $\widetilde{\mathbf{Y}}$. One pass of our EM algorithm allows us to obtain estimated values for $\widetilde{\mathbf{P}}$ and $\widetilde{\mathbf{Q}}$ as well as the predicted scores by using Eq. (ref) and (ref) on the imputed matrices, with the output denoted as $\widetilde{\mathbf{M}}$. Finally, given these values, we update our initial imputation of $\widetilde{\mathbf{X}}$ and $\widetilde{\mathbf{Y}}$ to a more accurate one by setting
By iteratively applying equations (ref), (ref), and (ref), our algorithm can adapt to relatively large contamination rates for both features and targets. This feature of our method is also explored through a Monte Carlo simulation in the next section.
Many economic time series do not follow common release schedules and information availability itself can change across time. This creates a similar missing-data problem, but not one where we can claim the data is missing at random given there is a clear pattern to the unobserved data points. As a running example, consider $\bm{x}$ to be a set of monthly indicators, where our targets $\bm{y}$ are economic indicators (such as GDP or inflation) available at a quarterly frequency. As one cannot observe the quarterly variables in the intermediate months for which one has feature data available, the pattern of missingness is clearly not random.
In order to provide real-time indices of economic activity under the problem of mixed-frequency data, many proposal have been introduced in the literature foroni_mixed-frequency_2014. We now show how the PTFA framework can be extended to account for mixed-frequency observations, particularly when the missing data is on the target. We can modify our initial model equations as follows by introducing a latent set of targets $\bm{y}^*$ that is available at the higher frequency of $\bm{x}$. By aggregating the latent variable into the lower-frequency analog, our model produces a likelihood that can directly be used to modify the EM algorithm introduced previously:
for $l = 1, \ldots, L$, where $L$ represents the amount of period of high-frequency observations with respect to the lower-frequency ones (e.g., if $\bm{x}$ is monthly and $\bm{y}$ is quarterly, then $L = 3$). This aggregation scheme is drawn from the parsimonious approach to mixed-frequency modeling advocated by giannone_nowcasting_2008. Combining equations (ref) and (ref), we can directly derive the following equation for the observable target in terms of the latent factors:
Assuming as before that $\bm{f}^{(l)} \sim \mathcal{N}_k(\bm{0}_k, \mathbf{V}_F)$ and $(\bm{e}_x, \bm{e}_y^{*(l)}) \sim \mathcal{N}_d(\bm{0}_d, \bm{\Sigma})$ independently across $l = 1, \ldots, L$, the conditional likelihood can be expressed as
Define the full set of features $\bm{x} = [\bm{x}^{(1)\top}, \ldots, \bm{x}^{(L)\top}]^{\top}$ and factors $\bm{f} = [\bm{f}^{(1)\top}, \ldots, \bm{f}^{(L)\top}]^{\top}$, of size $p \cdot L$ and $k \cdot L$, respectively. That is, these vectors collect the higher-frequency observations to match the lower frequency of the observed target $\bm{y}$. As before, we can use Bayes' rule to obtain a posterior for the factors as $\bm{f} \mid \bm{x}, \bm{y} \sim \mathcal{N}_{kL}(\bm{m}, \bm{\Omega})$, with posterior mean and covariance given by
We obtain a joint posterior for the high-frequency factors $\bm{f}^{(1)}, \ldots, \bm{f}^{(L)}$ where the aggregation equation (ref) results in a natural correlation structure within each low-frequency period (e.g., monthly factors are naturally correlated to predict a quarterly target).
We now present how the EM algorithm for PTFA can be adjusted to handle data observed at a mixed-frequency. Let $T$ denote the amount of lower-frequency observations available, such that there is a total of $\Bar{T} = L \cdot T$ high-frequency observations (e.g., $T$ quarters and $3T$ months of data). Collect all feature observations into an $T \times (pL)$ matrix $\bm{X}$ (e.g., each row has all monthly features associated to each quarter) and all targets into a $T \times q$ matrix $\bm{Y}$. Given the new posterior of the high-frequency factors, the E-step requires the construction of the expected likelihood over $\theta$. Collect the mean vectors $\bm{m}_1, \ldots, \bm{m}_T$ associated to (ref) into a $T \times (kL)$ matrix given by
Letting $\bm{V} \coloneqq \mathbb{E}_{F \mid Y, X; \theta}[\bm{F}^{\top} \bm{F}] = T \cdot \bm{\Omega}_\text{MF} + \bm{M}^{\top} \bm{M}$, this means the posterior expectation of the log-likelihood now takes the form
The M-step then simplifies to obtaining update rules for all components of $\theta$. As the features are now aggregated to a different scale compared to the target, the update steps for $\mathbf{P}$ and $\mathbf{Q}$ are no longer simplified if stacked using matrix operations. Therefore, we present updating steps whose computation will remain efficient even if computed separately.
To this end, write $\bm{X} = [\bm{X}^{(1)}, \ldots, \bm{X}^{(L)}]$ and $\bm{M} = [\bm{M}^{(1)}, \ldots, \bm{M}^{(L)}]$, where each $\bm{X}^{(\ell)}$ block is a $T \times p$ matrix and the $\bm{M}^{(\ell)}$ block is a $T \times k$ matrix, respectively for each $\ell \in \{1, \ldots, L\}$. Similarly, let each $k \times k$ block of $\bm{V}$ be denoted as $\bm{V}_{\ell, r}$ for $\ell, r \in \{1, \ldots, L\}$. The update rules for the loadings for the features and targets under a mixed-frequency setting can then be expressed as
As before, the first-order conditions for these updates allow us to obtain particularly simple and computationally efficient updates for $\sigma_x^2$ and $\sigma_y^2$. Additional details are presented in Appendix (ref).
Finally, we note that all previous derivations can be adapted to the case when $L$ itself changes with time, such that there are $L_t$ high-frequency observations per low-frequency period, which is also a common occurrence in practice. For example, in macroeconomic now-casting of monthly targets such as inflation and industrial production, due to lags and complex interaction between release schedules of useful high-frequency predictors. In finance, differences in trading cycle definitions and firm-specific factors can also cause high-frequency information to be available at differing lengths. By defining a sequence $(L_1, \ldots, L_T)$ of information availability at each time $t$, we can use summation notation instead of matrix operations to efficiently compute update rules without modifying the core derivations.
When working with economic or financial data, it is often unrealistic to assume that the volatility of the data is constant, which becomes a source of misspecification. Next, we show how to allow for stochastic volatility in the context of PTFA and the necessary changes to the EM algorithm to do so. Error covariance matrices in the context of multivariate time series models are usually modeled using multivariate stochastic volatility models, introducing significant computational costs primiceri_time_2005. However, note that in our model the Gaussian noise terms are assumed isotropic, thus depending on a single constant parameter. The computational burden can be further simplified by considering recursive, simulation-free variance matrix discounting methods as in quintana_time_1988. For $\sigma_x$ and $\sigma_y$ we use Exponential Weighted Moving Average (EWMA) estimators. These depend on decay factors $\lambda_x$ and $\lambda_y$ as follows
where $\widehat{\sigma}_{x}^2(t)$ and $\widehat{\sigma}_{y}^2(t)$ are the per-period estimates obtained from our model. In practice, the decay factors $\lambda_x$ and $\lambda_y$ are set to values close to 1, placing more weight on past volatility estimates, thereby making the process smoother and ensuring progressive learning from new data. Setting $\lambda_x = \lambda_y = 0$ mutes stochastic volatility completely and can be made equivalent to the static case by choosing the final estimate of the variances as the time-averages of $\sigma_{x}^2(t)$ and $\sigma_{y}^2(t)$.
These EWMA processes allow us to dynamically adjust the volatility estimates as the model iterates through time, capturing the time-varying nature of the volatility in both features and targets. The estimated volatilities, $\sigma_{x}(t)$ and $\sigma_{y}(t)$, are then used in the next iteration of the model, ensuring volatility is incorporated into parameter estimation. This iterative procedure ensures that the volatilities evolve over time, reflecting the dynamic nature of the system. We note that beyond being computationally trivial and very flexible, the EWMA provides an accurate approximation to an integrated GARCH process. The full EM implementation with this extension can be found in Algorithm (ref) of Appendix (ref).
Our final key extension to probabilistic targeted factor analysis is the introduction of dynamics in the factor process. In economic and financial applications, it is common to interpret the factors as common aggregate effects that drive co-movements between features and targets. It is then natural to assume that these aggregate effects show persistence, such that it is worthwhile to acknowledge dynamic relationships between the recovered factors from the features used to predict targets. These kind of model dynamics are collectively referred to as dynamic factor models stock_dynamic_2011.
We assume the following simple dynamic process where the factors are linearly related to their value in their previous period through a vector autoregressive (VAR) structure. Given initial condition $\bm{f}_0$, this results in the following modified model equations:
where we explicitly add a time index $t$ to emphasize the dynamics, $\mathbf{A}$ is a $k \times k$ matrix of coefficients, and $\bm{v}_t$ is assumed multivariate Gaussian white noise with covariance matrix $\mathbf{V}_F$ and uncorrelated to $\bm{e}_{x, t}$ and $\bm{e}_{y, t}$. Due to well-known identification restrictions forni_generalized_2000, stock_dynamic_2011, doz_dynamic_2020, the model defined by (ref)--(ref) is without loss of generality, as one cannot disentangle lagged factor effects on features and targets from higher-order dynamics (i.e., more lags) in the factor VAR.
Identification restrictions also prevent us from recovering the contemporaneous correlation structure of the factors unless one imposes structural constraints, similar to those required for Structural VAR analysis. Therefore, we do not pursue estimation of the variance of the errors in the dynamic factor equation and instead we let $\mathbf{V}_F$ represent a pre-specified prior variance as before, leaving $\mathbf{A}$ unrestricted. The $k$-dimensional initial condition $\bm{f}_0$ and coefficient matrix $\mathbf{A}$ then become additional parameters to estimate in this setting, such that we augment the model parameters to $\theta = (\mathbf{P}, \mathbf{Q}, \sigma_x^2, \sigma_y^2, \mathbf{A}, \bm{f}_0)$. As the conditional likelihood at each time $t$ does not change compared to the static model (ref), we focus on providing the modified EM update rules.
The autoregressive law of motion for the factors (ref) implies a joint conditional prior $p(\bm{f}_1, \ldots, \bm{f}_T \mid \bm{f}_0; \mathbf{A}, \mathbf{V}_F)$ that critically does not depend on the loadings or error variances. This leads to the following expression for the E-step likelihood
where we again collect all factors into a $T \times k$ matrix $\mathbf{F}$ and denote $\mathbf{M} \coloneqq \mathbb{E}_{\mathbf{F}|\mathbf{X},\mathbf{Y}; \theta}[\mathbf{F}]$ as well as $\mathbf{V} \coloneqq \mathbb{E}_{\mathbf{F}|\mathbf{X},\mathbf{Y}; \theta}[\mathbf{F}^{\top} \mathbf{F}] = \sum_{t=1}^{T} \mathbb{E}_{\mathbf{F}|\mathbf{X},\mathbf{Y}; \theta}[\bm{f}_t \bm{f}_t^{\top}]$. While these expectations had simple analytical expressions when assuming no factor dynamics (see Eqs. (ref) and (ref)), these are not available in the current scenario. Fortunately, there is a large body of literature that considers computationally efficient techniques to update these expectations across time, such as the Kalman filter welch_introduction_1995. As the parameters $\mathbf{A}$ and initial condition $\bm{f}_0$ defining the dynamics do not enter the conditional likelihood, the M-step update rules for the loadings are the same as in the static case, such that $\mathbf{L} = \mathbf{Z}^{\top} \mathbf{M} \mathbf{V}^{-1}$ regardless of the method used to compute the expectations $\mathbf{M}$ and $\mathbf{V}$, where we again stack loadings and data into $\mathbf{L}$ and $\mathbf{Z}$, respectively. The update rules for the variances are also the same as in (ref).
In the spirit of computational simplicity and efficiency of the paper, instead of the aforementioned Kalman filter we present the factor posterior and update rules for the remaining parameters using a banded matrix approach chan_bayesian_2019. Define the $TK$-dimensional prior mean vector as $\bm{\mu}_0$ and let $\mathbf{H}_A$ be a $Tk \times Tk$ block banded matrix, such that they can be written as
As $\mathbf{H}_A$ is lower-triangular also, its determinant equals the product of the determinants of the diagonal blocks, such that $\det (\mathbf{H}_A) = 1$ and it is therefore invertible. This allows us to express (ref) succinctly as $\operatorname{vec}(\mathbf{F}^{\top}) = \bm{\mu}_0 + \mathbf{H}_A^{-1} \operatorname{vec} (\mathbf{V}^{\top})$. The factor posterior in vectorized form is then $\operatorname{vec}(\mathbf{F}^{\top}) \mid \mathbf{X}, \mathbf{Y}; \theta \sim \mathcal{N}_{Tk}(\operatorname{vec}(\mathbf{M}^{\top}), \bm{\Omega}_\text{DFM})$, where the posterior mean and covariance are given by
Finally, define $\mathbf{V}_{s, t} \coloneqq \bm{\Omega}_{s, t} + \bm{m}_s \bm{m}_t^{\top}$ using the $st$-th $k \times k$ block of $\bm{\Omega}_\text{DFM}$ and the $t$-th $k$-dimensional row of $\mathbf{M}$, for $s, t \in \{1, \ldots, T\}$. We can then write the M-step update rules for the remaining parameters as
Computing the elements $\bm{\Omega}_{s, t}$ can be challenging, as even if $\bm{\Omega}^{-1}$ is a banded matrix, its inverse will be dense. However, note that we only ever require to compute the sums of elements along the bands of $\bm{\Omega}$. As the banded matrix structure of $\bm{\Omega}^{-1}$ is preserved in its Cholesky decomposition, one can obtain efficient algorithms to recover the band elements of the full inverse. For our EM method, we draw on the simple and general implementations of these operations provided by durrande_banded_2019. The full steps of pseudo-code can be found in Algorithm (ref) of Appendix (ref). Their derivative operations could also potentially be used for inference in our setting, though we leave this exploration for further research as it is outside the scope of the current paper.
Note that while we assumed a vector autoregressive (VAR) process with only one lag for the factors in (ref), we can easily extend the previous formulas to the case of more complex VAR$(\ell)$ dynamics by simply augmenting $\mathbf{H}_A$. Alternatively, one can use the fact that the VAR companion form implies any dynamical system of $n$ variables with $\ell$ lags can be expressed as a VAR$(1)$ system on a system with $n(\ell-1)$ variables. By selecting a larger number of components $k$, we could then also control for the presence of more persistent co-movements without needing to model these dynamics, at the expense of larger calculations. For this reason, one expects the number of static factors required to capture co-movements in time series to be larger than the number of dynamic factors required. Choosing the number of factors $k$ in supervised estimation methods for fitting and optimal forecasting performance remains challenging ahn_forecasting_2022. In our implementations and empirical application we either explore a range of values of $k$ and note their performance or obtain $k$ by cross-validation (CV).
In this section, we present several simulation exercises conducted to evaluate the performance of PTFA compared to traditional factor extraction techniques. We maintain PLS as the baseline for targeted decomposition, and include both standard PCA and PPCA as competing unsupervised decomposition techniques.\footnote{As the output from PCA and PPCA are non-targeted factors that maximize the explained variance of predictors $\mathbf{X}$; we fit a linear regression of $\mathbf{Y}$ on these factors to obtain the final predictions.} The goal is to assess the in-sample accuracy in predicting the targets under alternative data-generating processes for the noise in predictor and response variables. In the following section we then explore the out-of-sample performance of PTFA in three real-world applications to time series data. We focus the discussion of the simulation comparisons around PTFA and PLS as supervised decomposition methods, leaving the unsupervised PCA and PPCA results to Appendix (ref).
The first step of the simulation revolves around the factor structure of predictors and targets before adding noise. Throughout all exercises, we set $T = 200, p = 10, q = 3$, and $k = 2$. We first draw all entries in the loadings $\mathbf{P}$ and $\mathbf{Q}$ from a uniform distribution between 0 and 1. Then, we generate $\bm{f}_t \sim \mathcal{N}_k(\bm{0}_k, \mathbf{I}_k)$ independently for each time period $t \in \{1, \ldots, T\}$. Given the factors and loadings, we finally generate features and targets according to equations (ref) and (ref), respectively.
The main differences across each of the DGPs is the distribution of the errors $e_x$ and $e_y$. For the simplest (and correctly specified) DGP, we consider isotropic Gaussian errors
Figure (ref) summarizes the key finding on a single realization of simulated data from (ref), where we fix $\sigma_x = \sigma_y = 1$. Note how the predicted targets from PTFA align more closely to the true targets when compared to standard PLS. The $R^2$ score value for each is also higher resulting in an average score of 68.1% for PTFA compared to 56.5% in standard PLS on this single realization. Figure (ref) in the appendix presents the path taken by the values of $R^2$ of the fit as the EM iterations of our algorithm progress. Notice how the algorithm quickly adapts to a large level of explained variance in the targets and levels off once the estimates reach numerical convergence as measured by the $\ell_2$ distance between iterates.
As discussed in the main text and formally shown in Appendix (ref), PTFA uses the assumption of isotropic Gaussian noise terms simply to provide a probabilistic framework to targeted factor extraction. Performance of PTFA should therefore not depend on whether feature ($\mathbf{X}$) and target ($\mathbf{Y}$) variables are correlated or even normally distributed. Through this and the next simulation exercises, we show that the relative performance between PTFA and PLS does not depend on the assumed distribution of the variables being decomposed.
As a first extension, we dispose of the isotropic assumption and allow for the noises to be multivariate normal distributions with non-diagonal covariance matrices. For this example, we assume the following Toeplitz covariance structure for both features and targets:
such that $\bm{e}_x \sim \mathcal{N}(\bm{0}_p, \bm{\Omega}_x)$ and $\bm{e}_y \sim \mathcal{N}(\bm{0}_q, \bm{\Omega}_y)$, where $\rho_x, \rho_y \in [-1, 1]$ are correlation parameters. Figures (ref) and (ref) present the same statistics as before for a realization of data from (ref) using $\rho_x = \rho_y = 0.5$ (keeping the remaining values the same as in the previous exercise). Similar results as in the isotropic Gaussian case are obtained, with PTFA dominating PLS in terms of in-sample fit within only a small number of iterations of the EM algorithm.
We present evidence that similar results hold once we dispense the assumption of Gaussian noise altogether. Specifically, we consider the following setup design to produce heavy-tailed and asymmetric noise that results in clear deviations from Gaussian features and targets. Errors in features are drawn independently from a Student-$t$ distributions with 3 degrees-of-freedom and scale $\sigma_x$, while target noise is drawn from a $\chi^2$ distribution with 1 degree of freedom.
Figures (ref) and (ref) showcase that similar results to before arise from specification (ref). As long as the data is standardized prior to processing, it can be observed that PTFA will deliver targeted factors that are generally more accurate to summarize the information in the targets regardless of the distributions of the variables involved.
Crucially, Figure (ref) showcases that these performance gains do not depend on any given realization of data. By comparing the average (across targets) $R^2$ statistics over 1000 replications of the previous setting, we find that PTFA first-order stochastic-dominates PLS in generating better in-sample fit. PLS itself dominates the other two unsupervised techniques PCA and PPCA, while they are indistinguishable from each other under the correctly specified (ref).
Finally, given the critical role of noise in explaining the virtue of PLS, we additionally simulate data with differing levels of noise in both features $\mathbf{X}$ and targets $\mathbf{Y}$. That is, once again we simulate noisy data from (ref), adjusting both error scales $\sigma_x$ and $\sigma_y$ over a grid between 0.1 and 5, covering a wide range of signal-to-noise ratios.
Figure (ref) shows the median value of average $R^2$ statistics across targets over 1000 replications of this simulation setup. The superior performance of PTFA is evident from the heatmap. The gains in terms of goodness-of-fit of PTFA when compared to PLS are more salient when noise is increases, in particular when the noise is in the targets instead of the features. This is as seen in the PCA case, where perturbations to the data in the form of noise or outliers creates issues for consistently recovering the axes of maximal variance chen_spectral_2021.
Financial conditions indices (FCIs) offer a single, quantitative measure summarizing the state of financial conditions in the economy by combining information from a wide range of variables --- including credit spreads, asset prices, liquidity measures, and volatility indicators. These indices have become essential tools for policymakers, especially since the Global Financial Crisis, as they help to gauge the overall tightness or looseness of financial conditions and, in turn, provide early signals of potential risks to the real economy adrian_vulnerable_2019.
Traditional approaches to constructing FCIs typically rely on principal component analysis (PCA) or its variants applied to balanced panels of data hatzius_financial_2010, brave_monitoring_2011, koop_new_2014 and extract one single factor ($k = 1$) to match the idea of an index. However, high-frequency financial datasets are often incomplete due to mismatched reporting frequencies, data lags, or sporadic data availability. This incompleteness can lead to noisy and less reliable indices when using standard techniques. For instance, least-squares PCA methods may overfit or produce overly volatile indices in the presence of substantial missing data. Nevertheless, the main critique to current FCIs is related to their identification and interpretability. Since FCIs are constructed as unsupervised common components of a large number of financial time-series, their signal and relevance to future economic developments are unclear and hard to learn. Poor identification of the single factor that proxies an FCI has led some authors to extract two factors to fully characterize the state of financial conditions in the economy Lombardi2025.
To address these challenges, we propose constructing a targeted FCI. Unlike PCA-based FCIs, which summarize variation in financial variables without regard to a macro target, our targeted FCI extracts the component of financial conditions that is most informative about a chosen target $\textbf{Y}$. In essence, the targeted FCI aims to provide a cleaner and more timely measure of financial conditions by filtering out information irrelevant to the target of interest. This directly addresses the identification and interpretability critique of traditional FCIs. For instance, if we define $\textbf{Y}$ to be GDP growth, we are effectively looking for a FCI relevant in explaining economic growth. Whereas, if we target inflation, we would effectively be constructing a FCI that could signal inflation risk through time.
Therefore, our approach allows the econometrician to construct as many different FCIs as the number of targets they may wish to define, from the same pool of financial variables. This approach can significantly enhance the usefulness of the index for real-time monitoring and forecasting, ultimately supporting better-informed monetary policy decisions. To implement such an approach, we employ our Targeted Dynamic Factor Model (PTFA-DFM), which allows for mixed-frequency data and a more robust treatment of missing data in general, while enabling the weights to depend on the target variable.
We begin by collecting a large number of financial time-series, reported at monthly frequency as summarized below in Table (ref). The first three variables listed are candidate targets, which can be studied in isolation or jointly as multivariate objects. The $p = 16$ variables outlined in rows 4 to 19 are used to define $\textbf{X}$ and represent a cross-section of proxy financial conditions indicators such as the VIX which captures sentiment in financial markets, the return of the S&P 500 index, a series of credit spreads, which have been found to contain relevant information to predict real economic development, among others.
Next, we fit our targeted Dynamic Factor Model to this data to examine the host of potential FCIs one can construct. Figure (ref) plots different Targeted FCIs, which differ with respect to the variables they target. The main result is that the target choice matters: the GDP and unemployment targeted FCIs are the closest to the PCA counterpart and the most volatile, reacting strongly around growth slowdowns, while the inflation-targeted FCI is more muted and its dynamics are more pronounced in later stages of recessions rather than at mature stages of the business cycle.
An interesting result emerges from the inspection of the loadings of each financial variables onto the alternative FCIs we consider. Figure (ref) describes the most influential variables per FCI. We can observe that credit-related variables such as the treasury bill spread, debt-to-GDP ratio, or the mortgage spread load more heavily on the inflation targeted FCI. In contrast, market volatility (VIX) and equity returns (S&P500) weigh more heavily in the GDP and unemployment targeted indices. It is apparent that Financial Conditions Indices are sensitive to their targets and their dynamics vary accordingly.
Partial Least Squares may be valuable when the key motivation for dimensionality reduction is prediction. To demonstrate the practical relevance of PTFA, we conduct a simple macroeconomic forecasting exercise using the Federal Reserve Economic Data -- Monthly Database mccracken_fred-md_2016. The data includes $p = 126$ variables that track macroeconomic developments in the United States at a monthly frequency. In this section, we discuss the main results and compare PTFA with popular alternatives in the literature.
Table (ref) shows out-of-sample (OOS) Mean Squared Forecast Error (MSFE) of prediction for each target variable using $k = 7$ factors extracted with either PLS, PTFA, PTFA-DFM, PTFA-SV, PPCA or PCA. The targets in our exercise are Industrial Production, CPI inflation and the Unemployment Rate (such that $q=1$ for each panel). The number of factors $k$ is chosen to be consistent with the FRED-MD factors calculated according to mccracken_fred-md_2016, with code made available by the authors. The key message from Table (ref) is that the PTFA outperforms both PLS and PCA, adding value to forecasts across forecast horizons and macroeconomic variables considered. PTFA delivers systematic OOS gains by extracting factors that maximize predictive covariance with targets, down-weighting directions in predictors that inflate variance but do not aid prediction. This mechanism mirrors our Monte Carlo exercise, discussed in previous section, where improvements intensify as the noise in targets rise. We note that allowing for stochastic volatility (PTFA-SV) does not materially improve macro forecasts over the baseline PTFA except for the case of CPI inflation for $h=1$, suggesting the gains stem primarily from targeted factor extraction rather than time-varying error dynamics. We also note that the Dynamic Factor extension of PTFA performs quite well, outperforming competitors for many combinations of forecast horizons $h$ and target variables.
Attempts to predict stock returns or the equity premium are in no short supply in the Finance literature. welch_comprehensive_2007 and goyal_comprehensive_2024, provide a review and comprehensive assessment of the performance of 46 different variables that have been suggested by the academic literature to be good predictors of the equity premium. Following this large body of empirical work, our financial application studies the predictability of U.S. aggregate stock returns, using the goyal_comprehensive_2024 dataset.
The goal of this exercise is to predict the equity risk premia, using $p = 26$ signals that are available at a monthly frequency. Thus, in this case $\mathbf{Y}$ is the equity premium and $\mathbf{X}$ are the various predictors, lagged by one period, following standard practice. In this setting, it is less clear how many factors $k$ should be considered. Therefore, we estimate PTFA and all competing models for different values of $k$ and report the forecasting performance results in Table (ref). As for the forecasting horizon, we only consider 1 month ahead forecasts given the nature of the problem of forecasting stock returns, since information is priced-in quite fast.
The main message from Table (ref) is that PTFA and extensions add value as compared to PLS and PCA in predicting the equity risk premia. Similar to our application with macroeconomic data, our model with stochastic volatility (PTFA-SV) does not seem to outperform PTFA, that is relatively more parsimonious. We observe that MSFE associated to PTFA forecasts are quite competitive, regardless of how many factors one chooses. PTFA-DFM delivers uniformly lower MSFEs than PCA and PLS across different choices of $k$ with the tightest gains for $k=1$. As $k$ increases, OOS MSFE rise for all models, reflecting higher variance from estimating weakly identifiable directions in high-noise, short-sample settings. The one-factor PTFA solution concentrates on the most stable predictor–target covariance direction and avoids the overfitting that afflicts multi-component PLS, a clear advantage of PTFA.
We introduce a probabilistic framework for targeted factor extraction called PTFA and derive a fast Expectation-Maximization (EM) algorithm to estimate the model. PTFA is flexible and naturally handles parameter uncertainty, noise, and missing data in estimation. Through simulation exercises and three real-world applications in macroeconomic forecasting, equity premium prediction, and the construction of financial conditions indices, we demonstrate the superior performance of PTFA, especially in noisy and incomplete data environments. Along the way, we provide additional contributions to mixed-frequency data, stochastic volatility, time-series persistence in the latent factors, and give further theoretical insight to the probabilistic PLS solutions.
Our probabilistic foundation also opens many avenues for future research, including interesting methodological extensions using probabilistic (fully Bayesian) or variational inference. By providing an open-source implementation of the method, our hope is that practitioners of time-series forecasting and researchers alike will continue to expand and improve upon the technique.
\begingroup \setstretch{0.85} \endgroup