## Forecasting functional time series using weighted likelihood methodology

Ufuk Beyaztas Department of Statistics Bartin University

## Han Lin Shang

Research School of Finance, Actuarial Studies and Statistics Australian National University

August 2, 2019

## Abstract

Functional time series whose sample elements are recorded sequentially over time are frequently encountered with increasing technology. Recent studies have shown that analyzing and forecasting of functional time series can be performed easily using functional principal component analysis and existing univariate/multivariate time series models. However, the forecasting performance of such functional time series models may be affected by the presence of outlying observations which are very common in many scientific fields. Outliers may distort the functional time series model structure, and thus, the underlying model may produce high forecast errors. We introduce a robust forecasting technique based on weighted likelihood methodology to obtain point and interval forecasts in functional time series in the presence of outliers. The finite sample performance of the proposed method is illustrated by Monte Carlo simulations and four real-data examples. Numerical results reveal that the proposed method exhibits superior performance compared with the existing method(s).

Keywords: Bootstrap; Functional principal components; Functional time series; Weighted likelihood

## 1 Introduction

A functional time series arises when functional objects are collected sequentially over time. In other words, let { yt } ( t = 1, · · · , N ) denote a sequence of observed functional data; then, it is termed an functional time series if each functional element yt ( t ) is defined on a bounded interval with a continuous variable t ; t ∈ [ a , b ] . Denote by { yt ( xj ) } , for t = 1, · · · , N and j = 1, · · · , J , the noisy functional time series data observed at time points { x 1, · · · , xJ } . It is assumed that the functional time series is characterized by a continuous smooth function f and an error process # such that

$$y _ { t } ( x _ { j } ) = f _ { t } ( x _ { j } ) + \sigma _ { t } ( x _ { j } ) \varepsilon _ { t, j }$$

where # t j , is an independently and identically distributed (i.i.d.) Gaussian random variable with zero mean and unit variance and s t ( xj ) allowing for heteroskedasticity.

Forecasting unobservable future realizations of functional time series is of great interest. In practice, forecasts can be obtained in the form of point and/or interval forecasts. A point forecast corresponds to an estimate (conditionally on the available data) of the unknown future realization of the underlying process. However, this approach may not produce reliable inferences for future observations since it does not provide any information about the degree of uncertainty associated with the point forecasts. By contrast, interval forecasts, such as the prediction interval, provide better inferences taking into account the uncertainty associated with point forecasts; see for example Chatfield (1993), Kim (2001), and Jore et al. (2010). In the context of functional time series, Hyndman and Ullah (2007) propose a functional data approach to obtain point and interval forecasts for age-specific mortality and fertility rates observed over time. While doing so, they consider that the functional data { yt ( xj ) } consist of random functions separated by consecutive and continuous time intervals. Their approach is as follows: ( ) Approximate the i smooth functions { f t ( x ) } separately using a nonparametric smoothing technique on each function yt ( x ) ( t = 1, · · · , N ). ( ii ) Decompose the smoothed functional time series data into K orthonormal principal components f k ( k = 1, · · · , K ) and associated uncorrelated scores b k t , ( t = 1, · · · , N ) using a functional principal component model. ( iii ) Apply a univariate time series model to each principal component score to obtain their future values. ( iv ) Calculate the future realization of the functional time series by multiplying principal components with the forecasted part of the principal component scores. ( v ) Obtain the prediction intervals under the assumption of Gaussian distributed error terms. This approach (or a modified version) has received extensive attention in the literature and has successfully been used in many areas; see for example, Hyndman and Shang (2009), Shang and Hyndman (2011), Shang (2013), Kosiorowski (2014), Gao and Shang (2017), Wagner-Muns et al. (2019), Curceac et al. (2019) and references therein.

The aforementioned literature uses non-robust time series estimation methodologies when

forecasting future values of the principal component scores. However, this approach is affected by the outliers which are common in real data sets; see, for example, Shang (2019). An outlier is an observation that has been generated by a stochastic process with a distribution different from that of the vast majority of the remaining observations (Ra˜a n et al. , 2015). In the context of functional time series, three types of outliers are observed: (1) magnitude outlier , which is a point far from the bulk of the data; (2) shape outlier , which falls within the range of the data but differs in shape from the bulk of the data; and (3) the combination of both outlier types. For more information about the outliers in functional time series, see Febrero et al. (2007), Hyndman and Shang (2010), Sun and Genton (2011) and Ra˜ na et al. (2015).

In the case of outlying observation(s), non-robust techniques produce biased estimates, and high forecasting errors correspond to outlying observations. In such cases, the high forecasting errors may severely affect point forecasts as well as prediction intervals and could lead to unreliable inferences. We propose a robust functional time series forecasting method based on the minimum density power divergence estimator of Basu et al. (1998). The proposed method uses a modified version of the usual maximum likelihood score equations, called weighted score equations, to estimate the model parameters. The weighted score equations are defined as a function of the Pearson residuals for which large values are obtained when the observations diverge from the underlying model. Hence, this approach makes it possible to check whether the maximum likelihood estimators are affected by a set of observations that are inconsistent with the model, and provides robust estimates by downweighting such observations. It also provides weighted residuals, which are used to obtain point and/or interval forecasts. We use several Monte Carlo experiments and real-data examples to compare the finite sample performance of the proposed and existing methods. Our numerical records, which are discussed in Sections 3 and 4, reveal that the proposed method provides finite sample performance competitive with that of the existing methods when no outlier is present in the observed data. In addition, when outliers are present in the data, its performance is shown to be superior compared with that of available techniques.

The remaining of this paper is organized as follows. Section 2 provides an overview of the functional time series methods considered and the weighted likelihood estimation methodology. Several Monte Carlo experiments under different scenarios are conducted to evaluate the finite sample performance of the proposed method, and the results are presented in Section 3. Section 4 reports the findings obtained by applying the proposed method to some real-data examples. Section 5 concludes the paper.

