AR Model

by Professor Throckmorton
for Time Series Econometrics
W&M ECON 408/PUBP 616
Slides

AR(\(p\)) Model

• The AR(\(p\)) model is: \(y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_p y_{t-p} + \varepsilon_t\)

• \(\phi_1\), \(\phi_2\), \(\ldots\), \(\phi_p\) are the autoregressive parameters that quantify the influence of past values on the current value

• \(\varepsilon_t\) is white noise error term (or innovation) with mean zero and constant variance, \(\sigma^2\)

AR(\(1\)) Model

• An AR(\(1\)) Model is a special case of the AR(\(p\))

  \[ y_t = \phi y_{t-1} + \varepsilon_t \]

• An AR(\(1\)) is often used to model cyclical dynamics

• e.g., the fluctuations in an AR(1) series where \(\phi = 0.95\) appear much more persistent (i.e., smoother) than those where \(\phi = 0.4\)

AR(\(1\)) for any mean

• An AR(\(1\)) model is covariance stationary if \(|\phi| < 1\)

• Suppose \(E[y_t] = E[y_{t-1}] \equiv \bar{y}\)

• Then \(E[y_t] = 0\) since \(E[y_t] = \phi E[y_{t-1}] + E[\varepsilon_t] \rightarrow (1-\phi)\bar{y} = 0 \rightarrow \bar{y} = 0\)

• Instead consider the AR(1) model with a generic mean: \(y_t - \mu = \phi(y_{t-1} - \mu) + \varepsilon_t\)

• Then \(E[y_t] = \mu\) since

  \[\begin{align*} E[y_t] - \mu &= \phi E[y_{t-1}] -\phi\mu + E[\varepsilon_t] \\ \rightarrow (1-\phi)\bar{y} &= (1-\phi)\mu \\ \rightarrow \bar{y} &= \mu \end{align*}\]

• You can rearrange the (any mean) AR(1) model as \(y_t = (1-\phi)\mu + \phi y_{t-1} + \varepsilon_t\)

AR(1) as MA(\(\infty\))

• The AR(1) model structure is the same at all points in time

  \[\begin{align*} y_t &= \phi y_{t-1} + \varepsilon_t \\ y_{t-1} &= \phi y_{t-2} + \varepsilon_{t-1} \\ y_{t-2} &= \phi y_{t-3} + \varepsilon_{t-2} \end{align*}\]

• Combine, i.e., recursively substitute, to get \(y_t = \phi^3 y_{t-3} + \phi^2 \varepsilon_{t-2}+ \phi \varepsilon_{t-1} + \varepsilon_t\)

• Rinse and repeat to get \(y_t = \sum_{i=0}^\infty \phi^i \varepsilon_{t-i}\) (i.e., the MA(\(\infty\)) Model or Wold Representation)

• And an MA(\(\infty\)) is invertible back to an AR(\(1\)) so long as \(|\phi| < 1\).

• An AR model that can be written as an MA(\(\infty\)) is causal (or exhibits causality).

• That property differentiates it from an invertible MA model, which can be expressed as an AR(\(\infty\)) model.

Variance

• Suppose \(\mu=0\) without loss of generality.

• Given \(Var(\varepsilon_t) = \sigma^2\) and \(|\phi| < 1\), then \(\gamma(0) \equiv Var(y_t)\)

  \[\begin{align*} &= Var\left(\sum_{i=0}^\infty \phi^i \varepsilon_{t-i}\right) \\ &= Var(\varepsilon_t + \phi \varepsilon_{t-1} + \phi^2 \varepsilon_{t-2} + \phi^3 \varepsilon_{t-3} + \phi^4 \varepsilon_{t-4} + \cdots) \\ &= \sigma^2 + \phi^2 \sigma^2 + (\phi^2)^2 \sigma^2 + (\phi^3)^2 \sigma^2 + (\phi^4)^2 \sigma^2 + \cdots \\ &= \sigma^2(1 + \phi^2 + (\phi^2)^2 + (\phi^2)^3 + (\phi^2)^4 + \cdots ) \\ &= \sigma^2/(1-\phi^2) \end{align*}\]

  See Geometric Series.

• \(Std(y_t) = \sqrt{\sigma^2/(1-\phi^2)}\)

