Multilevel Functional Time Series Analysis

• Multilevel functional time series are models that jointly analyze dependent functional observations over time within hierarchical frameworks, capturing both within-function and cross-population dynamics.
• Key methodologies include Bayesian state-space formulations, dual-factor decompositions, dynamic FPCA, and frequency-domain approaches to separate common trends from unit-specific variations.
• These approaches enhance forecasting and inference by borrowing strength across groups, ensuring parsimonious representations and robust dynamic updates through advanced estimation techniques.

Multilevel functional time series refers to jointly modeling multiple, dependent functions observed over time within a hierarchical (multilevel) structure. In this setting, each observation is a curve, surface, or other function on a continuous domain, while dependence arises simultaneously from within-function structure, serial evolution, and cross-series or cross-population coupling. The literature includes Bayesian state-space formulations with common functional bases, multilevel decompositions into common and level-specific trends, dual-factor structures for multiple populations, frequency-domain marginalizations that integrate subject-specific spectra, and interpretable additive formulations with bivariate coefficient surfaces. Across these formulations, the central objective is to separate shared dynamic structure from unit-specific heterogeneity while preserving functional inference and forecastability (Kowal et al., 2014).

1. Observation schemes and multilevel structure

A standard multilevel functional time-series setup observes related functions over time. In one formulation, $y_{j,t}(s)$ denotes the $j$-th function at time $t$, observed over a continuous domain $s\in S$, with $j=1,\dots,J$ and $t=1,\dots,T$. The observation grid may vary across series and times, since $S_{j,t}\subset S$ can be irregular. This accommodates continuous-domain inference together with discrete observation (Kowal et al., 2014).

A multiple-population formulation writes $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$, where populations are indexed by $p=1,\dots,P$, time by $t=1,\dots,T$, and the functional domain by $j$0. The panel can be arranged as a $j$1 matrix of functions, with each column a functional time series and each row a cross-sectional snapshot across populations. In this representation, multilevel structure is represented by population-specific basis functions and loadings that capture cross-sectional heterogeneity, coupled with a set of common factors/loadings that model shared dynamics across populations (Tang et al., 2021).

A more explicitly integrative formulation considers $j$2 for subject or group $j$3, with multi-way dependence consisting of within-curve dependence, temporal dependence, and cross-subject interactions. When subjects are grouped, dependencies arise both within groups and across groups. The stated motivation for an integrative approach is to borrow strength across subjects to improve statistical efficiency, capture shared dynamic structure that drives multiple series, and produce parsimonious representations that enable robust downstream modeling and forecasting (Guo et al., 24 Mar 2026).

This combination of functional observations, serial dependence, and cross-sectional hierarchy distinguishes multilevel functional time series from single-series functional data analysis and from conventional multivariate time series. A plausible implication is that the modeling problem is not only dimension reduction, but also allocation of variation across shared, group-specific, and idiosyncratic functional subspaces.

2. Principal model architectures

Several model classes recur in the literature. They differ mainly in how they represent shared structure, how they encode temporal dynamics, and whether hierarchy is imposed in the time domain, frequency domain, or regression surface.

Framework	Core representation	Multilevel mechanism
Bayesian MFDLM	$j$4	Common basis, hierarchical priors, cross-series dependence
Dual-factor HDFTS	$j$5	Common front-loadings and population-specific back-loadings
MDKLE / spectral integration	$j$6	Shared frequency-dependent filters across subjects
Common/residual multilevel FPCA	$j$7	Common trend plus maturity- or size-specific residual trend
Additive surface model	$j$8	Cross-unit and cross-level coefficient surfaces

The Bayesian Multivariate Functional Dynamic Linear Model couples a time-invariant functional basis with time-varying coefficients that evolve stochastically: $j$9 In the common-basis variant, $t$0 for all $t$1, which pools information across series and makes coefficients directly comparable across $t$2. Stacking all series allows joint $t$3 and $t$4 with block structures capturing cross-series dependence (Kowal et al., 2014).