## 2 Methodology

Let us consider a sequence of stationary functional time series { yt ( x ) : t ∈ Z , x ∈ J } where J is a bounded interval. It is assumed that the functions yt are elements of the metric, semi-metric, Hilbert or the Banach space, in general. We assume that the functions are elements of a squareintegrable function y ∈ L 2 [ 0, 1 ] residing in Hilbert space H satisfying ∫ y 2 ( x dx ) &lt; ¥ with an inner product 〈 y z , 〉 = ∫ y x z x dx ( ) ( ) , ∀ y z , ∈ L 2 [ 0, 1 . Denote by ] ( W S , , P ) the probability space where W , S and P represent the sample space, s -algebra on W and the probability measure on S , respectively. Then the random functional variable y is defined as y : ( W S , , P ) → H so that y is assumed to be an element of L 2 and y -1 ( B ) ∈ S where B is a Borel set of the Borel s -algebra generated by L 2 . We further assume that the random variable y ∈ L 2 ( W ) with finite second-order moment is a second-order stochastic process so that E [ | y | 2 ] = ∫ W | y | 2 dP &lt; ¥ . The mean and covariance functions of the random variable y are defined as in (2) and (3), respectively; see Ramsay and Dalzell (1991).

$$\mu ( x ) \ = \ E \left [ y ( x ) \right ] = \int _ { \cap } y ( x ) d P$$

$$C ( x, s ) \ & = \ \stackrel { \dots } { E } \left [ \left ( y ( x ) - \stackrel { \dots } { \mu ( x ) } \right ) \left ( y ( s ) - \mu ( s ) \right ) \right ] \\ & = \ \int _ { \Omega } \left [ \left ( y ( x ) - \mu ( x ) \right ) \left ( y ( s ) - \mu ( s ) \right ) \right ] d P$$

$$& \cdot ( x ) ] = \int _ { \Omega } y ( x ) d P & & ( 2 ) \\ & y ( x ) - \mu ( x ) ) \left ( y ( s ) - \mu ( s ) \right ) \right ] \\ & \cdot ( y ( x ) - \mu ( x ) ) \left ( y ( s ) - \mu ( s ) \right ) \right ] d P$$

Let { yt ( x ) : t = 1, · · · , N x , ∈ J } be an observed functional time series of size N with the same distribution as y . Then, the sample mean and sample covariance functions are given by (4) and (5), respectively.

$$\bar { \mu } ( x ) \ = \ N ^ { - 1 } \sum _ { t = 1 } ^ { N } y _ { t } ( x ) _ { \quad, }$$

$$\hat { C } ( x, s ) \ = \ ( N - 1 ) ^ { - 1 } \sum _ { t = 1 } ^ { N } \left ( y _ { t } ( x ) - \bar { y } ( x ) \right ) \left ( y _ { t } ( s ) - \bar { y } ( s ) \right )$$

Functional principal component analysis is frequently used to analyse functional time series. Briefly, it represents the data by a linear combination of orthonormal principal components f k and their associated scores b k t , ( k = 1, 2, · · · ). In doing so, it decomposes the covariance operator given in (3) into orthogonal bases of eigenfunctions. Let y k and l k , respectively, denote the k th eigenfunction and eigenvalue. Then, the covariance operator is decomposed as follows:

$$C ( x, s ) = \sum _ { k = 1 } ^ { \infty } \lambda _ { k } \psi _ { k } ( x ) \psi _ { k } ( s ).$$

The k th principal component score is then defined as b k t , = ∫ yt ( x ) y k ( x dx ) . In what follows, the