Autocorrelation

• \(\gamma(1) = E[y_t y_{t-1}]\)

  \[\begin{align*} &= E[(\phi y_{t-1} + \varepsilon_t)y_{t-1}] \\ &= E[\phi y_{t-1}^2] + E[\varepsilon_t y_{t-1}] \\ &= \phi Var[y_{t-1}] + 0 \\ &= \phi \gamma(0) \\ &= \phi \sigma^2/(1-\phi^2) \end{align*}\]

• \(\rho(1) = \gamma(1)/\gamma(0) = \phi\)

• \(\gamma(2) = E[y_t y_{t-2}] = E[(\phi^2y_{t-2} + \varepsilon_t + \phi \varepsilon_{t-1})y_{t-2}] = \phi^2 E[y_{t-2}^2] = \phi^2 \gamma(0)\)

• \(\rho(2) \equiv \gamma(2)/\gamma(0) = \phi^2\)

• In general, \(\rho(k) = \phi^k\), which is a very nice looking ACF.

• Thus, for an AR(\(1\)) process/model we can refer to \(\phi\) as the autocorrelation (AC) parameter/coefficient.

AR Model as DGP

Let’s compare a couple different AR models to their underlying white noise.

1. White noise:

   \[ y_t = \mu + \varepsilon_t \]

2. AR(1), low AC:

   \[ y_t = (1-\phi_l)\mu + \phi_l y_{t-1} + \varepsilon_t \]

3. AR(1), high AC:

   \[ y_t = (1-\phi_h)\mu + \phi_h y_{t-1} + \varepsilon_t \]

# Libraries
import numpy as np
# Assign parameters
T = 501; mu = 2; phi_l = 0.5; phi_h = 0.9;
# Draw random numbers
rng = np.random.default_rng(seed=0)
eps = rng.standard_normal(T)
# Simulate time series
y = np.zeros([T,3])
for t in range(1,T):
    y[t,0] = mu + eps[t]
    y[t,1] = (1-phi_l)*mu + phi_l*y[t-1,1] + eps[t]
    y[t,2] = (1-phi_h)*mu + phi_h*y[t-1,2] + eps[t]

# Libraries
import pandas as pd
# Convert to DataFrame and plot
df = pd.DataFrame(y)
df.columns = ['White Noise','AR(1) low AC','AR(1) high AC']

# libraries
import matplotlib.pyplot as plt
# Plot variables
fig, axs = plt.subplots(len(df.columns),figsize=(6.5,3.5))
for idx, ax in enumerate(axs.flat):
    col = df.columns[idx]; y = df[col]
    ax.plot(y); ax.set_ylabel(col)
    ax.grid(); ax.autoscale(tight=True); ax.label_outer()
    print(f'{col}: \
E(y) = {y.mean():.2f}, Std(y) = {y.std():.2f}, \
AC(y) = {y.corr(y.shift(1)):.2f}')

White Noise: E(y) = 1.97, Std(y) = 1.02, AC(y) = -0.02
AR(1) low AC: E(y) = 1.94, Std(y) = 1.16, AC(y) = 0.48
AR(1) high AC: E(y) = 1.72, Std(y) = 2.18, AC(y) = 0.88

# Auto-correlation function
from statsmodels.graphics.tsaplots import plot_acf as acf
# Plot Autocorrelation Function
fig, ax = plt.subplots(figsize=(6.5,1.5))
acf(df['AR(1) high AC'],ax)
ax.grid(); ax.autoscale(tight=True);

# ARIMA model
from statsmodels.tsa.arima.model import ARIMA
# Define model
mod = ARIMA(df['AR(1) high AC'],order=(1,0,0))
# Estimate model
res = mod.fit()
summary = res.summary()
# Print summary tables
tab0 = summary.tables[0].as_html()
tab1 = summary.tables[1].as_html()
tab2 = summary.tables[2].as_html()
# print(tab0); print(tab1); print(tab2)

Estimation Results

Dep. Variable:	AR(1) high AC	No. Observations:	501
Model:	ARIMA(1, 0, 0)	Log Likelihood	-718.607

	coef	std err	z	P>|z|	[0.025	0.975]
