library(tsibble)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(fpp3)
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──
## ✔ tibble      3.3.0          ✔ ggplot2     4.0.0.9000
## ✔ dplyr       1.1.4          ✔ tsibbledata 0.4.1     
## ✔ tidyr       1.3.1          ✔ feasts      0.4.2     
## ✔ lubridate   1.9.4          ✔ fable       0.5.0
## Warning: package 'fable' was built under R version 4.5.2
## ── 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()
library(USgas)
library(dplyr)
library(ggplot2)
library(readxl)

#Question 3.1

The dataset global_economy explores the economic indicators as provided by the World Bank and covers the years from 1960 through 2017.

#exploring the data
?global_economy
head(global_economy)
#returns messy overlapping plot of GDP per capita over time for all countries
global_economy |>
  autoplot(GDP / Population) +
  labs(
    title = "GDP per Capita by Country Over Time",
    y = "GDP per Capita ($USD)",
    x = "Year"
  )

global_economy |>
  #filter for the country with the highest GDP per capita
  mutate(GDP_per_capita = GDP / Population) |>
  filter(GDP_per_capita == max(GDP_per_capita, na.rm = TRUE)) |>
  select(Country, Year, GDP_per_capita)
global_economy |>
  #shows which country has the maximum capita GDP per capita each year over time
  mutate(GDP_per_capita = GDP / Population) |>
  as_tibble() |>
  group_by(Year) |>
  filter(GDP_per_capita == max(GDP_per_capita, na.rm = TRUE)) |>
  select(Year, Country, GDP_per_capita) |>
  arrange(Year)

Monaco has the highest GDP per capita in the dataset at $185,829 in 2014. The maximum GDP per capita across all countries has shown significant growth over time, increasing from around $15,000 in 1960 to over $180,000 by 2014. This represents a large increase over the 54-year period. The growth shows a particular steepness since the 1990s.

# Top 10 countries by average GDP per capita
top_10 <- global_economy |>
  mutate(GDP_per_capita = GDP / Population) |>
  group_by(Country) |>
  summarise(avg_GDP = mean(GDP_per_capita, na.rm = TRUE)) |>
  arrange(desc(avg_GDP)) |>
  head(10)
## Warning: Current temporal ordering may yield unexpected results.
## ℹ Suggest to sort by `Country`, `Year` first.
## Current temporal ordering may yield unexpected results.
## ℹ Suggest to sort by `Country`, `Year` first.
top_10
## Warning: Current temporal ordering may yield unexpected results.
## ℹ Suggest to sort by `Country`, `Year` first.
# Plot these top countries over time
global_economy |>
  filter(Country %in% top_10$Country) |>
  mutate(GDP_per_capita = GDP / Population) |>
  autoplot(GDP_per_capita) +
  labs(
    title = "Top 10 Wealthiest Countries, GDP per Capita",
    x = "Year",
    y = "GDP per Capita ($USD)"
  )
## Warning: Removed 22 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Tells which countries held the record in different decades
global_economy |>
  mutate(GDP_per_capita = GDP / Population,
         Decade = floor(Year / 10) * 10) |>
  as_tibble() |>
  group_by(Decade) |>
  filter(GDP_per_capita == max(GDP_per_capita, na.rm = TRUE)) |>
  select(Decade, Year, Country, GDP_per_capita) |>
  distinct(Decade, .keep_all = TRUE)

#Question 3.2

US GDP from global economy

# Filter for United States
us_gdp <- global_economy |>
  filter(Country == "United States")

# Original data plot
us_gdp |>
  autoplot(GDP) +
  labs(
    title = "United States GDP - Original",
    y = "GDP (Billions USD)",
    x = "Year"
  )

# Apply log transformation
us_gdp |>
  autoplot(log(GDP)) +
  labs(
    title = "United States GDP - Log Transformed",
    y = "Log(GDP)",
    x = "Year"
  )

#Box-Cox transformation
us_gdp |>
  features(GDP, features = guerrero)
us_gdp |>
  autoplot(box_cox(GDP, lambda = 0)) +
  labs(
    title = "United States GDP - Box-Cox",
    y = "Transformed GDP",
    x = "Year"
  )

