Forecasting functional time series using weighted likelihood methodology

TL;DR

This paper introduces a robust weighted likelihood approach for forecasting functional time series, effectively handling outliers to improve prediction accuracy, demonstrated through simulations and real data.

Contribution

The paper proposes a novel weighted likelihood methodology for robust forecasting of functional time series, addressing outlier sensitivity in existing models.

Findings

Proposed method outperforms existing forecasting techniques in simulations.

Numerical results show improved accuracy in real-data applications.

Method provides reliable point and interval forecasts despite outliers.

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).

Tables

Table 1. 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
0%	$h=1$		WLE	0.0228	0.9365	0.7244
			ARIMA	0.0229	0.9119	0.7419
	$h=5$		WLE	0.0228	0.9360	0.7258
			ARIMA	0.0228	0.9122	0.7413
	$h=10$		WLE	0.0229	0.9350	0.7312
			ARIMA	0.0228	0.9106	0.7492
1%	$h=1$	MO	WLE	0.0229	0.9375	0.7299
			ARIMA	0.0229	0.9245	0.7507
		SO	WLE	0.0238	0.9153	0.7525
			ARIMA	0.0242	0.9086	0.7749
	$h=5$	MO	WLE	0.0231	0.9438	0.7303
			ARIMA	0.0246	0.9553	1.1374
		SO	WLE	0.0225	0.9366	0.7315
			ARIMA	0.0258	0.9413	1.3039
	$h=10$	MO	WLE	0.0242	0.9467	0.7404
			ARIMA	0.0296	0.9594	1.2357
		SO	WLE	0.0228	0.9373	0.7398
			ARIMA	0.0343	0.935	1.4226
5%	$h=1$	MO	WLE	0.0229	0.9378	0.7277
			ARIMA	0.0230	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.9420	0.7329
			ARIMA	0.0249	0.9557	1.1344
		SO	WLE	0.0297	0.9306	0.7375
			ARIMA	0.0263	0.9342	1.3158
	$h=10$	MO	WLE	0.0233	0.9507	0.7298
			ARIMA	0.0295	0.9592	1.2390
		SO	WLE	0.0247	0.9346	0.7583
			ARIMA	0.0361	0.9322	1.4349
10%	$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
	$h=5$	MO	WLE	0.0232	0.9437	0.7315
			ARIMA	0.0249	0.9546	1.1469
		SO	WLE	0.0231	0.9225	0.7652
			ARIMA	0.0276	0.9217	1.2838
	$h=10$	MO	WLE	0.0236	0.9491	0.7324
			ARIMA	0.0303	0.9582	1.2457
		SO	WLE	0.0232	0.9286	0.7769
			ARIMA	0.0374	0.9338	1.4522

Table 2. 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 ℎ 1 h=1 .

Contamination	Method	AMSE	Cp	Score
1%	WLE	0.0300	0.9352	0.7575
	ARIMA	0.0321	0.9100	0.8003
5%	WLE	0.0229	0.9533	0.7564
	ARIMA	0.0636	0.9710	4.1816
10%	WLE	0.0305	0.9612	0.9618
	ARIMA	0.1805	0.9780	4.4548

Equations

$$y_t(x_j)=f_t(x_j)+{\sigma }_t(x_j){\epsilon }_{t,j}$$

$$\mu (x)$$

$$C(x,s)$$

$$\bar{\mu }(x)$$

$$\hat{C}(x,s)$$

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

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

$$y_t(x_j)=f_t(x_j)+{\sigma }_t(x_j){\epsilon }_{t,j}$$

$$f_t(x)=\hat{\mu }(x)+\sum _{k=1}^K{\hat{\varphi }}_k(x){\hat{\beta }}_{k,t}+{ϵ}_t(x)$$

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

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

$${\beta }_{k,t}={\varphi }_1{\beta }_{t-1,k}+\cdots ,{\varphi }_p{\beta }_{t-p,k}+{\zeta }_t,t=1,\cdots ,N$$

$$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\}$$

$$u_{\Phi }\left({\zeta }_t(\Phi );{\sigma }_{\zeta }^2\right)$$

$$u_{{\sigma }_{\zeta }}\left({\zeta }_t(\Phi );{\sigma }_{\zeta }^2\right)$$

$$f^{\ast }\left({\zeta }_t(\Phi ),{\hat{F}}_N(\Phi )\right)$$

$$m^{\ast }\left({\zeta }_t(\Phi ),{\sigma }_{\zeta }^2\right)$$

