A pre-loaded example dataset in R

Main page: https://www.rdocumentation.org/packages/datasets/versions/3.6.2/topics/UKDriverDeaths

UKDriverDeaths
##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 1969 1687 1508 1507 1385 1632 1511 1559 1630 1579 1653 2152 2148
## 1970 1752 1765 1717 1558 1575 1520 1805 1800 1719 2008 2242 2478
## 1971 2030 1655 1693 1623 1805 1746 1795 1926 1619 1992 2233 2192
## 1972 2080 1768 1835 1569 1976 1853 1965 1689 1778 1976 2397 2654
## 1973 2097 1963 1677 1941 2003 1813 2012 1912 2084 2080 2118 2150
## 1974 1608 1503 1548 1382 1731 1798 1779 1887 2004 2077 2092 2051
## 1975 1577 1356 1652 1382 1519 1421 1442 1543 1656 1561 1905 2199
## 1976 1473 1655 1407 1395 1530 1309 1526 1327 1627 1748 1958 2274
## 1977 1648 1401 1411 1403 1394 1520 1528 1643 1515 1685 2000 2215
## 1978 1956 1462 1563 1459 1446 1622 1657 1638 1643 1683 2050 2262
## 1979 1813 1445 1762 1461 1556 1431 1427 1554 1645 1653 2016 2207
## 1980 1665 1361 1506 1360 1453 1522 1460 1552 1548 1827 1737 1941
## 1981 1474 1458 1542 1404 1522 1385 1641 1510 1681 1938 1868 1726
## 1982 1456 1445 1456 1365 1487 1558 1488 1684 1594 1850 1998 2079
## 1983 1494 1057 1218 1168 1236 1076 1174 1139 1427 1487 1483 1513
## 1984 1357 1165 1282 1110 1297 1185 1222 1284 1444 1575 1737 1763
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