const	1.6774	0.390	4.302	0.000	0.913	2.442
ar.L1	0.8843	0.021	41.211	0.000	0.842	0.926
sigma2	1.0282	0.059	17.442	0.000	0.913	1.144

Heteroskedasticity (H):	1.14	Skew:	-0.15
Prob(H) (two-sided):	0.41	Kurtosis:	3.45

• Parameter estimates are close to DGP and significantly different from \(0\)

• Residuals have skewness near \(0\) and kurtosis near \(3\) so they are probably normally distributed

Modeling UR as AR(\(1\))

• The U.S. unemployment rate is very persistent.

• What if we model it with an AR(1)?

# Libraries
from fredapi import Fred
import pandas as pd
# Read Data
fred_api_key = pd.read_csv('fred_api_key.txt', header=None)
fred = Fred(api_key=fred_api_key.iloc[0,0])
data = fred.get_series('UNRATE').to_frame(name='UR')
data.index = pd.DatetimeIndex(data.UR.index.values,freq=data.UR.index.inferred_freq)
# Plot
fig, ax = plt.subplots(figsize=(6.5,2));
ax.plot(data.UR); ax.set_title('U.S. Unemployment Rate');
ax.yaxis.set_major_formatter('{x:,.1f}%')
ax.grid(); ax.autoscale(tight=True)

# Plot Autocorrelation Function
fig, ax = plt.subplots(figsize=(6.5,2.5))
acf(data.UR,ax);
ax.grid(); ax.autoscale(tight=True)

• The ACF decays very slow \(\rightarrow\) UR might be non-stationary.

• Let’s assume it is stationary for a moment and estimate an AR(\(1\)) model

# ARIMA model
from statsmodels.tsa.arima.model import ARIMA
# Define model
mod = ARIMA(data.UR,order=(1,0,0))
# Estimate model
res = mod.fit()
summary = res.summary()
# Print summary tables
tab0 = summary.tables[0].as_html()
tab1 = summary.tables[1].as_html()
tab2 = summary.tables[2].as_html()
#print(tab0); print(tab1); print(tab2)

Estimation Results

Dep. Variable:	VALUE	No. Observations:	925
Model:	ARIMA(1, 0, 0)	Log Likelihood	-496.551

	coef	std err	z	P>|z|	[0.025	0.975]
const	5.5470	0.761	7.289	0.000	4.055	7.039
ar.L1	0.9704	0.008	122.023	0.000	0.955	0.986
sigma2	0.1708	0.001	159.798	0.000	0.169	0.173

Heteroskedasticity (H):	6.35	Skew:	16.87
Prob(H) (two-sided):	0.00	Kurtosis:	428.90

• Estimate of AR coeffecient is very high at \(0.9704\), but the data is not a random walk

• Residuals are skewed and have excess kurtosis \(\rightarrow\) model is missing nonlinearity

Actual vs. Simulated Data

• What if we simulated data from AR(1) model setting the parameters to the previous estimates?

• What does a simulation from this estimated model look like?

# Assign parameters
T = len(data.UR); mu = 5.547; rho = 0.9704; sigma = np.sqrt(0.1708);
# Draw random numbers
rng = np.random.default_rng(seed=0)
eps = rng.standard_normal(T)
# Simulate time series
URsim = mu*np.ones([T,1])
for t in range(1,T):
    URsim[t,0] = (1-rho)*mu + rho*URsim[t-1,0] + sigma*eps[t];

# Plot Autocorrelation Function
fig, ax = plt.subplots(figsize=(6.5,1.5))
acf(URsim,ax)
ax.grid(); ax.autoscale(tight=True);

• This ACF approaches zero at shorter lags compared to the actual data.

• We know an AR(\(1\)) model with \(\phi = 0.9704\) generates stationary time series.

• Maybe the actual UR isn’t stationary?

# Add simulation to DataFrame
data['URsim'] = URsim
#  Variables to plot and ylabel
plotme = [['UR','Actual UR'],['URsim','Simulated UR']]
# Plot variables
fig, axs = plt.subplots(len(plotme),figsize=(6.5,3))
for idx, ax in enumerate(axs.flat):
    y = data[plotme[idx][0]]; ylabel = plotme[idx][1]
    ax.plot(y); ax.set_ylabel(ylabel)
    ax.yaxis.set_major_formatter('{x:.0f}%')
    ax.grid(); ax.autoscale(tight=True); ax.label_outer()
    print(f'{ylabel}: \
E(y) = {y.mean():.2f}, Std(y) = {y.std():.2f}, \
AC(y) = {y.corr(y.shift(1)):.2f}')

