Adv Data Anal Classif DOI 10.1007/s11634-014-0170-x REGULAR ARTICLE

Modeling and forecasting interval time series with threshold models Paulo M. M. Rodrigues · Nazarii Salish

Received: 27 February 2013 / Revised: 15 March 2014 / Accepted: 20 March 2014 © Springer-Verlag Berlin Heidelberg 2014

Abstract This paper proposes threshold models to analyze and forecast intervalvalued time series. A relatively simple algorithm is proposed to obtain least square estimates of the threshold and slope parameters. The construction of forecasts based on the proposed model and methods for the analysis of their forecast performance are also introduced and discussed, as well as forecasting procedures based on the combination of different models. To illustrate the usefulness of the proposed methods, an empirical application on a weekly sample of S&P500 index returns is provided. The results obtained are encouraging and compare very favorably to available procedures. Keywords Interval-valued data · Time series · Nonlinearities · Threshold models · Combined forecasts · S&P500 index returns Mathematics Subject Classification (2000)

62F10 · 62P20

1 Introduction Modeling and forecasting interval-valued time series (ITS, hereafter) has received considerable attention in recent literature. The possibility to analyze large data sets (such as, e.g., high frequency data) based on the ITS framework provides a new domain for the analysis of statistically detectable patterns and has attracted the attention of a

P. M. M. Rodrigues Economic Research Department, Banco de Portugal, Av. Almirante Reis 71-6th Floor, 1150-012 Lisbon, Portugal e-mail: [email protected] N. Salish (B) Bonn Graduate School of Economics, University of Bonn, Kaiserstrasse 1, 53113 Bonn, Germany e-mail: [email protected]

123

P. M. M. Rodrigues, N. Salish

large number of researchers in economics and finance. The core of available literature on ITS is based on the extension of classical data analysis and statistical methods. In particular, exponential smoothing, pattern recognition and multivariate models are among the most influential methodologies used (see Sect. 3 for further references). The main interest of this paper lies in the possible nonlinear structure of ITS. Hence, we extend the available literature by proposing nonlinear models for the ITS context. The motivation for this contribution lies in the remarkable amount of empirical evidence in the conventional time series literature suggesting that nonlinear models do provide a richer understanding of the dynamics of variables of interest. For instance, Henry et al. (2001) present evidence of threshold nonlinearity in the Australian real exchange rate, Dueker et al. (2007) develop a contemporaneous threshold autoregressive model for the bonds market, and Guidolin et al. (2009) studied extensively nonlinear predictability of stock and bond returns. In the context of ITS, several studies have also stressed the importance of nonlinearities (see e.g., Muñoz et al. 2007; Maia and de Carvalho 2011). However, these procedures aim to produce forecasts without explicitly modeling the nonlinear characteristics of the data and to the best of our knowledge there have been no studies in the ITS context that attempt to model and explain nonlinear features. Furthermore, no empirical results on regime dependent ITS forecasts are available so far. In this paper two contributions are made. First, econometric methods for regime switching threshold models are adapted to ITS characterized by their center and radius and least squares techniques are used to estimate them. The estimation algorithm is simple to implement and involves a joint grid search over the vector of thresholds for center and radius. Second, we discuss forecasting techniques for nonlinear ITS. The main focus is on one step ahead forecasts from the proposed threshold model. The development of multi step ahead forecasting techniques for threshold models is far from being trivial (see, e.g., Clements and Smith 1999) and is beyond the scope of this paper. Due to the fact that data may have different time series patterns, not necessarily purely linear or nonlinear, we stress the importance of forecast combination methods for practical applications. Moreover, it has been found in empirical studies that forecast combination methods produce on average better forecasts than methods based on ex-ante best individual forecasting methods (see, e.g., Timmermann 2006). Maia and de Carvalho (2011) motivated combination methods in ITS by combining neural network models and Holt’s exponential smoothing using the hybrid approach of Zhang (2003). Therefore, we propose procedures for forecast combinations and corrections for the ITS context. Finally, in order to provide exhaustive information for an objective decision regarding the interval forecast performance we also discuss quality measures designed for the interval framework (such as center-radius or lowerupper bound); as well as additional forecast descriptive statistics (such as efficiency and coverage rates). This paper is structured as follows. Section 2 introduces the theoretical framework of the interval threshold model, the estimation procedure and discusses related computational issues. Section 3 provides a brief overview of available ITS forecasting methods, discusses accuracy measures and comparison procedures for forecast evaluation, and presents new forecasting approaches. Section 4 discusses an empirical application on the S&P500 index returns and Sect. 5 concludes the paper.

123

Modeling and forecasting interval time series

2 Threshold models for ITS 2.1 The model T , with Yt = 0 for t ≤ 0, such that Consider the observed interval time series {Yt }t=1 each observation Yt , t = 1, . . . , T has an interval structure, i.e. Yt = [at , bt ], and is represented by its center ct = (at + bt )/2 and its radius rt = (bt − at )/2; we briefly write Yt = ct , rt . In our analysis the center-radius representation of intervalvalued data is favored since it is known from literature that the range (or double radius of an interval) is a suitable measure of the variability of a random variable. Several empirical applications, in the context of financial data, use the range for the estimation of volatility (see, e.g., Beckers 1983; Chou 2005). Moreover, the choice of this representation is also based on the available evidence of the good predictive power of the interval radius (see, e.g., Arroyo et al. 2011) and its superiority when compared to the lower-upper bound representation of ITS.1 The model considered in this paper is a center-radius self-exciting threshold model (CR-SETAR, hereafter) for ITS. For the sake of simplicity we focus only on a two regime CR-SETAR model, but this is with no loss of generality since the methods presented can be extended to higher-order threshold models in a straightforward manner as shown in Sect. 2.4. The general two-regime CR-SETAR model that we propose takes the form,

