Sunspots — R (ar Yule-Walker + forecast), cyclicality¶

Companion to the Python/PyMC notebooks. R's classic ar() fits an autoregression by Yule-Walker, choosing the order by AIC; forecast() produces the out-of-sample forecast. We read the implied cycle period from the AR polynomial's complex roots, and plot with ggplot/autoplot.

In [1]:
suppressMessages({library(forecast); library(ggplot2)})
d <- read.csv('sunspots.csv'); y <- ts(d$sunspots, start=d$year[1], frequency=1)
H <- 40; train <- window(y, end=d$year[1]+length(y)-H-1); test <- window(y, start=d$year[1]+length(y)-H)
fit <- ar(train, order.max=15, aic=TRUE)                    # Yule-Walker, AIC order
cat(sprintf('ar() selected order = AR(%d)\n', fit$order))
rts <- polyroot(c(1, -fit$ar)); cz <- rts[abs(Im(rts)) > 1e-6]
dom <- cz[which.min(Mod(cz))]; cat(sprintf('implied cycle period = %.1f years\n', 2*pi/abs(Arg(dom))))
fc <- forecast(fit, h=H)
cat(sprintf('OUT-OF-SAMPLE RMSE = %.1f   (naive-mean RMSE %.1f)\n',
            sqrt(mean((as.numeric(fc$mean)-as.numeric(test))^2)), sqrt(mean((as.numeric(test)-mean(train))^2))))

ar() selected order = AR(9)

implied cycle period = 10.4 years

OUT-OF-SAMPLE RMSE = 39.6   (naive-mean RMSE 52.8)
In [2]:
options(repr.plot.width = 9, repr.plot.height = 5)
autoplot(fc) +
  autolayer(test, series='actual (held out)') +
  geom_hline(yintercept=mean(train), linetype='dotted', colour='grey40') +
  scale_colour_manual(values=c('actual (held out)'='red')) +
  labs(title=sprintf('ar() AR(%d) out-of-sample forecast — cycle damps to mean', fit$order),
       x='year', y='sunspots', colour='') +
  theme_minimal(base_size=13) + theme(legend.position='top')
No description has been provided for this image

Results¶

  • ar() (Yule-Walker) independently selects AR(9) and recovers the ~10–11-year cycle from the complex roots of the AR polynomial — matching statsmodels and the Bayesian fit.
  • forecast() shows the same damping toward the mean: a stochastic cycle is only predictable a short way out. RMSE beats the naive-mean benchmark.
  • Three engines agree on the cyclicality story. This contrasts cleanly with AirPassengers' fixed seasonality — the bridge from seasonality to the GARCH/volatility work to come.