Actual UR: E(y) = 5.67, Std(y) = 1.71, AC(y) = 0.97
Simulated UR: E(y) = 5.21, Std(y) = 1.35, AC(y) = 0.95

Does it fit?

• The actual UR is very smooth, but the AR(1) model is not.

• The simulated UR often falls below \(2\%\), whereas the actual UR does not.

• The actual UR meets or exceeds \(10\%\) a few times, whereas the simulated UR does not.

• The actual UR rises quickly and falls slowly, whereas the simulated UR seems to rise and fall at roughly the same rates.

• This all suggests that the AR(1) model is unable to reproduce some important features of the actual UR.

Partial ACF

• Most of the slides in this chapter are about the AR(\(1\)) model.

• Suppose the true DGP is an AR model, but we do not know the order.

• When the DGP is a MA model, we can use the ACF to determine the order since it drops off to zero after a certain lag.

• But for AR model, the autocorrelation at higher lags is influenced by the intermediate observations between \(t\) and \(t-k\), and even for an AR(\(1\)) model the ACF tapers slowly.

• Alternatively, Partial ACF (PACF) helps to identify the order of an AR model.

• The PACF at lag \(k\) is the correlation between \(y_t\) and \(y_{t-k}\) that is not explained by observations in between, \(\{y_{t-k+1}:y_{t-1}\}\)

• The PACF at lag \(k\) is denoted as \(\phi_{kk}\), and plotting \(\phi_{kk}\) against lag \(k\) forms the PACF plot.

AR(\(2\)) Model

• Let’s simulate the following AR(\(2\)) model

  \[ y_t = (1-\phi_1-\phi_2)\mu + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \sigma \varepsilon_t \]

  where \(\mu = 2\), \(\phi_1 = 0.8\), \(\phi_2 = 0.15\), and \(\sigma = 0.5\)

• We know the DGP and the true order of the AR model. Does the PACF suggest the same order?

# Assign parameters
T = 501; mu = 2; phi_1 = 0.8; phi_2 = 0.15; sigma = 0.5;
# Draw random numbers
rng = np.random.default_rng(seed=0)
eps = rng.standard_normal(T)
# Simulate time series
y = mu*np.ones([T,1])
for t in range(2,T):
    y[t,0] = (1-phi_1-phi_2)*mu + phi_1*y[t-1,0] + phi_2*y[t-2,0] + sigma*eps[t]

# Convert to DataFrame and plot
df = pd.DataFrame(y)
fig, ax = plt.subplots(figsize=(6.5,2.5))
ax.plot(df); ax.set_title('AR(2) Simulated Data')
ax.grid(); ax.autoscale(tight=True)

# Partial auto-correlation function
from statsmodels.graphics.tsaplots import plot_pacf as pacf
# Plot Autocorrelation Function
fig, ax = plt.subplots(figsize=(6.5,1.5))
pacf(df,ax)
ax.grid(); ax.autoscale(tight=True);

• The PACF is significantly different from zero for \(k = \{1,2\}\), which indicates the DGP is probably an AR(\(2\)).

• So the PACF is good at deterimining order so long as AR coefficients are large enough.

Modeling UR as AR(\(2\))

• Is a higher order AR model a good fit for the U.S. unemployment rate?

• Let’s first plot the PACF of the UR to see if we can determine an order.

# Plot Autocorrelation Function
fig, ax = plt.subplots(figsize=(6.5,1.5))
pacf(data.UR,ax);
ax.grid(); ax.autoscale(tight=True)

This suggests that an AR(1) might be the best model.

• What if we ignore the PACF and estimate an AR(2) anyway?

# Define model
mod = ARIMA(data.UR,order=(2,0,0))
# Estimate model
res = mod.fit()
summary = res.summary()
# Print summary tables
tab0 = summary.tables[0].as_html()
tab1 = summary.tables[1].as_html()
tab2 = summary.tables[2].as_html()
#print(tab0); print(tab1); print(tab2)

