Multi-Level Dynamic Factor Models

• Multi-Level Dynamic Factor Models are dynamic latent-variable models that incorporate a hierarchical loading structure to separate pervasive global factors from group-specific local factors.
• They utilize sequential least squares and state-space approaches to extract and orthogonalize factors, ensuring precise identification and efficient estimation.
• These models are applied in macro-financial scenario design, micro–macro integration, and decomposing common versus idiosyncratic movements in large, structured datasets.

Multi-Level Dynamic Factor Models (ML-DFMs) are dynamic latent-variable models that extend the standard dynamic factor model by imposing a hierarchical loading structure on a large panel of observed series. In the canonical formulation, a small set of pervasive factors affects all variables, while additional group-specific factors affect only variables inside their own block; some formulations also add semipervasive factors shared by subsets of blocks. Across recent work, ML-DFMs appear in at least three closely related forms: block-structured approximate factor models estimated by sequential least squares, multidimensional state-space systems linking macro and micro data, and matrix-structured models with global and local factors. At the same time, several nearby models—deep nonlinear DFMs, structured matrix DFMs, and static multilevel covariance models—are methodologically relevant but are not, strictly speaking, ML-DFMs (Bellocca et al., 16 Feb 2026, Bellocca et al., 14 Jul 2025, Barigozzi et al., 2023, Zhang et al., 2023).

1. Conceptual definition and scope

The basic ML-DFM enriches the standard approximate factor model by separating variation into common components that are pervasive across the full cross section and local components that are restricted to known groups. With groups indexed by $s=1,\dots,S$, variables indexed by $i=1,\dots,N_s$, and time indexed by $t$, one canonical observation equation is

$$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$

where $\mathbf{G}_t$ collects global factors, $\mathbf{L}_{s,t}$ collects group-specific factors, and the local loading blocks are zero outside their own group. In stacked form, the loading matrix is block-structured, with global loadings stacked across groups and local loadings arranged block diagonally (Bellocca et al., 16 Feb 2026).

A related block formulation partitions the observed vector as $X_t=(X_{1,t},\dots,X_{K,t})'$ and writes

$X_t=P^{*}F_t^{*}+\epsilon_t,$

where $F_t^{*}$ can include pervasive factors $G_t$, block-specific factors $i=1,\dots,N_s$0, and, in overlapping-block designs, pairwise semipervasive factors such as $i=1,\dots,N_s$1. This extension is useful when some latent sources are shared only across selected blocks rather than the entire panel (Bellocca et al., 14 Jul 2025).

This suggests that “multi-level” refers first to the loading architecture rather than to any single estimation routine. In the papers considered here, the hierarchy may be non-overlapping, overlapping, additive, or embedded in a state-space system, but the defining feature is the imposition of structured zero restrictions or equality constraints that distinguish common and group-level information from idiosyncratic noise (Bellocca et al., 16 Feb 2026, Bellocca et al., 14 Jul 2025, Barigozzi et al., 2023).

2. Core mathematical structure and identification

The compact stacked ML-DFM representation is

$i=1,\dots,N_s$2

with $i=1,\dots,N_s$3. The decisive structural feature is the block zero restriction: local factors load only within their own group (Bellocca et al., 16 Feb 2026).

In non-overlapping block systems, the measurement equation can be written explicitly as

$i=1,\dots,N_s$4

while overlapping-block models augment the latent vector with pairwise semipervasive factors. The hierarchy is therefore not merely “more factors”; it is a factorization with economically meaningful loading restrictions (Bellocca et al., 14 Jul 2025).

Identification follows the usual factor-model logic, but admissible rotations must preserve the zero restrictions. One set of restrictions normalizes local factors and global factors separately: $i=1,\dots,N_s$5

$i=1,\dots,N_s$6

and imposes orthogonality between global and local factors,

$i=1,\dots,N_s$7

A notable technical point is that group-specific factors from different groups need not be orthogonal for identification, although enough normalization and orthogonality is imposed to place estimated and true factors in the same identified space (Bellocca et al., 16 Feb 2026).

In sequential least-squares implementations without an explicit normalization, only the common component $i=1,\dots,N_s$8 is initially just-identified; the separate factors and loadings are identified only up to rotation. The normalization is then recovered through orthogonalization across levels and principal-components normalization of estimated common components, with the orthogonalization implemented recursively via regressions equivalent to Gram-Schmidt (Bellocca et al., 14 Jul 2025).

3. Estimation, factor extraction, and inference

A central estimation route is the sequential least squares (SLS) procedure associated with Breitung and Eickmeier. In the formulation studied for ML-DFM inference, estimation proceeds by: extracting $i=1,\dots,N_s$9 factors by PC within each group, using Canonical Correlation Analysis (CCA) across groups to estimate common global combinations, regressing the data on initial global factors to obtain residuals, extracting local factors within each group, then alternating least-squares updates for loadings and factors until the residual sum of squares converges. After convergence, the estimated factors are renormalized so that the factor covariance matrices equal identity and the loading cross-products are diagonal (Bellocca et al., 16 Feb 2026).

