Creating autoregressive models in Python

Autoregressive models use values from previous time steps to create a linear regression equation that will predict the value at the next time step. It is a simple but powerful algorithm 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.

You might have heard of autoregressive models like the well known Autoregressive Moving Average (ARMA), Autoregressive Integrated Moving Average ARIMA or Seasonal ARIMA (SARIMA), where they are integrated with a moving average model. Autoregressive models are also used in deep learning text processing algorithms like RNN or LSTM.

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 uses a regression of 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 dependent variable is simply the value from the previous time period. 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). Autoregressive models can contain a single variable (univariate) or have multiple variables (multivariate).

In this tutorial, we'll be using the sktime and statsmodels libraries to analyze time series data set of CO2 readings from the Mauna Loa observatory in Hawaii and to create autoregressives 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 several libraries for creating our Autoregressive models.

First, the sktime library, a Python library for time series analysis and learning tasks such as classification, regression, clustering, annotation, and forecasting. Second, seaborn which is a library for data visualization and the creation of charts. You'll use matplotlib to interact with additional plots. The pandas library is the canonical library for working with loaded data. The statsmodels library will be used to create our ARIMA models and provides utilities for visualization as well.

%pip install sktime seaborn matplotlib pandas statsmodels;

Now that we've installed all of the libraries, we can import them into our environment:

import sktime
from matplotlib import pyplot
import matplotlib as plt
import seaborn
import datetime
import pandas as pd
import statsmodels
import numpy as np

For this tutorial, you'll use the monthly CO2 data from the Mauna Loa Observatory in Hawaii. The data is available here and a direct link is available here. To create a univariate time series, we'll use only the 'average' field from the dataset and construct a date which can be used as the index for the time series.

url = "https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.csv"

# change for your uploaded URL
df = pd.read_csv(url, comment='#')
df['date'] = pd.to_datetime(dict(year=df.year, month=df.month, day=1))
df = df[['date', 'average']]
df = df.set_index(df['date']).drop('date', axis=1)
df.head()

The head command outputs:

date        average
1958-03-01  315.71
1958-04-01  317.45
1958-05-01  317.51
1958-06-01  317.27
1958-07-01  315.87

The statsmodels library requires that our dataframe have an explicit temporality, so resample the data with 'MS' to ensure that the dataframe knows that the data is using a Month Start interval.

df = df.resample("MS").last()

Step 3. Load the data and perform EDA

Before we start estimating what kind of ARIMA model, we should do some quick exploratory data analysis.

First, we can look at the statistics of our dataset:

df.describe().T

There's quite a wide standard deviation, which might indicate some non-stationarity in the time series. Plotting the time series can also help show trend and seasonality:

pyplot.plot(df)
pyplot.xticks(rotation=45)
pyplot.show()

Step 4. Check stationarity and variance

Now it's time to establish whether our time series model is stationary or not. Autoregressive models generally aren't able to model non-stationary data.

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.

from statsmodels.tsa.stattools import adfuller
adf_result = adfuller(df)
print(adf_result[1])

This will output:

    1.0

from statsmodels.tsa.stattools import kpss
kpss_result = kpss(df)
print(kpss_result[1])

This will output:

    0.01

Since the ADF test value is far above 0.05 and the KPSS value is far below 0.05, this indicates that our time series is not stationary. We can check trend and seasonality using the seasonal_decompose method:

from statsmodels.tsa.seasonal import seasonal_decompose
result = seasonal_decompose(df, model='additive')
result.plot()
pyplot.show()

To correct this, we can difference the time series by subtracting xt-1 from xt, which will give us a time series which consists of the changes between values in the time series. This is an easy way to remove a trend from a time series and create a more stationary data set.

The sktime library has a Differencer class which allows you to difference both using lagged values and seasonal components. In this case, we'll just use a lag of 1 which will subtract x_t-1 from x_t.

from sktime.transformations.series.difference import Differencer
transformer = Differencer(lags=1)
df['average'] = transformer.fit_transform(df['average'])

Now we can inspect the differenced data using the seasonal decomposition method.

from statsmodels.tsa.seasonal import seasonal_decompose
result = seasonal_decompose(df, model='additive')
result.plot()
pyplot.show()

This looks much better. We can check the ADF test results on our new time series:

adf_result = adfuller(df)
print('{:.10f}'.format(adf_result[1]))

This will output:

    0.0000056113

Since this value is below 0.05 we can conclude that we now have a stationary time series.

