Box-Jenkins analysis of Internet use (DK §8.4)#
ARMA(p, q) models, p, q = 0, …, 5, for the first differences of the number of users logged on to an Internet server each minute (Makridakis, Wheelwright and Hyndman 1998), in state space form (DK §3.4), compared by AIC; then again with 14 observations treated as missing, and forecasts with 50% intervals from the ARMA(1, 1) model.
Results#
The log likelihoods agree with statsmodels’ SARIMAX to four decimals, with and without the missing observations; for example -254.1497 and -225.7704 for the ARMA(1, 1). Higher-order models can have several local optima, where different optimizers stop at different points.
DK’s tables use a scale for the AIC that we could not identify: for the complete series they are close to \(n^{-1}[-\log L + 2(p + q)]\), but no simple formula fits the table with missing observations, and DK report optimizer failures for several models. This page uses DK’s own formula (DK §7.4), \(n^{-1}[-2\log L + 2w]\). By it, ARMA(3, 0), ARMA(1, 1) and some higher-order models lie within 0.03 of each other, the conclusion DK draw.
Program#
!> DK 8.4: ARMA models for the first differences of the number of users
!> logged on to an Internet server each minute (99 observations), fitted in
!> state space form (DK 3.4), with the AIC of DK 7.4,
!>
!> AIC = n^-1 [-2 log L + 2 w], w = p + q + 1,
!>
!> for p, q = 0..5 (Table 8.1); then again with 14 observations treated as
!> missing (Table 8.2); then forecasts from the ARMA(1, 1) model with 50%
!> intervals (Fig. 8.9).
!>
!> The log likelihoods agree with statsmodels' SARIMAX to 4 decimals, with
!> and without the missing observations (e.g. ARMA(1, 1): -254.1497 and
!> -225.7704). Higher-order models can have several local optima.
!> DK's tables use another scale: they are close to n^-1 [-log L + 2 (p + q)]
!> for the complete series, but no simple formula fits the table with missing
!> observations, and DK report optimizer failures for several models. With
!> the formula above, AIC chooses among the same few models: ARMA(3, 0),
!> ARMA(1, 1) and some higher-order models are within 0.03 of each other.
!>
!> Run from the repo root, after `.venv/bin/python data/fetch_dk_data.py`:
!> fpm run --example dk_8_4_internet
program dk_8_4_internet
use statespace
use csv_io, only: read_csv, column
use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
implicit none
integer, parameter :: missing(14) = [6, 16, 26, 36, 46, 56, 66, 72, 73, 74, 75, 76, &
86, 96]
character(len=32), allocatable :: names(:)
real(dp), allocatable :: data(:, :), users(:), dy(:, :)
call read_csv("data/wwwusage.csv", names, data)
users = column(names, data, "value")
allocate (dy(1, size(users) - 1))
dy(1, :) = users(2:) - users(:size(users) - 1)
print '(a)', "Table 8.1: AIC for ARMA(p, q), differenced series"
call aic_table(dy)
dy(1, missing) = ieee_value(1.0_dp, ieee_quiet_nan)
print '(/, a)', "Table 8.2: the same with 14 observations missing"
call aic_table(dy)
print '(/, a)', &
"ARMA(1, 1) forecasts with 50% intervals, series with missing observations"
call forecasts(dy, 10)
contains
function arma(y, p, q) result(model)
real(dp), intent(in) :: y(:, :)
integer, intent(in) :: p, q
type(structural_model_t) :: model
type(component_holder_t) :: comps(1)
integer :: info
comps(1)%c = arima_t(ar=p, ma=q)
model = structural_model(y, comps, info)
if (info /= SS_OK) error stop "model"
end function arma
subroutine aic_table(y)
real(dp), intent(in) :: y(:, :)
type(structural_model_t) :: model
type(fit_result_t) :: res
type(fit_options_t) :: opts
character(len=10) :: cell(0:5)
integer :: p, q, info
opts%compute_cov = .false.
opts%factr = 10.0_dp
opts%pgtol = 1.0e-8_dp
print '(4x, 6(i10))', [(q, q=0, 5)]
do p = 0, 5
cell = ""
do q = 0, 5
if (p == 0 .and. q == 0) cycle
model = arma(y, p, q)
call fit(model, res, options=opts, info=info)
if (info /= SS_OK) then
cell(q) = "failed"
else
write (cell(q), '(f10.3)') (-2 * res%llf + 2 * (p + q + 1)) / size(y, 2)
end if
end do
print '(i4, 6a10)', p, cell
end do
end subroutine aic_table
subroutine forecasts(y, h)
real(dp), intent(in) :: y(:, :)
integer, intent(in) :: h
type(structural_model_t) :: model
type(fit_result_t) :: res
type(filter_result_t) :: fres
type(forecast_result_t) :: fc
real(dp), parameter :: z50 = 0.6744897501960817_dp ! 75% normal quantile
real(dp) :: sd
integer :: j, info
model = arma(y, 1, 1)
call fit(model, res, info=info)
if (info /= SS_OK) error stop "fit"
call model%filter(res%params, fres, info)
call forecast(model%rep, fres, h, fc, info)
if (info /= SS_OK) error stop "forecast"
print '(4x, a, 3a10)', "t", "forecast", "lower", "upper"
do j = 1, h
sd = sqrt(fc%cov(1, 1, j))
print '(i5, 3f10.3)', size(y, 2) + j, fc%mean(1, j), fc%mean(1, j) - z50 * sd, &
fc%mean(1, j) + z50 * sd
end do
end subroutine forecasts
end program dk_8_4_internet
Output#
Table 8.1: AIC for ARMA(p, q), differenced series
0 1 2 3 4 5
0 5.554 5.251 5.255 5.246 5.241
1 5.346 5.195 5.215 5.198 5.203 5.215
2 5.275 5.215 5.226 5.210 5.184 5.199
3 5.172 5.191 5.208 5.224 5.200 5.220
4 5.191 5.210 5.230 5.243 5.200 5.234
5 5.211 5.229 5.187 5.206 5.163 5.179
Table 8.2: the same with 14 observations missing
0 1 2 3 4 5
0 4.939 4.687 4.686 4.673 4.683
1 4.704 4.622 4.640 4.633 4.638 4.657
2 4.671 4.641 4.651 4.641 4.639 4.651
3 4.605 4.625 4.643 4.661 4.630 4.679
4 4.625 4.642 4.659 4.677 4.640 4.656
5 4.641 4.659 4.652 4.694 4.657 4.676
ARMA(1, 1) forecasts with 50% intervals, series with missing observations
t forecast lower upper
100 -0.714 -2.890 1.462
101 -0.469 -3.766 2.829
102 -0.308 -3.984 3.369
103 -0.202 -4.030 3.627
104 -0.132 -4.024 3.759
105 -0.087 -4.006 3.832
106 -0.057 -3.988 3.873
107 -0.037 -3.973 3.898
108 -0.025 -3.962 3.913
109 -0.016 -3.955 3.922