⎛ ct = ⎝μ1 + ⎛

pc 

α1i ct−i +

i=1

⎞

pr 

β1 j rt− j ⎠ × 1{zC,t−1 ≤γC }

j=1

+ ⎝μ2 +

pc 

α2i ct−i +

i=1

⎛ rt = ⎝λ1 + ⎛

qc 

j=1

ϕ1i ct−i +

i=1

+ ⎝λ2 +

pr 

qr 

i=1

β2 j rt− j ⎠ × 1{zC,t−1 >γC } + εt ,

(1)

⎞

φ1 j rt− j ⎠ × 1{z R,t−1 ≤γ R }

j=1 qc 

⎞

ϕ2i ct−i +

qr 

⎞ φ2 j rt− j ⎠ × 1{z R,t−1 >γ R } + ξt ,

(2)

j=1

where 1{·} is an indicator function, and z C,t−1 and z R,t−1 are known functions of the past history of Yt . That is, z C,t−1 and z R,t−1 are measurable w.r.t. Ft−1 , where Ft−1 is the σ -algebra generated by (Yt−1 , . . . , Y1 ). As an evident example one can consider for instance lagged values of the center ct−d and radius rt−d as threshold variables. The parameters pc , pr , qc , and qr represent the different lag orders of the components entering (1) and (2), and γC and γ R denote the threshold parameters 1 Note that it is not clear if the results obtained in the following sections will hold for the lower-upper bound representation of ITS, and we do not pursue this extension in this paper. However, this presents an interesting line of research for further investigation.

123

P. M. M. Rodrigues, N. Salish

which take values in bounded subsets of R: C = [γ C , γ C ] and  R = [γ R , γ R ], respectively. The errors εt and ξt are assumed to have finite second moments, to be independent identically distributed over time and contemporaneously correlated. Therefore, we assume that the  respective form  of the covariance matrix of the distur  σC2 σC R ⊗ IT = u ⊗ IT , where u = ε , ξ  , bances εt and ξt is V ar [u] = 2 σC R σ R IT denotes an identity matrix of dimension T , and ⊗ is the Kronecker product. For notational convenience we consider θ1C := (μ1 , α11 , . . . , α1 pc , β11 , . . . , β1 pr ), θ1R := (λ1 , ϕ11 , . . . , ϕ1qc , φ11 , . . . , φ1qr ), θ2C := (μ2 , α21 , . . . , α2 pc , β21 , . . . , β2 pr ) and θ2R := (λ2 , ϕ21 , . . . , ϕ2qc , φ21 , . . . , φ2qr ). Notice that the threshold model in (1)– (2) has two regimes in each equation, which are governed by θ1C and θ2C for Eq. (1) and θ1R and θ2R for Eq. (2). This model allows all coefficients to switch between regimes, depending on wether z C,t−1 and z R,t−1 exceed the thresholds γC and γ R , respectively. The vector representation of (1)–(2) be considered for the derivation of the esti   will  , θ  , θ   , y := (c , . . . , c , r , . . . , r ) , , θ2C mators. Thus, considering θ := θ1C 1 T 1 T 1R 2R X C,t−1 := (1, ct−1 , . . . , ct− pc , rt−1 , . . . , rt− pr ) , X R,t−1 := (1, ct−1 , . . . , ct−qr,

  rt−1 , . . . , rt−qr ) , X C,t−1 (γC ) := X C,t−1 × 1{zC,t−1 ≤γC } ; X C,t−1 × 1{zC,t−1 >γC } , 

   X R,t−1 (γ R ) := X R,t−1 × 1{z R,t−1 ≤γ R } ; X R,t−1 × 1{z R,t−1 >γ R } , X C (γC ) := (X C,0

     it fol(γC ) , . . . , X C,T −1 (γC )) , and X R (γ R ) := X R,0 (γ R ) , . . . , X R,T −1 (γ R ) lows that (1)–(2) can equivalently be written as

y = X (γC , γ R ) θ + u,

(3)

where X (γC , γ R ) is a block diagonal matrix with X C (γ R ) and X R (γ R ) on the main diagonal. 2.2 Estimation of CR-SETAR models The parameters of interest are θ and γ := (γC , γ R ). Since model (3) is nonlinear in parameters and discontinuous, the sequential conditional least squares principle (see, e.g., Hansen 1996) is used to obtain the parameter estimates. For any given value of γ the multivariate least-squares estimator of θ is the one that minimizes u (γ ) ( u (γ ) ⊗ IT )−1 u (γ ). Hence,  −1      u (γ ) ⊗ IT −1 X (γ ) u (γ ) ⊗ IT −1 y . θ (γ ) = X (γ ) X (γ )

(4)

From the classical VAR theory it is well known that a plausible estimator of u (γ ) is u (γ ) =

 T  1  ξt (γ )  εt (γ )   εt2 (γ ) ,  ξt (γ ) ξt2 (γ )  εt (γ )  T t=1

where  εt (γ ) and  ξt (γ ) are obtained from Eqs. (1) and (2), respectively.

123

Modeling and forecasting interval time series

However, typically γ is unknown and in order to obtain the model parameter estimates one requires an estimate of this vector of threshold parameters. Considering that γ is the value that optimizes the fit of the model, its estimation can be carried out through the minimization of the quadratic form,

u (γ ) ⊗ IT −1  u(γ ) u(γ ), γ = arg min   γ ∈

(5)

