US GNP growth — the business cycle with ARMA, out-of-sample

statsmodels (frequentist) — when the "cycle" is irregular and barely predictable

Quarterly US real GNP growth, 1951–1984 (Hamilton's 1989 series). The business cycle — expansions punctuated by recessions (shaded) — is the macroeconomic example of cyclicality: the swings have no fixed period (recessions arrive irregularly). Three things make this case different from the others in the arc:

• growth is already stationary (it's a difference of log GNP), so we model it directly with an ARMA (d=0);
• the autocorrelation is weak — GNP growth is close to unforecastable, the empirical basis for "GDP growth is hard to predict";
• the linear ARMA is symmetric, but real recessions are sharper than expansions — which is exactly why Hamilton introduced the Markov-switching model (see the companion regime-switching project). ARMA is the linear benchmark that motivates it.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
import pmdarima as pm

d = pd.read_csv('hamilton_gnp.csv', parse_dates=['date']).dropna(subset=['gnp_growth']).reset_index(drop=True)
g = pd.Series(d['gnp_growth'].values, index=d['date'].values); rec = d['recession'].values
H = 16; train, test = g.iloc[:-H], g.iloc[-H:]                  # hold out last 4 years
print('%d quarters %d..%d; test %d..%d' % (len(g), g.index[0].year, g.index[-1].year, test.index[0].year, test.index[-1].year))
print('mean growth %.2f%%/qtr; ADF p = %.4f (stationary => d=0); lag-1 autocorr = %.2f' %
      (g.mean(), adfuller(train)[1], pd.Series(train).autocorr(1)))
am = pm.auto_arima(train, seasonal=False, d=0, stepwise=True, suppress_warnings=True, error_action='ignore')
print('auto_arima -> ARMA%s   AIC %.1f' % (am.order, am.aic()))
res = ARIMA(train, order=am.order).fit()

135 quarters 1951..1984; test 1981..1984
mean growth 0.74%/qtr; ADF p = 0.0000 (stationary => d=0); lag-1 autocorr = 0.31

auto_arima -> ARMA(1, 0, 0)   AIC 343.1

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.
  self._init_dates(dates, freq)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.
  self._init_dates(dates, freq)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.
  self._init_dates(dates, freq)

def shade(ax):                                                  # NBER recession bands
    inr=False
    for i,r in enumerate(rec):
        if r and not inr: x0=g.index[i]; inr=True
        if inr and (not r or i==len(rec)-1): ax.axvspan(x0, g.index[i], color='gray', alpha=.2); inr=False
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].plot(g.index, g.values, color='steelblue', lw=.9); ax[0].axhline(g.mean(), color='firebrick', ls=':', lw=1)
shade(ax[0]); ax[0].set_title('Real GNP growth (%/qtr), NBER recessions shaded'); ax[0].set_ylabel('growth %'); ax[0].set_xlabel('year')
plot_acf(train, lags=16, ax=ax[1]); ax[1].set_title('ACF — weak, short-lived persistence'); ax[1].set_xlabel('lag (quarters)')
plt.tight_layout(); plt.savefig('gnp_py1.png', dpi=120, bbox_inches='tight'); plt.show()

ar = res.arroots if hasattr(res,'arroots') else np.array([])
cplx = ar[np.abs(ar.imag)>1e-6] if len(ar) else np.array([])
if len(cplx): print('implied business-cycle period from complex AR roots ~ %.1f quarters (%.1f years)' %
                     (2*np.pi/np.abs(np.angle(cplx[np.argmin(np.abs(cplx))])), 2*np.pi/np.abs(np.angle(cplx[np.argmin(np.abs(cplx))]))/4))
else: print('AR roots are real (no clean pseudo-cycle) — persistence without a sharp period')

AR roots are real (no clean pseudo-cycle) — persistence without a sharp period

# out-of-sample forecast: growth reverts to its mean within a couple of quarters (low predictability)
fc = res.get_forecast(H); fm = fc.predicted_mean; ci = fc.conf_int(alpha=0.05)
rmse = np.sqrt(np.mean((fm.values-test.values)**2)); rmse0 = np.sqrt(np.mean((test.values-train.mean())**2))
print('OUT-OF-SAMPLE (16 qtr):  ARMA RMSE = %.2f   naive-mean RMSE = %.2f  (barely better => growth is hard to forecast)' % (rmse, rmse0))

fig, ax = plt.subplots(figsize=(11, 4.4))
ax.plot(train.index, train.values, color='steelblue', lw=.9, label='train')
ax.plot(test.index, test.values, 'o-', color='black', ms=3, lw=.8, label='actual (held out)')
ax.plot(fm.index, fm.values, '--', color='firebrick', lw=1.8, label='forecast')
ax.fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], color='firebrick', alpha=.2, label='95% PI')
ax.axhline(train.mean(), color='gray', ls=':', lw=1, label='mean growth'); shade(ax)
ax.set_xlabel('year'); ax.set_ylabel('growth %'); ax.set_title('ARMA forecast (RMSE %.2f vs naive %.2f) — reverts to mean fast'%(rmse,rmse0)); ax.legend(fontsize=8, ncol=2)
plt.tight_layout(); plt.savefig('gnp_py2.png', dpi=120, bbox_inches='tight'); plt.show()

