Time-Series Foundations: Stationarity, ARMA, and Forecasting

1. Stationarity, ergodicity and dependence

A process \(\{y_t\}\) is strictly stationary if the joint distribution of \((y_{t},\dots,y_{t+h})\) does not depend on \(t\), for every \(h\). It is weakly (covariance) stationary if \(E[y_t]=\mu\) and \(\operatorname{Cov}(y_t,y_{t-j})=\gamma_j\) do not depend on \(t\). The autocorrelation function (ACF) is \(\rho_j=\gamma_j/\gamma_0\).

Stationarity alone does not justify time averages: \(y_t=Z\) for all \(t\), with \(Z\) random, is stationary, but \(\bar y_T=Z\) never learns \(E[Z]\). Ergodicity rules this out. Loosely, it requires that observations far apart are nearly independent.

These two results replace the WLLN and CLT of Lecture 2. The MDS condition is exactly what correct specification of a dynamic conditional mean delivers: if \(E[u_t\mid\text{past}]=0\), then \(x_tu_t\) is an MDS whenever \(x_t\) contains only past information. If the scores are serially correlated, the long-run variance and HAC standard errors of Lecture 4 are needed.

2. Moving-average and autoregressive processes

Let \(\varepsilon_t\) be white noise: mean zero, variance \(\sigma^2\), uncorrelated over time.

MA(q). \(y_t=\mu+\varepsilon_t+\theta_1\varepsilon_{t-1}+\dots+\theta_q\varepsilon_{t-q}\). It is always stationary, and \(\gamma_j=0\) for \(j>q\): the ACF cuts off after lag \(q\). It is invertible (expressible as an AR(\(\infty\)) in current and past \(y\)) if the roots of \(\theta(z)=1+\theta_1z+\dots+\theta_qz^q\) lie outside the unit circle.

AR(1). \(y_t=c+\phi y_{t-1}+\varepsilon_t\). By repeated substitution, if \(|\phi|<1\),

The ACF decays geometrically. The effect of a shock on \(y_{t+h}\) is \(\phi^h\), and the half-life of a shock is \(\ln0.5/\ln|\phi|\). With \(\phi=0.9\) per quarter, it is about 6.6 quarters.

AR(p). \(\phi(L)y_t=c+\varepsilon_t\), with \(\phi(L)=1-\phi_1L-\dots-\phi_pL^p\) and \(L\) the lag operator. It is stationary if all roots of \(\phi(z)=0\) lie outside the unit circle. Equivalently, all eigenvalues of the companion matrix

lie inside the unit circle. The autocovariances satisfy the Yule–Walker equations \(\gamma_j=\phi_1\gamma_{j-1}+\dots+\phi_p\gamma_{j-p}\) for \(j\ge1\). Complex roots generate damped cycles in the ACF and in impulse responses.

ARMA(p,q). \(\phi(L)y_t=c+\theta(L)\varepsilon_t\): a parsimonious representation of long autocorrelation patterns.

The Wold theorem justifies linear time-series models as approximations: any stationary process has an MA(\(\infty\)) representation, and finite ARMA models are parsimonious approximations to it. But the \(\varepsilon_t\) are forecast errors, not structural shocks. That distinction drives Lecture 11.

3. Model identification

The partial autocorrelation at lag \(j\) is the coefficient on \(y_{t-j}\) in a regression of \(y_t\) on \(y_{t-1},\dots,y_{t-j}\). Under white noise, sample autocorrelations are approximately \(N(0,1/T)\), giving the familiar \(\pm1.96/\sqrt T\) bands. The Ljung–Box statistic \(Q=T(T+2)\sum_{j=1}^m\widehat\rho_j^2/(T-j)\sim\chi^2_{m-p-q}\) tests jointly for remaining residual autocorrelation.

Information criteria. Choose the order minimizing

BIC is consistent: it selects the true finite order with probability approaching one. AIC tends to over-fit but is asymptotically efficient for forecasting when the true model is infinite-order. Compare criteria on a common estimation sample.

4. Estimation and inference

