knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(cli.unicode = FALSE)
library(fpp3)
library(seasonal)

Exercise 3.1

Consider the GDP information in global_economy. Plot the GDP per capita for each country over time. Which country has the highest GDP per capita? How has this changed over time?

gdp_pc <- global_economy |>
  mutate(GDP_per_capita = GDP / Population)

gdp_pc |>
  autoplot(GDP_per_capita, show.legend = FALSE) +
  labs(title = "GDP per capita by country", y = "US$")

Country with the highest GDP per capita:

gdp_pc |>
  as_tibble() |>
  filter(!is.na(GDP_per_capita)) |>
  slice_max(GDP_per_capita, n = 5) |>
  select(Country, Year, GDP_per_capita)
## # A tibble: 5 x 3
##   Country        Year GDP_per_capita
##   <fct>         <dbl>          <dbl>
## 1 Monaco         2014        185153.
## 2 Monaco         2008        180640.
## 3 Liechtenstein  2014        179308.
## 4 Liechtenstein  2013        173528.
## 5 Monaco         2013        172589.

How the leader has changed over time:

leaders <- gdp_pc |>
  as_tibble() |>
  filter(!is.na(GDP_per_capita)) |>
  group_by(Year) |>
  slice_max(GDP_per_capita, n = 1) |>
  ungroup()

leaders |> count(Country, sort = TRUE)
## # A tibble: 6 x 2
##   Country                  n
##   <fct>                <int>
## 1 Monaco                  43
## 2 United States            8
## 3 Kuwait                   2
## 4 Liechtenstein            2
## 5 United Arab Emirates     2
## 6 Luxembourg               1
gdp_pc |>
  filter(Country %in% c("Monaco", "Liechtenstein", "Luxembourg",
                        "United States", "Kuwait", "United Arab Emirates")) |>
  autoplot(GDP_per_capita) +
  labs(title = "GDP per capita of the countries that have ranked first", y = "US$")

Monaco has the highest GDP per capita, peaking at about $185,000 in 2014. It ranked first in 43 of the 58 years.

Exercise 3.2

For each of the following series, make a graph of the data. If transforming seems appropriate, do so and describe the effect.

United States GDP from global_economy

us_gdp <- global_economy |> filter(Country == "United States")
us_gdp |> autoplot(GDP) + labs(title = "US GDP")

lambda <- us_gdp |> features(GDP, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] 0.2819443
us_gdp |>
  autoplot(box_cox(GDP, lambda)) +
  labs(title = paste("US GDP, Box-Cox transformed, lambda =", round(lambda, 2)))

US GDP grows exponentially, so the curve bends upward. A Box-Cox transformation with lambda of about 0.28 straightens it into a nearly linear trend, which will be easier to model and forecast.

Victorian “Bulls, bullocks and steers” from aus_livestock

bulls <- aus_livestock |>
  filter(State == "Victoria", Animal == "Bulls, bullocks and steers")
bulls |> autoplot(Count) + labs(title = "Victorian bulls, bullocks and steers slaughtered")

lambda <- bulls |> features(Count, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] -0.04461887
bulls |>
  autoplot(log(Count)) +
  labs(title = "Log of Victorian bulls, bullocks and steers slaughtered")

The variation is larger in the late 1970s, when the counts were highest, and smaller when counts are lower. The Guerrero lambda is close to 0, so a log transformation is appropriate. It makes the size of the fluctuations more even across the series.

Victorian electricity demand from vic_elec

vic_elec |> autoplot(Demand) + labs(title = "Victorian electricity demand")

lambda <- vic_elec |> features(Demand, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] 0.09993089
vic_elec |>
  autoplot(box_cox(Demand, lambda)) +
  labs(title = paste("Victorian electricity demand, Box-Cox, lambda =", round(lambda, 2)))

The demand series has no trend. Its variation comes from daily, weekly and yearly seasonality and from temperature spikes (heat waves in summer), not from the level of the series. The transformed plot looks almost the same as the original, so a transformation is not really needed here.

Gas production from aus_production

aus_production |> autoplot(Gas) + labs(title = "Australian gas production")

lambda <- aus_production |> features(Gas, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] 0.1095171
aus_production |>
  autoplot(box_cox(Gas, lambda)) +
  labs(title = paste("Australian gas production, Box-Cox, lambda =", round(lambda, 2)))