The dual-factor model for high-dimensional functional time series decomposes centered curves as

$t$5

where $t$6 is a vector of common functional loadings, $t$7 is an $t$8 matrix-valued factor at time $t$9, and $s\in S$0 is a vector of population-specific functional loadings capturing cross-sectional heterogeneity. The paper also rewrites this in two stages by defining a common functional FTS $s\in S$1, followed by a further reduction of $s\in S$2 into common loadings plus matrix time series (Tang et al., 2021).

A frequency-domain alternative begins from subject-specific spectral density operators and integrates them into a marginal spectral operator,

$s\in S$3

with equal weights $s\in S$4 in the paper. The eigenfunctions of $s\in S$5 define optimal functional filters, and filtering yields structured multilevel score vectors $s\in S$6, one per dynamic component (Guo et al., 24 Mar 2026).

The multilevel FPCA decomposition used for foreign-exchange implied volatility surfaces writes each maturity-specific smile as

$s\in S$7

where $s\in S$8 is the common dynamic trend shared across maturities and $s\in S$9 is the maturity-specific residual trend. An analogous two-component decomposition is used for intraday particle number size distribution,

$j=1,\dots,J$0

with the common trend computed by averaging across sizes and the residual process modeled separately for each size (Shang et al., 2021, Shang et al., 2 Oct 2025).

A different direction keeps the model in the original function space and lets lagged curves from one unit affect another through bivariate coefficient surfaces: $j=1,\dots,J$1 For multilevel settings, the same paper writes

$j=1,\dots,J$2

so that cross-level effects are encoded directly by $j=1,\dots,J$3 (Wang et al., 28 Apr 2025).

3. Basis construction, identifiability, and estimation

The main technical difficulty is that multilevel functional time series require simultaneous control of smoothness, identifiability, and dependence.

In the MFDLM, factor loading curves are represented by cubic B-splines with knots $j=1,\dots,J$4,

$j=1,\dots,J$5

and smoothness is enforced through the roughness penalty

$j=1,\dots,J$6

The Bayesian formulation places

$j=1,\dots,J$7

so constant and linear parts are unpenalized and nonlinear components are shrunk by $j=1,\dots,J$8. Ordering is imposed through $j=1,\dots,J$9, which orders $t=1,\dots,T$0 by decreasing smoothness. Identifiability is enforced by orthonormality,

$t=1,\dots,T$1

or equivalently $t=1,\dots,T$2. Orthogonality is handled by conditioning Gaussian posteriors on linear constraints, and unit norm by post-sample normalization (Kowal et al., 2014).

Dynamic FPCA-based multilevel methods replace static covariance with long-run covariance. For stationary $t=1,\dots,T$3, long-run covariance is defined by

$t=1,\dots,T$4

and estimated by a kernel sandwich estimator with bandwidth $t=1,\dots,T$5 and kernel $t=1,\dots,T$6. The paper on foreign-exchange implied volatility surfaces states that dynamic FPCA properties used there are: it minimizes mean integrated squared reconstruction error over the whole functional sample, extracts long-run covariance, and produces uncorrelated dynamic principal component scores (Shang et al., 2021).

In the dual-factor model, estimation begins with dynamic FPCA on each population’s long-run covariance to obtain $t=1,\dots,T$7, chooses $t=1,\dots,T$8 by cumulative percentage of variance with $t=1,\dots,T$9, solves a functional concurrent regression for $S_{j,t}\subset S$0, and then estimates the common loadings $S_{j,t}\subset S$1 from eigenfunctions of an aggregated nonnegative operator $S_{j,t}\subset S$2. The resulting $S_{j,t}\subset S$3 is the low-dimensional matrix-valued time series that carries most temporal dynamics (Tang et al., 2021).

The frequency-domain approach estimates subject-specific spectra by a Bartlett lag-window estimator, integrates them into $S_{j,t}\subset S$4, computes eigenfunctions $S_{j,t}\subset S$5, aligns phase through $S_{j,t}\subset S$6, and forms time-domain filters