The original GDP shows strong exponential growth with increasing variance over time and the log transformation stabilizes the variance and converts the exponential growth to have a more linear trend, thus making the series easier to model and forecast. Since the variance is increasing with the level of the series and there is exponential growth, this data should be transformed.

# Filter for Victorian bulls, bullocks and steers
vic_bulls <- aus_livestock |>
  filter(
    State == "Victoria",
    Animal == "Bulls, bullocks and steers"
  )

# Plot original data
vic_bulls |>
  autoplot(Count) +
  labs(
    title = "Victorian Bulls, Bullocks and Steers Slaughter (Original)",
    y = "Count (thousands of head)",
    x = "Month"
  )

# Apply log transformation
vic_bulls |>
  autoplot(log(Count)) +
  labs(
    title = "Victorian Bulls, Bullocks and Steers Slaughter (Log Transformed)",
    y = "Log(Count)",
    x = "Month"
  )

# Box-Cox transformation
lambda_bulls <- vic_bulls |>
  features(Count, features = guerrero) |>
  pull(lambda_guerrero)

vic_bulls |>
  autoplot(box_cox(Count, lambda = lambda_bulls)) +
  labs(
    title = "Victorian Bulls, Bullocks and Steers Slaughter (Box-Cox Transformed)",
    y = "Transformed Count",
    x = "Month"
  )

The Victorian Bulls, Bullocks and Steers slaughter data shows relatively stable variance over time, with consistent seasonal fluctuations throughout the series.

Transformation isn’t necessary for this series because both the log and Box-Cox transformations produce plots that are visually nearly identical to the original data. This is because the variance remains relatively constant over time and the magnitude of seasonal fluctuations does not increase or decrease with the level of the series.

# Plot original data
vic_elec |>
  autoplot(Demand) +
  labs(
    title = "Victorian Electricity Demand (Original)",
    y = "Demand (MWh)",
    x = "Time"
  )

# Apply log transformation
vic_elec |>
  autoplot(log(Demand)) +
  labs(
    title = "Victorian Electricity Demand (Log Transformed)",
    y = "Log(Demand)",
    x = "Time"
  )

# Box-Cox transformation
lambda_elec <- vic_elec |>
  features(Demand, features = guerrero) |>
  pull(lambda_guerrero)

vic_elec |>
  autoplot(box_cox(Demand, lambda = lambda_elec)) +
  labs(
    title = "Victorian Electricity Demand (Box-Cox Transformed)",
    y = "Transformed Demand",
    x = "Time"
  )

The Victorian Electricity Demand data is a high-frequency series of half-hourly observations ranging from 2012-2015, with demand ranging from approximately 2,500 to 9,000 MWh. The series shows strong seasonal patterns with annual cycles, reflecting higher electricity demand during summer and winter months. Several extreme demand spikes are visible, particularly around 2014, likely corresponding to heat waves or extreme weather events.

Transformation isn’t necessary for this series as the log and Box-Cox transformations produce plots that are visually nearly identical to the original data. The variance remains relatively constant over time. While extreme spikes exist, they represent genuine demand events rather than variance heterogeneity. The magnitude of regular daily and seasonal fluctuations does not systematically increase or decrease with the level of the series.

# Plot original data
aus_production |>
  autoplot(Gas) +
  labs(
    title = "Australian Gas Production (Original)",
    y = "Gas Production (Petajoules)",
    x = "Quarter"
  )

# Apply log transformation
aus_production |>
  autoplot(log(Gas)) +
  labs(
    title = "Australian Gas Production (Log Transformed)",
    y = "Log(Gas)",
    x = "Quarter"
  )

# Box-Cox transformation
lambda_gas <- aus_production |>
  features(Gas, features = guerrero) |>
  pull(lambda_guerrero)

aus_production |>
  autoplot(box_cox(Gas, lambda = lambda_gas)) +
  labs(
    title = "Australian Gas Production (Box-Cox Transformed)",
    y = "Transformed Gas Production",
    x = "Quarter"
  )