The seasonal swings in gas production grow as the level of the series grows. A Box-Cox transformation with lambda of about 0.11 (close to a log) makes the seasonal variation roughly the same size throughout the series.

Exercise 3.3

Why is a Box-Cox transformation unhelpful for the canadian_gas data?

canadian_gas |> autoplot(Volume) + labs(title = "Canadian gas production")

canadian_gas |> gg_subseries(Volume)

A Box-Cox transformation works when the variation increases (or decreases) steadily with the level of the series. In canadian_gas that is not the case. The seasonal variation is small in the 1960s, becomes much larger from about 1975 to 1990, and then becomes smaller again in the 1990s and 2000s, even though production keeps rising. Because the size of the seasonal swings is not tied to the level, no single power transformation can make the variation constant.

Exercise 3.4

What Box-Cox transformation would you select for your retail data (from Exercise 7 in Section 2.10)?

set.seed(624)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
distinct(as_tibble(myseries), State, Industry)
## # A tibble: 1 x 2
##   State           Industry              
##   <chr>           <chr>                 
## 1 New South Wales Takeaway food services
myseries |> autoplot(Turnover) + labs(title = "Retail turnover")

lambda <- myseries |> features(Turnover, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] 0.002144737
myseries |> autoplot(log(Turnover)) + labs(title = "Log of retail turnover")

My series is takeaway food services turnover in New South Wales. The seasonal swings get bigger as turnover grows. The Guerrero method gives a lambda of about 0.002, which is essentially 0, so I would use a log transformation. After taking logs, the seasonal variation is about the same size across the whole series.

Exercise 3.5

For the following series, find an appropriate Box-Cox transformation in order to stabilise the variance.

Tobacco from aus_production

aus_production |> autoplot(Tobacco) + labs(title = "Tobacco production")

aus_production |> features(Tobacco, features = guerrero)
## # A tibble: 1 x 1
##   lambda_guerrero
##             <dbl>
## 1           0.926

The Guerrero lambda is about 0.93, which is very close to 1. A lambda of 1 means no transformation, so the tobacco series does not need one. Its variation is already fairly constant.

Economy class passengers between Melbourne and Sydney from ansett

mel_syd <- ansett |>
  filter(Airports == "MEL-SYD", Class == "Economy")
mel_syd |> autoplot(Passengers) + labs(title = "Economy passengers, Melbourne-Sydney")

lambda <- mel_syd |> features(Passengers, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] 1.999927
mel_syd |>
  autoplot(box_cox(Passengers, lambda)) +
  labs(title = paste("Economy passengers, Box-Cox, lambda =", round(lambda, 2)))

The Guerrero method suggests a lambda of about 2. This value is driven by unusual periods rather than by variance that changes with the level. Passengers drop to zero during the 1989 pilots’ strike, and there is another dip in 1992. Apart from those periods, the variation is fairly steady, so the transformation does not change much. A lambda of 1 (no transformation) would also be reasonable.

Pedestrian counts at Southern Cross Station from pedestrian

southern_cross <- pedestrian |>
  filter(Sensor == "Southern Cross Station")
southern_cross |> autoplot(Count) + labs(title = "Hourly pedestrian counts, Southern Cross Station")

lambda <- southern_cross |> features(Count + 1, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] -0.2566589
southern_cross |>
  autoplot(box_cox(Count + 1, lambda)) +
  labs(title = paste("Pedestrian counts, Box-Cox, lambda =", round(lambda, 2)))

The hourly counts include some zeros (late at night), so I added 1 before transforming. The Guerrero lambda is about -0.25, a strong transformation similar to a log. It compresses the very high peak-hour counts and spreads out the low overnight counts, which makes the variation more even.

Exercise 3.7

Consider the last five years of the Gas data from aus_production.

gas <- tail(aus_production, 5 * 4) |> select(Gas)

a. Plot the time series. Can you identify seasonal fluctuations and/or a trend-cycle?

gas |> autoplot(Gas) + labs(title = "Gas production, last five years")

There is strong seasonality. Production is highest in Q3 (winter) and lowest in Q1 every year. There is also a mild upward trend.

b. Use classical_decomposition with type = multiplicative to calculate the trend-cycle and seasonal indices.

gas_dcmp <- gas |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components()

gas_dcmp |> autoplot()

gas_dcmp |>
  as_tibble() |>
  distinct(quarter = quarter(Quarter), seasonal) |>
  arrange(quarter)
