# Question 1
library(tsibble)

## 
## Attaching package: 'tsibble'

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union

library(tsibbledata)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(feasts)

## Loading required package: fabletools

# Question 2
library(fpp3)

## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.3 ──

## ✔ tibble    3.3.1     ✔ ggtime    1.0.0
## ✔ tidyr     1.3.2     ✔ fable     0.5.0
## ✔ lubridate 1.9.5

## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()     masks base::date()
## ✖ dplyr::filter()       masks stats::filter()
## ✖ tsibble::intersect()  masks base::intersect()
## ✖ lubridate::interval() masks tsibble::interval()
## ✖ dplyr::lag()          masks stats::lag()
## ✖ tsibble::setdiff()    masks base::setdiff()
## ✖ tsibble::union()      masks base::union()

3.1 Exercise

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?

1. Plot GDP per capita for each country over time

global_economy %>%
  autoplot(GDP / Population, show.legend = FALSE) +
  labs(
    title = "GDP Per Capita for Each Country Over Time",
    y = "$US per capita",
    x = "Year"
  )

2. Find the country and year with the highest GDP per capita

highest_gdp_pc <- global_economy %>%
  mutate(GDP_per_capita = GDP / Population) %>%
  arrange(desc(GDP_per_capita))

Printing the top record

head(highest_gdp_pc, 1)

## # A tsibble: 1 x 10 [1Y]
## # Key:       Country [1]
##   Country Code   Year         GDP Growth   CPI Imports Exports Population
##   <fct>   <fct> <dbl>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>      <dbl>
## 1 Monaco  MCO    2014 7060236168.   7.18    NA      NA      NA      38132
## # ℹ 1 more variable: GDP_per_capita <dbl>

Analysis of Results for Exercise 3.1

Highest GDP per capita: Monaco (or Liechtenstein/Luxembourg depending on the specific recorded year) holds the highest peak GDP per capita in the dataset.
Changes over time: Overall, global GDP per capita shows a consistent upward trend for most nations, punctuated by temporary economic contractions (such as the 2009 global financial crisis). Small nations with high incomes and low populations frequently dominate the absolute highest per capita rankings.

3.2 Exercise

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.
Slaughter of Victorian “Bulls, bullocks and steers” in aus_livestock.
Victorian Electricity Demand from vic_elec.
Gas production from aus_production.

United States GDP from global_economy.

us_gdp <- global_economy %>% 
  filter(Country == "United States")

us_gdp %>% autoplot(GDP) + 
  labs(title = "United States GDP", y = "$US")

# Optional Box-Cox / Transformation check for US GDP
us_gdp %>% features(GDP, features = guerrero) # Finds optimal lambda

## # A tibble: 1 × 2
##   Country       lambda_guerrero
##   <fct>                   <dbl>
## 1 United States           0.282

us_gdp %>% autoplot(box_cox(GDP, 0.282)) + 
  labs(title = "Transformed US GDP (Box-Cox)")

Slaughter of Victorian “Bulls, bullocks and steers” in aus_livestock.

vic_bulls <- aus_livestock %>% 
  filter(State == "Victoria", Animal == "Bulls, bullocks and steers")

vic_bulls %>% autoplot(Count) + 
  labs(title = "Slaughter of Victorian Bulls, Bullocks and Steers", y = "Count")

Victorian Electricity Demand from vic_elec.

vic_elec %>% autoplot(Demand) + 
  labs(title = "Victorian Electricity Demand (Half-hourly)", y = "Demand")

# Aggregating to daily/monthly is often preferred since half-hourly is too dense
vic_elec %>% 
  mutate(Date = as_date(Time)) %>% 
  group_by(Date) %>% 
  summarise(Demand = sum(Demand)) %>% 
  autoplot(Demand) + 
  labs(title = "Daily Victorian Electricity Demand", y = "Demand (MWh)")

Gas production from aus_production.

aus_production %>% autoplot(Gas) + 
  labs(title = "Australian Gas Production", y = "Petajoules")

# Box-Cox transformation using Guerrero's optimal lambda
lambda_gas <- aus_production %>% 
  features(Gas, features = guerrero) %>% 
  pull(lambda_guerrero)

aus_production %>% autoplot(box_cox(Gas, lambda_gas)) + 
  labs(title = "Transformed Australian Gas Production (Box-Cox)", subtitle = paste("Lambda =", round(lambda_gas, 2)))

Effects of Transformations for Exercise 3.2

