Creating autoregressive models in R

Autoregression is the creation of an model that uses values from previous time steps to create a regression equation that can predict the value at the next time step. It is a simple, but very powerful, technique for time series analysis that can give you highly interpretable and very effective predictions if your time series data contains data that is strongly correlated across the time steps.

A time series is a sequence of measurements of the same variable or group of variables made over time. Usually the measurements are made at evenly spaced times, hourly, monthly, or yearly. As an example, we might have values that measure the number of airline passengers in a country, with measurements observed each month. In this case, the passenger counts would be the y that we are measuring and to emphasize that we have measured values over time, we use t as a subscript rather than the usual i, so that y_t means the value of y at time period t.

An autoregressive model is when we regress a value from a time series on previous values from that same time series. For example, y_t regressed on y_t-1. In this simple regression model, the response variable in the previous time period has become the predictor. The errors represent all the usual assumptions about errors in a simple linear regression model. We often look at the order of an autoregression is the number of preceding values in the series that are used to predict the value at the present time. So, y_t regressed on y_t-1 is a first-order autoregression, which is usually written as AR(1).

In this tutorial, we'll be using the forecast library both to load a sample time series data set of Australian Gas production and to create the autoregressive model based on that data set. Let's get started.

Prerequisites

You need an IBM Cloud account to create a watsonx.ai project.

Steps

Step 1. Set up your environment

While you can choose from several tools, this tutorial walks you through how to set up an IBM account to use a Jupyter Notebook. Jupyter Notebooks are widely used within data science to combine code, text, images, and data visualizations to formulate a well-formed analysis.

1. Log in to watsonx.ai using your IBM Cloud account.

2. Create a watsonx.ai project.

3. Create a Jupyter Notebook using the R Kernel.

This step will open a notebook environment where we can import a data set and create an autoregressive model.

Step 2. Install and import libraries

We’ll use two libraries for creating our autoregressive models (AR models). First, the forecast package, which is a library that contains methods and tools for displaying and analyzing univariate time series forecasts including exponential smoothing via state space models and automatic ARIMA modeling. Second, the tseries library which is a library containing utilities for series analysis and computational finance.

library(forecast)
library(tseries)

Step 3. Load the data and do EDA

Now, we'll load a library from the forecast package:

df <- forecast::gas

AR models assume that the time series is stationary, which means that its statistical properties do not change over time. In practice, many real-world time series are non-stationary, which require preprocessing steps like differencing.

You can check the descriptive statistics of the time series with the summary method:

summary(df)

This will output:

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   1646    2675   16788   21415   38628   66600

It's also helpful to plot a time series to see visible trends and variance:

autoplot(df)

We can see an obvious trend and a change in variance throughout the time series as well.

Step 4. Check the stationarity and variance

Now, we need to establish whether our time series model is stationary or not. ARIMA models can handle cases where the non-stationarity is due to a unit-root but might not work well at all when non-stationarity is of another form.

There are three basic criteria for a series to be classified as stationary series:

• The mean of the series should not be a function of time, which means that the mean should be roughly constant though some variance can be modeled.
• The variance of the series should not be a function of time. This property is known as homoscedasticity.
• The covariance of the terms, for instance the k^th term and the (k + n)^th term, should not be a function of time.

Essentially a stationary time series is a flat looking series, without trend, constant variance over time, a constant autocorrelation structure over time, and no periodic fluctuations. It should resemble white noise as closely as possible. It’s fairly evident that our time series is non-stationary so our model will likely need to reflect that.

We can confirm that our time series is non-stationary by using two tests: the Augmented Dickey-Fuller test (ADF) and the Kwiatkowski-Phillips-Schmidt-Shin test (KPSS).

First, the ADF test. This is a type of statistical test called a unit root test. The ADF test is conducted with the following assumptions:

• The null hypothesis for the test is that the series is non-stationary, or series has a unit root.
• The alternative hypothesis for the test is that the series is stationary, or series has no unit root.

If the null hypothesis cannot be rejected then this test might provide evidence that the series is non-stationary. This means that when you run the ADF test, if the Test statistic is less than the Critical Value and p-value is less than 0.05, then we reject the null hypothesis, which means the time series is stationary.

adf.test(df)

This will output:

    Augmented Dickey-Fuller Test

data:  df
Dickey-Fuller = -2.7131, Lag order = 7, p-value = 0.2764
alternative hypothesis: stationary

The ADF test here returns a p-value of 0.2764 which indicates that the time series is not stationary and has a unit root.

We also can run the KPSS test. The KPSS test is also a unit root test and it figures out if a time series is stationary around a mean or linear trend or is non-stationary due to a unit root. The KPSS results are interpreted in the opposite way from the ADF test:

• The null hypothesis for the test is that the data is stationary.
• The alternate hypothesis for the test is that the data is not stationary.

So unlike the ADF test, if the null hypothesis cannot be rejected then this test provides evidence that the series is stationary. This means that when you run the ADF test, if the Test statistic is less than the Critical Value and p-value is less than 0.05, then we reject the Null Hypothesis(HO); that is, the time series is not stationary. So, ADF with a p-value of less than 0.05 indicates stationarity and the KPSS with a p-value of more than 0.05 indicates stationarity.

