Robust time series models with trend and seasonal components

Caivano, Michele; Harvey, Andrew; Luati, Alessandra Article Robust time series models with trend and seasonal components SERIEs - Journal of the Spanish Economic Association Provided in Cooperation with: Spanish Economic Association Suggested Citation: Caivano, Michele; Harvey, Andrew; Luati, Alessandra (2016) : Robust time series models with trend and seasonal components, SERIEs - Journal of the Spanish Economic Association, ISSN 1869-4195, Springer, Heidelberg, Vol. 7, Iss. 1, pp. 99-120, https://doi.org/10.1007/s13209-015-0134-1 This Version is available at: https://hdl.handle.net/10419/158552 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/4.0/ SERIEs (2016) 7:99–120 DOI 10.1007/s13209-015-0134-1 ORIGINAL ARTICLE Robust time series models with trend and seasonal components Michele Caivano1·Andrew Harvey2· Alessandra Luati3 Received: 13 January 2015 / Accepted: 3 November 2015 / Published online: 10 December 2015 © The Author(s) 2015. This article is published with open access at Springerlink.com Abstract We describe observation driven time series models for Student-t and EGB2 conditional distributions in which the signal is a linear function of past values of the score of the conditional distribution. These specifications produce models that are easy to implement and deal with outliers by what amounts to a soft form of trimming in the case of t and a soft form of Winsorizing in the case of EGB2. We show how a model with trend and seasonal components can be used as the basis for a seasonal adjustment procedure. The methods are illustrated with US and Spanish data. Keywords Fat tails ·EGB2 ·Score ·Robustness ·Student’s t ·Trimming · Winsorizing JEL Classification C22 ·G17 1 Introduction Time series are often subject to observations that, when judged by the Gaussian yardstick, are outliers. This is a very real issue for many economic time series. Agustin Maravall’s seasonal adjustment program, TRAMO-SEATS, which he developed jointly with Victor Gomez, tackles the problem by identifying outliers and, where appropriate, replacing them by dummy variables. Here we take a different approach BAndrew Harvey [email protected] 1Bank of Italy, Rome, Italy 2University of Cambridge, Cambridge, UK 3University of Bologna, Bologna, Italy 123 100 SERIEs (2016) 7:99–120 in which we employ a new class of robust models where the dynamics of the level, or location, are driven by the score of the conditional distribution of the observations. These called dynamic conditional score (DCS) models have recently been developed by Creal et al. (2011,2013) and Harvey (2013). They are relatively easy to implement and their form facilitates the development of a comprehensive and relatively straightforward theory for the asymptotic distribution of the maximum likelihood estimator. The changing level of a Gaussian time series is usually obtained from an ARIMA process or explicitly modeled as an unobserved component. The statistical treatment of linear Gaussian unobserved components models is straightforward, with the Kalman filter playing a key role. Additive outliers may be captured by dummy variables. A different way forward is to let the noise have a Student t-distribution, thereby accommodating the outliers. However, the treatment of such a model requires computationally intensive procedures, as described in Durbin and Koopman (2012). The DCS-t model proposed by Harvey and Luati (2014) provides an alternative approach which is observation-driven in that the conditional distribution of the observations is specified. A model of this kind may be compared and contrasted with the methods in the robustness literature; see Maronna, Martin and Yohai (2006, ch 8) and McDonald and Newey (1988), where a parametric approach is called ‘partially adaptive’. Robust procedures for guarding against additive outliers typically respond to large observations in one of two ways: either the response function converges to a positive (negative) constant for observations tending to plus (or minus) infinity or it goes to zero. These two procedures are usually classified as Winsorizing or as trimming. The score for a t-distribution converges to zero and so can be regarded as a parametric form of trimming. Similarly a parametric form of Winsorizing is given by the exponential generalized beta distribution of the second kind (EGB2) distribution. The article by Caivano and Harvey (2014) sets out the theory for the DCS location model with an EGB2 distribution and illustrates its practical value. The article is organized as follows. Section 2reviews the idea behind the DCS location model and expands on the reason for using the conditional score. The various classifications of distributions in terms of their tails are then set out and this is followed by a discussion of the score as it relates to location and scale and the link with robust estimation. This material leads on to the contrast between the DCS t and EGB2 models in Sect. 3. Section 4compares the fit of EGB2 and t-distributions for two US macroeconomic time series. A comparison between the way in which Gaussian and DCS models adapt to structural breaks is made in Sect. 5. This issue is important, because it may be thought that the price paid by DCS models for their robustness is a slow response to structural breaks. DCS models with trend and seasonal components are described in Sect. 6. These models can be regarded as robust counterparts to the unobserved components ‘basic structural model’ (BSM). Seasonal adjustment can be carried out with the BSM by extracting smoothed components using the standard Kalman filter and smoother. Because DCS models only give filtered components, it is necessary to devise a method for smoothing. In Sect. 7we apply this method to a monthly series of tourists arriving in Spain and compare the extracted trend and seasonal components with those obtained by fitting a BSM with the outliers handled by dummy variables. 123 SERIEs (2016) 7:99–120 101 2 Filters, heavy tails and robust estimation The first sub-section below sets out a simple unobserved components model and shows how the innovations form of the Kalman filter may be adapted to form a DCS model. The second sub-section provides a rationale for the use of the conditional score. The way in which tails of distributions may be classified is reviewed in the third subsection and in the fourth, tail behaviour is related to the considerations of robustness. The treatment of these topics is more general and integrated than in Harvey (2013). 2.1 Unobserved components and filters A simple Gaussian signal plus noise model is yt=μt+εt,ε t∼NID0,σ2 ε,t=1,...,T, μt+1=φμt+ηt,η t∼NID0,σ2 η,(1) where the irregular and level disturbances, εtand ηtrespectively, are mutually independent and the notation NID0,σ2denotes normally and independently distributed with mean zero and variance σ2. The autoregressive parameter is φ, while the signalnoise ratio, q=σ2 η/σ2 ε,plays the key role in determining how observations should be weighted for prediction and signal extraction. The reduced form (RF) of (1)isan ARMA(1,1) process yt=φyt−1+ξt−θξt−1,ξ t∼NID0,σ2,t=1, ..., T(2) but with restrictions on θ. For example, when φ=1, 0 ≤θ≤1. The UC model in (1) is effectively in state space form and, as such, it may be handled by the Kalman filter (KF); see Harvey (1989). The parameters φand qmaybeestimated by maximum likelihood (ML), with the likelihood function constructed from the onestep ahead prediction errors. The KF can be expressed as a single equation which combines μt|t−1,the optimal estimator of μtbased on information at time t−1,with ytin order to produce the best estimator of μt+1. Writing this equation together with an equation that defines the one-step ahead prediction error, vt,gives the innovations form of the KF: yt=μt|t−1+vt, μt+1|t=φμt|t−1+ktvt.(3) The Kalman gain, kt,depends on φand q. In the steady-state, ktis constant. Setting it equal to κin (3) and re-arranging gives the ARMA model (2) with ξt=vtand φ−κ=θ. When the noise in (1) comes from a heavy-tailed distribution such as Student’s tit can give rise to observations which, when judged against the yardstick of a Gaussian 123 102 SERIEs (2016) 7:99–120 distribution, are additive outliers. As a result fitting a Gaussian model is inefficient and may even yield estimators which are inconsistent. Simulation methods, such as Markov chain Monte Carlo (MCMC) and particle filtering, provide the basis for a direct attack on such non-Gaussian models; see Durbin and Koopman (2012). However, simulationbased estimation can be time-consuming and subject to a degree of uncertainty. In addition the statistical properties of the estimators are not easy to establish. The DCS approach begins by writing down ft(yt), the distribution of the t-th observation conditional on past observations. The time-varying parameter is then updated by a suitably defined filter. Such a model is said to be observation driven. In a linear Gaussian UC model, the KF depends on the one step-ahead prediction error. The main ingredient in the DCS filter for non-Gaussian distributions is the replacement of vtin the KF equation by a variable, ut,that is proportional to the score of the conditional distribution; compare Maronna, Martin and Yohai (2006, p. 272–4) and the references therein. Thus the second equation in (3) becomes μt+1|t=φμt|t−1+κut,(4) where ut=∂ln ft(yt)/∂μt|t−1and κis treated as an unknown parameter. This filter could be regarded as an approximation to the computer intensive solution for the parameter driven unobserved components model. The attraction of regarding it as a model in its own right is that it becomes possible to derive the asymptotic distribution of the maximum likelihood estimator and generalize in various directions. 2.2 Why the Score? Suppose that at time t−1wehave θt−1the ML estimate of a parameter, θ. Then the score is zero, that is ∂ln Lt−1(θ) ∂θ = t−1  j=1 ∂ln j(θ) ∂θ =0atθ= θt−1,(5) where j(θ;yj)=f(yj;θ). When a new observation becomes available, a single iteration of the method of scoring gives  θt= θt−1+1 It θt−1∂ln Lt(θ) ∂θ = θt−1+1 It θt−1∂ln t(θ) ∂θ , where It( θt−1)=t.I( θt−1)is the information matrix for tobservations and the last expression follows because of (5). For a Gaussian distribution a single update goes straight to the ML estimate at time t(recursive least squares). Remark 1 We may sometimes be able to choose a link function so that the information quantity does not depend on θ. 123 SERIEs (2016) 7:99–120 103 As t→∞,It( θt−1)→∞so the recursion becomes closed to new information. If it is thought that θchanges over time, the filter needs to be opened up. This may be done by replacing 1/tby a constant, which may be denoted as κ. Thus  θt= θt−1+κ1 I θt−1∂ln t(θ) ∂θ . With no information about how θmight evolve, the above equation might be converted to the predictive form by letting  θt+1|t= θt.Thus  θt+1|t= θt|t−1+κ1 I θt|t−1∂ln t(θ) ∂θ .(6) For a Gaussian distribution in which θis the mean and the variance is known to be σ2, 1 I θt|t−1∂ln t(θ) ∂θ =yt− θt|t−1 and (6) is an exponentially weighted moving average (EWMA). If there is reason to think that the parameter tends to revert to an underlying level, ω, the updating scheme might become  θt+1|t=ω(1−φ) +φ θt|t−1+κ1 I θt|t−1∂ln t(θ) ∂θ where |φ|<1.This scheme corresponds to a first-order autoregressive process. More generally we might introduce lags so as to smooth out the changes or allow for periodic effects. 2.3 Heavy tails The Gaussian distribution has kurtosis of three and a distribution is said to exhibit excess kurtosis, or to be leptokurtic, if its kurtosis is greater than three. Although some researchers take excess kurtosis as defining heavy tails, it is not, in itself, an ideal measure, particularly for asymmetric distributions. Most classifications in the insurance and finance literature begin with the behaviour of the upper tail for a nonnegative variable, or one that is only defined above a minimum value; see Embrechts et al. (1997). The two which are relevant here are as follows. A distribution is said to be heavy-tailed if lim y→∞exp(y/α)F(y)=∞ for all α>0,(7) 123 104 SERIEs (2016) 7:99–120 where F(y)=Pr(Y>y)=1−F(y)is the survival function. When yhas an exponential distribution, F(y)=exp(−y/α), so exp(y/α)F(y)=1 for all y.Thus the exponential distribution is not heavy-tailed. A distribution is said to be fat-tailed if, for a fixed positive value of η, F(y)=cL(y)y−η,η>0,(8) where cis a non-negative constant and L(y)is slowly varying,1that is lim y→∞ L(ky) L(y)=1. The parameter ηis the tail index. The implied PDF is a power law PDF f(y)∼cL(y)ηy−η−1,y→∞,η>0,(9) where ∼is defined such that a(x)∼b(x)as x→x0if limx→x0(a/b)→1.The m-th moment exists if m<η. The Pareto distribution is a simple case in which F(y)=y−η for y>1.If a distribution is fat-tailed then it must be heavy-tailed, but the converse is not true; see Embrechts, Kluppelberg and Mikosch (1997, p. 41–2). The above criteria are related to the behavior of the conditional score and whether or not it discounts large observations. This, in turn, connects to robustness, as shown in the sub-section following. More specifically, consider a power law PDF, (9), with ydivided by a scale parameter,ϕ, so that F(y/ϕ) =cL(y/ϕ)(y/ϕ)−ηand f(y)∼ cL(y)ϕ−1η(y/ϕ)−η−1.Then ∂ln f/∂ϕ ∼η/ϕ as y →∞ (10) and so the score is bounded. With the exponential link function, ϕ=exp(λ), ∂ln f/∂λ ∼ηas y→∞.Similarly as y→0, ∂ln f/∂λ ∼η. The logarithm of a variable with a fat-tailed distribution has exponential tails. Let xdenote a variable with a fat-tailed distribution in which the scale is written as ϕ=exp(μ) and let y=ln x.Then for large y f(y)∼cL(ey)ηe−η(y−μ),η>0,as y →∞, whereas as y→−∞,f(y)∼cL(ey)ηeη(y−μ),η>0.Thus yis not heavy-tailed, but it may exhibit excess kurtosis. The score with respect to location, μ, is the same as the original score with respect to the logarithm of scale and so tends to ηas y→∞. 1More generally regularly varying is limy→∞(L(ky)/L(y)) =kβ;see Embrechts, Kluppelberg and Mikosch (1997, p. 37, 564). Fat-tailed distributions are regularly varying with η=−β>0. 123 SERIEs (2016) 7:99–120 105 2.4 Robust estimation The location-dispersion model is yt=μ+ϕεt,t=1, ..., T,(11) where εtis a standardized variable with PDF denoted exp ρ(εt)and the scale, ϕ,is called the dispersion for yt;see Maronna, Martin and Yohai (2006, p. 37–8). The density for ytis f(yt;μ, ϕ, ξ) =ϕ−1exp ρ((yt−μ)/ϕ), where ξdenotes one or more shape parameters, and the scores for μand ϕare given by differentiating ln f(yt)=ρ((yt−μ)/ϕ) −ln ϕ. The score for location is ∂ln ft ∂μ =∂ρ(zt) ∂μ =ψL(zt), where zt=(yt−μ)/ϕ, whereas the score for scale is ∂ln ft ∂ϕ =∂ρ(zt) ∂ϕ −1 ϕ=ψS(zt)−1 ϕ. Note that ψS(zt)=ϕ−1ztψL(zt). When the scale is parameterized with an exponential link function, ϕ=exp λ, thescoreis ∂ln f ∂λ =∂ρ(zt) ∂λ −1=ztψL(zt)−1.(12) If yt=ln xt,where xt=ϕεtas in (11) and then the logarithm of the scale parameter for xt,that is ln ϕ, becomes the location for yt.Hence ψL(yt)=ψS(xt). The ML estimators are asymptotically efficient, assuming certain regularity conditions hold. More generally ρ(.) may be any function deemed to yield estimators with good statistical properties. In particular, the estimators should be robust to observations which would be considered to be outliers for a normal distribution. When normality is assumed, the ML estimators of the mean and variance are just the corresponding sample moments, but these can be subject to considerable distortion when outliers are present. Robust estimators, on the other hand, are resistant to outliers while retaining relatively high efficiency when the data are from a normal distribution. The M-estimator, which features prominently in the robustness literature, has a Gaussian response until a certain threshold, K,whereupon it is constant; see Maronna, Martin and Yohai (2006, p. 25–31). This is known as Winsorizing as opposed to trimming, where observations greater than Kin absolute value are given a weight of zero.2 2In both cases a (robust) estimate of scale needs to be pre-computed and the process of computing M-estimates is then often iterated to convergence. 123 106 SERIEs (2016) 7:99–120 3 DCS location models The stationary first-order DCS model corresponds to the Gaussian innovations form, (3), and is yt=μt|t−1+vt=μt|t−1+exp(λ)εt,t=1, ..., T, μt+1|t=δ+φμt|t−1+κut,|φ|<1,(13) where ω=δ/(1−φ) is the unconditional mean of μt|t−1,εtis a serially independent variate with unit scale3and utis proportional to the conditional score, that is ut= k.∂ ln f(yt|yt−1,yt−2, ...)/∂μt|t−1,where kis a constant. More generally, an ARMA-type model of order (p,r)is μt+1|t=δ+φ1μt|t−1+... +φpμt−p+1|t−p+κ0ut+κ1ut−1+... +κrut−r. (14) In the Gaussian case ut=yt−μt|t−1and if qis defined as max(p,r+1), then ytis an ARMA(p,q)process with MA coefficients θi=φi−κi−1,i=1, .., q.Nonstationary ARIMA-type models may also be constructed as may structural times series models with trend and seasonal components. Explanatory variables can be introduced into DCS models, as described in Harvey and Luati (2014). Maronna, Martin and Yohai (2006, Sect 8.6 and 8.8) give a robust algorithm for AR and ARMA models with additive outliers. For a first-order model their filter is essentially the same as (13) except that their dynamic equation is driven by a robust ψ−function and they regard the model as an approximation to a UC model.4 3.1 Student tmodel The Student tνdistribution has fat tails for finite degrees of freedom, ν, with the tail index given by ν. Moments exist only up to and including ν−1. The excess kurtosis, that is the amount by which the normal distribution’s kurtosis of three is exceeded, is 6/(ν −4), provided that ν>4. When the location changes over time, it may be captured by a model in which, conditional on past observations, ythas a tν-distribution with μt|t−1generated by a linear function of ut=1+ν−1e−2λyt−μt|t−12−1 vt,t=1, ..., T,(15) 3The standard deviation is √ν/(ν −2)times the scale. 4Muler, Peña and Yohai (2009, p. 817) note two shortcomings of the estimates obtained in this way. They write: ‘First, these estimates are asymptotically biased. Second, there is not an asymptotic theory for these estimators, and therefore inference procedures like tests or confidence regions are not available.’ They then suggest a different approach and show that it allows an asymptotic theory to be developed. 123 SERIEs (2016) 7:99–120 113 μt=μt−1+βt−1+ηt,η t∼NID0,σ2 η, βt=βt−1+ζt,ζ t∼NID0,σ2 ζ, where ζtand ηtare independent of each other and of εt.The random walk plus noise or local level model is a special case. As regards the seasonal component, let γjt denote the effect of season jat time tand define γt=(γ1t, ..., γst). The full set of seasonals evolves as a multivariate random walk γt=γt−1+ωt,t=1, ...., T, where ωt=(ω1t, ..., ωst)is a zero mean disturbance with Var(ωt)=σ2 ωI−s−1ii, with σ2 ω>0.Although all sseasonal components are continually changing, only one affects the observations at any particular time, that is γt=γjt when season jis prevailing at time t.The requirement that the seasonal components always sum to zero is enforced by the restriction that the disturbances sum to zero at each t.This restriction is implemented by the correlation structure of ωt, where Variωt=0, coupled with initial conditions constraining the seasonals to sum to zero at t=0. The multiplicative seasonal ARIMA model known as the ‘airline model’ is syt=(1+θL)1+Lsξt,ξ t∼NID0,σ2 ξ Box and Jenkins (1976, pp. 305–6) gave a rationale for this model in terms of EWMAs at monthly and yearly intervals. The reduced form of the BSM is syt∼ MA(s+1).Maravall (1985), compared the autocorrelation functions of sytfor the BSM and airline model for some typical values of the parameters and finds them to be quite similar, particularly when the seasonal MA parameter, , is close to minus one. In the limiting case when is equal to minus one, the airline model is equivalent to a BSM with σ2 ζ=σ2 ω=0. 6.2 Stochastic trend and seasonal in the DCS model The DCS model for trends and seasonals, yt=μt|t−1+γt|t−1+vt,t=1, ..., T, has a structure which is similar to that of the innovations form of the Kalman filter for the BSM. The filter for the trend is μt+1|t=μt|t−1+βt|t−1+κ1ut βt+1|t=βt|t−1+κ2ut. An integrated random walk (IRW) trend in the UC local linear trend model implies the constraint κ2=κ2 1/(2−κ1), 0<κ 1<1,which may be found from Harvey 123 114 SERIEs (2016) 7:99–120 (1989, p. 177). The restriction can be imposed on the DCS-tmodel by treating κ1=κ as the unknown parameter, but without unity imposed as an upper bound. The filter for the seasonal is γt|t−1=z tγt|t−1,γ t+1|t=γt|t−1+κtut, where the s×1 vector ztpicks out the current season from the vector γt|t−1.If κjt, j=1, .., s,denotes the jth element of κt, then in season jwe set κjt =κs,where κsis a non-negative unknown parameter, whereas κit =−κs/(s−1),i= j,i=1, .., s. The amounts by which the seasonal effects change therefore sum to zero. The initial conditions at time t=0 are estimated by treating them as parameters. The above filter may be regarded as a robust version of the well-known Holt-Winters filter; see Harvey (1989, p. 31). However, it differs from Holt-Winters in the Gaussian case by enforcing the restriction that the seasonals sum to zero. This is an important advantage. 6.3 Seasonal adjustment In contrast to the Gaussian BSM, the DCS model has no exact solution for smoothing. Some possibilities are suggested in Harvey (2013, Sect 3.7), but these are difficult to generalize beyond the local level model. The best way to employ the DCS model for seasonal adjustment is to use it to mitigate the effects of outliers by modifying them rather than eliminating them by dummy variables. A dummy variable effectively means that the corresponding observation is treated as though it were missing; in other words it corresponds to hard trimming. We first fit the DCS model and use it to construct pseudo-observations from the signal, that is μs t|t−1=μt|t−1+γt|t−1.Thus ytμs t|t−1=μs t|t−1+ut,t=1, ..., T.(18) For the DCS-t model we can write ytμs t|t−1=(1−bt)yt+btμs t|t−1,t=1, ..., T, where btis as in (16) with μt|t−1replaced by μs t|t−1.In the EGB2, we work directly with (18). The pseudo-observations are now used to estimate the parameters in a BSM and the signal is estimated by smoothing. This new signal, denoted μs t|T,is then used to construct new pseudo-observations as ytμs t|T=μs t|T+ut,t=1, ..., T,(19) 123 SERIEs (2016) 7:99–120 115 with btreplaced by bt|T=bt|Tμs t|T=yt−μs t|T2 /ν exp(2λ) 1+yt−μs t|T2 /ν exp(2λ) for the t-distribution and by bt|T= exp yt−μs t|T/h/σ 1+exp yt−μs t|Th/σ for the EGB2. The BSM may be re-estimated and the whole process iterated until convergence. If desired, λand σcan be updated at each step using the sample variance of the u ts. The relationship between variance of the u tsand λfor the t-distribution is given in Harvey (2013, p. 62). For the EGB2 (with ξ=ς) it follows from Caivano and Harvey (2014) that σ2 u=σ2h2ξ2/(2ξ+1). The seasonally adjusted observations are constructed as ys t=yt−γt|T=μt|T+vt|T,t=1, ..., T,(20) where vt|T=yt−μs t|T.If it is felt that the outliers should be modified to make them less extreme, we could let ys t=μt|T+ut|T. The above procedure could be implemented with TRAMO-SEATS rather than the unobserved components BSM. 7 Tourists in Spain The logarithm of the number of tourists entering Spain from January 2000 to April 2014 (source: Frontur) is plotted in Fig. 5. Comparing a Gaussian unobserved components model with dummy variables with a DCS-t model provides a contrast between hard and soft trimming. Fitting a BSM gives a Bowman-Shenton statistic of 19.29 for the residuals so normality is clearly rejected. (The 1% critical value for a χ2 2is 9.21.) The automatic outlier detection option in STAMP finds five outliers, three of which are in 20012, along with a small structural break which, although not sudden, is attributed to April 2008. The corresponding (edited) output is shown below. Fitting a BSM with dummy variables included for the outliers and break gives the smoothed components shown in Fig. 5. Date Estimate SE ’t-stat’ Outlier 2002(3) 0.14295 0.02818 5.07283 Outlier 2002(4) -0.14438 0.02824 -5.11205 Outlier 2002(8) 0.10473 0.02775 3.77429 123 116 SERIEs (2016) 7:99–120 Ltourists Level+Intv 2000 2005 2010 14.5 15.0 15.5 16.0 Ltourists Level+Intv Ltourists-Level 2000 2005 2010 15.1 15.2 15.3 15.4 15.5 Ltourists-Level Ltourists-Seasonal 2000 2005 2010 -0.50 -0.25 0.00 0.25 0.50 Ltourists-Seasonal Ltourists-Irregular 2000 2005 2010 -0.01 0.00 0.01 0.02 Ltourists-Irregular Fig. 5 Smoothed estimates of components from BSM with outliers treated by dummy variables Outlier 2005(4) -0.10592 0.02766 -3.82885 Outlier 2010(4) -0.12736 0.02767 -4.60211 Level break 2008(4) -0.10926 0.02429 -4.49895 The ML estimates for the parameters of a DCS-t random walk model (the drift was not significant) with a seasonal component are as follows: κ=0.4906(0.0666)κs=1.0068 (0.1172) ν=6.006 (1.5678) λ=−3.2573 (0.0612) with initial values μ0=15.1018(0.0174)and γ0=[−0.5795(0.0270), −0.4761 (0.0269),−0.1780(0.0264), 0.1733(0.0244), 0.1475(0.0291), 0.2648(0.0261), 0.5791(0.0259),0.4566(0.0209),0.3157(0.0327), 0.1167(0.0304), −0.3862(0.0289), −0.4339], for the seasonal factors. The figures in parentheses are numerical standard errors and the initial value for the last seasonal factor is  γ0,12 =0.018 which is constructed by from the others and thus has no standard error. The filtered DCS-t trend shown in Fig. 6appears not to be affected by the outliers, which are downweighted, and the trend, shown in the two top panels, adjusts to the fall in 2008 as quickly as the Gaussian filter from the outlier-adjusted BSM which is fitted with a dummy for a break in 2008(4). The figure also shows the seasonal component and contrasts the prediction errors (irregular component) and scores; it is evident that extreme observations are downweighted by the score. The autocorrelation functions 123 SERIEs (2016) 7:99–120 117 Ltourists mu_nu_6 2000 2005 2010 2015 14.5 15.0 15.5 16.0 Ltourists mu_nu_6 mu_nu_6 2000 2005 2010 2015 15.1 15.2 15.3 15.4 mu_nu_6 seas_nu_6 2000 2005 2010 2015 -0.50 -0.25 0.00 0.25 0.50 seas_nu_6 u_nu_6 v_nu_6 2000 2005 2010 2015 -0.1 0.0 0.1 u_nu_6 v_nu_6 Fig. 6 Filtered estimates from DCS-t filter v_nu_6 2000 2005 2010 2015 -0.1 0.0 0.1 v_nu_6 u_nu_6 2000 2005 2010 2015 -0.1 0.0 0.1 0.2 u_nu_6 ACF-v_nu_6 0 5 10 15 20 -0.5 0.0 0.5 1.0 ACF-v_nu_6 ACF-u_nu_6 0 5 10 15 20 -0.5 0.0 0.5 1.0 ACF-u_nu_6 Fig. 7 Residuals and scores, together with ACFs of irregular and score are shown in Fig. 7. The autocorrelations of the scores tend to be slightly bigger in absolute value, reflecting the fact that they are not weakened by outliers in the same way as the raw predictions errors. Figure 8shows the filtered trend together with the smoothed trends for the first three iterations, obtained as described in Sect. 6.3. The trend at the second iteration is almost indistinguishable from the one obtained at the third. 123 118 SERIEs (2016) 7:99–120 mu_nu_6 Level_2 Level_1 Level_3 2000 2005 2010 2015 15.10 15.15 15.20 15.25 15.30 15.35 15.40 mu_nu_6 Level_2 Level_1 Level_3 Fig. 8 Smoothed trends from DCS-t model Level - DCS Level - TRAMO-SEATS 2000 2005 2010 2015 15.10 15.15 15.20 15.25 15.30 15.35 15.40 Level - DCS Level - TRAMO-SEATS Fig. 9 TRAMO-SEATS trend and TRAMO-SEATS-DCS trend The above procedure was repeated, but using TRAMO-SEATS to obtain a smoothed estimate of the trend from the DCS pseudo-observations. The result is shown in Fig. 9, together with the smoothed trend obtained from the full (automatic) TRAMO-SEATS procedure As can be seen, the two trends are very close. Figure 10 shows the corresponding smoothed irregular components. 123 SERIEs (2016) 7:99–120 119 Irregular - DCS Irregular - TRAMO-SEATS 2000 2005 2010 2015 -0.075 -0.050 -0.025 0.000 0.025 0.050 0.075 Irregular - DCS Irregular - TRAMO-SEATS Fig. 10 Irregular components for TRAMO-SEATS and TRAMO-SEATS-DCS We also fitted an EGB2 model, but the fit was not as good as for the DCS-t model. Generally the results were similar except that the seasonal pattern for EGB2 was found to be deterministic. The ML estimates for the asymmetric model were: κ=0.32(0.07)κs=0  ξ=0.74 (0.55)ς=0.66 (0.55) λ=−3.03 (0.07) The two shape parameters are very close and there is virtually no loss in imposing symmetry, that is ξ=ς. 8 Conclusions This article has shown how DCS models with changing location and/or scale can be successfully extended to cover EGB2 conditional distributions. Most of the theoretical results on the properties of DCS-t models, including the asymptotic distribution of ML estimators, carry over to EGB2 models. However, whereas the t-distribution has fattails, and hence subjects extreme observations to a form of soft trimming, the EGB2 distribution has light tails ( but excess kurtosis) and hence gives a gentle form of Winsorizing. The examples show that the EGB2 distribution can give a better fit to some macroeconomic series. The way in which DCS models respond to breaks was examined and it was shown that, contrary to what might be expected, they adjust almost as rapidly as Gaussian models. A seasonal adjustment procedure may be carried out with DCS models that include both trend and seasonal components. Having fitted the DCS model, the scores are used 123 120 SERIEs (2016) 7:99–120 to adjust the data before smoothing using a standard Gaussian model. Two or three iterations seem to be sufficient. The method was illustrated with data on tourists in Spain. There is a case for a structural break in April, 2008, but the DCS model quickly adjusts to it and indeed it could be reasonably argued that it is better to let the change take place over several months rather than assigning it to just one. In summary, our new DCS procedure, like TRAMO-SEATS, provides a practical approach to seasonal adjustment in the presence of outliers. Acknowledgments We are grateful to the Gabriele Fiorentini and a referee for helpful comments on the first draft. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References Caivano M, Harvey AC (2014) Time series models with an EGB2 conditional distribution. J Time Ser Anal 34:558–571 Creal D, Koopman SJ, Lucas A (2013) Generalized autoregressive score models with applications. J Appl Econ 28:777–795 Creal D, Koopman SJ, Lucas A (2011) A dynamic multivariate heavy-tailed model for time-varying volatilities and correlations. J Bus Econ Stat 29:552–563 Durbin J, Koopman SJ (2012) Time series analysis by state space methods, 2nd edn. Oxford Statistical Science Series, Oxford Embrechts P, Kluppelberg C, Mikosch T (1997) Modelling extremal events. Springer, Berlin Fernandez C, Steel MFJ (1998) On Bayesian modeling of fat tails and skewness. J Am Stat Assoc 99:359– 371 Harvey AC (1989) Forecasting, structural time series models and the kalman filter. Cambridge University Press, Cambridge Harvey AC (2013) Dynamic models for volatility and heavy tails., Econometric Society MonographCambridge University Press, Cambridge, New York Harvey AC, Luati A (2014) Filtering with heavy tails. J Am Stat Assoc 109:1112–1122 Koopman SJ, Harvey AC, Doornik JA, Shephard N (2009) STAMP 8.2 structural time series analysis modeller and predictor. Timberlake Consultants Ltd, London Maravall A (1985) On structural time series models and the characterization of components. J Bus Econ Stat 3:350–355 Maronna R, Martin D, Yohai V (2006) Robust statistics: theory and methods. Wiley, Chichester McDonald JB, Newey WK (1988) Partially adaptive estimation of regression models via the generalized t distribution. Econ Theory 4:428–457 Muler N, Pena D, Yohai VJ (2009) Robust estimation for ARMA models. Ann Stat 37:816–840 Zhu D, Galbraith JW (2011) Modelling and forecasting expected shortfall with the generalised asymmetric student-t and asymmetric exponential power distributions. Journal of Empir Financ 18:765–778 123