The FARS implementation generalizes the same sequential logic to non-overlapping and overlapping blocks. Initial global factors can be obtained by CCA or PCA; global components are filtered out block by block; lower-level factors are then extracted sequentially; block-specific factors are finally obtained by PC on residualized blocks; and the algorithm iterates between factor and loading updates until the residual sum of squares converges. The package exposes this structure through mldfm(), with method = 0 for CCA initialization and method = 1 for PCA initialization (Bellocca et al., 14 Jul 2025).

For inference, the key recent question is whether Bai’s asymptotic distribution for PC factors in standard DFMs can approximate the empirical distribution of SLS-estimated global and group-specific ML-DFM factors. The benchmark law is

$t$0

with

$t$1

Monte Carlo evidence shows that this approximation works well for moderate to large samples for both global and group-specific SLS factors, especially for the global factor and for local factors with stronger loadings. When local loadings are relatively weak, the asymptotic covariance can approximate empirical covariance better than full empirical mean squared error because finite-sample bias remains non-negligible (Bellocca et al., 16 Feb 2026).

The practical covariance estimator is

$t$2

Two alternatives for $t$3 are emphasized. The heteroscedasticity-robust Bai–Ng style estimator assumes no cross-sectional idiosyncratic correlation, whereas the FPR estimator uses a thresholded estimate of the idiosyncratic covariance matrix and therefore allows cross-sectional dependence. Across the simulations, the preferred uncertainty estimator is FPR with subsampling correction for loading-estimation uncertainty; when cross-sectional dependence is absent, HR and FPR perform similarly, but when dependence is present, FPR clearly dominates HR (Bellocca et al., 16 Feb 2026).

The same inferential logic underlies factor confidence regions in FARS. The corrected MSE combines the plug-in asymptotic term with a cross-sectional subsampling term,

$t$4

and yields confidence ellipsoids

$t$5

These ellipsoids are then used for stress design (Bellocca et al., 14 Jul 2025).

4. State-space, multidimensional, and matrix-structured generalizations

A major extension replaces a fixed $t$6 panel with multidimensional dependent data whose cross section changes over time. In the multidimensional dynamic factor model, the observed object at time $t$7 is a matrix $t$8, which is vectorized and embedded into a fixed-dimension series $t$9, $$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$0, with missingness treated as intrinsic rather than as zero padding. The generic model is

$$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$1

$$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$2

which yields a state-space ML-DFM capable of processing repeated cross-sections, missing values, and asynchronous releases (Barigozzi et al., 2023).

In the household–macro specialization, the observed vector is

$$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$3

with $$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$4 the vector of 8 macroeconomic aggregates and $$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$5 the vector of real income per head for households in demographic group $$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$6. The decomposition combines smooth trends $$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$7, a common business-cycle factor $$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$8, group- and series-specific idiosyncratic latent cycles $$y_{s,it}=\boldsymbol{\lambda}_{s,i}^{(g)\prime} \mathbf{G}_t+\boldsymbol{\lambda}_{s,i}^{(l)\prime} \mathbf{L}_{s,t}+\boldsymbol{\varepsilon}_{s,it},$$9, and measurement noise $\mathbf{G}_t$0. For micro incomes, all households in the same demographic group share the same trend and business-cycle loading through replication by $\mathbf{G}_t$1, yielding a hierarchical additive dynamic factor structure rather than a tensor factor model (Barigozzi et al., 2023).

Estimation in that framework uses penalized quasi maximum likelihood with an ECM algorithm built around a Kalman smoother. The E-step accommodates arbitrary missingness through time-varying observed sets $\mathbf{G}_t$2 and observation matrices $\mathbf{G}_t$3; the CM-step updates initial conditions, transition coefficients, innovation variances, and loadings under structural equality constraints. The model can therefore process rotating survey participation, ragged edges, and publication lags while remaining operational for real-time updating, localized predictions, counterfactuals, and impulse response functions (Barigozzi et al., 2023).

Matrix-valued observations lead to a different but related generalization. The multilevel matrix factor model observes $\mathbf{G}_t$4 and writes

$\mathbf{G}_t$5

where $\mathbf{G}_t$6 is global and $\mathbf{G}_t$7 is local to group $\mathbf{G}_t$8. In vectorized form,

$\mathbf{G}_t$9

so the model is a multilevel factor model with Kronecker-structured loadings. Estimation is sequential: identify global row and column loading spaces from cross-group covariance, project them out, estimate local loading spaces from residual matrix dynamics, then recover global and local signal parts. The identified objects are loading spaces rather than unique loading matrices, reflecting the usual rotational indeterminacy (Zhang et al., 2023).

5. Methodological boundaries and neighboring models

Several recent models are closely related to ML-DFMs but should not be conflated with them. The most explicit boundary case is the deep dynamic factor model $\mathbf{L}_{s,t}$0, which replaces the linear factor-observation map with a deep autoencoder while preserving dynamic latent states and idiosyncratic AR components. Its core equations are

$\mathbf{L}_{s,t}$1