kpss.test(df)

This will output:

    KPSS Test for Level Stationarity

data:  df
KPSS Level = 7.7018, Truncation lag parameter = 5, p-value = 0.01

Here, the KPSS test also indicates that our time series is not stationary.

Having two tests which may or may not agree can be confusing. If KPSS and ADF agree that the series is stationary then we can consider it stationary and there’s no need to to use a difference term in the ARIMA model.

If the ADF test finds a unit root but KPSS finds that the series is stationary around a deterministic trend then the series is trend-stationary, and it needs to be detrended. You can either difference the time series by using the diff method or a Box-Cox transformation may remove the trend.

If the ADF does not find a unit root but KPSS claims that it is non-stationary then the series is difference stationary. We'll need to difference the data before passing data to our autoregressive model.

If KPSS and ADF agree that the series is non-stationary then we can consider it non-stationary. We'll need to difference the data before passing data to our autoregressive model.

We'll use the decompose method to determine what may need to be removed from the time series in order to fix it.

decomp <- decompose(df)
autoplot(decomp)

This plot shows that our data has a strong upward trend as well as a constant seasonal copmonent as well. In order to remove this trend, we can difference the data. In differencing the value of a time series at time t is substracted from the value at time t-1 so that the time series becomes simply the changes between the readings in the time series.

The diff method generates a differenced time series that we can use to train an autoregressive model:

diff <- diff(df)
autoplot(diff)

This time series still shows high variance, or heteroskedasticity, so it will still be difficult to fit an autoregressive model to it. When we encounter heteroskedasticity in a time series it's common to take a log transformation to manage it. One tool for doing this is the Box-Cox transformation. The BoxCox method takes a lambda parameter which is defined as follows:

This lambda value is passed to the BoxCox method to remove the excess variance from the time series through a log transform:

lambda <- BoxCox.lambda( df )
transform <- BoxCox(df, lambda)

Now, we can plot our time series with variance corrected:

plot(transform)

Since we've created a differenced and transformed data set in order to correct for trend and heteroskedasticity, to get values from our predictions, we'll need to inverse the transform using the InvBoxCox method and then sum all of the predictions. We'll cover this in the final section of this tutorial.

Step 5. Plot autocorrelations

Now, we need to check to see where our de-trended time series contains autocorrelations. We can see the degree to which a time series is correlated with its past values by calculating the autocorrelation. Calculating the autocorrelation can answer questions about whether the data exhibits randomness and how related one observation is to an immediately adjacent observation. This can give us a sense of what sort of model might best represent the data.

The autocorrelations are often plotted to see the correlation between the points up to and including the lag unit. There are two approaches to autocorrelations: the autocorrelation function (ACF) and the partial autocorrelation function (PACF).

In ACF, the correlation coefficient is in the x-axis whereas the number of lags (referred to as the lag order) is shown in the y-axis. In R, we can use the acf function to create an autocorrelation plot.

acf(transform)

The de-trended series still shows a strongly positive correlations throughout the series. The rate of decay in an ACF plot can signify the persistence of relationships. A rapid decline indicates a short-lived correlation, common in volatile or random time series data. In contrast, a slow decay suggests a long-lasting relationship which is more typical of cyclical patterns.

It is also a good idea to check the PACF as well. The PACF is helpful in determining the order of the AR part of the ARIMA model. It is also useful to determine or validate how many seasonal lags to include in the forecasting equation of a moving average based forecast model for a seasonal time series.

pacf(transform)

We can see that the values of the PACF plot drop off quickly after the 1.0 lag, indicating that lags after 1.0 are less significant.

A rough rubric for when to use AR terms in the model is when:

• ACF plots show autocorrelation decaying towards zero
• PACF plot cuts off quickly towards zero
• ACF of a stationary series shows positive at Lag 1

A rough rubric for when to use MA terms in the model is when:

• The ACF is Negatively Autocorrelated at Lag 1
• ACF that drops sharply after a few lags
• PACF decreases gradually rather than suddenly

We can see from our two plots that our ACF plot shows a positive value at Lag 1 and that the PACF cuts off somewhat quickly. This might indicate that we should use an AR term in our model.

Now we'll split our data set into a training set and a held out 36 values to use as a test set to compare our models performance:

train <- head(transform, length(transform) - 36)
h <- length(transform) - length(train)
test <- tail(transform, h)

Step 6. Training the model

An autoregressive model is when a value from a time series is regressed on previous values from that same time series. For example, y_t on y_t-1 the linear model would be:

The order of an autoregression is the number of immediately preceding values in the series that are used to predict the value at the present time. So, the preceding model is a first-order autoregression, written as AR(1). An AR(2) model might look like the following:

Here, we use the 2 preceding values in the regression to calculate the value of the series at time t. When we create a model in R, we either specify an order for the autoregression or we can allow methods to select the most appropriate order for the time series data. The order selection can be done by various methods, one of the most common is to use the Aikike Information Criteria or AIC of each possible order and compare them to find the lowest value for the order to help in model selection.

