Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
59,930 characters · 8 sections · 73 citation commands
Prediction in locally stationary time series
AMS subject classification: 62M10; 62M20
Keywords and phrases: locally stationary time series, high dimensional auto-covariance, matrices, prediction, local linear regression,
\setcounter{equation}{0}
An important problem in time series analysis is to predict or forecast future observations from a given a stretch of data, say $X_{1}, \ldots , X_{n}$, and numerous authors have worked on this problem. Meanwhile there is a well developed theory for prediction under the assumption of stationary processes [see for example brockwell2002introduction, bickel2011banded, mcmurry2015high among many others]. On the other hand, if data is obtained over a long stretch of time it may be unrealistic to assume that the stochastic structure of a time series is stable. Moreover, in many shorter time series non-stationarity can also be observed and prediction under the assumption of stationarity might be misleading.
A common approach to deal with this problem of non-stationarity is to assume a location scale model with a smoothly changing trend and variance but a stationary error process, say $X_n = \mu (n) + \sigma (n) \varepsilon_n$ [see, for example, van2004forecasting, stuaricua2005nonstationarities, zhao2008confidence, guillaumin2017analysis, das2017predictive]. In this case the trend and variance function can be estimated and prediction can be performed applying methods for stationary data to the standardized residuals. However, there appear also more sophisticated features of non-stationarity in the data, which are not captured by a a simple location scale model, such as time-changing kurtosis or skewness, and the standardized residuals obtained by this procedure may not be stationary.
To address this type of non-stationarity various mathematical concepts modeling a slowly-changing stochastic structure have been developed in the literature [see for example, priestley1988non, dahlhaus1997fitting, nason2000wavelet, zhou2009local or VO2012]. The corresponding stochastic processes are usually called {\it locally stationary} and the problem of predicting future observations in these models is a very challenging one. An early reference is fryzlewicz2003forecasting who considered centered {\it locally stationary wavelet processes}. In this model the sample covariance matrix in the prediction equation is not estimable and the authors proposed an approximation using the (uniquely defined) wavelet spectrum. van2004forecasting considered the prediction problem in a location scale model with a smoothly changing variance and stationary error process. More recent work on forecasting in centered locally stationary time series can be found in roueff2018prediction and kley2019predictive. The first named authors investigated a predictor based on auto-regression of a given order, while kley2019predictive considered predictors in stationary and locally stationary models for (possibly) non-stationary data and selected the “better” prediction among the two estimates. A common feature of most of these methods is that they are all based on auto-regressive fitting.
In the present paper we contribute to this literature and propose an alternative method for prediction in physically dependent locally stationary times series, which does not rely on auto-regressive fitting and is therefore more flexible. To be precise we consider the model
where $\mu$ is a deterministic and smooth mean or trend function on the interval $[0,1]$ and $\{ \epsilon_{i,n} : i =1, \ldots , n \}_{n\in \mathbb{N}} $ is a triangular array modelled by a locally stationary process in the sense of zhou2009local - see Section (ref) for mathematical details. We then estimate the regression function $\mu$ by local linear smoothing and define a banded estimator for the corresponding auto-covariance matrix
from the residuals of the nonparametric fit, where the width of the band increases with the sample size. Banded estimates of auto-covariance matrices have been considered by wu2009banding and mcmurry2010banded for {\bf centered} and {\bf stationary} processes using the fact that in this case the matrix $\Sigma_{n}$ in (ref) is a Toeplitz matrix. Neither of these results is applicable under the assumption of non-stationarity (even if the locally stationary process $\{X_{i,n}\}_{i=1,\ldots,n}$ in (ref) is centered).
In Section (ref) we establish consistency (with respect to the operator norm) of the new covariance operator for locally stationary processes with a time varying mean function. These results are then used in Section (ref) to develop new prediction methods, which - in contrast to the currently available literature - do not use autoregressive fitting. In Section (ref) we investigate the finite sample properties of the estimator of the covariance matrix and compare the new predictor with the currently available methodology. Finally, all proofs of our main theoretical results and technical details can be found in Section (ref).
\setcounter{equation}{0}
Consider the time series model (ref) where $\{ \epsilon_{i,n}: i=1,\ldots,n\}_{n \in \mathbb{N}}$ is an array of centered random variables, and $\mu:[0,1]\rightarrow \mathbb R$ is a smooth mean function. More precisely we assume
In order to model a local stationary error process we use a concept introduced by zhou2009local. To be precise, define for an $L^q$-integrable random variable $X$ its norm by $\| X \|_q=(\mathbb{E} [|X|^q])^{1/q} (q \geq 1)$, let $\{\varepsilon_i : i \in \mathbb{Z} \}$ denote a sequence of independent identically distributed observations and define $\mathcal{F}_i = ( \ldots ,\varepsilon_{i-2},\varepsilon_{i-1},\varepsilon_i )$. We assume that there exists a function $G : [0,1]\times \mathbb R^\mathbb{N} \rightarrow \mathbb R$ such that
is a well defined random variable. For arbitrary functions $G$ it is not guaranteed that the stochastic structure of $\{\epsilon_{i,n} \colon i \in \mathbb Z\} $ varies smoothly, but we can achieve this by the following assumptions.
Next we quantify the dependence structure. For this purpose let $\{\varepsilon_i' : i \in \mathbb{Z} \} $ denote an independent copy of $\{\varepsilon_i : i \in \mathbb{Z}\}$, define $\mathcal{F}_i^* =(\ldots ,\varepsilon_{-2},\varepsilon_{-1},\varepsilon_0',\varepsilon_{1}, \ldots ,\varepsilon_i)$ and
as a measure of dependence. We assume for the same $q\geq 2$ as in assumption (L1) that
If assumption (L1) holds the covariance matrix $\Sigma_n=(\sigma_{i,j,n})_{1\leq i,j\leq n} $ in (ref) is well defined, where
Throughout this paper we do not reflect the dependence on $n$ in the notation of the entries of a matrix, whenever it is clear from the context. For example we will use $\sigma_{i,j}$ instead {of} $\sigma_{i,j,n}$ and similarly a simplified notation for corresponding estimates. We also define the (time dependent) auto-covariances
of the stationary (for fixed $t \in [0,1]$) process $\{G(t,\mathcal{F}_i)\}_{i \in \mathbb{Z}} $. To estimate the covariances in (ref) we use a local linear regression estimate of the function $\gamma_k$. In order to prove consistency of this estimator we require a smoothness condition on the auto-covariances in (ref), which is formulated as follows.
An application of the Cauchy-Schwarz inequality and the dominated convergence theorem show that a sufficient condition for assumptions (L2) and (A1), is given by (L1) and
In the following section we will use the local linear estimates for the function $\gamma_k$ to define a banded estimate of the covariance matrix $\Sigma_n$ of a locally stationary process of the form (ref) and investigate its asymptotic properties for increasing sample size. We also discuss a corresponding estimator in the stationary case because usually estimators are studied under the assumption of a centered stationary process, that is $\mu \equiv 0$. In the subsequent Section (ref) we use these results for prediction in locally stationary processes with a non-vanishing trend.
\setcounter{equation}{0}
The estimation of the covariance matrix has attracted considerable attention in the literature. We refer among many others to the work of bickel2008covariance, bickel2008regularized for high-dimensional independent identically distributed data and anderson2003introduction, wu2009banding, chen2013covariance, box2015time, and mcmurry2015high who considered this problem for time series. Most authors consider the case of a vanishing trend, i.e. $\mu \equiv 0$, and assume that the error process $\{\epsilon_{i,n} : i=1,\ldots,n \}$ is a sequence of independent identical observations or a stationary series. For example, in the case of a stationary centered process wu2009banding proposed the banded estimator
of the matrix $\Sigma_n$, where $\mathbf 1(A) $ denotes the indicator function of the set $A$ and $$ \tilde \sigma_{i,j}=\frac{1}{n-|i-j|}\sum_{s=1}^{n-|i-j|}X_{s,n}X_{s+|i-j|,n}, $$ is the sample auto-covariance of $\{X_{1,n}, \ldots , X_{n,n} \}$ at lag $|i-j|$ and $l_n \in \mathbb{N}$ denotes a tuning parameter satisfying $l_n\rightarrow \infty$, $l_n=o(n)$ as $n \to \infty$. mcmurry2010banded modified this statistic such that the new estimator leaves the band intact, and then gradually down-weighs increasingly distant off-diagonal entries instead of setting them to zero as in the banded matrix case. Both estimators use the fact that for stationary processes the matrix $\Sigma_n$ is a Toeplitz matrix.
Note that the estimator (ref) is not consistent for the auto-covariance if the mean function is not constant. As there are many applications where time series have a smoothly changing mean function we begin our discussion analyzing a mean-corrected estimator of the matrix $\Sigma_n$ for a stationary {error} process of the form (ref), which avoids this problem.
Let $\hat \mu $ be the local linear estimator defined by
where $\tau_n$ denotes the bandwidth. For the kernel $K$ we make the following assumption:
We consider the residuals
obtained from the local linear fit and denote by $$ \hat\sigma^\dag_{i,j} =\frac{1}{n-|i-j|}\sum_{s=1}^{n-|i-j|}\hat \epsilon_{s,n}\hat \epsilon_{s+|i-j|,n} ~~~(i,j=1 , \ldots , n) $$ the sample auto-covariance of the residuals $\{\hat \epsilon_{1,n} , \ldots , \hat \epsilon_{n,n} \}$ at lag $|i-j|$. Finally, we define for $ l_{n} \in \mathbb{N} $ the banded matrix
as an estimator of the matrix $\Sigma_n$. It will be shown below that the estimator $\hat \Sigma^\dag_{n}$ is consistent for $\Sigma_n$ in the case of a strictly stationary error process. To measure the distance between two matrices (of increasing dimension) we introduce the operator norm
of a matrix $A$, where $|\cdot|$ denotes the Euclidean norm (note that $\rho^2(A)$ is the largest eigenvalue of the matrix $A^\top A$).
Theorem (ref) establishes consistency of the estimator of the covariance matrix in model (ref) in the operator norm under the assumption of a stationary error process. However, there also exist many time series exhibiting a non-stationary behaviour in the higher order moments and dependence structure [see stuaricua2005nonstationarities, elsner2008increasing, guillaumin2017analysis among others], and estimation under the assumption of a location model with a stationary error process might be misleading. In this case the estimator $\hat \Sigma^\dag_{n}$ in (ref) is not necessarily consistent since the unknown covariance matrix $\Sigma_n$ is not a Toeplitz matrix. To address this problem we propose an alternative approach which also yields a consistent estimator for non-stationary time series. Roughly speaking, we estimate the elements $\sigma_{i,j} $ in the matrix $\Sigma_n$ by
where $\hat \gamma_k (t) $ is a local linear estimate of the auto-covariance function (ref) of the process $\{G(t,\mathcal{F}_i)\}_{i \in \mathbb{Z}} $.
To be precise, we distinguish between a lag of odd or even order and define
if the lag $k$ is of even order, where $b_n$ is a bandwidth and the residuals $\hat \epsilon_{i,n} $ are defined in (ref). In (ref) we use the notation $\hat \epsilon_{i,n}=0$ if the index $i$ satisfies $i<0$ or $i>n$. Similarly, for an odd lag $k$ we define
where
The estimator of the element $\sigma_{i,j}$ in $\Sigma_n$ is finally defined by (ref) and for the covariance matrix we use again a banded estimator, that is
Our next result yields the consistency of this estimator in the operator norm.
\setcounter{equation}{0}
In this section we discuss some applications of the proposed estimators in the problem to perform predictions in locally stationary processes. For centered time series this problem has been recently investigated by roueff2018prediction, kley2019predictive who proposed to fit a locally stationary AR model and perform the prediction using an AR approximation. In this section, we suggest an alternative method which is not based on AR fitting. To be precise, assume that we observe a stretch of data $X_{1,n}, \ldots, X_{m,n}$ from the model (ref) and that we are interested in a prediction of the next observation $X_{m+1,n}$. To be precise, our aim is the construction of best linear predictor of $X_{m+1,n}$ based on $X_{1,n},\ldots ,X_{m,n}$. For this purpose we define
where $\mathbf X_{m,n}=(1,X_{1,n},...,X_{m,n})^\top$ and the prediction vector $\mathbf a_m=(a_{m+1,n},a_{m,n}...,a_{1,n})^\top:=(a_{m+1,n}, (\mathbf{a_m}^{*} )^{\top} )^{\top}$ is given by
In order to estimate the vector $ \mathbf a_m$ we define the local linear estimators from the sample $X_{1,n}, \ldots, X_{m,n}$ by
and denote by
the covariance matrix of the vector $(X_{1,n}, \ldots, X_{m,n})^T$. The residuals (ref) for estimating the auto-covariances are then replaced by residuals by $$ \hat \epsilon_{i,n}^{1:m} =X_{i,n}-\hat \mu^{(1:m)} (i/n) \quad (i =1, \ldots , m) $$ from the nonparametric fit from the data $X_{1,n}, \ldots, X_{m,n}$. Next, we define $\hat \gamma_k^{1:m}$ as the analogue of the estimator (ref) (if the lag $k$ is even) and (ref) (if the lag is odd), where the residual $\hat \epsilon_{\ell,n}$ is replaced by $\hat \epsilon_{\ell,n}^{1:m}$. We further define
as a banded estimator of the covariance matrix $\Sigma_{n,m}:=\text{Cov}(X_{i,n},X_{j,n})_{1\leq j\leq m}$ in (ref). It can be shown that, if the assumptions of Theorem (ref) are satisfied and $m\geq \lfloor cn\rfloor$ for some positive constant $c$,
where the sequence $r_n$ is defined in (ref). We shall construct a predictor based on $ \hat \Sigma_{n,m}^{-1}$ and for this purpose we show that the consistency of the estimator $\hat \Sigma_{n,m}$ in (ref) can be transferred to its inverse. \\ Throughout this paper we denote $\lambda_{min}(A)$ the minimum eigenvalue of a symmetric matrix $A$ and make the following assumption.
We can now define an estimate $\hat{ \mathbf {a}}_m=(\hat a_{m+1,n},(\hat {\mathbf{a}}_m^*)^\top)^\top$ of the vector $\mathbf a_m$ in (ref) by $$ \hat a_{m+1,n} =\hat \mu^{1:m}(m/n)-\sum_{s=1}^m\hat a_{m+1-s,n}\hat \mu^{1:m}(s/n),\label{hatam}\\ $$ and
where
The final predictor of $X_{m+1,n}$ is defined by
The rate $r_n^\circ$ in (ref) results from convergence rate of the nonparmetric estimate of the time-varying mean and does not appear if the trend is not estimated because it is known to be $0$. Conditions (P2) and (P3) can be verified by checking the coefficients of the MA representation of the locally stationary AR process (ref). They assure that for any $i,j$, the process $\{\mathbb{E}(G(t,\mathcal{F}_i)G(s,\mathcal{F}_j))\}_{t,s\in [0,1]}$ is sufficiently smooth on $[0,1]\times [0,1]$.
\setcounter{equation}{0}
To implement our method we need to choose several tuning parameters: the bandwidths $\tau_n$ and $b_n$ for the local linear estimators of the trend $\mu$ and auto-covariance function $\gamma_k$ and the width $l_n$ of the banded estimator of the covariance matrix $\Sigma_n$. For choosing $\tau_n$, we recommend the Generalized Cross Validation (GCV) method proposed in zhou2010simultaneous.
More precisely, let $\hat \mu^{1:m}(\cdot,\tau), 1\leq i\leq m$ be the local linear estimate of the mean trend defined in (ref) using bandwidth $\tau$, then we choose $\tau_n$ as
where $T^{1:m}_{\tau,ii}$ is the $i_{th}$ diagonal entry of the matrix $$ J^{1:m}_0 \big ( (X^{1:m}(i/n))^\top W^{1:m}_\tau(i/n)X^{1:m}(i/n) \big )^{-1}(X^{1:m} (i/n))^\top W^{1:m}_\tau(i/n), $$ $J^{1:m}_0$ and $X^{1:m}(i/n)$ are $m\times 2$ matrices defined by
respectively, and $W_\tau(x)$ is an $m \times m$ diagonal matrix with elements $\big\{ K \big (\frac{x-s/n}{\tau} \big ) \big \}_{s=1, \ldots m}$. The bandwidth $b_n$ for the estimation of the auto-covariance function $\gamma_k$ in (ref) is defined similarly. For example, if $k$ is even, we choose $b_n$ as
where $\hat {\gamma}_k^{1:m}(i/n,c)$ is the local linear estimator with bandwidth $c$ defined as in (ref) using $m$ observations and $T^{1:m}_{c,ii}$ is defined as in the previous paragraph.
To motivate the choice of the width $l_n$ in the banded estimator of the covariance matrix, note that
[see Section 4.3 in zhang2012inference], where $\tilde \sigma_k^2=\int_0^{\frac{m}{n}\wedge \frac{n-k}{n}}g^2(t)dt$, and the function $g^2$ is the long-run variance of the locally stationary process $\{\epsilon_{i,n}\epsilon_{i+k,n}\}_{i=1}^{n-k}$. For its estimation we use a statistic proposed by dette2018change, which is defined as follows. Consider the partial sum of lag $k$ $$ ^kS^{1:m}_{r_0,r_1}=\sum_{i=r_0}^{r_1}\hat \epsilon^{1:m}_{i,n}\hat \epsilon^{1:m}_{i+k,n}, $$ where we use the notation $\hat \epsilon^{1:m}_{i,n}=0$ if the index $i$ satisfies $i<1$ or $i>m$. For an integer $b\geq 2$ we introduce the quantities $$ ^k\Delta^{1:m}_{j,b}=\frac{^k S^{1:m}_{j-b+1,j}- {^k S^{1:m}_{j+1,j+b}}}{b}. $$ Finally, we define for $t\in[b/n,(m-b)/n]$
where $$ \omega(t,i)=K\Big(\frac{i/n-t}{b_n}\Big){\Big /}\sum_{ {i}=1}^nK\Big(\frac{i/n-t}{b_n}\Big) $$ and the bandwidth $b_n$ is given by (ref) with $\hat \epsilon^{1:m}_{i-k/2,n}\hat \epsilon^{1:m}_{i+k/2,n}$ there replaced by $\hat \epsilon^{1:m}_{i,n}\hat \epsilon^{1:m}_{i+k,n}$. For $t\in[0,b/n)$ and $t\in((m-b)/n,m/n]$ we define $\hat g^2 (t)=\hat g^2(b/n)$ and $\hat g^2 (t)=\hat g^2 ((m-b)/n)$, respectively. Finally, we propose
as a data-driven choice of the width $l_n$, where $ \kappa(\alpha)$ is the $\frac{1+(1-\alpha)^{1/(l_1-l_0+1)}}{2}$-quantile of the standard normal distribution and $l_{0}$ and $l_1$ are constants (if the set $\big\{ n^{-1/2}|\sum_{i=1}^n\hat \epsilon^{1:m}_{i,n}\hat \epsilon^{1:m}_{i+l,n}|\geq \kappa({0.01}) \hat \sigma_l, 1_0\leq l\leq l_1\big\}$ is empty we define $l_n=l_0-1$).
In this section we investigate the finite sample properties of the estimators (ref) and (ref) for the covariance matrix $\Sigma_n$ of a locally stationary process, where we consider
as mean functions. Recalling the notation $\mathcal{F}_i = ( \ldots , \varepsilon_{i-1},\varepsilon_i ) $ we investigate four different distributions for the errors in model (ref):
Note that model (a) defines a stationary process and model (b) defines a locally stationary AR(1) process. Model (c) defines a nonlinear $tvMA(1)$ process. Since the innovations $\varepsilon_i$ in model (c) have a symmetric distribution, the covariance matrix of model (c) is diagonal. Model (d) defines a $tvMA(2)$ process, where only the entries in the diagonal and the first two off diagonals of the covariance matrix do not vanish.
We examine the estimator for covariance matrix $\Sigma_n$ for sample sizes $n=250$, $500$ and $1000$ using $1000$ simulation runs. For the estimation of the width $\l_{n}$ of the band in (ref) we use (ref) with $l_0=1$, $l_1=6$. In each simulation run the tuning parameters ($\tau_n$, $b_n$) are determined as described at the beginning of this section. In Table (ref) and (ref) we display the simulated mean squared error of the spectral loss $\rho (\hat \Sigma_n - \Sigma_n)$ for different estimators $\hat \Sigma_n$, where different mean functions and error processes in model (ref) are considered. In particular we compare the mean corrected estimator (ref) for non-stationary error processes with the mean corrected estimator (ref) which assumes a stationary error process. The numbers in brackets show the standard error of the estimates. We observe that in the stationary model (a) the accuracy of both estimators improve with increasing sample size. Moreover, the estimator (ref) outperforms (ref) because this estimator is constructed for stationary processes. On the other hand, for the dependence structures (b) - (d) corresponding to locally stationary processes the stationary method in (ref) is not consistent and the estimator (ref) shows a substantially superior behaviour.
To illustrate the finite sample properties of the estimator proposed in Section (ref) for prediction we examine the mean trend (ref). As error process we consider a locally stationary AR(6) model defined by
where the functions $a_1(t), \ldots , a_6(t)$ are given by
$ \sigma(t)=(1+0.5\sin 2\pi t)^{0.5}$ and $\mathcal B$ is the lag operator on the filter $\mathcal{F}_i$, i.e., $\mathcal B G(t,\mathcal{F}_i)=G(t,\mathcal{F}_{i-1}).$ We consider a standard normal as well as a $\chi^2(6)$ distribution for the errors $\varepsilon_i$ (centered and standardized such that $\mathbb{E}[\varepsilon_i]=0$ Var$(\varepsilon_i)=1$) and examine the mean squared error of the prediction for sample sizes $n=250, n=500,n=1000$. We also compare the new predictor with the methods in roueff2018prediction, kley2019predictive and giraud2015aggregation which were theoretically investigated for centered data. In a first step we used these methods with the residuals $\hat \epsilon^{1:m}_{i,n}$ to obtain a prediction for the de-trended series. In a second step we add to this estimate the value $\hat \mu^{{1:m} }(m/n)$ to obtain the final prediction of $X_{m+1,n}$. Notice that these authors use time-varying AR$(d)$ processes to approximate the time series for prediction without knowing $d$. Since the error process (ref) is a locally AR$(6)$ process, we investigate the performance of the methods proposed by roueff2018prediction, kley2019predictive and giraud2015aggregation for $d=3$, $d=6$ and $d=9$ (note that in the predictor of kley2019predictive $d$ denotes the maximum lag that their algorithm allows). These cases represent the situation of underestimation, correct-estimation and overestimation of $d$. Note that in the cited references there are no rules how to select $d$. Moreover, for the method proposed by kley2019predictive we choose the parameter $\delta$ in their procedure as $0.05$, as a small parameter $\delta$ prefers the choices of a time-varying model to a stationary model.
In Table (ref) and (ref) we present the simulated mean squared error $$ \mathbb{E}[(\hat X^{pred}_{m+1,n} - X_{m+1,n})^2] $$ for the four different prediction methods and different distributions of the innovations. The columns denoted by $t_{pred}=0.5$ and $t_{pred}=1$ correspond to a prediction of $X_{\lfloor n/2 \rfloor +1}$ from on $X_{1,1}, \ldots, X_{\lfloor n/2 \rfloor }$ and a prediction of $X_{n,n}$ from $X_{1,1},\ldots, X_{n-1,n}$, respectively, where {we use $l_0= \lceil \log (m) \rceil $ and $l_1=5+\lceil \log (m) \rceil$ in (ref).} The first row shows the simulated mean squared error of the prediction (ref). With increasing sample size this mean squared error approximates $1$. This corresponds to our theoretical result in Theorem (ref), because we have for the model under consideration $\sigma(0.5)=\sigma(1)=1$. The rows denoted by R-S, G-R-S and K-P-F show the simulated mean squared error for predictors proposed by roueff2018prediction, giraud2015aggregation and kley2019predictive, respectively, with different time lags $d=3,6,9$. In general, the non-stationary predictor (ref) performs better or similar as the alternative methods with different time lag $d$ in all scenarios. Our simulation results also demonstrate that the performance of R-S, G-R-S and K-P-F predictors depend sensitively on the choice of $d$. Finally, the large numbers in R-S predictor is due to the singularity of estimated local covariance matrix. We expect that this can be corrected by using an eigenvalue corrected positive definite covariance matrix estimator similar to (ref).
We also examine the distribution of the prediction error as investigated in Theorem (ref). For this purpose we show in Figure (ref) the QQ plot of prediction errors of the predictors (ref) for standard normal distributed errors and centered and standardized $\chi^2(6)$-distributed errors {in} model (ref), respectively. The model is given by (ref) and the sample sizes is $n=1000$. These results confirm the theoretical findings in Theorem (ref).
Finally, we compare the new predictor (ref) with the methods proposed by roueff2018prediction, giraud2015aggregation and kley2019predictive in a locally stationary MA(6) model defined by
where the time varying coefficients $a_1, \ldots a_{6}$ and the function $\sigma$ are the same as those defined in the locally stationary AR$(6)$ model (ref), the mean function is given by (ref) and the random variables $\varepsilon_{i}$ are independent standard normal distributed. The results are presented in Table (ref) and we observe similar properties as in the locally stationary AR$(6)$ model (ref). A detailed discussion is omitted for the sake of brevity.
In this section we apply our method to predict market indices. Let $p_t$ be the adjusted daily closing value at day $t$, then the log return $r_t$ is defined as
As pointed out by stuaricua2005nonstationarities, the sign of $r_t$ is unpredictable. As a result, these authors proposed to model $r_t$ as
where $\mu $ and $\sigma $ are time varying functions and $\epsilon_t$ denotes a zero-mean noise process. stuaricua2005nonstationarities used model (ref) to study the non-stationarity of stock returns. In this section we apply the new method to predict $y_t:=\log(|r_t|)$ for the SP500, NASDAQ and Dow Jones Index. We consider data from Dec. $19$, $2016$ to Dec. $17$, $2019$. For SP500, NASDAQ and Dow Jones Index, we delete the log return of Jan. 10, 2017, Nov. 13, 2018 and Nov. 12, 2019 respectively due to their negative infinity values. Therefore the lengths of the series are $752$. We use the new method to predict the market indices at trading days between April. 8, 2019 and Dec. 17, 2019 for SP500 and NASDAQ and at trading days between April. 5, 2019 and Dec. 17, 2019 for Dow Jones Series, respectively, and calculate the empirical mean squared error for these predictions. For the sake of comparison we also apply the methods of roueff2018prediction (R-S), giraud2015aggregation (G-R-S) and kley2019predictive (K-P-F) to the same series. As in the simulation, for fair comparison we perform those algorithms on non-parametrically de-trended data and use the outcome plus $\hat \mu((T-1)/T)$ as the prediction of indices at day $T$. The corresponding results are listed in Table (ref), where we use the different lags $3,6,9$ in the procedures based on autoregressive fitting. We observe that the new prediction method (ref) shows the best performance for all three market indices. For NASDAQ index the method proposed by \ kley2019predictive with $d=9$ shows a similar performance. In general the parameter $d$ for the prediction method proposed by roueff2018prediction, giraud2015aggregation and kley2019predictive is difficult to select, while it has a complicated impact on the predictions when applying those approaches. In Figure (ref) we also plot the prediction error of the different methods for the three market indices. The left panels display $\log|r_t|$, while the right panels show absolute prediction errors of the prediction (ref) and of the predictors proposed by roueff2018prediction (R-S), giraud2015aggregation (G-R-S) and kley2019predictive (K-P-F) for the corresponding parameter $d \in \{ 3,6,9\} $, which achieves the smallest mean squared error.