United States GDP: Shows a strong upward exponential trend where absolute variability increases over time. A Box-Cox transformation (with \(\lambda \approx 0.28\)) or a log transformation helps stabilize the increasing variance and linearizes the multiplicative growth trend.
Victorian Bulls, Bullocks and Steers: Displays stable variance or additive patterns; transformations are generally unnecessary unless specific seasonal variance stabilization is required.
Victorian Electricity Demand: The raw series is recorded half-hourly, making the standard graph a solid block of ink. Transforming the temporal scale by aggregating to daily or monthly totals makes the underlying seasonal and calendar patterns readable.
Australian Gas Production: Seasonal fluctuations grow larger as the level of production increases over time (multiplicative seasonality). A Box-Cox transformation (using Guerrero’s optimal \(\lambda \approx 0.11\)) successfully stabilizes the seasonal variation so that swings are uniform across the entire historical period.

3.3 Exercise

The canadian_gas dataset from the fpp3 package consists of monthly Canadian gas production data.
A Box-Cox transformation is designed to fix a specific issue: heteroscedasticity, which occurs when the size of the seasonal fluctuations or variance grows or shrinks consistently alongside the overall level/trend of the time series.
For the canadian_gas data, a Box-Cox transformation is unhelpful because the variance does not change consistently over time. Instead, the absolute seasonal variation starts small, wildly increases in the middle portion of the dataset (the 1970s and 1980s), and then drops significantly and stabilizes later on—even as the overall production volume remains high. Because the variance behavior fluctuates independently of the overall data level, no single mathematical exponent (\[\lambda \]) can stabilize it across the entire timeline.

library(fpp3)
library(patchwork) # For aligning plots side-by-side.

Plotting the raw data.

p1 <- canadian_gas %>% 
  autoplot(Volume) +
  labs(
    title = "Original Canadian Gas Production",
    y = "Volume (Billion cubic metres)",
    x = "Month"
  )

Finding the optimal Guerrero lambda parameter.

lambda_val <- canadian_gas %>% 
  features(Volume, features = guerrero) %>% 
  pull(lambda_guerrero)

Plot the Box-Cox transformed data.

p2 <- canadian_gas %>% 
  autoplot(box_cox(Volume, lambda_val)) +
  labs(
    title = "Box-Cox Transformed Data",
    subtitle = paste("Using optimal Guerrero lambda =", round(lambda_val, 2)),
    y = "Transformed Volume",
    x = "Month"
  )

# Combine the plots to visibly inspect the variance
p1 / p2

Visual Breakdown of the Results for Exercise 3.3

The Core Issue: In the original plot, the variance (the width of the seasonal swings) is low up to the mid-1970s, spikes drastically from the late 1970s to the 1990s, and then narrows significantly after 2000.
The Failed Fix: In the transformed plot, the middle section is still compressed awkwardly while the end sections are stretched. Because the variation doesn’t increase strictly as the trend goes up, a uniform power transformation cannot “flatten” the unequal volatility.
Instead of a standard Box-Cox transformation, handling the complex patterns in canadian_gas requires models that deal with changing structural regimes or a seasonal decomposition method like STL decomposition to pull out the shifting seasonal component dynamically.

3.4 Exercise

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

Because the aus_retail dataset contains many unique time series (each with a different Series ID), the specific Box-Cox transformation parameter (\[\lambda \]) depends entirely on the random seed or specific retail series you selected in your earlier exercise.
The best R code to address this question uses Guerrero’s method via the features() function to dynamically find the mathematically optimal \[\lambda \] for your exact series, and then visualizes the transformation side-by-side with the original data using patchwork.

library(fpp3)
library(patchwork) # Used again to display plots side-by-side.

# 1. Use the EXACT same random seed and subset from your Section 2.10 Exercise
set.seed(624) # Replace with the seed you used in Exercise 7
myseries <- aus_retail %>% 
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# 2. Automatically find the optimal Guerrero lambda 
optimal_lambda <- myseries %>% 
  features(Turnover, features = guerrero) %>% 
  pull(lambda_guerrero)

# 3. Create the original data plot
p1 <- myseries %>% 
  autoplot(Turnover) + 
  labs(title = "Original Retail Turnover", y = "Turnover")

# 4. Create the Box-Cox transformed plot using the optimal lambda
p2 <- myseries %>% 
  autoplot(box_cox(Turnover, optimal_lambda)) + 
  labs(
    title = "Box-Cox Transformed Data", 
    subtitle = paste("Guerrero Lambda =", round(optimal_lambda, 3)),
    y = "Transformed Turnover"
  )

# 5. Display plots stacked vertically to evaluate variance stabilization
p1 / p2

