key: cord-0155582-dnnke7fh authors: Dokumentov, Alexander; Hyndman, Rob J. title: STR: Seasonal-Trend Decomposition Using Regression date: 2020-09-13 journal: nan DOI: nan sha: ec659020395e69828400518f88baec5b2bfc7778 doc_id: 155582 cord_uid: dnnke7fh We propose a new method for decomposing seasonal data: STR (a Seasonal-Trend decomposition using Regression). Unlike other decomposition methods, STR allows for multiple seasonal and cyclic components, covariates, seasonal patterns that may have non-integer periods, and seasonality with complex topology. It can be used for time series with any regular time index including hourly, daily, weekly, monthly or quarterly data. It is competitive with existing methods when they exist, but tackles many more decomposition problem than other methods allow. STR is based on a regularized optimization, and so is somewhat related to ridge regression. Because it is based on a statistical model, we can easily compute confidence intervals for components, something that is not possible with most existing decomposition methods (such as STL, X-12-ARIMA, SEATS-TRAMO, etc.). Our model is implemented in the R package stR, so can be applied by anyone to their own data. Author: STR: Seasonal-Trend decomposition using Regression Article accepted for INFORMS Journal on Data Science 3 A different approach was followed by Cleveland et al. (1990) , who developed STL (Seasonal-Trend decomposition using Loess), which has come to be widely used outside the national statistics agencies, partly because of its availability in R (R Core Team 2020). This method uses iterative Loess smoothing to obtain an estimate of the trend, and then Loess smoothing again to extract a changing additive seasonal component. Several model-based methods for seasonal decomposition have been developed, including the TRAMO/SEATS procedure which was developed at the Bank of Spain (Gómez and Maravall 2001) , the TBATS model of De Livera et al. (2011) , and various structural time series model approaches (Harvey 1990 , Commandeur et al. 2011 ). One big advantage of using a model for seasonal decomposition and adjustment is that it provides a natural way to compute confidence and prediction intervals. Despite this long history, and the availability of many algorithms and models for time series decomposition, there are many time series characteristics that are not addressed in these approaches. We seek a time series decomposition method with the following attributes: • provides a meaningful and simple statistical model; • allows for computation of confidence intervals of components; • can take account of covariates; • allows for fractional seasonal periods (e.g., with weekly data); • allows for multiple seasonal periods (e.g., with daily or sub-daily data); • allows for complex seasonal topology (e.g., due to public holidays); • can incorporate covariates which interact with the seasonality. None of the existing methods satisfies all of these requirements. We aim to fill this gap with our new approach, which is clear, general, model-based, robust (if required), and simple. We show that the problem of seasonal decomposition can be re-cast in the framework of ordinary least squares or quantile regression. Moreover, our approach provides new features (such as covariates which affect data in a complex seasonal manner, and the ability to model complex seasonality) that have not been developed before. Our new STR method is the most general framework for the decomposition of seasonal data that is available at present. We suppose our time series y t can be decomposed additively as follows: where • T t is a smoothly changing trend; are smoothly changing seasonal components with possibly complex topology; • z p,t are covariates with coefficients φ p,t which may be time-varying and even seasonal; • R t is the "remainder". The seasonal component is assumed to have a repeating pattern which changes very slowly or is constant over time. The trend component describes the smooth underlying mean of the data (after accounting for seasonality). The remainder component (which we will assume to be uncorrelated) contains only noise and idiosyncratic patterns in the data. Where components behave multiplicatively, we can use a log or Box-Cox transformation (Box and Cox 1964) to made the additive assumption more reasonable. The total number of coefficients to be estimated is usually much larger than the number of observations. Consequently, we will impose some regularization on the coefficient estimates. In fact, without some regularization, the trend, seasonal components and time-changing coefficients will be unindentifiable. The innovation in our approach is to provide an estimation method that allows these regularizations to be imposed in an efficient manner using matrix differencing operators, which convert the estimation method to a linear model. In this way, the estimation becomes analogous to ridge regression. Although the matrices involved can become very large, they are also sparse, so sparse matrix algebra methods can be used to reduce the computational burden. Smoothness of a function implies that its second derivatives are small, since otherwise large values of the second derivatives would lead to high curvature of the function making it "wiggly" (see Wood 2006) . Thus, we consider a function smooth if its second derivatives are normally distributed with zero mean and where the variance controls the degree of smoothness. Because time is discrete, we approximate derivatives with differences. Here, we use the notation So a smooth trend is obtained by requiring that , where σ L controls the degree of smoothness. This can be expressed via the multivariate normal density is an n-vector of the trend component; D is the (n − 2) × n difference operator matrix such that D = ∆ 2 T t n t=3 ; and f L (σ L ) is a normalizing function. Note that the model subsumes the case of a linear trend, obtained when σ L = 0 (making the trend's second derivatives equal to zero). Let m i denote the number of "seasons" in the seasonal component S (i) t . For example, if y t is monthly data, there is only one seasonal component (I = 1) with m 1 = 12; if y t is daily data, there may be two seasonal components (I = 2), with a weekly pattern (m 1 = 7) and an annual pattern (m 2 = 365 ignoring leap years). Let us also define the function κ i (t) = t mod m i , which transforms time t into the corresponding season κ i (t) ∈ {1, . . . , m i }. At time t, we observe only one element of each seasonal component. It is reasonable to ask what the other elements of the component are at that moment t. For example, if we have daily data with weekly seasonality, it is valid to ask on a Wednesday, "What is the current value of the seasonal component corresponding to Friday?". In other words, we define latent components that are responsible for seasons other than κ i (t). In this way, we treat the ith seasonal component as a two-dimensional matrix, [S (i) k,t ] (k = 1, . . . , m i ; t = 1, . . . , n). This formulation of each seasonal component as a two-dimensional structure is one of the novelties and features of STR, and it allows for simple and natural smoothing constraints to be applied. The seasonal two-dimensional surface actually has the topology of a cylinder, where dimension t (time) is along the length of the cylinder, but dimension s (season) is "circular" around the cylinder. For example, in the case of a weekly seasonal component in daily data, the season dimension corresponds to days of the week. A time series corresponds to a spiral on the cylinder. The other locations on the cylinder (between the spirals) will represent "imaginary" seasonality. These correspond, for example, to the current value of Friday's seasonality when today is a Wednesday. To ensure they remain identifiable and interpretable, the seasonal terms must sum to zero at any time t, so they must satisfy the property k S (i) k,t = 0 for each t. Smoothness of the weekly seasonal pattern can be described as smoothness of S The seasonal smoothness can also be expressed using difference operator matrices. We define three such matrices, D tt,i , D ss,i and D st,i , corresponding to smoothness in the time, season and time-season directions respectively. For smoothness in the time direction, we require that the seasonal vectors ∆ 2 S Then the density for the second derivatives of the ith seasonal component is given by Smoothness in the season and time-season directions are defined analogously by restricting the derivatives ∂ 2 ∂s 2 and ∂ 2 ∂s∂t . It should be noted that the matrices D st and D ss take the cylinder topology into account when calculating proper "circular" differences. Again, by setting the smoothing parameters to zero, some simple special cases are obtained. Setting the variance of the seasonal component in the time-seasonal direction to zero causes the seasonal component to be periodic. By setting both trend and time-seasonal variances to zero we effectively force our model to have a linear trend with seasonal dummies. Setting the variance of the second derivatives of the seasonal component in the time direction to zero alone makes changes in the seasonal component linear over time, but does not preclude the seasonal pattern from changing. A seasonal two-dimensional surface can have a topology other than a cylinder. Suppose that we are going to model some social behaviour (e.g., electricity demand) during working days and holidays (including weekends). The topology modelling the human behaviour is shown in Figure 1 . The left circle represents a working day, the right circle represents a holiday. They are connected by lines representing transition periods. Points A and C represent hour 0 of a day, and points B and D represent hour 12. Every day has 24 hours, and the transition periods take 12 hours. According to the diagram, a working day can either be followed by another working day, or flow until hour 12 and then undergo a 12 hour transition (line B-C) into a holiday. Similarly, a holiday can either be followed by another holiday, or undergo a 12 hour transition (line D-A) into a working day. Equivalently, the topology shown in Figure 1 can be described as two connected cylinders. The differencing matrices D ss and D st are defined as above for all data points except A, B, C and D. At points A, B, C and D, the second derivatives can be restricted in various ways. One possibility is to define D ss and D st to regularise the derivatives twice, once for each path in the diagram. An example of a decomposition with complex topology can be found in Section 4.2. The remainder terms R t are assumed to be iid N (0, σ 2 R ). Let y = [y 1 , . . . , y n ] be an n-vector of observations, Z = [z t,p ] be a matrix of covariates with coefficients given by where f (σ R ) is a normalizing function. Covariates are often available to help explain some of the variation seen in a time series. For example, Findley et al. (2012) considers the effects of various moving holidays on human activities, and Bell and Martin (2004) considers covariates with time-varying coefficients. We consider time series that are affected by covariates of three types: (1) static covariates where we assume that the associated coefficients are constant over time; (2) flexible covariates with coefficients whose magnitudes change smoothly over time, but where the coefficients do not exhibit any seasonal pattern; and (3) seasonal covariates with coefficients whose magnitudes change smoothly over time in a seasonal pattern. As far as we know, this last type of time-varying coefficient is new in the decomposition of time series. To distinguish the three types of covariates, we write where: • a t,h are static covariates with constant coefficients; • b t,j are flexible covariates with time-varying but non-seasonal coefficients; • c t,k are seasonal covariates with time-varying coefficients, where the coefficients have seasonal patterns with corresponding seasonal periods p k ∈ {m 1 , . . . , m I }; • A = [a t,h ] is an n × H matrix of static covariates, where every covariate occupies a single column; • α is an H-vector of coefficients of the static covariates; • β j is the jth n-vector of changing coefficients for the jth flexible covariate; • γ k is the kth n-vector of changing coefficients for the kth seasonal covariate. To ensure the time-varying coefficients change slowly over time, we will impose smoothness via differencing operators, much as we did for the trend and seasonal components. For this purpose, we define matrices ∆ tt,k , ∆ st,k and ∆ ss,k that take second differences of the kth seasonal coefficients in the time, time-season and season dimensions, respectively. The parameters of the model are given by α, β 1 , . . . , β J , γ 1 , . . . , γ K , , s 1 , . . . , s I , and the standard deviations σ R , σ L , σ 1 , . . . , σ I along with those corresponding to the other smoothness constraints. The model is over-parametrized with more parameters than observations. However, the regularization via the smoothness constraints will allow it to be estimated. Combining the preceding results, and assuming that all terms are independent of each other and that the standard deviations are known, we find that the minus log likelihood function for this where the λ and θ coefficients are ratios of standard deviations. Then maximum likelihood estimates can be obtained by minimizing (2) over Note that (2) corresponds to the minus log likelihood function for the linear model where y + = [y , 0 ] is the vector y padded with zeros, η is a vector of unknown coefficients (with seasonal components adjusted by removing the last row of seasonal observations), ε ∼ N (0, σ 2 R I), and with fixed parameters Λ = {λ tt,1 , λ st,1 , λ ss,1 , . . . , λ tt,r , λ st,r , λ ss,r , λ , θ 1 , . . . , θ J , θ tt,1 , θ st,1 , θ ss,1 , . . . , θ tt,K , θ st,K , θ ss,K } . If some values of Λ are zeros, the corresponding rows of matrix X can (and should) be removed, as they have no effect and removing them improves the computation time. The total number of coefficients (the length of η) that we need to estimate is usually much larger than the number of observations (the length of y). This does not cause computational problems because the coefficients are regularised via smoothness constraints, and the estimation is performed using (4), where y + is longer than η. Using standard linear regression results, the maximum likelihood solution is given bŷ with covariance matrix The trend componentT t and seasonal componentsŜ (1) t , . . . ,Ŝ t , and their corresponding confidence intervals can be obtained directly fromη and Cov(η). We will use cross-validation for estimating the λ and θ parameters of (5). Since the model in (4) is a linear model, the leave-one-out cross-validation residuals can be calculated (see Seber and Lee 2003) using where y i is the ith element of the vector y,ŷ i is the ith element of the vectorŷ = Hy, and h ii is the ith diagonal element of the hat matrix H = X(X X) −1 X . Therefore, we can use the well-known formula for cross-validation in linear regression (see, for example, Ruppert et al. 2003, p.45) : STR finds the optimal λ and θ smoothing parameters by minimising CV. The problem of minimising CV can be complex, for example if there are many local minima, and we have no method that guarantees the finding of the global minimum. However, rather simple methods often work well in practice. We perform such optimisation using the Nelder-Mead method as implemented in the optim() function from the stats package in R (R Core Team 2020). In cases where numerical difficulties may arise in using leave-one-out cross-validation, we resort to K-fold cross-validation where we split the data set into K subsets such that the observation at time t belongs to subset i if and only if (t − 1) mod (Kg) ∈ [ig, . . . , (i + 1)g − 1]. When g = 1, no consecutive observations lie in the same subset of data. This gives a reasonable sparsity of the subsets, and we speculate that the resulting K-fold cross-validation will not differ much from the result of leave-one-out cross-validation. However, in this case, the trend component usually absorbs any serial correlation that may be present in the data. For g > 1, blocks of g consecutive observations are selected within each subset, similar to a blockbootstrap method, but with a greater separation of neighbouring blocks. The value of g determines the amount of correlation that is allowed to remain in the residuals, and therefore controls the trade-off between correlated residuals and trend flexibility. To demonstrate that our method estimates the appropriate time series components, we use simulations where we know the true underlying components. We will compare our results against two other decomposition methods that handle multiple seasonal components, STL and TBATS. STL (Cleveland et al. 1990 ) is a well known decomposition method, and therefore a natural choice to use as a method to compare. The original methodology and implementation was designed for only one seasonal period. However, the implementation in the forecast package for R (Hyndman et al. 2021) will compute an STL decomposition iteratively with multiple seasonal periods. TBATS (De Livera et al. 2011 ) is a forecasting method designed for multiple seasonalities; it provides a decomposition of the time series into components, and so can provide an alternative comparison. Our simulated daily data consists of four components: where T t is a trend, S W t is a weekly seasonal component, S Y t is a yearly seasonal component, R t is the remainder, and α, β, and γ are parameters which control the contribution of the components to y t . We use two data generating processes (DGPs) to simulate data. For the first DGP, we use a quadratic trend function with random coefficients, such that T * t = N 1 (t + n/2(N 2 − 1)) 2 where N 1 and N 2 are independent N(0,1) random variables. Then T * t is normalized to give T t with mean zero and unit variance. The seasonal components consist of five pairs of Fourier terms with random N(0,1) coefficients, which are then also normalized. Because this DGP has deterministic components, we refer to it as the "Deterministic" GDP. For the second "Stochastic" DGP, the trend is given by an ARIMA(0,2,0) model with standard normal errors, which is then normalized. The seasonal components are also smooth and obtained using a similar process constrained to be periodic. Specifically, for each seasonal component we start with a vector of iid N(0,1) values of the length of one season, normalized and replicated to be of length n. This is then integrated and normalized twice to make it smooth. So each seasonal component is as smooth as the trend, is centred on zero and is periodic. For comparison of the three decomposition methods, we use parameter sets (α, β, γ) = We will illustrate our proposed method using two data sets. The first is a simple application with monthly data, that allows us to compare our results with existing methods. The second is more complicated and demonstrates the full range of capabilities of our method including multiple seasonal patterns, complex seasonal topology due to public holidays and interacting seasonalities, with a nonlinear covariate. The time series shown in the top panel of Figure 3 concerns logarithms of monthly supermarket and grocery store revenue in New South Wales, Australia, from 2000 to 2009. We apply the model in (4) to the logged data with no covariates and one seasonal component (m 1 = 12), with the smoothing parameters selected using leave-one-out cross-validation. The resulting decomposition is shown in Figure 3 , along with 95% confidence intervals shown as blue bands. Here, the trend component measures the overall health and size of the economy; the seasonal component allows the study of how human behaviour changes over time; the remainder component allows us to identify unusual features of the data. We illustrate the power of the STR approach by presenting a more complicated example of time series decomposition using half-hourly electricity consumption in the state of Victoria, Australia, during the 115 days starting on 10 January 2000. For each 30-minute period, we also have the air temperature at the Melbourne weather station (near the centre of the largest city in Victoria). We use concurrent temperatures and their squared values as predictors. Because the effect of temperature can change over time (a mild day in summer will have a different effect from a day of the same temperature in winter), we allow for changing coefficients over time. This data set has two seasonal patterns. The first pattern is a weekly seasonal pattern (m 1 = 48 × 7 = 336) that represents the specific demand features that are attributable to a particular day of the week. The second pattern is a daily seasonal pattern (m 2 = 48) with the topology in Figure 1 , which allows the model to distinguish between working days and holidays/weekends, and to make transitions between them. The pattern reflects the tendency to have a higher electricity demand during standard working hours and a lower demand at night. It also reflects the tendency to have different demand patterns on working days and holidays/weekends. A longer series would also have an annual seasonal pattern, but with only 115 days, this cannot be distinguished from the trend component. Figure 4 shows the time series decomposed using STR. The λ coefficients were chosen semiautomatically (the starting point for the minimization procedure was chosen according to experiments involving the minimization of the same problem with fewer predictors). Five-fold crossvalidation with a gap length of g = 336 (one week) for the optimal smoothing parameters yielded the optimal mean squared error. Two seasonal patterns and two regressors are used for the decomposition. Thus, the data are represented as the sum of six components: trend, weekly seasonality, daily seasonality with a complex topology (work-day and non-work-day seasonality including transition periods), temperature and squared temperature (which are time-varying but non-seasonal), and the remainder. It has to be emphasized that the second seasonal component represents work-day, non-work-day seasonal patterns and the patterns of the transition periods between them as a single structure. To the best of our knowledge, no other seasonal decomposition method has this powerful ability. Another unique STR feature, illustrated in this decomposition, is the use of flexible regressors (components four and five). We could also have used a more complex approach with seasonal regressors (which is another unique feature of STR). It is interesting to look carefully at the remainder term, as it will contain evidence of events that affect (or are correlated with) the electricity consumption, but which have not been captured in the model. The ten residuals that are largest in absolute value are shown in Table 2, corresponding to the grey vertical lines in Figure 4 . February, 2000, probably because it was one of the hottest days in Melbourne (40 • C), but with a sudden drop in temperature (see the first five red lines in Figure 5 ). We can speculate that although the temperature dropped at around 2:30 pm, the buildings stayed hot because they had heated up during the very hot previous three days and two nights. Thus, the electricity consumption was higher than that expected by the model, which only accounts for the contemporaneous temperature. A negative outlier appeared at 4:00 pm on the following day, which was also very hot (39.8 • C) (see the separate single red line in Figure 5 ). It is very unusual to have a hot day following a week later than the first four outliers. The explanation is similar to the previous case: it was a very hot day (38.1 • C) after a series of rather cool days and nights, and the model overestimated electricity consumption. The model could be modified to better handle these situations by using lagged as well as concurrent temperatures, along with the change in temperatures over recent periods. Because we have used a linear modelling framework, it is very easy to extend our approach in many directions, exploiting the vast literature on linear regression models. We highlight a few possible extensions here. We can reduce the dimension of the model by using a linear combination of smooth basis functions (such as splines) for the trend: so that = Ψν, where ν contains coefficients to be estimated. Then no difference operators are needed to ensure smoothness. Instead, the value of q and shape of the basis functions ψ j will determine the smoothness. Similarly, seasonal functional components can be obtained, using Fourier terms for example (reminiscent of De Livera et al. 2011) . This approach leads, with no additional effort, to handling seasonality with a fractional or varying period. If we replace the L 2 norm used in (2), with the L 1 norm, we obtain a robust variant of STR. Equivalently, we replace the normality assumption with a double exponential distribution in the expression for the log-likelihood. This can be written as a quantile regression (Koenker 2005) : where y + , X and η are defined as in (4). This can be solved numerically using a quantile regression algorithm (Koenker 2020) . The L 1 and L 2 norms can also be mixed according to different assumptions regarding the distributions of the model components. In such cases, the minimization problem can be reduced to the LASSO minimization problem. A mixture of norms can be useful, for example, in cases where the trend and seasonal patterns are smooth but the noise component has outliers. It can also be useful when the trend changes abruptly, but the seasonal components are smooth, and the noise is distributed normally (or at least has no outliers). A variation of (4) allows the errors to be correlated with covariance matrix Σ ε . Then the maximum likelihood estimates are given by the generalized least squares solution with The STR approach allows us to deal easily with level shifts, spikes, and other intervention patterns as defined by Box and Tiao (1975) . Where these are linear terms, they can simply be added to the X matrix in the usual way. For nonlinear intervention effects such as "shift and decay" patterns, non-linear least squares will need to be used. The model can be used for forecasting by simply treating future observations as missing, and then estimating them. The smooth trend T t will be continued linearly, much like a natural spline. If the remainder term R t is autocorrelated, the forecasts can be adjusted to allow for the shortterm dynamics represented by the autocorrelations (e.g., by fitting an ARMA model to R t ). Prediction intervals can be obtained in the standard way for linear regression models. When the STR model is being used for forecasting, better results will probably be obtained by selecting the smoothing parameters using K-fold cross-validation where the value of g should be comparable to the forecasting horizon h in order to avoid the trend being too flexible. The STR model provides a way to identify new features in time series data, such as weekly or daily patterns. These can then be used within other algorithms. For example, they could be used to classify time series (Lubba et al. 2019) , to find anomalous time series (as distinct from anomalous observations) (Hyndman et al. 2015) , to select a forecasting model (Talagala et al. 2018) , or to create a weighted forecast combination (Montero-Manso et al. 2020). We have looked at two different business applications in monthly supermarket revenue and halfhourly electricity demand. It is easy to imagine many others, and we conclude with just three to illustrate the range of possibilities that could be used with this method. 1. A retail outlet wishes to measure the effect of a promotion by studying the increased number of customers in their stores. The promotion took place during a week that includes a public holiday. They can use STR with seasonal, trend, and holiday effects, along with a categorical variable indicating the periods before, during and after the promotion. This gives a direct measure of the number of additional customers during the promotion, and the number of additional customers in the weeks following the promotion. 3. An energy company needs to forecast gas usage in a region for the next month, and they have ten years of hourly data. However, COVID-19 has caused changes in the usage profile as many people are working from home. STR can be used to estimate the trend, regular seasonal patterns with public holiday effects, and new seasonal patterns since people started working from home (including possibly new public holiday effects). The model can include temperature variables to allow for heating effects, and these can interact with the COVID-19 variables. This article has introduced a new flexible approach to seasonal-trend decomposition using a linear regression model which can be used for a wide range of time series. The STR method allows for multiple seasonal periods and complex seasonality, missing values in input data, provides confidence intervals, finds smoothing parameters, and allows regressors to be taken into account with coefficients that may be time-varying and seasonal. The main novelties of the STR approach are as follows. First, the seasonal components are represented as two dimensional structures of optionally complex topology. This is a completely new way of viewing seasonal components in the field of seasonal-trend decomposition. Second, the influence of covariates is allowed to be time-varying and can be seasonal, cyclic, etc. In all cases, the seasonal, cyclic or other components can have complex constraints which are difficult (if not impossible) to represent using other methods. Third, the approach proposes a regression perspective on time series decomposition, rather than the filtering or stochastic process approach which underpins most other methods. Fourth, beyond the methodological contributions mentioned above, the main theoretical contribution of the new approach is that it allows us to build a bridge between some stochastic processes and linear regression problems, which can give new mathematical interpretations and new computational methods for solving stochastic process problems. Finally, the new approach provides a unified framework for handling a very wide variety of seasonal-trend decomposition problems, which is easy to adapt in order to take into account further time series features. An R package, stR (Dokumentov and Hyndman 2018) , implementing our model is available on CRAN. The package also contains a vignette, which describes various practical aspects of the implementation in detail. This paper was produced using Rmarkdown (Allaire et al. 2020) . The source files to reproduce the paper, including all the examples, are available from github.com/robjhyndman/STR paper. 2020) rmarkdown: Dynamic Documents for R The elimination of spurious correlation due to position in time or space Modeling time-varying trading-day effects in monthly time series An analysis of transformations Intervention analysis with applications to economic and environmental problems STL: A seasonal-trend decomposition procedure based on Loess Statistical software for state space methods Author: STR: Seasonal-Trend decomposition using Regression 22 Article accepted for The X11ARIMA/88 Seasonal Adjustment Method: Foundations and User's Manual (Statistics Canada Seasonal adjustment methods and real time trend-cycle estimation Forecasting time series with complex seasonal patterns using exponential smoothing stR: STR Decomposition Some recent developments and directions in seasonal adjustment New capabilities and methods of the X-12-ARIMA seasonal-adjustment program Stock series holiday regressors generated by flow series holiday regressors Seasonal adjustment and signal extraction in economic time series Forecasting, structural time series models and the Kalman filter The suspension of the Berlin produce exchange and its effect upon corn prices forecast: Forecasting Functions for Time Series and Linear Models Large-scale unusual time series detection Quantile regression 2020) quantreg: Quantile Regression Seasonal Adjustment with the X-11 Method catch22: Canonical time-series characteristics Author: STR: Seasonal-Trend decomposition using Regression Article accepted for FFORMA: Feature-based forecast model averaging 1884) A comparison of the fluctuations in the price of wheat and in the cotton and silk imports into great britain R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing Semiparametric Regression Linear regression analysis The X-11 variant of the Census II method seasonal adjustment program Electronic computers and business indicators On the graduation of the rates of sickness and mortality presented by the experience of the Manchester Unity of Oddfellows during the period 1893-97 Meta-learning how to forecast time series Generalized Additive Models: an introduction with R Rob Hyndman gratefully acknowledges the support of the Australian Centre of Excellence for Mathematical and Statistical Frontiers.