In this chapter, we consider stationary models that may be suitable for residual series that contain no obvious trends or seasonal cycles. The fitted stationary models may then be combined with the fitted regression model to improve forecasts.
Book code
rho <- function(k, beta) {
q <- length(beta) - 1
if (k > q) ACF <- 0 else {
s1 <- 0; s2 <- 0
for (i in 1:(q-k+1)) s1 <- s1 + beta[i] * beta[i+k]
for (i in 1:(q+1)) s2 <- s2 + beta[i]^2
ACF <- s1 / s2}
ACF
}
beta <- c(1, 0.7, 0.5, 0.2)
rho.k <- rep(1, 10)
for (k in 1:10) rho.k[k] <- rho(k, beta)
plot(0:10, c(1, rho.k), pch = 4, ylab = expression(rho[k]))
abline(0, 0)
Tidyverts code
pacman::p_load("tsibble", "fable", "feasts",
"tsibbledata", "fable.prophet", "tidyverse",
"patchwork", "slider", "urca")
rho <- function(k, beta) {
q <- length(beta) - 1
if (k > q) ACF <- 0 else {
s1 <- 0; s2 <- 0
for (i in 1:(q-k+1)) s1 <- s1 + beta[i] * beta[i+k]
for (i in 1:(q+1)) s2 <- s2 + beta[i]^2
ACF <- s1 / s2}
ACF
}
dat <- tibble(
order = 0:10,
betas = list(c(1, 0.7, 0.5, 0.2)),
betas2 = list(c(1, -0.7, 0.5, -0.2)),
rho.k = map2_dbl(order, betas, ~rho(.x, .y)),
rho.k2 = map2_dbl(order, betas2, ~rho(.x, .y)))
dat |>
ggplot(aes(x = order, y = rho.k)) +
geom_hline(yintercept = 0, color = "darkgrey") +
geom_point() +
labs(y = expression(rho[k]), x = "lag k") +
theme_bw()
beta2 <- c(1, -0.7, 0.5, -0.2)
rho.k <- rep(1, 10)
for (k in 1:10) rho.k[k] <- rho(k, beta2)
plot(0:10, c(1, rho.k), pch = 4, ylab = expression(rho[k]))
abline(0, 0)
dat |>
ggplot(aes(x = order, y = rho.k2)) +
geom_hline(yintercept = 0, color = "darkgrey") +
geom_point() +
labs(y = expression(rho[k]), x = "lag k") +
theme_bw()
set.seed(1)
b <- c(0.8, 0.6, 0.4)
x <- w <- rnorm(1000)
for (t in 4:1000) {
for (j in 1:3) x[t] <- x[t] + b[j] * w[t - j]
# 0.9267779, 1.036964, 0.7863824
}
plot(x, type = "l")
set.seed(1)
dat <- tibble(
w = rnorm(1000),
betas = list(c(0.8, 0.6, 0.4))) |>
mutate(
w_lag = slide(w, ~.x, .before = 3, .after = -1),
w_lag = map(w_lag, ~rev(.x)),
t = 1:n()) |>
slice(-c(1:3)) |>
mutate(
lag_betas = map2_dbl(
w_lag,
betas,
\(.x, .y) sum(.x *.y)),
x = w + lag_betas) |>
tsibble::as_tsibble(index = t)
autoplot(dat, .var = x)
acf(x)
ACF(dat, var = x) |>
autoplot()
rho <- function(k, beta) {
q <- length(beta) - 1
if (k > q) ACF <- 0 else {
s1 <- 0; s2 <- 0
for (i in 1:(q-k+1)) s1 <- s1 + beta[i] * beta[i+k]
for (i in 1:(q+1)) s2 <- s2 + beta[i]^2
ACF <- s1 / s2}
ACF
}
beta <- c(1, 0.7, 0.5, 0.2)
rho.k <- rep(1, 10)
for (k in 1:10) rho.k[k] <- rho(k, beta)
plot(0:10, c(1, rho.k), pch = 4, ylab = expression(rho[k]))
abline(0, 0)
beta2 <- c(-0.7, 0.5, -0.2)
rho.k <- rep(1, 10)
for (k in 1:10) rho.k[k] <- rho(k, beta2)
plot(0:10, c(1, rho.k), pch = 4, ylab = expression(rho[k]))
abline(0, 0)
set.seed(1)
b <- c(0.8, 0.6, 0.4)
x <- w <- rnorm(1000)
for (t in 4:1000) {
for (j in 1:3) x[t] <- x[t] + b[j] * w[t - j]
}
plot(x, type = "l")
acf(x)pacman::p_load("tsibble", "fable", "feasts",
"tsibbledata", "fable.prophet", "tidyverse",
"patchwork", "slider", "urca")
rho <- function(k, beta) {
q <- length(beta) - 1
if (k > q) ACF <- 0 else {
s1 <- 0; s2 <- 0
for (i in 1:(q-k+1)) s1 <- s1 + beta[i] * beta[i+k]
for (i in 1:(q+1)) s2 <- s2 + beta[i]^2
ACF <- s1 / s2}
ACF
}
dat <- tibble(
order = 0:10,
betas = list(c(1, 0.7, 0.5, 0.2)),
betas2 = list(c(1, -0.7, 0.5, -0.2)),
rho.k = map2_dbl(order, betas, ~rho(.x, .y)),
rho.k2 = map2_dbl(order, betas2, ~rho(.x, .y)))
dat |>
ggplot(aes(x = order, y = rho.k)) +
geom_hline(yintercept = 0, color = "darkgrey") +
geom_point() +
labs(y = expression(rho[k]), x = "lag k") +
theme_bw()
dat |>
ggplot(aes(x = order, y = rho.k2)) +
geom_hline(yintercept = 0, color = "darkgrey") +
geom_point() +
labs(y = expression(rho[k]), x = "lag k") +
theme_bw()
set.seed(1)
dat <- tibble(
w = rnorm(1000),
betas = list(c(0.8, 0.6, 0.4))) |>
mutate(
w_lag = slide(w, ~.x, .before = 3, .after = -1),
w_lag = map(w_lag, ~rev(.x)),
t = 1:n()) |>
slice(-c(1:3)) |>
mutate(
lag_betas = map2_dbl(
w_lag,
betas,
\(.x, .y) sum(.x *.y)),
x = w + lag_betas) |>
tsibble::as_tsibble(index = t)
autoplot(dat, .var = x)
ACF(dat, var = x) |>
autoplot()Book code
x.ma <- arima(x, order = c(0, 0, 3))
x.ma
Call:
arima(x = x, order = c(0, 0, 3))
Coefficients:
ma1 ma2 ma3 intercept
0.7898 0.5665 0.3959 -0.0322
s.e. 0.0307 0.0351 0.0320 0.0898
sigma^2 estimated as 1.068: log likelihood = -1452.41, aic = 2914.83Tidyverts code
dat |>
model(arima = ARIMA(x ~ 1 + pdq(0, 0, 3) + PDQ(0, 0, 0))) |>
tidy()# A tibble: 4 × 6
.model term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 arima ma1 0.790 0.0307 25.7 1.91e-112
2 arima ma2 0.569 0.0351 16.2 1.62e- 52
3 arima ma3 0.394 0.0320 12.3 1.57e- 32
4 arima constant -0.0296 0.0900 -0.329 7.42e- 1Book code
www <- "data/pounds_nz.dat"
x <- read.table(www, header = T)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ma
Call:
arima(x = x.ts, order = c(0, 0, 1))
Coefficients:
ma1 intercept
1.000 2.8329
s.e. 0.072 0.0646
sigma^2 estimated as 0.04172: log likelihood = 4.76, aic = -3.53Tidyverts code
dat <- read_table("data/pounds_nz.dat") |>
mutate(
date = seq(
ymd("1991-01-01"),
by = "3 months",
length.out = n()) |>
yearquarter()) |>
as_tsibble()
dat_ma <- dat |>
model(arima = ARIMA(xrate ~ 1 + pdq(0,0,1) + PDQ(0, 0, 0)))
tidy(dat_ma)# A tibble: 2 × 6
.model term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 arima ma1 1.00 0.0720 13.9 1.12e-16
2 arima constant 2.83 0.0646 43.9 8.38e-35acf(x.ma$res[-1])
dat_ma |>
residuals() |>
ACF() |>
autoplot()
Book code
set.seed(1)
x <- arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))
coef(arima(x, order = c(1, 0, 1))) ar1 ma1 intercept
-0.596966371 0.502703368 -0.006571345 Tidyverts code
set.seed(1)
dat <- tibble(
x = arima.sim(
n = 10000,
list(ar = -0.6, ma = 0.5))) |>
mutate(index = 1:n()) |>
as_tsibble(index = index)
dat |>
model(arima = ARIMA(x ~ 1 + pdq(1,0,1) + PDQ(0, 0, 0))) |>
tidy()# A tibble: 3 × 6
.model term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 arima ar1 -0.597 0.0494 -12.1 2.20e-33
2 arima ma1 0.503 0.0530 9.49 2.87e-21
3 arima constant -0.0105 0.0152 -0.690 4.90e- 1Book code
www <- "data/pounds_nz.dat"
x <- read.table(www, header = T)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ar <- arima(x.ts, order = c(1, 0, 0))
x.arma <- arima(x.ts, order = c(1, 0, 1))
c("ma" = AIC(x.ma), "ar" = AIC(x.ar), "arma" = AIC(x.arma)) ma ar arma
-3.526895 -37.404165 -42.273571 Tidyverts code
dat <- read_table("data/pounds_nz.dat") |>
mutate(
date = seq(
ymd("1991-01-01"),
by = "3 months",
length.out = n()) |>
yearquarter()) |>
as_tsibble()
dat_fit <- dat |>
model(
ma = ARIMA(x ~ 1 + pdq(0,0,1) + PDQ(0, 0, 0)),
ar = ARIMA(x ~ 1 + pdq(1,0,0) + PDQ(0, 0, 0)),
arma = ARIMA(x ~ 1 + pdq(1,0,1) + PDQ(0, 0, 0)))
glance(dat_fit)# A tibble: 3 × 8
.model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
1 ma 0.0440 4.76 -3.53 -2.84 1.46 <cpl [0]> <cpl [1]>
2 ar 0.0192 21.7 -37.4 -36.7 -32.4 <cpl [1]> <cpl [0]>
3 arma 0.0163 25.1 -42.3 -41.1 -35.6 <cpl [1]> <cpl [1]>x.arma
Call:
arima(x = x.ts, order = c(1, 0, 1))
Coefficients:
ar1 ma1 intercept
0.8925 0.5319 2.9597
s.e. 0.0759 0.2021 0.2435
sigma^2 estimated as 0.01505: log likelihood = 25.14, aic = -42.27dat_fit |>
select(arma) |>
tidy()# A tibble: 3 × 6
.model term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 arma ar1 0.892 0.0759 11.8 2.19e-14
2 arma ma1 0.532 0.202 2.63 1.21e- 2
3 arma constant 0.318 0.0262 12.2 7.78e-15acf(resid(x.arma))
dat_fit |>
select(arma) |>
residuals() |>
ACF() |>
autoplot()
www <- "data/pounds_nz.dat"
x <- read.table(www, header = T)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ar <- arima(x.ts, order = c(1, 0, 0))
x.arma <- arima(x.ts, order = c(1, 0, 1))
c("ma" = AIC(x.ma), "ar" = AIC(x.ar), "arma" = AIC(x.arma))
x.arma
acf(resid(x.arma))dat <- read_table("data/pounds_nz.dat") |>
mutate(
date = seq(
ymd("1991-01-01"),
by = "3 months",
length.out = n()) |>
yearquarter()) |>
as_tsibble()
dat_fit <- dat |>
model(
ma = ARIMA(x ~ 1 + pdq(0,0,1) + PDQ(0, 0, 0)),
ar = ARIMA(x ~ 1 + pdq(1,0,0) + PDQ(0, 0, 0)),
arma = ARIMA(x ~ 1 + pdq(1,0,1) + PDQ(0, 0, 0)))
glance(dat_fit)
dat_fit |>
select(arma) |>
tidy()
dat_fit |>
select(arma) |>
residuals() |>
ACF() |>
autoplot()Book code
www <- "data/cbe.dat"
CBE <- read.table(www, header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
Time <- 1:length(Elec.ts)
Imth <- cycle(Elec.ts)
Elec.lm <- lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth))
acf(resid(Elec.lm))
Tidyverts code
cbe_ts <- read_table("data/cbe.dat") |>
select(elec) |>
mutate(
log_elec = log(elec),
date = seq(ymd("1958-01-01"),
ymd("1958-01-01") + months(n()-1),
by = "1 months"),
month = tsibble::yearmonth(date),
time = 1:n(),
imth = month(date)) |>
as_tsibble(index = month)
elec_lm <- cbe_ts |>
model("lm" = TSLM(log_elec ~ time + I(time^2) + factor(imth)))
elec_lm |>
residuals() |>
ACF() |>
autoplot()
best.order <- c(0, 0, 0)
aic.fits <- NULL
best.aic <- Inf
for (i in 0:2) for (j in 0:2) {
fit.aic <- AIC(arima(resid(Elec.lm), order = c(i, 0, j)))
if (fit.aic < best.aic) {
best.order <- c(i, 0, j)
best.arma <- arima(resid(Elec.lm), order = best.order)
best.aic <- fit.aic
}
names(fit.aic) <- paste0(i, 0, j)
aic.fits <- c(aic.fits, fit.aic)
}
# best.order of 2, 0, 0 is the minimum
aic.fits 000 001 002 100 101 102 200 201
-1617.759 -1766.314 -1825.612 -1888.537 -1913.830 -1913.173 -1913.910 -1912.896
202
-1911.186 model_resid <- elec_lm |>
residuals() |>
select(-.model) |>
model(
auto = ARIMA(.resid ~ 1 + pdq(0:2,0,0) + PDQ(0, 0, 0)),
a000 = ARIMA(.resid ~ 1 + pdq(0,0,0) + PDQ(0, 0, 0)),
a001 = ARIMA(.resid ~ 1 + pdq(0,0,1) + PDQ(0, 0, 0)),
a002 = ARIMA(.resid ~ 1 + pdq(0,0,2) + PDQ(0, 0, 0)),
a100 = ARIMA(.resid ~ 1 + pdq(1,0,0) + PDQ(0, 0, 0)),
a101 = ARIMA(.resid ~ 1 + pdq(1,0,1) + PDQ(0, 0, 0)),
a102 = ARIMA(.resid ~ 1 + pdq(1,0,2) + PDQ(0, 0, 0)),
a200 = ARIMA(.resid ~ 1 + pdq(2,0,0) + PDQ(0, 0, 0)),
a201 = ARIMA(.resid ~ 1 + pdq(2,0,1) + PDQ(0, 0, 0)),
a202 = ARIMA(.resid ~ 1 + pdq(2,0,2) + PDQ(0, 0, 0)))
# best order is 2, 0, 0
# auto would pick it without running each
model_best <- select(model_resid, "a200")
model_resid |>
glance() |>
filter(AIC == min(AIC))# A tibble: 2 × 8
.model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
1 auto 0.000459 961. -1914. -1914. -1898. <cpl [2]> <cpl [0]>
2 a200 0.000459 961. -1914. -1914. -1898. <cpl [2]> <cpl [0]>acf(resid(best.arma))
model_best |>
residuals() |>
ACF() |>
autoplot()
new.time <- seq(length(Elec.ts), length = 36)
new.data <- data.frame(Time = new.time, Imth = rep(1:12, 3))
predict.lm <- predict(Elec.lm, new.data)
predict.arma <- predict(best.arma, n.ahead = 36)
elec.pred <- ts(exp(predict.lm + predict.arma$pred), start = 1991, freq = 12)
ts.plot(cbind(Elec.ts, elec.pred), lty = 1:2)
new_data <- tibble(
date = seq(
max(cbe_ts$date) + months(1),
max(cbe_ts$date) + months(36),
by = "1 months"),
month = tsibble::yearmonth(date),
time = seq(nrow(cbe_ts), length = 36),
imth = rep(1:12, 3)) |>
as_tsibble(index = month)
plot_dat <- new_data |>
mutate(
residuals = forecast(
model_best,
h = "3 years") |>
pull(.mean),
expected_log = forecast(
elec_lm,
new_data = new_data) |>
pull(.mean),
expected = exp(expected_log + residuals)
)
ggplot() +
geom_line(data = cbe_ts, aes(x = date, y = elec)) +
geom_line(data = plot_dat, aes(x = date, y = expected), linetype = 2)
www <- "data/cbe.dat"
CBE <- read.table(www, header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
Time <- 1:length(Elec.ts)
Imth <- cycle(Elec.ts)
Elec.lm <- lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth))
acf(resid(Elec.lm))
best.order <- c(0, 0, 0)
best.aic <- Inf
for (i in 0:2) for (j in 0:2) {
fit.aic <- AIC(arima(resid(Elec.lm), order = c(i, 0,
j)))
if (fit.aic < best.aic) {
best.order <- c(i, 0, j)
best.arma <- arima(resid(Elec.lm), order = best.order)
best.aic <- fit.aic
}
best.order
acf(resid(best.arma))
new.time <- seq(length(Elec.ts), length = 36)
new.data <- data.frame(Time = new.time, Imth = rep(1:12, 3))
predict.lm <- predict(Elec.lm, new.data)
predict.arma <- predict(best.arma, n.ahead = 36)
elec.pred <- ts(exp(predict.lm + predict.arma$pred), start = 1991, freq = 12)
ts.plot(cbind(Elec.ts, elec.pred), lty = 1:2)cbe_ts <- read_table("data/cbe.dat") |>
select(elec) |>
mutate(
date = seq(ymd("1958-01-01"),
ymd("1958-01-01") + months(n()-1),
by = "1 months"),
month = tsibble::yearmonth(date),
time = 1:n(),
imth = month(date)) |>
as_tsibble(index = month)
elec_lm <- cbe_ts |>
model("lm" = TSLM(log(elec) ~ time + I(time^2) + factor(imth)))
elec_lm |>
residuals() |>
ACF() |>
autoplot()
model_resid <- elec_lm |>
residuals() |>
select(-.model) |>
model(
auto = ARIMA(x ~ 1 + pdq(0:2,0,0) + PDQ(0, 0, 0)),
a000 = ARIMA(x ~ 1 + pdq(0,0,0) + PDQ(0, 0, 0)),
a001 = ARIMA(x ~ 1 + pdq(0,0,1) + PDQ(0, 0, 0)),
a002 = ARIMA(x ~ 1 + pdq(0,0,2) + PDQ(0, 0, 0)),
a100 = ARIMA(x ~ 1 + pdq(1,0,0) + PDQ(0, 0, 0)),
a101 = ARIMA(x ~ 1 + pdq(1,0,1) + PDQ(0, 0, 0)),
a102 = ARIMA(x ~ 1 + pdq(1,0,2) + PDQ(0, 0, 0)),
a200 = ARIMA(x ~ 1 + pdq(2,0,0) + PDQ(0, 0, 0)),
a201 = ARIMA(x ~ 1 + pdq(2,0,1) + PDQ(0, 0, 0)),
a202 = ARIMA(x ~ 1 + pdq(2,0,2) + PDQ(0, 0, 0)))
# best order is 2, 0, 1
model_best <- select(model_resid, "a201")
model_resid |>
glance() |>
filter(AIC == min(AIC))
model_best |>
residuals() |>
ACF() |>
autoplot()
new_data <- tibble(
date = seq(
max(cbe_ts$date) + months(1),
max(cbe_ts$date) + months(36),
by = "1 months"),
month = tsibble::yearmonth(date),
time = seq(nrow(cbe_ts), length = 36),
imth = rep(1:12, 3)) |>
as_tsibble(index = month)
plot_dat <- new_data |>
mutate(
residuals = forecast(
model_best,
h = "3 years") |>
pull(.mean),
expected_log = forecast(
elec_lm,
new_data = new_data) |>
pull(.mean),
expected = exp(expected_log + residuals)
)
ggplot() +
geom_line(data = cbe_ts, aes(x = date, y = elec)) +
geom_line(data = plot_dat, aes(x = date, y = expected), linetype = 2)Book code
www <- "data/wave.dat"
wave.dat <- read.table(www, header = T)
# attach (wave.dat)
layout(1:3)
plot(as.ts(wave.dat$waveht), ylab = 'Wave height')
acf(wave.dat$waveht)
pacf(wave.dat$waveht)
Tidyverts code
wave_dat <- read_table("data/wave.dat") |>
mutate(index = 1:n()) |>
as_tsibble(index = index)
p1 <- autoplot(wave_dat) +
labs(y = 'Wave height')
p2 <- ACF(wave_dat) |>
autoplot()
p3 <- PACF(wave_dat) |>
autoplot()
p1 / p2 / p3
layout(1:3)
wave.arma <- arima(wave.dat$waveht, order = c(4,0,4))
acf(wave.arma$res[-(1:4)])
pacf (wave.arma$res[-(1:4)])
hist(wave.arma$res[-(1:4)], xlab='height / mm', main='')
wave_arma <- wave_dat |>
model(arma = ARIMA(waveht ~ 1 + pdq(4,0,4) + PDQ(0,0,0)))
p1 <- wave_arma |>
residuals() |>
ACF() |>
autoplot()
p2 <- wave_arma |>
residuals() |>
PACF() |>
autoplot()
p3 <- wave_arma |>
residuals() |>
ggplot(aes(x = .resid)) +
geom_histogram()
p1 / p2 / p3
www <- "data/wave.dat"
wave.dat <- read.table(www, header = T)
# attach (wave.dat)
layout(1:3)
plot(as.ts(wave.dat$waveht), ylab = 'Wave height')
acf(wave.dat$waveht)
pacf(wave.dat$waveht)
wave.arma <- arima(wave.dat$waveht, order = c(4,0,4))
acf(wave.arma$res[-(1:4)])
pacf (wave.arma$res[-(1:4)])
hist(wave.arma$res[-(1:4)], xlab='height / mm', main='')R commands used in examplesfeasts::PACF(): Function PACF computes an estimate of the partial autocorrelation function of a (possibly multivariate) tsibble.feasts::ACF(): The function ACF computes an estimate of the autocorrelation function of a (possibly multivariate) tsibble.purrr:map2_dbl(): The map functions transform their input by applying a function to each element of a list or atomic vector and returning an object of the same length as the input.ggplot2:geom_hline(): These geoms add reference lines (sometimes called rules) to a plotbase::expression(): If the text argument to one of the text-drawing functions in R is an expression, the argument is interpreted as a mathematical expression and the output will be formatted according to TeX-like rules.slider::slide(): iterates through .x using a sliding window, applying .f to each sub-window of .x.base::rev(): provides a reversed version of its argument.fable::ARIMA(): Searches through the model space specified in the specials to identify the best ARIMA model, with the lowest AIC, AICc or BIC value.pdq(): The pdq special is used to specify non-seasonal components of the model.PDQ(): The PDQ special is used to specify seasonal components of the model. To force a non-seasonal fit, specify PDQ(0, 0, 0) in the RHS of the model formula.