Rpubs URL: https://rpubs.com/umaisabdullah/DATA624_HW2

Objective

Complete and submit exercises 3.1, 3.2, 3.3, 3.4, 3.5, 3.7, 3.8 and 3.9 for Homework 2 from the Hyndman online Forecasting book. We have to submit both our Rpubs link as well as attach the .pdf file with our code.

Overview

Homework 1 was about looking at time series and naming the patterns we could see. This homework is about two things we do to a series before modelling it.

The first is transformation. Many series get noisier as they get bigger, so the size of the wiggle depends on the level of the series. That is a problem for models that assume the noise stays about the same size throughout. A Box-Cox transformation squashes the large values more than the small ones and evens that variation out.

The second is decomposition. This splits a series into three parts: a trend cycle for the slow movement, a seasonal part for the pattern that repeats each year, and a remainder for whatever is left over. Once the parts are separated we can look at each one on its own, which makes it much easier to see what is really driving the series.

Pre-Requisite

All of the data for this homework comes from packages, so nothing needs to be downloaded into the working directory. Two packages are required:

  1. fpp3 loads the core forecasting functions and all of the datasets used below.
  2. seasonal provides the X-11 method used in Exercise 3.8.

Exercise 3.1

Question. 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?

My thought process. The dataset gives GDP and population as separate columns, so GDP per capita is not in the data and I have to build it by dividing one by the other. That division is the whole point of the exercise. Comparing raw GDP would only tell me which countries are large, since a big country produces more in total no matter how well off its people are. Dividing by population puts every country on a per person basis, which is what makes them comparable.

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

gdppc |>
  autoplot(GDP_per_capita, show.legend = FALSE) +
  labs(title = "GDP per capita for all countries",
       x = "Year", y = "GDP per capita (USD)")

That plot draws one line per country with no legend, so it is far too crowded to read any single country from. It is still worth showing, because it makes clear that most countries sit in a low band while a few rise well above the rest. To actually answer the question I need to ask the data directly instead of reading the plot.

# Highest GDP per capita ever recorded, and the year it happened.
gdppc |>
  as_tibble() |>
  filter(!is.na(GDP_per_capita)) |>
  slice_max(GDP_per_capita, n = 5) |>
  select(Country, Year, GDP_per_capita)
## # A tibble: 5 × 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.
# Which country leads in each of the most recent years.
gdppc |>
  as_tibble() |>
  filter(!is.na(GDP_per_capita)) |>
  group_by(Year) |>
  slice_max(GDP_per_capita, n = 1) |>
  ungroup() |>
  filter(Year >= 2008) |>
  select(Year, Country, GDP_per_capita)
## # A tibble: 10 × 3
##     Year Country       GDP_per_capita
##    <dbl> <fct>                  <dbl>
##  1  2008 Monaco               180640.
##  2  2009 Monaco               149221.
##  3  2010 Monaco               144569.
##  4  2011 Monaco               162155.
##  5  2012 Monaco               152000.
##  6  2013 Liechtenstein        173528.
##  7  2014 Monaco               185153.
##  8  2015 Liechtenstein        167591.
##  9  2016 Monaco               168011.
## 10  2017 Luxembourg           104103.
# Plot only the handful of countries that have led at some point.
leaders <- c("Monaco", "Liechtenstein", "Luxembourg", "Norway")

gdppc |>
  filter(Country %in% leaders) |>
  autoplot(GDP_per_capita) +
  labs(title = "GDP per capita for the usual leaders",
       x = "Year", y = "GDP per capita (USD)")

What I found. Monaco has the highest GDP per capita in the data, reaching about 185,000 US dollars in 2014. Liechtenstein is close behind, and Luxembourg and Norway are the other two that appear near the top.

How it has changed over time. Two things stand out. The first is that all of these countries trend upward, so GDP per capita grows over the period rather than staying flat. The second is that the leader is not always the same country. Monaco and Liechtenstein trade the top spot back and forth depending on the year, and the data for these very small countries is missing in some years, which is part of why the leader jumps around.

The countries at the top are worth a comment. Monaco, Liechtenstein and Luxembourg are all very small. That is not a coincidence. A small country with a large banking or finance sector produces a lot of GDP relative to a very small population, so the ratio comes out high. This is a good reminder that GDP per capita measures output per person, not how much money an ordinary person actually takes home.

Exercise 3.2

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

