Sunspots — cyclicality with AR models (complex roots & a damped forecast)¶
statsmodels (frequentist) — when the period is variable, not fixed¶
Annual sunspot counts, 1700–2008 (Wolf/SIDC numbers). The famous ~11-year solar cycle is the textbook example of cyclicality — and it is not seasonality:
| Seasonality (AirPassengers) | Cyclicality (sunspots) | |
|---|---|---|
| period | fixed & known (12 months) | variable & unknown (~9–13 yr, wandering) |
| modelled by | seasonal differencing $(1-B^{12})$, dummies, Fourier | AR terms with complex roots / a stochastic-cycle component |
| ACF | spikes at the seasonal lag | damped sine wave |
| long forecast | repeats forever | damps toward the mean |
So there is no seasonal differencing here. A plain AR(p) generates a pseudo-cycle whenever its characteristic polynomial has complex roots; the angle of those roots sets the cycle's period.
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.ar_model import ar_select_order
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.tsa.stattools import adfuller
d = pd.read_csv('sunspots.csv'); s = pd.Series(d['sunspots'].values, index=d['year'].values)
H = 40; train, test = s.iloc[:-H], s.iloc[-H:] # hold out the last 40 years
print('%d years %d-%d; train %d (%d-%d), test %d (%d-%d)' %
(len(s), s.index[0], s.index[-1], len(train), train.index[0], train.index[-1], len(test), test.index[0], test.index[-1]))
print('ADF p = %.4f (small => stationary, i.e. mean-reverting cycle, NOT a trend => d=0)' % adfuller(train)[1])
sel = ar_select_order(train, maxlag=15, ic='aic') # data-driven AR order
p = max(sel.ar_lags); print('AIC selects AR(%d) (Yule famously used AR(2); Box-Jenkins/Tong AR(9))' % p)
res = ARIMA(train, order=(p, 0, 0)).fit()
309 years 1700-2008; train 269 (1700-1968), test 40 (1969-2008) ADF p = 0.0622 (small => stationary, i.e. mean-reverting cycle, NOT a trend => d=0) AIC selects AR(9) (Yule famously used AR(2); Box-Jenkins/Tong AR(9))
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: An unsupported index was provided. As a result, forecasts cannot be generated. To use the model for forecasting, use one of the supported classes of index. self._init_dates(dates, freq) C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: An unsupported index was provided. As a result, forecasts cannot be generated. To use the model for forecasting, use one of the supported classes of index. self._init_dates(dates, freq) C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: An unsupported index was provided. As a result, forecasts cannot be generated. To use the model for forecasting, use one of the supported classes of index. self._init_dates(dates, freq) C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: An unsupported index was provided. As a result, forecasts cannot be generated. To use the model for forecasting, use one of the supported classes of index. self._init_dates(dates, freq) C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: An unsupported index was provided. As a result, forecasts cannot be generated. To use the model for forecasting, use one of the supported classes of index. self._init_dates(dates, freq)
In [2]:
# complex AR roots -> implied pseudo-cycle period (2*pi / angle of the root nearest the unit circle)
roots = res.arroots # |root|>1 for a stationary AR
cplx = roots[np.abs(roots.imag) > 1e-6]
dom = cplx[np.argmin(np.abs(cplx))] # closest to unit circle = least damped = dominant
period = 2*np.pi/np.abs(np.angle(dom))
print('dominant complex root modulus %.3f, implied cycle PERIOD = %.1f years (the ~11-yr solar cycle)' % (np.abs(dom), period))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ax[0].plot(s.index, s.values, color='steelblue', lw=.9)
ax[0].axvline(test.index[0], color='gray', ls=':', lw=1); ax[0].set_title('Annual sunspot number 1700-2008 (variable-amplitude ~11-yr cycle)')
ax[0].set_xlabel('year'); ax[0].set_ylabel('sunspots')
plot_acf(train, lags=40, ax=ax[1]) # damped sine-wave ACF = cyclicality signature
ax[1].set_title('ACF: a damped oscillation (the cyclicality fingerprint)'); ax[1].set_xlabel('lag (years)')
plt.tight_layout(); plt.savefig('sunspots_py1.png', dpi=120, bbox_inches='tight'); plt.show()
dominant complex root modulus 1.030, implied cycle PERIOD = 10.5 years (the ~11-yr solar cycle)
In [3]:
# out-of-sample forecast: an AR cycle DAMPS toward the mean as the horizon grows
# (integer-year index has no freq -> statsmodels returns a positional forecast index; plot vs test.index = the real years)
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)); mae = np.mean(np.abs(fm.values - test.values))
print('OUT-OF-SAMPLE (40 yr): RMSE = %.1f MAE = %.1f (naive mean RMSE %.1f)' %
(rmse, mae, np.sqrt(np.mean((test.values - train.mean())**2))))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
# (A) AR characteristic roots in the complex plane (complex pair => pseudo-cycle)
th = np.linspace(0, 2*np.pi, 200); ax[0].plot(np.cos(th), np.sin(th), 'k-', lw=.8)
inv = 1/roots # inverse roots: inside unit circle
ax[0].scatter(inv.real, inv.imag, color='firebrick', zorder=3)
ax[0].axhline(0, color='gray', lw=.4); ax[0].axvline(0, color='gray', lw=.4); ax[0].set_aspect('equal')
ax[0].set_title('Inverse AR roots (complex pair → %.1f-yr cycle)' % period); ax[0].set_xlabel('Re'); ax[0].set_ylabel('Im')
# (B) forecast damping
ax[1].plot(train.index, train.values, color='steelblue', lw=.8, label='train')
ax[1].plot(test.index, test.values, 'o-', color='black', ms=3, lw=.8, label='actual (held out)')
ax[1].plot(test.index, fm.values, '--', color='firebrick', lw=1.6, label='forecast (damps to mean)')
ax[1].fill_between(test.index, ci.iloc[:,0], ci.iloc[:,1], color='firebrick', alpha=.2, label='95% PI')
ax[1].axhline(train.mean(), color='gray', ls=':', lw=1); ax[1].set_xlabel('year'); ax[1].set_ylabel('sunspots')
ax[1].set_title('Forecast (RMSE %.0f) — the cycle amplitude decays' % rmse); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('sunspots_py2.png', dpi=120, bbox_inches='tight'); plt.show()
OUT-OF-SAMPLE (40 yr): RMSE = 37.4 MAE = 29.4 (naive mean RMSE 52.8)
Spectral view — the frequency-domain fingerprint of the cycle¶
The complex AR roots and the ACF both implied an ≈11-year cycle; the spectrum shows it directly. The AR(9) spectral density (a parametric smoothing of the noisy raw periodogram) has a broad peak at ≈1/10.5 yr — exactly the frequency implied by the dominant complex root ( $f=\angle\,\text{root}/2\pi$ ) — plus a weaker harmonic near 5 yr. Crucially it is a bump, not a spike: a stochastic cycle spreads power over a band of frequencies, whereas deterministic seasonality (e.g. AirPassengers) concentrates power in sharp lines at $k/12$. The raw periodogram is an inconsistent (noisy) estimator — Welch or the parametric AR spectrum give the clean peak.
In [4]:
# spectral view: periodogram + AR(p) spectral density = the frequency-domain of the AR cycle
from scipy.signal import periodogram, welch
phi = res.arparams; sig2 = res.params[list(res.param_names).index("sigma2")]
f_pg, Pxx = periodogram(train.values - train.mean(), fs=1.0, detrend="constant")
f_w, Pw = welch(train.values - train.mean(), fs=1.0, nperseg=64, 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
fpk = np.abs(np.angle(dom)) / (2*np.pi)
fig, ax = plt.subplots(figsize=(8.6, 4.7))
ax.semilogy(f_pg, Pxx, color="0.75", lw=.7, label="raw periodogram (noisy)")
ax.semilogy(f_w, Pw, color="goldenrod", lw=1.3, label="Welch (smoothed)")
ax.semilogy(ff, Sar, color="firebrick", lw=2.2, label="AR(%d) spectral density" % p)
ax.axvline(fpk, color="steelblue", ls="--", lw=1.5, label="complex-root peak: %.1f-yr cycle" % (1/fpk))
ax.set_xlim(0, 0.5); ax.set_ylim(3e2, 4e4)
ax.set_xlabel("frequency (cycles / year)"); ax.set_ylabel("power (log scale)")
ax.set_title("Sunspots spectrum: a BROAD peak at ~1/%.1f yr = a stochastic CYCLE\n(a fixed seasonal period would be a sharp spike, not a bump)" % (1/fpk))
ax.legend(fontsize=8, loc="upper right")
sec = ax.secondary_xaxis("top", functions=(lambda x: 1/np.maximum(x, 1e-9), lambda pr: 1/np.maximum(pr, 1e-9)))
sec.set_xlabel("period (years)"); sec.set_xticks([2, 3, 5, 10, 20, 50])
plt.tight_layout(); plt.savefig("sunspots_spectrum.png", dpi=120, bbox_inches="tight"); plt.show()
Results¶
- Stationary, no differencing. Unlike AirPassengers, sunspots are mean-reverting — the series is bounded and clearly cyclical, and the ADF test is borderline (p≈0.06, rejecting a unit root at the 10% but not the 5% level) — so
d=0— we model the level with autoregression, not differences. - The ACF is a damped sine wave — the signature of cyclicality (a seasonal series instead shows sharp spikes at the fixed seasonal lag).
- Complex AR roots → the cycle. The fitted AR has a complex-conjugate root pair near the unit circle; its angle implies a ~10–11-year period, recovering the solar cycle directly from the autoregression. (Yule's 1927 AR(2) on this very series was the original demonstration; AIC here prefers a higher order that fits the asymmetric rise/fall better.)
- Forecast damping. The out-of-sample forecast tracks the first cycle or two and then decays toward the mean, with the prediction band swelling to the unconditional variance. This is fundamental: a stochastic cycle is only predictable a short way ahead (its phase drifts), whereas deterministic seasonality would repeat indefinitely. RMSE still beats the naive-mean benchmark.
- Cross-checks: PyMC (Bayesian AR, with a posterior over the cycle period) and R (
ar/forecast) follow.