Before we begin plotting autocorrelations and training a model, we'll create a test and train split for our data set so that we can test the accuracy of our trained model on validation data.

from sktime.split import temporal_train_test_split
train, test = temporal_train_test_split(df, test_size=36)

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 auto-correlation. Calculating the autocorrelation can answer questions about whether the data exhibit 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 and the Partial Autocorrelation Function.

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 python we can use the autocorrelation_plot function to create an auto-correlation plot.

from pandas.plotting import autocorrelation_plot
#just use the first 48 to make the plot easier to read
autocorrelation_plot(train[-48:])

We're looking for values that are outside of the critical range of 0.3 to -0.3. Across 48 time units, the time series shows a significant correlation at the 12th and 24th lag and a weakly negative correlation at lags 6 through 9. This would indicate that our time series still has a seasonality component to it and that a model of an order of 12 or 13 may fit our data well. The autocorrelation drifts to zero across the lags which is a good sign that our data may be compatible with Autoregressive modeling.

It is also a good idea to check the Partial Autocorrelation Function (PACF) as well. An ACF plot shows the the relationship between yt and yt−k for different values of k. If y_t and y_t−1 are correlated, then y_t−1 and y_t−2 will also be correlated. But it's also possible for y_t and y_t−2 to be correlated because they are both connected to y_t−1, rather than because of any new information contained in y_t−2 that could be used in forecasting y_t. To overcome this problem, we can use partial autocorrelations to remove a number of lag observations. These measure the relationship between y_t and y_t−k after removing the effects of lags 1 to k. So the first partial autocorrelation is identical to the first autocorrelation, because there is nothing between them to remove. This can often give a more precise estimate of which lags might contain indications of seasonality and unit roots.

from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(train)
pyplot.show()

Again, we can see positive values at lags 1 and 12 and negative values at lags 6 through 9 because of the seasonality of our data set. After the 12th lag the values decrease and are mostly within the critical range so this indicates that an autoregressive model may fit our data well.

A rough rubric for when to use Autoregressive models 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

We can see all three of these in our time series data, so it will be a good candidate for modeling with an autoregressive model.

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, Yt on Yt-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 Python 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.

First, we'll specify an autoregressive model using the lags at 12 and 24 to account for the seasonality. We create the model by passing the training data and then specifying all of the lags that we would like to include. This can take two forms: if an integer is passed then it's interpreted as the number of lags to include in the model, if it's a list then this is interpreted as the lag indices to include. In this case, we pass an array with 12 and 24 to indicate that we want to only use those lags.

from statsmodels.tsa.ar_model import AutoReg
params = AutoReg(train, lags = [12, 24])
res = params.fit()
res.summary()

The summary of the autoregressive model indicates the fitted parameters for the autoregressive model. If you've worked with linear regressions in the past then these statistics should be familiar to you and their interpretation is much the same. The most important statistics are the coefficients for the lag and the P-value of that lag. The lags with the strongest coefficient and the lowest P-value are likely to be the best fitting order for your time series data. The confidence intervals for the coefficients can also be indicative of the best fitting lags.

For the sake of comparison we can utilize a convenience function built into the statsmodel library called ar_select_order which will select an number of lags to optimize for the metric selected. The method signature has the following parameters:

• endog: this is the endogenous data
• maxlag: the maximum number of lags that can be used
• ic: the information criterion used to select the number of lags

• trend: The trend to include in the model. It has

  • ‘n’ - No trend.
  • ‘c’ - Constant only.
  • ‘t’ - Time trend only.
  • ‘ct’ - Constant and time trend.

• seasonal: Flag indicating whether to include seasonal dummies in the model

• period: The period of the data. Only used if seasonal is True.

In this case, we can specify AIC as the metric that we would like to use. We can also specify that we believe the data to be seasonal and that a constant trend and a time trend may be present.

from statsmodels.tsa.ar_model import ar_select_order
select = ar_select_order(train, 36, "aic", trend="ct", seasonal=True)
print(select.ar_lags) # show all selected lags

We can see that the ar_select_order has taken 2 years worth of monthly values:

    [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]

Let's check the model fit:

res = select.model.fit()
print("AIC: {0} \nBIC: {1} \nHQIC: {2}".format(res.aic, res.bic, res.hqic))

This will output:

    AIC: 364.2537584764123
    BIC: 543.5424408499916
    HQIC: 433.41145119060553