My thought process. The question says to transform only if it seems appropriate, so the real task is deciding when a transformation helps and when it does not. My rule is to look at whether the size of the variation changes with the level of the series. If the swings get bigger as the series climbs, a transformation will even them out and is worth doing. If the swings stay about the same size throughout, transforming adds a step and makes the units harder to explain for no real gain. So I plot each series first and decide afterward, rather than transforming everything out of habit.

United States GDP

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

us_gdp |>
  autoplot(GDP) +
  labs(title = "United States GDP", x = "Year", y = "GDP (USD)")

This curves upward rather than rising in a straight line, which is what steady percentage growth looks like. A transformation can straighten that out.

lambda_us <- us_gdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

lambda_us
## [1] 0.2819443
us_gdp |>
  autoplot(box_cox(GDP, lambda_us)) +
  labs(title = "United States GDP after Box-Cox transformation",
       x = "Year",
       y = paste0("Transformed GDP (lambda = ", round(lambda_us, 3), ")"))

Effect. The Guerrero method picks a lambda near 0.28. After transforming, the curve becomes close to a straight line. That is the useful part. A straight line is easier to model and to extend forward than a curve, because a constant slope is simpler to describe than a slope that keeps increasing.

Bulls, bullocks and steers in Victoria

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

bulls |>
  autoplot(Count) +
  labs(title = "Slaughter of bulls, bullocks and steers in Victoria",
       x = "Month", y = "Count")

bulls |>
  features(Count, features = guerrero) |>
  pull(lambda_guerrero)
## [1] -0.04461887

Effect. Here I would leave the series alone. The level drifts downward over the years, but the size of the month to month variation stays roughly similar instead of shrinking along with the level. Guerrero returns a lambda of about -0.04, which is close to zero and so amounts to a log transformation, but applying it does not visibly improve the plot. Since a transformation makes the numbers harder to explain and there is no clear payoff, the honest choice is not to transform.

This is the case the question is really testing. Not every series needs a transformation, and knowing when to leave one alone matters as much as knowing how to transform it.

Victorian electricity demand

vic_elec |>
  autoplot(Demand) +
  labs(title = "Half-hourly electricity demand in Victoria",
       x = "Time", y = "Demand (MWh)")

The half-hourly data is so dense that it prints as a solid block, so I aggregate to daily totals to actually see the shape.

vic_daily <- vic_elec |>
  index_by(Date) |>
  summarise(Demand = sum(Demand))

vic_daily |>
  autoplot(Demand) +
  labs(title = "Daily electricity demand in Victoria",
       x = "Date", y = "Demand (MWh)")

Effect. I would not transform this one either, but for a different reason than the cattle series. The variation here is not growing with the level. It is driven by temperature, with demand spiking in summer for air conditioning and rising again in winter for heating. That produces sharp peaks at both extremes, and a Box-Cox transformation does not fix that, because the problem is not unstable variance. It is a relationship with weather. The right fix would be to bring temperature into the model later on, not to reshape the series now.

Gas production

aus_production |>
  autoplot(Gas) +
  labs(title = "Australian gas production",
       x = "Quarter", y = "Gas (petajoules)")

This is the clearest case for transforming in the whole exercise. The seasonal swings are tiny at the start of the series and very large by the end, so the size of the pattern grows right along with the level.

lambda_gas <- aus_production |>
  features(Gas, features = guerrero) |>
  pull(lambda_guerrero)

lambda_gas
## [1] 0.1095171
aus_production |>
  autoplot(box_cox(Gas, lambda_gas)) +
  labs(title = "Australian gas production after Box-Cox transformation",
       x = "Quarter",
       y = paste0("Transformed gas (lambda = ", round(lambda_gas, 3), ")"))

Effect. Guerrero picks a lambda near 0.11, which is close to a log transformation. After transforming, the seasonal swings are about the same size across the whole series instead of fanning out at the end. That is exactly what a variance stabilising transformation is for, and this series is the textbook example of when to use one.

Exercise 3.3

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

My thought process. The way the question is phrased tells me the expected answer is no, so my job is to explain why rather than to hunt for the best lambda. A Box-Cox transformation fixes one specific problem, which is variance that grows or shrinks along with the level of the series. So the way to answer is to check whether this series actually has that problem. If the variance moves in some other way, then Box-Cox is the wrong tool and no choice of lambda will help.

canadian_gas |>
  autoplot(Volume) +
  labs(title = "Monthly Canadian gas production",
       x = "Month", y = "Volume (billions of cubic metres)")