For this series, working with transformed data would be strongly preferred for modeling purposes as the variance increases dramatically with the level of the series. Early periods show minimal seasonal fluctuations of less than 5 petajoules, while recent periods show large seasonal swings exceeding 50 petajoules.

Both the log and Box-Cox transformations successfully linearize the exponential growth trend and stabilize the variance, making seasonal fluctuations consistent across all time periods. The transformed data is much more amenable to standard forecasting methods.

#Question 3.3

#taking a look at the data
?canadian_gas
head(canadian_gas)
#plotting Canadian gas
canadian_gas |>
  autoplot(Volume) +
  labs(
    title = "Canadian Monthly Gas Production (Original)",
    y = "Gas Production (Billion cubic meters)",
    x = "Month"
  )

# Find optimal lambda using Guerrero method
lambda_canadian <- canadian_gas |>
  features(Volume, features = guerrero) |>
  pull(lambda_guerrero)

print(paste("Optimal lambda:", lambda_canadian))
## [1] "Optimal lambda: 0.576764791258982"
# Apply Box-Cox transformation
canadian_gas |>
  autoplot(box_cox(Volume, lambda = lambda_canadian)) +
  labs(
    title = "Canadian Monthly Gas Production (Box-Cox Transformed)",
    y = "Transformed Gas Production",
    x = "Month"
  )

Box-Cox transformation is unhelpful for the canadian_gas data because the variance remains relatively stable over time. The optimal lambda of 0.58 (close to 1, meaning minimal transformation) and the nearly identical appearance of both plots confirm this. Seasonal fluctuations stay consistent at 1-3 billion cubic meters throughout the series regardless of production level. Since transformation provides no variance stabilization benefit, the original scale is preferred for its interpretability.

#Question 3.4

# Load retail series
set.seed(12345678)  
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# Check series
myseries |>
  distinct(State, Industry)
# Plotting the original data
myseries |>
  autoplot(Turnover) +
  labs(
    title = "Retail Data (Original)",
    y = "Turnover ($Million AUD)",
    x = "Month"
  )

# Find optimal Box-Cox lambda using Guerrero method
lambda_retail <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

print(paste("Optimal lambda:", lambda_retail))
## [1] "Optimal lambda: 0.0830363119999255"
# Apply Box-Cox transformation with optimal lambda
myseries |>
  autoplot(box_cox(Turnover, lambda = lambda_retail)) +
  labs(
    title = "Retail Data (Box-Cox Transformed)",
    y = "Transformed Turnover",
    x = "Month"
  )

For my retail series, the Guerrero method selected an optimal Box-Cox transformation with lambda equal to 0.08. This lambda value close to 0 suggests a log transformation is appropriate. The transformation stabilizes the variance and makes seasonal fluctuations more consistent across time. Comparing the original and transformed plots, the Box-Cox transformation successfully addresses the increasing variance issue present in the original data and would be appropriate for modeling purposes.

#Question 3.5

# Explore tobacco data
aus_production |>
  autoplot(Tobacco) +
  labs(
    title = "Australian Tobacco Production (Original)",
    y = "Tobacco Production (tonnes)",
    x = "Quarter"
  )
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Find optimal lambda
lambda_tobacco <- aus_production |>
  features(Tobacco, features = guerrero) |>
  pull(lambda_guerrero)

print(paste("Optimal lambda for Tobacco:", lambda_tobacco))
## [1] "Optimal lambda for Tobacco: 0.926463585274383"
# Apply Box-Cox transformation
aus_production |>
  autoplot(box_cox(Tobacco, lambda = lambda_tobacco)) +
  labs(
    title = "Australian Tobacco Production (Box-Cox Transformed)",
    y = "Transformed Tobacco Production",
    x = "Quarter"
  )
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Filter for Economy class Melbourne-Sydney
economy_melb_syd <- ansett |>
  filter(
    Class == "Economy",
    Airports == "MEL-SYD"
  )

# Plot original
economy_melb_syd |>
  autoplot(Passengers) +
  labs(
    title = "Economy Passengers MEL-SYD (Original)",
    y = "Passengers",
    x = "Week"
  )