With an optimal Guerrero lambda of 0.002, the transformation is indeed very close to a log transformation (\(\lambda = 0\)). This value indicates that your specific retail series has a strong multiplicative seasonal component that expands rapidly as the overall turnover grows.

3.5 Exercise

For the following series, find an appropriate Box-Cox transformation in order to stabilise the variance. Tobacco from aus_production, Economy class passengers between Melbourne and Sydney from ansett, and Pedestrian counts at Southern Cross Station from pedestrian.

Here is my workflow: each series requires filtering before analysis, with the pedestrian data needing extra scrutiny to handle an edge case upfront. I then use the Guerrero method to identify an optimal \(\lambda\), always sanity-checking the resulting transformation instead of blindly relying on the automated value.

Tobacco Production

lambda_tobacco <- aus_production |>
  features(Tobacco, features = guerrero) |>
  pull(lambda_guerrero)

lambda_tobacco

## [1] 0.9264636

aus_production |>
  filter(!is.na(Tobacco)) |>
  autoplot(Tobacco) +
  labs(title = "Tobacco Production before Transformation",
       x = "Quarter", y = "Tobacco (tonnes)")

aus_production |>
  filter(!is.na(Tobacco)) |>
  autoplot(box_cox(Tobacco, lambda_tobacco)) +
  labs(title = "Tobacco Production after Box-Cox Transformation",
       x = "Quarter",
       y = paste0("Transformed tobacco (lambda = ", round(lambda_tobacco, 3), ")"))

The Guerrero method yields an optimal parameter of \(\lambda \approx 0.93\). Because a value this close to 1 corresponds to a simple linear shift rather than a meaningful power transformation, the algorithm indicates that the series requires no stabilization. This is corroborated by comparing the pre- and post-transformation plots, which show virtually identical profiles due to the data’s already stable variance. Rather than applying an unnecessary transformation, I retain the series in its original scale; a \(\lambda \approx 1\) confirms that the variance is well-behaved as is.

Melbourne–Sydney Economy Route Analysis

ansett_econ <- ansett |>
  filter(Class == "Economy", Airports == "MEL-SYD")

ansett_econ |>
  autoplot(Passengers) +
  labs(title = "Melbourne–Sydney Economy Route Analysis",
       x = "Week", y = "Passengers")

lambda_ansett <- ansett_econ |>
  features(Passengers, features = guerrero) |>
  pull(lambda_guerrero)

lambda_ansett

## [1] 1.999927

ansett_econ |>
  autoplot(box_cox(Passengers, lambda_ansett)) +
  labs(title = "Melbourne–Sydney Economy Route after Box-Cox Transformation",
       x = "Week",
       y = paste0("Transformed Passengers (lambda = ", round(lambda_ansett, 3), ")"))

The Guerrero method yields an optimal parameter of \(\lambda \approx 2\), which squares the series rather than compressing it. This is an counterintuitive result that merits closer inspection: values of \(\lambda > 1\) amplify large values relative to small ones, running counter to the typical purpose of a variance-stabilizing transformation. The root cause is a prolonged collapse to zero passenger traffic in 1989 caused by an industrial dispute. This extended run of zeros heavily distorts the optimization, causing Guerrero to react to an external operational shock rather than the series’ underlying seasonal or structural variance. Consequently, a Box-Cox transformation is ill-suited here; this structural anomaly requires targeted intervention—such as intervention modeling or filtering—rather than reliance on a distorted parameter estimate.

Southern Cross Station Pedestrian Traffic Analysis

ped <- pedestrian |>
  filter(Sensor == "Southern Cross Station")

# There are two issues for examining prior to transforming. Zero-count hours.
sum(ped$Count == 0)

## [1] 159

# Does the time series contain temporal gaps or missing hourly timestamps?
nrow(has_gaps(ped) |> filter(.gaps))

## [1] 1

The hourly series presents two distinct data-quality issues, either of which compromises downstream transformations. First, it includes 159 observations of exactly zero—primarily during overnight periods of inactivity. Second, the series exhibits implicit temporal gaps where specific hourly timestamps are omitted entirely. Because power transformations such as the Box-Cox family are undefined at zero when \(\lambda \le 0\), attempting a log-like transformation yields infinite values. Aggregating the observations into daily counts resolves both challenges concurrently: daily foot traffic is strictly positive, eliminating the zero-count issue, while temporal aggregation constructs a continuous, regular time index.

ped_daily <- ped |>
  index_by(Date) |>
  summarise(Count = sum(Count))

sum(ped_daily$Count == 0)

## [1] 0

nrow(has_gaps(ped_daily) |> filter(.gaps))

## [1] 0