$$(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)$$

$$=$$

Full text

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 $\left\{y_{t}\right\}$ ($t=1,\cdots,N$) denote a sequence of observed functional data; then, it is termed an functional time series if each functional element $y_{t}(\tau)$ is defined on a bounded interval with a continuous variable $\tau$; $\tau\in[a,b]$. Denote by $\left\{y_{t}(x_{j})\right\}$, for $t=1,\cdots,N$ and $j=1,\cdots,J$, the noisy functional time series data observed at time points $\left\{x_{1},\cdots,x_{J}\right\}$. It is assumed that the functional time series is characterized by a continuous smooth function $f$ and an error process $\varepsilon$ such that

[TABLE]

where $\varepsilon_{t,j}$ is an independently and identically distributed (i.i.d.) Gaussian random variable with zero mean and unit variance and $\sigma_{t}(x_{j})$ 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 $\left\{y_{t}(x_{j})\right\}$ consist of random functions separated by consecutive and continuous time intervals. Their approach is as follows:

($i$) Approximate the smooth functions $\left\{f_{t}(x)\right\}$ separately using a nonparametric smoothing technique on each function $y_{t}(x)$ ($t=1,\cdots,N$).

($ii$) Decompose the smoothed functional time series data into $K$ orthonormal principal components $\phi_{k}$ ($k=1,\cdots,K$) and associated uncorrelated scores $\beta_{k,t}$ ($t=1,\cdots,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 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ña 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 $\left\{y_{t}(x):t\in\mathcal{Z},x\in J\right\}$ where $J$ is a bounded interval. It is assumed that the functions $y_{t}$ are elements of the metric, semi-metric, Hilbert or the Banach space, in general. We assume that the functions are elements of a square-integrable function $y\in L^{2}[0,1]$ residing in Hilbert space $\mathcal{H}$ satisfying $\int y^{2}(x)dx<\infty$ with an inner product $\langle y,z\rangle=\int y(x)z(x)dx$, $\forall y,z\in L^{2}[0,1]$. Denote by $\left(\Omega,\Sigma,P\right)$ the probability space where $\Omega$, $\Sigma$ and $P$ represent the sample space, $\sigma$-algebra on $\Omega$ and the probability measure on $\Sigma$, respectively. Then the random functional variable $y$ is defined as $y:\left(\Omega,\Sigma,P\right)\rightarrow\mathcal{H}$ so that $y$ is assumed to be an element of $L^{2}$ and $y^{-1}\left(\mathcal{B}\right)\in\Sigma$ where $\mathcal{B}$ is a Borel set of the Borel $\sigma$-algebra generated by $L^{2}$. We further assume that the random variable $y\in L^{2}\left(\Omega\right)$ with finite second-order moment is a second-order stochastic process so that $E\left[|y|^{2}\right]=\int_{\Omega}|y|^{2}dP<\infty$. The mean and covariance functions of the random variable $y$ are defined as in (2) and (3), respectively; see Ramsay and Dalzell (1991).

[TABLE]

Let $\left\{y_{t}(x):t=1,\cdots,N,x\in J\right\}$ 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.

[TABLE]

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 $\phi_{k}$ and their associated scores $\beta_{k,t}$ ($k=1,2,\cdots$). In doing so, it decomposes the covariance operator given in (3) into orthogonal bases of eigenfunctions. Let $\psi_{k}$ and $\lambda_{k}$, respectively, denote the $k$th eigenfunction and eigenvalue. Then, the covariance operator is decomposed as follows:

[TABLE]

The $k$th principal component score is then defined as $\beta_{k,t}=\int y_{t}(x)\psi_{k}(x)dx$. In what follows, the random functions are expressed using Karhunen-Loève expansion as:

[TABLE]

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 $\left\{y_{t}(x_{j})\right\}$, for $t=1,\cdots,N$ and $j=1,\cdots,J$ be the observed functional time series (with noise) characterized by a continuous smooth function $f$ and an error process $\varepsilon$:

[TABLE]

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

[TABLE]

where $\hat{\mu}(x)=N^{-1}\sum_{t=1}^{N}f_{t}(x)$ is the estimated mean of functions, $\hat{\phi}_{k}(x)$ and $\hat{\beta}_{k,t}$ denote the estimated functional principal components and their associated principal component scores, respectively, $\epsilon_{t}(x)$ is the error function with mean zero and $K$ is the number of basis functions with $K<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 $\hat{\beta}_{k,N+h|N}$ the $h$-step-ahead forecast of $\beta_{k,N+h}$ conditionally on the available data up to time $N$. Then, the point predictor of $y_{N+h}(x)$ is obtained by multiplying the forecasted principal component scores with the estimated functional principal components, as follows:

[TABLE]

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-\alpha)100\%$ prediction interval for $y_{N+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:

($i$) smoothing error $\varepsilon_{t,j}$ given in (1),

($ii$) the error caused by the functional principal component decomposition $\epsilon_{t}(x)$ in (6), and ($iii$) the error occurring owing to forecasting principal component scores $\beta_{k,t}$. Let $\xi_{k,h,t}$ denote the $h$-step-ahead forecast error $\hat{\xi}_{k,h,t}=\hat{\beta}_{k,t}-\hat{\beta}_{k,t|t-h}$ ($t=h+1,\cdots,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 $y_{t}(x)=f_{t}(x)+\sigma_{t}(x)\varepsilon_{t}$ for $t=1,\cdots,N$ to obtain smooth functions $f_{t}(x)$, smoothing error $\hat{\varepsilon}_{t}$ and variance component $\hat{\sigma}_{t}(x)$.

• Step 2.

Decompose the smooth functions into $k=1,\cdots,K$ orthonormal functional principal components and associated principal component scores and obtain the fitted functions $\hat{f}_{t}(x)$. Then, calculate the error functions $\hat{\epsilon}_{t}(x)=f_{t}(x)-\hat{f}_{t}(x)$.

• Step 3.

Obtain $h$-step-ahead forecasts of the principal component scores $\beta_{k,t}$, $\hat{\beta}_{k,N+h|N}$ for $k=1,\cdots,K$ and $t=N+1,\cdots,N+h$ using a univariate time series method. In addition, obtain the $h$-step-ahead forecast errors $\hat{\xi}_{k,h,t}$.

• Step 4.

Calculate the future bootstrap values as follows:

[TABLE]

where $\hat{\beta}^{*}_{k,N+h|N}=\hat{\beta}_{k,N+h|N}+\xi_{k,h,*}$, $\xi_{k,h,*}$ is an i.i.d. random sample from $\left\{\xi_{k,h,t}\right\}$, $\hat{\epsilon}^{*}_{N+h|N}(x)$, $\hat{\varepsilon}^{*}_{N+h|N,j}$ and $\hat{\sigma}^{*}_{N+h}(x)$ are random samples with replacement from $\left\{\hat{\epsilon}_{t}(x)\right\}$, $\left\{\hat{\varepsilon}_{t}\right\}$ and $\left\{\hat{\sigma}_{t}(x)\right\}$, respectively.

• Step 5.

Repeat Step 4. $B$ times by drawing random samples of the error terms to obtain $B$ sets of bootstrap replicates of $\hat{y}^{*}_{N+h|N}(x)$, $\left\{\hat{y}^{*,1}_{N+h|N}(x),\cdots,\hat{y}^{*,B}_{N+h|N}(x)\right\}$ for each $h$, where $B$ denotes the number of bootstrap simulations.

The $(1-\alpha)100\%$ bootstrap prediction intervals for $y_{N+h}$ are then obtained by the $\alpha/2$th and $(1-\alpha/2)$th quantiles of the bootstrap replicates $\left\{\hat{y}^{*,1}_{N+h|N}(x),\cdots,\hat{y}^{*,B}_{N+h|N}(x)\right\}$. 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 $\beta_{k,t}$ ($t=1,\cdots,N$) is characterized by a zero mean AR($p$) process, as follows:

[TABLE]

where $\Phi=(\phi_{1},\cdots,\phi_{p})$ denotes the parameter vector and $\zeta_{t}$ is an i.i.d. white noise sequence with mean zero and variance $\sigma^{2}_{\zeta}$. Under the assumption of normality, the probability density function for the model (7) is given by

[TABLE]

where $\zeta_{t}(\Phi)=\beta_{k,t}-\sum_{i=1}^{p}\phi_{i}\beta_{k,t-i}$ and $\zeta_{N}=\left\langle\zeta_{p+1},\cdots,\zeta_{N}\right\}$. Let $\ell\left(\zeta_{t}(\Phi);\sigma^{2}_{\zeta}\right)=\log\mathcal{L}\left(\zeta_{t}(\Phi);\sigma^{2}_{\zeta}\right)$ denote the conditional log-likelihood function. Then, the maximum likelihood estimators of $\Phi$ and $\sigma^{2}_{\zeta}$ conditionally on the first $p$ observations, respectively, are obtained by the solution of the score functions $u_{\Phi}$ and $u_{\sigma_{\zeta}}$, as follows:

[TABLE]

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^{*}\left(\zeta_{t}(\Phi),\hat{F}_{N}(\Phi)\right)$ and $m^{*}\left(\zeta_{t}(\Phi),\sigma^{2}_{\zeta}\right)$ denote a kernel density estimator with bandwidth $h$ based on the empirical distribution function $\hat{F}_{N}(\Phi)$ and the kernel smoothed model density, as follows:

[TABLE]

where $M(\sigma^{2}_{\zeta})$ and $k\left(\zeta_{t}(\Phi);r,g\right)$ represent the normal distribution function with mean zero and variance $\sigma^{2}_{\zeta}$ and a kernel density with bandwidth $g$, respectively. Denote by $\delta_{t}\left(\zeta_{t}(\Phi);M(\sigma^{2}_{\zeta}),\hat{F}_{N}(\Phi)\right)=\frac{f^{*}\left(\zeta_{t}(\Phi),\hat{F}_{N}(\Phi)\right)}{m^{*}\left(\zeta_{t}(\Phi),\sigma^{2}_{\zeta}\right)}$ and $\omega(\delta_{t})=\min\left\{1,\frac{[A(\delta_{t})+1]^{+}}{\delta_{t}+1}\right\}$ the Pearson residual and the weight function, respectively, where $[\cdot]^{+}$ denotes the positive part and $A(\cdot)$ is the residual adjustment function. Note that we consider the Hellinger residual adjustment function of Lindsay (1994), $A(\delta)=2[(\delta+1)^{0.5}-1]$. Based on these definitions, the conditional weighted likelihood of the parameters $\Phi$ and $\sigma_{\zeta}$ are obtained by solving the following estimating equations:

[TABLE]

Now let $\hat{\Phi}^{\omega}=\left(\hat{\phi}^{\omega}_{1},\cdots,\hat{\phi}^{\omega}_{p}\right)$ denote the estimated parameter vector of AR($p$) process using the weighted likelihood methodology. Denote by $\hat{\beta}^{\omega}_{k,N+h|N}=\sum_{i=1}^{p}\hat{\phi}^{\omega}_{i}\hat{\beta}^{\omega}_{k,N+h-i|N}$ the weighted likelihood version of $h$-step ahead forecast of $\beta_{k,N+h}$. Then, the weighted likelihood-based $h$-step-ahead 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, $\hat{\omega}_{t}$, when the model is accurately specified and no outlying observation is present in the data; $\sup_{t}|\hat{\omega}_{t}-1|\xrightarrow{p}0$, see Agostinelli (1998). Following by Agostinelli (2002b), it can be shown that $\left|\sum_{t=1}^{N}\hat{\omega}_{t}\zeta_{t}(\hat{\Phi}_{\omega})^{2}-\sum_{t=1}^{N}\zeta_{t}(\hat{\Phi})^{2}\right|=o_{p}(N)$. This result indicates that

($i$) the weighted likelihood method tends to perform similar to the maximum likelihood 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) $\gamma=[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; $y_{t}(x_{j})=15+\cos(\pi j/4)+N(0,0.15^{2})$ ($t=1,\cdots,N=100$ and $j=1,\cdots,J=12$). The magnitude outlier(s) are generated by contaminating $N\gamma$ randomly selected function(s) by a random function $y^{c}_{t}(x_{j})=|N(0.75,0.15^{2})|$, that is $y_{t}(x_{j})=y_{t}(x_{j})+y^{c}_{t}(x_{j})$. Conversely, the shape outliers are generated from the following process: $15+\sin(\pi j/4)+N(0,0.15^{2})$. Examples of the simulated data (with magnitude and shape outliers) are presented in Figure 1.

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 $\alpha$ 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:

[TABLE]

where $\mathds{1}\left\{\cdot\right\}$ denotes the indicator function and $\hat{y}^{*,\alpha}_{N+h|N}(x_{j})$ is the $\alpha$th quantile of the bootstrap replicates.

The results obtained from the Monte Carlo experiments are reported in Table LABEL:tab:sim. 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.

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}(x_{j})=|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 LABEL:tab:sim2. Compared to the results reported in Table LABEL:tab:sim, Table LABEL:tab:sim2 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.

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 constructed based on the first 91 functions. This figure supports the results presented in Figure 3.

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, 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.