# Find optimal lambda
lambda_ansett <- economy_melb_syd |>
  features(Passengers, features = guerrero) |>
  pull(lambda_guerrero)

print(paste("Optimal lambda for Ansett:", lambda_ansett))
## [1] "Optimal lambda for Ansett: 1.9999267732242"
# Apply Box-Cox transformation
economy_melb_syd |>
  autoplot(box_cox(Passengers, lambda = lambda_ansett)) +
  labs(
    title = "Economy Passengers MEL-SYD (Box-Cox Transformed)",
    y = "Transformed Passengers",
    x = "Week"
  )

# Filter for Southern Cross Station
southern_cross <- pedestrian |>
  filter(Sensor == "Southern Cross Station")

# Plot original
southern_cross |>
  autoplot(Count) +
  labs(
    title = "Pedestrian Count at Southern Cross Station (Original)",
    y = "Pedestrian Count",
    x = "Date-Time"
  )

# Find optimal lambda
lambda_pedestrian <- southern_cross |>
  features(Count, features = guerrero) |>
  pull(lambda_guerrero)

print(paste("Optimal lambda for Pedestrian:", lambda_pedestrian))
## [1] "Optimal lambda for Pedestrian: -0.250161562391101"
# Apply Box-Cox transformation
southern_cross |>
  autoplot(box_cox(Count, lambda = lambda_pedestrian)) +
  labs(
    title = "Pedestrian Count at Southern Cross Station (Box-Cox Transformed)",
    y = "Transformed Count",
    x = "Date-Time"
  )

For the three series, the Guerrero method selected optimal lambda values of 0.93 for tobacco production, 2.0 for economy MEL-SYD passengers, and -0.25 for pedestrian counts at Southern Cross Station. Tobacco’s lambda near 1 indicates minimal transformation is needed as variance is already stable. Ansett’s squaring transformation provides limited benefit due to a structural break around 1989-1990 where counts drop to zero. Pedestrian’s inverse transformation helps compress the upper range but the series remains challenging due to regular nighttime zeros. None of these series exhibit the classic fan-shaped variance pattern where Box-Cox transformation would be most effective.

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

# Part a: Plot the time series
gas |>
  autoplot(Gas) +
  labs(
    title = "Australian Gas Production - Last 5 Years",
    y = "Gas Production (Petajoules)",
    x = "Quarter"
  )

The time series shows a clear upward trend over the five-year period, with gas production increasing from approximately 450 to 520 petajoules. Strong seasonal fluctuations are visible, with regular peaks and troughs occurring each year. The seasonal pattern appears consistent, with higher production in certain quarters (likely Q2 and Q3) and lower production in others.

# Part b: Classical decomposition with multiplicative model
gas_decomp <- gas |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components()

# View the components
gas_decomp
# Plot the decomposition
gas_decomp |>
  autoplot() +
  labs(title = "Classical Multiplicative Decomposition of Gas Production")
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

The classical decomposition produces four components: the original data, trend-cycle, seasonal indices, and random component. The trend-cycle shows smooth upward growth. The seasonal indices reveal a consistent quarterly pattern with values above and below 1, indicating multiplicative seasonal effects.

The decomposition results support the graphical interpretation from part (a). The smooth upward trend-cycle confirms the overall increasing pattern observed in the original plot. The seasonal component clearly shows the regular quarterly fluctuations, with peaks occurring in Q3 and troughs in Q1 matching the visual pattern of high summer production and low winter production. The consistency of the seasonal pattern across all five years validates the observation that seasonal amplitude remains stable.

gas |>
  autoplot(Gas) +
  autolayer(gas_decomp, season_adjust, color = "tomato") +
  labs(
    title = "Original vs Seasonally Adjusted Gas Production",
    y = "Gas Production (Petajoules)",
    x = "Quarter"
  )

The seasonally adjusted data (shown in red) removes the quarterly fluctuations, revealing the underlying trend-cycle more clearly. The adjusted series shows a smooth upward trend without the regular peaks and troughs present in the original data, making it easier to identify the long-term growth pattern in gas production.