ped_daily |>
  autoplot(Count) +
  labs(title = "Southern Cross Station Pedestrian Traffic Analysis",
       x = "Date", y = "Count")

lambda_ped <- ped_daily |>
  features(Count, features = guerrero) |>
  pull(lambda_guerrero)

lambda_ped

## [1] 0.2726316

ped_daily |>
  autoplot(box_cox(Count, lambda_ped)) +
  labs(title = "Southern Cross Station Pedestrian post Box-Cox Transformation",
       x = "Date",
       y = paste0("Transformed count (lambda = ", round(lambda_ped, 3), ")"))

Applied to the aggregated daily series, the Guerrero method identifies an optimal parameter of \(\lambda \approx 0.27\), falling between a logarithmic and a square-root transformation. This transformation compresses the scale, pulling lower-traffic days closer to peak periods and achieving noticeably more uniform variance throughout the series. The broader takeaway is that exploratory data screening is far more critical than automated parameter estimation: running a Box-Cox transformation directly on the raw hourly data would have silently failed or produced invalid values due to the unhandled zero counts.

3.7 Exercise

Consider the last five years of the Gas data from aus_production. gas <- tail(aus_production, 5*4) |> select(Gas)

Plot the time series. Can you identify seasonal fluctuations and/or a trend-cycle?
Use classical_decomposition with type=multiplicative to calculate the trend-cycle and seasonal indices.
Do the results support the graphical interpretation from part a?
Compute and plot the seasonally adjusted data.
Compute and plot the seasonally adjusted data.
Does it make any difference if the outlier is near the end rather than in the middle of the time series?

library(fpp3)
library(patchwork)

# Filtering the last 5 years (5 years * 4 quarters = 20 observations).

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


# a. Rendering the time series.

p_raw <- gas |> autoplot(Gas) +
  labs(title = "Last 5 Years of Australian Gas Production", y = "Petajoules")
print(p_raw)

# b. Classical Multiplicative Decomposition.

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

# Extract and plot components.
components(classical_fit) |> autoplot() +
  labs(title = "Classical Multiplicative Decomposition")

# d. Seasonally Adjusted Data.

p_adj <- components(classical_fit) |>
  autoplot(season_adjust) +
  labs(title = "Seasonally Adjusted Gas Data (Original)", y = "Petajoules")
print(p_adj)

# e. Introducing an Outlier in the Middle.

gas_mid_outlier <- gas
# Row 10 represents a point roughly in the middle of the 20 observations.
gas_mid_outlier$Gas[10] <- gas_mid_outlier$Gas[10] + 300

p_mid <- gas_mid_outlier |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components() |>
  autoplot(season_adjust) +
  labs(title = "Seasonally Adjusted (Outlier in Middle)")


# f. Introducing an Outlier at the End.

gas_end_outlier <- gas
# Row 20 represents the final observation.
gas_end_outlier$Gas[20] <- gas_end_outlier$Gas[20] + 300

p_end <- gas_end_outlier |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components() |>
  autoplot(season_adjust) +
  labs(title = "Seasonally Adjusted (Outlier at End)")

# Comparing outlier positions side-by-side.
p_mid / p_end

Exercise Answers Point-by-Point

Identification: The raw time series plot shows a clear upward trend-cycle because production levels increase over the 5-year span. Strong seasonal fluctuations are visible: production consistently peaks during Q3 (winter months with high heating demand) and bottoms out during Q1.
Classical Decomposition: The code utilizes classical_decomposition(..., type = "multiplicative") to break the series into trend, seasonal, and random components.
Comparison: Yes, the results fully support the interpretation from part a. The seasonal panel shows a rigid, recurring pattern peaking in Q3, while the trend panel isolates a smooth, continuous upward slope.
Seasonally Adjusted Data: Plotting season_adjust reveals a much smoother line. By removing the predictable quarterly spikes, you can track the underlying economic trajectory directly.
Effect of a Mid-Series Outlier: Because classical decomposition relies heavily on symmetric moving averages to compute the trend, adding a massive value spreads the distortion across nearby dates. It creates an artificial bump in the surrounding trend values and a sharp spike in the season_adjust line at that specific index.
End-of-Series Outlier vs. Middle: Yes, the location makes a significant difference.

When an outlier is in the middle, the moving average formulas smooth it partially into the trend of surrounding periods.
When the outlier is at the very end, the symmetric moving average cannot be fully computed due to missing future data points. As a result, the outlier aggressively forces the end of the calculated trend line upward or downward, heavily contaminating final-quarter indicators.

3.8 Exercise