Autoregressions by OLS. For a stationary AR(p) with MDS errors, OLS of \(y_t\) on \((1,y_{t-1},\dots,y_{t-p})\) is consistent and asymptotically normal, by Theorems 9.1–9.2. It is biased in finite samples, because \(E[\varepsilon_t\mid y_{t+1}]\neq0\): strict exogeneity fails although contemporaneous exogeneity holds. In the AR(1), the Kendall bias is approximately \(-(1+3\phi)/T\) (with an intercept). This is non-trivial for persistent series in short samples.

MA terms require nonlinear least squares or maximum likelihood, using the Kalman filter or conditional sums of squares.

A lagged dependent variable with serially correlated errors. Suppose \(y_t=\beta y_{t-1}+u_t\) with \(u_t=\rho u_{t-1}+e_t\). Then \(y_{t-1}\) depends on \(u_{t-1}\), which is correlated with \(u_t\), so OLS is inconsistent. HAC standard errors do not fix this: the problem is in the conditional mean. The remedy is to model the dynamics fully: here \(y_t=(\beta+\rho)y_{t-1}-\beta\rho y_{t-2}+e_t\) is an AR(2) with white-noise errors. Test for remaining serial correlation with the Breusch–Godfrey LM test, which, unlike Durbin–Watson, is valid with lagged dependent variables.

5. Forecasting

Under quadratic loss, the optimal forecast is the conditional mean \(\widehat y_{T+h\mid T}=E[y_{T+h}\mid\mathcal F_T]\). For the AR(1),

Forecasts revert to the mean, and forecast uncertainty grows toward the unconditional variance. Iterating a one-step model ("iterated") and regressing \(y_{t+h}\) on information at \(t\) ("direct") are both common. Iterated forecasts are more efficient if the model is correct; direct forecasts are more robust to misspecification. Reported intervals usually ignore parameter uncertainty, which matters in small samples.

Evaluation. Use pseudo out-of-sample forecasts: estimate on data up to \(t\), forecast \(t+h\), roll forward, and compute the root mean squared forecast error (RMSFE). Compare two forecasts with the Diebold–Mariano test: with loss differential \(d_t=e_{1t}^2-e_{2t}^2\), the statistic \(\bar d/\widehat{se}(\bar d)\to_dN(0,1)\), where the standard error is HAC because \(h\)-step errors are MA(\(h-1\)). For nested models, the DM test is non-standard. Use Clark and West (2007).

import numpy as np
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import acf, pacf

rng = np.random.default_rng(9)
proc = ArmaProcess(ar=[1, -1.2, 0.35], ma=[1])          # AR(2) from problem 1
y = proc.generate_sample(400, distrvs=rng.standard_normal)
print("ACF :", acf(y, nlags=5).round(2))
print("PACF:", pacf(y, nlags=5).round(2))

bic = {p: ARIMA(y[4:], order=(p, 0, 0)).fit().bic for p in range(5)}
p_hat = min(bic, key=bic.get)
print("BIC chooses p =", p_hat)

errs = []
for t in range(300, 399):                               # expanding-window 1-step forecasts
    fc = ARIMA(y[:t], order=(p_hat, 0, 0)).fit().forecast(1)[0]
    errs.append(y[t] - fc)
print("pseudo out-of-sample RMSFE:", np.sqrt(np.mean(np.square(errs))).round(3))

import numpy as np, statsmodels.api as sm
from statsmodels.tsa.ar_model import AutoReg, ar_select_order
y = sm.datasets.macrodata.load_pandas().data["infl"].to_numpy()
# observations 4..103 are 1960Q1..1984Q4; 104..202 are 1985Q1..2009Q3
# TODO: ar_select_order(y[4:104], maxlag=8, ic="bic"); loop over t, fit AutoReg(y[:t], p), forecast y[t]