## # A tibble: 4 x 2
##   quarter seasonal
##     <int>    <dbl>
## 1       1    0.875
## 2       2    1.07 
## 3       3    1.13 
## 4       4    0.925

c. Do the results support the graphical interpretation from part a?

Yes. The seasonal indices show Q3 about 13% above average and Q2 about 7% above, while Q4 is about 7% below and Q1 about 12% below. This matches the pattern in the plot. The trend-cycle rises from about 200 to 226, which confirms the upward trend. The trend levels off during 2008.

d. Compute and plot the seasonally adjusted data.

gas_dcmp |>
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "Data")) +
  geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
  labs(title = "Gas production and seasonally adjusted gas production",
       y = "Petajoules", colour = NULL)

e. Change one observation to be an outlier and recompute the seasonally adjusted data. What is the effect of the outlier?

I added 300 to the 10th observation (2007 Q4), which is in the middle of the series.

gas_mid <- gas
gas_mid$Gas[10] <- gas_mid$Gas[10] + 300

gas_mid_dcmp <- gas_mid |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components()

gas_mid_dcmp |>
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "Data")) +
  geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
  labs(title = "Outlier in the middle (2007 Q4)", y = "Petajoules", colour = NULL)

The outlier has a big effect:

  • The seasonally adjusted series has a huge spike at 2007 Q4, because the outlier is not seasonal and passes straight through to the adjusted data.
  • The trend-cycle is pulled up in the quarters around the outlier.
  • The seasonal indices change as well. The Q4 index jumps from about 0.93 to 1.12, so Q4 now looks like the highest quarter, which is wrong. Because of this, the seasonally adjusted values for every Q4 in the series are pushed down, and the other quarters are pushed up. The whole seasonally adjusted series becomes more jagged, not just the point with the outlier.

f. Does it make any difference if the outlier is near the end rather than in the middle of the time series?

gas_end <- gas
gas_end$Gas[20] <- gas_end$Gas[20] + 300

gas_end_dcmp <- gas_end |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components()

gas_end_dcmp |>
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "Data")) +
  geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
  labs(title = "Outlier at the end (2010 Q2)", y = "Petajoules", colour = NULL)

Yes. When the outlier is at the end, the seasonally adjusted series has a spike at that last point, but the rest of the series is almost unchanged. Classical decomposition cannot estimate the trend for the last two quarters (the moving average needs values on both sides), so the outlier barely enters the seasonal index calculations. The seasonal indices change only slightly. An outlier in the middle distorts the trend and the seasonal indices, and so it affects the whole seasonally adjusted series.

Exercise 3.8

Recall your retail time series data (from Exercise 7 in Section 2.10). Decompose the series using X-11. Does it reveal any outliers, or unusual features that you had not noticed previously?

x11_dcmp <- myseries |>
  model(x11 = X_13ARIMA_SEATS(Turnover ~ x11())) |>
  components()

autoplot(x11_dcmp) +
  labs(title = "X-11 decomposition of NSW takeaway food turnover")

x11_dcmp |> gg_subseries(seasonal)

The X-11 decomposition shows several things that were not obvious from the time plot:

Exercise 3.9

Figures 3.19 and 3.20 show the result of decomposing the number of persons in the civilian labour force in Australia each month from February 1978 to August 1995.

a. Write about 3–5 sentences describing the results of the decomposition. Pay particular attention to the scales of the graphs in making your interpretation.

The civilian labour force grew steadily from about 6,400 thousand people in 1978 to about 9,000 thousand in 1995, and the trend accounts for almost all of the movement in the data. The seasonal component is small by comparison: it only ranges from about -100 to +100 thousand, while the trend changes by about 2,600 thousand over the period. The grey bars on the left show this: the bar on the seasonal panel is much longer than the one on the trend panel, which means the seasonal panel is drawn on a much smaller scale. The seasonal pattern also changes over time. For example, March and April rose until the mid-1980s and then fell, August and November became more negative, and September and December became more positive. Most remainder values are small (within about ±50 thousand), except for a few large negative values in 1991 and 1992.

b. Is the recession of 1991/1992 visible in the estimated components?

Yes. The trend flattens out from about 1991 to 1993, after growing steadily before that. The recession is clearest in the remainder, which has very large negative values in 1991 (down to about -400 thousand) and more negative values in 1992 and 1993. These are far bigger than the remainder values anywhere else in the series. The decomposition could not absorb the sudden drop in the labour force into the smooth trend, so it shows up in the remainder.