$S_{j,t}\subset S$7

This explicitly addresses the fact that frequency-domain eigenfunctions are determined up to unit-modulus multipliers (Guo et al., 24 Mar 2026).

The trend-and-seasonality framework with variable phase treats each observed function as

$S_{j,t}\subset S$8

with $S_{j,t}\subset S$9. Identifiability is obtained by the Karcher mean constraint on inverse warps, $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$0, and $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$1, $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$2 for orthogonal subspaces $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$3 and $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$4 (Tai et al., 2017).

The additive surface model represents $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$5 with bivariate splines on a triangulation $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$6, imposes $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$7 continuity with a constraint matrix $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$8, adds a roughness penalty $X_{p,t}(u)\in \mathcal H:=\mathcal L^2(\mathcal I)$9, and induces sparsity by a group bridge penalty

$p=1,\dots,P$0

This yields both global removal of entire $p=1,\dots,P$1 surfaces and local removal of specific triangles within a surface (Wang et al., 28 Apr 2025).

4. Dynamic inference and computation

Multilevel functional time-series computation is dominated by latent-state updating, basis or filter estimation, and hyperparameter learning.

The MFDLM uses a Gibbs sampler that partitions parameters into blocks. Spline coefficients are sampled from constrained Gaussian full conditionals; latent states are updated with forward-filtering backward-sampling; evolution matrices $p=1,\dots,P$2 can be sampled under Gaussian priors; and variance components such as $p=1,\dots,P$3 and $p=1,\dots,P$4 use conjugate Gamma and Wishart updates. The paper states that spline updates scale with $p=1,\dots,P$5 per iteration, while FFBS scales with $p=1,\dots,P$6 in the worst case, though block-diagonal structures and diagonal $p=1,\dots,P$7 greatly reduce this burden. Common-basis estimation shrinks functional parameter dimension to $p=1,\dots,P$8 irrespective of $p=1,\dots,P$9 (Kowal et al., 2014).

A more general Bayesian function-space formulation places the dynamic linear model directly on separable Banach spaces: $t=1,\dots,T$0 with Gaussian-process induced covariance operators. Under consistent discretizations, the paper derives Gaussian filtering updates and employs standard MCMC under finite-grid approximations. This makes Kalman filtering and smoothing on function spaces available in hierarchical settings after basis expansion or discretization (Petris, 2013).

In the dual-factor model, temporal dependence is concentrated in $t=1,\dots,T$1, and forecasting is implemented by fitting a VAR model to the stacked factor matrices, with lag order selected by information criteria such as AIC. In the foreign-exchange multilevel FPCA framework, by contrast, the principal component scores are forecast with univariate ARIMA$t=1,\dots,T$2, chosen automatically via the Hyndman–Khandakar algorithm; exponential smoothing is also reported to yield qualitatively similar results (Tang et al., 2021, Shang et al., 2021).

The frequency-domain MDKLE framework estimates scores through a Bayesian MAP step with Whittle prior,

$t=1,\dots,T$3

solved by gradient-based optimization, with per-iteration complexity $t=1,\dots,T$4. The end-to-end procedure proceeds from mean removal and autocovariance estimation to spectral smoothing, marginalization, eigen-decomposition, phase optimization, filter construction, prior construction, score extraction, and diagnostics (Guo et al., 24 Mar 2026).

The particle-number application forecasts each score sequence independently using robust univariate methods such as exponential smoothing, while allowing VAR$t=1,\dots,T$5 as an extension. Dynamic updating with newly observed partial curves is implemented through OLS, ridge, or penalized least squares regression on FPC bases; empirically, ridge is reported as the most accurate and stable (Shang et al., 2 Oct 2025).

The additive surface model is fitted by iterative backfitting with blockwise convex subproblems after introducing auxiliary parameters $t=1,\dots,T$6. Predictor blocks are updated by weighted $t=1,\dots,T$7 proximal steps; zero groups are dropped; and a post-selection refit without sparsity penalty is used to debias the active surfaces (Wang et al., 28 Apr 2025).