canadian_gas |>
  gg_season(Volume) +
  labs(title = "Seasonal plot: Canadian gas production",
       y = "Volume (billions of cubic metres)")

lambda_cg <- canadian_gas |>
  features(Volume, features = guerrero) |>
  pull(lambda_guerrero)

lambda_cg
## [1] 0.5767648
canadian_gas |>
  autoplot(box_cox(Volume, lambda_cg)) +
  labs(title = "Canadian gas production after Box-Cox transformation",
       x = "Month",
       y = paste0("Transformed volume (lambda = ", round(lambda_cg, 3), ")"))

Why Box-Cox does not help here. Look at how the seasonal variation behaves over time. It is small in the 1960s, grows to its largest through the middle of the series around the 1970s and 1980s, and then gets smaller again in the 1990s. The variation does not simply increase with the level of the series. It increases and then decreases while the level keeps climbing the whole time.

That is the heart of the problem. A Box-Cox transformation applies one fixed rule to the entire series, so it can only correct variance that moves in one direction along with the level. Here a lambda small enough to calm the large middle years would over-squash the later ones, and a lambda that suits the ends would leave the middle untouched. There is no single value that works for all three periods at once, which is why the transformed plot still shows uneven variation.

The deeper issue is that the shape of the seasonal pattern itself changes over time, which is a different problem than unstable variance. Fixing that needs a method that lets the seasonal component evolve, such as an STL decomposition with a flexible seasonal window, rather than one transformation applied to everything at once.

Exercise 3.4

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

My thought process. This one continues from Homework 1, so I use the same seed and get the same series. That matters because the answer depends on which series was drawn, and using a different one would make this exercise inconsistent with the earlier homework. Then I let the Guerrero method choose lambda rather than picking a value by eye, and I plot the result to confirm the choice actually did something useful.

set.seed(624)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# Confirm this is the same series as Homework 1.
myseries |>
  as_tibble() |>
  distinct(State, Industry, `Series ID`)
## # A tibble: 1 × 3
##   State           Industry               `Series ID`
##   <chr>           <chr>                  <chr>      
## 1 New South Wales Takeaway food services A3349792X
myseries |>
  autoplot(Turnover) +
  labs(title = "Retail turnover before transformation",
       x = "Month", y = "Turnover ($ million)")

The seasonal swings clearly get bigger as the series rises, so a transformation is worth doing here.

lambda_retail <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

lambda_retail
## [1] 0.002144737
myseries |>
  autoplot(box_cox(Turnover, lambda_retail)) +
  labs(title = "Retail turnover after Box-Cox transformation",
       x = "Month",
       y = paste0("Transformed turnover (lambda = ", round(lambda_retail, 4), ")"))

What I would select. Guerrero picks a lambda of about 0.002, which is so close to zero that it is effectively a log transformation. In practice I would just use lambda equal to 0 and take logs.

Why I would round it to zero. A lambda of exactly 0 is defined as the log, and logs have a meaning people can actually explain. A change on the log scale is roughly a percentage change in the original units, so saying turnover grew by a steady percentage each year is something a reader can follow. A lambda of 0.002 gives nearly identical numbers but has no plain English meaning at all. When two choices fit equally well, the one that is easier to explain is the better choice.

Effect. After the transformation the seasonal swings are about the same size from the start of the series to the end, rather than fanning out as turnover grows. The December spike is still clearly there, which is what I want. The goal is to even out the size of the variation, not to erase the seasonal pattern.

Exercise 3.5

Question. 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.

My thought process. Each of these three needs some filtering before I can touch it, and the pedestrian one has a trap in it that I checked for before writing any code. I let Guerrero pick lambda each time, then look at whether the answer actually makes sense rather than just reporting the number it gives me.

Tobacco

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), ")"))

What I selected. Guerrero picks a lambda near 0.93, which is very close to 1.

A lambda of 1 means no transformation at all apart from shifting the numbers. So the method is essentially telling me this series does not need transforming. The before and after plots look almost the same, which confirms it. The variance here is already fairly steady, so there is nothing for Box-Cox to fix.

I am reporting this as it is rather than pretending the transformation did something useful. A lambda near 1 is a real answer, and it means leave the series as it is.

Economy class passengers between Melbourne and Sydney

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

