Chapter 1: Time Series Data

In most branches of science, engineering, and commerce, there are variables measured sequentially in time. Reserve banks record interest rates and exchange rates each day. The government statistics department will compute the country’s gross domestic product on a yearly basis. Newspapers publish yesterday’s noon temperatures for capital cities from around the world. Meteorological offices record rainfall at many different sites with differing resolutions. When a variable is measured sequentially in time over or at a fixed interval, known as the sampling interval, the resulting data form a time series.

We can use varied methods on the AP R object which is of the class ts which is an abbreviation for matrix of a ‘time series’ format. It is not a data.frame or a tibble as you have seen in your previous R courses. The base plot() function has a plotting method for ts data objects. Within the tidyverts we can use autoplot() to create default plots on our tsibble objects.

1.4.1 ‘A flying start: Air passenger bookings’ Modern Look

Book code

data(AirPassengers)
AP <- AirPassengers
AP
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 178 163 172 178 199 199 184 162 146 166
1952 171 180 193 181 183 218 230 242 209 191 172 194
1953 196 196 236 235 229 243 264 272 237 211 180 201
1954 204 188 235 227 234 264 302 293 259 229 203 229
1955 242 233 267 269 270 315 364 347 312 274 237 278
1956 284 277 317 313 318 374 413 405 355 306 271 306
1957 315 301 356 348 355 422 465 467 404 347 305 336
1958 340 318 362 348 363 435 491 505 404 359 310 337
1959 360 342 406 396 420 472 548 559 463 407 362 405
1960 417 391 419 461 472 535 622 606 508 461 390 432

Tidyverts code

# install.packages("pacman")
pacman::p_load("tsibble", "fable", "feasts",
    "tsibbledata", "fable.prophet", "tidyverse", "patchwork")

ap <- as_tsibble(AirPassengers)  |>
    mutate(
        Year = year(index),
        Month = month(index)
        )
ap
# A tsibble: 144 x 4 [1M]
      index value  Year Month
      <mth> <dbl> <dbl> <dbl>
 1 1949 Jan   112  1949     1
 2 1949 Feb   118  1949     2
 3 1949 Mar   132  1949     3
 4 1949 Apr   129  1949     4
 5 1949 May   121  1949     5
 6 1949 Jun   135  1949     6
 7 1949 Jul   148  1949     7
 8 1949 Aug   148  1949     8
 9 1949 Sep   136  1949     9
10 1949 Oct   119  1949    10
# ℹ 134 more rows
class(AP)
[1] "ts"
class(ap)
[1] "tbl_ts"     "tbl_df"     "tbl"        "data.frame"
start(AP)
[1] 1949    1
head(ap, 1)
# A tsibble: 1 x 4 [1M]
     index value  Year Month
     <mth> <dbl> <dbl> <dbl>
1 1949 Jan   112  1949     1
end(AP)
[1] 1960   12
tail(ap, 1)
# A tsibble: 1 x 4 [1M]
     index value  Year Month
     <mth> <dbl> <dbl> <dbl>
1 1960 Dec   432  1960    12
frequency(AP)
[1] 12
guess_frequency(ap$index)
[1] 12
plot(AP, ylab = "Passengers (1000's)",
    main = "Monthly totals of international airline passengers, 1949 to 1960")

autoplot(ap) +
    labs(
        y = "Passengers (1000's)",
        title = "Monthly totals of international airline passengers, 1949 to 1960"
    )

layout(1:2)
plot(aggregate(AP))
boxplot(AP ~ cycle(AP))

# pacman::p_load("patchwork")
psum <- autoplot(summarise(index_by(ap, Year), value = sum(value)))
pbox <- ggplot(ap, aes(x = factor(Month), y = value)) +
    geom_boxplot()
psum / pbox

Show the book code in full 
data(AirPassengers)
AP <- AirPassengers
AP
class(AP)
start(AP)
end(AP)
frequency(AP)
plot(AP, ylab = "Passengers (1000's)",
    main = "Monthly totals of international airline passengers, 1949 to 1960")