5. Forecasting, uncertainty quantification, and diagnostics

Forecasting in multilevel functional time series typically reconstructs future curves from forecast scores, forecast states, or forecast multivariate score vectors.

For the MFDLM, under linear-Gaussian assumptions,

$t=1,\dots,T$8

and therefore

$t=1,\dots,T$9

Credible bands are formed from posterior samples of $j$00 and latent states, and cross-series dependencies can be summarized through common-trend parameters $j$01 or off-diagonal entries of $j$02. Standardized residuals

$j$03

are used for outlier detection and regime diagnostics. Model checking includes DIC or marginal likelihood, residual checks, and posterior predictive checks (Kowal et al., 2014).

The dual-factor model reconstructs the $j$04-th population by

$j$05

Interval forecasts are obtained by nonparametric bootstrap on in-sample forecast errors and evaluated by the interval score

$j$06

The foreign-exchange study instead reports mean absolute forecast error, mean squared forecast error, and mean mixed error, and evaluates predictive superiority by the Model Confidence Set procedure with $j$07 and $j$08 statistics (Tang et al., 2021, Shang et al., 2021).

The spectral integration framework reconstructs or forecasts functions through the truncated filter expansion

$j$09

and then models $j$10 by VAR$j$11 or VARMA. The MAP step naturally imputes scores under irregular sampling and noise, so missing functional values are reconstructed through the same filter expansion (Guo et al., 24 Mar 2026).

The particle-number study develops two distribution-free interval procedures. The first calibrates a width parameter $j$12 in pointwise intervals $j$13 by minimizing the absolute gap between empirical and nominal coverage on a validation set. The second applies split conformal prediction using absolute residual conformity scores and the empirical $j$14 quantile. The paper states that split conformal prediction provides finite-sample marginal coverage guarantees under exchangeability (Shang et al., 2 Oct 2025).

A recurrent distinction concerns coherence. The foreign-exchange multilevel FPCA framework explicitly states that there is no reconciliation; coherence across maturities is achieved by the common/residual decomposition and joint estimation, not by summing constraints (Shang et al., 2021). By contrast, the frequency-domain paper proposes that filtered score representations could be used to perform hierarchical reconciliation and notes that integrating spectra first may reduce reconciliation burden (Guo et al., 24 Mar 2026). This suggests that “coherence” in the multilevel functional literature may refer either to structural sharing inside the model or to explicit aggregation constraints imposed after modeling.

6. Applications, empirical findings, and methodological boundaries

Applications span macro-finance, neuroscience, mortality, environmental monitoring, and option-implied volatility.

In multi-economy yield curves, the MFDLM uses maturity $j$15 months as the functional domain and weekly changes in yield curves for $j$16 economies. A common basis yields interpretable components resembling level, slope, curvature; common-trend dependence identifies direct coupling to U.S. factors; and posterior means show BOE and BOC closely track Fed on level and slope, while ECB does so less strongly. Standardized residuals reveal weeks, including late 2008, when other economies diverge from Fed-driven trends. In local field potential analysis, the same framework models log-spectra and squared coherence over $j$17 Hz; feature-binding trials exhibit increased squared coherence and spectral power in Theta, Alpha, and Beta bands around event time, indicating increased synchronization and activity in PFC/PPC (Kowal et al., 2014).

In Japanese subnational mortality, the dual-factor model reports more accurate point and interval forecasts than several alternatives. The reported mean RMSFE values are: females, FDFM $j$18 versus TSHDFTS $j$19, HDFTS $j$20, HDFFM $j$21, MFM $j$22, CPD $j$23, Tucker $j$24; males, FDFM $j$25 versus TSHDFTS $j$26, HDFTS $j$27, HDFFM $j$28, MFM $j$29, CPD $j$30, Tucker $j$31. The same study translates improved mortality forecasts into a life annuity pricing scheme and reports, for a female aged $j$32 in $j$33, true $j$34, FDFM $j$35 with error $j$36 million under the example portfolio (Tang et al., 2021).