ansett_econ |>
  autoplot(Passengers) +
  labs(title = "Economy class passengers, Melbourne to Sydney",
       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 = "Economy passengers after Box-Cox transformation",
       x = "Week",
       y = paste0("Transformed passengers (lambda = ", round(lambda_ansett, 3), ")"))

What I selected. Guerrero picks a lambda near 2, which squares the values rather than squashing them.

That is an unusual answer and it is worth stopping on. A lambda above 1 spreads the large values further apart instead of pulling them in, which is the opposite of what a variance stabilising transformation normally does. What is really going on is that this series has a stretch in 1989 where passenger numbers fall to zero, because an industrial dispute grounded the airline. That flat run of zeros dominates the calculation, so the lambda is responding to the dispute rather than to any natural pattern in the variance.

So my honest answer is that Box-Cox is not doing useful work on this series. The unusual period is the real feature of the data, and I would handle that stretch directly rather than trusting a lambda that was driven by it.

Pedestrian counts at Southern Cross Station

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

# Two problems to check for before transforming.
sum(ped$Count == 0)                             # hours with a count of zero
## [1] 159
nrow(has_gaps(ped) |> filter(.gaps))            # is the series missing hours?
## [1] 1

There are two problems here and either one would cause trouble. The hourly data contains 159 readings of exactly zero, which happens overnight when nobody walks past the sensor. It also has gaps where hours are missing altogether. A log transformation is undefined at zero, so a lambda near 0 would produce infinite values.

Aggregating to daily totals fixes both problems at once. A whole day always has some foot traffic, so the zeros disappear, and summing by date gives a complete daily series.

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

# Confirm both problems are gone.
sum(ped_daily$Count == 0)
## [1] 0
nrow(has_gaps(ped_daily) |> filter(.gaps))
## [1] 0
ped_daily |>
  autoplot(Count) +
  labs(title = "Daily pedestrian counts at Southern Cross Station",
       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 = "Daily pedestrian counts after Box-Cox transformation",
       x = "Date",
       y = paste0("Transformed count (lambda = ", round(lambda_ped, 3), ")"))

What I selected. On the daily series Guerrero picks a lambda near 0.27, which sits between a log and a square root. After the transformation the low days are pulled up closer to the busy ones, so the variation is more even across the series.

The lesson from this one is that checking the data first mattered more than the lambda itself. If I had run Box-Cox straight on the hourly counts without looking, the zeros would have caused a problem that is easy to miss.

Exercise 3.7

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

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

My thought process. This exercise builds on itself, so each part uses the one before it. The interesting part is the last two, where I deliberately damage the data and watch what happens. That is a good way to learn how a method behaves, because I know exactly what I changed and can put any difference down to that one change.

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 time series

gas |>
  autoplot(Gas) +
  labs(title = "Australian gas production, last five years",
       x = "Quarter", y = "Gas (petajoules)")

What I see. Both patterns are present. There is a clear seasonal pattern that repeats every year, with production peaking in Q3 and dropping to its lowest in Q1. Q3 is the Australian winter, so more gas is being used for heating. On top of that there is a gentle upward trend across the five years, though it is mild compared with the size of the seasonal swing.

(b) Classical multiplicative decomposition

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

gas_decomp |>
  autoplot() +
  labs(title = "Classical multiplicative decomposition of gas production")

# The seasonal index for each quarter.
gas_decomp |>
  as_tibble() |>
  mutate(Quarter_num = quarter(Quarter)) |>
  group_by(Quarter_num) |>
  summarise(seasonal_index = round(mean(seasonal, na.rm = TRUE), 4))
## # A tibble: 4 × 2
##   Quarter_num seasonal_index
##         <int>          <dbl>
## 1           1          0.875
## 2           2          1.07 
## 3           3          1.13 
## 4           4          0.925

(c) Do the results support part (a)?

Yes, and the numbers make it precise. The seasonal indices come out at about 0.875 for Q1, 1.074 for Q2, 1.126 for Q3, and 0.925 for Q4.

Because this is a multiplicative decomposition, these are multipliers rather than amounts. Q3 at 1.126 means production in Q3 runs about 12.6 percent above the trend, and Q1 at 0.875 means production runs about 12.5 percent below it. That matches exactly what I said from the plot, that Q3 is the peak and Q1 is the low point, but now I have numbers instead of an impression.

The trend panel also confirms the mild upward movement I described in part a.

(d) Seasonally adjusted data

