Performance of Machine Learning Models on Time Series Data

In the past few years there has been an increased interest among economists in machine learning. For more discussion see here, here, here, here, here, here, here,  and here.  See also Mindy Mallory's recent post here.

While some folks like Susan Athey are beginning to develop the theory to understand how machine learning can contribute to causal inference, it has carved out a niche in the area of prediction. But what about times series analysis and forecasting?

That is a question taken up by authors this past March in an interesting paper (Statistical and Machine Learning forecasting methods: Concerns and ways forward). They took a good look at the performance of popular machine learning algorithms relative to traditional statistical time series approaches. The authors found that traditional approaches including exponential smoothing and econometric time series approaches out performed algorithmic approaches from machine learning across a number of model specifications, algorithms, and time series data sources.

Below are some interesting excerpts and takeaways from the paper:

When I think of time series methods, I think of things like cointegration, stationarity, autocorrelation, seasonality, auto-regressive conditional heteroskedasticity etc. (I recommend Mindy Mallory's posts on time series here)

Hearing so much about the ability of some machine learning approaches (like deep learning) to mimick feature engineering, I wondered how well algorithmic approaches would handle these issues in time series applications. The authors looked at some of the previous literature in relation to this:

"In contrast to sophisticated time series forecasting methods, where achieving stationarity in both the mean and variance is considered essential, the literature of ML is divided with some studies claiming that ML methods are capable of effectively modelling any type of data pattern and can therefore be applied to the original data [62]. Other studies however, have concluded the opposite, claiming that without appropriate preprocessing, ML methods may become unstable and yield suboptimal results [28]."

One thing about this paper, as I read it, is that it does not take an adversarial or luddite tone toward machine learning methods in favor of more traditional approaches. While they found challenges related to predictive accuracy, they seemed to proactively look deeper to understand why ML algorithms performed the way they did and how to make ML approaches better at time series.

One of the challenges with ML, even with crossvalidation was overfitting and confusion of signals, patterns, and noise in the data:

"An additional concern could be the extent of randomness in the series and the ability of ML models to distinguish the patterns from the noise of the data, avoiding over-fitting....A possible reason for the improved accuracy of the ARIMA models is that their parameterization is done through the minimization of the AIC criterion, which avoids over-fitting by considering both goodness of fit and model complexity."

They also recommend instances where ML methods may offer advantages:

"even though M3 might be representative of the reality when it comes to business applications, the findings may be different if nonlinear components are present, or if the data is being dominated by other factors. In such cases, the highly flexible ML methods could offer significant advantage over statistical ones"

It was interesting that basic exponential smoothing approaches outperformed much more complicated ML methods:

"the only thing exponential smoothing methods do is smoothen the most recent errors exponentially and then extrapolate the latest pattern in order to forecast. Given their ability to learn, ML methods should do better than simple benchmarks, like exponential smoothing."

However the authors note it is often the case that smoothing methods can offer advantages over more complex econometric time series as well (i.e. ARIMA, VAR, GARCH etc.)

Toward the end of the paper the authors go on to discuss in detail the differences in the domains where we have seen a lot of success in machine learning (speech and image recognition, games, self driving cars etc. ) vs. time series and forecasting applications.

In table 10 of the paper, they drill into some of these specific differences and discuss structural instabilities related to time series data, how the 'rules' change and how forecasts themselves can influence future values, and how this kind of noise might be hard for ML algorithms to capture.

This paper is definitely worth going through again and one to keep in mind if you are about to embark on an applied forecasting project.

Reference: 

Makridakis S, Spiliotis E, Assimakopoulos V (2018) Statistical and Machine Learning forecasting methods: Concerns and ways forward. PLoS ONE 13(3): e0194889. https://doi.org/10.1371/journal.pone.0194889

See also Paul Cuckoo's LinkedIn post on this paper: https://www.linkedin.com/pulse/traditional-statistical-methods-often-out-perform-machine-paul-cuckoo/ 

Vector Autoregressions and Bayesian Vector Autoregressions

In a previous post I highlighted the basic algebra associated with ARIMA models:

  

Ex ARMA(2,1,1) : ( Y*_t = b₀ + b₁Y*_t-1  + b₂Y*_t-2  + b₃ e_t-1  + e_t

The terms in this model consist only of a single observed time series and its associated error terms. Vector Autoregressions (VAR's) are class of models that allow for a multivariate analysis of time series. For example 2 series Y₁  and Y₂ (perhaps GNP and treasury yields). We could specify the two series as follows:

Y_1,t  =  c₁ + A_1,1Y_1,t-1 + A_1,2Y_2,t-1 + e_1,t

Y_2,t  =  c₂ + A_2,1Y_1,t-1 + A_2,2Y_2,t-1 + e_2,t    

 

More generally, writing one equation for each series:

Y_t  =  c + A₁Y_t-1 + …+ A_pY_t-p + e

  

Or in matrix form: Y = AZ + U

Y =  series of K different endogenous variables

A=  coefficient matrix

Z =  matrix of lagged variables

U = matrix of error terms

Exogenous variables ‘x’ can be incorporated as follows:

  Y_t  =  c + A₁Y_t-1 + …+ A_pY_t-p + B₁X_t-1 + …+ B_pX_t-p + u

Problems with VAR’s

Given the large number of parameters in VAR’s, overparamaterization and over fitting becomes an issue. Some coefficients may turn out to be different from zero by coincidence, leading to spurious relationships in the model specification.

Bayesian Vector Auto regression (BVAR)

Assume that the model parameters are random variable. A prior distribution is specified based on prior information and combined with objective information from observed data to obtain the posterior distribution using Bayes Theorem.  The prior distribution specification acts as a barrier, or provides ‘shrinkage’ preventing the estimated parameters from depicting what are only spurious correlations.

Estimation

Recall , Bayes Theorem, for parameter θ  and data X,  in its simplest form:

p(θ |X) = p(X|θ)p(θ ) / p(X)  or  p(X|θ)p(θ ) / ∑_i p(X| θ_i) p(θ_i)

The posterior distribution can be defined based on the prior and the likelihood function as follows:

p(θ|X) = p(θ ) L(θ|X) / ∫_θ p(θ) L(θ |X) dθ α p(θ) L(θ |X)

This should look familiar- see also the post on Bayesian Econometrics.

p(θ|X) is then estimated based on the above formulation using MCMC methods. Various prior distributions have been discussed in the literature, including the Minnessota Prior, Conjugate Priors, Flat Priors, and the Sims-Zha Prior.

References:

The Usage of Different Prior Distributions in Bayesian Vector Autoregressive Models. Volkan Sevinc and Gul Ergun. Hacettepe Journal of Mathematics and Statistics Volume 38(1) 2009.

Advances in Bayesian Time Series Modeling and the Study of Politics: Theory, Testing, Forecasting, and Policy Analysis. Patrick T. Brandt and John R. Freeman. Political Analysis. (2006) 14:1-36.

Vector Autoregression. (presentation) Jamie Monogan. Washington University. St. Louis. November 15 & 17, 2010. Link: http://artsci.wustl.edu/~jmonogan/teaching/ts/14var.pdf

ARIMA Model Forecast of Annual Fertilizer Use in KY

I recently found some data on annual fertilizer use (measured by an index) in KY. This looked like a good series to try out  time series analysis in R (see program documentation for data source)as well as provide a basic demo for my earlier post on time series analysis.

The series is fertilizer use in Kentucky from 1960-2004, and is plotted below:

A look at the autocorrelation plot indicates a sinusoidal decay toward zero, which indicates an autoregressive process and partial autocorrelations show a cut off after the first lag, representing a possible AR(1) process(Newbold & Bos, 1990).

Fitting the ARMA(0,0,1) model in R based on values through 2003, yields a forecast of 2.447325 for 2004 vs. the actual observed value of 2.49. This is a much closer estimate than the naive estimate based soley on the previous value for 2003, which was 2.63. The forecast for 2005 is 2.298228, for which I have no data to verify. It would be interesting to build a more complex model that incorporates additional data, perhaps a vector auto-regression.

R code for for this post:

# *------------------------------------------------------------------
# |  forecast for annual fertilizer use in KY (based on index)           
# *-----------------------------------------------------------------
 
# data source: http://www.ers.usda.gov/Data/AgProductivity/
 
# read in year/time variable
yr <- c(1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004)
# read in fertilizer index
fertilizer <- c(1.4279, 1.4114, 1.5069, 1.6599, 1.8183, 1.8229, 2.1135, 1.9576, 1.3736, 1.5381, 1.5041, 1.5090, 1.4749, 1.6329, 1.9789, 1.8462, 2.0444, 1.9997, 1.9384, 1.8392, 2.1652, 1.9528, 1.4552, 1.3719, 1.3323, 1.3283, 1.5702, 1.5339, 1.4294, 1.4010, 1.4534, 1.4372, 1.4505, 1.5637, 1.4135, 1.6958, 1.7138, 1.7090, 2.3354, 2.4032, 2.7095, 1.9193, 2.3007, 2.6395, 2.4940)
 
kydat <- data.frame(yr,fertilizer) # create data frame, necessary for subsetting by yr
names(kydat) # list var names in data frame
dim(kydat)  # rows and columns
kydat2 <- kydat[kydat$yr<2004,] # subset to remove last 5 values for forecasting
dim(kydat2) # rows and columns
ky_ftlzr <- ts(kydat2) # convert to a time series object
plot(ky_ftlzr[,2]) # plot the series objects
acf(ky_ftlzr[,2], 40,xlim=c(1,40)) # autocorrelation plot- sinusoidal decay to zero -indicates autoregressive process 
pacf(ky_ftlzr[,2],40, xlim=c(1,40)) # consistent with a AR(1) model
 
# model 
(x.fit = arima(ky_ftlzr[,2], order = c(1,0,0))) # fit a ARIMA(0,0,1)) AIC = 5084.73
tsdiag(x.fit, gof.lag = 20) # ARIMA model diagnostics
x.fore = predict(x.fit, n.ahead=2) # forecast the next 5 values
print(x.fore) # 2.447325 2.298228
 
# compare these values to the observed value
print(kydat[kydat$yr>=2004,]) # 2004 = 2.49

Created by Pretty R at inside-R.org

Forecasting and Time Series Analysis

Time Series Analysis

Exponential Smoothing:

S = α Y _t +  (1-α) Y_t-1              
S =  portion of current average + portion of previous average = smooth value

Classical Equation:

Forecasted Value = T x S x C                        
FV =trend*seasonal index*cyclical index} 

Regression Equation:

T + S + C    
Y = trend regression with seasonal and cyclical dummies

Seasonal Indices: 

AVG_{season i}  /  AVG_total  = seasonal index for season i.

Econometric Forecasting Models

Distributed Lag Models

Y_t = a₀ + b₀ X _t +  λY_t-1  +  u_t

where λ is on the interval (0,1)  
 

e_t = r e_t-1  +  u_t    (serially correlated error terms)

Serially correlated error terms can be corrected using GLS, iterated maximum likelihood, instrumental variables.

Non-stationary: upward or downward trend in data, this can be corrected by take first differences:

           
         Y* = Y­_t – Y_t-1

Application: To predict  Y_t+1 use the current value of Y.  Predicting Y_t+2 can be predicted as a conditional forecast, given a forecast of Y_t+1 

ARIMA Models

Combines moving average (MA) and autoregressive (AR) models as well as differencing.

ARIMA(p,d,q)  where p = # of autoregressive terms
                                    d = # of 1st differences
                                    q = # moving average terms

Application:  ARIMA(2,1,1) :    Y*_t = b₀ + b₁Y*_t-1  + b₂Y*_t-2  + b₃ e_t-1  + e_t

2 autoregressive terms,  1 1^st difference, 1 moving average term

Seasonal effects can also be incorporated:

ARIMA(p,d,q)(P,D,Q)


where P = # of seasonally differenced autoregressive terms
          D = # of seasonal differences
          Q = # of seasonally differenced moving average terms


This model can also be represented as such:

[θ(B) θ(B^s) / Φ(B) Φ(B^s) Δ^dΔ^s  ] *   a_t 

where  

B = backshift operater such that BY = Yt-1

θ(B) = MA(q)

θ(B^s) = seasonal MA(Q)

Φ(B) = AR(p)

Φ(B^s) = seasonal AR(P)

Δ^d=  differencing

Δ^s = seasonal differencing

a_t = white noise term

References:

Using Econometrics: A Practical Guide (4th Edition)

A.H. Studenmund