layout(1:2)
plot(aggregate(AP))
boxplot(AP ~ cycle(AP))
Show the Tidyverts code in full 
data(AirPassengers)
install.packages("pacman")
pacman::p_load("tsibble", "fable", "feasts",
    "tsibbledata", "fable.prophet", "tidyverse", "patchwork")
ap <- as_tsibble(AirPassengers)  |>
    mutate(
        Year = year(index),
        Month = month(index)
        )
ap
class(ap)
head(ap, 1)
tail(ap, 1)
guess_frequency(ap$index)
autoplot(ap) +
    labs(
        y = "Passengers (1000's)",
        title = "Monthly totals of international airline passengers, 1949 to 1960"
    )
psum <- autoplot(summarise(index_by(ap, Year), value = sum(value)))
pbox <- ggplot(ap, aes(x = factor(Month), y = value)) +
    geom_boxplot()
psum / pbox

1.4.2 ‘Unemployment: Maine’ Modern Look

Note: The book uses attach() in their examples, a terrible practice (tp) in modern R programming. We have opted not to recreate their code exactly when using attach().

Book code

 Maine.month <- read.table("data/Maine.dat", header = TRUE)
 # attach(Maine.month) tp
 Maine.month.ts <- ts(Maine.month, start = c(1996, 1), freq = 12)
 Maine.annual.ts <- aggregate(Maine.month.ts)/12
 Maine.annual.ts
Time Series:
Start = 1996 
End = 2005 
Frequency = 1 
      unemploy
 [1,] 5.258333
 [2,] 5.125000
 [3,] 4.508333
 [4,] 3.950000
 [5,] 3.275000
 [6,] 3.733333
 [7,] 4.341667
 [8,] 4.991667
 [9,] 4.616667
[10,] 4.841667

Tidyverts code

# install.packages("pacman")
pacman::p_load("tsibble", "fable", 
    "feasts", "tsibbledata",
    "fable.prophet", "tidyverse",
    "patchwork")

maine_month <- read_table("data/Maine.dat")
start_date <- lubridate::ymd("1996-01-01")
date_seq <- seq(start_date,
    start_date + months(nrow(maine_month)-1),
    by = "1 months")
maine <- tibble(
    dates = date_seq,
    year = lubridate::year(date_seq),
    month = lubridate::month(date_seq),
    value = pull(maine_month, unemploy)
)
maine_month_ts <- maine |>
    mutate(index = tsibble::yearmonth(dates)) |>
    as_tsibble(index = index)

maine_annual_ts <- summarise(index_by(maine_month_ts, year), value = mean(value))
maine_annual_ts 
# A tsibble: 11 x 2 [1Y]
    year value
   <dbl> <dbl>
 1  1996  5.26
 2  1997  5.12
 3  1998  4.51
 4  1999  3.95
 5  2000  3.28
 6  2001  3.73
 7  2002  4.34
 8  2003  4.99
 9  2004  4.62
10  2005  4.84
11  2006  4.66
layout(1:2)
plot(Maine.month.ts, ylab = "unemployed (%)")
plot(Maine.annual.ts, ylab = "unemployed (%)")

mp <- autoplot(maine_month_ts, .vars = value) +
    labs(y = "unemployed (%)")
yp <- autoplot(maine_annual_ts) +
    labs(y = "unemployed (%)") +
    scale_x_continuous(breaks = seq(1900, 2010, by = 2))
mp / yp

