Auto Regressive Models
In an Autoregressive model, or AR model, today’s value equals a mean plus a fraction phi of yesterday’s value, plus noise. Since there is only one lagged value on the right hand side, this is called an AR model of order 1, or simply an
AR(1) model. If the AR parameter, phi, is one, then the process is a random walk. If phi is zero, then the process is white noise. In order for the process to be stable and stationary, phi has to be between -1 and +1.
As an example, suppose
R(t) is a time series of stock returns. If
phi is negative, then a positive return last period, at time
t-1, implies that this period’s return is more likely to be negative. This is referred to as “mean reversion”. On the other hand, if
phi is positive, then a positive return last period implies that this period’s return is expected to be positive. This is referred to as “momentum”.
AR(1) model parameters
Here are four simulated time series with different AR parameters. When phi equals
0.9, it looks close to a random walk. When phi equals
-0.9, the process looks more erratic - a large positive value is usually followed by a large negative one. The bottom two are similar, but are less exaggerated and closer to white noise.
Here are four autocorrelation functions for different AR parameters. The autocorrelation decays exponentially at a rate of phi. Therefore if phi is
lag-1 autocorrelation is
lag-2 autocorrelation is
lag-3 autocorrelation is
(0.9)^3, etc. When phi is negative, the autocorrelation function still decays exponentially, but the signs of the autocorrelation function reverse at each lag.
Higher order AR models
So far, we’ve been only looking at AR(1) models. The model can be extended to include more lagged values and more phi parameters. Here we show an AR(1), an AR(2), and an AR(3).
Simulating an AR process
Often, if you want to study and understand a pure AR process, it is useful to work with simulated data.
statsmodels provides modules for simulating AR processes. First, import the class,
ArmaProcess. Then define the
parameters of the AR process.
Simulate AR(1) Time Series using
We will simulate and plot a few
AR(1) time series, each with a different parameter,
phi, using the
arima_process module in
statsmodels. First we will look at an AR(1) model with a large positive and a large negative
There are a few conventions when using the
arima_process module that require some explanation. First, these routines were made very generally to handle both AR and MA models.
Second, when inputting the coefficients we must include the zero-lag coefficient of 1, and the sign of the other coefficient is the opposite of what we have been using. For example, for an
AR(1) process with
phi = 0.9, the second element of the
ar array should be the opposite sign,
-0.9. This is consistent with the time series literature in the field of signal processing. We also have to input the MA parameters. We will learn about MA models later, so for now, just ignore the MA part. Then, we create an instance of the class
ArmaProcess. To simulate data, use the method
generate_sample, with the number of simulated samples as an argument.
So, for example, for an
AR(1) process with
phi = 0.9, the array representing the AR parameters would be
ar = np.array([1, -0.9]).
import matplotlib.pyplot as plt import numpy as np # import the module for simulating data from statsmodels.tsa.arima_process import ArmaProcess plt.figure(figsize=(10, 10)) # Plot 1: AR parameter = +0.9 plt.subplot(2,1,1) ar1 = np.array([1, -0.9]) ma1 = np.array() AR_object1 = ArmaProcess(ar1, ma1) simulated_data_1 = AR_object1.generate_sample(nsample=1000) plt.plot(simulated_data_1) # Plot 2: AR parameter = -0.9 plt.subplot(2,1,2) ar2 = np.array([1, 0.9]) ma2 = np.array() AR_object2 = ArmaProcess(ar2, ma2) simulated_data_2 = AR_object2.generate_sample(nsample=1000) plt.plot(simulated_data_2) plt.show()
The two AR parameters produce very different looking time series plots, but next we’ll see how to really be able to distinguish the time series.
Compare the ACF for Several AR Time Series
The autocorrelation function decays exponentially for an AR time series at a rate of the AR parameter. For example, if the AR parameter,
phi = 0.9 , the first-lag autocorrelation will be
0.9, the second-lag will be
(0.9)^2, the third-lag will be
(0.9)^3, etc. A smaller AR parameter will have a steeper decay. For a negative AR parameter, say
-0.9, the decay will flip signs, so the first-lag autocorrelation will be
-0.9, the second-lag will be positive and third will be negative, etc.
Let’s compute the autocorrelation function for each of the three simulated datasets using the
plot_acf function with 20 lags (and suppress the confidence intervals by setting alpha=1).
# Import the plot_acf module from statsmodels from statsmodels.graphics.tsaplots import plot_acf # Plot 1: AR parameter = +0.9 plot_acf(simulated_data_1, alpha=1, lags=20) plt.show() # Plot 2: AR parameter = -0.9 plot_acf(simulated_data_2, alpha=1, lags=20) plt.show() ar3 = np.array([1, -0.3]) ma3 = np.array() AR_object3 = ArmaProcess(ar3, ma3) simulated_data_3 = AR_object3.generate_sample(nsample=1000) # Plot 3: AR parameter = +0.3 plot_acf(simulated_data_3, alpha=1, lags=20) plt.show()
The ACF plots match what we predicted.
Estimating and forecasting an AR Model
statsmodels has another module for estimating the parameters of a given AR model.
import ARMA, which is a class, and create an instance of that class called mod, with the arguments being the data that you’re trying to fit, and the order of the model. The order (1,0) means you’re fitting the data to an AR(1) model. An order (2,0) would mean you’re fitting the data to an AR(2) model. The second part of the order is the MA part, which will be discussed in the next chapter. Once you instantiate the class, you can use the method fit to estimate the model, and store the results in result.
To see the full output, use the summary method on result. The coefficients for the mean
phi are highlighted in red. In the simulated data,
mu was zero and
phi = 0.9, and you can see that the estimated parameters are very close to the true parameters.
If you just want to see the coefficients rather than the entire regression output, you can use the
.params property, which returns an array of the fitted coefficients,
phi in this case.
Example: Estimate the AR parameter using simulated data
In this example we will estimate the
phi, of one of the simulated series that we generated earlier. Since the parameters are known for a simulated series, it is a good way to understand the estimation routines before applying it to real data.
simulated_data_1 with a true
phi of 0.9, we will print out the estimate of
phi. In addition, we will also print out the entire output that is produced when we fit a time series, so we can get an idea of what other tests and summary statistics are available in statsmodels.
Create an instance of the
modusing the simulated data
order (p,q)of the model (in this case, for an AR(1)), is order=(1,0).
import warnings warnings.filterwarnings('ignore') # Import the ARMA module from statsmodels from statsmodels.tsa.arima_model import ARMA # Fit an AR(1) model to the first simulated data mod = ARMA(simulated_data_1, order=(1, 0)) res = mod.fit() # Print out summary information on the fit print(res.summary()) # Print out the estimate for the constant and for phi print("When the true phi=0.9,true mu=" + str(np.round(np.mean(simulated_data_1), 2)) + " the estimate of phi and mu are:") print(res.params)
ARMA Model Results ============================================================================== Dep. Variable: y No. Observations: 1000 Model: ARMA(1, 0) Log Likelihood -1429.721 Method: css-mle S.D. of innovations 1.010 Date: Sun, 28 Feb 2021 AIC 2865.442 Time: 18:42:08 BIC 2880.166 Sample: 0 HQIC 2871.038 ============================================================================== coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ const 0.2416 0.321 0.753 0.451 -0.387 0.870 ar.L1.y 0.9013 0.014 66.041 0.000 0.875 0.928 Roots ============================================================================= Real Imaginary Modulus Frequency ----------------------------------------------------------------------------- AR.1 1.1095 +0.0000j 1.1095 0.0000 ----------------------------------------------------------------------------- When the true phi=0.9,true mu=0.25 the estimate of phi and mu are: [0.24164318 0.90133825]
Notice how close the estimated parameter is to the true parameter.
Forecasting with an AR Model
In addition to estimating the parameters of a model that we did in the last section, we can also do forecasting, both in-sample and out-of-sample using
in-sampleis a forecast of the next data point using the data up to that point, and
out-of-sampleforecasts any number of data points in the future. These forecasts can be made using either the
predict()method if you want the forecasts in the form of a series of data, or using the
plot_predict()method if we want a plot of the forecasted data. We supply the starting point for forecasting and the ending point, which can be any number of data points after the data set ends.
For the simulated series
phi=0.9, we will plot
- Create an instance of the ARMA class called mod using the simulated data simulated_data_1 and the order (p,q) of the model (in this case, for an AR(1) order=(1,0)
- Fit the model mod using the method .fit() and save it in a results object called res
- Plot the in-sample and out-of-sample forecasts of the data using the plot_predict() method
- Start the forecast 10 data points before the end of the 1000 point series at 990, and end the forecast 10 data points after the end of the series at point 1010
# Import the ARMA module from statsmodels from statsmodels.tsa.arima_model import ARMA fig, ax = plt.subplots(figsize=(10, 5)) # Forecast the first AR(1) model mod = ARMA(simulated_data_1, order=(1, 0)) res = mod.fit() res.plot_predict(start=990, end=1010, ax=ax) plt.show()
Notice how, when phi is high like here, the forecast gradually moves to the long term mean of zero, but if phi were low, it would move much quicker to the long term mean. As you can see below for
phi = 0.3
fig, ax = plt.subplots(figsize=(10, 5)) # Forecast the first AR(1) model mod = ARMA(simulated_data_3, order=(1, 0)) res = mod.fit() res.plot_predict(start=990, end=1010, ax=ax) plt.show()
Example: Forecast Interest Rates
We will now use the forecasting techniques and apply it to real data rather than simulated data. The annual data of 10-year interest rates going back 56 years, which is in a Series called
annual_interest_rate_data. Being able to forecast interest rates is of enormous importance, not only for bond investors but also for individuals like new homeowners who must decide between fixed and floating rate mortgages.
In an earlier post, we saw that there is some mean reversion in interest rates over long horizons. In other words, when interest rates are high, they tend to drop and when they are low, they tend to rise over time. Currently they are below long-term rates, so they are expected to rise, but an AR model attempts to quantify how much they are expected to rise.
import pandas as pd daily_rates = pd.read_csv('data/daily_interest_rates_us.csv', index_col=['DATE'], parse_dates=True) daily_rates.columns = ['APR'] daily_rates.head()
annual_interest_rate_data = daily_rates.resample('A').mean() annual_interest_rate_data.head()
# Import the ARMA module from statsmodels from statsmodels.tsa.arima_model import ARMA fig, ax = plt.subplots(figsize=(10, 5)) # Forecast interest rates using an AR(1) model mod = ARMA(annual_interest_rate_data, order=(1, 0)) res = mod.fit() # Plot the original series and the forecasted series res.plot_predict(start='2000', end='2026', ax=ax) plt.legend(fontsize=8) plt.show()
According to an AR(1) model, 10-year interest rates are forecasted to rise from 0.56%, towards the end of 2026 to ~2.00% in five years.
Note: This is simulated data. I wish interest rates were really this low (0.56%) at the time of writing this post.
Is AR Model just a Random Walk?
Sometimes it is difficult to distinguish between a time series that is slightly mean reverting and a time series that does not mean revert at all, like a random walk. We will compare the
ACF for the slightly mean-reverting interest rate series of the last section with a simulated random walk with the same number of observations.
We should notice when plotting the autocorrelation of these two series side-by-side that they look very similar if the AR model is Random Walk.
# Generate random steps with mean=0 and standard deviation=1 simulated_data = np.random.normal(loc=0, scale=1, size=annual_interest_rate_data.shape)
# Import the plot_acf module from statsmodels from statsmodels.graphics.tsaplots import plot_acf # Plot the interest rate series and the simulated random walk series side-by-side fig, axes = plt.subplots(1,2, figsize=(10,5)) # Plot the autocorrelation of the interest rate series in the top plot fig = plot_acf(annual_interest_rate_data, alpha=1, lags=12, ax=axes) # Plot the autocorrelation of the simulated random walk series in the bottom plot fig = plot_acf(simulated_data, alpha=1, lags=12, ax=axes) # Label axes axes.set_title("Interest Rate Data") axes.set_title("Simulated Random Walk Data") plt.show()
It appears that the interest rates are indeed mean-reverting and not just a random walk.
Choosing the right model
In practice, we will ordinarily not be told the order of the model that we’re trying to estimate.
Identifying the order of an AR Model
There are two techniques that can help determine the order of the AR model:
- The Partial Autocorrelation Function, and
- The Information Criteria
The Partial Autocorrelation Function measures the incremental benefit of adding another lag.
Imagine running several regressions, where you regress returns on more and more lagged values. The coefficients in the red boxes represent the values of the partial autocorrelation function for different lags. For example, in the bottom row, the coefficient in the red box, phi 4-4, is the lag-4 value of the Partial Autocorrelation Function, and it represents how significant adding a fourth lag is when you already have three lags.
plot_pacf is the
statsmodels function for plotting the partial autocorrelation function. The arguments are the same as that of the
plot_acf module that we saw earlier. The input
x is a
series or array. The argument
lags indicates how many lags of the parital autocorrelation function will be plotted. And the
alpha argument sets the width of the confidence interval.
Comparison of PACF for different AR Models
These plots show the Partial Autocorrelation Function for AR models of different orders.
In the upper left, for an
AR(1) model, only the lag-1 PACF is significantly different from zero. Similarly, for an
AR(2) model, two lags are different from zero, and for and
AR(3), three lags are different from zero. Finally, for White Noise, there are no lags that are significantly different from zero.
The need for
The more parameters in a model, the better the model will fit the data. But this can lead to overfitting of the data. This is where information criteria comes into play. The information criteria adjusts the goodness-of-fit of a model by imposing a penalty based on the number of parameters used. Two common adjusted goodness-of-fit measures are called the
Akaike Information Criterion and the
Bayesian Information Criterion.
AIC and BIC
This is the full output from estimating an
AR(2) model. The
BIC are highlighted in the red box.
To get the AIC and BIC statistics, you follow the same procedure from the last section to fit the data to a model. In the last section, we learned how to get the full output using summary or just the AR parameters using the params attribute. We can also get the AIC or BIC using those attributes.
Where to use information criteria?
In practice, the way to use the information criteria is to fit several models, each with a different number of parameters, and choose the one with the lowest Bayesian information criterion.
Suppose we are given a time series of data, and unknown to us, it was simulated from an AR(3) model. Here is a plot of the BIC when we fit the data to an AR(1) up to an AR(8) model. You can see that the lowest BIC occurs for an AR(3).
Example using PACF to identify the order of AR Model
One useful tool to identify the order of an AR model is to look at the Partial Autocorrelation Function (PACF). In this example, we will simulate two time series, an
AR(1) and an
AR(2), and calculate the sample PACF for each. We will notice that for an AR(1), the PACF should have a significant lag-1 value, and roughly zeros after that. And for an AR(2), the sample PACF should have significant lag-1 and lag-2 values, and zeros after that.
# Import the modules for simulating data and for plotting the PACF from statsmodels.tsa.arima_process import ArmaProcess from statsmodels.graphics.tsaplots import plot_pacf # Simulate AR(1) with phi=+0.6 ma = np.array() ar = np.array([1, -0.6]) AR_object = ArmaProcess(ar, ma) simulated_data_1 = AR_object.generate_sample(nsample=5000) # Plot PACF for AR(1) plot_pacf(simulated_data_1, lags=20) plt.show() # Simulate AR(2) with phi1=+0.6, phi2=+0.3 ma = np.array() ar = np.array([1, -0.6, -0.3]) AR_object = ArmaProcess(ar, ma) simulated_data_2 = AR_object.generate_sample(nsample=5000) # Plot PACF for AR(2) plot_pacf(simulated_data_2, lags=20) plt.show()
Notice that the number of significant lags for the PACF indicate the order of the AR model
Estimate the Order of Model using BIC
Another tool to identify the order of a model is to look at the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). These measures compute the goodness of fit with the estimated parameters, but apply a penalty function on the number of parameters in the model. We will take the
AR(2) simulated data from above section, saved as
simulated_data_2, and compute the BIC as you vary the order,
p, in an
0 to 6.
# Import the module for estimating an ARMA model from statsmodels.tsa.arima_model import ARMA # Fit the data to an AR(p) for p = 0,...,6 , and save the BIC BIC = np.zeros(7) for p in range(7): mod = ARMA(simulated_data_2, order=(p,0)) res = mod.fit() # Save BIC for AR(p) BIC[p] = res.bic # Plot the BIC as a function of p plt.plot(range(1,7), BIC[1:7], marker='o') plt.xlabel('Order of AR Model') plt.ylabel('Bayesian Information Criterion') plt.show()
AR(2), the BIC achieves its minimum at
p=2, which is what we expect.