We can see that the AIC is roughly half of our previous model. The BIC and HQIC are also significantly lower as well.

When we check the estimated parameters and the errors associated with them, not all of the parameters seem highly significant. This could be a sign that we might want to recreate our model using just the significant parameters and removing the nuisance parameters but since the model is still fairly small and the time series not too long we can leave it as is:

res_stats = pd.DataFrame({'parameter estimations':res.params, 'errors':res.bse})
print(res_stats)

This will show us the parameter estimations for all of our lags:

                 parameter estimations    errors
    const                     0.487935  0.195844
    trend                     0.000495  0.000089
    s(2,12)                   0.549128  0.138676
    s(3,12)                   0.415967  0.237492
    s(4,12)                  -0.236964  0.301942
    s(5,12)                  -0.929309  0.332358
    s(6,12)                  -1.536317  0.338571
    s(7,12)                  -1.657907  0.339382
    s(8,12)                  -0.821790  0.343709
    s(9,12)                   0.100541  0.338144
    s(10,12)                  0.122370  0.306786
    s(11,12)                  0.002686  0.239431
    s(12,12)                 -0.149991  0.139858
    average.L1               -0.336318  0.036770
    average.L2               -0.180672  0.038767
    average.L3               -0.170620  0.039369
    average.L4               -0.133652  0.039802
    average.L5               -0.082281  0.039929
    average.L6               -0.066625  0.039841
    average.L7               -0.069048  0.039491
    average.L8               -0.041535  0.039413
    average.L9               -0.022176  0.039236
    average.L10              -0.065015  0.039154
    average.L11               0.037157  0.039230
    average.L12               0.102359  0.039222
    average.L13               0.039471  0.039316
    average.L14               0.025184  0.039187
    average.L15              -0.012151  0.039183
    average.L16              -0.067055  0.039187
    average.L17              -0.109872  0.039233
    average.L18              -0.081415  0.039247
    average.L19              -0.140589  0.039251
    average.L20              -0.102440  0.039321
    average.L21              -0.098774  0.039345
    average.L22              -0.049864  0.039241
    average.L23              -0.005906  0.038709
    average.L24               0.047031  0.038281
    average.L25               0.092464  0.036224

You can see that not all of these are highly significant, so you could remove some of the parameters to make the model more parsimonious.

Now we check the diagnostic plots of our model. We're looking for more or less steady residuals across the time series, an error density centered around 0, and a correlogram which doesn't show evidence of strong seasonality in the residuals.

fig = pyplot.figure(figsize=(16, 9))
fig = res.plot_diagnostics(fig=fig, lags=24)

We see steady residuals across the time series, an error density centered around 0, and a correlogram which doesn't show evidence of strong seasonality, so this is an indication that our model fits well.

Step 7. Forecast using the autoregressive model

Now that you've fit an AR model you can use it for time series forecasting by predicting values and then comparing those forecasted future values to the held-out test data set which we created in Step 4.

We'll update the model with the previous values each time we forecast the current value:

testing = pd.concat([train,test])
ar24_predictions = []

for i in range(len(train), len(train) + len(test)):
    updated_res = res.apply(testing[:i])
    ar24_predictions.append(updated_res.forecast(1).iloc[0])

Now we can plot our predictions and the actual values for the data set.

predictions = pd.DataFrame({
    "AR": ar24_predictions,
    "Actual": test['average'],
})
_, ax = pyplot.subplots()
ax = predictions.plot(ax=ax)

We can see that our regression model did quite well. To test the efficacy of our model we can even forecast without providing updated values to the model as it forecasts the test set:

predictions = pd.DataFrame({
    "AR": res.forecast(36),
    "Actual": test['average'],
})
_, ax = pyplot.subplots()
ax = predictions.plot(ax=ax)

You can check the Root Mean Squared Error for the predicted values and the test using the mean_squared_error function from the sklearn.metrics library:

from sklearn.metrics import mean_squared_error
rmse = np.sqrt(mean_squared_error(predictions['AR'], test['average']))
print('Test RMSE: %.3f' % rmse)

This will output:

Test RMSE: 0.462

Your autoregressive model is able to model the original data points with good accuracy.

Summary and next steps

In this tutorial, you learned how autoregressive models can be used to predict time series data. We explored the Mauna Loa CO2 dataset and transformed the time series using differencing so that an autoregressive model could be more easily fit. You then developed two different Autoregressive models using the ar and ar_select_order methods and compared the two of them.

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.