Maine.Feb <- window(Maine.month.ts, start = c(1996,2), freq = TRUE)
Maine.Aug <- window(Maine.month.ts, start = c(1996,8), freq = TRUE)
(Feb.ratio <- mean(Maine.Feb) / mean(Maine.month.ts))
[1] 1.222529
(Aug.ratio <- mean(Maine.Aug) / mean(Maine.month.ts))
[1] 0.8163732
maine_feb = filter(maine, month == 2)
maine_aug = filter(maine, month == 8)
(feb_ratio <- mean(pull(maine_feb, value)) / mean(pull(maine_month_ts, value)))
[1] 1.222529
(aug_ratio <- mean(pull(maine_aug, value)) / mean(pull(maine_month_ts, value)))
[1] 0.8163732
Show the Book code in full 
Maine.month <- read.table("data/Maine.dat", header = TRUE)
# attach(Maine.month) this is a bad deal
class(Maine.month)
Maine.month.ts <- ts(Maine.month, start = c(1996, 1), freq = 12)
Maine.annual.ts <- aggregate(Maine.month.ts)/12
layout(1:2)
plot(Maine.month.ts, ylab = "unemployed (%)")
plot(Maine.annual.ts, ylab = "unemployed (%)")
Maine.Feb <- window(Maine.month.ts, start = c(1996,2), freq = TRUE)
Maine.Aug <- window(Maine.month.ts, start = c(1996,8), freq = TRUE)
(Feb.ratio <- mean(Maine.Feb) / mean(Maine.month.ts))
(Aug.ratio <- mean(Maine.Aug) / mean(Maine.month.ts))
Show the Tidyverts code in full 
install.packages("pacman")
pacman::p_load("tsibble", "fable", "feasts", "tsibbledata",
    "fable.prophet", "tidyverse", "patchwork")

maine_month <- read_table("data/Maine.dat")
start_date <- lubridate::ymd("1996-01-01")
date_seq <- seq(start_date,
    start_date + months(nrow(maine_month)-1),
    by = "1 months")
maine <- tibble(
    dates = date_seq,
    year = lubridate::year(date_seq),
    month = lubridate::month(date_seq),
    value = pull(maine_month, unemploy)
)
maine_month_ts <- maine |>
    mutate(index = tsibble::yearmonth(dates)) |>
    as_tsibble(index = index)

maine_annual_ts <- summarise(index_by(maine_month_ts, year), value = mean(value))
mp <- autoplot(maine_month_ts, .vars = value) +
    labs(y = "unemployed (%)")
yp <- autoplot(maine_annual_ts) +
    labs(y = "unemployed (%)") +
    scale_x_continuous(breaks = seq(1900, 2010, by = 2))
mp / yp
maine_feb = filter(maine, month == 2)
maine_aug = filter(maine, month == 8)
(feb_ratio <- mean(pull(maine_feb, value)) / mean(pull(maine_month_ts, value)))
(aug_ratio <- mean(pull(maine_aug, value)) / mean(pull(maine_month_ts, value)))

1.4.2 ‘US unemployment’ Modern Look

Book code

US.month <- read.table("data/USunemp.dat", header = T)
# attach(US.month) tp
US.month.ts <- ts(US.month$USun, start=c(1996,1), end=c(2006,10), freq = 12)
plot(US.month.ts, ylab = "unemployed (%)")

Tidyverts code

us_month <- read_table("data/USunemp.dat")
start_date <- lubridate::ymd("1996-01-01")
date_seq <- seq(start_date,
    start_date + months(nrow(us_month)-1),
    by = "1 months")
us_month_ts <- tibble(
    dates = date_seq,
    month = tsibble::yearmonth(dates),
    value = pull(us_month, USun)) |>
    as_tsibble(index = month)
autoplot(us_month_ts) + labs(y = "unemployed (%)")

Show the Book code in full 
US.month <- read.table("data/USunemp.dat", header = T)
attach(US.month)
US.month.ts <- ts(USun, start=c(1996,1), end=c(2006,10), freq = 12)
plot(US.month.ts, ylab = "unemployed (%)")
Show the Tidyverts code in full 
us_month <- read_table("data/USunemp.dat")
start_date <- lubridate::ymd("1996-01-01")
date_seq <- seq(start_date,
    start_date + months(nrow(us_month)-1),
    by = "1 months")
us_month_ts <- tibble(
    dates = date_seq,
    month = tsibble::yearmonth(dates),
    value = pull(us_month, USun)) |>
    as_tsibble(index = month)
autoplot(us_month_ts) + labs(y = "unemployed (%)")

1.4.3 ‘electricity, beer, and chocolate’ Modern Look

Book code