The foreign-exchange implied-volatility study finds that dynamic FPCA generally improves out-of-sample forecast accuracy, and more specifically that the dynamic univariate functional time-series method shows the greatest improvement. Using $j$37, DFTS is in the superior set in $j$38 cases, DMLFTS in $j$39, and static FTS in $j$40; using $j$41, DFTS is superior in $j$42, DMLFTS in $j$43, and static FTS in $j$44. A stylized trading strategy based on daily implied-volatility forecasts yields statistically significant positive mean daily returns for EUR-GBP $j$45, modest positive for EUR-USD $j$46, and negative for EUR-JPY $j$47 over the reported period (Shang et al., 2021).

In hourly PM$j$48 trajectories in Baoding, the spectral integration method reports, for forecasting with $j$49, Spectral MPCA NMSPE $j$50 versus DFPCA $j$51, PADA $j$52, LRFPCA $j$53, and VMFPCA $j$54. For imputation in the same case study, Spectral MPCA NMSE is $j$55 versus PADA $j$56, DFPCA $j$57, VMFPCA $j$58, and LRFPCA $j$59. The paper interprets this as superior forecasting and competitive imputation under integrated spectra across subjects (Guo et al., 24 Mar 2026).

In London intraday particle number size distributions, the MLFTS framework uses $j$60 particle size categories and $j$61 hourly time points per day, with approximately $j$62 weeks of data. Weekly curve construction with $j$63-point functions tends to minimize MAPE relative to per-day $j$64-point curves; Factor+MLFTS slightly improves over MLFTS alone; and dynamic updating with partially observed intraday data yields substantial gains, with regression-based ridge updates outperforming block-moving and OLS/PLS alternatives. At $j$65 nominal coverage, conformal and sd calibration perform similarly; at $j$66, sd calibration is slightly tighter with better calibration in several cases (Shang et al., 2 Oct 2025).

The interpretable additive model for Japanese mortality reports, at lag $j$67, FBM MAFE $j$68, MSFE $j$69, compared with UFTS $j$70, MFTS $j$71, and MLFTS $j$72. For Tokyo, global sparsity selects Aomori, Iwate, and Wakayama as influencing prefectures, while local sparsity identifies age-specific effects such as ages $j$73–$j$74 and $j$75–$j$76 in Aomori affecting Tokyo’s mortality surface (Wang et al., 28 Apr 2025).

Across the literature, limitations are stated clearly. Common examples are time-invariant basis functions or loadings, linear-Gaussian dynamics, white-noise or white-in-time idiosyncratic errors, and stationarity assumptions for long-run covariance or spectral operators [(Kowal et al., 2014); (Tang et al., 2021); (Guo et al., 24 Mar 2026)]. Proposed extensions include non-Gaussian observation errors, time-varying bases $j$77, covariate-driven dynamics in $j$78 or $j$79, nonlinear state-space or factor models, correlated $j$80, hierarchical shrinkage across levels, and integrated reconciliation procedures [(Kowal et al., 2014); (Tang et al., 2021); (Wang et al., 28 Apr 2025)].

A common misconception is that multilevel functional time series is merely multivariate FPCA applied to stacked curves. The cited formulations instead distinguish common bases from population-specific back-loadings, common trends from residual trends, shared spectral filters from subject-specific scores, and cross-level coefficient surfaces from post hoc reconciliation. Another misconception is that hierarchical coherence must always be imposed after forecasting; some frameworks achieve coherence through shared latent structure, while others explicitly leave reconciliation outside the model (Shang et al., 2021, Guo et al., 24 Mar 2026).

Taken together, the literature portrays multilevel functional time series as a family of models for dependent curves observed across time and hierarchy, with common themes of functional dimension reduction, cross-sectional borrowing of strength, and explicit dynamic modeling. This suggests that the field is best understood not as a single estimator, but as a set of closely related inferential strategies for allocating variation across shared and level-specific functional dynamics.