gas_decomp |>
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "Original")) +
  geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
  labs(title = "Gas production, original and seasonally adjusted",
       x = "Quarter", y = "Gas (petajoules)", colour = "Series")

What seasonally adjusting does. It divides out the seasonal index, which removes the regular yearly pattern and leaves the trend plus whatever is left over. The adjusted line is much smoother and the upward drift is easier to see, because the seasonal swing is no longer hiding it.

(e) Effect of an outlier in the middle

gas_outlier_mid <- gas |>
  mutate(Gas = if_else(row_number() == 10, Gas + 300, Gas))

decomp_mid <- gas_outlier_mid |>
  model(classical_decomposition(Gas, type = "multiplicative")) |>
  components()

decomp_mid |>
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "With outlier")) +
  geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
  labs(title = "Effect of an outlier in the middle of the series",
       x = "Quarter", y = "Gas (petajoules)", colour = "Series")

# How the seasonal indices changed compared with the clean data.
decomp_mid |>
  as_tibble() |>
  mutate(Quarter_num = quarter(Quarter)) |>
  group_by(Quarter_num) |>
  summarise(seasonal_index = round(mean(seasonal, na.rm = TRUE), 4))
## # A tibble: 4 × 2
##   Quarter_num seasonal_index
##         <int>          <dbl>
## 1           1          0.821
## 2           2          0.998
## 3           3          1.06 
## 4           4          1.12

What the outlier does. Adding 300 to one observation causes damage in two separate places, and the second one is the more interesting of the two.

The obvious effect is a large spike in the seasonally adjusted series at that point. That is expected, since the outlier is a genuine jump in the data and seasonal adjustment does not remove it.

The less obvious effect is that the seasonal indices themselves move. Comparing the table above with the clean one from part b, the Q1 index drops from about 0.875 to 0.821 and the Q2 index falls from about 1.074 to 0.998. Those quarters are not where I put the outlier, yet their indices changed anyway.

That happens because classical decomposition builds the index for each quarter by averaging across all the years, so one huge value pulls the averages around. The result is that a single bad observation contaminates the estimate for other quarters and other years, not just the point where it appeared. That is a real weakness of classical decomposition, and it is one of the reasons more robust methods like STL are usually preferred.

(f) Outlier near the end instead of the middle

gas_outlier_end <- gas |>
  mutate(Gas = if_else(row_number() == n() - 1, Gas + 300, Gas))

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

decomp_end |>
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Gas, colour = "With outlier")) +
  geom_line(aes(y = season_adjust, colour = "Seasonally adjusted")) +
  labs(title = "Effect of an outlier near the end of the series",
       x = "Quarter", y = "Gas (petajoules)", colour = "Series")

# Compare the trend estimates from all three versions.
mid_trend <- decomp_mid |> as_tibble() |> select(Quarter, trend_mid = trend)
end_trend <- decomp_end |> as_tibble() |> select(Quarter, trend_end = trend)

gas_decomp |>
  as_tibble() |>
  select(Quarter, trend_clean = trend) |>
  left_join(mid_trend, by = "Quarter") |>
  left_join(end_trend, by = "Quarter") |>
  mutate(across(where(is.numeric), \(x) round(x, 1)))
## # A tibble: 20 × 4
##    Quarter trend_clean trend_mid trend_end
##      <qtr>       <dbl>     <dbl>     <dbl>
##  1 2005 Q3         NA        NA        NA 
##  2 2005 Q4         NA        NA        NA 
##  3 2006 Q1        200.      200.      200.
##  4 2006 Q2        204.      204.      204.
##  5 2006 Q3        207       207       207 
##  6 2006 Q4        210.      210.      210.
##  7 2007 Q1        213       213       213 
##  8 2007 Q2        216.      254.      216.
##  9 2007 Q3        219.      294.      219.
## 10 2007 Q4        219.      294.      219.
## 11 2008 Q1        219.      294.      219.
## 12 2008 Q2        219       256.      219 
## 13 2008 Q3        219       219       219 
## 14 2008 Q4        220.      220.      220.
## 15 2009 Q1        222.      222.      222.
## 16 2009 Q2        223.      223.      223.
## 17 2009 Q3        225.      225.      263.
## 18 2009 Q4        226       226       301 
## 19 2010 Q1         NA        NA        NA 
## 20 2010 Q2         NA        NA        NA

Does the position matter? Yes. The difference comes down to how the trend is calculated.