random functions are expressed using Karhunen-Lo` eve expansion as:

$$y _ { t } ( x ) = \sum _ { k = 1 } ^ { \infty } \beta _ { k, t } \psi _ { k } ( x )$$

See Ramsay and Silverman (2002), Ramsay et al. (2009) and Shang (2014) for more details about functional principal component analysis and its practical demonstration.

Let { yt ( xj ) } , for t = 1, · · · , N and j = 1, · · · , J be the observed functional time series (with noise) characterized by a continuous smooth function f and an error process # :

$$y _ { t } ( x _ { j } ) = f _ { t } ( x _ { j } ) + \sigma _ { t } ( x _ { j } ) \varepsilon _ { t, j }$$

Each function can be decomposed using a basis function expansion, as follows:

$$f _ { t } ( x ) = \hat { \mu } ( x ) + \sum _ { k = 1 } ^ { K } \hat { \phi } _ { k } ( x ) \hat { \beta } _ { k, t } + \epsilon _ { t } ( x )$$

where ˆ m ( x ) = N -1 GLYPH&lt;229&gt; N t = 1 f t ( x ) is the estimated mean of functions, ˆ f k ( x ) and ˆ b k t , denote the estimated functional principal components and their associated principal component scores, respectively, e t ( x ) is the error function with mean zero and K is the number of basis functions with K &lt; N . The choice of the basis function in (6) is arbitrary under orthogonality restriction. As Hyndman and Ullah (2007) highlight, each principal component score can be forecasted independently using univariate time series methods since they are uncorrelated to each other. They also note that there may be cross-correlations at non-zero lags but these can be negligible. However, multivariate time series methods, such as vector autoregression, can be used to take into account the cross-correlations; see Aue et al. (2015). Let h represent the forecast horizon. Denote by ˆ b k N , + h N | the h -step-ahead forecast of b k N , + h conditionally on the available data up to time N . Then, the point predictor of yN + ( h x ) is obtained by multiplying the forecasted principal component scores with the estimated functional principal components, as follows:

$$\hat { y } _ { N + h | N } = \hat { \mu } ( x ) + \sum _ { k = 1 } ^ { K } \hat { \phi } _ { k } ( x ) \hat { \beta } _ { k, N + h | N }$$

As noted in Section 1, point forecasts do not provide any information about the uncertainty of future realization of the functional time series. Conversely, the prediction interval is capable of producing valid inferences taking into account the uncertainty of each forecast. Hyndman and Ullah (2007) propose a ( 1 -a ) 100% prediction interval for yN + h under the assumption of normality. However, this approach may seriously be affected by any departure from the normality assumption, which is not known in practice. In such cases, the bootstrap method is commonly used method

to overcome this issue since it does not require full knowledge of the underlying distributional assumption. Therefore, we consider only the bootstrap prediction interval. Hyndman and Shang (2009) introduce a bootstrap approach to construct prediction intervals for the future values of the functional time series. Based on the functional time series model discussed above, we incorporate three error sources to construct the bootstrap prediction interval: ( ) smoothing error i # t j , given in (1), ( ii ) the error caused by the functional principal component decomposition e t ( x ) in (6), and ( iii ) the error occurring owing to forecasting principal component scores b k t , . Let x k h t , , denote the h -step-ahead forecast error ˆ x k h t , , = ˆ b k t , -ˆ b k t t , | -h ( t = h + 1, · · · , N ). Then, the algorithm of the bootstrap procedure proposed by Hyndman and Shang (2009) is as follows.

Step 1. Smooth the entire functional time series yt ( x ) = f t ( x ) + s t ( x ) # t for t = 1, · · · , N to obtain smooth functions f t ( x ) , smoothing error ˆ # t and variance component ˆ s t ( x ) .

Step 2. Decompose the smooth functions into k = 1, · · · , K orthonormal functional principal components and associated principal component scores and obtain the fitted functions ˆ f t ( x ) . Then, calculate the error functions ˆ e t ( x ) = f t ( x ) -ˆ f t ( x ) .

Step 3. Obtain h -step-ahead forecasts of the principal component scores b k t , , ˆ b k N , + h N | for k = 1, · · · , K and t = N + 1, · · · , N + h using a univariate time series method. In addition, obtain the h -step-ahead forecast errors ˆ x k h t , , .

## Step 4. Calculate the future bootstrap values as follows:

$$\hat { y } _ { N + h | N } ^ { * } ( x ) = \hat { \mu } ( x ) + \sum _ { k = 1 } ^ { K } \hat { \beta } _ { k, N + h | N } ^ { * } + \hat { \epsilon } _ { N + h | N } ^ { * } ( x ) + \hat { \sigma } _ { N + h } ^ { * } ( x ) \hat { \varepsilon } _ { N + h | N, j } ^ { * }$$

where ˆ b ∗ k N , + h N | = ˆ b k N , + h N | + x k h , , ∗ , x k h , , ∗ is an i.i.d. random sample from { x k h t , , } , ˆ e ∗ N + h N | ( x ) , ˆ # ∗ N + h N j | , and ˆ s ∗ N + h ( x ) are random samples with replacement from { ˆ e t ( x ) } , { ˆ # t } and { ˆ s t ( x ) } , respectively.

Step 5. Repeat Step 4. B times by drawing random samples of the error terms to obtain B sets of bootstrap replicates of ˆ y ∗ N + h N | ( x ) , { ˆ y ∗ ,1 N + h N | ( x ) , · · · , ˆ y ∗ , B N + h N | ( x ) } for each h , where B denotes the number of bootstrap simulations.

The ( 1 -a ) 100% bootstrap prediction intervals for yN + h are then obtained by the a /2 th and ( 1 -a /2 ) th quantiles of the bootstrap replicates { ˆ y ∗ ,1 N + h N | ( x ) , · · · , ˆ y ∗ , B N + h N | ( x ) } . This approach works well when the functional time series is relatively smooth, and no outlier is present in the data. It has successfully been used in a wide range of applications; see, for example, Hyndman et al. (2013), Aneiros et al. (2013), Husin et al. (2016), Aneiros et al. (2016), Canale and Vantini (2016), Shang (2017), Kearney et al. (2018) and Shang et al. (2018). By contrast, the traditional univariate

time series methods used for modelling the functional principal component scores may produce biased estimates as well as high forecasting errors, leading to unreliable results when outliers are present in the data. To overcome this issue, we propose a robust forecasting functional time series method. It is based on replacing the traditional estimators and residuals by the minimum power divergence estimator of Basu et al. (1998) and weighted residuals, respectively.

## 2.1 The weighted likelihood methodology

To employ the weighted likelihood methodology in the forecasting of functional time series, we consider the stationary autoregressive model of order p (AR( p )). Let us suppose that the k th principal component score b k t , ( t = 1, · · · , N ) is characterized by a zero mean AR( p ) process, as follows:

$$\beta _ { k, t } = \phi _ { 1 } \beta _ { t - 1, k } + \cdots, \phi _ { p } \beta _ { t - p, k } + \zeta _ { t }, \ t = 1, \cdots, N$$

where F = ( f 1, · · · , f p ) denotes the parameter vector and z t is an i.i.d. white noise sequence with mean zero and variance s z 2 . Under the assumption of normality, the probability density function for the model (7) is given by

$$P ( \zeta _ { N } | \Phi ; \sigma _ { \zeta } ^ { 2 } ) = \left ( \frac { 1 } { 2 \pi \sigma _ { \zeta } ^ { 2 } } \right ) ^ { - \frac { N - p } { 2 } } \exp \left \{ - \frac { 1 } { 2 \sigma _ { \zeta } ^ { 2 } } \sum _ { t = p + 1 } ^ { N } \zeta _ { t } ^ { 2 } ( \Phi ) \right \}$$

where z t ( F ) = b k t , -GLYPH&lt;229&gt; p i = 1 f b i k t , -i and z N = 〈 z p + 1, · · · , z N } . Let glyph[lscript] ( z t ( F ) ; s z 2 ) = log L ( z t ( F ) ; s z 2 ) denote the conditional log-likelihood function. Then, the maximum likelihood estimators of F and s z 2 conditionally on the first p observations, respectively, are obtained by the solution of the score functions u F and u sz , as follows:

$$u _ { \Phi } \left ( \zeta _ { t } ( \Phi ) ; \sigma _ { \zeta } ^ { 2 } \right ) \ & = \ \frac { \partial } { \partial \Phi } \ell \left ( \zeta _ { t } ( \Phi ) ; \sigma _ { \zeta } ^ { 2 } \right ) \\ u _ { \sigma _ { \zeta } } \left ( \zeta _ { t } ( \Phi ) ; \sigma _ { \zeta } ^ { 2 } \right ) \ & = \ \frac { \partial } { \partial \sigma _ { \zeta } } \ell \left ( \zeta _ { t } ( \Phi ) ; \sigma _ { \zeta } ^ { 2 } \right )$$

The weighted likelihood methodology is proposed by Markatou (1996) and Basu et al. (1998) to construct efficient and robust estimators by replacing the usual score functions with weighted score equations that measure the discrepancy between the estimated and hypothesized model densities. It has been extended to a wide variety of statistical inference problems, see, for example, Agostinelli and Markatou (2001), Agostinelli (2002a,b), Agostinelli (2003) and Agostinelli and Bisaglia (2010). Let us consider the probability density function given by (8) and let f ∗ ( z t ( F ) , ˆ FN ( F ) ) and m ∗ ( z t ( F ) , s z 2 ) denote a kernel density estimator with bandwidth h based on the empirical

distribution function ˆ FN ( F ) and the kernel smoothed model density, as follows:

$$f ^ { * } \left ( \zeta _ { t } ( \Phi ), \hat { F } _ { N } ( \Phi ) \right ) \ & = \ \int k \left ( \zeta _ { t } ( \Phi ) ; r, g \right ) d \hat { F } _ { N } ( r ; \Phi ) \\ m ^ { * } \left ( \zeta _ { t } ( \Phi ), \sigma _ { \zeta } ^ { 2 } \right ) \ & = \ \int k \left ( \zeta _ { t } ( \Phi ) ; r, g \right ) d M ( r ; \sigma _ { \zeta } ^ { 2 } )$$

where M ( s z 2 ) and k ( z t ( F ) ; r , g ) represent the normal distribution function with mean zero and variance s z 2 and a kernel density with bandwidth g , respectively. Denote by d t ( z t ( F ) ; M ( s z 2 ) , ˆ FN ( F ) ) = f ∗ ( z t ( F ) , ˆ FN ( F ) ) m ∗ ( z t ( F ) , s z 2 ) and w d ( t ) = min { 1, [ A ( d t )+ ] + 1 d t + 1 } the Pearson residual and the weight function, respectively, where [ · ] + denotes the positive part and A ( ) · is the residual adjustment function. Note that we consider the Hellinger residual adjustment function of Lindsay (1994), A ( d ) = 2 [( d + ) 1 0.5 -1 . Based on these definitions, the conditional weighted likelihood of the parameters ] F and sz are obtained by solving the following estimating equations:

$$& ( N - p ) ^ { - 1 } \sum _ { t = p + 1 } ^ { N } \omega ( \delta _ { t } ) u _ { \Phi } \left ( \zeta _ { t } ( \Phi ) ; \sigma _ { \zeta } ^ { 2 } \right ) \\ & ( N - p ) ^ { - 1 } \sum _ { t = p + 1 } ^ { N } \omega ( \delta _ { t } ) u _ { \sigma _ { \zeta } } \left ( \zeta _ { t } ( \Phi ) ; \sigma _ { \zeta } ^ { 2 } \right )$$

Now let ˆ F w = ( ˆ f w 1 , · · · , ˆ f w p ) denote the estimated parameter vector of AR( p ) process using the weighted likelihood methodology. Denote by ˆ b w k N , + h N | = GLYPH&lt;229&gt; p i = 1 ˆ f w i ˆ b w k N , + h - | i N the weighted likelihood version of h -step ahead forecast of b k N , + h . Then, the weighted likelihood-based h -stepahead point forecasts and bootstrap prediction intervals for y N + h N | are obtained similarly as in the algorithm given in the previous subsection. In the proposed forecasting strategy, the following holds for the final weights obtained from the full model, ˆ w t , when the model is accurately specified and no outlying observation is present in the data; sup t | ˆ w t - | -1 p → 0, see Agostinelli (1998). Following by Agostinelli (2002b), it can be shown that ∣ ∣ ∣ GLYPH&lt;229&gt; N t = 1 ˆ w z t t ( ˆ F w ) 2 -GLYPH&lt;229&gt; N t = 1 z t ( ˆ F ) 2 ∣ ∣ ∣ = op ( N ) . This result indicates that ( ) the weighted likelihood method tends to perform similar to the i maximumlikelihood method when no outlier is present in the data, and ( ii ) the weighted likelihood is expected to have better performance than the maximum likelihood method when the data are contaminated by the outlier(s) (since it downweights the 'bad' forecast errors caused by the contamination).

## 3 Numerical results

This section reports the finite sample performance of the proposed functional time series forecasting method via several Monte Carlo experiments. Throughout the experiments, two different

simulation scenarios are considered: (1) The data are generated from a relatively smooth process with no outlying observations, and (2) g = [ 1%,5%,10% of the generated data is contaminated by ] the deliberately inserted magnitude and shape outlier(s). The following process is used to generate the data; yt ( xj ) = 15 + cos ( p j /4 ) + N ( 0, 0.15 2 ) ( t = 1, · · · , N = 100 and j = 1, · · · , J = 12). The magnitude outlier(s) are generated by contaminating N g randomly selected function(s) by a random function y c t ( xj ) = | N ( 0.75, 0.15 2 ) | , that is yt ( xj ) = yt ( xj ) + y c t ( xj ) . Conversely, the shape outliers are generated from the following process: 15 + sin ( p j /4 ) + N ( 0, 0.15 2 ) . Examples of the simulated data (with magnitude and shape outliers) are presented in Figure 1.

Figure 1 Example of the simulated functional time series (grey lines) with outliers (black lines): magnitude outlier (first row) and shape outlier (second row).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00336v1-figure-80b58b87-dd1e-4001-87a4-dadc3c3d3652.jpg)

We divide the generated data into two parts so that we can construct the model using the first N -h functions to obtain point forecasts and the last h functions for the bootstrap prediction intervals, where h = [ 1, 5, 10 . ] To construct an functional time series model, first, the noisy functional time series is converted to a smooth function by the smoothing spline method. Then, the hybrid principal component model of Hyndman and Ullah (2007) is used to decompose

the smooth functions. The nominal level a is set to 0.05 to calculate 95% bootstrap prediction intervals. For each scenario, MC = 1, 000 Monte Carlo simulations with B = 999 bootstrap resamples are performed, and only the first K = 3 principal components and their scores are used to obtain forecasts. The performance of the proposed method is compared with that of the autoregressive integrated moving average (ARIMA) model as used by Hyndman and Ullah (2007) and most of the references cited. To evaluate the forecasting performance of the proposed method, we calculate several performance metrics, including the average mean squared error (AMSE), coverage probability (Cp) and the average interval scores (Score), as follows:

$$\begin{array} { l l l } \text{:average probability (Cp) and the average interval scores (Score), as follows: } \\ \\ \text{AMSE} & = & \frac { 1 } { h \times 1 2 } \sum _ { i = 1 } ^ { h } \sum _ { j = 1 } ^ { 1 2 } \left \{ y _ { N + h } ( x _ { j } ) - \hat { y } _ { N + h | N } ( x _ { j } ) \right \} ^ { 2 } \\ \\ \text{CP} & = & \frac { 1 } { h \times 1 2 } \sum _ { i = 1 } ^ { h } \sum _ { j = 1 } ^ { 1 2 } 1 \left \{ \hat { y } _ { N + h | N } ^ { *, 0. 0 2 5 } ( x _ { j } ) \leq y _ { N + h } ( x _ { j } ) \leq \hat { y } _ { N + h | N } ^ { *, 0. 9 7 5 } ( x _ { j } ) \right \} \\ \text{Score} & = & \frac { 1 } { h \times 1 2 } \sum _ { i = 1 } ^ { h } \sum _ { j = 1 } ^ { 1 2 } \left \{ \left ( \hat { y } _ { N + h | N } ^ { *, 0. 9 7 5 } ( x _ { j } ) - \hat { y } _ { N + h | N } ^ { *, 0. 0 2 5 } ( x _ { j } ) \right ) \\ \\ & + & \frac { 2 } { \alpha } \left ( \hat { y } _ { N + h | N } ^ { *, 0. 0 2 5 } ( x _ { j } ) - y _ { N + h } ( x _ { j } ) \right ) 1 \left \{ y _ { N + h } ( x _ { j } ) < \hat { y } _ { N + h | N } ^ { *, 0. 0 2 5 } ( x _ { j } ) \right \} \\ \\ & + & \frac { 2 } { \alpha } \left ( y _ { N + h } ( x _ { j } ) - \hat { y } _ { N + h | N } ^ { *, 0. 9 7 5 } ( x _ { j } ) \right ) 1 \left \{ y _ { N + h } ( x _ { j } ) > \hat { y } _ { N + h | N } ^ { *, 0. 9 7 5 } ( x _ { j } ) \right \} \right \} \\ \\ \text{where} \, 1 \left \{ \cdot \right \} \text{ denotes the indicator function and } \hat { y } _ { N ^ { \alpha } _ { \perp H | N } } ^ { *, \alpha } ( x _ { j } ) \text{ is the } \alpha ^ { \text{th} } \text{ quantile of the bootstra} \text{,} \end{array}$$

where 1 {·} denotes the indicator function and ˆ y ∗ , a N + h N | ( xj ) is the a th quantile of the bootstrap replicates.

The results obtained from the Monte Carlo experiments are reported in Table 1. Our findings show that, regardless of the forecast horizon, both weighted likelihood and ARIMA-based forecasting models tend to have similar AMSE values for the point forecasts when no outliers are present in the data. However, the proposed method produces better coverage probabilities and interval scores than those obtained from ARIMA. These results indicate that compared with ARIMA, we can obtain more reliable and more precise prediction intervals using the weighted likelihood-based bootstrap method even if the outliers do not contaminate the data. When the data has outlier(s) (magnitude and/or shape), both methods still have similar AMSE values for the point forecasts when h = 1 but the proposed method has better coverage probabilities and interval scores compared with those of ARIMA. By contrast, for long-term forecast horizons ( h = 5 and h = 10), the proposed method has smaller AMSE values than ARIMA, and the difference becomes more prominent as the forecast horizon increases. For the prediction intervals, while it seems that both methods have similar coverage performance, the weighted likelihood-based bootstrap method produces significantly narrower prediction intervals than those of ARIMA. This is because the bootstrap method based on ARIMA is considerably affected by the large forecast errors produced by the outliers. However, the proposed method downweights the effects of

outliers, and the structure of the bootstrap prediction intervals are not distorted.

Table 1 Simulation results: Estimated AMSE, Cp and Score values for the proposed (weighted likelihood estimator (WLE)) and ARIMA models when the data have magnitude (MO) and shape (SO) outliers.

| Contamination   | h      | Outlier   | Method   |   AMSE |     Cp |   Score |
|-----------------|--------|-----------|----------|--------|--------|---------|
|                 | h      |           | WLE      | 0.0228 | 0.9365 |  0.7244 |
|                 |        |           | ARIMA    | 0.0229 | 0.9119 |  0.7419 |
|                 | h = 5  |           | WLE      | 0.0228 | 0.936  |  0.7258 |
|                 | h = 5  |           | ARIMA    | 0.0228 | 0.9122 |  0.7413 |
|                 | h = 10 |           | WLE      | 0.0229 | 0.935  |  0.7312 |
| 0%              | h = 10 |           | ARIMA    | 0.0228 | 0.9106 |  0.7492 |
| 1%              | h = 1  | MO        | WLE      | 0.0229 | 0.9375 |  0.7299 |
| 1%              | h = 1  |           | ARIMA    | 0.0229 | 0.9245 |  0.7507 |
| 1%              | h = 1  | SO        | WLE      | 0.0238 | 0.9153 |  0.7525 |
| 1%              | h = 1  |           | ARIMA    | 0.0242 | 0.9086 |  0.7749 |
|                 | h = 5  | MO        | WLE      | 0.0231 | 0.9438 |  0.7303 |
|                 | h = 5  |           | ARIMA    | 0.0246 | 0.9553 |  1.1374 |
|                 | h = 5  | SO        | WLE      | 0.0225 | 0.9366 |  0.7315 |
|                 | h = 5  |           | ARIMA    | 0.0258 | 0.9413 |  1.3039 |
|                 | h = 10 | MO        | WLE      | 0.0242 | 0.9467 |  0.7404 |
|                 | h = 10 |           | ARIMA    | 0.0296 | 0.9594 |  1.2357 |
|                 |        | SO        | WLE      | 0.0228 | 0.9373 |  0.7398 |
|                 |        |           | ARIMA    | 0.0343 | 0.935  |  1.4226 |
|                 | h = 1  | MO        | WLE      | 0.0229 | 0.9378 |  0.7277 |
|                 | h = 1  |           | ARIMA    | 0.023  | 0.9246 |  0.7415 |
|                 |        | SO        | WLE      | 0.0233 | 0.9167 |  0.7537 |
|                 |        |           | ARIMA    | 0.0242 | 0.9089 |  0.7847 |
|                 | h = 5  | MO        | WLE      | 0.0233 | 0.942  |  0.7329 |
| 5%              | h = 5  |           | ARIMA    | 0.0249 | 0.9557 |  1.1344 |
|                 | h = 5  | SO        | WLE      | 0.0297 | 0.9306 |  0.7375 |
|                 | h = 5  |           | ARIMA    | 0.0263 | 0.9342 |  1.3158 |

Continued on next page

| Contamination   | h      | Outlier   | Method   |   AMSE |     Cp |   Score |
|-----------------|--------|-----------|----------|--------|--------|---------|
|                 | h = 10 | MO        | WLE      | 0.0233 | 0.9507 |  0.7298 |
|                 |        |           | ARIMA    | 0.0295 | 0.9592 |  1.239  |
|                 |        | SO        | WLE      | 0.0247 | 0.9346 |  0.7583 |
|                 |        |           | ARIMA    | 0.0361 | 0.9322 |  1.4349 |
|                 | h = 1  | MO        | WLE      | 0.0227 | 0.9377 |  0.7265 |
|                 |        |           | ARIMA    | 0.0228 | 0.9244 |  0.7476 |
|                 |        | SO        | WLE      | 0.0228 | 0.9197 |  0.7601 |
|                 |        |           | ARIMA    | 0.0238 | 0.9071 |  0.8843 |
|                 |        | MO        | WLE      | 0.0232 | 0.9437 |  0.7315 |
| 10%             | h = 5  |           | ARIMA    | 0.0249 | 0.9546 |  1.1469 |
|                 |        | SO        | WLE      | 0.0231 | 0.9225 |  0.7652 |
|                 |        |           | ARIMA    | 0.0276 | 0.9217 |  1.2838 |
|                 |        | MO        | WLE      | 0.0236 | 0.9491 |  0.7324 |
|                 | h = 10 |           | ARIMA    | 0.0303 | 0.9582 |  1.2457 |
|                 |        | SO        | WLE      | 0.0232 | 0.9286 |  0.7769 |
|                 |        |           | ARIMA    | 0.0374 | 0.9338 |  1.4522 |

The finite-sample performance of a maximum likelihood based forecasting method, such as ARIMA, depends on the magnitude of outliers. Via several Monte-Carlo experiments, we examine the effects of outlier size on forecasting accuracy of the ARIMA and weighted likelihood based functional time series methods. Our results indicate that the performance of the ARIMA becomes worse as the magnitude of the outliers increases while the proposed weighted likelihood based functional time series method produces consistent results. For instance, when the magnitude outliers are generated from y c t ( xj ) = | N ( 3.75, 0.15 2 ) | , the one-step-ahead forecasting performances of the ARIMA and proposed weighted likelihood based functional time series methods are presented in Table 2. Compared to the results reported in Table 1, Table 2 shows that when [ 1%,5%,10% of the data are contaminated with larger magnitude outliers, the ARIMA method ] produced about [ 1.40, 2.76, 7.91 ] times larger AMSE than the ones obtained when the data are contaminated with small magnitude outliers and [ 1.06, 5.64, 5.96 ] times larger score values. On the other hand, the proposed method produced almost the similar AMSE and score values in both cases. Similar results can be obtained from the corresponding author upon request for the case when the data have large shape outliers.

Table 2 Simulation results: Estimated AMSE, Cp and Score values for the proposed (weighted likelihood estimator (WLE)) and ARIMA models when the data have large magnitude outliers and forecast horizon h = 1.

| Contamination   | Method   |   AMSE |     Cp |   Score |
|-----------------|----------|--------|--------|---------|
| 1%              | WLE      | 0.03   | 0.9352 |  0.7575 |
| 1%              | ARIMA    | 0.0321 | 0.91   |  0.8003 |
| 5%              | WLE      | 0.0229 | 0.9533 |  0.7564 |
| 5%              | ARIMA    | 0.0636 | 0.971  |  4.1816 |
| 10%             | WLE      | 0.0305 | 0.9612 |  0.9618 |
| 10%             | ARIMA    | 0.1805 | 0.978  |  4.4548 |

## 4 Real-data examples

This section evaluates the finite sample performance of the proposed method using four environmental datasets-hourly bare soil temperature, wind speed, solar radiation, and wind chill from 1 May 2017 to 31 July 2017 (92 days in total)-which are collected from the Michigan weather station (the data are obtained from North Dakota Agricultural Weather Network Center https://ndawn.ndsu.nodak.edu/ ). The functional time series representation of the datasets is presented in Figure 2. It is clear from this figure that the bare soil temperature functions are reasonably smooth and have no clear outlier. Conversely, the functions of the functional time series of wind speed, solar radiation, and wind chill are noisy, and all three datasets have several outlying functions. For all four datasets, we obtain only one-step-ahead forecasts. The forecasting performance of the ARIMA and of the proposed method are compared using the rolling holdout testing samples, in line with Shang (2019) and Hyndman and Ullah (2007); thus, 80% of the datasets are used as training samples and the remaining 20% are used for validation. In the modelling step, the number of principal components is also determined based on rolling holdout validation samples. The optimal value of K is found as 4 for the bare soil temperature dataset and 6 for the other three datasets.

The obtained AMSE and Score values are presented in Figure 3. The results indicate that compared with ARIMA, the proposed method produces slightly better AMSE values for the one-step-ahead point forecast. The results also show that the weighted likelihood-based bootstrap method produces better prediction intervals than those obtained from the ARIMA-based bootstrap procedure especially for the datasets contaminated by outliers. Figure 4 presents a graphical representation of the one-step-ahead forecasts for which the functional time series models are

Figure 2 Functional time series representations of the environmental variables.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00336v1-figure-b4df4665-db0a-4d6e-8188-03a4bead14dc.jpg)

constructed based on the first 91 functions. This figure supports the results presented in Figure 3.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00336v1-figure-08f74c0b-ac85-4318-89a0-7474f72dbecf.jpg)

Solar radiation

Wind chill

Figure 3 Obtained AMSE and Score values of the ARIMA and proposed (WLE) models for the real-data examples.

## 5 Conclusion

functional time series are frequently observed in many scientific fields owing to novel data collection tools. Consequently, several functional time series models have been developed to analyse such datasets and to obtain forecasts of their unobserved realizations. Recent studies have shown that the univariate time series models together with functional principal component regression can be used to obtain valid point forecasts and bootstrap prediction intervals. However,

Figure 4 Obtained one-step-ahead point forecasts: ARIMA (blue solid lines) and weighted likelihood (red solid lines), and bootstrap 95% prediction intervals: ARIMA (blue dashed lines) and weighted likelihood (red solid lines), together with the observed function (black solid lines).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00336v1-figure-27d8815d-155b-4475-91f7-9a424c52eb50.jpg)

the traditional univariate time series models such as ARIMA, which are commonly used in the analyses, may severely be affected in the presence of outliers, leading to poor forecasting results. Therefore, we propose a robust functional time series forecasting approach based on the weighted likelihood methodology. We evaluate the finite sample performance of the proposed method using several Monte Carlo experiments and four environmental datasets. The numerical results produced showed that the proposed method is a good competitor for ARIMA and would be widely adopted because of its narrower prediction intervals.

## References

Agostinelli, C., 1998. Inferenza statistica robusta basata sulla funzione di verosimiglianza pesata: alcuni sviluppi , Phd thesis, Dipartimento di Scienze Statistiche, Universita di Padova.

Agostinelli, C., 2002a. Robust model selection in regression via weighted likelihood methodology, Statistics &amp; Probability Letters , 56 (3), 289-300.

Agostinelli, C., 2002b. Robust stepwise regression, Journal of Applied Statistics , 29 (6), 825-840.

Agostinelli, C., 2003. Robust time series estimation via weighted likelihood, in: R. Dutter, P. Filzmoser, U. Gather, and P.J. Rousseeuw, eds., Developments in Robust Statistics , Heidelberg: Springer.

Agostinelli, C. and Bisaglia, L., 2010. ARFIMA processes and outliers: A weighted likelihood approach, Journal of Applied Statistics , 37 (9), 1569-1584.

Agostinelli, C. and Markatou, M., 2001. Test of hypotheses based on the weighted likelihood methodology, Statistica Sinica , 11 (2), 499-514.

Aneiros, G., Vilar, J., and Ra˜a, P., 2016. Short-term forecast of daily curves of electricity demand n and price, International Journal of Electrical Power &amp; Energy Systems , 80, 96-108.

Aneiros, G., Vilar, J.M., Cao, R., and Roque, A.M.S., 2013. Functional prediction for the residual demand in electricity spot markets, IEEE Transactions on Power Systems , 28 (4), 4201-4208.

Aue, A., Norinho, D.D., and Hormann, S., 2015. On the prediction of stationary functional time series, Journal of the American Statistical Association , 110, 378-392.

Basu, A., Harris, I.R., Hjort, N.L., and Jones, M.C., 1998. Robust and efficient estimation by minimising a density power divergence, Biometrika , 85 (3), 549-559.

Canale, A. and Vantini, S., 2016. Constrained functional time series: Applications to the Italian gas market, International Journal of Forecasting , 32 (4), 1340-1351.

Chatfield, C., 1993. Calculating interval forecasts, Journal of Business and Economics Statistics , 11 (2), 121-135.

Curceac, S., Ternynck, C., Ouarda, T.B.M.J., Chebana, F., and Niang, S.D., 2019. Short-term air temperature forecasting using nonparametric functional data analysis and SARMA models, Environmental Modelling &amp; Software , 111, 394-408.

| Febrero, M., Galeano, P., and Gonzalez-Manteiga, W., 2007. Outlier detection in functional data by depth measures, with application to identify abnormal NOx levels, Environmetrics , 19 (4), 331-345.   |
|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Gao, Y. and Shang, H.L., 2017. Multivariate functional time series forecasting: Application to age-specific martality rates, Risks , 5 (21), 1-18.                                                       |
| Husin, W.Z.W., Ramli, N.M., Zainol, M.S., and Ghani, N.Z., 2016. Forecasting Malaysian mortality rates: Afunctional data analysis approach, Advanced Science Letters , 2 (12), 4259-4262.                |
| Hyndman, R.J., Booth, H., and Yasmeen, F., 2013. Coherent mortality forecasting: The product-ratio method with functional time series models, Demography , 50 (1), 261-283.                              |
| Hyndman, R.J. and Shang, H.L., 2009. Forecasting functional time series, Journal of the Korean Statistical Society , 38 (3), 199-211.                                                                    |
| Hyndman, R.J. and Shang, H.L., 2010. Rainbow plots, bagplots, and boxplots for functional data, Journal of Computational and Graphical Statistics , 19 (1), 29-45.                                       |
| Hyndman, R.J. and Ullah, M.S., 2007. Robust forecasting of mortality and fertility rates: Afunc- tional data approach, Computational Statistics &Data Analysis , 51 (10), 4942-4956.                     |
| Jore, A.S., Mitchell, J., and Vahey, S.P., 2010. Combining forecast densities from vars with uncertain instabilities, Journal of Applied Econometrics , 25 (4), 621-634.                                 |
| Kearney, F., Cummins, M., and Murphy, F., 2018. Forecasting implied volatility in foreign exchange markets: Afunctional time series approach, The European Journal of Finance , 24, 1-8.                 |
| Kim, J.H., 2001. Bootstrap after bootstrap prediction intervals for autoregressive models, Journal of Business and Economic Statistics , 19 (1), 117-128.                                                |
| Kosiorowski, D., 2014. Functional regression in short-term prediction of economic time series, Statistics in Transition , 15 (4), 611-626.                                                               |
| Lindsay, B.G., 1994. Efficiency versus robustness: the case for minimum hellinger distance and related methods, Annals of Statistics , 22 (2), 1018-1114.                                                |
| Markatou, M., 1996. Robust statistical inference: weighted likelihoods or usual mestimation?, Communications in Statistics - Theory and Methods , 25 (11), 2597-2613.                                    |
| Ra˜a, n P., Aneiros, G., and Vilar, J.M., 2015. Detection of outliers in functional time series, Environ- metrics , 26 (3), 178-191.                                                                     |

| Ramsay, J.O. and Dalzell, C., 1991. Some tools for functional data analysis, Journal of the Royal Statistical Society, Series B , 53 (3), 539-572.                                                                          |
|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Ramsay, J.O., Hooker, G., and Graves, S., 2009. Functional Data Analysis with R and MATLAB , New York: Springer.                                                                                                            |
| Ramsay, J.O. and Silverman, B.W., 2002. Applied Functional Data Analysis: Methods and Case Studies , New York: Springer.                                                                                                    |
| Shang, H.L., 2013. Functional time series approach for forecasting very short-term electricity demand, Journal of Applied Statistics , 40 (1), 152-168.                                                                     |
| Shang, H.L., 2014. Asurvey of functional principal component analysis, AStA Advances in Statistical Analysis , 98 (2), 121-142.                                                                                             |
| Shang, H.L., 2017. Forecasting intraday S&P500 index returns: Afunctional time series approach, Journal of forecasting , 36 (7), 741-455.                                                                                   |
| Shang, H.L., 2019. Arobust functional time series forecasting method, Journal of Statistical Computa- tion and Simulation , 89 (5), 795-814.                                                                                |
| Shang, H.L. and Hyndman, R.J., 2011. Nonparametric time series forecasting with dynamic updating, Mathematics and Computers in Simulation , 81 (7), 1310-1324.                                                              |
| Shang, H.L., Yang, Y., and Kearney, F., 2018. Intraday forecasts of a volatility index: Functional time series method with dynamic updating, Annals of Operations Research , in press, doi:10.1007/s10479- 018-3108-4.      |
| Sun, Y. and Genton, M.G., 2011. Functional boxplot, Journal of Computational and Graphical Statistics , 20 (2), 316-334.                                                                                                    |
| Wagner-Muns, I.M., Guardiola, I.G., Samaranayke, V.A., and Kayani, W.I., 2019. Afunctional data analysis approach to traffic volume forecasting, IEEE Transactions on Intelligent Transportation Systems , 19 (3), 878-888. |
Ufuk Beyaztas Han Lin Shang
Forecasting functional time series using weighted likelihood methodology
Functional time series whose sample elements are recorded sequentially over
time are frequently encountered with increasing technology. Recent studies have
shown that analyzing and forecasting of functional time series can be performed
easily using functional principal component analysis and existing
univariate/multivariate time series models. However, the forecasting
performance of such functional time series models may be affected by the
presence of outlying observations which are very common in many scientific
fields. Outliers may distort the functional time series model structure, and
thus, the underlying model may produce high forecast errors. We introduce a
robust forecasting technique based on weighted likelihood methodology to obtain
point and interval forecasts in functional time series in the presence of
outliers. The finite sample performance of the proposed method is illustrated
by Monte Carlo simulations and four real-data examples. Numerical results
reveal that the proposed method exhibits superior performance compared with the
existing method(s).