To create an autoregressive model in R, we pass our differenced train data set to the AR method to create the fitted model. There are different ways to call the method:

• x: a univariate or multivariate time series.
• aic: If TRUE then the Akaike Information Criterion is used to choose the order of the autoregressive model.

• method: character string specifying the method to fit the model. Must be one of the strings in the default argument (the first few characters are sufficient). Available methods include:

  • Ordinary Least Squares indicated by 'ols'
  • Maximum Likelihood Estimation, indicated by 'mle'
  • Yule-Walker, indicated by 'yw'
  • Burg method, indicated by 'burg'.

This is optional as the Yule-Walker method is chosen by default. These can also be called by appending the method to ar, like ar.mle or ar.yw.

In this case, we'll use the OLS method and use AIC to select the order of the model:

fit <- ar.ols(train, aic = TRUE)

We can check the order of the model by checking the fit$order property. In this case, an order of 26 has been selected which reflects the monthly seasonality of the time series.

fit$order

This will output:

[1] 26

To check the residuals of our model we can graph the fit$resid property:

autoplot(fit$resid)

The residuals increase as the variance of the data increases since autoregressive models can't account for heteroskedasticity as well as other model types. Now that we've fitted a model we can use it to predict our test dataset and plot the predictions:

prediction <- predict(fit, n.ahead = 36)
autoplot(prediction$pred, series="Prediction") + autolayer(test, series="Actual")

Our model overpredicts the troughs in the time series but does mostly approximate the target time series.

Another way to create an autoregression model is to simply use the ARIMA method and specify only an AR component to the ARIMA model. Instead of an AR model with an order of 26, we'll create a model with an order of 13 for comparison. The ARIMA model allows for both autoregressive and moving average components but to have a better comparison with the AR model, we pass order=c(13, 1, 0) to the ARIMA method. This indicates an autoregressive order of 13, a difference of 1, and no moving average component. Using the ARIMA method is identical to the AR method but additionally returns both log-likelihood and sigma² metrics for additional model diagnostics, which can make it easier to work with.

fit_arima <- arima(train, order=c(13, 1, 0), include.mean = TRUE)
prediction_arima <- predict(fit_arima, n.ahead = 36)
autoplot(prediction_arima$pred, series="Prediction") + autolayer(test, series="Actual")

The differencing order will remove one from the number of observations that our second model will observe.

Step 7. Diagnostics

We can see that the values look fairly similar for the predictions. Since both models use standard errors we comparing the residuals from the two models may help us understand the differences between them.

compare <- ts.union(fit_arima$residuals, fit$resid)
autoplot(compare, facets = TRUE)

We can see that the residuals in each case mirror one another closely. To see how the two different models and their compare different accuracy metrics of the two predictions. First the generated ARIMA model:

forecast::accuracy(prediction_arima$pred, test)

                 ME     RMSE      MAE        MPE      MAPE      ACF1 Theil's U
Test set -0.1460751 0.215147 0.164628 -0.8541112 0.9582771 0.2391624 0.7748088

Second, the hand-selected model:

forecast::accuracy(prediction$pred, test)

                 ME      RMSE       MAE        MPE      MAPE      ACF1 Theil's U
Test set -0.1159915 0.1959326 0.1391823 -0.6770549 0.8091223 0.1263917 0.6976378

We can see here that the AR(13) with differencing model performs better than the AR(26) model without. This is common with data that exhibit a strong trend even after transformation like the data set that we're working with in this tutorial. One thing you'll note is that the model is predicting the transformed values rather than the original values of the data set. This is because we used the Box Cox Transformation on the time series before training our AR models. While this leads to better predictive power in our models, it can be difficult to map our predictions back to the original values of the data set. To convert the predictions to the scale of the original values we'll reverse the Box-Cox transformation:

inv <- InvBoxCox(arima$pred, BoxCox.lambda( df ))
autoplot(tail(df, 36), series="Actual") + autolayer(inv, series="Prediction")

This shows us the untransformed predictions from the AR(13) model.

Summary and next steps

In this tutorial, you learned how autoregressive models can be used to predict time series data. We look at the Australian Gas Monthly data set and transformed the time series using a Box-Cox transform so that an autoregressive model could be more easily fit. You then developed two different autoregressive models using the AR and ARIMA methods and compared the two of them. After including differencing to create a better fit, we inversed the Box-Cox transform to get untransformed values of the prediction.

Try watsonx for free

Build an AI strategy for your business on one collaborative AI and data platform called IBM watsonx, which brings together new generative AI capabilities, powered by foundation models, and traditional machine learning into a powerful platform spanning the AI lifecycle. With watsonx.ai, you can train, validate, tune, and deploy models with ease and build AI applications in a fraction of the time with a fraction of the data.

Try watsonx.ai, the next-generation studio for AI builders.

Next steps

Explore more articles and tutorials about watsonx on IBM Developer.