Estimation Results

Dep. Variable:	VALUE	No. Observations:	925
Model:	ARIMA(2, 0, 0)	Log Likelihood	-495.314

	coef	std err	z	P>|z|	[0.025	0.975]
const	5.5611	0.740	7.510	0.000	4.110	7.012
ar.L1	1.0208	0.009	109.007	0.000	1.002	1.039
ar.L2	-0.0517	0.010	-5.057	0.000	-0.072	-0.032
sigma2	0.1703	0.001	143.540	0.000	0.168	0.173

Heteroskedasticity (H):	6.37	Skew:	16.65
Prob(H) (two-sided):	0.00	Kurtosis:	423.48

• Both AR coefficients are significantly different from \(0\)

• But the AR(2) is still linear, i.e., normal and symmetric, so the residuals are still skewed and fat-tailed

• PACF did not help to determine that second lag might improve fit

UR First Difference

# Year-over-year growth rate
dUR = data.UR.diff(1); dUR = dUR.dropna()
# Plot
fig, ax = plt.subplots(figsize=(6.5,1.5));
ax.plot(dUR); ax.set_title('Change in U.S. Unemployment Rate (month to month)');
ax.yaxis.set_major_formatter('{x:,.0f}pp'); ax.grid(); ax.autoscale(tight=True)

Is it stationary?

• Yes, mean looks constant around \(0\).

• Yes, volatility looks constant, except for COVID.

• Maybe not, it’s impossible to see if the autocovariance is constant.

# Plot Autocorrelation Function
fig, ax = plt.subplots(figsize=(6.5,1.5))
acf(dUR,ax)
ax.grid(); ax.autoscale(tight=True);

• The ACF indicates that the series is probably stationary.

• What does the PACF say?

# Plot Autocorrelation Function
fig, ax = plt.subplots(figsize=(6.5,1.5))
pacf(dUR,ax);
ax.grid(); ax.autoscale(tight=True)

• The PACF shows that the autocorrelation is not significantly different from zero for all lags.

• What is the best model to estimate? Let’s try MA(\(1\)) and AR(\(1\)) and compare.

Estimation Results

# Define model
mod0 = ARIMA(dUR,order=(0,0,1))
mod1 = ARIMA(dUR,order=(1,0,0))
# Estimate model
res = mod1.fit(); summary = res.summary()
# Print summary tables
tab0 = summary.tables[0].as_html()
tab1 = summary.tables[1].as_html()
tab2 = summary.tables[2].as_html()
#print(tab0); print(tab1); print(tab2)

Dep. Variable:	VALUE	No. Observations:	924
Model:	ARIMA(0, 0, 1)	Log Likelihood	-500.699

	coef	std err	z	P>|z|	[0.025	0.975]
const	0.0007	0.025	0.027	0.979	-0.047	0.049
ma.L1	0.0398	0.008	4.893	0.000	0.024	0.056
sigma2	0.1731	0.001	185.366	0.000	0.171	0.175

Heteroskedasticity (H):	6.54	Skew:	16.40
Prob(H) (two-sided):	0.00	Kurtosis:	419.10

• The MA coefficient estimate is statistically significantly different from \(0\).

• However, the residuals are still probably heteroskedastic and not normally distributed.

• Thus, a linear/normal model is probably not the best choice.

Dep. Variable:	VALUE	No. Observations:	924
Model:	ARIMA(1, 0, 0)	Log Likelihood	-500.757

	coef	std err	z	P>|z|	[0.025	0.975]
const	0.0007	0.025	0.026	0.979	-0.048	0.049
ar.L1	0.0363	0.008	4.344	0.000	0.020	0.053
sigma2	0.1731	0.001	185.413	0.000	0.171	0.175

Heteroskedasticity (H):	6.55	Skew:	16.42
Prob(H) (two-sided):	0.00	Kurtosis:	419.51

• The AR coefficient is small/weak but estimate is statistically significantly different from \(0\).

• Like the MA(\(1\)) model, the residuals are still probably heteroskedastic and not normally distributed.

• This again confirms that a linear/normal model is probably not the best choice.