Since aus_retail contains multiple time series, the code utilizes a random seed to draw a specific series (reproducing the structure of Section 2.10, Exercise 7), extracts the X-11 components, and plots them.

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

x11_decomp |>
  autoplot() +
  labs(title = "X-11 Decomposition of Retail Turnover")

x11_decomp |>
  ggplot(aes(x = Month)) +
  geom_line(aes(y = Turnover, colour = "Original")) +
  geom_line(aes(y = season_adjust, colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(title = "Retail Turnover, Original, Seasonally Adjusted, and Trend",
       x = "Month", y = "Turnover ($ million)", colour = "Series")

irr_mean <- mean(x11_decomp$irregular, na.rm = TRUE)
irr_sd   <- sd(x11_decomp$irregular, na.rm = TRUE)

x11_decomp |>
  as_tibble() |>
  filter(abs(irregular - irr_mean) > 3 * irr_sd) |>
  select(Month, Turnover, irregular) |>
  arrange(Month)

## # A tibble: 10 × 3
##       Month Turnover irregular
##       <mth>    <dbl>     <dbl>
##  1 1983 Oct     80.1     0.917
##  2 1983 Dec     89.6     0.914
##  3 1986 Dec    136.      1.07 
##  4 1990 Sep    135.      0.914
##  5 1993 Mar    122.      0.893
##  6 1994 Mar    171.      1.07 
##  7 1999 Jul    152.      0.930
##  8 1999 Nov    165.      1.07 
##  9 2002 Jan    229.      1.09 
## 10 2003 Mar    202.      0.930

The X-11 decomposition identifies ten months where the irregular component deviates by more than three standard deviations from its mean, several of which highlight clear structural phenomena. Notably, March 1993 and March 1994 exhibit sharp opposing extremes—the former unusually depressed and the latter substantially elevated. This swing illustrates the classic ‘Easter effect’: because Easter alternates between March and April, holiday-driven demand shifts across calendar boundaries. Since standard monthly seasonal filters assume static, fixed periodic patterns, these calendar-shift dynamics spill directly into the irregular residual. Additionally, January 2002 registers the series’ largest positive anomaly, while October and December 1983 display pronounced negative outliers near the beginning of the series, likely reflecting early-stage reporting noise or lower historical data reliability.

Building upon the exploratory analysis in Homework 1—which correctly identified a dominant upward trend, a pronounced December seasonal peak, and a lack of cyclicality—this decomposition provides a more granular perspective. In the raw time plot, the sheer magnitude of the trend and seasonality masks localized anomalies. By stripping away these dominant structural components, decomposition isolates the irregular residual. This process allows subtle, single-month deviations to emerge—most notably the moving Easter effect, which was entirely obscured in the original data but becomes clearly identifiable in the decomposed residual.

3.9 Exercise

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.

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 decomposition is heavily dominated by the secular trend, as the total labour force expands steadily from roughly 6.4 million to 9.0 million individuals over the seventeen-year sample period. Consequently, the trend component exhibits the widest dynamic range across all panels. By contrast, the seasonal variation is modest, spanning approximately \(\pm 100{,}000\) persons—or roughly one percent of the aggregate series. Interpreting the vertical scales is critical here: while the seasonal plot resembles a pronounced oscillation, this visual amplitude is an artifact of vertical axis magnification. The gray scale bars on the right margins confirm this disparity; despite having similar physical heights, they denote value ranges differing by two orders of magnitude. Furthermore, the seasonal profile exhibits slow structural evolution rather than strict periodicity—notably, the March peak climbs from 65 to nearly 90 thousand in the mid-1980s before receding to roughly 40 thousand by 1995. This time-varying behavior highlights why flexible filtering via STL is far superior to classical decomposition with fixed seasonal components. Finally, the remainder is largely white noise, punctuated only by a sharp negative shock in the early 1990s that sits far outside its typical distribution.

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

The recessionary shock is distinctly visible across two of the three decomposed components. First, the trend series shows a clear structural deceleration, flattening substantially between 1990 and 1992 before resuming its long-term upward trajectory. Second—and more prominently—the remainder component exhibits a steep double-dip trough reaching approximately \(-380\), markedly exceeding its standard stationary noise band of \(\pm 100\). This negative deviation signifies that realized labor force numbers dropped sharply below what the combined trend and seasonal estimates anticipated. Conversely, the seasonal component remains unaffected throughout this period. Because an economic downturn represents a macroeconomic shock rather than a calendar-driven phenomenon, it disrupts the underlying trend and inflates irregular residuals without altering baseline month-to-month employment seasonality.

Puerto Rico

Homework 02 Puerto Rico

Dariusz

2026-09-20