Nile river flow — R break detection + intervention (Arima)¶

Companion to the Python/PyMC notebooks. We detect the level break with a self-contained least-squares scan (the single-break case of strucchange), then fit forecast::Arima with a step-dummy xreg (Box-Tiao intervention) and forecast the held-out years.

In [1]:
suppressMessages({library(forecast); library(ggplot2)})
d <- read.csv('nile.csv'); y <- ts(d$flow, start=d$year[1])
H <- 20; ntr <- length(y)-H; train <- window(y, end=d$year[ntr])   # keep as ts so the forecast carries proper years
# least-squares single-break detection: pick k minimising within-segment RSS
n <- length(train); cand <- 4:(n-4)
rss <- sapply(cand, function(k) sum((train[1:k]-mean(train[1:k]))^2) + sum((train[(k+1):n]-mean(train[(k+1):n]))^2))
kstar <- cand[which.min(rss)]; bk_year <- d$year[kstar]
cat(sprintf('detected break at year %d  (pre mean %.0f, post mean %.0f)\n', bk_year, mean(train[1:kstar]), mean(train[(kstar+1):n])))
step <- as.numeric(d$year >= bk_year)
fit <- Arima(train, order=c(1,0,0), xreg=step[1:ntr], include.mean=TRUE)
cat(sprintf('level-shift coefficient (xreg) = %.0f units   AIC %.1f\n', coef(fit)['xreg'], fit$aic))
fit0 <- Arima(train, order=c(1,1,1)); cat(sprintf('no-intervention ARIMA(1,1,1) AIC %.1f (higher)\n', fit0$aic))
fc <- forecast(fit, h=H, xreg=step[(ntr+1):length(y)])
cat(sprintf('OUT-OF-SAMPLE RMSE = %.0f\n', sqrt(mean((as.numeric(fc$mean)-as.numeric(window(y, start=d$year[ntr+1])))^2))))

detected break at year 1898  (pre mean 1098, post mean 840)

level-shift coefficient (xreg) = -251 units   AIC 1012.5

no-intervention ARIMA(1,1,1) AIC 1016.6 (higher)

OUT-OF-SAMPLE RMSE = 127
In [2]:
options(repr.plot.width = 9, repr.plot.height = 5)
test <- window(y, start=d$year[ntr+1])
autoplot(fc) +
  autolayer(test, series='actual (held out)') +
  geom_vline(xintercept=bk_year, linetype='dashed', colour='firebrick') +
  annotate('text', x=bk_year+1, y=max(y), label=paste('break', bk_year), hjust=0, colour='firebrick') +
  scale_colour_manual(values=c('actual (held out)'='black')) +
  labs(title='Nile: intervention model (AR(1) + step dummy) forecast', x='year', y='flow (10^8 m^3)', colour='') +
  theme_minimal(base_size=13) + theme(legend.position='top')
No description has been provided for this image

A stochastic local level — the break without a dummy¶

StructTS(y, type="level") fits the classic local-level (random-walk-plus-noise) model; its smoothed level adapts through the ~1898–1899 drop on its own — the state-space counterpart to the step dummy, and the R twin of the Python UnobservedComponents fit. Deterministic step (dashed) vs stochastic level (red): two views of one break.

In [3]:
# stochastic local level (StructTS) -- adapts to the 1899 drop, no dummy (cf. the deterministic step)
options(repr.plot.width = 9, repr.plot.height = 5)
ss  <- StructTS(y, type = "level")
lvl <- tsSmooth(ss)[, "level"]
tt  <- as.numeric(time(y)); yy <- as.numeric(y)
pre <- mean(yy[tt < bk_year]); post <- mean(yy[tt >= bk_year])
plot(tt, yy, pch = 20, col = "grey55", xlab = "year", ylab = "flow (10^8 m^3)",
     main = "Nile: a stochastic local level (StructTS) adapts to the 1899 drop")
lines(tt, lvl, col = "firebrick", lwd = 2.5)
lines(c(min(tt), bk_year), c(pre, pre), col = "steelblue", lwd = 2, lty = 2)
lines(c(bk_year, max(tt)), c(post, post), col = "steelblue", lwd = 2, lty = 2)
abline(v = bk_year, lty = 3, col = "grey40")
legend("topright", c("flow", "local level (smoothed)", "deterministic step"),
       col = c("grey55", "firebrick", "steelblue"), pch = c(20, NA, NA),
       lty = c(NA, 1, 2), lwd = c(NA, 2.5, 2), bty = "n", cex = .8)
No description has been provided for this image

Results¶

  • The least-squares scan detects the break at ~1898-1899 — agreeing with the Bayesian changepoint posterior — confirming the Aswan-dam shift directly from the data.
  • Arima with the step xreg estimates the same ~-250-unit level drop, with a lower AIC than a differenced ARIMA(1,1,1).
  • Three engines agree on the break location and its size. This intervention/changepoint toolkit (step dummies, break detection, Bayesian changepoints) is the ARIMA complement to seasonality and cyclicality — and the changepoint idea reappears in regime-switching and volatility models.