gdp_pc <- global_economy |>
mutate(GDP_per_capita = GDP / Population)
gdp_pc |>
autoplot(GDP_per_capita, show.legend = FALSE) +
labs(title = "GDP per capita, all countries", y = "US$")
# Highest single-year GDP per capita on record
gdp_pc |>
filter(!is.na(GDP_per_capita)) |>
filter(GDP_per_capita == max(GDP_per_capita)) |>
select(Country, Year, GDP_per_capita)
## # A tsibble: 1 x 3 [1Y]
## # Key: Country [1]
## Country Year GDP_per_capita
## <fct> <dbl> <dbl>
## 1 Monaco 2014 185153.
# Highest in the most recent year (2017)
gdp_pc |>
filter(Year == 2017, !is.na(GDP_per_capita)) |>
slice_max(GDP_per_capita, n = 5) |>
select(Country, GDP_per_capita)
## # A tsibble: 5 x 3 [1Y]
## # Key: Country [5]
## Country GDP_per_capita Year
## <fct> <dbl> <dbl>
## 1 Luxembourg 104103. 2017
## 2 Macao SAR, China 80893. 2017
## 3 Switzerland 80190. 2017
## 4 Norway 75505. 2017
## 5 Iceland 70057. 2017
Monaco has the highest GDP per capita for most of the record, with Liechtenstein also among the top. GDP per capita has risen over time for nearly all countries, and the spread between rich and poor countries has widened.
global_economy |>
filter(Country == "United States") |>
autoplot(GDP) +
labs(title = "US GDP", y = "US$")
US GDP grows smoothly and roughly exponentially with stable relative variation, so a variance-stabilising transformation is not needed (dividing by population to get per-capita GDP would be a reasonable adjustment).
aus_livestock |>
filter(Animal == "Bulls, bullocks and steers", State == "Victoria") |>
autoplot(Count) +
labs(title = "Victorian bulls, bullocks and steers slaughtered")
The variability is roughly constant over the series, so no transformation is needed.
vic_elec |>
autoplot(Demand) +
labs(title = "Victorian electricity demand", y = "MWh")
The magnitude of fluctuations is fairly stable, so a transformation is not helpful here — the modelling challenge is the multiple seasonal periods, not variance.
aus_production |>
autoplot(Gas) +
labs(title = "Australian gas production (untransformed)")
lambda_gas <- aus_production |>
features(Gas, features = guerrero) |>
pull(lambda_guerrero)
lambda_gas
## [1] 0.1095171
aus_production |>
autoplot(box_cox(Gas, lambda_gas)) +
labs(title = paste0("Box-Cox transformed gas (lambda = ",
round(lambda_gas, 2), ")"))
Gas production shows seasonal swings that grow with the level of the series, so a Box-Cox transformation is appropriate. The Guerrero method chooses a small lambda (close to a log transform), which evens out the seasonal variation.
canadian_gas |>
autoplot(Volume) +
labs(title = "Canadian gas production", y = "billion cubic metres")
canadian_gas |> gg_season(Volume)
A Box-Cox transformation is unhelpful because the seasonal variation is not monotonically related to the level of the series. The seasonal fluctuations are smallest when volume is low (early years) and when it is high (later years), and largest in the middle of the series. A Box-Cox transformation can only stabilise variance that increases or decreases steadily with the level, so no single value of lambda fixes this non-monotonic pattern.
set.seed(12345678)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
lambda_retail <- myseries |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero)
lambda_retail
## [1] 0.08303631
myseries |>
autoplot(box_cox(Turnover, lambda_retail)) +
labs(title = paste0("Box-Cox transformed retail turnover (lambda = ",
round(lambda_retail, 2), ")"))
The Guerrero method selects the lambda shown above. It is well below 1, meaning a strong transformation is needed because the seasonal variation in this retail series grows with the level; the transformed series has much more even variance.
# Tobacco (drop the missing tail values)
lambda_tobacco <- aus_production |>
filter(!is.na(Tobacco)) |>
features(Tobacco, features = guerrero) |>
pull(lambda_guerrero)
# Economy class passengers, Melbourne - Sydney
lambda_ansett <- ansett |>
filter(Class == "Economy", Airports == "MEL-SYD") |>
features(Passengers, features = guerrero) |>
pull(lambda_guerrero)
# Pedestrian counts at Southern Cross Station.
# The hourly series contains zeros (overnight), for which a Box-Cox transform
# is undefined, so we aggregate to daily totals first.
ped_daily <- pedestrian |>
filter(Sensor == "Southern Cross Station") |>
index_by(Date) |>
summarise(Count = sum(Count))
lambda_ped <- ped_daily |>
features(Count, features = guerrero) |>
pull(lambda_guerrero)
tibble(
series = c("Tobacco", "ansett Economy MEL-SYD", "Southern Cross pedestrians (daily)"),
lambda = c(lambda_tobacco, lambda_ansett, lambda_ped)
)
## # A tibble: 3 × 2
## series lambda
## <chr> <dbl>
## 1 Tobacco 0.926
## 2 ansett Economy MEL-SYD 2.00
## 3 Southern Cross pedestrians (daily) 0.273
aus_production |> filter(!is.na(Tobacco)) |>
autoplot(box_cox(Tobacco, lambda_tobacco)) +
labs(title = paste0("Tobacco, Box-Cox (lambda = ", round(lambda_tobacco, 2), ")"))
ansett |> filter(Class == "Economy", Airports == "MEL-SYD") |>
autoplot(box_cox(Passengers, lambda_ansett)) +
labs(title = paste0("Economy MEL-SYD, Box-Cox (lambda = ",
round(lambda_ansett, 2), ")"))
ped_daily |>
autoplot(box_cox(Count, lambda_ped)) +
labs(title = paste0("Southern Cross pedestrians (daily), Box-Cox (lambda = ",
round(lambda_ped, 2), ")"))
The Guerrero-selected lambdas differ a lot by series. Tobacco returns a lambda near 1, meaning almost no transformation is needed — its variance is already fairly stable. Economy MEL-SYD hits the upper boundary (lambda near 2) because that series contains a structural break (the 1989 pilots’ strike drives passengers to near zero), so Box-Cox is not really the right tool there. For the daily pedestrian counts, a moderate lambda evens out the variance well once the series is aggregated to daily totals.
gas <- tail(aus_production, 5 * 4) |> select(Gas)
autoplot(gas, Gas) +
labs(title = "Australian gas production (last 5 years)")
(a) There is a clear upward trend and strong seasonality — production peaks in Q3 (winter) and troughs in Q1 each year.
(b) Classical multiplicative decomposition:
dcmp <- gas |>
model(classical_decomposition(Gas, type = "multiplicative")) |>
components()
autoplot(dcmp) +
labs(title = "Classical multiplicative decomposition of gas")
(c) Yes — the decomposition confirms part (a): the trend-cycle rises steadily and the seasonal indices are above 1 in Q3 and below 1 in Q1.
(d) Seasonally adjusted series:
plot_seasadj <- function(data, title) {
data |>
model(classical_decomposition(Gas, type = "multiplicative")) |>
components() |>
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Gas, colour = "Data")) +
geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
scale_colour_manual(values = c("Data" = "grey60",
"Seasonally adjusted" = "#0072B2"),
name = NULL) +
labs(title = title, y = "Gas")
}
plot_seasadj(gas, "Seasonally adjusted gas")
# Outlier in the middle
gas_mid <- gas |>
mutate(Gas = if_else(row_number() == 10, Gas + 300, Gas))
plot_seasadj(gas_mid, "Outlier in the middle (+300 at obs 10)")
(e) The outlier passes almost entirely into the remainder, so the seasonally adjusted series shows a large spike at that point. It also distorts the trend-cycle and seasonal indices near the outlier.
# Outlier near the end
gas_end <- gas |>
mutate(Gas = if_else(row_number() == 19, Gas + 300, Gas))
plot_seasadj(gas_end, "Outlier near the end (+300 at obs 19)")
(f) Yes. Classical decomposition cannot estimate the trend-cycle for the first and last two quarters, so an outlier near the end falls where the trend is missing and contaminates fewer seasonal estimates, whereas a middle outlier distorts both the trend-cycle and the seasonal component around it. Either way a spike appears in the seasonally adjusted series.
library(seasonal)
x11_dcmp <- myseries |>
model(x11 = X_13ARIMA_SEATS(Turnover ~ x11())) |>
components()
autoplot(x11_dcmp) +
labs(title = "X-11 decomposition of retail turnover")
x11_dcmp |>
ggplot(aes(x = Month, y = irregular)) +
geom_line() +
labs(title = "X-11 irregular (remainder) component")
The X-11 irregular component makes outliers easy to spot as spikes that stand out from the otherwise flat remainder. It also shows that the seasonal component evolves over time (its shape changes across the series), which the single fixed seasonal shape of a classical decomposition would hide.
(a) The decomposition is dominated by the trend component, which shows the Australian civilian labour force growing steadily from 1978 to 1995. The seasonal component is small — the grey scale bar on its panel is much larger than on the trend panel, indicating the seasonal swings are tiny relative to the overall level. The remainder is also small except around 1991-1992. The seasonal pattern is regular and changes only slowly over the period.
(b) Yes. The 1991/1992 recession is visible mainly in the remainder component, which shows a pronounced dip at that time, and as a flattening of the otherwise steadily rising trend. It is not visible in the seasonal component.
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin25.4.0
## Running under: macOS Tahoe 26.6.2
##
## Matrix products: default
## BLAS: /opt/homebrew/Cellar/openblas/0.3.34/lib/libopenblasp-r0.3.34.dylib
## LAPACK: /opt/homebrew/Cellar/r/4.6.1/lib/R/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] seasonal_1.11.0 fable_0.5.0 feasts_0.5.0 fabletools_0.8.0
## [5] ggtime_1.0.0 tsibbledata_0.4.1 tsibble_1.2.0 ggplot2_4.0.3
## [9] lubridate_1.9.5 tidyr_1.3.2 dplyr_1.2.1 tibble_3.3.1
## [13] fpp3_1.0.3
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 rappdirs_0.3.4 utf8_1.2.6
## [4] generics_0.1.4 anytime_0.3.13 digest_0.6.39
## [7] magrittr_2.0.5 evaluate_1.0.5 grid_4.6.1
## [10] timechange_0.4.0 RColorBrewer_1.1-3 mixtime_0.3.0
## [13] fastmap_1.2.0 jsonlite_2.0.0 tinytex_0.60
## [16] purrr_1.2.2 scales_1.4.0 jquerylib_0.1.4
## [19] cli_3.6.6 rlang_1.3.0 x13binary_1.1.61.2
## [22] crayon_1.5.3 vecvec_1.3.0 cachem_1.1.0
## [25] withr_3.0.3 yaml_2.3.12 tools_4.6.1
## [28] tzdb_0.5.0 vctrs_0.7.3 R6_2.6.1
## [31] lifecycle_1.0.5 pkgconfig_2.0.3 bslib_0.12.0
## [34] progressr_1.0.0 pillar_1.11.1 gtable_0.3.6
## [37] glue_1.8.1 Rcpp_1.1.2 xfun_0.60
## [40] tidyselect_1.2.1 rstudioapi_0.19.0 knitr_1.52
## [43] farver_2.1.2 htmltools_0.5.9 rmarkdown_2.32
## [46] labeling_0.4.3 compiler_4.6.1 S7_0.2.2
## [49] distributional_0.9.0