Classical decomposition estimates the trend with a moving average, and a moving average needs values on both sides of a point. That means the trend cannot be computed for the first two and last two quarters of the series, which is why those entries are blank in the table above.

The table makes the difference easy to count. With the outlier in the middle, the trend is wrong for five quarters in a row, from 2007 Q2 through 2008 Q2, where it jumps to around 294 instead of the true 219. With the outlier near the end, only two quarters are affected, 2009 Q3 and Q4.

The reason is that an outlier in the middle sits fully inside the region where the trend is estimated, so every moving average window that contains it gets pulled upward, and that is several windows. An outlier near the end falls partly outside that region, because the last two quarters have no trend estimate at all. Fewer windows include it, so it does less damage.

The practical lesson is that the same bad value does different amounts of damage depending on where it lands, and a point in the middle of a series is the one to worry about most.

Exercise 3.8

Question. 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?

My thought process. X-11 is a more sophisticated decomposition than the classical one I used in Exercise 3.7. It estimates the trend all the way to both ends of the series instead of leaving gaps, and it is more resistant to outliers. The question asks what it reveals that I had not already noticed, so rather than eyeballing the remainder panel I will flag values that sit unusually far from the rest and let the numbers point me at 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")

# Flag irregular values more than three standard deviations from the mean.
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

What X-11 revealed. The test flags ten months where the irregular component sits more than three standard deviations away from its mean. A few of them are worth a comment.

The March values in 1993 and 1994 stand out, one unusually low and the next unusually high in consecutive years. That pattern is a hint about Easter. Easter moves between March and April from one year to the next, so a takeaway food series can get a boost in March one year and in April the next. A fixed monthly seasonal pattern cannot follow a holiday that moves around, so the leftover effect ends up in the irregular component.

January 2002 is the largest positive value in the list. October and December 1983 are both unusually low and sit near the very start of the series, where the data may simply be less reliable than it is later on.

What this adds over Homework 1. In Homework 1 I described this series as having a clear December peak, a strong upward trend, and no obvious cycle. That was all correct, but at that level of detail individual odd months were invisible, because the trend and the seasonal pattern dominate the plot and everything else is small next to them.

Decomposition helps precisely because it takes those two dominant parts out. Once the trend and season are removed, what is left is small enough that a single unusual month becomes easy to see. The moving Easter effect is a good example, since it is completely hidden in the original plot but shows up clearly once the series is decomposed.

Exercise 3.9

Question. 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.

My thought process. This exercise has no code, because the figures are printed in the textbook rather than produced from a dataset I have. The question puts weight on the scales of the graphs, so that is what I focus on. Reading the scales is the whole trick here, since a panel can look dramatic while actually showing a very small effect once you check the numbers on its axis.

(a) Describing the decomposition

The decomposition is dominated by the trend. The labour force grows steadily from roughly 6.4 million to 9 million people across the seventeen years, and the trend panel covers by far the widest range of values of any component. The seasonal component is small next to that it spans only about ±100,000 people, close to one percent of the total labour force. This is where reading the scales matters: the seasonal panel looks like a large repeating wave, but that’s only because it’s drawn on a far more magnified scale than the trend panel. The grey reference bars on the right of each panel make this explicit they’re roughly the same physical height, yet represent ranges that differ by two orders of magnitude. The seasonal pattern also changes shape gradually over the period rather than staying fixed (for example, March’s seasonal peak grows from around 65 in the early 1980s to nearly 90 by the mid-80s, then falls back to around 40 by 1995), which is why a method like STL that lets the season evolve suits this data better than a classical decomposition would. The remainder is mostly small and without pattern, apart from a short, sharp episode of unusually large negative values in the early 1990s that stands well outside its normal range.

(b) Is the 1991/1992 recession visible?

Yes, and it shows up in two of the three components. In the trend panel, the steady climb pauses and flattens out around 1990-1992 rather than continuing to rise, before resuming its climb afterward. The clearer sign is in the remainder panel, which shows a sharp double-dip down to around -380 far larger than the normal noise band of roughly ±100 seen elsewhere in that panel. That means the actual labour force fell well below what the trend and seasonal components together would predict, and that gap is what a recession looks like in a decomposition. The seasonal component is essentially unchanged through the recession, which makes sense: a recession is an economic event rather than a calendar one, so it interrupts the trend and produces a spike in the remainder, but it doesn’t change which months of the year are typically busy.