where  :=  C ×  R is the

 set of all possible values of threshold parameters (γC , γ R ) γ is obtained, and  u(γ ) =  ε(γ ),  ξ (γ ) are the residuals obtained from (3). Once  the slope estimates  θ :=  θ ( γ ), the least-squares residuals  u :=  u ( γ ) and the sample u ( u = variance matrix γ ) can also be computed. Minimization of (5) can be performed by direct search. Since the quadratic form given in (5) depends on γ through the indicator functions 1{zC,t−1 ≤γC } and 1{z R,t−1 ≤γ R } , it is a step function with at most T 2 steps occurring at distinct values of the observed threshold variables z C,t−1 and z R,t−1 . Thus, the minimization problem in (5) can over values of γ equalling the distinct values of z t−1 := be reduced to searching z C,t−1 , z R,t−1 in the sample, i.e.,

u (z t−1 ) ⊗ IT −1  γ = arg min   u(z t−1 ) u(z t−1 ). z t−1 ∈

(6)

It is useful to note two related issues with the computation of (6). First, thresholds γ which allocate a small number of observations into one regime are undesirable. This possibility can be excluded by restricting the search in (6) to values of γ such that a minimum number of observations in each regime is ensured (see, e.g., Hansen 1996, for more details). Second, to make computation less intensive for applications with large T , instead of searching over all values of z C,t−1 and z R,t−1 the search may be limited to specific quantities which will greatly reduce the search and yield nearly identical results. 2.3 Model diagnostics In this section we stress issues related with CR-SETAR model diagnostics. In particular, lag selection, residual serial correlation and testing for threshold effects are discussed. Model selection criteria for conventional SETAR models, such as AIC, BIC or HQ, asymptotically tend to select an incorrect lag length as indicated, for instance, by Pitarakis (2006). Interesting contributions have already been provided for conventional threshold models by Gonzalo and Pitarakis (2002) and Nieto (2005) among others. However, it is still an open challenge for future research on nonlinear parametric ITS models which requires a careful theoretical treatment and development of new methods. For an empirical application, given the lack of available lag order selection tools in the symbolic data literature, we suggest the use of the autocorrelation and partial autocorrelation functions.

123

P. M. M. Rodrigues, N. Salish