# Part e: Create outlier in the middle
gas_outlier <- gas |>
  mutate(Gas = if_else(row_number() == 10, Gas + 300, Gas))

# Decompose with outlier
gas_outlier_decomp <- gas_outlier |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components()

# Plot comparison
gas_outlier |>
  autoplot(Gas) +
  autolayer(gas_outlier_decomp, season_adjust, color = "red") +
  labs(
    title = "Seasonally Adjusted with Outlier (Middle)",
    y = "Gas Production",
    x = "Quarter"
  )

Adding an outlier of 300 petajoules in the middle of the series significantly distorts the decomposition. The outlier creates a spike in both the trend-cycle component and the seasonally adjusted data. The moving average smoothing spreads the outlier’s effect to surrounding observations, creating an artificial hump in the trend around that period.

# Part f: Outlier at the end
gas_outlier_end <- gas |>
  mutate(Gas = if_else(row_number() == n(), Gas + 300, Gas))

# Decompose
gas_outlier_end_decomp <- gas_outlier_end |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components()

# Plot
gas_outlier_end |>
  autoplot(Gas) +
  autolayer(gas_outlier_end_decomp, season_adjust, color = "red") +
  labs(
    title = "Seasonally Adjusted with Outlier (End)",
    y = "Gas Production",
    x = "Quarter"
  )

Yes, the location matters. An outlier at the end has less impact on the trend-cycle estimation because there are fewer observations for the moving average to smooth over. The end outlier creates a sharp spike in the seasonally adjusted data but affects fewer surrounding observations compared to a middle outlier, which distorts the trend over a wider range due to the symmetric moving average.

#Question 3.6 To show that a 3×5 MA is equivalent to a 7-term weighted moving average, we apply the moving averages sequentially.

5-MA

The 5-term moving average at time t is: M_t^(5) = (1/5)(y_{t-2} + y_{t-1} + y_t + y_{t+1} + y_{t+2})

3-MA

The 3-term moving average of M_t^(5) is: M_t^(3×5) = (1/3)(M_{t-1}^(5) + M_t^(5) + M_{t+1}^(5))

Substitute and expand

Substituting the 5-MA expressions: M_t^(3×5) = (1/3)[(1/5)(y_{t-3} + y_{t-2} + y_{t-1} + y_t + y_{t+1}) + (1/5)(y_{t-2} + y_{t-1} + y_t + y_{t+1} + y_{t+2}) + (1/5)(y_{t-1} + y_t + y_{t+1} + y_{t+2} + y_{t+3})]

** terms**

M_t^(3×5) = (1/15)[y_{t-3} + 2y_{t-2} + 3y_{t-1} + 3y_t + 3y_{t+1} + 2y_{t+2} + y_{t+3}]

weighted average

M_t^(3×5) = (1/15)y_{t-3} + (2/15)y_{t-2} + (3/15)y_{t-1} + (3/15)y_t + (3/15)y_{t+1} + (2/15)y_{t+2} + (1/15)y_{t+3}

Converting to decimals: - 1/15 -> 0.067 - 2/15 -> 0.133 - 3/15 -> 0.200

Therefore, the 3×5 MA equals: M_t = 0.067y_{t-3} + 0.133y_{t-2} + 0.200y_{t-1} + 0.200y_t + 0.200y_{t+1} + 0.133y_{t+2} + 0.067y_{t+3}

This confirms the given weights: 0.067, 0.133, 0.200, 0.200, 0.200, 0.133, 0.067.

Question 3.7 - Your Complete Answer:

Part (a): The time series shows a clear upward trend over the five-year period, with gas production increasing from approximately 450 to 520 petajoules. Strong seasonal fluctuations are visible, with regular peaks and troughs occurring each year. The seasonal pattern appears consistent, with higher production in certain quarters (likely Q2 and Q3) and lower production in others.

Part (b): The classical decomposition produces four components: the original data, trend-cycle, seasonal indices, and random component. The trend-cycle shows smooth upward growth. The seasonal indices reveal a consistent quarterly pattern with values above and below 1, indicating multiplicative seasonal effects.