CBE <- read.table("data/cbe.dat", header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
Beer.ts <- ts(CBE[, 2], start = 1958, freq = 12)
Choc.ts <- ts(CBE[, 1], start = 1958, freq = 12)
plot(cbind(Elec.ts, Beer.ts, Choc.ts))

Tidyverts code

cbe <- read_table("data/cbe.dat")
start_date <- lubridate::ymd("1958-01-01")
date_seq <- seq(start_date,
    start_date + months(nrow(cbe)-1),
    by = "1 months")
cbe_ts <- cbe |>
    mutate(
        date = date_seq,
        month = tsibble::yearmonth(date)) |>
    pivot_longer(choc:elec) |>
    mutate(name = factor(name,
        levels = c("elec", "beer", "choc"))) |>
    as_tsibble(index = month, key = name)
autoplot(cbe_ts) +
    facet_wrap(~name, ncol = 1, scales = "free_y")

data(AirPassengers)
AP <- AirPassengers
AP.elec <- ts.intersect(AP, Elec.ts)
AP <- AP.elec[,1]
Elec <- AP.elec[,2]
layout(1:2)
plot(AP, main = "", ylab = "Air passengers / 1000's")
plot(Elec, main = "", ylab = "Electricity production / MkWh")

data(AirPassengers)
ap <- as_tsibble(AirPassengers)  |>
        mutate(
            Year = year(index),
            Month = month(index)
            )
ap_elec <- cbe_ts |>
    rename(elec = value) |>
    filter(name == "elec") |>
    inner_join(
        rename(ap, airp = value),
        by = join_by(month == index))            
ap <- autoplot(ap_elec, .vars = airp) +
    labs(y = "Air passengers / 1000's")
ep <- autoplot(ap_elec, .vars = elec) +
    labs(y = "Electricity production / MkWh")
ap / ep

plot(as.vector(AP), as.vector(Elec),
    xlab = "Air passengers / 1000's",
    ylab = "Electricity production / MWh")
abline(reg = lm(Elec ~ AP))

ap_elec |>
    ggplot(aes(x = airp, y = elec)) +
    geom_point() +
    geom_smooth(method = "lm") +
    labs(
        x = "Air passengers / 1000's",
        y = "Electricity production / MWh"
    )

cor(AP, Elec)
[1] 0.8841668
cor(pull(ap_elec, airp), pull(ap_elec, elec))
[1] 0.8841668
Show the Book code in full 
CBE <- read.table("data/cbe.dat", header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
Beer.ts <- ts(CBE[, 2], start = 1958, freq = 12)
Choc.ts <- ts(CBE[, 1], start = 1958, freq = 12)
plot(cbind(Elec.ts, Beer.ts, Choc.ts))
data(AirPassengers)
AP <- AirPassengers
AP.elec <- ts.intersect(AP, Elec.ts)
AP <- AP.elec[,1]
Elec <- AP.elec[,2]
layout(1:2)
plot(AP, main = "", ylab = "Air passengers / 1000's")
plot(Elec, main = "", ylab = "Electricity production / MkWh")
plot(as.vector(AP), as.vector(Elec),
    xlab = "Air passengers / 1000's",
    ylab = "Electricity production / MWh")
abline(reg = lm(Elec ~ AP))
cor(AP, Elec)
Show the Tidyverts code in full 
cbe <- read_table("data/cbe.dat")
start_date <- lubridate::ymd("1958-01-01")
date_seq <- seq(start_date,
    start_date + months(nrow(cbe)-1),
    by = "1 months")
cbe_ts <- cbe |>
    mutate(
        date = date_seq,
        month = tsibble::yearmonth(date)) |>
    pivot_longer(choc:elec) |>
    mutate(name = factor(name,
        levels = c("elec", "beer", "choc"))) |>
    as_tsibble(index = month, key = name)
autoplot(cbe_ts) +
    facet_wrap(~name, ncol = 1, scales = "free_y")
data(AirPassengers)
ap <- as_tsibble(AirPassengers)  |>
        mutate(
            Year = year(index),
            Month = month(index)
            )
ap_elec <- cbe_ts |>
    rename(elec = value) |>
    filter(name == "elec") |>
    inner_join(
        rename(ap, airp = value),
        by = join_by(month == index))            
ap <- autoplot(ap_elec, .vars = airp) +
    labs(y = "Air passengers / 1000's")
ep <- autoplot(ap_elec, .vars = elec) +
    labs(y = "Electricity production / MkWh")
ap / ep
ap_elec |>
    ggplot(aes(x = airp, y = elec)) +
    geom_point() +
    geom_smooth(method = "lm") +
    labs(
        x = "Air passengers / 1000's",
        y = "Electricity production / MWh"
    )
cor(pull(ap_elec, airp), pull(ap_elec, elec))

1.4.4 ‘Quarterly exchnage rate’ Modern Look

Book code

Z <- read.table("data/pounds_nz.dat", header = T)
Z.ts <- ts(Z, st = 1991, fr = 4)
plot(Z.ts,
    xlab = "time / years",
    ylab = "Quarterly exchange rate in $NZ / pound")

Tidyverts code

z <- read_table("data/pounds_nz.dat")
date_seq <- seq(
    lubridate::ymd("1991-01-01"),
    by = "3 months",
    length.out = nrow(z))
z_ts <- z |>
    mutate(
        date = date_seq,
        quarter = tsibble::yearquarter(date)) |>
    as_tsibble(index = quarter)
autoplot(z_ts) +
    labs(
        x = "time / years",
        y = "Quarterly exchange rate in $NZ / pound"
    )

Z.92.96 <- window(Z.ts, start = c(1992, 1), end = c(1996, 1))
Z.96.98 <- window(Z.ts, start = c(1996, 1), end = c(1998, 1))
layout (1:2)
plot(Z.92.96,
    ylab = "Exchange rate in $NZ/pound",
    xlab = "Time (years)" )
 plot(Z.96.98,
    ylab = "Exchange rate in $NZ/pound",
    xlab = "Time (years)" )

z_ts |>
    mutate(interval = ifelse(
            date < ymd("1996/01/01"),
            "92-95",
            "96-98")) |>
    ggplot(aes(x = quarter, y = xrate)) +
    geom_line() +
    facet_wrap(~interval, ncol = 1,
        scales = "free_x")

Show the Book code in full 
Z <- read.table("data/pounds_nz.dat", header = T)
Z.ts <- ts(Z, st = 1991, fr = 4)
plot(Z.ts,
    xlab = "time / years",
    ylab = "Quarterly exchange rate in $NZ / pound")
Z.92.96 <- window(Z.ts, start = c(1992, 1), end = c(1996, 1))
Z.96.98 <- window(Z.ts, start = c(1996, 1), end = c(1998, 1))
layout (1:2)
plot(Z.92.96,
    ylab = "Exchange rate in $NZ/pound",
    xlab = "Time (years)" )
 plot(Z.96.98,
    ylab = "Exchange rate in $NZ/pound",
    xlab = "Time (years)" )
Show the Tidyverts code in full 
z <- read_table("data/pounds_nz.dat")
date_seq <- seq(
    lubridate::ymd("1991-01-01"),
    by = "3 months",
    length.out = nrow(z))
z_ts <- z |>
    mutate(
        date = date_seq,
        quarter = tsibble::yearquarter(date)) |>
    as_tsibble(index = quarter)
autoplot(z_ts) +
    labs(
        x = "time / years",
        y = "Quarterly exchange rate in $NZ / pound"
    )
z_ts |>
    mutate(interval = ifelse(
            date < ymd("1996/01/01"),
            "92-95",
            "96-98")) |>
    ggplot(aes(x = quarter, y = xrate)) +
    geom_line() +
    facet_wrap(~interval, ncol = 1,
        scales = "free_x")

1.4.5 ‘Global temperature series’ Modern Look

Book code

Global <- scan("data/global.dat")
Global.ts <- ts(Global,
    st = c(1856, 1),
    end = c(2005, 12),
    fr = 12)
Global.annual <- aggregate(Global.ts, FUN = mean)
layout (1:2)
plot(Global.ts)
plot(Global.annual)

Tidyverts code

global <- scan("data/global.dat")
date_seq <- seq(
    lubridate::ymd("1856-01-01"),
    by = "1 months",
    length.out = length(global))
global_df <- tibble(
    date = date_seq,
    value = global,
    month = tsibble::yearmonth(date_seq),
    year = lubridate::year(date_seq))

global_ts <-  as_tsibble(global_df, index = month)
global_annual_ts <-  global_df |>
    summarise(mean(value), .by = year) |>
    as_tsibble(index = year)
gm <- autoplot(global_ts)
ga <- autoplot(global_annual_ts)
gm / ga

New.series <- window(Global.ts, start=c(1970, 1), end=c(2005, 12))
New.time <- time(New.series)
plot(New.series); abline(reg=lm(New.series ~ New.time))

global_ts |>
    filter(date >= ymd("1970-01-01"), date <= ymd("2005-12-31")) |>
    ggplot(aes(x = date, y = value)) +
    geom_point(color = "lightgrey") +
    geom_line() +
    geom_smooth(method = "lm") +
    theme_bw()

Show the Book code in full 
Global <- scan("data/global.dat")
Global.ts <- ts(Global,
    st = c(1856, 1),
    end = c(2005, 12),
    fr = 12)
Global.annual <- aggregate(Global.ts, FUN = mean)
plot(Global.ts)
plot(Global.annual)
New.series <- window(Global.ts, start=c(1970, 1), end=c(2005, 12))
New.time <- time(New.series)
plot(New.series); abline(reg=lm(New.series ~ New.time))
Show the Tidyverts code in full 
global <- scan("data/global.dat")
date_seq <- seq(
    lubridate::ymd("1856-01-01"),
    by = "1 months",
    length.out = length(global))
global_df <- tibble(
    date = date_seq,
    value = global,
    month = tsibble::yearmonth(date_seq),
    year = lubridate::year(date_seq))

global_ts <-  as_tsibble(global_df, index = month)
global_annual_ts <-  global_df |>
    summarise(mean(value), .by = year) |>
    as_tsibble(index = year)
gm <- autoplot(global_ts)
ga <- autoplot(global_annual_ts)
gm / ga

global_ts |>
    filter(date >= ymd("1970-01-01"), date <= ymd("2005-12-31")) |>
    ggplot(aes(x = date, y = value)) +
    geom_point(color = "lightgrey") +
    geom_line() +
    geom_smooth(method = "lm") +
    theme_bw()

1.5.5 ‘Decomposition in R’ Modern Look

Book code

plot(decompose(Elec.ts))

Tidyverts code

elec_decompose <- filter(cbe_ts, name == "elec") |>
    model(feasts::classical_decomposition(value,
        type = "add"))  |>
    components()
autoplot(elec_decompose)

Elec.decom <- decompose(Elec.ts, type = "mult")
plot(Elec.decom)

elec_decompose <- filter(cbe_ts, name == "elec") |>
    model(feasts::classical_decomposition(value,
        type = "mult"))  |>
    components()
autoplot(elec_decompose)

Trend <- Elec.decom$trend
Seasonal <- Elec.decom$seasonal
ts.plot(cbind(Trend, Trend * Seasonal), lty = 1:2)

elec_decompose |>
    ggplot(aes(x = month)) +
    geom_line(aes(y = trend)) +
    geom_line(aes(y = trend * seasonal), linetype = 2)

Show the Book code in full 
plot(decompose(Elec.ts))
Elec.decom <- decompose(Elec.ts, type = "mult")
plot(Elec.decom)
Trend <- Elec.decom$trend
Seasonal <- Elec.decom$seasonal
ts.plot(cbind(Trend, Trend * Seasonal), lty = 1:2)
Show the Tidyverts code in full 
elec_decompose <- filter(cbe_ts, name == "elec") |>
    model(feasts::classical_decomposition(value,
        type = "add"))  |>
    components()
autoplot(elec_decompose)
elec_decompose <- filter(cbe_ts, name == "elec") |>
    model(feasts::classical_decomposition(value,
        type = "mult"))  |>
    components()
autoplot(elec_decompose)
elec_decompose |>
    ggplot(aes(x = month)) +
    geom_line(aes(y = trend)) +
    geom_line(aes(y = trend * seasonal), linetype = 2)

1.6 Summary of modern R commands used in examples