Forecasting in the Presence of in and Out of Sample Breaks

Forecasting in the presence of in and out of sample breaks Abstract We present a frequentist-based approach to forecast time series in the presence of in-sample and out-of-sample breaks in the parameters of the forecasting model. We first model the parameters as following a random level shift process, with the occurrence of a shift governed by a Bernoulli process. In order to have a structure so that changes in the parameters be forecastable, we introduce two modifications. The first models the probability of shifts according to some covariates that can be forecasted. The second incorporates a built-in mean reversion mechanism to the time path of the parameters. Similar modifications can also be made to model changes in the variance of the error process. Our full model can be cast into a conditional linear and Gaussian state space framework. To estimate it, we use the mixture Kalman filter and a Monte Carlo expectation maximization algorithm. Simulation results show that our proposed forecasting model provides improved forecasts over standard forecasting models that are robust to model misspecifications. We provide two empirical applications and compare the forecasting performance of our approach with a variety of alternative methods. These show that substantial gains in forecasting accuracy are obtained. 1 Introduction Forecasting is obviously of paramount importance in time series analyses. The theory of constructing and evaluating forecasting models is well established in the case of stable relationships. However, there is growing evidence that forecasting models are subject to instabilities, leading to imprecise and unreliable forecasts. This is so in a variety of fields including macroeconomics and finance; e.g., Stock and Watson (1996) who documented widespread prevalence of instabilities in macroeconomic time series relationships. Prominent examples are forecasting inflation (Stock and Watson, 2007) and the equity premium (Pastor and Stambaugh, 2001). Also, Paye and Timmermann (2006) examined model instability in the coefficients of ex post predictable components of stock returns. See also Pesaran and Timmermann (2002), Rapach and Wohar (2006) and Pettenuzzo and Timmermann (2011). There is a vast literature on testing for and estimating structural changes within a given sample of data; see, e.g., Andrews (1993), Bai and Perron (1998, 2003) and Perron (2006) for a survey. Much of the literature does not model the breaks as being stochastic. Hence, the scope for improving forecasts is limited. There can be improvements by relying on the estimates of the last regime (or at least putting more weights on them) but even then such improvements are possible if there are no out-of-sample breaks. In the presence of out-ofsample breaks the limitation imposed by treating the breaks as deterministic mitigates the forecasting ability of models corrected for in-sample breaks. This renders forecasting in the presence of structural breaks quite a challenge; see, e.g., Clements and Hendry (2006). Some Bayesian models have been proposed to address this problem; see, e.g., Pesaran et al. (2006), Koop and Porter (2007), Maheu and Gordon (2008), Maheu and McCurdy (2009) and Hauwe et al. (2011). The advantage of the Bayesian approach steams from the fact that it treats the parameters as random and by imposing a prior (or meta-prior) distribution one can model the breaks and allow them to occur out-of-sample with some probability. Such methods can, however, be sensitive to the exact prior distributions used. We propose a frequentist-type approach with a forecasting model in which the changes in the parameters have a probabilistic structure so that the estimates can help forecast future out-of-sample breaks. Our approach is best suited to the case for which breaks occur both in and out-of-sample, which in particular avoids the problematic use of a trimming window assumed to have a stable structure. The method will work best indeed if there are many in-sample breaks, so that a long span of data is beneficial. This is unavoidable since good outof-sample forecasts of breaks require in-sample information about the process generating such 1 breaks, the more so the more efficient the forecasts will be. The same applies to previously proposed Bayesian methods, though the use of tight priors can partially substitute for the lack of precise in-sample information. Having said that, our method still yields considerable improvements even if relatively few breaks are present in-sample. Our approach is similar in spirit to unobserved components models in which the parameters are modeled as random walk processes. There are, however, important departures. Most importantly, a shift need not occur every period. It does so with some probability dictated by a Bernoulli process for the occurrence of shifts and a normal random variable for its magnitude. This leads to a specification in which the parameters evolve according to a random level shift process. Some or all of the parameters of the model can be allowed to change and the latent variables that dictate the changes can be common or different for each parameters. Also, the variance of the errors may change in a similar manner. The basic random level shift model has been used previously to model changes in the mean of a time series, whether stationary or long-memory, in particular to try to assess whether a seemingly long-memory model is actually a random level shift process or a genuine long-memory one; see Ray and Tsay (2002), Perron and Qu (2010), Lu and Perron (2010), Qu and Perron (2013), Xu and Perron (2014), Li et al. (2016) and Varneskov and Perron (2016). It has been shown to provide improved forecasts over commonly used short or longmemory models. Our basic framework is a generalization in which any or all parameters of a forecasting model are modeled as random level shift processes. To improve the forecasting performance we augment the basic model in two directions. First, we model the probability of shifts as a function of some covariates which can be forecasted. Second, we allow a meanreversion mechanism such that the parameters tend to revert back to the pre-forecast average. This last feature is especially influential in providing forecasting at long horizons. Functional forms for these two modifications are suggested for which the parameters can be estimated and incorporated in the forecast scheme to model the future path of the parameters. Modeling parameters as random level shifts has been suggested previously but, to our knowledge, only in a Bayesian framework. McCulloch and Tsay (1993) considered an autoregression in which the intercept is subject to random level shifts, though the autoregressive parameters are held fixed. They also allow the probability of shifts to depend on some covariates and changes in the variance of the errors though with a different specification. Gerlach, Carter and Kohn (2000) consider a class of conditionally linear Gaussian state-space models with a vector of latent variables indicating the occurrence of changes in the coefficients that follow a Markov process. Pesaran et al. (2006) extend the Markovian structure of Chib 2 (1998) with a fixed number of regimes by adopting a hierarchical prior with a constant transition probability matrix out of sample, thereby allowing breaks to occur at each date in the post-sample period. Koop and Potter (2007) consider models with a random number of regimes with the transitions from one regime to another being dictated by a Markov process and the durations of the regimes following a Poisson distribution. Giordani and Kohn (2008) extend their analysis, and that of Gerlach, Carter and Kohn (2000) to allow an arbitrary number of shifts occurring independently for the coefficients and error variance using a random level shift process with constant probability of shifts. Giordani, Kohn and van Dijk (2007) consider a class of conditionally linear and Gaussian state-space models which allows nonlinearity, structural change and outliers that can accommodate a fixed number of regimes with Markov transitions probabilities or random level shift processes, though in the applications they restrict the magnitudes of change and impose restrictive structures on the latent variables indicating the occurrence of changes. Groen, Paap and Ravazzolo (2013) use a model with random level shifts in the coefficients and error variance with constant probabilities to model and forecast inflation. Smith (2012) consider a Markov breaks regression model akin to a random coefficient model with all parameters changing at the same time and the probability of shifts being Markovian. As noted in some of the applications, the results can be quite sensitive to the prior used. Our approach is closest to that of McCulloch and Tsay (1993) except that we consider a general forecasting linear model with the same type of changes in coefficients and variance of the errors, allowing the probabilities of shifts to depend on some covariates. We also incorporate a mean-reversion mechanism. More importantly, we do not adopt a Bayesian approach and thereby bypass the need to specify priors and have the results influenced by them. Also, our focus is explicitly on providing improved forecasts. As stated above, the basic ingredient of the structure adopted has been considered previously, though not advanced as a widely applicable forecasting framework. Our aim is to generalized it and provide a ‘general purpose’ forecasting model that performs well for diverse scenarios with or without breaks. We also propose a novel algorithm to construct in-sample and out-of-sample confidence bands of the estimates and the forecasts. Our model can be cast into a non-linear non-Gaussian state space framework for which standard Kalman filter type algorithms cannot be used. The state space representation of our model is a linear dynamic mixture model in the sense that it is linear and Gaussian conditional on some latent random variables. Chen and Liu (2000) propose a special sequential Monte Carlo method, the mixture Kalman filter, which uses a random mixture of Gaussian distributions to approximate a target distribution. Giordani et al. (2007) discuss 3 the advantages of the class of conditionally linear and Gaussian state space models. The EM (Expectation Maximization) algorithm is used to obtain the maximum likelihood estimates of the parameters. This allows treating the latent state variables as missing data (Bilmes, 1998) and using a complete or data-augmented likelihood function which is easier to evaluate than the original likelihood. Since the missing information is random, the complete-data likelihood function is a random variable and we end up maximizing the expectation of the complete-data log-likelihood with respect to the missing data. Wei and Tanner (1990) introduce the Monte Carlo EM algorithm where the evaluation step is executed by Monte Carlo methods. Random samples from the conditional distribution of the missing data (state variables) can be obtained via a particle smoothing algorithm. The forecasting procedure is then relatively simple and straightforward once the model has been estimated. Simulations show that the estimation method provides reliable results in finite samples. The parameters are estimated precisely and the filtered estimates of the time path of the parameters follow closely the true process. To show the robustness of our forecasting model, we design simulations comparing the forecasting performances of various popular models (various form of the RLS models, historical average, rolling average, ARMA, ARIMA, Markov Switching, Time Varying Parameters (TVP)) when the Data Generating Process (DGP) is one of the forecasting models considered. The results show that our random level shift model with built-in mean reversion always performs nearly as well as the model corresponding to the true DGP, and can even be better (e.g., when the DGP is ARIMA). All other forecasting methods perform poorly in one or more of the cases considered. Hence, our method provides reliable results robust to a wide range of processes. We apply our forecasting model to two series which have been the object of considerable attention from a forecasting point of view. The emphasis is on the equity premium. We compare the forecast accuracy of our model relative to the most important forecasting methods applicable for this variable, as well as different forecasting sub-samples or periods. The results show clear gains in forecasting accuracy, sometimes by a very wide margin; e.g., over 90% reduction in mean squared forecast error relative to popular contenders. For this particular series,the TVP performs quite well being a close second best. To show the robustness of our forecasting model, we also consider the Treasury bill rate. Our method continues to provide the best forecasts overall, while the TVP model lead to very poor forecasts in most samples considered. Other applications can be found in Xu (2017). Finally, note that given the availability of the proper code for estimation and forecasting, the method is very flexible and easy to implement. For a given forecasting model, all that 4 is required by the users are: 1) which parameters (including the variance of the errors if desired) are subject to change; 2) whether the same or different latent Bernoulli processes dictates the timing of the changes in each parameters; 3) which covariates are potential explanatory variables to model the probability of shifts. The rest of the paper is organized as follows. Section 2 describes the basic model with random level shifts in the parameters. Section 3 discusses the modifications introduced to improve forecasting: the modeling of the probability of shifts and the allowance for a mean-reverting mechanism. Section 4 presents the estimation methodology: the mixture Kalman filtering algorithm in Section 4.1, the particle smoothing algorithm in Section 4.2, the Monte Carlo Expectation Maximization method to evaluate the likelihood function in Section 4.3. Section 5 introduces the construction of in-sample confidence bands and out-ofsample forecast bands. Section 6 provides forecasting simulations of various models to show the reliability and robustness of our proposed method. Section 7 contains the applications and comparisons with other forecasting methods. Section 8 offers brief concluding remarks. Detailed estimation algorithms are included in an appendix. 2 Model setup Let yt is a scalar variable to be forecasted with the following basic model: yt = Xt β t + et ; t = 1, ..., T, (1) where Xt a k-vector of covariates, T is the samples size and, in the base case, et ∼ i.i.d. N(0, σ 2e ). It is assumed that some or all of the parameters are time-varying and exhibit structural changes at some unknown time. The specification adopted for the time-variation in the β β , . . . , Kk,t ) and ηt = (η1,t , . . . , η k,t )0 ∼ parameters is β t = β t−1 + Ktβ η t , where Ktβ = diag(K1,t β ∼ Ber(p(j) ) and are independent across j. Hence, i.i.d. N(0, Σ). The latent variables Kj,t each parameter evolves according to a Random Level Shift (RLS) process such that the shifts β β . When Kj,t = 1, a shift are dictated by the outcomes of the Bernoulli random variables Kj,t β 2 ηj,t occurs drawn from a N(0, σ η,j ) distribution, otherwise when Kj,t = 0 the parameter does not change. The shifts can be rare (small values of p(j) ) or frequent (larger values of p(j) ). This specification is ideally suited to model changes in the parameters occurring at unknown dates. Many specifications are possible depending on the assumptions imposed on β β = . . . = Kk,t , we can interpret the model as one in which all paKtβ and Σ. First, when K1,t rameters are subject to change at the same times, akin to the pure structural change model of Bai and Perron (1998). A partial structural model, can be obtained by setting p(j) = 0 5 for the parameters not allowed to change, or equivalently by setting the corresponding rows β β = . . . = Kk,t is arguably the most interesting for and columns of Σ to 0. The case with K1,t a variety of applications. However, it is also possible not to impose equality for the different β . This allows the timing of the changes in the different parameters to be governed by difKj,t ferent independent latent processes. This may be desirable in some cases. For instance, it is reasonable to expect changes in the constant to be related to low frequency variations of the random level shifts type, while changes in the coefficients associated with random regressors to be related to business-cycle type variations. Of course, many different specifications are possible, and the exact structure needs to be tailored to the specific application under study. β are independent across j may The assumption that the latent Bernoulli processes Kj,t seem strong. It implies that the timing of the changes are independent across parameters. As stated above, this can be relaxed by imposing a perfect correlation, i.e., setting some latent variables to be the same. Ideally, one may wish to have a more flexible structure that would allow imperfect though non-zero correlation. This generalization is not feasible in our framework. In many cases, it may also be sensible to impose that Σ is a diagonal matrix. This implies that the magnitudes of the changes in the various parameters are independent. In our applications, we follow this approach as it appears the most relevant case in practice and also considerably reduces the complexity of the estimation algorithm to be discussed in Section 4. Hence, for the j th parameter β j (j = 1, . . . , k), we have β η j,t β j,t = β j,t−1 + Kj,t (2) β ∼ Ber(p(j) ). In some cases, it may also be of interest to where ηj,t ∼ N(0, σ 2η,j ) and Kj,t allow for changes in the variance of the errors, in which case et = σ ²,t ²t with lnσ 2²,t = φlnσ 2²,t−1 + Ktσ v²,t (3) where ²t ∼ N(0, 1), Ktσ ∼ Ber(pσ ) and v²,t ∼ N(0, σ 2v ). Note that when p(j) = pσ = 0 for all j, the model reduces to the classic regression model with time invariant parameters. When p(j) = 1 for all j and pσ = 0, it becomes the standard time varying parameter model; e.g., Rosenberg (1973), Chow (1984), Nicholls and Pagan (1985) and Harvey (2006). If φ = 1 in (3), we have a random level shift model for volatility. Adding that Ktσ = 1, we have the stochastic volatility modeled as a random walk and if |φ| < 1 and Ktσ = 1, the commonly used stochastic volatility (e.g., Hull and White, 1987). Stock and Watson (2007) used a similar unobserved component stochastic volatility (UC-SV) model to forecast inflation, in which the stochastic volatility equation is specified with φ = 1 and Ktσ = 1. 6 3 Modifications useful for forecasting improvements So far the framework is well tailored to model in-sample breaks in the parameters but does not allow future breaks to play a role in forecasting. In order to do so, we incorporate some modifications. Two features likely to improve the fit and the forecasting performance are to allow for changes in the probability of shifts and model explicitly a mean-reverting mechanism (j) for the level shift component. First, we specify the jump probability to be pt = f (γ, wt ), where γ is a m-vector of parameters, wt are m covariates that would allow to better predict the probability of shifts and f is a function that ensures pt ∈ [0, 1]. Note that wt needs to be in the information set at time t in order for the model to be useful for forecasting. We shall adopt a linear specification with the standard normal cumulative distribution function Φ(·), (j) (j) β ∼ Ber(pt ) with pt = Φ(r0 + r10 wt ), where r0 is a scalar and r1 and m-vector of so that Kj,t parameters. As similar specification can be made for the probability of the Bernoulli random variable Ktσ affecting the shifts in the variance of the errors. The second step allows a mean reverting mechanism. The motivation for doing so is that we often observe evidence that parameters do not jump arbitrarily and that large upward movements tend to be followed by a decrease. One example is when Xt ≡ 1, then β t stands for the unconditional mean of yt . As shown in Xu and Perron (2014), mean reversion properties are commonly detected in volatility series of equity returns. This feature can be beneficial to improve the forecasting performance if explicitly modeled. The specification adopted is: (t−1) ηj,t ∼ N(μη,j,t , σ 2η,j ), and μη,j,t = ρ(β j,t|t−1 − β̄ j ), where β j,t|t−1 is the filtered estimate (t−1) is the mean of all the filtered of the parameter subject to change at time t − 1 and β̄ j estimates of the jump component from the beginning of the sample up to time t − 1. This implies a mean-reverting mechanism provided ρ < 0. The magnitude of ρ then dictates the speed of reversion. If ρ = 0, there is no mean reversion. Note that only data up to time t − 1 is used in order to be useful for forecasting purposes. Also, it will have an impact on forecasts since being in a high (low) values state implies that in future periods the values will be lower (higher), and more so as the forecasting horizon increases. Hence, this specification has an effect on the forecasts of both the sign and size of future jumps in the parameters. Similar specifications can be made to pσ and v²,t for the changes in the variance of the errors. 4 Estimation methodology The model is within the class of conditional linear Gaussian State Space models of the form yt = Xt β t + et 7 (4) β t = β t−1 + Ktβ ηt , lnσ 2²,t = φlnσ 2²,t−1 + v²,t (5) where yt is the variable to be forecasted and (β t , Ktβ , lnσ 2²,t ) is the state vector. The measurement and transition equations are (4) and (5), respectively. Conditional on (Ktβ , lnσ 2²,t ), the system is a linear and Gaussian state space model and p(β t |Ktβ , lnσ 2²,t , Y ; Θ), where Y = (y1 , ..., yT ) and Θ is the vector of parameters, can be evaluated by the Kalman filter. The method of estimation is based on the mixture Kalman filtering using particle filters due to Chen and Liu (2000). Note that in (5), random level shifts in the stochastic volatility process can be added as in (3). See the appendix for details. 4.1 Mixture Kalman filtering Here xt denotes the state variable and λt ≡ (Ktβ , lnσ 2²,t ) the latent variables. Let y t = (y1 , ..., yt ), Λt = (λ1 , ..., λt ), and λt be realizations of Λt . The filtering distribution of xt is ¡ ¢ ¡ ¢ R p(xt |y t ) = p(xt |y t , λt )p(λt |y t )dλt , where p(xt |y t , λt ) ∼ N μt (λt ), Σt (λt ) , and μt (λt ), Σt (λt ) can be obtained via the Kalman filter for a given trajectory λt . The mixture Kalman fil(1) (M) ter uses a weighted sample of the indicators St = {(λt,(1) , wt ), ..., (λt,(M) , wt )} to rep(i) resent the distribution p(λt |y t ), where wt are weights defined below and λt,(i) are simβ,(i) β,(i) ulated latent variables; e.g., in the basic model λt,(i) ≡ (K1 , ..., Kt ), so that given a jump probability p, λt,(i) can be generated as random draws from the Bernoulli distribution with probability p. One then uses a random mixture of Gaussian distributions P P (i) (i) t,(i) ), Σt (λt,(i) )), where Wt = M Wt−1 M i=1 wt N(μt (λ i=1 wt , to approximate the target distribution p(xt |y t ). The detailed mixture Kalman filtering algorithm is provided for the basic model (1)-(2), and the one with stochastic volatility (1)-(3) in the appendix. To illustrate the adequacy of this method, we present simple illustrative examples. First, β t is generated using (1) and (2) with mean reversion and time varying probability with (r0 , r1 , σ e , σ η , ρ) = (−1.96, 4, 0.2, 0.2, −0.1) and T = 1000. The covariate wt is set to be 1 every 50 observations, 0 otherwise. Figure 1 presents the true path of β t along with the filtered estimates of β t showing a close agreement between the two. Figure 2 considers the case with stochastic volatility, as specified by (1)-(3) with mean reversion and time varying probability and (r0 , r1 , φ, σ v , σ η , ρ) = (−1.96, 4, 0.95, 0.2, 0.2, −0.1). A plot the true path of β t along with the filtered estimates of β t are presented in Panel A. The corresponding values for the volatility process are presented in Panel B. Again, the filtered values closely follow the true paths in both cases. While obviously limited, the cases reported are representative of what one can expect in most cases (from unreported additional simulations performed), showing the adequacy of the filtering method adopted. 8 4.2 Particle smoothing (i) The particle smoothing algorithm is designed to obtain particle smothers {st }M i=1 with (i) M certain weights {wt }i=1 from p(xt |y T ). Godsill et al. (2004) provide a forward-filtering and backward-simulation smoothing procedure. It allows drawing random samples from the joint density p(x0 , x1 , . . . xT |y T ), not only the individual marginal smoothing densities p(xt |y T ). The smoothing algorithm relies on a pre-filtering procedure and a previously obtained set of (i) (i) filters {wt , xt }M i=1 for each time period. The main ingredients are the relations: p(x1 , . . . , xT |y T ) = p(xT |y T ) QT −1 t=1 p(xt |xt+1 , . . . , xT , y T ) p(xt |xt+1 , . . . , xT , y T ) = p(xt |xt+1 , y t ) = p(xt |y t )p(xt+1 |xt )/p(xt+1 |y t ) ∝ p(xt |y t )p(xt+1 |xt ) The first equality follows from the Markov property of the model and the second from (i) t Bayes’ rule. Since random samples {xt }M i=1 from p(xt |y ) can be obtained from the mixture P (i) (xt ) Kalman filtering algorithm, p(xt |xt+1 , . . . , xT , y T ) can be approximated as M i=1 wt|t+1 δ x(i) t PM (i) (i) (i) (i) (i) with modified weights wt|t+1 = wt p(xt+1 |xt )/ i=1 wt p(xt+1 |xt ), where δ x(i) (xt ) is the t Dirac delta function. This procedure is performed in a reverse-time direction conditioning on future states. Given a random sample {st+1 , . . . , sT } drawn from p(xt+1 , . . . , xT |y T ), we take one step back and sample st from p(xt |st+1 , . . . , sT , y T ). The smoothing algorithm is summarized in the appendix in the context of the various versions of our model. 4.3 MCEM algorithm Frequentist likelihood-based parameter estimation of conditional linear and Gaussian state space models using the mixture Kalman filters and smoothers is not straightforward. The gradient-based optimizer suffers from a discontinuity problem caused by the resampling. Here, we follow the Monte Carlo Expectation Maximization (MCEM) method proposed by Olsson et al. (2008). The Basic EM algorithm is a general method to obtain the maximumlikelihood estimates of the parameters of an underlying distribution from a given data set with missing values. Suppose the complete data set is Z = (Y, X), in which Y is observed but X is unobserved, and Θ is the parameter vector. For the joint density p(z|Θ) = p(y, x|Θ) = p(y|Θ)p(x|y, Θ), we define the complete-data likelihood function by L(Θ|Y, X) = p(Y, X|Θ). The original likelihood L(Θ|Y ) is the incomplete-data likelihood. Since X is unobserved and may be generated from an underlying distribution, e.g., the transition equation in a state space model, L(Θ|Y, X) is indeed a random variable. Therefore, we maximize the expectation of logL(Θ|Y, X) with respect to X, with the expectation, conditional on Y and 9 some input value for the parameters Θ(k−1) , defined by: R Q(Θ, Θ(k−1) ) = E[logL(Θ|Y, X)|Y, Θ(k−1) ] = logp(Y, x|Θ)p(x|Y, Θ(k−1) )dx which permits an iterative procedure to update the values of the parameters Θ. The difference between the MCEM and the basic EM algorithms is that when evaluating Q(Θ, Θ(k−1) ), the MCEM uses a Monte-Carlo based sample average to approximate the expectation. The P (i) Monte Carlo Expectation or E-step is Q∗ (Θ, Θ(k−1) ) = M −1 M i=1 log(p(Y, x |Θ)), where (k−1) ). Given current parameter estimates, ran{x(i) }M i=1 are random samples from p(x|Y, Θ (i) dom samples from p(x|Y, Θ(k−1) ) are simply the particle smoothers {st }M i=1 obtained as (k) (k−1) ). These two described above. The Maximization or M-step is Θ = arg maxΘ Q(Θ, Θ steps are repeated until Θ(k) converges. The rate of convergence has been studied by many researchers; e.g., Dempster et al. (1977), Wu (1983) and Xu and Jordan (1996). In the context of the simple version of our model, the specifics of the algorithm are in the appendix. Overall, the estimation procedure is the following. Let Θ(0) be a vector of initial (i) parameter values: 1) (Mixture Kalman filtering): obtain mixture Kalman filters {xt }M i=1 from p(xt |y t , Θ(k−1) ), i = 1, ..., M, t = 1, ..., T ; 2) (Particle smoothing): obtain particle (i) T (k−1) ), i = 1, ..., M, t = 1, ..., T ; 3) (Estimation): evaluate smoothers {st }M i=1 from p(xt |y , Θ (i) M (k−1) ) using {st }i=1 and maximize it to obtain updated parameter estimates Θ(k) ; 4) Q(Θ, Θ Repeat steps 1-3 with k updated to k + 1 until the parameter estimates converge. 5 In-sample confidence bands and out-of-sample forecast bands We now propose a simulation based method to construct in-sample confidence bands and outof-sample forecast bands following a modification of the method of Blasques et al. (2016) who dealt with observation-driven time varying parameter models, for which the observations are given by yt ∼ p(yt |β t ; θ) with β t following the updating equation β t+1 = φ(β t , yt ; θ), where φ(.) is a differentiable recurrence function and θ is the static parameter. Their framework does not aply to ours since it is a parameter-driven time varying model. Perron and Xu (2016) pointed out that the updating equation (process) for the time-varying parameters in parameter-driven models can be written as β t+1 = h(β t , η t ; ψ), where ηt ∼ gη (ψ) is the idiosyncratic innovation and ψ is the static parameter. The parameter β t then follows a recurrence process with its own innovations. Therefore, the in-sample confidence bands need to incorporate both parameter uncertainty and innovation uncertainty. The parameter estimate ψ̂ is constructed via Monte Carlo maximum likelihood estimation. Let the estimate of the asymptotic covariance matrix of ψ̂ be defined by Σ̂ = −{∂ 2 log L̂(ψ̂)/∂ψ∂ψ 0 }−1 , where 10 L̂(ψ̂) is the Monte Carlo estimate of the likelihood function evaluated at ψ̂. The estimate Σ̂ can be computed numerically. Once an estimate of the asymptotic distribution of ψ̂ is obtained, the in-sample confidence bands for β̂ t+1 can be constructed using simulation methods similar to the filtering forecast band method proposed in Blasques et al. (2016). (i) The procedure can be described as follows: 1) Draw M parameter values ψ̂ from the (i) (i) asymptotic distribution ψ̂ ∼ N(ψ̂, T −1 Σ̂); see Olson and Ryden (2008); 2) Given ψ̂ , and (i) (1) (S) (s) for each time t, draw S sequences η t , ..., ηt from the estimated density ηt ∼ gη (ψ̂ ) for (s) (1) (S) s = 1, ..., S and t = 1, ..., T ; 3) Given the observations η t , ..., ηt , the filtered sequence β̂ 1 , (s) (s) (s) (s) (i) ..., β̂ T is obatined via the updating function β̂ t+1 = h(β̂ t , ηt ; ψ̂ ); 4) Repeat steps 2-3 for (i),(s) ; 5) Calculate the appropriate percentiles i = 1, ..., M to obtain M × S filtered paths of β̂ t (i),(s) to obtain the in-sample confidence bands for β̂ t . for each t over the M × S draws of β̂ t The procedure to construct the out-of-sample forecast bands for β̂ t+h is actually the (i),(s) same as described above. We simply need to obtain M × S extrapolated paths of β̂ t+h to compute the percentiles. To illustrate, we again use a simple example. The true process for β t is generated using (1)-(2) with mean reversion and time varying probability with parameters (r0 , r1 , σ e , σ η , ρ) = (−1.96, 4, 0.2, 1, −0.5). The covariate wt is set to be 1 every 20 observations, 0 otherwise. The computation of the in-sample and out-of-sample bands are based on M = 1000 and S = 1000 simulations. The number of observations is 500. Panel A of Figure 3 presents the true process β t , the filtered estimates and the 2.5% ad 97.5% quantiles of the simulated distribution. The confidence bands are quite narrow around the true process showing precisely estimated parameters. Panel B presents the results for the out-of-sample confidence bands. The forecasting horizon is set to be 100 steps. The figures show the forecasts and the 2.5% ad 97.5% quantiles of the simulated distribution. In order to show the relevance of the mean reversion scheme and modelling the changes in the probability of shifts, we consider two different illustrations using the realizations from Panel A, one starting from a high value (top figure) and the other starting from a low value (bottom figure). The results show that in both cases, the forecasts revert to mean level as the forecasting horizon increases, which leads to improved forecasts over simply using the value of the last regime, say. When constructing the forecasts, we model the probability of shifts using the same sequence of wt . This is done for simplicity of illustration (as it is a deterministic sequence) and allows us to gauge important features. It does not imply the same shifts, though they occur at the same time as the shifts in the out-of-sample data. The following features are noteworthy: 1) The size of the shifts is always in a direction that induces mean reversion when −1 < ρ < 0 (this can be shown theoretically to be a 11 feature of our forecasting model); 2) the magnitudes of the shifts decrease as the forecast horizon increases, which mirrors the fact that the mean-reverting effect is itself decreasing. Overall, the confidence bands are reasonable given the wide fluctuations in the out-of-sample values. They are narrow for short-horizons and increase as the forecasting horizon increases, as expected. In both cases, the true values of the series are within the 95% confidence band. 6 Simulations We now show, via simulations, the reliability and robustness of RLS type models in forecasting even with misspecified models. We consider nine Data Generating Processes (DGPs): 1) RLS basic model (infrequent large breaks): yt = β t + et , with β t = β t−1 + Kt η t , where Kt ∼ Ber(p), et ∼ N(0, σ 2e ), η t ∼ N (0, σ 2η ), θ = (p, σ e , σ η ) = (0.025, 2, 2). 2) RLS basic model (frequent small breaks): Same as in DGP-1 with (p, σ e , σ η ) = (0.5, 2, 0.2). 3) RLS with mean reversion: Same as in DGP-1, except that the probability of shifts is now a function of some covariate wt and β t follows a mean reverting process; i.e., pt = Φ(r0 + r1 wt ), (t−1) ); θ = (r0 , r1 , σ e , σ η , ρ) = (−1.96, 4, 0.2, 0.2, −0.1). ηt ∼ N(μηt , σ 2η ), μηt = ρ(β t|t−1 − β̄ The covariate wt is set to be 1 every 50 observations, 0 otherwise so that the probability of level shifts is small most of the time and close to 1 every 50 periods. 4) RLS without mean reversion: Same as DGP-3 with ρ = 0. 5) RLS_SV: Same as DGP-3, with added stochastic volatility to the error term of the form et = σ ε,t εt , εt ∼ N(0, 1), with lnσ 2ε,t = φlnσ 2ε,t−1 +vε,t , where vε,t ∼ N(0, σ 2v ) and independent of εt and et ; θ = (r0 , r1 , φ, σ v , σ η , ρ) = (−1.96, 4, 0.95, 0.2, 0.2, −0.1). 6) ARMA(1,1) (Autoregressive and Moving Average process): (1 − φL)yt = (1 + θL)εt , φ = 0.95, θ = −0.5 and εt ∼ N(0, 1). 7) ARIMA(1,1,1) (Autoregressive Integrated and Moving Average process): (1 − L)(1 − φL)yt = (1 + θL)εt , φ = 0.1, θ = −0.5 and εt ∼ N(0, 1). 8) TVP (Time Varying Parameter Model): yt = β t + et with β t = β t−1 + η t , where et ∼ N(0, σ 2e ) and η t ∼ N(0, σ 2η ) independent of each other; (σ e , σ η ) = (0.2, 0.2). 9) Markov Switching (MS): A two states regime switching model (e.g., Hamilton, 1994): yt = μSt + et , where et ∼ N(0, σ 2St ), St = 1, 2, with [μ1 , μ2 ] = [0.5, −0.5], [σ 21 , σ 22 ] = [1, 2] and the elements of the transition matrix from state i to state j for i, j = 1, 2 given by: p11 = 0.95, p12 = 0.1, p21 = 0.05, p22 = 0.9. In each case, we generate 100 true data paths and 1000 observations for each path. We use the first 800 observations for in-sample estimation and the rest to evaluate out-of-sample 12 forecasting accuracy. The forecasting horizon is up to 24 periods. The forecasting models considered are: the ‘RLS_m’: the RLS model with mean reversion and time varying probability; ‘RLS_SV’: the RLS model with mean reversion, time varying probability and stochastic volatility; ‘Average’: the historical average, namely the average over all observations in the expanding in-sample period; ‘Rolling’: the average of the last 50 observations of the in-sample period; ‘ARMA’: an ARMA(1,1) model; ‘ARIMA’: an ARIMA(1,1,1) model; ‘TVP’: a Time Varying Parameter Model as specified in DGP 8; ‘MS’: a Markov switching model as described in DGP 9. For DGP (1)-(2), we also consider the basic RLS model without mean reversion nor time varying probability, which acts as the benchmark model. For each DGP, we report the relative MSFEs of some other misspecified models with respect to the benchmark model, which is, in all cases, the true model with estimated parameters. The results are summarized in Table 1. Numbers smaller than 1 indicate a better forecasting performance than that obtained with the corresponding true model. Bold numbers indicate the smallest relative cumulative MSFEs for a given DGP and forecast horizon. Consider first the results in Panels 1-5, for which some type of RLS model is the true DGP. With few exceptions, the best performing forecasting model is the ‘RLS_m’. In the few cases for which it is not the best, the preferred one is the ‘RLS_SV’ for forecast horizons equal to 4 and 12 for DGP-2. The differences are, however, minor between the two. From DGP 1 and 2, the performance is shown to be robust to the frequency and the size of the breaks. The ’RLS_m’ even have relative RMSEs smaller than 1 for the large infrequent breaks model, which can be viewed as basically a structural change model. What is especially interesting is that introducing a mean reverting component even when not present leads to better forecasts, see Panel 4. The ‘TVP’ and ‘Markov Switching’ models perform poorly, especially at long-horizons. The ‘ARMA’ and ‘ARIMA’ models perform quite well but still produce inferior forecasts compared to the ‘RLS_m’. The ‘historical average’ is prone to severe deficiencies; e.g. Panel 1. The ‘rolling average’ has better performance compared to ‘historical average’ in most cases but is still inferior to the ’RLS_m’. We now turn to Panels 6-9, when the true DGP is not RLS. The results show that the RLS type models still have robust or even better performance compared to the benchmark model. The ‘RLS_m’ or ‘RLS_SV’ are second best (relative to the benchmark model) in most cases. In cases of model misspecifications, the performances of the various alternative models considered can be very poor; e.g. DGPs 7 and 8 for the ‘historical average’ and the ‘rolling window average’, DGP 8 for ‘Markov Switching’ and DGP 9 for ‘TVP’. As for the ‘ARMA’ and ‘ARIMA’ models, the performances are robust but still inferior to the 13 ’RLS_SV’ even under an ARIMA model specification for longer horizons. The results show that our random level shift model with built-in mean reversion always performs nearly as well as the model corresponding to the true DGP, and can even be better (e.g., when the true DGP is RLS without mean reversion or ARIMA). All other forecasting methods perform very poorly in one or more of the cases considered. Hence, our method provides reliable results that are robust to a wide range of processes. 7 Forecasting applications We consider two applications pertaining to variables which have been the object of sustained attention in the literature: the equity premium and the Treasury Bill rates, with the emphasis on the former. We compare the forecast accuracy of our model relative to the most important forecasting methods applicable for each variable. For the equity premium, the Time Varying Parameter Model (TVP) performs well being a close second best. As shown in the simulations, it is not robust to a variety of DGPs, while our method is. To illustrate this, we also consider the Treasury Bill rate for which our method still provides the best forecasts overall, while the TVP model leads to very poor forecasts in most samples. The out-of-sample forecasts are constructed in two steps. The first involves forecasting the covariates wt using a preliminary model; e.g., using an AR(k) or the random level shift model with a fixed probability of shift. The h-step ahead forecast of the jump probability is then pt+h|t = Φ(r̂0 + r̂1 wt+h|t ) where wt+h|t is the h-step ahead forecast of wt+h at time t and (r̂0 , r̂1 ) are the parameter estimates. We also use forecasted values for Xt+h and wt+h using an AR(p) model with p selected using the Akaike Information Criterion (AIC) with a maximal value of 4. The second step is to forecast {β t+s }hs=1 . The 1-step-ahead forecast P (i) (i) (i) is calculated as β t+1|t = E[β t+1 |It ] = M i=1 wt f (β t+1|t ), where β t+1|t is obtained via the Kalman filtering steps. For h step-ahead forecasts, β t+h|t = E[β t+h |It ] can be calculated recursively by repeating the filtering algorithm from time t + 1 to t + h, and treating the observations {yt+s }hs=1 as missing values. We can continue to apply the above algorithm setting vt = 0, Kt = 0 for t = t + 1, ..., t + h. Throughout, the out-of-sample forecasting experiments aim at evaluating the experience of a real-time forecaster by performing all model specifications and estimations using data through date t, making a h-step ahead forecast for date t + h, then moving forward to date t + 1 and repeating this through the sub-sample used to construct the forecasts. Unless otherwise indicated, the estimation of each model is recursive, using an increasing data window starting with the same initial observations. In all cases, we allow mean reversion in the parameters when constructing forecasts using our 14 RLS model. The forecasting performance is evaluated using the mean square forecast error −1 PTout 2 (MSFE) criterion defined as MSF E(h) = Tout t=1 (ȳt,h − ȳt+h|t ) , where Tout is the number P P of forecasts produced, h is the forecasting horizon, ȳt,h = hk=1 yt+k and ȳt+h|t = hk=1 yt+k|t with yt+k the actual observation at time t + k and yt+k|t its forecast conditional at time t. To ease presentation, the MSFE are reported relative to some benchmark model, usually the most popular forecasting model in the literature. The cumulative MSFE defined above gives the same relative measure of forecast performance as the root mean squared errors. Since our interest is not in the absolute level, it makes no difference which one is used. 7.1 Equity premium Forecasts of excess returns at both short and long-horizons are important for many economic decisions. Much of the existing literature has focused on the conditional return dynamics and studied the implications of structural breaks in regression coefficients including the lagged dividend yield, short-term interest rate, term spread and the default premium. However, most of the research has focused on modeling the equity premium assuming a certain number of structural breaks in-sample while ignoring potential out-of-sample structural breaks. Recently, Maheu and McCurdy (2009) studied the effect of structural breaks on forecasts of the unconditional distribution of returns, focusing on the long-run unconditional distribution in order to avoid model misspecification problems. Their empirical evidence strongly argue against ignoring structural breaks for out-of-sample forecasting. We consider using our forecasting model with different specifications. One models the unconditional mean of excess returns incorporating random level shifts in mean, with the time varying jump probabilities influenced by the lagged value of the absolute rate of growth in the earning price (EP) ratio. We also consider a conditional mean model using the dividend yield as the covariate. Following Jagannathan et al. (2000), we approximate the equity premium of S&P 500 returns as the difference between stock yield and bond yield. The data were obtained from Robert Shiller’s website (http://www.econ.yale.edu/~shiller/data.htm). According to Gordon’s valuation model, stock returns are the sum of the dividend yields and the expected future growth rate in stock dividends. We use the average dividend growth rate (over the pre-forecasting sample) to proxy for the expected future growth rate. The data consist of monthly series and cover the period from 1871:1 to 2012.5. High quality monthly data are available after 1927, before 1927 the monthly data are interpolated from lower frequency data. We use the 10-years Treasury constant maturity rate (GS10) as the risk free rate. 15 We start with a simple random level shift model without explanatory variables given by: yt = β t + et , β t = β t−1 + Ktβ η t (6) (t−1) ), Ktβ ∼ Ber(pt ) with where et ∼ i.i.d.N(0, σ 2e ), ηt ∼ i.i.d.N(μηt , σ 2η ), μηt = ρ(β t|t−1 − β̄ pt = Φ(r0 + r1 wt ). The covariate wt used to model the time variation in the probability of shifts is the lagged absolute value of the rate of change in the EP ratio, given that large fluctuations in the EP ratio are likely to induce a higher probability that excess stock returns will experience level shifts. We also consider a conditional forecasting model that uses the lagged dividend price (DP) ratio as the regressor. The specifications are yt = β 1t + β 2t dpt−1 + et (7) where, with β t = (β 1t , β 2t ), β t = β t−1 + Ktβ ηt , and dpt is the DP ratio. Lettau and van Nieuwerburgh (2008) analyzed the implications of structural breaks in the mean of the DP ratio for conditional return predictability. Xia (2001) studied model instability using a continuous time model relating excess stock returns to dividend yields. They specify β t to follow an Ornstein—Uhlenbeck process and the ensuing estimates of the time varying coefficient β 2t revealed instability of the forecasting relationship. Hence, instabilities have been shown to be of concern when using this conditional forecasting model, which motivates the use of our forecasting model. Besides the addition of the lagged DP ratio as regressors, the specifications are the same as for the unconditional mean model (6). We consider various versions depending on which coefficients are allowed to change and if so whether they change at the same time. These are: 1) the unconditional mean model (6) with level shifts, 2) the conditional mean model (7) with the constant allowed to change β β 6= 0, K2t = 0), 3) the conditional mean model (7) with the coefficient on the lagged DP (K1t β β = 0, K2t 6= 0). We compare our forecasting model with the most ratio allowed to change (K1t popular forecasting models used in the literature. These are: 1) the historical average (used as the benchmark model); 2) a rolling ten-years average; 3) the conditional model with the lagged DP ratio as the regressor without changes in the parameter; 4) a rolling version over ten years of the model previously stated in 3); 5) a TVP model with the unconditional mean following a random walk; 6) a two-states regime switching model. We first consider 1998-2012 as the forecasting period, with forecasting horizons 1, 6, 12, 18, 24, 30 and 36 months. The results are presented in Table 2.1. All three versions involving random level shifts perform very well and are comparable. The best model for short horizons less than 6 months is the conditional mean model (7) with the constant allowed to 16 β β change (K1t 6= 0, K2t = 0), though the difference are quite minor. For longer horizons, the conditional mean model (7) with the coefficient on the lagged DP ratio allowed to change β β = 0, K2t 6= 0) is the best. Of note is that our model performs much better than any (K1t competing forecasting models except the TVP model. This is especially the case at shorthorizons, for which the gain in forecasting accuracy translates into a reduction in MSFE of up to 90% when compared to the conditional model with no breaks (and even more so when compared to the rolling 10 years average or the historical average, the latter performing especially badly). At longer horizons, the conditional mean model (7) with level shifts still perform better than the conditional model with constant coefficients but to a lesser extent. The rolling version of the DP ratio model performs better than the one using the full sample for short horizons but less so at long horizons. In no case is it better than any of the versions with random level shifts. We also provide p-values from the Model Confidence Set (MCS) of Hansen et al (2011) with p-values greater than 0.1 indicating that the model belongs to the 10% MCS. The conditional mean model with the coefficient on the lagged DP ratio allowed β β = 0, K2t 6= 0) belongs to the MCS for all forecasting horizons. Other RLS to change (K1t type models and the TVP model belong to the MCS for 1-step-ahead forecasts. This can be viewed as strong evidence that the performance of our RLS model is good at forecasting the equity premium compared to popular candidates. To assess the robustness of the results we also consider the forecasting period 1988-1996, given that it offers an historical episode with different features; see Table 2.2. Note that the conditional mean model with constant parameters now performs very poorly with MSFEs more than four times those of the rolling 10 years average. The benchmark historical average performs even worse during this time period. On the other hand, the models with random level shifts continue to perform very well, with MSFEs around 0.2% of the historical average at short horizons, and around 2.5% at longer horizons up to 60 months. All models with random level shifts have comparable performance at short horizons, but the unconditional mean model (6) with level shifts is best at longer horizons. Meanwhile, the TVP model is also a strong candidate being the best at very short horizons and remaining in the 10% confidence set for all horizons. The conditional mean model with the coefficient on the lagged DP-ratio allowed to change also belongs to the 10% MCS for h > 1. Given that the results show impressive improvements in forecast accuracy using our forecasting method and the importance of forecasting the equity premium, we performed further sensitivity analyses using a different data set. The data are now those used in Welch and Goyal (2008); downloaded from Amit Goyal’s website (http://www.hec.unil.ch/agoyal/). It 17 is also the same dataset used by Pettenuzzo et al. (2014). As we will show, the improvements in forecast accuracy using our framework continue to be large whatever the sampling intervals used (yearly or monthly). We follow the common practice of simply calculating the equity premium as the historical average difference between returns on stocks and on risk-free assets. Annualized equity premiums are calculated from monthly data as compounded excess returns. Goyal’s data set goes back to 1871 and includes monthly and annual data. The most recent update ends in 2014. An added advantage of using Goyal’s dataset is that many other economic variables available for a long span of time. Here, we use the book-to-market ratio as the regressor to help predict the equity premium, which is one of the three factors in Fama and French (1993) three factors model. The Book-to-Market Ratio (BM) is the ratio of book value to market value for the Dow Jones Industrial Average. The covariate used to forecast the level shift probability is the absolute change in the EP ratio. In all cases, we use the historical average as the benchmark model. To analyze the effect of the oil shock, the forecast period is from 1975 to 2009 for both annual and monthly data. The results are presented in Tables 3.1 (annual) and 3.2 (monthly). They show that all three models we propose beat the benchmark for all sampling intervals. We obtain 80% reduction in cumulative MSFE for one-step-ahead forecasts and 58% for longer horizon forecasts. The improvements becomes larger when higher frequency data are used. For monthly data, we get 98% reduction in MSFE for short horizon forecasts and almost 74% reduction for long horizon forecasts. Note that from Table 3.1 (annual data) almost all competing models belong to the 10% model confidence set except for the rolling 10 years average and the constant parameter model with the BM ratio as a regressor for longer horizons forecasts. On the other hand, from Table 3.2 (monthly data) the TVP model for short horizons forecasts and the constant parameter model with a rolling estimation window for longer horizons forecasts are the only competing models that remain in the 10% MCS. With annual data, the total number of observations for out-of-sample forecasting is 34, considerably less than for the monthly case. The small samples available for MCS testing is likely the reason for the difficulty in selecting the “competing models”. Nevertheless, the relative MSFEs indicates relative performances similar to those in Tables 2.1-2.2. We also looked for other regressors such as financial variables with the purpose of further improving forecasting accuracy. Due to the short span of the financial data, we use monthly data from 1990/01/31 to 2008/12/31 for in-sample estimation and forecast from 2009 to 2015 at horizons up to 24 months. The financial variables are the VIX index from the Chicago Board Options Exchange (CBOE) and the returns on the monthly S&P 500 index 18 option. The results are reported in Table 4, which indicate a gradually improving forecasting performance for horizons longer than 12 months with these two additional regressors. For a 24 months forecasting horizon, the cumulative MSFE of the conditional mean model including financial variables is only half that with only the DP ratio. In summary, the results provide strong evidence that our forecasting model offers marked improvements in forecast accuracy. It does so at all horizons with results robust to different forecasting periods and data sets. It remains that the TVP model is a close second best forecasting model for the equity premium series analyzed. According to the simulations, the TVP model is much less robust to model mispecification than the RLS-type models. To illustrate this issue, we next consider the issue of forecasting the interest rate. 7.2 Interest rate Another variable of interest, which has attracted attention from a forecasting perspective, is the U.S. T-bill rate. Various studies have shown that it exhibits structural instability in both mean and variance (e.g., Garcia and Perron, 1996, Gray, 1996, Ang and Bekaert, 2002, and Pesaran et al., 2006). We use monthly data on the 3-months Treasury Bill rate from 1947:07-2002:12, obtained from the Federal Bank of St. Louis database; as in Pesaran et al. (2006). The period prior to 1968:12 is used for in-sample estimation, and we consider forecasting horizons of 12, 24, 36, 48 and 60 months. The basic model adopted is a simple AR(1) process given by: yt = β 1t + β 2t yt−1 + ²t . In all cases, we allow mean-reversion in the parameters and the covariate wt used to model the time-varying probabilities of shifts is the lagged value of the growth rate of GDP when a single latent Bernoulli variable is present. When two are present, the additional covariate is the lagged value of the absolute change in β = 0); stock returns (S&P 500). We consider two possible specifications: 1) ‘AR_K2t ’ (K1t β β 2) ‘AR_K1t , K2t ’ with K1t and K2t allowed to be different latent Bernoulli processes. The performance of the models is assessed relative to four competing forecasting methods1 : 1) ‘Recursive OLS’: a recursively estimated first-order autoregression with fixed parameters, used as the benchmark model; 2) ‘Rolling 5 years’: a first-order autoregression with fixed parameters estimated using a 5-year rolling window; 3) ‘Rolling 10 years’: a first-order autoregression with fixed parameters estimated using a 10-year rolling window; 4) ‘TVP’: a 1 We also compared the forecasting performance using Pesaran et al. (2006)’s composite and last regime model. Our RLS model performs better in most cases for forecasts computed every 12 months as in their paper. We do not include their model in our comparisons set due to computational constraints. Their method is highly computationally intensive and we could not apply it to forecasts computed every months. The sample obtained using forecasts computed every 12 months makes it too small for a valid MCS testing. 19 time-varying probability model in which β t = (β 1t , β 2t ) is modelled as a random walk. The results are presented in Table 5 for various forecast periods and forecast horizons h = 12, 24, 36, 48, 60 months. Consider first the results for the longest forecasting period 1968-2002. Here, the best forecasting model for all horizons is the ‘AR_K1t , K2t ’ with both the constant term and the AR coefficient allowed to follow a random level shift process. The gains in forecast accuracy vary between 3% and 8% and increase as the forecasting horizon increases. We then separate the forecasting period into three decades: the 70s, the 80s and the 90s. In the 70s, the 5 years rolling-average is overall the best predictor, though the ‘AR_K1t , K2t ’ model catches up and and is superior at h = 60. For the 80s and 90s, the best forecasting model is again the ‘AR_K1t , K2t ’ with both the constant term and the AR coefficient in the AR regression allowed to follow a random level shift process. The improvements in forecast accuracy are on average 5% reduction in MSFE. The ‘AR_K1t , K2t ’ performs the best in 16 out of 20 cases and is in the 10% MCS for 17 out of 20 cases. Note, however, that all models with random level shifts in parameters perform much better than the TVP model for the different time periods and horizons. The gains in forecast accuracy when using the RLS-type models over the TVP model are substantial. For the full sample, they range from a (roughly) 50% reduction at the shortest horizon to a 90% reduction at the longest one. For the 90s sub-sample, the corresponding reductions are of the order of 65% to 96%. The only period in which the TVP model does not perform badly is the 70’s, but it is still inferior to the RLS models and not within the 10% MCS. These results are consistent with the simulations reported in Section 6, which show the RLS-type model to be robust to a wide range of DGPs, while the other forecasting methods are not and can produce very poor forecasts for some DGPs. Overall, the evidence again indicates that important gains in forecast accuracy can be obtained using our forecasting models and that they are robust since in no case do they perform substantially worse than popular forecasting methods. 8 Conclusion We propose a forecasting framework based on modeling the parameters as random level shift processes dictated by a Bernoulli process for the occurrence of shifts and a normal random variable for its magnitude. Some or all of the parameters of the model can be allowed to change and the latent variables dictating the changes can be common or different for each parameters. Also, the variance of the errors may change in a similar manner. To improve the forecasting performance, we allow the probability of shifts to be a function of some covariates which can be forecasted and incorporate a mean-reversion mechanism. Our model 20 can be cast into a conditional linear and Gaussian state space framework and we provide a computationally efficient method of estimation relying on recent developments on mixture Kalman filtering methods. Simulations show that the proposed model has robust forecasting performance even under model misspecification. We applied our forecasting model to the equity premium and the Treasury bill interest rates and consider different forecasting subsamples or periods. The results show clear gains in forecasting accuracy, sometimes by a very wide margin. Finally, note that given the availability of the proper code for estimation and forecasting, the method is very flexible and easy to implement. For a given forecasting model, all that is required by the users are: 1) which parameters (including the variance of the errors if desired) are subject to change; 2) whether the same or different latent Bernoulli processes dictate the timing of the changes in each parameters; 3) which covariates are potential explanatory variables to model the probability of shifts. 21