import numpy as np, statsmodels.api as sm, statsmodels.formula.api as smf
from statsmodels.tsa.ar_model import AutoReg, ar_select_order
from statsmodels.tsa.stattools import acf, pacf
y = sm.datasets.macrodata.load_pandas().data["infl"].to_numpy()
print("ACF  lags 1-6:", acf(y[4:104], nlags=6)[1:].round(2))
print("PACF lags 1-6:", pacf(y[4:104], nlags=6)[1:].round(2))
p = len(ar_select_order(y[4:104], maxlag=8, ic="bic").ar_lags)
print("BIC-selected AR order on 1960-1984:", p)

e_ar, e_rw = [], []
for t in range(104, len(y)):
    fit = AutoReg(y[:t], lags=p, trend="c").fit()
    e_ar.append(y[t] - fit.predict(start=t, end=t)[0])
    e_rw.append(y[t] - y[t - 1])
e_ar, e_rw = np.array(e_ar), np.array(e_rw)
print(f"RMSFE 1985Q1-2009Q3:  AR({p}) {np.sqrt(np.mean(e_ar**2)):.3f}   random walk {np.sqrt(np.mean(e_rw**2)):.3f}")
d = e_ar**2 - e_rw**2
dm = smf.ols("d ~ 1", data={"d": d}).fit(cov_type="HAC", cov_kwds={"maxlags": 4})
print(f"Diebold-Mariano: mean loss differential {d.mean():.3f}, HAC t = {dm.tvalues.iloc[0]:.2f}, p = {dm.pvalues.iloc[0]:.3f}")

ACF  lags 1-6: [0.75 0.72 0.72 0.62 0.52 0.48]
PACF lags 1-6: [ 0.76  0.36  0.31 -0.06 -0.2  -0.1 ]
BIC-selected AR order on 1960-1984: 3
RMSFE 1985Q1-2009Q3:  AR(3) 2.457   random walk 2.913
Diebold-Mariano: mean loss differential -2.450, HAC t = -2.01, p = 0.044

import numpy as np, wooldridge as wd
from statsmodels.tsa.ar_model import AutoReg
u = wd.data("phillips")["unem"].to_numpy()
fit = AutoReg(u, lags=1, trend="c").fit()
# TODO: half-life = ln(0.5)/ln(rho); bias-corrected rho_bc = rho + (1+3 rho)/T

import numpy as np, wooldridge as wd
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.stats.diagnostic import acorr_ljungbox
ph = wd.data("phillips")
u = ph["unem"].to_numpy(); T = len(u) - 1
fit = AutoReg(u, lags=1, trend="c").fit()
c, rho = fit.params
hl = lambda r: np.log(0.5) / np.log(r)
rho_bc = rho + (1 + 3 * rho) / T
print(f"rho = {rho:.3f} (se {fit.bse[1]:.3f}), mean = {c/(1-rho):.2f}%, half-life = {hl(rho):.2f} years")
print(f"bias-corrected rho = {rho_bc:.3f}, half-life = {hl(rho_bc):.2f} years")
lb = acorr_ljungbox(fit.resid, lags=[4], return_df=True)
print(f"Ljung-Box(4) on residuals: Q = {lb['lb_stat'].iloc[0]:.2f}, p = {lb['lb_pvalue'].iloc[0]:.3f}")
s2 = fit.sigma2; last = u[-1]; mu = c / (1 - rho)
for h in range(1, 6):
    f = mu + rho**h * (last - mu); v = s2 * sum(rho**(2*j) for j in range(h))
    print(f"  {int(ph['year'].iloc[-1]) + h}: {f:.2f}%  [{f - 1.96*np.sqrt(v):.2f}, {f + 1.96*np.sqrt(v):.2f}]")

rho = 0.742 (se 0.088), mean = 5.78%, half-life = 2.33 years
bias-corrected rho = 0.801, half-life = 3.12 years
Ljung-Box(4) on residuals: Q = 3.06, p = 0.548
  2004: 5.94%  [4.02, 7.87]
  2005: 5.90%  [3.51, 8.30]
  2006: 5.87%  [3.25, 8.49]
  2007: 5.85%  [3.12, 8.58]
  2008: 5.83%  [3.04, 8.63]