with AR laws for factors and idiosyncratic terms. However, it has one shared latent bottleneck $\mathbf{L}_{s,t}$2, no grouped or nested decomposition, no global versus local blocks, and no identification of global versus local shocks. Its relevance to ML-DFMs is therefore methodological rather than structural (Andreini et al., 2020).

Dynamic matrix factor models are another neighboring class. They model matrix-valued time series through

$\mathbf{L}_{s,t}$3

or, after vectorization, through Kronecker-structured loadings and transitions. This yields a structured dynamic factor model with row and column loading spaces, but not an explicit hierarchy of global and group-specific factors. The structure is two-way rather than multilevel in the standard nested sense (Yu et al., 2024).

Parshakova, Hastie, and Boyd study a static Gaussian multilevel factor model with covariance $\mathbf{L}_{s,t}$4 constrained to have multilevel low-rank structure. The model is

$\mathbf{L}_{s,t}$5

with no time index and no state equation. Its main contribution is computational: a fast EM algorithm, recursive Sherman–Morrison–Woodbury inversion, and Cholesky/determinant machinery with linear time and storage complexities per iteration under fixed hierarchy and rank allocation. This is directly relevant to ML-DFM computation on the observation side, but it is not itself dynamic (Parshakova et al., 2024).

A further source of confusion is the phrase “level” in models for means and volatilities. In the joint level–volatility dynamic factor model,

$\mathbf{L}_{s,t}$6

the baseline contribution is a joint dynamic factor model for first moments and second moments, not a hierarchical ML-DFM over blocks or regions. The international inflation application does introduce a global plus regional factor decomposition for both means and volatilities, making the connection to ML-DFMs direct there, but “level” in the title refers to the conditional mean rather than to hierarchical level (Mumtaz et al., 4 Apr 2026).

These distinctions matter because the recent literature uses related language—multidimensional, multilevel, deep, matrix, level, volatility—for models that solve different problems. A common misconception is therefore to treat any structured latent dynamic model as an ML-DFM. The papers surveyed here instead support a narrower usage: ML-DFMs require an explicit hierarchical or grouped factor architecture, not merely nonlinear encoding, matrix geometry, or multilevel covariance structure (Andreini et al., 2020, Yu et al., 2024, Parshakova et al., 2024, Mumtaz et al., 4 Apr 2026).

6. Empirical roles, scenario design, and substantive applications

One applied role of ML-DFMs is macro-financial scenario design. FARS uses ML-DFM factors as inputs to factor-augmented quantile regressions,

$\mathbf{L}_{s,t}$7

then recovers predictive densities by fitting a skew-$\mathbf{L}_{s,t}$8 distribution to the implied quantiles. Stress scenarios are defined by optimizing a target quantile over the boundary of the factor confidence ellipsoid; the stressed factors are then propagated through the quantile regression and density-recovery steps. In the U.S. GDP illustration with three blocks—63 global macro variables, 248 domestic macro variables, and 208 global financial variables—the specification uses one global factor, one pairwise factor for blocks 1 and 3, and one local factor for each block. The application reports that GiS is more negative than GaR and that this gap persists across horizons $\mathbf{L}_{s,t}$9 and becomes more adverse at higher stress levels $X_t=(X_{1,t},\dots,X_{K,t})'$0 (Bellocca et al., 14 Jul 2025).

A second role is micro–macro integration. The multidimensional state-space formulation jointly models 8 macro aggregates from ALFRED and quarterly household income from the Consumer Expenditure Public Use Microdata over 1990–2020, with about 87,000 households restricted to prime-age urban households. The empirical hierarchy has four levels in practice: time, macro series, four household groups defined by education and ethnicity, and households within groups. The model finds persistent income differentials by education and ethnicity, common cyclical co-movement across demographics driven by the business cycle, additional group-specific idiosyncratic movements in income, and the ability to nowcast delayed household survey outcomes using timely macro releases (Barigozzi et al., 2023).

A third role is the decomposition of commonality across cross sections and moments. In the international inflation application of the level–volatility model, inflation is decomposed into a world factor, an AE regional factor, and an EMDE regional factor for both levels and volatilities: $X_t=(X_{1,t},\dots,X_{K,t})'$1

$X_t=(X_{1,t},\dots,X_{K,t})'$2

The reported substantive pattern is a dominant global level component in advanced economies and stronger regional and volatility contributions in emerging and developing economies, indicating that a standard mean-only hierarchy can miss important heterogeneity in common uncertainty (Mumtaz et al., 4 Apr 2026).

Across these applications, a consistent theme is that ML-DFMs are valuable when the panel has known structure that a single pervasive factor block would treat too coarsely. This suggests a unifying interpretation: the central object is not merely a low-dimensional latent state, but a low-dimensional latent state with imposed structure that determines which parts of the panel are allowed to move together. In current work, that structure may describe blocks of variables, subsets of blocks, households within demographic groups, or, in matrix formulations, groups of matrix time series with shared global and local dynamics (Bellocca et al., 14 Jul 2025, Barigozzi et al., 2023, Zhang et al., 2023).