Part (c): The decomposition results support the graphical interpretation from part (a). The smooth upward trend-cycle confirms the overall increasing pattern observed in the original plot. The seasonal component clearly shows the regular quarterly fluctuations, with peaks occurring in Q3 and troughs in Q1 matching the visual pattern of high summer production and low winter production. The consistency of the seasonal pattern across all five years validates the observation that seasonal amplitude remains stable.

Part (d): The seasonally adjusted data (shown in red) removes the quarterly fluctuations, revealing the underlying trend-cycle more clearly. The adjusted series shows a smooth upward trend without the regular peaks and troughs present in the original data, making it easier to identify the long-term growth pattern in gas production.

Part (e): Adding an outlier of 300 petajoules in the middle of the series significantly distorts the decomposition. The outlier creates a spike in both the trend-cycle component and the seasonally adjusted data. The moving average smoothing spreads the outlier’s effect to surrounding observations, creating an artificial hump in the trend around that period.

Part (f): Yes, the location matters. An outlier at the end has less impact on the trend-cycle estimation because there are fewer observations for the moving average to smooth over. The end outlier creates a sharp spike in the seasonally adjusted data but affects fewer surrounding observations compared to a middle outlier, which distorts the trend over a wider range due to the symmetric moving average.

#Question3.8

# Load the same retail series from Question 3.4
set.seed(12345678)  
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# Check which series
myseries |>
  distinct(State, Industry)
myseries |>
  autoplot(Turnover) +
  labs(
    title = "Retail Turnover - Time Series",
    y = "Turnover ($Million AUD)",
    x = "Month"
  )

# Seasonal plot
myseries |>
  gg_season(Turnover) +
  labs(title = "Seasonal Plot of Retail Turnover")
## Warning: `gg_season()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_season()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Subseries plot
myseries |>
  gg_subseries(Turnover) +
  labs(title = "Subseries Plot of Retail Turnover")
## Warning: `gg_subseries()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_subseries()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Lag plot
myseries |>
  gg_lag(Turnover, geom = "point") +
  labs(title = "Lag Plots of Retail Turnover")
## Warning: `gg_lag()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_lag()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# ACF plot
myseries |>
  ACF(Turnover) |>
  autoplot() +
  labs(title = "ACF of Retail Turnover")

The retail series (Northern Territory) shows a strong upward trend with clear monthly seasonality. December consistently peaks due to holiday shopping, followed by January troughs. The seasonal plot confirms this pattern, with increasing seasonal variation over time. The subseries plot shows upward trends across all months while maintaining seasonal differences.

The lag plots display strong positive correlations at all lags, especially at lag 12, indicating pronounced seasonality. The ACF shows slow decay characteristic of trend, with a scalloped pattern peaking at lags 12, 24, and 36, confirming strong annual seasonality. Both trend and seasonality are clearly present, with no distinct cyclicity beyond the regular seasonal pattern.

#Question 3.9

#Question 3.9

# STL decomposition of retail series
stl_decomp <- myseries |>
  model(STL(Turnover ~ trend(window = 21) + season(window = 13))) |>
  components()

# Plot decomposition
stl_decomp |>
  autoplot() +
  labs(title = "STL Decomposition of Retail Turnover")

The STL decomposition reveals the retail turnover separated into trend, seasonal, and remainder components. The trend shows smooth long-term growth from approximately $3 million to $15 million, with some acceleration in recent years and increased volatility after 2015. The seasonal component displays consistent monthly patterns throughout, with the amplitude of seasonal fluctuations growing proportionally with the trend level.

The remainder component reveals several notable features not immediately apparent in the original series. A large spike appears around 1990, likely corresponding to an economic event or data anomaly. Another significant spike occurs around 2000, and increased irregular variation is visible in the post-2015 period, suggesting structural changes in the retail environment or increased volatility. The remainder is generally small relative to the trend and seasonal components, indicating the series is well-captured by these two elements, but the identified spikes represent genuine outliers or unusual periods worth investigating further.