DATA 624 Homework: Chapter 3 Exercises
Exercises 3.1, 3.2, 3.3, 3.4, 3.5, 3.7, 3.8 and 3.9
Jonathan Cruz
September 20, 2026
library(fpp3) # loads tsibble, tsibbledata, feasts, fable, dplyr, ggplot2, ...
#below is a small helper that will be used multiple times, it takes a tsbibble which takes
# a tstible, a column and a Box-cox lambda then it draws the original vs transformed series
#stacked so both can be compared in the same axis
plot_boxcox <- function(data, y, lambda, title = NULL) {
data |>
# box_cox() is a fabletools function: (y^lambda - 1)/lambda, or log(y) if lambda = 0
mutate(Transformed = box_cox({{ y }}, lambda),
Original = {{ y }}) |>
# pivot_longer stacks both columns so ggplot can facet them
pivot_longer(c(Original, Transformed), names_to = "Series", values_to = "Value") |>
# keep the panel order original then it transform
mutate(Series = factor(Series, levels = c("Original", "Transformed"))) |>
ggplot(aes(x = !!index(data), y = Value)) +
geom_line() +
# scales = "free_y" because the two panels are on completely different scales
facet_wrap(~ Series, ncol = 1, scales = "free_y") +
labs(title = title,
subtitle = paste0("Box cox transformation with lambda = ", round(lambda, 3)),
x = NULL)
}Exercise 3.1 — GDP per capita (global_economy)
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?
# mutate() adds a column to the dataset because the object is a tsibble
# it has the index Year and the key Counter that is preserves
# so autoplot() knows to draw one line per country.
gdp_pc <- global_economy |>
mutate(GDP_per_capita = GDP / Population)
gdp_pc |>
autoplot(GDP_per_capita) + # one line per key
theme(legend.position = "none") + # 263 legend entries would make it unreadable
labs(title = "GDP per capita, all the countries in global_economy",
y = "US dollars (prices current)", x = "Year")
top_by_year <- gdp_pc |>
as_tibble() |> # just tibble we no longer need time-series(tstibble) rules
filter(!is.na(GDP_per_capita)) |> # some countries have missing GDP or population
group_by(Year) |>
slice_max(GDP_per_capita, n = 1) |> # the single highest country in each year
ungroup() |>
select(Year, Country, GDP_per_capita)
# How many years did each country hold the number 1 spot for ?
top_by_year |> count(Country, sort = TRUE)## # A tibble: 6 × 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# Show the most recent 15 years
top_by_year |> filter(Year >= 2003)## # A tibble: 15 × 3
## Year Country GDP_per_capita
## <dbl> <fct> <dbl>
## 1 2003 Monaco 108978.
## 2 2004 Monaco 123382.
## 3 2005 Monaco 124374.
## 4 2006 Monaco 133195.
## 5 2007 Monaco 167125.
## 6 2008 Monaco 180640.
## 7 2009 Monaco 149221.
## 8 2010 Monaco 144569.
## 9 2011 Monaco 162155.
## 10 2012 Monaco 152000.
## 11 2013 Liechtenstein 173528.
## 12 2014 Monaco 185153.
## 13 2015 Liechtenstein 167591.
## 14 2016 Monaco 168011.
## 15 2017 Luxembourg 104103.# The highest country year in the whole dataset
gdp_pc |>
as_tibble() |>
slice_max(GDP_per_capita, n = 3) |>
select(Year, Country, GDP_per_capita)## # A tibble: 3 × 3
## Year Country GDP_per_capita
## <dbl> <fct> <dbl>
## 1 2014 Monaco 185153.
## 2 2008 Monaco 180640.
## 3 2014 Liechtenstein 179308.leaders <- c("Monaco", "Liechtenstein", "Luxembourg", "Switzerland", "Norway",
"United States", "Kuwait", "United Arab Emirates")
gdp_pc |>
filter(Country %in% leaders) |> # holding on to only countries of interest
autoplot(GDP_per_capita) +
labs(title = "GDP per capita: countries that have led at some point",
y = "US dollars (current prices)", x = "Year", colour = "Country")
Answer. Monaco has had the highest GDP per capita for the most part, it’s been number 1 for 43 years out of all 58 and it holds the all time maximum. Overall GDP per capita has risen by large margins since 1960. United States led for the early part of 1960s and Kuwait actually lead during 1965 and 1966 and basically from there Monaco has dominated.
Exercise 3.2 — Graph each series; transform if appropriate
For each of the following series, make a graph of the data. If transforming seems appropriate, do so and describe the effect.
(a) United States GDP
us_gdp <- global_economy |>
filter(Code == "USA") # Code is the 3 letter code
us_gdp |>
autoplot(GDP / 1e12) + # divide by 1e12 so the axis reads in trillions
labs(title = "United States GDP", y = "Trillions of US dollars", x = "Year")
# Guerrero's method picks the lambda that stabilises the variance
(lambda_us <- us_gdp |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero))## [1] 0.2819443plot_boxcox(us_gdp, GDP, lambda_us, "United States GDP")
Effect. US GDP is growth is sort of exponential, the curve bends upward more and more as the year to year goes. After the transformation we see more of a straight line which suggest the growth being more additive rather than multiplicative and since the data is anual there’s no seasonality so were looking more to linearise the trend rather then stabilizing it for seasonal swings.
(b) Slaughter of Victorian “Bulls, bullocks and steers”
vic_bulls <- aus_livestock |>
filter(Animal == "Bulls, bullocks and steers", # one animal category
State == "Victoria") # one state
vic_bulls |>
autoplot(Count / 1e3) +
labs(title = "Victoria -slaughter of bulls, bullocks and steers",
y = "Thousands of head", x = "Month")
(lambda_bulls <- vic_bulls |>
features(Count, features = guerrero) |>
pull(lambda_guerrero))## [1] -0.04461887plot_boxcox(vic_bulls, Count, lambda_bulls, "Victoria - bulls, bullocks and steers")
Effect. This series takes a dive in the late 1970s and then it steps down around the 2000s again and then it has stabalized ear the 50,000s ever since. Swings are larger when level are highers and more stable on the lower levels. Since Guerrero is lambda was ~ 0 this would suggest a log transformation making the variation look constant as we see in our transformed chart That means that the effect is real but minimal so I would just use the original.
(c) Victorian electricity demand (vic_elec)
vic_elec |>
autoplot(Demand) +
labs(title = "Victorian half-hourly electricity demand", y = "MWh", x = "Time")
# index_by() is the tsibble version of group_by() in tibble but for the time index
# here we collapse the half-hourly index down to a daily one.
vic_daily <- vic_elec |>
index_by(Date = as_date(Time)) |> # new index = calendar date
summarise(Demand = sum(Demand)) # total daily demand
vic_daily |>
autoplot(Demand / 1e3) +
labs(title = "Victorian daily electricity demand", y = "GWh per day", x = "Date")
(lambda_elec <- vic_elec |>
features(Demand, features = guerrero) |>
pull(lambda_guerrero))## [1] 0.09993089plot_boxcox(vic_daily, Demand, lambda_elec, "Victorian daily electricity demand")
Effect. The demand for electricity has many seasonal patters, which are daily, weekly and year which we can assume summer spikes are due to air conditioning and the winder might be due to heat pumps. We see no trend here. After transformation, the plot looks virtually the same so I would just leave this series untransformed.
(d) Gas production (aus_production)
aus_production |>
autoplot(Gas) +
labs(title = "Australian quarterly gas production", y = "Petajoules", x = "Quarter")
(lambda_gas <- aus_production |>
features(Gas, features = guerrero) |>
pull(lambda_guerrero))## [1] 0.1095171plot_boxcox(aus_production, Gas, lambda_gas, "Australian gas production")
Effect.Now this is textbook transformation, look at that 20 fold from 1956 to 2010 and seasonal swings are almost completely proportional. We see that Guerrero is lambda yeilds a near 0 so were going with log for this one. The Box Cox transformation shows uo a constant seasonal swing during the whole series and the trend is closer to linear. The variance has been stabalised which is what we are looking for when trying to fit an additive model
Exercise 3.3 — Why is Box-Cox unhelpful for
canadian_gas?
canadian_gas |>
autoplot(Volume) +
labs(title = "Monthly Canadian gas production", y = "Billions of cubic metres", x = "Month")
(lambda_cg <- canadian_gas |>
features(Volume, features = guerrero) |>
pull(lambda_guerrero))## [1] 0.5767648plot_boxcox(canadian_gas, Volume, lambda_cg, "Canadian gas production")
Answer. A Box-Cox transformation helps with the sizing of variations however they are fixed. Smaller levels yeild smaller swings and higher levels yeild larger swings so the issue becomes that the middle it’s in between large swings in the beggining and smaller swings in the end sort of fighting between both to find it’s sizing between two swing sizes for levels.
Exercise 3.4 — Box-Cox for my retail series
What Box-Cox transformation would you select for your retail data (from Exercise 7 in Section 2.10)?
set.seed(12345678) #iill use same see as example from HW #1
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`, 1)) # one random Series ID
# Which series did we get?
myseries |> as_tibble() |> distinct(State, Industry, `Series ID`)## # A tibble: 1 × 3
## State Industry `Series ID`
## <chr> <chr> <chr>
## 1 Northern Territory Clothing, footwear and personal accessory reta… A3349767Wmyseries |>
autoplot(Turnover) +
labs(title = paste(unique(myseries$State), "-", unique(myseries$Industry)),
y = "Turnover ($ million AUD)", x = "Month")
(lambda_retail <- myseries |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero))## [1] 0.08303631plot_boxcox(myseries, Turnover, lambda_retail,
paste(unique(myseries$State), "-", unique(myseries$Industry)))
Answer. For the Northern Terrirorry series, the Box-Cox transformation I choose is log and here’s why. Guerrero recommends lambda as 0.08 which is very close to zero. The seasonal swings are the same across all 30 year and plus it’s easier to explain as logs are just changes in percentages.
Exercise 3.5 — Find a Box-Cox transformation to stabilise the variance
Tobacco from
aus_production, Economy class passengers between Melbourne and Sydney fromansett, and Pedestrian counts at Southern Cross Station frompedestrian.
(a) Tobacco (aus_production)
tobacco <- aus_production |>
filter(!is.na(Tobacco)) |> # drop the trailing NAs because Tobacco has missing values after 2004 Q2
select(Quarter, Tobacco)
tobacco |>
autoplot(Tobacco) +
labs(title = "Australian quarterly tobacco and cigarette production",
y = "Tonnes", x = "Quarter")
(lambda_tob <- tobacco |>
features(Tobacco, features = guerrero) |>
pull(lambda_guerrero))## [1] 0.9264636plot_boxcox(tobacco, Tobacco, lambda_tob, "Tobacco production")
Answer. Guerrero’s lambda practically gives us 1 and it makes complete sense. The seasonal swings are essentially the same wether series is rising or falling so the variance is stable already. The transformed panel is very identical to the original so I would leave it unstransformed.
(b) Economy class passengers, Melbourne–Sydney
(ansett)
mel_syd <- ansett |>
filter(Class == "Economy", Airports == "MEL-SYD")
mel_syd |>
autoplot(Passengers / 1e3) +
labs(title = "Ansett economy-class passengers, Melbourne–Sydney (weekly)",
y = "Thousands of passengers", x = "Week")
# How many weeks have zero passengers?
mel_syd |> as_tibble() |> summarise(zero_weeks = sum(Passengers == 0))## # A tibble: 1 × 1
## zero_weeks
## <int>
## 1 7(lambda_ansett <- mel_syd |>
features(Passengers, features = guerrero) |>
pull(lambda_guerrero))## [1] 1.999927plot_boxcox(mel_syd, Passengers, lambda_ansett, "Ansett MEL-SYD economy passengers")
Answer. Guerrero’s lambda is practically a 2 which is in the upper bounds so it’s telling us that the variance is not due to the levels. The problem is more structural than anything, you have seven weeks of zero passengers, nose-first dive on holidays . Since lambda is 2, the transformation occurs by squaring the data which exagerates peaks and does nothing to the zeros. I would stay away from transforming in this case as it brings more confusion than anything.
(c) Pedestrian counts at Southern Cross Station
(pedestrian)
The pedestrian data are hourly counts,
and a train station is empty in the middle of the night, so the series
contains genuine zeros. That matters for Box-Cox: with \(\lambda < 0\),
box_cox(0, lambda) is undefined (it produces
NaN), and \(\lambda = 0\)
(log) gives \(-\infty\).
# Pedestrian data are hourly_counts while the train is empty in the middle of the night so this series contains real zeros.
# this matters for Box-Cox because ifi lambda < 0, box_cox(0, lambda) will be undefined producing NAN which then goes to negative infinity
sth_cross <- pedestrian |>
filter(Sensor == "Southern Cross Station")
sth_cross |>
autoplot(Count) +
labs(title = "Hourly pedestrian counts, Southern Cross Station (2015–2016)",
y = "Pedestrians per hour", x = "Time")
# how many hours of zero?
sth_cross |> as_tibble() |> summarise(hours = n(), zero_hours = sum(Count == 0))## # A tibble: 1 × 2
## hours zero_hours
## <int> <int>
## 1 17539 159(lambda_ped <- sth_cross |>
features(Count, features = guerrero) |>
pull(lambda_guerrero))## [1] -0.2501616# with a negative lambda the zero-count hours become NaN and are simply
# dropped from the plot and this is a warning ggplot normally shows.
plot_boxcox(sth_cross, Count, lambda_ped, "Southern Cross Station: hourly counts")
Because the hourly plot is so compact, we will also look at daily totals, which removes the zeros and shows the weekly pattern and any variance change more clearly.
sth_daily <- sth_cross |>
index_by(Date) |> # Date already exists as a column in pedestrian
summarise(Count = sum(Count)) # total pedestrians per day
sth_daily |>
autoplot(Count / 1e3) +
labs(title = "Daily pedestrian counts, Southern Cross Station",
y = "Thousands per day", x = "Date")
(lambda_ped_daily <- sth_daily |>
features(Count, features = guerrero) |>
pull(lambda_guerrero))## [1] 0.2726316plot_boxcox(sth_daily, Count, lambda_ped_daily, "Southern Cross Station: daily counts")
Answer. Guerrero picks a negative lambda for the hourly data which shows the strong right skew, most hours are quiet and then you have rush hour and the values are enormous during those times.The transformation makes daily make daily cycles look more even by compressing the rush hour peaks but this negative lambda made it arwarkd because NAN can’t be transformed so the better choice is the lambda Guerrero chooses for the daily totals which stabalized variance with no problem, anyhow the goal was reached by transformation. Weekdays and weekends swings are siilar sizes and not dominated by weekday peaks
Exercise 3.7 — Classical decomposition of the last five years of Gas
# tail(x, 20) keeps the last 20 quarters which are 4 quarters in 5 years so 5 * 4,
gas <- tail(aus_production, 5 * 4) |> select(Gas)
gas## # A tsibble: 20 x 2 [1Q]
## Gas Quarter
## <dbl> <qtr>
## 1 221 2005 Q3
## 2 180 2005 Q4
## 3 171 2006 Q1
## 4 224 2006 Q2
## 5 233 2006 Q3
## 6 192 2006 Q4
## 7 187 2007 Q1
## 8 234 2007 Q2
## 9 245 2007 Q3
## 10 205 2007 Q4
## 11 194 2008 Q1
## 12 229 2008 Q2
## 13 249 2008 Q3
## 14 203 2008 Q4
## 15 196 2009 Q1
## 16 238 2009 Q2
## 17 252 2009 Q3
## 18 210 2009 Q4
## 19 205 2010 Q1
## 20 236 2010 Q2(a) Plot the series
gas |>
autoplot(Gas) +
labs(title = "Australian gas production, last five years (2005 Q3 – 2010 Q2)",
y = "Petajoules", x = "Quarter")
Answer. Yes to both. Clear seasonal pattern that repeats every four quarters and you see the upward trend-cycle.
(b) Classical multiplicative decomposition
# model() fits a "decomposition model" and components() extracts the pieces.
# type = "multiplicative" means Gas = trend * seasonal * random.
gas_dcmp <- gas |>
model(classical_decomposition(Gas, type = "multiplicative")) |>
components()
gas_dcmp # trend cycle, seasonal indices, random all seasonally adjusted## # A dable: 20 x 7 [1Q]
## # Key: .model [1]
## # : Gas = trend * seasonal * random
## .model Quarter Gas trend seasonal random season_adjust
## <chr> <qtr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "classical_decomposition(G… 2005 Q3 221 NA 1.13 NA 196.
## 2 "classical_decomposition(G… 2005 Q4 180 NA 0.925 NA 195.
## 3 "classical_decomposition(G… 2006 Q1 171 200. 0.875 0.974 195.
## 4 "classical_decomposition(G… 2006 Q2 224 204. 1.07 1.02 209.
## 5 "classical_decomposition(G… 2006 Q3 233 207 1.13 1.000 207.
## 6 "classical_decomposition(G… 2006 Q4 192 210. 0.925 0.987 208.
## 7 "classical_decomposition(G… 2007 Q1 187 213 0.875 1.00 214.
## 8 "classical_decomposition(G… 2007 Q2 234 216. 1.07 1.01 218.
## 9 "classical_decomposition(G… 2007 Q3 245 219. 1.13 0.996 218.
## 10 "classical_decomposition(G… 2007 Q4 205 219. 0.925 1.01 222.
## 11 "classical_decomposition(G… 2008 Q1 194 219. 0.875 1.01 222.
## 12 "classical_decomposition(G… 2008 Q2 229 219 1.07 0.974 213.
## 13 "classical_decomposition(G… 2008 Q3 249 219 1.13 1.01 221.
## 14 "classical_decomposition(G… 2008 Q4 203 220. 0.925 0.996 219.
## 15 "classical_decomposition(G… 2009 Q1 196 222. 0.875 1.01 224.
## 16 "classical_decomposition(G… 2009 Q2 238 223. 1.07 0.993 222.
## 17 "classical_decomposition(G… 2009 Q3 252 225. 1.13 0.994 224.
## 18 "classical_decomposition(G… 2009 Q4 210 226 0.925 1.00 227.
## 19 "classical_decomposition(G… 2010 Q1 205 NA 0.875 NA 234.
## 20 "classical_decomposition(G… 2010 Q2 236 NA 1.07 NA 220.# The four seasonal indices one per quarter. They multiply trend, so
# 1.13 means "Q3 is typically 13% above the trend-cycle".
gas_dcmp |>
as_tibble() |>
mutate(Qtr = quarter(Quarter)) |>
distinct(Qtr, seasonal)## # A tibble: 4 × 2
## Qtr seasonal
## <int> <dbl>
## 1 3 1.13
## 2 4 0.925
## 3 1 0.875
## 4 2 1.07autoplot(gas_dcmp) +
labs(title = "Classical multiplicative decomposition of Gas (last 5 years)")
(c) Do the results support the graphical interpretation from part (a)?
Answer. Yes. you the trend clearly going upwards. For seasonal you indeed to see similar seasonal fluctuations and with random I definitely see that trend * seasonal explains almost all of the movement.
(d) Compute and plot the seasonally adjusted data
# for multiplicative decomposition, a seosonally adjusted series is is y/s to divide out seasonal index
gas_dcmp |>
mutate(by_hand = Gas / seasonal) |> # same thing as season_adjust
as_tibble() |>
select(Quarter, Gas, seasonal, season_adjust, by_hand) |>
head(6)## # A tibble: 6 × 5
## Quarter Gas seasonal season_adjust by_hand
## <qtr> <dbl> <dbl> <dbl> <dbl>
## 1 2005 Q3 221 1.13 196. 196.
## 2 2005 Q4 180 0.925 195. 195.
## 3 2006 Q1 171 0.875 195. 195.
## 4 2006 Q2 224 1.07 209. 209.
## 5 2006 Q3 233 1.13 207. 207.
## 6 2006 Q4 192 0.925 208. 208.gas_dcmp |>
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Gas, colour = "Original")) +
geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
geom_line(aes(y = trend, colour = "Trend-cycle")) +
scale_colour_manual(values = c("Original" = "grey60",
"Seasonally adjusted" = "#0072B2",
"Trend-cycle" = "#D55E00")) +
labs(title = "Gas production: original vs. seasonally adjusted",
y = "Petajoules", x = "Quarter", colour = NULL)
The seasonally adjusted series removes the consistent up and down fluctuations and leaves a slowly rising line with a little bit of noise around it.
(e) Add an outlier and recompute
We add 300 to one observation in the middle of the series (2007 Q4, the 10th of 20 quarters) and redo everything.
# I'm adding 300 to an observation in the middle of the series which is after the 10th observation
# if_else() replaces just the one quarter while everything else stays the same
gas_mid <- gas |>
mutate(Gas = if_else(Quarter == yearquarter("2007 Q4"), Gas + 300, Gas))
gas_mid_dcmp <- gas_mid |>
model(classical_decomposition(Gas, type = "multiplicative")) |>
components()
autoplot(gas_mid_dcmp) +
labs(title = "Decomposition with an outlier in the middle")
# Compare the seasonal indices before and after the outlier
bind_rows(
gas_dcmp |> as_tibble() |> mutate(Series = "Original"),
gas_mid_dcmp |> as_tibble() |> mutate(Series = "Outlier in middle")
) |>
mutate(Qtr = quarter(Quarter)) |>
distinct(Series, Qtr, seasonal) |>
pivot_wider(names_from = Series, values_from = seasonal)## # A tibble: 4 × 3
## Qtr Original `Outlier in middle`
## <int> <dbl> <dbl>
## 1 3 1.13 1.06
## 2 4 0.925 1.12
## 3 1 0.875 0.821
## 4 2 1.07 0.998# Seasonally adjusted series with original vs. with the outlier
bind_rows(
gas_dcmp |> as_tibble() |> mutate(Series = "Original"),
gas_mid_dcmp |> as_tibble() |> mutate(Series = "Outlier in middle")
) |>
ggplot(aes(x = Quarter, y = season_adjust, colour = Series)) +
geom_line() +
labs(title = "Seasonally adjusted Gas: effect of an outlier in the middle",
y = "Petajoules", x = "Quarter", colour = NULL)
(f) Does it matter if the outlier is near the end instead of the middle?
gas_end <- gas |>
mutate(Gas = if_else(Quarter == yearquarter("2010 Q1"), Gas + 300, Gas)) # 19th of 20
gas_end_dcmp <- gas_end |>
model(classical_decomposition(Gas, type = "multiplicative")) |>
components()
autoplot(gas_end_dcmp) +
labs(title = "Decomposition with an outlier near the end (2010 Q1 + 300)")
# Seasonal indices - original vs middle-outlier vs end-outlier
bind_rows(
gas_dcmp |> as_tibble() |> mutate(Series = "Original"),
gas_mid_dcmp |> as_tibble() |> mutate(Series = "Outlier in middle"),
gas_end_dcmp |> as_tibble() |> mutate(Series = "Outlier near end")
) |>
mutate(Qtr = quarter(Quarter)) |>
distinct(Series, Qtr, seasonal) |>
pivot_wider(names_from = Series, values_from = seasonal)## # A tibble: 4 × 4
## Qtr Original `Outlier in middle` `Outlier near end`
## <int> <dbl> <dbl> <dbl>
## 1 3 1.13 1.06 1.11
## 2 4 0.925 1.12 0.889
## 3 1 0.875 0.821 0.897
## 4 2 1.07 0.998 1.10# All three seasonally adjusted series on one plot
bind_rows(
gas_dcmp |> as_tibble() |> mutate(Series = "Original"),
gas_mid_dcmp |> as_tibble() |> mutate(Series = "Outlier in middle"),
gas_end_dcmp |> as_tibble() |> mutate(Series = "Outlier near end")
) |>
ggplot(aes(x = Quarter, y = season_adjust, colour = Series)) +
geom_line() +
facet_wrap(~ Series, ncol = 1) +
theme(legend.position = "none") +
labs(title = "Seasonally adjusted Gas: outlier in the middle vs. near the end",
y = "Petajoules", x = "Quarter")
Answer. Yes it makes a diffirence because trend cycle is derrives from moving averages and this case since were dealing with quarters you get a 2 * 4-MA which leaves the first two and last two values as NA so an outlier in the middle dissort the trend and seasonal indices across the entire series while an outlier at end doesn’t necesarily interupts with the decomposition since it’s a place where decomposition can’t estimate it.
Exercise 3.8 — X-11 decomposition of my retail series
Recall your retail time series data. Decompose the series using X-11. Does it reveal any outliers, or unusual features that you had not noticed previously?
X_13ARIMA_SEATS() from feasts calls the Census Bureau’s
X-13 program through the seasonal package. Turnover ~ x11()
asks for the X-11 variant, which doesn’t function like classical
decomposition, the estimates trends all the way to the ends and lets the
seasonal pattern change slowly over time, and is robust to outliers
x11_dcmp <- myseries |>
model(x11 = X_13ARIMA_SEATS(Turnover ~ x11())) |>
components()
autoplot(x11_dcmp) +
labs(title = "X-11 decomposition: Northern Territory clothing & footwear turnover")
A second view that has the data with the trend and the seasonally adjusted series overlaid on top, which makes the outliers in the adjusted series easier to see.
x11_dcmp |>
ggplot(aes(x = Month)) +
geom_line(aes(y = Turnover, colour = "Data")) +
geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
geom_line(aes(y = trend, colour = "Trend")) +
scale_colour_manual(values = c("Data" = "grey60",
"Seasonally adjusted" = "#0072B2",
"Trend" = "#D55E00")) +
labs(title = "Retail turnover with X-11 trend and seasonally adjusted series",
y = "$ million AUD", x = "Month", colour = NULL)
Because X-11 is multiplicative here, the irregular component is a ratio that’s around 1. We flag months where the irregular is more than 2.5 from the standard deviations and it’s its mean as candidate outliers.
x11_dcmp |>
as_tibble() |>
mutate(z = (irregular - mean(irregular)) / sd(irregular)) |> # standardise
filter(abs(z) > 2.5) |> # unusually large/small
select(Month, Turnover, trend, seasonal, irregular, z) |>
arrange(Month)## # A tibble: 13 × 6
## Month Turnover trend seasonal irregular z
## <mth> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1988 May 2.9 2.63 0.958 1.15 3.92
## 2 1988 Jul 2.8 2.67 1.20 0.874 -3.35
## 3 1990 Feb 2.8 4.20 0.752 0.887 -3.02
## 4 1992 Mar 4.1 5.86 0.825 0.849 -4.03
## 5 1993 Dec 7.9 6.52 1.38 0.881 -3.18
## 6 1994 Feb 5.2 6.15 0.741 1.14 3.72
## 7 1997 Nov 6.4 5.61 1.03 1.11 2.79
## 8 1998 Dec 7.5 6.40 1.37 0.857 -3.81
## 9 1999 Jan 5.9 6.44 0.815 1.12 3.27
## 10 2000 Jun 8.5 7.46 0.997 1.14 3.74
## 11 2010 Jul 16.1 12.5 1.14 1.13 3.40
## 12 2013 Jan 11.8 12.7 0.821 1.14 3.56
## 13 2016 Feb 12.3 14.3 0.747 1.15 3.95# How the seasonal pattern itself has drifted over 30 years
x11_dcmp |>
gg_subseries(seasonal) +
labs(title = "X-11 seasonal component by month (note the changing December effect)",
y = "Seasonal index")
Answer. Yes, X-11 reveals a couple of things the plot hides, the irregular component showed us a few outliers that were invisible in the seasonal. It also hides that the trend was not smooth, there was some bumps, stair step growth and no growth since 2016 that I had not noticed before. Also in seasonal, classical decomposition would be completely wrong as December’s rises and falls are inconsistent so a fixed seasonal pattern hides such details.
Exercise 3.9 — Australian civilian labour force decomposition
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.
To read the scales with more care I had to reproduce the figures from
the same data (fma::labour, the dataset the book uses
fma which is on CRAN and accesible by
install.packages("fma")), using STL with a 13 month
seasonal window.
# fma::labour is the classic 'ts' object, as_tsibble() converts it to a tsibble.
labour <- fma::labour |>
as_tsibble() |>
rename(Month = index, Persons = value) # nicer and neater names than index/value
labour_dcmp <- labour |>
model(STL(Persons ~ season(window = 13))) |>
components()
autoplot(labour_dcmp) +
labs(title = "STL decomposition - Australian civilian labour force (thousands)")
labour_dcmp |>
gg_subseries(season_year) +
labs(title = "Seasonal component by month", y = "Seasonal effect (thousands)")
# Numbers behind the description below
labour_dcmp |>
as_tibble() |>
summarise(data_min = min(Persons), data_max = max(Persons),
seas_min = min(season_year), seas_max = max(season_year),
rem_min = min(remainder), rem_max = max(remainder))## # A tibble: 1 × 6
## data_min data_max seas_min seas_max rem_min rem_max
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 6404. 9003. -93.9 122. -260. 147.# The largest remainders — where do they fall?
labour_dcmp |>
as_tibble() |>
slice_min(remainder, n = 5) |>
select(Month, Persons, trend, remainder)## # A tibble: 5 × 4
## Month Persons trend remainder
## <mth> <dbl> <dbl> <dbl>
## 1 1991 Mar 8199. 8430. -260.
## 2 1991 Jun 8199. 8448. -207.
## 3 1991 Apr 8323. 8432. -106.
## 4 1991 Jan 8263 8431. -82.2
## 5 1993 Apr 8520. 8607. -79.9# How much did the trend grow each year around the recession?
labour_dcmp |>
as_tibble() |>
mutate(Year = year(Month)) |>
filter(Year %in% 1989:1994) |>
group_by(Year) |>
summarise(trend_growth = last(trend) - first(trend)) # Dec trend minus Jan trend## # A tibble: 6 × 2
## Year trend_growth
## <dbl> <dbl>
## 1 1989 242.
## 2 1990 70.4
## 3 1991 125.
## 4 1992 36.9
## 5 1993 84.8
## 6 1994 171.(a) Describe the results (3–5 sentences, paying attention to the scales)
The decomposition is mainly influence by the trend, during observation I see the labour force grew from 6.4 M to 9.0 M persons from 1978 - 1995, basically a rise ~ 2.6 M. The trend panel looks almost exactly like the data panel but that’s mainly due to the other components being smaller. If we look at the seasonal component, it only ranges from -100 to 120 with the labor force highest in December and lowest during January and August. The gray bars we see in the left of each panel makes the scale difference visible from a glance. The remainder is usually small and looks like noise and is clearly reflected as we look at 1991, it affects the data panel slightly but since is small but not. much.
(b) Is the 1991/1992 recession visible in the estimated components?
Yes but mainly in the remainder as discussed at the end of part a of this question. Those large negative falls we see in remainder perfectly reflects this, this is the most obvious however it is also slightly observable in the trend panel as you see a small dip during that same time.
Session info
# Package versions for reproducibility
sessionInfo()## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Sequoia 15.7.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/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] fable_0.5.0 feasts_0.5.0 fabletools_0.8.0 ggtime_1.0.0
## [5] tsibbledata_0.4.1 tsibble_1.2.0 ggplot2_4.0.3 lubridate_1.9.5
## [9] tidyr_1.3.2 dplyr_1.2.1 tibble_3.3.1 fpp3_1.0.3
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.6 rappdirs_0.3.4 sass_0.4.10
## [4] generics_0.1.4 anytime_0.3.13 lattice_0.22-9
## [7] digest_0.6.39 magrittr_2.0.5 evaluate_1.0.5
## [10] grid_4.6.1 timechange_0.4.0 RColorBrewer_1.1-3
## [13] mixtime_0.3.0 fastmap_1.2.0 jsonlite_2.0.0
## [16] seasonal_1.11.0 forecast_9.0.2 purrr_1.2.2
## [19] scales_1.4.0 jquerylib_0.1.4 cli_3.6.6
## [22] x13binary_1.1.61.2 rlang_1.3.0 crayon_1.5.3
## [25] vecvec_1.3.0 withr_3.0.3 cachem_1.1.0
## [28] yaml_2.3.12 otel_0.2.0 parallel_4.6.1
## [31] tools_4.6.1 tzdb_0.5.0 colorspace_2.1-3
## [34] vctrs_0.7.3 R6_2.6.1 zoo_1.9-0
## [37] fma_2.5 lifecycle_1.0.5 urca_1.3-4
## [40] pkgconfig_2.0.3 progressr_1.0.0 pillar_1.11.1
## [43] bslib_0.12.0 gtable_0.3.6 glue_1.8.1
## [46] Rcpp_1.1.2 xfun_0.60 tidyselect_1.2.1
## [49] rstudioapi_0.19.0 knitr_1.51 farver_2.1.2
## [52] nlme_3.1-169 htmltools_0.5.9 labeling_0.4.3
## [55] rmarkdown_2.32 timeDate_4052.112 fracdiff_1.5-4
## [58] compiler_4.6.1 S7_0.2.2 distributional_0.8.1