Main page: https://www.rdocumentation.org/packages/datasets/versions/3.6.2/topics/Seatbelts
head(Seatbelts)
## DriversKilled drivers front rear kms PetrolPrice VanKilled law
## [1,] 107 1687 867 269 9059 0.1029718 12 0
## [2,] 97 1508 825 265 7685 0.1023630 6 0
## [3,] 102 1507 806 319 9963 0.1020625 12 0
## [4,] 87 1385 814 407 10955 0.1008733 8 0
## [5,] 119 1632 991 454 11823 0.1010197 10 0
## [6,] 106 1511 945 427 12391 0.1005812 13 0
require(stats); require(graphics)
## work with pre-seatbelt period to identify a model, use logs
work <- window(log10(UKDriverDeaths), end = 1982+11/12)
par(mfrow = c(3, 1))
plot(work); acf(work); pacf(work)
par(mfrow = c(1, 1))
(fit <- arima(work, c(1, 0, 0), seasonal = list(order = c(1, 0, 0))))
##
## Call:
## arima(x = work, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0)))
##
## Coefficients:
## ar1 sar1 intercept
## 0.4378 0.6281 3.2274
## s.e. 0.0764 0.0637 0.0131
##
## sigma^2 estimated as 0.00157: log likelihood = 300.85, aic = -593.7
z <- predict(fit, n.ahead = 24)
ts.plot(log10(UKDriverDeaths), z$pred, z$pred+2*z$se, z$pred-2*z$se,
lty = c(1, 3, 2, 2), col = c("black", "red", "blue", "blue"))
## now see the effect of the explanatory variables
X <- Seatbelts[, c("kms", "PetrolPrice", "law")]
X[, 1] <- log10(X[, 1]) - 4
arima(log10(Seatbelts[, "drivers"]), c(1, 0, 0),
seasonal = list(order = c(1, 0, 0)), xreg = X)
##
## Call:
## arima(x = log10(Seatbelts[, "drivers"]), order = c(1, 0, 0), seasonal = list(order = c(1,
## 0, 0)), xreg = X)
##
## Coefficients:
## ar1 sar1 intercept kms PetrolPrice law
## 0.3348 0.6672 3.3539 0.0082 -1.2224 -0.0963
## s.e. 0.0775 0.0612 0.0441 0.0902 0.3839 0.0166
##
## sigma^2 estimated as 0.001476: log likelihood = 349.73, aic = -685.46