OUT-OF-SAMPLE (16 qtr):  ARMA RMSE = 1.20   naive-mean RMSE = 1.20  (barely better => growth is hard to forecast)

Why a linear ARMA isn't the whole story — business-cycle asymmetry

Hamilton's motivation for regime-switching: recessions are not merely lower growth, they are deeper and sharper. Split by the NBER indicator, recession growth averages ≈ −0.6%/qtr vs ≈ +1.1% in expansions, and the overall distribution is left-skewed (skew ≈ −0.46) — rare, abrupt contractions that a symmetric, constant-parameter ARMA cannot represent. This is exactly the nonlinearity the Markov-switching notebooks below capture.

# business-cycle ASYMMETRY: recessions are deeper (left skew) -> a symmetric ARMA can't capture it
from scipy.stats import skew, norm
gv = g.values; rc = rec.astype(bool)
exp_g, rec_g = gv[~rc], gv[rc]; sk = skew(gv)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4)); bins = np.linspace(gv.min(), gv.max(), 22)
ax[0].hist(exp_g, bins=bins, color="steelblue", alpha=.6, density=True, label="expansion (mean %.2f)" % exp_g.mean())
ax[0].hist(rec_g, bins=bins, color="firebrick", alpha=.6, density=True, label="recession (mean %.2f)" % rec_g.mean())
ax[0].axvline(0, color="k", lw=.7); ax[0].set_xlabel("GNP growth (%/qtr)"); ax[0].set_ylabel("density")
ax[0].set_title("Asymmetry: recessions are deeper, not just lower\n(a symmetric ARMA cannot represent this)"); ax[0].legend(fontsize=8)
xs = np.linspace(gv.min(), gv.max(), 200)
ax[1].hist(gv, bins=22, color="slateblue", alpha=.6, density=True, label="growth (skew %.2f)" % sk)
ax[1].plot(xs, norm.pdf(xs, gv.mean(), gv.std()), "k--", lw=1.4, label="matched normal (skew 0)")
ax[1].axvline(gv.mean(), color="firebrick", ls=":", lw=1); ax[1].set_xlabel("GNP growth (%/qtr)"); ax[1].set_ylabel("density")
ax[1].set_title("Left-skew: rare sharp contractions\n(the nonlinearity that motivates Markov switching)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig("gnp_asymmetry.png", dpi=120, bbox_inches="tight"); plt.show()

Spectral view — cyclicality as a broad band, not a line

The GNP-growth spectrum is broad and gently declining from low frequency (growth is close to white noise, echoing the weak ACF and poor forecastability): about half the power sits in the business-cycle band (periods 6–32 quarters). There is no sharp spectral line — the business cycle has no fixed period, the frequency-domain signature of cyclicality (contrast the AirPassengers / CO₂ seasonal combs).

# spectral view: GNP growth is near white noise with a gentle business-cycle tilt (no sharp line)
from scipy.signal import periodogram, welch
gv = g.values
r4 = ARIMA(gv, order=(1, 0, 0), trend="c").fit(); phi = r4.arparams; sig2 = r4.params[list(r4.param_names).index("sigma2")]
f_pg, Pxx = periodogram(gv, fs=1.0, detrend="constant"); f_w, Pw = welch(gv, fs=1.0, nperseg=32, detrend="constant")
ff = np.linspace(1e-3, 0.5, 600); kk = np.arange(1, len(phi) + 1)
Sar = sig2 / np.abs(1 - (phi[None, :] * np.exp(-2j*np.pi*ff[:, None]*kk[None, :])).sum(1))**2
lo, hi = 1/32, 1/6; mass = Pxx[(f_pg >= lo) & (f_pg <= hi)].sum() / Pxx[f_pg > 0].sum()
fig, ax = plt.subplots(figsize=(8.8, 4.7))
ax.axvspan(lo, hi, color="goldenrod", alpha=.18, label="business-cycle band (6-32 qtr)")
ax.semilogy(f_pg, Pxx, color="0.78", lw=.6, label="raw periodogram"); ax.semilogy(f_w, Pw, color="goldenrod", lw=1.3, label="Welch")
ax.semilogy(ff, Sar, color="firebrick", lw=2.2, label="AR(1) spectral density")
ax.set_xlim(0, 0.5); ax.set_ylim(2e-1, 3e0)
ax.set_xlabel("frequency (cycles / quarter)"); ax.set_ylabel("power (log scale)")
ax.set_title("GNP-growth spectrum: gently declining from low frequency (business cycle), no sharp line\nnear white noise - %.0f%% of power in the 6-32-qtr band" % (100*mass))
ax.legend(fontsize=8, loc="upper right")
sec = ax.secondary_xaxis("top", functions=(lambda x: 1/np.maximum(x, 1e-9), lambda p: 1/np.maximum(p, 1e-9)))
sec.set_xlabel("period (quarters)"); sec.set_xticks([2, 3, 4, 6, 8, 16, 32])
plt.tight_layout(); plt.savefig("gnp_spectrum.png", dpi=120, bbox_inches="tight"); plt.show()

Results

• Stationary growth, weak persistence. GNP growth is stationary (ADF rejects a unit root) with only mild positive autocorrelation; auto_arima picks a low-order ARMA. The ACF dies out within a few quarters.
• Recessions are visible but irregular — the business "cycle" has no fixed length, which is why we model it as autoregression rather than seasonality.
• Growth is barely forecastable. Out of sample, the ARMA's RMSE only slightly beats simply predicting the mean, and the forecast reverts to the mean within ~2 quarters. This is the honest empirical finding behind "nobody can forecast GDP growth" — the opposite extreme from CO₂.
• Why regime-switching exists. A linear ARMA treats expansions and recessions symmetrically. Real recessions are sharper and shorter than expansions; capturing that asymmetry is exactly what the Markov-switching AR model adds on this very series. ARMA is the linear benchmark it improves on.
• Cross-checks: PyMC (Bayesian ARMA, posterior persistence) and R (forecast::auto.arima) follow.

Connection — from this linear AR to Markov-switching

This AR/ARMA is the linear, single-regime benchmark for the same Hamilton GNP series modelled with regimes:

• markov_switching_python.ipynb / markov_switching_R.ipynb — frequentist MS-AR (statsmodels / MSwM)
• bayesian_msar_python.ipynb — Bayesian MS-AR (NumPyro)
• mts_python.ipynb — the McCulloch–Tsay AR-shift model (a close relative)

The ladder — each step relaxes "parameters are constant over time":

model	regime structure	latent state
ARIMA (this notebook)	one regime, constant forever	none
changepoint / intervention (nile_pymc.ipynb)	one permanent switch	a single break time τ
Markov-switching (MS-AR)	regimes recur	a hidden state $s_t$ on a Markov chain

A symmetric ARMA can't capture the sharp, recurring recessions visible in the shaded plot above — that asymmetry is exactly what MS-AR adds.