Interpretable models for forecasting high-dimensional functional time series

Interpr etable models for for ecasting high-dimensional functional time series Han Lin Shang Department of Actuarial Studies and Business Analytics Macquarie University Cristian F . Jim ´ enez-V ar ´ on * Department of Mathematics University of Y ork Abstract W e study the modeling and for ecasting of high-dimensional functional time series, which can be temporally dependent and cr oss-sectionally correlated. W e implement a functional analysis of variance (F ANOV A) to decompose high-dimensional functional time series, such as subnational age- and sex-specific mortality observed over years, into two distinct components: a deterministic mean structur e and a residual process varying over time. Unlike purely statistical dimensionality- reduction techniques, the F ANOV A decomposition provides a direct and interpr etable framework by partitioning the series into effects attributable to data-specific factors, such as r egional and sex-level variations, and a grand functional mean. From the residual process, we implement a functional factor model to capture the remaining stochastic trends. By combining the forecasts of the residual component with the estimated deterministic structur e, we obtain the forecasted curves for high-dimensional functional time series. Illustrated by the age-specific Japanese subna- tional mortality rates fr om 1975 to 2023, we evaluate and compare the accuracy of the point and interval forecasts acr oss various forecast horizons. The results demonstrate that leveraging these interpretable components not only clarifies the underlying drivers of the data but also impr oves forecast accuracy , providing mor e transpar ent insights for evidence-based policy decisions. Keywords: one-way functional analysis of variance; two-way functional analysis of variance; functional factor model; functional panel data; subnational age-specific mortality rates * Department of Mathematics, Ian W and Building, University of Y ork, Deramore Lane, Y ork YO10 5GH, United Kingdom; Email: cristian.jimenezvar on@york.ac.uk 1 1 Introduction Recent advances in computer storage and r ecording facilitate the pr esence of functional data acr oss several scientific fields, from cosmic demographics (see, e.g., Loredo et al. , 2025 ) to small-area official estimation (see, e.g., Maiti et al. , 2016 ). When a functional variable is measured over time, it gives rise to functional time series analysis (see Horv ´ ath and Kokoszka , 2012 ; Kokoszka and Reimherr , 2017 , for a compr ehensive review), wher e the crux is to captur e temporal dependence. When one observes multiple functional variables, each has its own temporal dependence; together , they also exhibit cross-sectional dependence. That gives rise to high-dimensional functional time series (HDFTS) (see, e.g., Gao et al. , 2019 ; Chang et al. , 2025 ), also known as functional panel data (see, e.g., Jim ´ enez-V ar ´ on et al. , 2025 ). Applications include temperature curves in hundreds of weather stations (see, e.g., Delaigle et al. , 2019 ); r eturn curves in finance are typically available for hundreds of stocks (see, e.g., Leng et al. , 2026 ). HDFTS can be expressed as Y t = ( Y 1 t , . . . , Y Nt ) ⊤ , i = 1, . . . , N , t = 1, . . . , T , with Y it = ( Y it ( u ) , u ∈ C i ) ∈ H i , where H i is a Hilbert space defined as a set of measurable and squar e-integrable functions on a compact set C i . A distinct feature of HDFTS is that the number of cross sections can exceed the number of curves, i.e., N > T . In the HDFTS literature, Zhou and Dette ( 2023 ) derived Gaussian and multiplier bootstrap approximations for the sums of HDFTS. Using these appr oximations, they construct joint simulta- neous confidence bands for the mean functions and develop a hypothesis test to assess whether the mean functions in the panel dimension exhibit parallel behavior . Hallin et al. ( 2023 ) investigated the r epresentation of HDFTS using a factor model, identifying conditions on the eigenvalues of the covariance operator crucial for establishing the existence and uniqueness of the factor model. T ang et al. ( 2022 ) studied clustering HDFTS, while Li et al. ( 2024 ) proposed hypothesis tests for change-point detection and estimation, and further clustering of common change points in HDFTS using an information criterion. In estimating HDFTS models, T avakoli et al. ( 2023 ) developed a functional factor model with functional factor loadings and a matrix of real-valued factors. W ith real-valued factor loadings and functional factors, Guo et al. ( 2026 ) studied a different factor r epr esentation for decomposing HDFTS. By unifying both models, Leng et al. ( 2026 ) considered functional factor loadings and functional factors. W ith a two-stage repr esentation, Happ and Greven ( 2018 ) applied a univari- 2 ate functional principal component analysis to extract a set of principal components and their associated scores; by stacking the scores, one can further reduce dimensionality through prin- cipal component analysis. In the same vein, Gao et al. ( 2019 ) stacked the principal component scores acr oss all populations by their eigenvalue orders; using a factor model, the dimensional- ity can be further reduced. While these factor-based approaches are computationally efficient for dimensionality reduction, a significant challenge in their widespread application is that the resulting factors and loading matrices are often difficult to interpr et ( Liu and W ang , 2026 ). In many high-dimensional settings, the latent factors ar e purely mathematical constructs that do not directly corr espond to observable categorical drivers, such as geographic r egions or demographic cohorts. This ”black-box” nature can limit the utility of the model for practitioners who require an understanding of the specific sour ces of variation, such as sex-specific tr ends or regional disparities, that drive the overall series. T o address this limitation, we consider a structural appr oach that prioritizes interpretability without sacrificing the flexibility of factor modeling. Following earlier work by Jim ´ enez-V ar ´ on et al. ( 2024 ), we consider a two-way functional analysis of variance (ANOV A) to decompose an HDFTS. Functional ANOV A models evaluate the functional effects of categorical variables (known as factors) by determining how functions differ across their levels. Differing from Jim ´ enez-V ar ´ on et al. ( 2024 ), our two-way functional ANOV A incorporates an interaction term, modeled via a one-way functional ANOV A, ther eby providing a more flexible repr esentation of cross-group heterogeneity than existing appr oaches. Functional ANOV A models are not new , since they are useful for analyzing data across a range of applications, such as human tactile perception ( Spitzner et al. , 2003 ), menstrual cycle data ( Brumback and Rice , 1998 ), and circadian rhythms with random ef fects and smoothing-spline ANOV A decomposition ( W ang et al. , 2003 ). Kaufman and Sain ( 2010 ) established a Bayesian framework for functional ANOV A modeling to estimate the effect of geographical regions on Canadian temperatur es. Sun and Genton ( 2012 ) pr oposed a functional median polish modeling as an extension of the univariate median polish of T ukey ( 1977 ), and a functional rank test was developed to determine the significance of the ef fects of the functional main factor . While Shang ( 2025 ) consider ed a one-way functional ANOV A for modeling and forecasting a time series of Lorenz curves, Jim ´ enez-V ar ´ on et al. ( 2024 ) proposed a two-way functional ANOV A to model and 3 forecast subnational age-specific mortality observed over time. Regarding our age- and sex-specific subnational mortality rates, our modeling framework consists of three steps: (I) we begin with the two-way functional ANOV A to model mean patterns associated with region and sex; (II) the interaction term in the two-way functional ANOV A can be captured by one-way functional ANOV A; (III) the sex-specific residuals acr oss regions over time from the one-way functional ANOV A are then modeled thr ough a functional factor model. As a by-product, we want to examine whether or not the interaction term in step II can enhance point and interval forecast accuracy . This paper contributes to the literature on high-dimensional functional time series (HDFTS) by developing an interpretable and structur ed modeling framework that separates deterministic and stochastic variation. W e begin by decomposing HDFTS into data-driven functional components using a two-way functional analysis of variance (F ANOV A), which attributes systematic variation directly to observable factors, such as region and sex. W e then extend this decomposition by incorporating an interaction structur e, modeled via a one-way F ANOV A, to capture heteroge- neous effects across gr oups. The remaining variation is modeled using a functional factor model, allowing us to account for time-varying stochastic dynamics in the residual component. This unified framework combines interpr etability with flexibility . Using Japanese age and sex-specific subnational mortality data, we demonstrate that the proposed approach clarifies the sour ces of cross-r egional heter ogeneity and improves point and interval forecast accuracy r elative to existing methods. The remainder of the paper is organized as follows. Section 2 introduces the data and provides a preliminary analysis. Section 3 presents the pr oposed decomposition framework, including the two-way and one-way F ANOV A and a functional factor model. Section 4 outlines the construction of prediction intervals. Section 5 discusses model fitting, while Section 6 evaluates forecasting performance. Section 7 provides concluding r emarks and discusses future r esearch dir ections. 2 Japanese subnational age-specific mortality rates In many developed countries, such as Japan, the incr ease in longevity and the aging population have raised concerns about the sustainability of pension, health and aged car e systems (see, e.g., 4 Coulmas , 2007 ). These concerns have led to increased attention among government policymakers to accurately model and for ecast age-specific mortality . Although much ef fort has been devoted to national mortality forecasting, subnational age-specific mortality forecasts are particularly helpful for informing policy within local r egions and for allocating curr ent and futur e resour ces at the national and subnational levels. Despite their usefulness for understanding cross-r egion heterogeneity , subnational mortality data are generally of poor data quality . 0 20 40 60 80 −8 −6 −4 −2 Okinawa female data Age Log death rate (a) Raw data 0 20 40 60 80 −8 −6 −4 −2 Okinawa male data Age Log death rate (b) Raw data 0 20 40 60 80 −8 −6 −4 −2 Age Log death rate (c) Smoothed data 0 20 40 60 80 −8 −6 −4 −2 Age Log death rate (d) Smoothed data Figure 1: Age-specific raw and smoothed log mortality rates from ages 0 to 100+ between 1975 and 2023 in Okinawa. Curves are or dered chr onologically accor ding to the colors of the rainbow . The left vertical axis measures log mortality rates. Mortality rates dip in early childhood, climb in the teen years, stabilize in the early 20s, and then steadily increase with age. Obtained from Japanese Mortality Database ( 2026 ), we investigate the subnational age- and 5 sex-specific mortality rates in Japan fr om 1975 to 2023. Mortality rates are the ratios of register ed death counts to population exposur e in the r elevant year for the given age (based on a one-year age group). W e examine age groups ranging from 0 to 99 in a single year of age, with the last age group including all ages at and above 100. In Figures 1a and 1b , we present a rainbow plot of age- and sex-specific mortality in Okinawa, the most southern pr efecture in Japan. Although a decreasing tr end can be observed, it is evident that the noise level of such a series is quite high. These multiple functions ar e observed with err or at times t = 1, . . . , T and we aim to for ecast the functions for times t = T + 1, . . . , T + h , wher e h denotes a for ecast horizon. Let {Z it ( u j ) } denote the log of mortality rates, wher e ages ( u 1 , . . . , u p ) are observed discretely . W e assume that ther e are underlying L 2 continuous and smooth functions {Y it ( u j ) } such that ln Z it ( u j ) = ln Y it ( u j ) + σ t ( u j ) ε it , j , j = 1, . . . , p , where { ε it , j } are independent and identically distributed variables with zer o mean and unit vari- ance, and σ t ( u j ) allows heter oskedasticity . Because of monotonic increasing beyond age 65, we con- sider a P -spline with monotonic constraint ( W ood , 1994 ). Computationally , the smooth.demogdata function in the demography package ( Hyndman , 2026 ) is used. In Figures 1c and 1d , we display smoothed log mortality rates for Okinawa. Not only do smoothing techniques help addr ess missing values and mortality rates that exceed one, but they also enhance forecast accuracy . 3 Decomposition of high-dimensional functional time series 3.1 T wo-way functional analysis of variance Through a two-way functional ANOV A, smoothed mortality rates Y g it ( u ) can be decomposed as Y g it ( u ) = µ ( u ) + α i ( u ) + δ g ( u ) + γ ig ( u ) + X g it ( u ) , (1) where Y g it ( u ) denotes the smoothed mortality rate at age u for gender g in prefectur e i = 1, . . . , N and year t = 1, . . . , T . In this model, µ ( u ) repr esents the functional grand effect, α i ( u ) denotes 6 the functional row (pr efecture) ef fect, δ g ( u ) signifies the functional column (gender) effect, γ ig ( u ) denotes the interaction term, and X g it ( u ) repr esents the r esidual component. By grouping the interaction term and the residual component, we denote X g , † it ( u ) = γ ig ( u ) + X g it ( u ) . The functional grand, row , and column effects can be estimated by their sample means (see also Chapter 13 Ramsay and Silverman , 2006 ). They are given by b µ ( u ) = 1 N × 2 × T N X i = 1 2 X g = 1 T X t = 1 Y g it ( u ) b α i ( u ) = 1 2 × T 2 X g = 1 T X t = 1 Y g it ( u ) − b µ ( u ) b δ g ( u ) = 1 N × T N X i = 1 T X t = 1 Y g it ( u ) − b µ ( u ) . T o ensure identifiability , several constraints are in place: N X i = 1 α i ( u ) = 2 X g = 1 δ g ( u ) = 0, N X i = 1 X g , † it ( u ) = 2 X g = 1 X g , † it ( u ) = 0, ∀ t . The estimates b µ ( u ) , b α i ( u ) , and b δ g ( u ) capture the overall mean, the regional mean, and the gender mean, X g , † it ( u ) is time-varying and itself an HDFTS. Since females and males have differ ent mortality profiles, we adopt a one-way F ANOV A in Section 3.2 or a functional factor model of Leng et al. ( 2026 ) in Section 3.3 to capture the r egional pattern. 3.2 One-way functional analysis of variance The interaction term fr om the two-way F ANOV A decomposition states that the prefectur e ef fect influences each gender dif fer ently . T o estimate the interaction term, we employ a one-way func- tional ANOV A decomposition (see, e.g., Shang , 2025 ). In this framework, the observations are decomposed as follows: X g , † it ( u ) = θ g ( u ) + η g i ( u ) + R g it ( u ) for g ∈ { F, M } , (2) 7 where θ g ( u ) repr esents the functional grand effect, η g i ( u ) repr esents the row ef fect corr esponding to the prefectur e mean, and R g it ( u ) denotes the remaining r esidual component. Empirically , these first two terms in ( 2 ) can be approximated by b θ g ( u ) = 1 N × T N X i = 1 T X t = 1 X g , † it ( u ) b η g i ( u ) = 1 T T X t = 1 X g , † it ( u ) − b θ g ( u ) . 3.3 A functional factor model It is r ecognized that HDFTS ar e influenced by common functions over the temporal dimension, leading to strong cr oss-sectional dependence. For example, age- and sex-specific mortality curves collected across differ ent r egions may be influenced by common patterns in each r egion. Hence, the functional factor model plays an important r ole in modeling HDFTS. In the statistical literature, Hallin et al. ( 2023 ) and T avakoli et al. ( 2023 ) present a factor model with functional factor loadings and scalar-valued factors. Guo et al. ( 2026 ) propose a factor model with scalar-valued factor loadings and functional factors. By unifying both factor models, Leng et al. ( 2026 ) introduce a factor model with functional factor loadings and functional factors, which we adopt. For each gender g , we consider the following functional factor model: R it ( u ) = k ∗ X k = 1 Z B g ik ( u , v ) F g tk ( v ) dv + ε g it ( u ) . (3) where F g tk ( v ) repr esents time-varying factors, B g ik ( u , v ) is a continuous two-dimensional linear operator , k ∗ is the number of factors, and ε g it ( u ) is the error term of the model. By imposing a low-dimensional functional factor condition on the latent factor F g tk ( v ) , we obtain the series approximation F g tk ( v ) = Φ g k ( v ) ⊤ G g t + η g tk ( v ) , (4) where Φ g k ( v ) is a q -dimensional vector of basis functions, G g t is a q -dimensional vector of stationary random variables, and η g tk ( v ) is the approximation err or . 8 By plugging ( 4 ) into ( 3 ), we obtain R g it ( u ) = k ∗ X k = 1 Z C ∗ k B g ik ( u , v )  Φ g k ( v ) ⊤ G g t + η g tk ( v )  dv + ε g it ( u ) = k ∗ X k = 1 Z B g ik ( u , v ) Φ g k ( v ) ⊤ dv G g t + k ∗ X k = 1 Z B g ik ( u , v ) η g tk ( v ) dv + ε g it ( u ) = Λ g i ( u ) ⊤ G g t + ε ∗ g it ( u ) , where Λ g i ( u ) = P k ∗ k = 1 R B g ik ( u , v ) Φ g k ( v ) ⊤ dv and ε ∗ g it ( u ) = P k ∗ k = 1 R B g ik ( u , v ) η g tk ( v ) dv + ε g it ( u ) . T o estimate b G g t , we first estimate the covariance of R g it ( u ) by ∆ g = ( ∆ g tt ′ ) T × T with ∆ g tt ′ = 1 N N X i = 1 Z u ∈ C s R g it ( u ) R g it ′ ( u ) du . By eigenanalysis of the T × T matrix ∆ g , we obtain e G g = ( e G g 1 , . . . , e G g T ) ⊤ as a T × q matrix with columns being the eigenvectors (multiplied by √ T ) corresponding to the q largest eigenvalues of ∆ g . The factor loading functions are estimated as e Λ g i ( u ) = 1 T T X t = 1 R it ( u ) e G g t , i = 1, . . . , N , via least squares, using the normalization r estriction 1 T P T t = 1 e G g t e G g ⊤ t = I g q . T o estimate the number of factors q , we r esort to an information criterion introduced in Leng et al. ( 2026 ), which can slowly diverge to infinity . Let ν ℓ ( ∆ g /T ) be the ℓ th largest eigenvalue of ∆ g /T and define b q = arg min 1 ⩽ ℓ ⩽ q max [ ν ℓ ( ∆ g /T ) + ℓϕ NT ] − 1, where ϕ NT is the penalty term and q max is a user-specified positive integer . In practice, we set ϕ NT = [ max ( T , N )] − 1 2 and q max = T . 4 Construction of prediction intervals In statistics, for ecast uncertainty is often quantified using a statistical model, although such approaches may be vulnerable to model misspecification, selection bias, and limited finite-sample 9 validity . While bootstrapping can potentially mitigate some of these concerns (see, e.g., Papar oditis and Shang , 2023 ), it is often computationally demanding. Her e, we take a model agnostic and distribution-free approach, namely conformal prediction, to construct prediction intervals in HDFTS. Among a rich family of conformal prediction methods, we consider split and sequential conformal predictions. Specifically , for the TW A + OW A + FFM framework in ( 5 ) , we first jointly model all prefectur es and genders using two-way functional ANOV A to extract the deterministic grand and main effects. W e then addr ess the interaction term by applying a gender -specific one-way functional ANOV A, which decomposes the remaining variation into deterministic prefecture-specific means and a stochastic residual pr ocess R g it ( u ) . The integrated point forecast on the original scale is thus: b Y g i , T + h ( u ) = b µ ( u ) + b α i ( u ) + b δ g ( u ) | {z } Deterministic (TW A) + b θ g ( u ) + b η g i ( u ) | {z } Deterministic (OW A) + b R g i , T + h ( u ) | {z } Stochastic (FFM) . (5) For the TW A+FFM model, we assume that there is no deterministic interaction component ( θ g ( u ) = η g i ( u ) = 0). Consequently , the construction of the pr ediction interval utilizes the r esiduals X g , † it ( u ) obtained directly after the two-way functional ANOV A. In both modeling frameworks, we implement split and sequential conformal pr ediction meth- ods. In the split approach, we use a validation set to calibrate a tuning parameter ξ g α , i to scale pointwise summary measures (such as the standar d deviation or absolute quantiles) of the residu- als as described below . In the sequential approach, we avoid a fixed validation set by using an autoregr essive model on the absolute r esiduals to update predictive quantiles as new data arrive. Finally , these calibrated stochastic components are added back to the deterministic structure to produce the complete interval for ecasts. 4.1 Split conformal prediction The split conformal prediction method uses a validation set to calibrate the empirical coverage probability to match the nominal coverage probability closely . W e divide the 49-year sample (1975- 2023) into training, validation, and test sets with pr oportions of 60%, 20%, and 20%, respectively . Using the initial training data from 1975 to 2002, we implement an expanding-window for ecasting 10 scheme to generate h -step-ahead forecasts for the validation period fr om 2003 to 2013, for h = 1, 2, . . . , 10. Under this forecasting scheme, the training sample is pr ogressively enlar ged at each iteration. The number of curves in the validation set or test set varies with the forecast horizon h . For instance, when h = 1, there are 11 years to compute residual functions, corr esponding to the dif ferences between the observed curves in the validation set and their for ecasts, whereas when h = 10, there are only two years. From these residual functions, we compute pointwise summary measur es γ g i ( u ) , such as the pointwise standar d deviation, for each r egion i and gender g . Alternatively , for quantile-based intervals, we take the absolute residuals and calculate the 100 ( 1 − α ) % empirical quantiles, where α is the significance level, typically α = 0.05. For a given forecast horizon h , let us denote the residual functions on the original scale b ϵ g im ( u ) = Z g im ( u ) − b Z g im ( u ) for m = 1, . . . , M , when M denotes the number of years in the validation set. W ith the objective of minimizing the absolute distance between the empirical and nominal coverage pr obabilities, we seek to determine a tuning parameter ξ g α , i such that 100 ( 1 − α ) % of the residual functions satisfy − ξ g α , i γ g i ( u ) ⩽ b ϵ g i ( u ) ⩽ ξ g α , i γ g i ( u ) , where b ϵ g i ( u ) = [ b ϵ g i ,1 ( u ) , . . . , b ϵ g i , M ( u )] . By the law of large numbers, when M is r easonably large, one could achieve P h − ξ g α , i γ g i ( u ) ⩽ Z g iT + h ( u ) − b Z g iT + h ( u ) ⩽ ξ g α , i γ g i ( u ) i ≈ 1 M M X m = 1 1  − ξ g α , i γ g i ( u ) ⩽ b ϵ g im ( u ) ⩽ ξ g α , i γ g i ( u )  , where Z g iT + h ( u ) denotes the holdout age-specific mortality rates in prefectur e i and gender g in year T + h , b Z g iT + h ( u ) denotes the corresponding forecasts in the original scale (after taking exponential back-transformation), and 1 { · } repr esents the binary indicator function. 4.2 Sequential conformal prediction The split conformal pr ediction method r equires a validation data set to calibrate a tuning parameter . Not only can it be time-consuming, but it also reduces the samples used to predict the test set. W ithout the need for a validation set, this sequential conformal prediction can automatically 11 tune the pr edictive quantiles of the absolute residual functions as new data arrive. For instance, with the last 10 years as the test set, we use the other years to compute the absolute r esiduals {| e g i 3 ( u j ) | , . . . , | e g iι ( u j ) |} for a given age u j . W e requir e at least the first two curves to pr oduce a forecast, so the r esiduals begin with the 3 rd curve. At the quantile of 1 − α , we fit a quantile regr ession on lagged residuals, where the order of autoregression, denoted by p, is determined by an information criterion, such as the Akaike information criterion ( Akaike , 1969 ). Conditional on the most recent p number of absolute residuals as input, we predict a one-step-ahead quantile, denoted b q g iι + 1, α ( u j ) , where ι repr esents the data in the end of the training period. The prediction intervals ar e then given by b Z g i ( ι + 1 ) ( u j ) ± b q g i ( ι + 1 ) , α ( u j ) . Once the actual curve Z g i ( ι + 1 ) ( u ) arrives, we can update the absolute residual | e g i ( ι + 1 ) ( u j ) | and refit. 5 Model fitting 5.1 Model fitting via two-way functional analysis of variance Using subnational Japanese age- and sex-specific log mortality rates, we apply a two-way func- tional ANOV A to decompose HDFTS into differ ent mean effects. In Figure 2 , we display the grand mean effect, the r ow effect (prefectur e), the column effect (sex) and the time-varying residuals that can contain the interaction term between prefectur e and sex. 0 20 40 60 80 −8 −6 −4 −2 Overall mean Age Ground effect (a) Functional grand effect −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0 20 40 60 80 10 20 30 40 Mean by prefecture Age Prefecture (b) Functional row ef fect 0 20 40 60 80 −1.0 −0.5 0.0 0.5 1.0 Mean by sex Age Column effect Male Female (c) Functional column effect Figure 2: T wo-way functional analysis of variance decomposition: Functional grand effect, capturing the overall age profile; functional r ow effect, capturing the means across pr efectures; and functional column effect, capturing the means acr oss sexes. The functional grand effect r eveals an overall mean trend. The functional row effect r eveals 12 cross-r egional heter ogeneity , particularly for younger ages. The functional column effect shows a contrast between the male and female data. 5.2 Model fitting via one-way functional analysis of variance Like many developed nations, females generally have lower mortality rates than males at all ages in Japan ( Ikeda et al. , 2011 ). Since females and males are biologically differ ent ( Rogers et al. , 2010 ), we treat each of them separately . Thus, we apply one-way functional ANOV A to decompose the residual functions into a grand ef fect, a r ow effect, and the r emainder term. 0 20 40 60 80 −2e−14 0e+00 2e−14 4e−14 Age Grand eff ect (a) Functional grand effect of the female data −0.10 −0.05 0.00 0.05 0.10 0.15 0 20 40 60 80 10 20 30 40 Age Prefecture (b) Functional row ef fect of the female data 0 20 40 60 80 −4e−14 −2e−14 0e+00 2e−14 4e−14 Age Grand eff ect (c) Functional grand effect of the male data −0.15 −0.10 −0.05 0.00 0.05 0.10 0 20 40 60 80 10 20 30 40 Age Prefecture (d) Functional row ef fect of the male data Figure 3: From the interaction term and time-varying r esiduals, we implement the one-way func- tional ANOV A for each gender to capture the functional grand ef fect and row ef fect, respectively . 13 In Figur e 3 , we pr esent the decomposition terms obtained fr om the one-way functional ANOV A. For either females or males, the functional grand effects are very close to zero. In contrast, the functional row ef fects r eveal cross-r egion heterogeneity , more so for the male data. 5.3 Model fitting via functional factor model From the r eminder term of the one-way functional ANOV A, we apply the functional factor model to extract the first set of factor loadings and its associated factors for the female and male data in Figure 4 . While the two-way and one-way functional ANOV A models extract deterministic mean effects, the functional factor model aims to model the information associated with variance. −0.6 −0.4 −0.2 0.0 0 20 40 60 80 10 20 30 40 Factor loading (F) Age Prefecture (a) First functional factor loading of the female data −0.8 −0.6 −0.4 −0.2 0.0 0 20 40 60 80 10 20 30 40 Factor loading (M) Age Prefecture (b) First set of the factors of the male data 1980 1990 2000 2010 2020 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 Y ear F actor (c) First functional factor loading of the female data 1980 1990 2000 2010 2020 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 Y ear F actor (d) First set of the factors of the male data Figure 4: V ia the functional factor model, we display the first set of factor loadings and factors, based on the estimated variance of the HDFTS residuals. 14 6 Results 6.1 Expanding-window scheme An expanding window analysis is a robust method for evaluating the stability of model parameters and the accuracy of temporal predictions. This approach assesses structural consistency by iteratively computing parameter estimates and corresponding for ecasts as the sample size grows (see Zivot and W ang , 2006 , pp. 313–314). In this study , we utilize a training set consisting of the first 39 years (1975–2013) of Japanese subnational age- and sex-specific mortality data to generate one- to 10-step-ahead for ecasts. W e then r e-estimate the time-series model parameters by incr ementing the sample size by one year at each iteration until the full data period ending in 2023 is exhausted. This recursive pr ocedure yields a total of 10 one-step-ahead forecasts, 9 two-step-ahead for ecasts, and continues down to a single 10-step-ahead for ecast. These values ar e compared with holdout observations to quantify out-of-sample accuracy . Figure 5 illustrates this expanding window scheme for the horizon h = 1, though horizons up to h = 10 are evaluated. T ime T rain F T rain F T rain F T rain F 1975:2013 1975:2014 . . . 1975:2022 T raining W indow Forecast (F) Figure 5: A diagram of the expanding-window forecast scheme. The data begin in 1975 and end in 2023. 6.2 Point forecast error metrics T o measure the point for ecast errors, we use the r oot mean square for ecast error (RMSFE) and the mean absolute forecast err or (MAFE). For each r egion i and gender g , the RMSFE and MAFE are 15 RMSFE g i ( h ) = v u u t 1 ( 11 − h ) × p 11 − h X ξ = h p X j = 1 h Z g η + ξ , i ( u j ) − b Z g η + ξ , i ( u j ) i 2 , MAFE g i ( h ) = 1 ( 11 − h ) × p 11 − h X ξ = h p X j = 1    Z g η + ξ , i ( u j ) − b Z g η + ξ , i ( u j )    , where Z g η + ξ , i ( u j ) repr esents the holdout sample for age u j and gender g in region i , and b Z g η + ξ , i ( u j ) repr esents the corr esponding point forecasts. 6.3 Comparison of point forecast accuracy In T able 1 , we report the one-step-ahead to ten-step-ahead point forecast errors for the holdout samples of the female and male data. Compared to the work of Jim ´ enez-V ar ´ on et al. ( 2024 ), our proposed methods based on functional ANOV A and the functional factor model impr ove the point forecast accuracy . Between the two frameworks, there is a slight prefer ence for including the one-way functional ANOV A to capture the interaction term between prefectur e and sex. T able 1: A veraged across 47 pr efectures, we evaluate and compare the point forecast accuracy measur ed by RMSE and MAE. Forecasting method is exponential smoothing. The method with the smallest overall error is bolded. TW A + OW A + FFM represents our method combining two-way functional ANOV A, one-way functional ANOV A and a functional factor model. TW A + FFM represents our method combining two-way functional ANOV A and a functional factor model. TW A + MFTS is a benchmark method combining two-way functional ANOV A and multivariate functional time series method in Jim ´ enez-V ar ´ on et al. ( 2024 ) Female Male Metric h TW A+ OW A+FFM TW A+ FFM TW A+ MFTS TW A+ OW A+FFM TW A+ FFM TW A+ MFTS RMSFE 1 0.0252 0.0284 0.0309 0.0448 0.0493 0.0540 2 0.0282 0.0299 0.0323 0.0440 0.0499 0.0556 3 0.0306 0.0305 0.0336 0.0423 0.0499 0.0563 4 0.0296 0.0316 0.0353 0.0452 0.0522 0.0580 5 0.0263 0.0330 0.0368 0.0528 0.0571 0.0599 6 0.0267 0.0350 0.0388 0.0612 0.0633 0.0626 7 0.0325 0.0375 0.0417 0.0694 0.0697 0.0656 Continued on next page 16 Female Male Metric h TW A+ OW A+FFM TW A+ FFM TW A+ MFTS TW A+ OW A+FFM TW A+ FFM TW A+ MFTS 8 0.0395 0.0389 0.0477 0.0626 0.0653 0.0659 9 0.0501 0.0444 0.0571 0.0427 0.0545 0.0682 10 0.0504 0.0430 0.0575 0.0427 0.0482 0.0630 Mean 0.0339 0.0352 0.0412 0.0508 0.0559 0.0609 Median 0.0301 0.0340 0.0378 0.0450 0.0533 0.0613 MAFE 1 0.0103 0.0109 0.0117 0.0172 0.0188 0.0212 2 0.0113 0.0115 0.0121 0.0172 0.0193 0.0222 3 0.0120 0.0118 0.0126 0.0170 0.0195 0.0228 4 0.0117 0.0121 0.0131 0.0185 0.0208 0.0237 5 0.0107 0.0127 0.0137 0.0218 0.0234 0.0248 6 0.0106 0.0134 0.0144 0.0251 0.0262 0.0260 7 0.0124 0.0141 0.0152 0.0282 0.0287 0.0271 8 0.0149 0.0140 0.0167 0.0251 0.0264 0.0265 9 0.0189 0.0150 0.0191 0.0173 0.0206 0.0264 10 0.0195 0.0147 0.0195 0.0163 0.0185 0.0246 Mean 0.0132 0.0130 0.0148 0.0204 0.0222 0.0245 Median 0.0118 0.0131 0.0140 0.0179 0.0207 0.0247 In Appendix A , we present an additional comparison of point for ecast errors based on the ARIMA forecasting method. 6.4 Interval forecast error metrics T o evaluate and compare interval forecast accuracy , we compute the empirical coverage pr obability (ECP), the coverage pr obability differ ence (CPD), and the interval scor e of Gneiting and Raftery ( 2007 ). For each year in the test set, the h -step-ahead pr ediction intervals are calculated at 1 − α nominal coverage probability , the lower and upper bounds, denoted by [ b Z lb, g η + ξ , i ( u j ) , b Z ub, g η + ξ , i ( u j )] , are not r equir ed to be centered ar ound the point forecasts ( Jim ´ enez-V ar ´ on et al. , 2024 ; Shang and Haberman , 2026 ). For region i and gender g , the ECP and CPD are defined as 17 ECP g α , h , i = 1 96 × ( 11 − h ) 10 X ξ = h 96 X j = 1 1  b Z lb, g η + ξ , i ( u j ) ⩽ Z g η + ξ , i ( u j ) ⩽ b Z ub, g η + ξ , i ( u j )  CPD g α , h , i =      1 96 × ( 11 − h ) 10 X ξ = h 96 X j = 1 h 1  Z g η + ξ , i ( u j ) > b Z ub, g η + ξ , i ( u j )  + 1  Z g η + ξ , i ( u j ) < b Z ub, g η + ξ , i ( u j )  i − α      , where η repr esents the year before the test data set. The ECP assesses coverage without evaluating the sharpness of the prediction interval. By combining coverage and sharpness, we use a scoring rule for the prediction interval at age u j , denoted as S α , ξ h b Z lb, g η + ξ , i ( u j ) , b Z ub, g η + ξ , i ( u j ) , Z g η + ξ , i ( u j ) i = h b Z ub, g η + ξ , i ( u j ) − b Z lb, g η + ξ , i ( u j ) i + 2 α h b Z lb, g η + ξ , i ( u j ) − Z g η + ξ , i ( u j ) i 1  Z g η + ξ , i ( u j ) < b Z lb, g η + ξ , i ( u j )  + 2 α h Z g η + ξ , i ( u j ) − b Z ub, g η + ξ , i ( u j ) i 1  Z g η + ξ , i ( u j ) > b Z ub, g η + ξ , i ( u j )  , where α is a level of significance. A veraging the number of ages and the number of years in the testing set, the mean interval score (IS) is given as S g α , h , i = 1 96 × ( 11 − h ) 10 X ξ = h 96 X j = 1 S α , ξ h b Z lb, g η + ξ , i ( u j ) , b Z ub, g η + ξ , i ( u j ) , Z g η + ξ , i ( u j ) i . Given the same ECP , the mean interval score r ewar ds narrower pr ediction intervals. 6.5 Comparison of interval forecast accuracy In T able 2 , we pr esent the one-step-ahead to ten-step-ahead interval for ecast errors based on the split and sequential conformal pr ediction methods. Between the split and sequential conformal prediction methods, the former one leads to under -estimation, in which the empirical coverage probability is smaller than the nominal one. In contrast, the sequential conformal prediction method r esults in over-estimation, in which the empirical coverage pr obability is larger than the nominal one. In the split conformal prediction method, the sd is a much better summary statistic than the quantile. 18 T able 2: Using the combination of two-way functional ANOV A, one-way functional ANOV A and functional factor model, we present the one-step-ahead to ten-step-ahead interval forecast accuracy for the female and male data at the 95% nominal coverage probability . Female Male Split (sd) Split (quantile) Sequential Split (sd) Split (quantile) Sequential h ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ARIMA 1 0.942 0.017 0.008 0.908 0.046 0.007 0.971 0.031 0.013 0.942 0.016 0.013 0.912 0.042 0.012 0.980 0.033 0.021 2 0.930 0.027 0.009 0.891 0.063 0.008 0.972 0.032 0.014 0.938 0.019 0.014 0.902 0.052 0.013 0.981 0.034 0.021 3 0.932 0.023 0.010 0.881 0.071 0.009 0.978 0.032 0.014 0.937 0.019 0.014 0.888 0.064 0.014 0.983 0.035 0.022 4 0.930 0.025 0.011 0.871 0.081 0.009 0.981 0.035 0.016 0.933 0.023 0.015 0.878 0.075 0.015 0.987 0.037 0.023 5 0.922 0.032 0.013 0.833 0.119 0.011 0.982 0.034 0.018 0.930 0.026 0.018 0.857 0.095 0.017 0.988 0.039 0.025 6 0.921 0.034 0.015 0.807 0.145 0.013 0.980 0.032 0.019 0.926 0.029 0.018 0.838 0.114 0.019 0.988 0.039 0.026 7 0.922 0.032 0.018 0.779 0.173 0.014 0.975 0.029 0.020 0.922 0.034 0.021 0.815 0.137 0.022 0.990 0.040 0.027 8 0.939 0.018 0.019 0.765 0.187 0.014 0.975 0.029 0.022 0.924 0.032 0.025 0.781 0.171 0.025 0.991 0.041 0.028 9 0.952 0.017 0.026 0.741 0.211 0.012 0.971 0.029 0.024 0.929 0.027 0.033 0.718 0.234 0.031 0.991 0.041 0.029 10 0.940 0.024 0.085 0.673 0.277 0.012 0.969 0.028 0.025 0.935 0.023 0.094 0.618 0.332 0.043 0.993 0.043 0.030 Mean 0.933 0.025 0.021 0.815 0.137 0.011 0.975 0.031 0.019 0.932 0.025 0.027 0.821 0.132 0.021 0.987 0.038 0.025 Median 0.931 0.024 0.014 0.820 0.132 0.011 0.975 0.032 0.018 0.931 0.025 0.018 0.848 0.104 0.018 0.988 0.039 0.025 ETS 1 0.927 0.028 0.008 0.882 0.072 0.007 0.963 0.030 0.013 0.940 0.017 0.013 0.910 0.044 0.012 0.979 0.033 0.020 2 0.926 0.029 0.008 0.873 0.080 0.008 0.966 0.030 0.013 0.939 0.018 0.014 0.902 0.051 0.012 0.981 0.035 0.021 3 0.928 0.027 0.009 0.862 0.090 0.008 0.968 0.030 0.015 0.940 0.018 0.014 0.889 0.063 0.013 0.983 0.035 0.022 4 0.934 0.024 0.009 0.855 0.097 0.008 0.971 0.033 0.016 0.940 0.018 0.015 0.882 0.070 0.013 0.986 0.037 0.023 5 0.932 0.026 0.009 0.838 0.114 0.008 0.971 0.033 0.018 0.938 0.021 0.016 0.863 0.089 0.015 0.986 0.038 0.025 6 0.935 0.026 0.009 0.827 0.126 0.008 0.969 0.033 0.020 0.942 0.019 0.017 0.846 0.105 0.017 0.987 0.039 0.027 7 0.941 0.022 0.012 0.815 0.137 0.009 0.968 0.032 0.020 0.943 0.018 0.021 0.831 0.121 0.019 0.985 0.037 0.026 8 0.939 0.018 0.014 0.768 0.184 0.011 0.965 0.035 0.022 0.957 0.017 0.022 0.836 0.116 0.015 0.985 0.039 0.028 9 0.931 0.025 0.022 0.676 0.276 0.017 0.959 0.036 0.023 0.962 0.019 0.032 0.814 0.138 0.017 0.987 0.039 0.029 10 0.918 0.040 0.056 0.555 0.395 0.026 0.955 0.038 0.024 0.956 0.018 0.124 0.742 0.208 0.022 0.988 0.040 0.030 Mean 0.931 0.027 0.016 0.795 0.157 0.011 0.966 0.033 0.018 0.946 0.018 0.029 0.851 0.101 0.015 0.985 0.037 0.025 Median 0.931 0.026 0.009 0.832 0.120 0.008 0.967 0.033 0.019 0.941 0.018 0.017 0.854 0.097 0.015 0.986 0.038 0.025 Using the combination of two-way functional ANOV A and functional factor model, we present 19 the one-step-ahead to ten-step-ahead interval forecast accuracy for the female and male data at the 95% nominal coverage pr obability in Appendix B . Extending the analysis of Jim ´ enez-V ar ´ on et al. ( 2024 ), we consider forecast horizons ranging from h = 1 to h = 10. The r esults of the corresponding pr ediction interval ar e reported in Appendix C . Overall, the split conformal appr oach calibrated via the standard deviation delivers the most competitive performance among the benchmark methods. Nevertheless, the proposed framework, combining a functional ANOV A decomposition with a functional factor model to capture and forecast time-varying r esidual dynamics, achieves a clear improvement in the accuracy of the interval forecast. In particular , this approach consistently yields smaller interval forecast errors compared to the method pr oposed in Jim ´ enez-V ar ´ on et al. ( 2024 ). 7 Conclusion W e present a novel statistical method for extracting patterns in high-dimensional functional time series and demonstrate its usage using Japanese subnational age- and sex-specific log mortality rates fr om 1975 to 2023. Using a two-way functional analysis of variance, we decompose the HDFTS into a functional grand effect, a row effect, a column effect, an interaction term, and time-varying residuals. The interaction term and time-varying residuals ar e further modeled via a one-way functional ANOV A; and the remainder term is then modeled by a functional factor model. T o facilitate r eproducibility , the code for computing the point and interval forecast errors is available at the following Github Repository . As a means of quantifying forecast uncertainty , we consider conformal prediction methods. Split conformal prediction requires sample splitting, which can lead to inferior interval for ecast accuracy at longer for ecast horizons. As a model-agnostic and tuning-parameter-fr ee approach, sequential conformal prediction gradually updates the predictive quantiles when new data arrive. Because it does not requir e calibration using a validation set, this conservative approach is recommended for quantifying finite-sample prediction uncertainty . There are several ways in which the current paper may be extended, and we briefly outline three: 1) In the sequential conformal pr ediction, we model the temporal dependence of the absolute residuals via an autoregr essive process in a quantile regr ession. Other time-series models can 20 also be applied. 2) While we consider modeling subnational age- and sex-specific mortality rates, one could also model life-table death counts, which themselves resemble probability density functions. 3) W e utilize a hierarchical decomposition to incorporate interaction effects, which enhances empirical forecast accuracy . While this paper focuses on the predictive benefits of such terms, developing formal hypothesis tests for interaction significance in an HDFTS framework, accounting for complex spatio-temporal dependencies, presents a valuable direction for future resear ch. Acknowledgment This resear ch is financially supported by the Australian Research Council Discovery Project DP230102250 and the Australian Resear ch Council Future Fellowship FT240100338. Jim ´ enez-V ar ´ on gratefully acknowledges support from the EPSRC NeST Pr ogramme Grant EP/X002195/1. References Akaike, H. (1969), ‘Power spectrum estimation through autor egressive model fitting’, Annals of the Institute of Statistical Mathematics 21 (1), 407–419. Brumback, B. A. and Rice, J. A. (1998), ‘Smoothing spline models for the analysis of nested and crossed samples of curves’, Journal of the American Statistical Association: Theory and Methods 93 , 961–976. Chang, J., Fang, Q., Qiao, X. and Y ao, Q. (2025), ‘On the modeling and prediction of high- dimensional functional time series’, Journal of the American Statistical Association: Theory and Methods 120 (552), 2181–2195. Coulmas, F . (2007), Population Decline and Ageing in Japan – the Social Consequences , Routledge, New Y ork. Delaigle, A., Hall, P . and Pham, T . (2019), ‘Clustering functional data into groups by using projec- tions’, Journal of the Royal Statistical Society Series B: Statistical Methodology 81 (2), 271–304. 21 Gao, Y ., Shang, H. L. and Y ang, Y . (2019), ‘High-dimensional functional time series forecasting: An application to age-specific mortality rates’, Journal of Multivariate Analysis 170 , 232–243. Gneiting, T . and Raftery , A. E. (2007), ‘Strictly proper scoring rules, prediction, and estimation’, Journal of the American Statistical Association: Review Article 102 (477), 359–378. Guo, S., Qiao, X., W ang, Q. and W ang, Z. (2026), ‘Factor modeling for high-dimensional functional time series’, Journal of Business and Economic Statistics 44 (1), 106–119. Hallin, M., Nisol, G. and T avakoli, S. (2023), ‘Factor models for high-dimensional functional time series I: Representation r esults’, Journal of T ime Series Analysis 44 , 578–600. Happ, C. and Greven, S. (2018), ‘Multivariate functional principal component analysis for data observed on different (dimensional) domains’, Journal of the American Statistical Association: Theory and Methods 113 (522), 649–659. Horv ´ ath, L. and Kokoszka, P . (2012), Inference for Functional Data with Applications , Springer , New Y ork. Hyndman, R. (2026), demography: Forecasting Mortality , Fertility , Migration and Population Data . R package version 2.0.1. URL: https://CRAN.R- project.org/package=demography Ikeda, N., Saito, E., Kondo, N., Inoue, M., Ikeda, S., Satoh, T ., W ada, K., Stickley , A., Katanoda, K., Mizoue, T ., Noda, M., Iso, H., Fujino, Y ., Sobue, T ., T sugane, S., Naghavi, M., Ezzati, M. and Shibuya, K. (2011), ‘What has made the population of japan healthy?’, The Lancet 378 (9796), 1094– 1105. Japanese Mortality Database (2026), National Institute of Population and Social Security Research . A vail- able at https://www.ipss.go.jp/p- toukei/JMD/index- en.asp (data downloaded on January 24, 2026). Jim ´ enez-V ar ´ on, C. F ., Sun, Y . and Shang, H. L. (2024), ‘Forecasting high-dimensional functional time series: Application to sub-national age-specific mortality’, Journal of Computational and Graphical Statistics 33 (4), 1160–1174. 22 Jim ´ enez-V ar ´ on, C. F ., Sun, Y . and Shang, H. L. (2025), ‘Forecasting density-valued functional panel data’, Australian & New Zealand Journal of Statistics 67 (3), 401–415. Kaufman, C. G. and Sain, S. R. (2010), ‘Bayesian functional ANOV A modeling using Gaussian process prior distributions’, Bayesian Analysis 5 , 123–149. Kokoszka, P . and Reimherr , M. (2017), Introduction to Functional Data Analysis , Chapman and Hall/CRC, New Y ork. Leng, C., Li, D., Shang, H. L. and Xia, Y . (2026), ‘Covariance function estimation for high- dimensional functional time series with dual factor structur es’, Journal of Business and Economic Statistics in press . Li, D., Li, R. and Shang, H. L. (2024), ‘Detection and estimation of structural breaks in high- dimensional functional time series’, The Annals of Statistics 52 (4), 1716–1740. Liu, X. and W ang, X. (2026), ‘Regularized estimation of the loading matrix in factor models for high-dimensional time series’. URL: https://arxiv .org/abs/2506.11232 Loredo, T ., Budav ´ ari, T ., Kent, D. and Ruppert, D. (2025), ‘Bayesian functional data analysis in astronomy’, Physical Sciences Forum 12 . Maiti, T ., Sinha, S. and Zhong, P .-S. (2016), ‘Functional mixed ef fects model for small area estima- tion’, Scandinavian Journal of Statistics 43 (3), 886–903. Paparoditis, E. and Shang, H. L. (2023), ‘Bootstrap prediction bands for functional time series’, Journal of the American Statistical Association: Theory and Methods 118 (542), 972–986. Ramsay , J. and Silverman, B. (2006), Functional Data Analysis , Springer Series in Statistics, Springer , New Y ork. Rogers, R. G., Everett, B. G., Onge, J. M. S. and Krueger , P . M. (2010), ‘Social, behavioral, and biological factors, and sex differ ences in mortality’, Demography 47 (3), 555–578. Shang, H. L. (2025), ‘For ecasting a time series of Lor enz curves: One-way functional analysis of variance’, Journal of Applied Statistics 52 (15), 2924–2940. 23 Shang, H. L. and Haberman, S. (2026), ‘Constructing pr ediction intervals for the age distribution of deaths’, Scandinavian Actuarial Journal in press . Spitzner , D. J., Marr on, J. S. and Essick, G. K. (2003), ‘Mixed-model functional ANOV A for studying human tactile perception’, Journal of the American Statistical Association: Applications & Case Studies 98 , 263–272. Sun, Y . and Genton, M. G. (2012), ‘Functional median polish’, Journal of Agricultural, Biological, and Environmental Statistics 17 , 354–376. T ang, C., Shang, H. L. and Y ang, Y . (2022), ‘Clustering and forecasting multiple functional time series’, The Annals of Applied Statistics 16 , 2523–2553. T avakoli, S., Nisol, G. and Hallin, M. (2023), ‘Factor models for high-dimensional functional time series II: Estimation and forecasting’, Journal of T ime Series Analysis 44 , 601–621. T ukey , J. (1977), Exploratory Data Analysis , Addison-W esley , Reading. W ang, Y ., Ke, C. and Brown, M. B. (2003), ‘Shape-invariant modeling of cir cadian rhythms with random effects and smoothing spline ANOV A decomposition’, Biometrics 59 , 804–812. W ood, S. N. (1994), ‘Monotonic smoothing splines fitted by cross validation’, SIAM Journal on Scientific Computing 15 (5), 1126–1133. Zhou, Z. and Dette, H. (2023), ‘Statistical inference for high-dimensional panel functional time series’, Journal of the Royal Statistical Society: Series B 85 (2), 523–549. Zivot, E. and W ang, J. (2006), Modeling Financial T ime Series with S-PLUS , Springer , New Y ork. 24 A Point forecast accuracy based on ARIMA In T able 3 , we present an additional comparison of point forecast errors based on the ARIMA forecasting method. Between the ETS and ARIMA forecasting methods, there is an advantage of using the ETS method, as it consistently yields lower mean and median RMSFE and MAFE values across the majority of for ecast horizons for both female and male datasets. T able 3: A veraged across 47 pr efectures, we evaluate and compare the point forecast accuracy measur ed by RMSE and MAE. The number of components is determined by the EVR criterion. Forecasting method is ARIMA. The method with the smallest overall error is bolded. Female Male Metric h TW A+ OW A+FFM TW A+ FFM TW A+ MFTS TW A+ OW A+FFM TW A+ FFM TW A+ MFTS RMSFE 1 0.0025 0.0028 0.0034 0.0045 0.0049 0.0055 2 0.0028 0.0030 0.0039 0.0044 0.0050 0.0057 3 0.0031 0.0031 0.0043 0.0042 0.0050 0.0059 4 0.0030 0.0032 0.0048 0.0045 0.0052 0.0061 5 0.0026 0.0033 0.0052 0.0053 0.0057 0.0063 6 0.0027 0.0035 0.0057 0.0061 0.0063 0.0066 7 0.0032 0.0037 0.0063 0.0069 0.0070 0.0071 8 0.0040 0.0039 0.0076 0.0063 0.0065 0.0080 9 0.0050 0.0044 0.0091 0.0043 0.0054 0.0094 10 0.0050 0.0043 0.0095 0.0043 0.0048 0.0094 Mean 0.0034 0.0035 0.0060 0.0051 0.0056 0.0070 Median 0.0030 0.0034 0.0055 0.0045 0.0053 0.0064 MAFE 1 0.0010 0.0011 0.0013 0.0017 0.0019 0.0021 2 0.0011 0.0011 0.0014 0.0017 0.0019 0.0022 3 0.0012 0.0012 0.0015 0.0017 0.0020 0.0023 4 0.0012 0.0012 0.0017 0.0019 0.0021 0.0024 5 0.0011 0.0013 0.0019 0.0022 0.0023 0.0025 6 0.0011 0.0013 0.0021 0.0025 0.0026 0.0026 7 0.0012 0.0014 0.0023 0.0028 0.0029 0.0028 8 0.0015 0.0014 0.0028 0.0025 0.0026 0.0031 9 0.0019 0.0015 0.0033 0.0017 0.0021 0.0036 10 0.0019 0.0015 0.0035 0.0016 0.0019 0.0037 Mean 0.0013 0.0013 0.0022 0.0020 0.0022 0.0027 Median 0.0012 0.0013 0.0020 0.0018 0.0021 0.0025 25 B Interval forecast accuracy from two-way functional ANOV A and functional factor model Using the combination of two-way functional ANOV A and functional factor model, in T able 4 , we present the one-step-ahead to ten-step-ahead interval for ecast accuracy for the female and male data at the 95% nominal coverage probability . T able 4: Interval forecast Accuracy for the female and male data: one-step-ahead to ten-step- ahead. TW A+FFM represents two-way functional ANOV A and functional factor model for the 95% nominal coverage probability . Female Male Split(sd) Split(q) Seq Split(sd) Split(q) Seq h ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ARIMA 1 0.939 0.018 0.008 0.899 0.054 0.008 0.975 0.031 0.014 0.946 0.012 0.014 0.918 0.036 0.013 0.983 0.034 0.021 2 0.932 0.023 0.009 0.888 0.066 0.009 0.976 0.032 0.014 0.938 0.019 0.014 0.905 0.049 0.014 0.984 0.034 0.021 3 0.925 0.030 0.010 0.870 0.083 0.011 0.976 0.032 0.015 0.931 0.025 0.015 0.888 0.066 0.015 0.985 0.036 0.022 4 0.921 0.033 0.011 0.855 0.097 0.013 0.979 0.035 0.016 0.923 0.033 0.016 0.871 0.083 0.016 0.989 0.039 0.023 5 0.912 0.041 0.013 0.833 0.119 0.016 0.979 0.034 0.018 0.912 0.044 0.018 0.849 0.105 0.019 0.990 0.040 0.024 6 0.909 0.044 0.017 0.815 0.137 0.018 0.977 0.033 0.019 0.905 0.050 0.019 0.823 0.130 0.022 0.990 0.040 0.025 7 0.911 0.042 0.018 0.793 0.159 0.019 0.977 0.031 0.020 0.900 0.055 0.022 0.796 0.158 0.027 0.990 0.040 0.026 8 0.933 0.022 0.018 0.787 0.165 0.016 0.975 0.030 0.021 0.898 0.057 0.025 0.764 0.189 0.030 0.991 0.041 0.027 9 0.950 0.024 0.028 0.765 0.187 0.013 0.975 0.030 0.023 0.911 0.045 0.034 0.701 0.252 0.039 0.990 0.040 0.028 10 0.908 0.048 0.088 0.695 0.255 0.012 0.975 0.029 0.024 0.896 0.057 0.092 0.610 0.340 0.048 0.992 0.042 0.029 Mean 0.924 0.033 0.022 0.820 0.132 0.013 0.976 0.032 0.018 0.916 0.040 0.027 0.812 0.141 0.024 0.988 0.039 0.024 Median 0.923 0.032 0.015 0.824 0.128 0.013 0.976 0.031 0.018 0.911 0.044 0.018 0.836 0.118 0.020 0.990 0.040 0.025 ETS 1 0.938 0.018 0.008 0.898 0.055 0.007 0.973 0.030 0.013 0.947 0.012 0.014 0.918 0.036 0.013 0.981 0.032 0.020 2 0.932 0.023 0.009 0.886 0.067 0.009 0.974 0.031 0.013 0.944 0.014 0.015 0.908 0.046 0.013 0.981 0.032 0.020 3 0.926 0.028 0.009 0.869 0.084 0.010 0.975 0.032 0.015 0.939 0.018 0.015 0.893 0.061 0.014 0.983 0.033 0.021 4 0.923 0.032 0.010 0.857 0.096 0.010 0.977 0.033 0.015 0.932 0.023 0.016 0.873 0.081 0.016 0.984 0.035 0.021 5 0.914 0.040 0.011 0.831 0.121 0.012 0.979 0.034 0.016 0.927 0.030 0.018 0.844 0.109 0.020 0.985 0.035 0.022 6 0.907 0.046 0.012 0.810 0.142 0.014 0.978 0.033 0.017 0.929 0.029 0.018 0.830 0.122 0.021 0.985 0.036 0.022 7 0.906 0.048 0.013 0.785 0.167 0.015 0.977 0.032 0.017 0.940 0.022 0.020 0.832 0.121 0.020 0.984 0.035 0.022 8 0.912 0.042 0.015 0.763 0.189 0.016 0.975 0.030 0.018 0.942 0.023 0.020 0.810 0.142 0.022 0.986 0.036 0.022 9 0.920 0.034 0.019 0.725 0.227 0.017 0.974 0.029 0.018 0.945 0.024 0.031 0.784 0.169 0.024 0.986 0.036 0.022 10 0.888 0.065 0.047 0.631 0.319 0.020 0.971 0.029 0.019 0.921 0.038 0.117 0.713 0.237 0.026 0.988 0.038 0.023 Mean 0.917 0.038 0.015 0.806 0.147 0.013 0.975 0.031 0.016 0.937 0.023 0.028 0.841 0.112 0.019 0.984 0.035 0.021 Median 0.917 0.037 0.011 0.820 0.132 0.013 0.975 0.031 0.016 0.940 0.023 0.018 0.838 0.115 0.020 0.985 0.035 0.022 26 C Interval forecast accuracy from two-way functional ANOV A of Jim ´ enez-V ar ´ on et al. ( 2024 ) and conformal prediction Extending the analysis of Jim ´ enez-V ar ´ on et al. ( 2024 ), we consider forecast horizons ranging from h = 1 to h = 10. The corr esponding prediction interval r esults are r eported in T able 5 . T able 5: Interval forecast Accuracy for the female and male data: one-step-ahead to ten-step- ahead. TW A+MFTS r epresents two-way functional ANOV A and multivariate functional time series method for the 95% nominal coverage probability . Female Male Split(sd) Split(q) Seq Split(sd) Split(q) Seq h ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ECP CPD IS ARIMA 1 0.930 0.028 0.010 0.880 0.075 0.010 0.982 0.036 0.024 0.938 0.019 0.015 0.912 0.043 0.014 0.990 0.040 0.028 2 0.919 0.037 0.026 0.860 0.095 0.013 0.981 0.035 0.023 0.932 0.023 0.016 0.901 0.053 0.015 0.990 0.040 0.027 3 0.916 0.039 0.011 0.842 0.110 0.012 0.981 0.035 0.023 0.927 0.025 0.016 0.884 0.066 0.016 0.989 0.041 0.027 4 0.913 0.041 0.013 0.824 0.128 0.013 0.982 0.035 0.024 0.920 0.033 0.016 0.870 0.081 0.017 0.990 0.041 0.028 5 0.903 0.050 0.016 0.800 0.152 0.015 0.983 0.036 0.022 0.912 0.040 0.017 0.850 0.101 0.018 0.991 0.042 0.029 6 0.904 0.051 0.016 0.780 0.173 0.017 0.982 0.035 0.021 0.910 0.043 0.019 0.824 0.127 0.020 0.989 0.041 0.027 7 0.904 0.051 0.019 0.750 0.203 0.020 0.981 0.032 0.021 0.909 0.042 0.021 0.800 0.150 0.023 0.988 0.042 0.028 8 0.915 0.041 0.024 0.717 0.236 0.023 0.979 0.031 0.021 0.915 0.038 0.024 0.768 0.182 0.026 0.988 0.042 0.028 9 0.931 0.029 0.033 0.658 0.294 0.026 0.977 0.030 0.020 0.917 0.039 0.035 0.710 0.241 0.033 0.989 0.043 0.027 10 0.913 0.041 0.075 0.564 0.386 0.032 0.976 0.028 0.020 0.916 0.038 0.089 0.611 0.339 0.045 0.987 0.043 0.027 Mean 0.915 0.041 0.024 0.767 0.185 0.018 0.981 0.033 0.022 0.920 0.034 0.027 0.813 0.138 0.023 0.989 0.041 0.028 Median 0.914 0.041 0.018 0.790 0.162 0.016 0.981 0.035 0.022 0.917 0.038 0.018 0.837 0.114 0.019 0.989 0.042 0.028 ETS 1 0.935 0.022 0.008 0.889 0.065 0.008 0.980 0.034 0.023 0.941 0.016 0.015 0.911 0.044 0.014 0.989 0.040 0.028 2 0.929 0.027 0.009 0.877 0.077 0.008 0.981 0.033 0.023 0.938 0.018 0.016 0.899 0.055 0.015 0.989 0.039 0.027 3 0.923 0.032 0.010 0.860 0.092 0.009 0.981 0.035 0.023 0.931 0.024 0.016 0.883 0.067 0.016 0.989 0.039 0.027 4 0.918 0.036 0.010 0.848 0.104 0.010 0.981 0.035 0.024 0.928 0.026 0.017 0.871 0.079 0.017 0.989 0.039 0.028 5 0.915 0.040 0.011 0.828 0.124 0.011 0.983 0.036 0.023 0.928 0.027 0.018 0.855 0.095 0.019 0.989 0.039 0.028 6 0.912 0.043 0.012 0.812 0.140 0.012 0.981 0.035 0.022 0.925 0.031 0.019 0.836 0.114 0.021 0.988 0.038 0.027 7 0.909 0.046 0.013 0.788 0.164 0.013 0.978 0.033 0.022 0.928 0.029 0.021 0.823 0.127 0.022 0.987 0.038 0.027 8 0.911 0.042 0.016 0.757 0.195 0.016 0.978 0.032 0.021 0.935 0.028 0.022 0.805 0.145 0.023 0.986 0.037 0.027 9 0.918 0.037 0.022 0.711 0.240 0.020 0.976 0.031 0.021 0.940 0.025 0.034 0.776 0.174 0.026 0.988 0.038 0.027 10 0.898 0.054 0.049 0.607 0.343 0.029 0.973 0.032 0.020 0.937 0.030 0.101 0.704 0.246 0.029 0.989 0.038 0.027 Mean 0.917 0.038 0.016 0.798 0.154 0.014 0.979 0.034 0.022 0.933 0.025 0.028 0.836 0.115 0.020 0.988 0.038 0.027 Median 0.916 0.039 0.011 0.820 0.132 0.011 0.981 0.034 0.022 0.933 0.026 0.019 0.846 0.105 0.020 0.989 0.039 0.027 27