ETF3231/5231: Business forecasting
Ch3. Time series decomposition OTexts.org/fpp3/
Outline
1 Transformations and adjustments
2 Time series components
3 Moving averages
4 Classical decomposition
5 History of time series decomposition
6 STL decomposition
2
Outline
1 Transformations and adjustments
2 Time series components
3 Moving averages
4 Classical decomposition
5 History of time series decomposition
6 STL decomposition
3
Per capita adjustments
global_economy %>%
filter(Country == “Australia”) %>%
autoplot(GDP)
1.5e+12
1.0e+12
5.0e+11
0.0e+00
1960 1980
2000
Year [1Y]
4
GDP
Per capita adjustments
global_economy %>%
filter(Country == “Australia”) %>%
autoplot(GDP / Population)
60000
40000
20000
0
1960 1980
2000
Year [1Y]
5
GDP/Population
Your turn
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?
6
Inflation adjustments
print_retail <- aus_retail %>%
filter(Industry == “Newspaper and book retailing”) %>%
group_by(Industry) %>%
index_by(Year = year(Month)) %>%
summarise(Turnover = sum(Turnover))
aus_economy <- filter(global_economy, Code == "AUS")
print_retail %>%
left_join(aus_economy, by = “Year”) %>%
mutate(Adj_turnover = Turnover / CPI) %>%
pivot_longer(c(Turnover, Adj_turnover),
names_to = “Type”, values_to = “Turnover”
) %>%
ggplot(aes(x = Year, y = Turnover)) +
geom_line() +
facet_grid(vars(Type), scales = “free_y”) +
xlab(“Years”) + ylab(NULL) +
ggtitle(“Turnover: Australian print media industry”)
7
Inflation adjustments
Turnover: Australian print media industry
50
45
40
35
30
4000
3000
2000
1990 2000 2010
Years
8
Adj_turnover Turnover
Mathematical transformations
If the data show different variation at different levels of the series, then a transformation can be useful.
9
Mathematical transformations
If the data show different variation at different levels of the series, then a transformation can be useful.
Denote original observations as y1, . . . , yT and transformed observations as w1, . . . , wT.
9
Mathematical transformations
If the data show different variation at different levels of the series, then a transformation can be useful.
Denote original observations as y1, . . . , yT and transformed observations as w1, . . . , wT.
Mathematical transformations for stabilizing variation
wt = √yt wt = √3 yt wt = log(yt)
Square root Cube root Logarithm
↓
Increasing strength
9
Mathematical transformations
If the data show different variation at different levels of the series, then a transformation can be useful.
Denote original observations as y1, . . . , yT and transformed observations as w1, . . . , wT.
Mathematical transformations for stabilizing variation
Square root Cube root Logarithm
wt = √yt wt = √3 yt wt = log(yt)
↓
Increasing strength
Logarithms, in particular, are useful because they are more interpretable: changes in a log value are relative (percent) changes on the original scale.
9
Mathematical transformations
food <- aus_retail %>%
filter(Industry == “Food retailing”) %>%
summarise(Turnover = sum(Turnover))
10000
5000
1990 2000 2010 2020
Month [1M]
10
Turnover ($AUD)
Mathematical transformations
food %>% autoplot(sqrt(Turnover)) +
labs(y = “Square root turnover”)
11
100
75
50
1990 2000 2010 2020
Month [1M]
Square root turnover
Mathematical transformations
food %>% autoplot(Turnoverˆ(1/3)) +
labs(y = “Cube root turnover”)
12
20
15
10
1990 2000 2010 2020
Month [1M]
Cube root turnover
Mathematical transformations
food %>% autoplot(log(Turnover)) +
labs(y = “Log turnover”)
13
9.5
9.0
8.5
8.0
7.5
7.0
1990
2000 2010 2020
Month [1M]
Log turnover
Mathematical transformations
food %>% autoplot(-1/Turnover) +
labs(y = “Inverse turnover”)
14
−0.00025
−0.00050
−0.00075
1990 2000 2010 2020
Month [1M]
Inverse turnover
Box-Cox transformations
Each of these transformations is close to a member of the family of Box-Cox transformations:
log(yt), λ = 0; (ytλ−1)/λ, λ̸=0.
wt =
15
Box-Cox transformations
Each of these transformations is close to a member of the family of Box-Cox transformations:
log(yt), λ = 0; (ytλ−1)/λ, λ̸=0.
wt =
λ = 1: (No substantive transformation)
λ = 1 : (Square root plus linear transformation) 2
λ = 0: (Natural logarithm) λ = −1: (Inverse plus 1)
15
Box-Cox transformations
16
Box-Cox transformations
food %>%
features(Turnover, features = guerrero)
## # A tibble: 1 x 1
## ## ## 1
lambda_guerrero
0.0524
17
Box-Cox transformations
food %>%
features(Turnover, features = guerrero)
## # A tibble: 1 x 1
## ## ## 1
lambda_guerrero
0.0524
This attempts to balance the seasonal fluctuations and random variation across the series.
Always check the results.
A low value of λ can give extremely large prediction intervals.
17
Box-Cox transformations
food %>% autoplot(box_cox(Turnover, 0.0524)) +
labs(y = “Box-Cox transformed turnover”)
18
12
11
10
9
1990
2000 2010 2020
Month [1M]
Box−Cox transformed turnover
Transformations
Often no transformation needed.
Simple transformations are easier to explain and work well enough.
19
Transformations
Often no transformation needed.
Simple transformations are easier to explain and work well enough.
Transformations can have very large effect on PI.
19
Transformations
Often no transformation needed.
Simple transformations are easier to explain and work well enough.
Transformations can have very large effect on PI.
If the data contains zeros, then can’t take logs. log1p() can be useful.
If some data are negative, no power transformation is possible unless a constant is added to all values.
19
Transformations
Often no transformation needed.
Simple transformations are easier to explain and work well enough.
Transformations can have very large effect on PI.
If the data contains zeros, then can’t take logs. log1p() can be useful.
If some data are negative, no power transformation is possible unless a constant is added to all values. Choosing logs is a simple way to force forecasts to be positive
Transformations must be reversed to obtain forecasts on the original scale. (Handled automatically by fable.)
19
Your turn
1 For the following series, find an appropriate transformation in order to stabilise the variance.
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
2 Why is a Box-Cox transformation unhelpful for the canadian_gas data?
20
Outline
1 Transformations and adjustments
2 Time series components
3 Moving averages
4 Classical decomposition
5 History of time series decomposition
6 STL decomposition
21
Time series patterns
Recall
Trend Cyclic
pattern exists when there is a long-term increase or decrease in the data.
pattern exists when data exhibit rises and falls that are not of fixed period (duration usually of at least 2 years).
pattern exists when a series is influenced by seasonal factors (e.g., the quarter of the year, the month, or day of the week).
Seasonal
22
Time series decomposition
where
yt = Tt = St = Rt =
yt =f(St,Tt,Rt)
data at period t
trend-cycle component at period t seasonal component at period t remainder component at period t
23
Time series decomposition
where yt = Tt = St = Rt =
yt =f(St,Tt,Rt)
data at period t
trend-cycle component at period t seasonal component at period t remainder component at period t
Additive decomposition: yt = St + Tt + Rt. Multiplicative decomposition: yt = St × Tt × Rt.
23
Time series decomposition
Additive model appropriate if magnitude of seasonal fluctuations does not vary with level.
If seasonal are proportional to level of series, then multiplicative model appropriate.
Multiplicative decomposition more prevalent with economic series
Alternative: use a Box-Cox transformation, and then use additive decomposition.
Logs turn multiplicative relationship into an additive relationship:
yt =St ×Tt ×Et ⇒ logyt =logSt +logTt +logRt.
24
US Retail Employment
us_retail_employment <- us_employment %>%
filter(year(Month) >= 1990, Title == “Retail Trade”) %>%
select(-Series_ID)
us_retail_employment
## # A tsibble: 357 x 3 [1M]
## Month Title
##
## 1 1990 Jan Retail Trade
## 2 1990 Feb Retail Trade
## 3 1990 Mar Retail Trade
## 4 1990 Apr Retail Trade
## 5 1990 May Retail Trade
## 6 1990 Jun Retail Trade
## 7 1990 Jul Retail Trade
## 8 1990 Aug Retail Trade
## 9 1990 Sep Retail Trade
## 10 1990 Oct Retail Trade
Employed
13256.
12966.
12938.
13012.
13108.
13183.
13170.
13160.
13113.
13185.
25
US Retail Employment
us_retail_employment %>%
autoplot(Employed) +
xlab(“Year”) + ylab(“Persons (thousands)”) +
ggtitle(“Total employment in US retail”)
16000
15000
14000
13000
Total employment in US retail
1990 2000
2010 2020
Year
26
Persons (thousands)
US Retail Employment
us_retail_employment %>%
model(stl = STL(Employed))
## # A mable: 1 x 1
## stl
##
## 1
27
US Retail Employment
dcmp <- us_retail_employment %>%
model(stl = STL(Employed))
components(dcmp)
## # A dable: 357 x 7 [1M]
## # Key: .model [1]
## # STL Decomposition: Employed = trend + season_year + remainder
## .model Month Employed trend season_year remainder
##
-38.1
-261.
-291.
-221.
-115.
-44.2 -23.0
0.0892
## 1 stl
## 2 stl
## 3 stl
## 4 stl
## 5 stl
## 6 stl
## 7 stl
## 8 stl
## 9 stl
1990 Jan 13256. 13291.
1990 Feb 12966. 13272.
1990 Mar 12938. 13252.
1990 Apr 13012. 13233.
1990 May 13108. 13213.
1990 Jun 13183. 13193.
1990 Jul 13170. 13173.
1990 Aug 13160. 13152.
1990 Sep 13113. 13131.
9.98
-25.6 15.7
-24.4 22.0
-11.8 19.5
-43.4 25.7
28
US Retail Employment
us_retail_employment %>%
autoplot(Employed, color=’gray’) +
autolayer(components(dcmp), trend, color=’red’) +
xlab(“Year”) + ylab(“Persons (thousands)”) +
ggtitle(“Total employment in US retail”)
29
Total employment in US retail
16000
15000
14000
13000
1990 2000
2010 2020
Year
Persons (thousands)
US Retail Employment
components(dcmp) %>% autoplot() + xlab(“Year”)
STL decomposition
Employed = trend + season_year + remainder
16000 15000 14000 13000 16000
15000 14000 13000
500 250 0 −250 100 50 0
−50 −100
1990 2000
2010 2020
Year
30
Employed trend season_year remainder
US Retail Employment
components(dcmp) %>% gg_subseries(season_year)
Jan Feb Mar
Apr May Jun Jul Aug Sep Oct
Nov Dec
500
250
0
−250
Month
31
stl
1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020 1990 2000 2010 2020
season_year
Seasonal adjustment
Useful by-product of decomposition: an easy way to calculate seasonally adjusted data. Additive decomposition: seasonally adjusted data given by
yt − St = Tt + Rt
Multiplicative decomposition: seasonally adjusted data given by
yt/St = Tt × Rt
32
US Retail Employment
us_retail_employment %>%
autoplot(Employed, color=’gray’) +
autolayer(components(dcmp), season_adjust, color=’blue’) +
xlab(“Year”) + ylab(“Persons (thousands)”) +
ggtitle(“Total employment in US retail”)
33
Total employment in US retail
16000
15000
14000
13000
1990 2000
2010 2020
Year
Persons (thousands)
Seasonal adjustment
We use estimates of S based on past values to seasonally adjust a current value.
Seasonally adjusted series reflect remainders as well as trend. Therefore they are not “smooth”” and “downturns”” or “upturns” can be misleading.
It is better to use the trend-cycle component to look for turning points.
34
Outline
1 Transformations and adjustments
2 Time series components
3 Moving averages
4 Classical decomposition
5 History of time series decomposition
6 STL decomposition
35
Moving averages
The simplest estimate of the trend-cycle uses moving averages.
mMA
ˆ 1 k
Tt = m
yt+j
where k = m−1 . 2
j=−k
36
Moving averages
global_economy %>% filter(Country == “Australia”) %>%
autoplot(Exports) + xlab(“Year”) + ylab(“% of GDP”) +
ggtitle(“Total Australian exports”)
Total Australian exports
21
18
15
12
1960 1980
2000
Year
37
% of GDP
Moving average smoothing
Country
Australia Australia Australia Australia Australia … Australia Australia Australia Australia Australia Australia
Year Exports
1960.00 12.99445 1961.00 12.40310 1962.00 13.94301 1963.00 13.00589 1964.00 14.93825 … … 2012.00 21.51897 2013.00 19.98772 2014.00 21.07577 2015.00 20.01296 2016.00 19.25303 2017.00 21.27035
5-MA
13.45694 13.50208 13.60794 … 20.77956 20.81365 20.36969 20.31997
38
Moving average smoothing
21
18
15
12
3−MA
1960 1980
2000
Year
39
Exports
Moving average smoothing
21
18
15
12
5−MA
1960 1980
2000
Year
40
Exports
Moving average smoothing
21
18
15
12
7−MA
1960 1980
2000
Year
41
Exports
Moving average smoothing
21
18
15
12
9−MA
1960 1980
2000
Year
42
Exports
Moving average smoothing
21
18
15
12
11−MA
1960 1980
2000
Year
43
Exports
Moving average smoothing
21
18
15
12
13−MA
1960 1980
2000
Year
44
Exports
Moving average smoothing
21
18
15
12
15−MA
1960 1980
2000
Year
45
Exports
Moving average smoothing
So a moving average is an average of nearby points
observations nearby in time are also likely to be close in value.
average eliminates some randomness in the data, leaving a smooth trend-cycle component.
3-MA: Tˆt = (yt−1 + yt + yt+1)/3
5-MA: Tˆt = (yt−2 + yt−1 + yt + yt+1 + yt+2)/5
each average computed by dropping oldest observation and including next observation. averaging moves through time series until trend-cycle computed at each observation possible.
46
Endpoints
Why is there no estimate at ends?
For a 3 MA, there cannot be estimates at time 1 or time T because the observations at time 0 and T + 1 are not available.
Generally: there cannot be estimates at times near the endpoints.
The order of the MA
larger order means smoother, flatter curve
larger order means more points lost at ends order = length of season or cycle removes pattern But so far odd orders?
47
Centered MA
4 MA:
1(yt−2 + yt−1 + yt + yt+1) 41
or 4(yt−1 + yt + yt+1 + yt+2)
48
Centered MA
4 MA:
1(yt−2 + yt−1 + yt + yt+1) 41
or 4(yt−1 + yt + yt+1 + yt+2) Solution: take a further 2-MA to “centre” result.
Tt = 11(yt−2 + yt−1 + yt + yt+1) 24
+ 1(yt−1 + yt + yt+1 + yt+2) 4
= 1yt−2 + 1yt−1 + 1yt + 1yt+1 + 1yt+2 84448 48
Centered MA
Year Data
1992 Q1 443.00
4-MA
2 × 4-MA
1992 Q2 410.00
451.25
1992 Q3 420.00
1992 Q4 532.00
1993 Q1 433.00
1993 Q2 421.00 1993 Q3 410.00 1993 Q4 512.00
450.00 450.12 450.25 446.50 446.00 443.00 .
.
448.75 451.50 449.00 444.00 448.00
438.00 . .
49
Moving average trend-cycle
A moving average of the same length as the season removes the seasonal pattern.
For quarterly data: use a 2 × 4 MA For monthly data: use a 2 × 12 MA
Tˆt = 1yt−6+ 1yt−5+···+ 1yt+5+ 1yt+6 24 12 12 24
50
Moving average trend-cycle
us_retail_employment_ma <- us_retail_employment %>%
mutate(
`12-MA` = slide_dbl(Employed, mean, .size = 12, .align = “cr”),
`2×12-MA` = slide_dbl(`12-MA`, mean, .size = 2, .align = “cl”) )
us_retail_employment_ma %>%
autoplot(Employed, color=’gray’) + autolayer(us_retail_employment_ma, vars(`2×12-MA`), color=’red’) + xlab(“Year”) + ylab(“Persons (thousands)”) +
ggtitle(“Total employment in US retail”)
51
Moving average trend-cycle
16000
15000
14000
13000
Total employment in US retail
1990 2000
2010 2020
Year
52
Persons (thousands)
Outline
1 Transformations and adjustments
2 Time series components
3 Moving averages
4 Classical decomposition
5 History of time series decomposition
6 STL decomposition
53
Trend-cycle
Additive decomposition: yt = Tt + St + Rt = Tˆt + Sˆt + Rˆt Multiplicative decomposition: yt = Tt × St × Rt = Tˆt × Sˆt × Rˆt
54
Trend-cycle
Additive decomposition: yt = Tt + St + Rt = Tˆt + Sˆt + Rˆt Multiplicative decomposition: yt = Tt × St × Rt = Tˆt × Sˆt × Rˆt
Estimate Tˆ using (2 × m)-MA if m is even. Otherwise, estimate Tˆ using m-MA
Compute de-trended series
Additive decomposition: yt − Tˆt
Multiplicative decomposition: yt/Tˆt De-trending
Remove smoothed series Tˆt from yt to leave St and Rt. Additivemodel:yt −Tˆt =(Tˆt +Sˆt +Rˆt)−Tˆt =Sˆt +Rˆt
Multiplicative model: yt = Tˆt×Sˆt×Rˆt = Sˆt × Rˆt 54 Tˆt Tˆt
Estimating seasonal component
Seasonal index for each quarter is estimated as an average of the detrended series for that quarter of successive years.
E.g., take averages across all Januaries to get S(1) if your data is monthly.
If necessary, adjust the seasonal indices so that:
foradditive:S(1) +S(2) +…+S(12) =0
for multiplicative: S(1) + S(2) + . . . + S(12) = m
The seasonal component Sˆt simply consists of replications of the seasonal indices.
55
Remainder component
Additive decomposition: Rˆt = yt − Tˆt − Sˆt Multiplicative decomposition: Rˆt = yt/(TˆtSˆt)
56
Remainder component
Additive decomposition: Rˆt = yt − Tˆt − Sˆt Multiplicative decomposition: Rˆt = yt/(TˆtSˆt)
Classical decomposition
Choose additive or multiplicative depending on which gives the most stable components.
For multiplicative model, this method of estimation is known as ratio-to-moving-average method.
56
US Retail Employment
Classical additive decomposition of total US retail employment Employed = trend + seasonal + random
16000 15000 14000 13000 16000
15000 14000 13000
600 400 200
0 −200
100 50 0 −50 −100
1990 2000
2010 2020
Year
57
Employed trend seasonal random
US Retail Employment
Classical multiplicative decomposition of total US retail employment Employed = trend * seasonal * random
16000 15000 14000 13000 16000
15000 14000 13000
1.04 1.02 1.00 0.98
1.005 1.000 0.995
1990 2000
2010 2020
Year
58
Employed trend seasonal random
US Retail Employment
us_retail_employment %>%
model(classical_decomposition(Employed, type = “additive”))
components() %>%
autoplot() + xlab(“Year”) +
ggtitle(“Classical additive decomposition of total
US retail employment”)
%>%
59
Comments on classical decomposition
Estimate of trend is unavailable for first few and last few observations.
Seasonal component repeats from year to year. May not be realistic.
Not robust to outliers.
Newer methods designed to overcome these problems.
60
Outline
1 Transformations and adjustments
2 Time series components
3 Moving averages
4 Classical decomposition
5 History of time series decomposition
6 STL decomposition
61
History of time series decomposition
Classical method originated in 1920s.
Census II method introduced in 1957. Basis for X-11 method and variants (including X-12-ARIMA, X-13-ARIMA)
STL method introduced in 1983
TRAMO/SEATS introduced in 1990s.
62
History of time series decomposition
Classical method originated in 1920s.
Census II method introduced in 1957. Basis for X-11 method and variants (including X-12-ARIMA, X-13-ARIMA)
STL method introduced in 1983
TRAMO/SEATS introduced in 1990s.
National Statistics Offices
ABS uses X-12-ARIMA
US Census Bureau uses X-13ARIMA-SEATS Statistics Canada uses X-12-ARIMA
ONS (UK) uses X-12-ARIMA
EuroStat use X-13ARIMA-SEATS
62
X-11 decomposition
Advantages
Relatively robust to outliers
Allows for slowly changing seasonal component Completely automated choices for trend and seasonal changes
Very widely tested on economic data over a long period of time.
63
X-11 decomposition
Advantages
Relatively robust to outliers
Allows for slowly changing seasonal component Completely automated choices for trend and seasonal changes
Very widely tested on economic data over a long period of time.
Disadvantages
No prediction/confidence intervals
Ad hoc method with no underlying model
Only developed for quarterly and monthly data 63
Extensions: X-12-ARIMA and X-13-ARIMA
The X-11, X-12-ARIMA and X-13-ARIMA methods are based on Census II decomposition.
These allow adjustments for trading days and other explanatory variables.
Known outliers can be omitted.
Level shifts and ramp effects can be modelled. Missing values estimated and replaced. Holiday factors (e.g., Easter, Labour Day) can be estimated.
64
X-13ARIMA-SEATS
Advantages
Model-based
Smooth trend estimate
Allows estimates at end points Allows changing seasonality Developed for economic data
65
X-13ARIMA-SEATS
Advantages
Model-based
Smooth trend estimate
Allows estimates at end points Allows changing seasonality Developed for economic data
Disadvantages
Only developed for quarterly and monthly data
65
Outline
1 Transformations and adjustments
2 Time series components
3 Moving averages
4 Classical decomposition
5 History of time series decomposition
6 STL decomposition
66
STL decomposition
STL: “Seasonal and Trend decomposition using Loess” Very versatile and robust.
Unlike X-12-ARIMA, STL will handle any type of seasonality.
Seasonal component allowed to change over time, and rate of change controlled by user. Smoothness of trend-cycle also controlled by user. Robust to outliers
Not trading day or calendar adjustments. Only additive.
Take logs to get multiplicative decomposition. Use Box-Cox transformations to get other decompositions.
67
STL decomposition
us_retail_employment %>%
model(STL(Employed ~ season(window=9), robust=TRUE)) %>%
components() %>% autoplot() +
ggtitle(“STL decomposition: US retail employment”)
STL decomposition: US retail employment
Employed = trend + season_year + remainder
16000 15000 14000 13000 16000 15000 14000 13000
500 250 0 −250
100 0 −100
1990 2000
2010 2020
Month
68
Employed trend season_year remainder
STL decomposition
69
STL decomposition
us_retail_employment %>%
model(STL(Employed ~ season(window=5))) %>%
components()
us_retail_employment %>%
model(STL(Employed ~ trend(window=15) +
season(window=”periodic”),
robust = TRUE)
) %>% components()
trend(window = ?)controlswigglinessoftrend component.
season(window = ?)controlsvariationonseasonal component.
season(window = ‘periodic’)isequivalenttoan infinite window.
70
STL decomposition
us_retail_employment %>%
model(STL(Employed)) %>%
components() %>%
autoplot()
STL decomposition
Employed = trend + season_year + remainder
16000 15000 14000 13000 16000 15000 14000 13000
500 250 0 −250 100 50 0
−50 −100
1990 2000
2010 2020
Month
71
Employed trend season_year remainder
STL decomposition
Algorithm that updates trend and seasonal components iteratively.
Starts with Tˆt = 0
Uses a mixture of loess and moving averages to successively refine the trend and seasonal estimates. The trend window controls loess bandwidth applied to deasonalised values.
The season window controls loess bandwidth applied to detrended subseries.
Robustness weights based on remainder. Defaultseasonwindow = 13
Defaulttrendwindow = nextodd( ceiling((1.5*period)/(1-(1.5/s.window))) 72