A further important step in the statistical verification of models is to perform diagnostic checks of the residuals. Some tests that are used in the traditional linear framework can also be applied for testing in the nonlinear framework. For instance, Jarque– Bera tests can be used for normality testing in both frameworks. On the other hand, the common Ljung–Box test is not usable because its asymptotic null distribution is unknown when the test is based on estimated residuals from TAR type models. Regarding testing for serial correlation of the errors, convenient results for the threshold context are obtained by Eitrheim and Teräsvirta (1996). They developed a Lagrange multiplier type test for nonlinear time series models (in particular, TAR-type models), which resembles closely the Breusch–Godfrey test for serial correlation. This procedure is employed in the empirical part of the paper. Notice that our main aim is in the investigation of forecasting abilities of the proposed threshold model. However, if one is also interested in the statistical justification of nonlinear versus linear ITS models, appropriate test statistics are required. For the completeness of the discussion we suggest using an adaptation of the well-established techniques for classical multivariate threshold models, such as for instance by Tsay (1998) or more recently by Nieto (2005). The former is also illustrated in Sect. 4.2. For a review of related topics on threshold type models see also Tong (2011). 2.4 Multiple regimes The CR-SETAR model proposed in (3) has a single threshold vector given by γ . However, in some applications multiple thresholds can be considered. For example, model (3) with three regimes will have the following regressor structure    × 1{zC,t−1 ≤γ1C } ; X C,t−1 × 1{γ1C
(7)

(8)

where the thresholds are ordered such that γ1C < γ2C and γ1R < γ2R . For the sake of space and without loss of generality we will focus on this three regime model since the method extends in a straightforward manner to higher-order threshold models. θ (γ ) is easy to compute For a given γ := (γ1C , γ2C , γ1R , γ2R ) the estimator of  using (4), where the regressors have the three regime form given in (7) and (8). The threshold parameter estimates are by definition the values that jointly minimize

u (γ ) ⊗ IT −1  u(γ ). Thus, this quadratic form can be estimated as in the sin u(γ ) gle threshold model. Note that grid search over γ may become quite cumbersome to implement in practice. For a general m-regime model, where m ≥ 2 it requires the estimation of approximately T 2m regressions. However, to avoid this computational burden, in addition to a search over specific quantities of threshold variables, sequential estimation can be applied to the threshold models as suggested, for instance, by Bai (1997) for change point models.

123

Modeling and forecasting interval time series

3 Forecasting ITS A number of methods have been proposed for forecasting in ITS. For instance, smoothing techniques such as exponential and Holt’s smoothing methods (see e.g., Arroyo et al. 2011; Maia and de Carvalho 2011); nonparametric pattern recognition methods, such as the k-nearest neighbors (k-NN) and the multi-layer perceptions (MLP) (see e.g., Arroyo et al. 2010; Muñoz et al. 2007); and multivariate models such as ARIMA, VAR and VEC models (see, e.g., Cheung 2007; Hu and He 2007) have been considered. However, since real world systems are often nonlinear (see, e.g., Granger and Terasvirta 1993) this paper introduces threshold models as a mechanism capable of explaining nonlinearities in interval valued data and to forecast these series. In this section, we focus our attention on one step ahead forecasts constructed from CR-SETAR models. Moreover, two related issues with forecasting ITS are stressed. First, given that one of our objectives is to discuss potential benefits of ITS threshold models in forecasting when compared to existing models, we discuss a wide array of alternative performance measures and procedures to evaluate predictive accuracy across pairs of different ITS models. Second, it is consensual in the forecasting literature that no single method is superior in every situation (see, e.g., Zou and Yang 2004). Moreover, it is observed that in many cases a combination of time series models achieves better results than individual models. Thus, we also implement a mechanism to forecast ITS based on the combination of different models using conventional forecast combination techniques and the hybrid approach of Zhang (2003). 3.1 Measuring forecast accuracy For comparative and accuracy analysis of interval forecasts one can proceed with standard difference measures of forecast accuracy by applying them separately to the descriptive components of interval-valued data (i.e., center and radius or lower and upper bounds). These measures are the root mean square forecast error (RMSFE), the mean absolute forecast error (MAFE), the forecast error bias (FEB), the forecast error variance (FEV) and the mean absolute percent forecast error (MAPFE). However, the ranking of forecasting models for ITS based on these measures will not be exhaustive. The fact that data has an interval structure implies that both characteristics (center and radius) that describe intervals have to be taken into consideration jointly. Hence, in order to quantify the overall accuracy of the fitted or forecasted ITS, the mean distance error (MDE) of intervals is used, i.e.,  t } = MDEd {Yt }, {Y

H



t=1 d(Yt , Yt )

H

,

(9)

t } is its forecast with t = 1, . . . , H . The assesswhere {Yt } is the observed ITS and {Y ment of dissimilarities between the ITS proceeds through the choice of a particular t ). Arroyo and Maté (2006) offer a set of accuracy measures based loss function d(Yt , Y on distances for interval data. In what follows, two alternative loss functions are used. The first is based on a generalization of the Minkowski-type distance of order ρ ≥ 1 defined as,

123

P. M. M. Rodrigues, N. Salish

t ) = d L ρ (Yt , Y

 ρ |ct −  ct |ρ + |rt −  rt |ρ ,

(10)

t =  where Yt = ct , rt  and Y ct , rt ; and the second on the normalized symmetric difference (NSD hereafter) of intervals, t ) = d N S D (Yt , Y

  t − w Yt ∩ Y t w Yt ∪ Y , t ) w(Yt ∪ Y

(11)

where w(·) indicates the width of the interval. Several modifications of (11) appear in related literature, each of which with advantages and disadvantages. For instance, Hsu and Wu (2008) considered the d N S D loss function normalized to the width of the realized interval Yt ; Ichino and Yaguchi (1994) proposed a loss function for a multidimensional space of mixed variables, where the numerator of (11) has an additional term that controls the effect of the inner-bound and outer-bound nearness between t and depends on the supplementary parameter. For convenience of handling Yt and Y ) we suggest a modification which presents useful features. It aims at d N S D (Yt , Y correcting the width of the symmetric difference in (11) with a term that provides information distance between two intervals with empty intersection, on the shortest  t 1Y ∩Y =∅ , where ∪ denotes the interval hull. Thus, using the t / Yt ∪ Y w Yt ∪ Y t t properties of set operations the new NSD distance for intervals can be rewritten as t ) = 2 − d N∗ S D (Yt , Y

t ) w(Yt ) + w(Y . t ) w(Yt ∪ Y

(12)

The usefulness of the metric proposed in (12) results from its scale-independence and simplicity of application. From (12) we observe that the range of this statistic t ) = 2 if and only if the considered is [0, 2]. In particular, note that (i) d N∗ S D (Yt , Y t ) < 2 if the intervals are degenerated (i.e., rt = 0) and not equal; (ii) 1 < d N∗ S D (Yt , Y ∗ t ) = 1 for any intersection of the non degenerated intervals is empty; (iii) d N S D (Yt , Y pair of tangent intervals or if one of the intervals is degenerated and is contained in t ) < 1 if and only if Yt ∩ Y t = ∅. the second one; and (iv) d N∗ S D (Yt , Y Furthermore, we also introduce a set of descriptive statistics defined for ITS. One is the coverage rate H t ) 1  w(Yt ∩ Y , (13) RC = H w(Yt ) t=1

and the other the efficiency rate RE =

H t ) 1  w(Yt ∩ Y . t ) H w(Y

(14)

t=1

The main purpose of considering RC , and R E is to provide additional information on what part of the realized ITS is covered by its forecast (i.e., coverage rate) and what part of the forecast covers the realized ITS (i.e., efficiency rate). However, it should be noted that these statistics need to be considered jointly, given that if analyzed

123

Modeling and forecasting interval time series

individually wrong conclusions may be drawn. For instance, if the realized intervals are fully enclosed in the predicted intervals then the coverage rate will be 100 %, but the efficiency rate may be less than 100 % and reveal the fact that the forecasted ITS is actually wider than the realized ITS. Note that the reverse may also be observed. Hence, only when the coverage and efficiency rates are reasonably high and the difference between them is small will these provide indication of a good forecast. 3.2 Forecast combinations and corrections Consider the problem of forecasting at time T the future value of the target ITS T {Yt }t=1 . The information set available at time T will be denoted by FT . The one T +1 ≡ E[YT +1 |FT ], step ahead forecasts from the CR-SETAR model, defined as Y T +1 =<  cT +1 , r T +1 >=< can be obtained from (3) in a straightforward manner, Y γ ) θC , X R,t ( γ ) θ R >. If forecasts from other N models are available at time T let X C,t ( T +1,1 , . . . , Y T +1,N ) denote an N × 2 vector of interval forecasts. vector  YT +1 := (Y The general forecast combination problem uses the information in a potentially N  high-dimensional vector of forecasts  YT +1 ∈ R × R+ to construct a lower dimen  K  sional summary measure M  YT +1 , w ∈ R × R+ , where K < N and w are the parameters associated with the combination. The vast majority of studies on forecasting has dealt with point interval forecasts. In our case, one dimensional aggregation will suffice and therefore we consider,  T +1 = M  YT +1 , w ∈ R × R+ , Y

(15)

which is the combined interval forecast as a function of the underlying forecasts  YT +1 and the combination parameters w ∈ W , where W is often assumed to be a compact subset in R N . The variables w are thought of as combination weights (although the parameters need not always have this interpretation). For example, equal weights would give N  1  (16) M  YT +1 , w = YT +1,k . N k=1

In general, the optimal combination parameters, w, are not known but can be obtained from the minimization of the combination forecast error using one of the loss functions presented in Sect. 3.1. For an extensive review of the theory underlying the general forecast combination problem see, for instance, Timmermann (2006). An example with equal weight is provided in the following section. An interesting approach to forecast combination is introduced by Zhang (2003). The main idea of the given method is to explore different aspects of the linear and nonlinear patterns in the data using a hybrid methodology. This methodology has been successfully applied to forecast ITS using Holt’s model and neural networks. However, this method differs from the one presented in (16). The algorithm aims to model first the linear structure and then to model the resulting residuals with a nonlinear model. Finally, forecasts are computed from the linear predictions and corrected with the forecasts of the residuals obtained from the nonlinear model. We will refer to

123

P. M. M. Rodrigues, N. Salish

this method as the forecast correction approach while (16) will be referred to as the forecast combination method. Moreover, we suggest using the hybrid procedure in Zhang (2003) not only as a forecast correction procedure for nonlinear models but also as a general procedure for forecast correction of different methods.

4 Empirical application on S&P500 index returns To illustrate the methods previously introduced and their potential usefulness, in this section we provide an empirical application to S&P500 index returns. The objective of the analysis given below is twofold. One is to illustrate the potential of the CRSETAR model in explaining and forecasting ITS and the other to compare the proposed procedures with other methods available in the literature. The last aspect is covered by examining the comparative forecast performance of the CR-SETAR with that of linear models for ITS, the k-NN approach and their combination and correction as suggested in Sect. 3.2. We also supplement our comparative analysis with a number of standard benchmarks commonly employed in empirical finance and forecasting literature (see, for instance, Guidolin et al. 2009) which are adapted for interval-valued data. The first benchmark considers that the S&P500 Index follows a random walk with drift so that, ct+h = μ + εt+h rt+h = λ + ξt+h ,

(17)

where interval is the sample mean of the returns computed at time t, i.e., the forecast

t+h = E Yt+h |Yt , ... = [μ − λ, μ + λ]. The second is a naïve forecast given as Y Yt+h−1 , and the third benchmark is a simple autoregressive model, i.e., ct+h = α + βct + εt+h rt+h = ϕ + φrt + ξt+h .

(18)

Our comparison is based on the methods presented in Sect. 3.1. It includes ITS accuracy measures, and standard accuracy measures used separately for the descriptive components of ITS (i.e., center and radius). Note that in the analysis, the Minkowskitype loss function is used with ρ = 1 and ρ = 2 due to its intuitive and mathematical appeal.

4.1 The data Our sample consists of weekly S&P 500 interval returns observed from January 3, 1997 to November 12, 2010. This sample is divided into two parts: the first sub-sample, from January 3, 1997 to December 28, 2007 (574 weeks), is used for estimation purposes and the second sub-sample, from January 4, 2008 to November 12, 2010 (150 weeks), to compute the ITS forecasts and analyze their quality. Note that the forecast period includes the financial crisis that started in 2008.

123

Modeling and forecasting interval time series

Fig. 1 Lower-upper limits of expanded and centered returns

Since the objective of this analysis is to characterize the evolution of returns, the starting point for building our data set is the adequate definition of interval returns. In what follows, two types of interval returns are considered: (i) The Expanded returns which consist of intervals that contain all possible differences between index price levels that occur in two neighboring trading periods, i.e., (19) YtE = {yt − yt−1 : yt ∈ [lt , u t ] ∧ yt−1 ∈ [lt−1 , u t−1 ]}, where [lt , u t ] is an interval of logarithmic prices. Note that, (19) satisfies the definition of interval subtraction, thus YtE = Yt − Yt−1 . (ii) The Centered returns which are intervals of returns for the current period based on the center of the price index levels for the previous period, i.e. YtC = {yt − yC,t−1 : yt ∈ [lt , u t ]},

(20)

where yC,t−1 is the center of Yt−1 . Figure 1 illustrates the different representations of the S&P 500 returns. Here, the expanded returns describe the range of all possible returns based on two trading periods, while centered returns represent the range of returns assuming that the reference is the centre of its range. Note that the following relation between YtE and YtC can be considered. Considering YtE =< ctE , rtE > and YtC =< ctC , rtC >, then intervals of expanded and centered returns have the same center (i.e., ctE = ctC ) and the radius of the expanded returns is the aggregation C ). of the contemporaneous and lagged centered returns radius (i.e., rtE = rtC + rt−1 4.2 The CR-SETAR model Due to the relation between the two types of interval returns defined in (19) and (20) in what follows we will only consider the CR-SETAR model for centered returns since for the expanded returns it will follow along similar lines. Note that four lags of ct and rt were considered in both equations of model (1)–(2) as this appeared, based on the behaviour of the autocorrelation function, to be the minimum necessary to adequately describe the short-run dynamics. In order to analyze the presence of nonlinearities and threshold effects in the S&P500 returns we considered a standard lag of the radius as the threshold vari-

123

P. M. M. Rodrigues, N. Salish Table 1 Estimation results for the center Eq. (1) Regime 1

Regime 2

Coeff.

Stand. Error

Coeff.

Abs. Changes Stand. Error

Constant

0.0010

0.0029

−0.0230

0.0043

0.0239

ct−1

0.1209

0.0658

0.2078

0.0583

0.0869

ct−2

0.0556

0.0602

0.1039

0.0707

0.0484

ct−3

−0.0618

0.0603

0.1643

0.0638

0.2261

ct−4

0.0094

0.0562

−0.0766

0.0606

0.0860

rt−1

0.1561

0.2494

0.7236

0.1369

0.5675

rt−2

−0.2089

0.1211

0.7808

0.1856

0.9897

rt−3

−0.2504

0.1566

−0.0509

0.1313

0.1995

rt−4

0.3041

0.1569

−0.5157

0.1315

0.8199

γC

0.0204

N

419

155

Regime 1 and Regime 2 is specified as rt−1 ≤ 0.0204 and rt−1 > 0.0204, respectively

able, i.e., z C,t−1 = rt−d and z R,t−1 = rt−l for 1 ≤ d, l ≤ 4, but we also investigated alternative threshold variable candidates, such as lags of the centered returns, ct−d , and different stationary transformations of the YtC (e.g., combinations such as z C,t−1 = rt−d and z R,t−1 = ct−l ). Selection of the optimal threshold variable is done by using the procedure suggested in Hansen (1997) for conventional SETAR models. We report that the results when z C,t−1 = rt−1 and z R,t−1 = ct−1 are used as threshold variables provide the best model fit in terms of the sum of squared residuals obtained from the minimization procedure described in (4) and (5). As suggested in Sect. 2.3 the presence of the threshold effects in the given setup can be formally justified by using the two step conditional least squares procedure proposed by Tsay (1998). Since the structure of the CR-SETAR models allows for different threshold variables in both equations the Tsay-test is performed for each equation. Results indicate a clear rejection in the center equation with related p-value equal to 6.44 × 10−8 and in the radius equation with p-value equal to 0.022. For the selected CR-SETAR model setup Tables 1 and 2 present the estimation outputs. We report parameter estimates, standard errors, absolute parameter changes over regimes, the number of observations in each regime and the threshold estimates. In addition, we found that the residuals for the given CR-SETAR process are not serially correlated, if we look at the LM test for serial correlation (see, Eitrheim and Ter¨asvirta, 1996) with up to 4 lags (producing p-values equal to 0.993 for the residuals from the center equation and 0.823 for the residuals from the radius equation, respectively).

4.3 Forecasting performance Table 3 displays the performance of the forecasts obtained from the different approaches. We compare the forecasts of the random walk with drift (M1), the naïve

123

Modeling and forecasting interval time series Table 2 Estimation results for the radius Eq. (2) Regime 1

Regime 2

Abs. Changes

Coeff.

Stand. Error

Coeff.

Stand. Error

Constant

0.0063

0.0065

0.0038

0.0009

0.0025

ct−1

−0.1749

0.1479

−0.0659

0.0236

0.1090

ct−2

−0.0306

0.0946

−0.0200

0.0215

0.0106

ct−3

−0.0485

0.0885

−0.0160

0.0210

0.0325

ct−4

−0.0095

0.0868

−0.0054

0.0194

0.0040

rt−1

−0.0977

0.0906

0.3894

0.0603

0.4871

rt−2

0.0928

0.1763

0.1487

0.0477

0.0559

rt−3

0.4052

0.2912

0.2543

0.0470

0.1509

rt−4

0.1859

0.2629

0.0038

0.0472

0.1821

γR N

−0.0277 540

34

Regime 1 and Regime 2 is specified as ct−1 ≤ −0.0277 and ct−1 > −0.0277, respectively Table 3 Overview of forecasting performance M1

M2

M3

M4

M5

M6

M7

M8

M9

M10

Centered returns MDE (L1)

0.032 0.035

0.028 0.027

0.029

0.029

0.024

0.027

MDE (L2)

0.036 0.041

0.032 0.032

0.033

0.033

0.028

0.032

0.024 0.025 0.028 0.029

MDE (NSD)

0.611 0.624

0.556 0.544

0.557

0.570

0.502

0.550

0.504 0.516

Efficiency

0.600 0.514

0.611 0.602

0.600

0.620

0.653

0.596

0.652 0.639

Coverage

0.512 0.525

0.572 0.595

0.581

0.537

0.638

0.586

0.636 0.619

MDE (L1)

0.042 0.035

0.033 0.031

0.035

0.036

0.027

0.031

0.027 0.028

MDE (L2)

0.046 0.041

0.036 0.035

0.039

0.039

0.031

0.035

0.031 0.032

MDE (NSD)

0.434 0.387

0.354 0.351

0.371

0.380

0.313

0.349

0.314 0.320

Expanded returns

Efficiency

0.817 0.750

0.817 0.798

0.794

0.829

0.829

0.799

0.828 0.826

Coverage

0.679 0.751

0.756 0.777

0.763

0.717

0.801

0.777

0.800 0.790

Radius of ITS RMSFE

0.020 0.014

0.014 0.012

0.015

0.016

0.010

0.012

0.010 0.010

Bias

0.008 0.000

0.004 0.002

0.004

0.006

0.001

0.002

0.001 0.001

FEV

0.000 0.000

0.000 0.000

0.000

0.000

0.000

0.000

0.000 0.000

MAFE

0.012 0.009

0.008 0.007

0.008

0.009

0.006

0.007

0.006 0.007

MAPFE

0.392 0.367

0.312 0.293

0.324

0.302

0.262

0.297

0.265 0.275

0.030 0.038

0.030 0.029

0.031

0.030

0.027

0.030

0.027 0.028

Center of ITS RMSFE Bias

−0.003 0.000 −0.002 0.023 −0.002 −0.002

−0.001 −0.003 −0.001 0.000

FEV

0.001 0.001

0.001 0.001

0.001

0.001

0.001

0.001

0.001 0.001

MAFE

0.021 0.026

0.020 0.023

0.021

0.020

0.018

0.020

0.018 0.019

MAPFE

1.445 6.879

1.795 0.886

3.022

3.260

1.891

2.778

1.940 3.226

123

P. M. M. Rodrigues, N. Salish

model (M2), the simple autoregressive framework (M3), a linear VAR model for center and radius (M4), the k-NN method based on NSD loss function (M5), the k-NN method based on L 2 (M6), the proposed CR-SETAR (M7), the VAR model corrected with the k-NN approach (M8), the CR-SETAR model corrected using a VAR (M9), and the CR-SETAR model corrected using the k-NN approach (M10). The upper part of Table 3 corresponds to the forecasting performance of each model for centered and expanded S&P500 returns. The bottom part of the table shows the results of the different approaches to forecast the ITS based on the performance of individual components (i.e., center and radius). We also investigate results for model combination in pairs using equal weights as in (16). However, for the sake of space we do not report these results in detail, but the majority of findings are consistent with the model correction procedures’ outcomes presented below. According to all MDE measures (d L 1 , d L 2 and d N∗ S D ), the CR-SETAR model, the CR-SETAR corrected by a VAR and the CR-SETAR model corrected with the k-NN approach show better performance indicating that the contribution of nonlinear models to a good forecast performance is significant.2 Moreover, these methods also present the highest forecast coverage and efficiency rates. As previously discussed, the closeness of the results of these two statistics can be taken as an indicator of the quality of the forecasts. Turning to the performance of individual attributes of forecasted ITS an interesting observation that can be made is that linear VAR and SETAR based models have only marginal improvements in the forecasting of center time series when compared to the rest of the models. However, the RMFSE and MAFE measures indicate significant improvements in radius time series forecasts when SETAR based models are employed. It means that the main predictive power of CR-SETAR is driven by the radius forecasts. This confirms previous results on forecasting stock returns with ITS techniques (see, e.g., Arroyo et al. 2011). Since the center of ITS can be associated with the conditional mean of the returns series and the radius can be used to estimate its volatility, the results obtained are also partially consistent with well-known results in the empirical finance literature. In many occasions the conditional mean of returns have close to white noise behavior and the main forecast power lies in the conditional variance. Finally, it is interesting to fully explore the advantages of the parametric estimates of CR-SETAR. As mentioned in the previous section the threshold variable z C,t−1 in (1) was set to be rt−1 based on the best model fit. Since the radius of ITS returns is associated with volatility, the estimated structure of the CR-SETAR model allows us to map the regimes of Eq. (1) into high- and low-volatility periods of the returns. Moreover, the estimate of the threshold variable γC provides a reference value for identification of both regimes. Figure 2 shows real and forecasted radius (blue and read colour, respectively), and the regime splitting for the estimation period (1997–2007) and for the forecast period (2008–2010). According to the estimate of γC most of the high volatility periods, as

2 Following the suggestion of one referee, to compare predictive accuracy we have also conducted Diebold–

Mariano tests (see Appendix) and the results obtained corroborate these conclusions.

123

Modeling and forecasting interval time series

Fig. 2 Regime splitting for forecasts and fitted values of radius

expected, have been found during recessions (the “dot-com” aftermath and the 2008 crisis).

5 Conclusion This paper develops new methods for the estimation and forecasting of interval-valued time series. We have introduced a threshold-type model for the center-radius representation of ITS and shown that the model is comparatively straightforward to estimate using sequential conditional LS (as in the conventional time series context). We also discussed combination and correction approaches to ITS forecasting. These methods may find constructive use in empirical applications as presented in the empirical section, in which they were applied on the S&P500 index returns. The empirical results obtained clearly suggest that nonlinear threshold-type models are able to improve forecast performance. Another important point to highlight from the results obtained is that the model is potentially useful in the analysis of volatility dynamics. In particular by distinguishing high volatile periods based on the behavior of the returns’ radius. Several extensions of the methods presented will be pursued in future work. First, it will be interesting to compare our results with those of other nonlinear models (such as, e.g., Markov switching, smooth transition, and ANN). Second, it will also be important to investigate how well the radius component of ITS does in forecasting volatility. Acknowledgments We are grateful to the participants of the ISF 2010 conference in San Diego and the SDA Workshop 2012 in Madrid for useful comments and suggestions. We would also like to thank Guest Editor Dr. Javier Arroyo, Coordinating Editor Professor Maurizio Vichi and three anonymous referees for helpful and constructive comments on an earlier version of this paper.

6 Appendix: Diebold–Mariano predictive accuracy tests Table 4 provides Diebold and Mariano (1995) test results for the comparison of different forecast methods. The results in this table are for center and radius forecasts and are based on loss differentials which are defined as the difference of squared forecast errors.

123

P. M. M. Rodrigues, N. Salish Table 4 DM test results for center / radius forecasts M1 M2 M1 × −1.811**

2.794***

M3

M4

M5

M6 −2.114**

M7

M8

M9

1.549*

0.039

1.059

2.262**

0.340

2.280**

3.464***

3.572***

4.057***

3.060***

2.971***

3.479***

3.188***

3.104***

1.984**

1.877**

2.076**

0.847

2.221**

1.986**

2.641***

2.515***

−1.816**

2.300**

1.839**

−2.504***

1.658** −0.299

M2 ×

×

M3 ×

×

×

M4 ×

×

×

×

M5 ×

×

×

×

×

M6 ×

×

×

×

×

−0.365

2.036** −1.556* −0.843

0.454

2.582*** −2.107** −2.827*** 0.927

−1.991**

−2.890*** −2.915***

3.003***

2.771**

2.168**

1.410*

1.732**

2.248**

2.053**

1.920**

0.459

2.006**

1.440*

1.220

−0.838

1.547*

1.367*

−2.780***

1.073

−0.674

2.005**

1.323*

−1.646*

2.193**

2.729***

2.459***

2.349***

×

2.946***

2.253**

4.457***

4.252***

2.277**

2.795***

2.513***

2.423***

1.103

0.495

−0.203

−0.599

×

×

×

×

×

×

M8 ×

×

×

×

×

×

×

−1.261 −1.299* ×

1.971** 1.650**

×

×

×

×

1.582*

1.823**

M7 ×

M9 ×

M10

×

×

×

×

1.362* 1.478* −1.419* −1.547*

***, ** and * denote 1 %, 5 % and 10 % significance level. The results above the main diagonal refer to Diebold–Mariano test results for the forecasts obtained from the center (top of each cell) and the radius (bottom of each cell) equation. M1 refers to the random walk with drift, M2 to the naïve model, M3 to the simple autoregressive framework, M4 to a linear VAR model for center and radius, M5 to the k-NN method based on NSD loss function, M6 the k-NN method based on L 2 , M7 to the proposed CR-SETAR model, M8 to the VAR model corrected with the k-NN approach, M9 to the CR-SETAR model corrected using a VAR, and M10 to the CR-SETAR model corrected using the k-NN approach

References Arroyo J, Maté C (2006) Introducing interval time series: accuracy measures. COMPSTAT 2006. Proceeding in Computational statistics, Heidelberg, pp 1139–1146 Arroyo J, Espínola R, Maté C (2011) Different approaches to forecast interval time series: a comparison in finance. Comput Econ 37:169–191 Arroyo J, Gonzalez-Rivera G, Maté C (2010) Forecasting with Interval and Histogram Data. Some Financial Applications. In: Ullah A, Giles D (eds.) Handbook of Empirical Economics and Finance, Chapman and Hall, pp 247–280 Bai J (1997) Estimating multiple breaks one at a time. Econ Theory 13:315–352 Beckers S (1983) Variance of security price return based on high, low and closing prices. J Bus 56:97–112 Cheung YW (2007) An empirical model of daily highs and lows. Int J Financ Econ 12:1–20 Chou RY (2005) Forecasting financial volatilities with extreme values: the conditional autoregressive range (CARR) model. J Money Credit Bak 37(3):561–582 Clements MP, Smith J (1999) A monte carlo study of the forecasting performance of empirical SETAR models. J Appl Econ 14(2):123–141 Diebold F, Mariano R (1995) Comparing predictive accuracy. J Bus Econ Statistics 13(3):253–263 Dueker M, Martin S, Spangnolo F (2007) Contemporaneous threshold autoregressive models: estimation, testing and forecasting. J Econ 141:517–547 Eitrheim Ã, Teräsvirta T (1996) Testing the adequacy of smooth transition autoregressive models. J Econ 74:59–75

123

Modeling and forecasting interval time series Gonzalo J, Pitarakis J-Y (2002) Estimation and model selection based inference in single and multiple threshold models. J Econ 110:319–352 Granger CWJ, Terasvirta T (1993) Modelling non-linear economic relationships. OUP Catalogue. Oxford University Press. ISBN 9780198773207 Guidolin M, Hyde S, McMillan D, Ono S (2009) Non-linear predictability in stock and bond returns: when and where is it exploitable? Int J Forecast 25:373–399 Hansen B (1996) Inference when a nuisance parameter is not identified under the null hypothesis. Econometrica 64:413–430 Hansen B (1997) Inference in TAR models. Studies in Nonlinear Dynamics and Econometrics 2 Henry Ó, Olekaln N, Summers PM (2001) Exchange rate instability: a threshold autoregressive approach. Econ Record 77:160–166 Hu C, He LT (2007) An application of interval methods to stock market forecasting. Reliab Comput 13(5):423–434 Hsu HL, Wu B (2008) Evaluating forecasting performance for interval data. Comput Math Appl 56:2155– 2163 Ichino M, Yaguchi H (1994) Generalized Minkowski metrics for mixed and feature-type data analysis. IEEE Trans Systems Man Cybern 24(1):698–708 Maia ALS, de Carvalho FdAT (2011) Holt’s exponential smoothing and neural network models for forecasting interval-valued time series. Int J Forecast 27:740–759 Muñoz SR, Maté AC, Arroyo J, Sarabia Á (2007) iMLP: Applying multi-layer perceptrons to interval-valued data. Neural Process Lett 25:157–169 Nieto FH (2005) Modelling bivariate threshold autoregressive processes in the presence of missing data. Commun Stat Theory Methods 34:905–930 Pitarakis J-Y (2006) Model selection uncertainty and detection of threshold effects. Studies in Nonlinear Dynamics and Econometrics 10 Timmermann A (2006) Forecast Combinations. Handbook of Economic Forecasting, Elsevier, Amsterdam Tong H (2011) Threshold models in time series analysis â 30 years on (with discussions by P. Whittle, M. Rosenblatt, B. E. Hansen, P. Brockwell, N. I. Samia and F. Battaglia). Stat Interface 4:107–136 Tsay RS (1998) Testing and modelling multivariate threshold models. J Am Stat Assoc 93:1188–1202 Zhang G (2003) Time series forecasting using a hybrid ARIMA and neural network model. Neurocomputing 50:159–175 Zou H, Yang Y (2004) Combining time series models for forecasting. Int J Forecast 20:69–84

123