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?
# A tsibble: 15,150 x 9 [1Y]
# Key: Country [263]
Country Code Year GDP Growth CPI Imports Exports Population
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Afghanistan AFG 1960 537777811. NA NA 7.02 4.13 8996351
2 Afghanistan AFG 1961 548888896. NA NA 8.10 4.45 9166764
3 Afghanistan AFG 1962 546666678. NA NA 9.35 4.88 9345868
4 Afghanistan AFG 1963 751111191. NA NA 16.9 9.17 9533954
5 Afghanistan AFG 1964 800000044. NA NA 18.1 8.89 9731361
6 Afghanistan AFG 1965 1006666638. NA NA 21.4 11.3 9938414
7 Afghanistan AFG 1966 1399999967. NA NA 18.6 8.57 10152331
8 Afghanistan AFG 1967 1673333418. NA NA 14.2 6.77 10372630
9 Afghanistan AFG 1968 1373333367. NA NA 15.2 8.90 10604346
10 Afghanistan AFG 1969 1408888922. NA NA 15.0 10.1 10854428
# ℹ 15,140 more rows
Given that plotting a total of 263 countries at once would be unreadable, I’ll proceed to create the GDP per capita column first and then plot everything with no legends
gdp_per_capita <- global_economy |>mutate(GDP_per_capita = GDP / Population)gdp_per_capita |>autoplot(GDP_per_capita) +labs(title ="GDP per capita", y ="US dollars") +theme(legend.position ="none")
Warning: `autoplot.tbl_ts()` was deprecated in fabletools 0.6.0.
ℹ Please use `ggtime::autoplot.tbl_ts()` instead.
ℹ Graphics functions have been moved to the {ggtime} package. Please use
`library(ggtime)` instead.
Warning: Removed 3242 rows containing missing values or values outside the scale range
(`geom_line()`).
The plot shows one gdp pulling far away from the rest and this is a strong visual hint about which country it is, therefore this is a clear indicator of the single row with the GDP per capita
The US GDP shows a clear upward trend that increases over time in dollar terms which is much higher in the last decade that in early years. This was a great example to use since the guerrero method selected a lambda of ~0.28. When I applied this transformation it shows a series that increases in a much closer linear trend macros the years.
Slaughter of Victorian “Bulls, bullocks and steers” in aus_livestock.
aus_livestock |>filter(State =="Victoria", Animal =="Bulls, bullocks and steers") |>autoplot(Count) +labs(title ="Bulls, bullocks and steers")
aus_livestock |>filter(State =="Victoria", Animal =="Bulls, bullocks and steers") |>features(Count, features = guerrero) |>pull(lambda_guerrero)
[1] -0.04461887
lambda_livestock <- aus_livestock |>filter(State =="Victoria", Animal =="Bulls, bullocks and steers") |>features(Count, features = guerrero) |>pull(lambda_guerrero)aus_livestock |>filter(State =="Victoria", Animal =="Bulls, bullocks and steers") |>autoplot(box_cox(Count, lambda_livestock)) +labs(title ="Transformed Victoria slaughter counts",y =paste("lambda =", round(lambda_livestock, 2)))
Bullocks and steers show a sharp spike to about 130,000 in the 0’s followed by a long declining series with noisy, cycle like fluctuations that persists at both high and low levels up until 2020. The guerrero method selected a lambda of ~0.04 indicating that a transformation close to a natural log is a good fit.
Victorian Electricity Demand from vic_elec.
vic_elec |>autoplot(Demand) +labs(title ="Victorian 1/2 hour electricity demand", y ="MWh")
vic_elec |>features(Demand, features = guerrero) |>pull(lambda_guerrero)
This example show a very dense graph. with a broader seasonal pattern of taller demand spikes cluster around summer and winter months. The guerrero method selected a lambda of ~0.1. Because the graph is dense visually we cannot confirm a change in spread. On the raw scale, the maximum demand is about 2.02 times the median demand, while on the Box Cox transformed scale that ratio drops to about 1.13 times the median.
aus_production |>autoplot(Gas) +labs(title ="Australian quarterly gas production", y ="Petajoules")
lambda_gas <- aus_production |>features(Gas, features = guerrero) |>pull(lambda_guerrero)lambda_gas
[1] 0.1095171
The aus_production data set shows seasonal growths with seasonal changes. The guerrero method selected a lambda of ~0.11 matching what we had for reading on this chapter.
3.3
Why is a Box-Cox transformation unhelpful for the canadian_gas data?
A Box Cox transformation won’t really help with the canadian_gas data. The issue here isn’t that seasonal variation grows proportionally with the level, which is what a Box Cox transformation is actually designed to fix. Instead, the strength of the seasonal pattern changes over time somewhat independently of the trend. It’s stronger in the middle decades of the series and weaker toward the beginning and end. Since a Box Cox transformation applies the same fixed function at every time point, it can only stabilize variance that scales consistently with the series level. It just doesn’t have a mechanism to correct a seasonal pattern whose shape and strength evolve for reasons unrelated to the current level. Because of that, no choice of lambda is going to make the seasonality look consistent across the full series.
3.4
What Box-Cox transformation would you select for your retail data (from Exercise 7 in Section 2.10)?
set.seed(123456) #Using a random number given that I don't want my ID posted in Rpubs.myseries <- aus_retail |>filter(`Series ID`==sample(aus_retail$`Series ID`, 1))myseries
# A tsibble: 441 x 5 [1M]
# Key: State, Industry [1]
State Industry `Series ID` Month Turnover
<chr> <chr> <chr> <mth> <dbl>
1 Victoria Electrical and electronic goods retai… A3349564W 1982 Apr 55.8
2 Victoria Electrical and electronic goods retai… A3349564W 1982 May 58.4
3 Victoria Electrical and electronic goods retai… A3349564W 1982 Jun 53.7
4 Victoria Electrical and electronic goods retai… A3349564W 1982 Jul 56.9
5 Victoria Electrical and electronic goods retai… A3349564W 1982 Aug 56.7
6 Victoria Electrical and electronic goods retai… A3349564W 1982 Sep 58.9
7 Victoria Electrical and electronic goods retai… A3349564W 1982 Oct 59.6
8 Victoria Electrical and electronic goods retai… A3349564W 1982 Nov 63.2
9 Victoria Electrical and electronic goods retai… A3349564W 1982 Dec 90.3
10 Victoria Electrical and electronic goods retai… A3349564W 1983 Jan 55.5
# ℹ 431 more rows
myseries |>autoplot(Turnover) +labs(title ="Electrical and electronic goods retailing", y ="Turnover")
lambda_myseries <- myseries |>features(Turnover, features = guerrero) |>pull(lambda_guerrero)lambda_myseries
[1] -0.03391672
myseries |>autoplot(box_cox(Turnover, lambda_myseries)) +labs(title ="Transformed Victoria electrical and electronic goods retailing",y =paste("lambda =", round(lambda_myseries, 2)))
If we look at the retail series for Victoria, specifically the turnover for electrical and electronic goods, the raw data shows a clear December seasonal spike. This spike gets much larger in absolute size as the overall trend rises. It goes from a modest peak in the middle of the 1980s to a jump of several hundred million dollars by the late 2010s. The Guerrero method selected a lambda of ~-0.03, which is very close to a natural log transformation. If you compare individual years directly, the raw seasonal spread grows from 1985 to2018, which is roughly a six fold increase. On the transformed scale, however, the spread is nearly identical.
3.5
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.
aus_production |>autoplot(Tobacco) +labs(title ="Australian tobacco and cigarette production", y ="Tonnes")
Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_line()`).
lambda_tobacco <- aus_production |>features(Tobacco, features = guerrero) |>pull(lambda_guerrero)lambda_tobacco
Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_line()`).
For Tobacco production in aus_production, the guerrero method selected a lambda of ~0.93, very close to 1. Since a lambda of exactly 1 leaves the shape of a series unchanged aside from a shift, this near unity value indicates that a Box Cox transformation offers little practical benefit here.
ansett |>filter(Airports =="MEL-SYD", Class =="Economy") |>autoplot(Passengers) +labs(title ="Economy class passengers, Melbourne-Sydney")
lambda_ansett <- ansett |>filter(Airports =="MEL-SYD", Class =="Economy") |>features(Passengers, features = guerrero) |>pull(lambda_guerrero)lambda_ansett
[1] 1.999927
ansett |>filter(Airports =="MEL-SYD", Class =="Economy") |>autoplot(box_cox(Passengers, lambda_ansett)) +labs(title ="Transformed Economy class passengers",y =paste("lambda =", round(lambda_ansett, 2)))
ansett |>filter(Airports =="MEL-SYD", Class =="Economy", Passengers >0) |>features(Passengers, features = guerrero) |>pull(lambda_guerrero)
[1] 0.2199838
For Economy class passengers on the Melbourne to Sydney route in ansett, applying guerrero directly to the full series produces a lambda of approximately 2, which would stretch large values apart rather than compress them, an unusual result for this kind of transformation. This is driven by an extended stretch in late 1980s where passenger counts drop to exactly 0. Recalculating guerrero after excluding those zero weeks produces a much more typical lambda of approximately 0.22, close to a log transformation and consistent with the general purpose of Box Cox transformations
Hourly pedestrian counts at Southern Cross Station show a repeating daily commuting pattern, with counts spiking to around 3000 to 3800 during peak hours and dropping near 0 overnight, along with a visible data gap spanning late 2015 into early 2016. The guerrero method selected a lambda of ~-0.25, a fairly strong transformation. This shows the transformation aggressively compresses the upper range of the series, which may be a reasonable way to stabilize variance across the huge range between overnight and peak hour counts, but is also worth noting as a tradeoff, since it makes it harder to distinguish between a strong commuting peak and an unusually strong one after transformation
3.7
a. Plot the time series. Can you identify seasonal fluctuations and/or a trend-cycle?
gas <-tail(aus_production, 5*4) |>select(Gas)
gas |>autoplot(Gas) +labs(title ="Australian gas production, last 5 years", y ="Petajoules")
Plotting the last five years of gas production reveals a clear repeating seasonal pattern, with production dipping to a low point each Q1, rising to a peak each Q3, then falling again through Q4, forming a consistent shape across all five years.
b. Use classical_decomposition with type=multiplicative to calculate the trend-cycle and seasonal indices.
gas |>model(classical_decomposition(Gas, type ="multiplicative")) |>components() |>autoplot() +labs(title ="Classical multiplicative decomposition of gas production")
Warning: `autoplot.dcmp_ts()` was deprecated in fabletools 0.6.0.
ℹ Please use `ggtime::autoplot.dcmp_ts()` instead.
ℹ Graphics functions have been moved to the {ggtime} package. Please use
`library(ggtime)` instead.
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_line()`).
Using a multiplicative classical decomposition separates the gas series into three distinct components. The trend cycle rises steadily from about 200 at the start of the five year window to about 230 by the end, which confirms the mild upward drift we identified visually. The seasonal component is a fixed pattern that repeats identically every year, with the first quarter running about 12% below trend and the Q3 running about 12% above trend. The remainder component fluctuates close to one with no strong remaining pattern, though it is undefined for the first two and last two quarters of the series. This happens because the trend cycle itself cannot be estimated at the very ends of a short series using this method.
c. Do the results support the graphical interpretation from part a?
The results do support the graphical interpretation from part a. The decomposition formalizes and confirms the pattern that was already visible by eye in the original plot.
d. Compute and plot the seasonally adjusted data.
gas |>model(classical_decomposition(Gas, type ="multiplicative")) |>components() |>autoplot(season_adjust) +labs(title ="Seasonally adjusted gas production", y ="Petajoules")
Removing the seasonal component from the multiplicative decomposition produces a seasonally adjusted series that climbs from ~200 petajoules at the start of the 5 year window to ~230 petajoules near the end, following a much smoother path than the raw data. The sharp Q1 to Q3 sawtooth pattern visible in the original series is no longer present, leaving primarily the underlying trend along with the smaller remainder fluctuations.
e. Change one observation to be an outlier (e.g., add 300 to one observation), and recompute the seasonally adjusted data. What is the effect of the outlier?
gas_outlier_mid <- gas |>mutate(Gas =if_else(Quarter ==yearquarter("2007 Q4"), Gas +300, Gas))gas_outlier_mid |>autoplot(Gas) +labs(title ="Gas production with outlier added at 2007 Q4", y ="Petajoules")
gas_outlier_mid |>model(classical_decomposition(Gas, type ="multiplicative")) |>components() |>autoplot(season_adjust) +labs(title ="Seasonally adjusted data with outlier at 2007 Q4", y ="Petajoules")
Adding 300 to the 2007 Q4 observation creates a clear outlier in the raw series, rising to about 505 petajoules against a typical range
f. Does it make any difference if the outlier is near the end rather than in the middle of the time series?
gas_outlier_end <- gas |>mutate(Gas =if_else(Quarter ==yearquarter("2010 Q2"), Gas +300, Gas))gas_outlier_end |>model(classical_decomposition(Gas, type ="multiplicative")) |>components() |>autoplot(season_adjust) +labs(title ="Seasonally adjusted data with outlier at 2010 Q2", y ="Petajoules")
Placing the outlier at the final quarter of 2010 produces a different effect than placing it in the middle of the series. Because the moving average trend in classical decomposition cannot be estimated for the last two quarters of any series, the outlier at the final observation never enters the calculation that builds the seasonal index for each quarter. As a result, the seasonal component here remains just as clean and undistorted as it was in the original decomposition. However, this also means that no seasonal adjustment is effectively applied to that final point, so the outlier passes through into the seasonally adjusted series almost entirely intact, reaching close to 500.
3.8
Recall your retail time series data (from Exercise 7 in Section 2.10). Decompose the series using X-11. Does it reveal any outliers, or unusual features that you had not noticed previously?
library(seasonal)
x11_dcmp <- myseries |>model(x11 =X_13ARIMA_SEATS(Turnover ~x11())) |>components()x11_dcmp |>autoplot() +labs(title ="X-11 decomposition of Victoria electrical and electronic goods retailing")
x11_dcmp |>as_tibble() |>slice_max(abs(irregular -1), n =5) |>select(Month, Turnover, irregular)
# A tibble: 5 × 3
Month Turnover irregular
<mth> <dbl> <dbl>
1 2008 Dec 652. 1.14
2 2000 Apr 171 0.881
3 2011 Oct 368. 1.11
4 1986 Mar 72.9 0.891
5 2000 Jun 252. 1.10
The X-11 decomposition reveals two things the raw plot did not show. The trend shows a plateau and slight dip around 2000 rather than a smooth climb. The irregular component flags two outlier months, December 2008 running about 14% above expectation during the financial crisis, and April 2000 about 12% below.
3.9
Figures 3.19 and 3.20 show the result of decomposing the number of persons in the civilian labour force in Australia each month from February 1978 to August 1995.
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 trend and data share a scale of roughly 7000 to 9000 and rise steadily throughout, while the seasonal component is much smaller, only about negative 100 to 100. The remainder is dominated by one sharp dip near negative 350 around 1991, far larger than the usual noise elsewhere in that panel.
Is the recession of 1991/1992 visible in the estimated components?
Yes, the 1991 to 1992 recession is clearly visible, mainly as a sharp negative spike in the remainder reaching ~-350 around early 1990s, far outside the normal remainder range of ~-50 to 50.
