library(feasts) library(tsibble) library(dplyr)
The feasts package provides a smorgasbord of tools for analysing tidy temporal data in the tsibble format. The package name is an acronym of its primary functionality: Feature Extraction and Statistics for Time Series.
Suppose we wanted to analyse seasonal patterns in the number of
domestic travellers to Melbourne, Australia. In the
tsibble::tourism data set, this can be further broken down
into 4 reasons of travel: “business”, “holiday”, “visiting friends and
relatives” and “other reasons”. The first observation from each series
are shown below.
<- tourism %>% tourism_melb filter(Region == "Melbourne") %>% tourism_melb group_by(Purpose) %>% slice(1) #> # A tsibble: 4 x 5 [1Q] #> # Key: Region, State, Purpose  #> # Groups: Purpose  #> Quarter Region State Purpose Trips #> <qtr> <chr> <chr> <chr> <dbl> #> 1 1998 Q1 Melbourne Victoria Business 405. #> 2 1998 Q1 Melbourne Victoria Holiday 428. #> 3 1998 Q1 Melbourne Victoria Other 79.9 #> 4 1998 Q1 Melbourne Victoria Visiting 666.
A useful first look at a time series is with a standard time series plot.
%>% tourism_melb autoplot(Trips)
This plot reveals that each of the reasons for visiting Melbourne follow a similar pattern, and that in recent years the amount of tourism is increasing. It is also clear that most people travel to Melbourne for one of three reasons (Business, Holiday or Visiting), and fewer travellers report other travel purposes.
While this plot reveals many useful large scale features of the data, it can be difficult to identify specifics about the seasonal pattern (other than its existence). Which quarter has the most tourism travel in Melbourne? To find out we will need to investigate other plot types.
%>% tourism_melb gg_season(Trips)
The seasonal plot (
gg_season()) wraps a seasonal period
(in this case, years) over the x axis, allowing you to see how each
quarter varies. In particular, it is apparent that Q3 is a low point for
people visiting friends and relatives, which noticeably increases in Q4.
Similarly Q2 and Q3 are the time periods with the most business
The trend of recent years can also be seen in the spread between the lines. The more recent years (purple/pink) are higher than the previous years. This also reveals that the drop in visiting tourism in Q3 was far less extreme in the most recent two years.
%>% tourism_melb gg_subseries(Trips)
An alternative visualisation of seasonal patterns is the subseries
gg_subseries()), which isolates seasonal periods into
separate plots. The blue lines indicate the average number of trips in
each quarter, making the increase in visiting tourism from Q3 to Q4 more
obvious. This plot style is especially useful in seeing how seasonality
changes over time. Focusing on the visiting tourism (last row of
facets), the number of tourists in Q3 and Q4 are increasing much more
than in Q1 and Q2 (suggesting that the trend may vary between
A look at the correlations in each series could reveal structures which are difficult to identify in the above plots.
%>% tourism_melb ACF(Trips) #> # A tsibble: 76 x 5 [1Q] #> # Key: Region, State, Purpose  #> Region State Purpose lag acf #> <chr> <chr> <chr> <cf_lag> <dbl> #> 1 Melbourne Victoria Business 1Q 0.487 #> 2 Melbourne Victoria Business 2Q 0.193 #> 3 Melbourne Victoria Business 3Q 0.331 #> 4 Melbourne Victoria Business 4Q 0.470 #> 5 Melbourne Victoria Business 5Q 0.219 #> 6 Melbourne Victoria Business 6Q 0.0750 #> 7 Melbourne Victoria Business 7Q 0.168 #> 8 Melbourne Victoria Business 8Q 0.321 #> 9 Melbourne Victoria Business 9Q 0.166 #> 10 Melbourne Victoria Business 10Q -0.0183 #> # ... with 66 more rows
The above code computes autocorrelations (
however it is also possible to compute partial autocorrelations
PACF()) and cross-correlations (
The tables given from these correlation functions also have a nice
autoplot() method, which will show the correlations along
with a threshold for significance (controllable with the
%>% tourism_melb ACF(Trips) %>% autoplot()
Another helpful strategy in investigating the patterns in a time series is to decompose it into components of interest. A useful decomposition for this is the STL decomposition, which allows you to extract multiple seasonal patterns with any seasonal period.
%>% tourism_melb model(STL(Trips ~ season(window = "periodic"))) %>% components() #> # A dable: 320 x 10 [1Q] #> # Key: Region, State, Purpose, .model  #> # : Trips = trend + season_year + remainder #> Region State Purpose .model Quarter Trips trend seaso~1 remai~2 seaso~3 #> <chr> <chr> <chr> <chr> <qtr> <dbl> <dbl> <dbl> <dbl> <dbl> #> 1 Melbourne Victoria Busine~ "STL(~ 1998 Q1 405. 437. -58.0 25.8 463. #> 2 Melbourne Victoria Busine~ "STL(~ 1998 Q2 408. 434. 18.1 -44.1 390. #> 3 Melbourne Victoria Busine~ "STL(~ 1998 Q3 486. 430. 37.9 18.3 448. #> 4 Melbourne Victoria Busine~ "STL(~ 1998 Q4 429. 431. 1.99 -3.48 427. #> 5 Melbourne Victoria Busine~ "STL(~ 1999 Q1 361. 427. -58.0 -7.82 419. #> 6 Melbourne Victoria Busine~ "STL(~ 1999 Q2 486. 409. 18.1 58.6 468. #> 7 Melbourne Victoria Busine~ "STL(~ 1999 Q3 359. 418. 37.9 -97.0 321. #> 8 Melbourne Victoria Busine~ "STL(~ 1999 Q4 426. 440. 1.99 -16.7 424. #> 9 Melbourne Victoria Busine~ "STL(~ 2000 Q1 495. 468. -58.0 84.1 552. #> 10 Melbourne Victoria Busine~ "STL(~ 2000 Q2 499. 501. 18.1 -19.8 481. #> # ... with 310 more rows, and abbreviated variable names 1: season_year, #> # 2: remainder, 3: season_adjust
The above call to
STL() has decomposed the
Trips variable into three components such that
Trips = trend + season_year + remainder. By setting
season(window = "periodic"), we have set the seasonal
pattern to be unchanging — you can control how quickly the seasonal
pattern can change by setting this to some number (smaller numbers
correspond to more rapid change).
%>% tourism_melb model(STL(Trips ~ season(window = 9))) %>% components() %>% autoplot()
Much like the table from
ACF(), decompositions can also
be plotted with
autoplot(). This gives the expected faceted
plot of the components extracted from the measured variable. The plot
shows that each purpose of travel has a different seasonal pattern, and
that the strength and structure of this pattern has changed over time.
As these components are often on different scales, this plot includes a
set of scale bars which are of equal scaled size across all plots.
The above plots and analysis are useful for if you’re looking at a
few series, but what can be done if you needed to look at and compare
more time series? Extracting features from a collection of time series
is a scalable approach to analysing many data sets. Each feature is a
numerical summary of the data set’s defining characteristics, and a set
of features can be created using
<- tourism_melb %>% tourism_melb_features features(Trips, feature_set(tags = "stl")) tourism_melb_features#> # A tibble: 4 x 12 #> Region State Purpose trend~1 seaso~2 seaso~3 seaso~4 spiki~5 linea~6 curva~7 #> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> #> 1 Melbour~ Vict~ Busine~ 0.782 0.559 3 1 515. 351. 326. #> 2 Melbour~ Vict~ Holiday 0.864 0.308 1 0 626. 793. 121. #> 3 Melbour~ Vict~ Other 0.826 0.183 3 0 11.8 200. 104. #> 4 Melbour~ Vict~ Visiti~ 0.759 0.475 0 3 932. 457. 320. #> # ... with 2 more variables: stl_e_acf1 <dbl>, stl_e_acf10 <dbl>, and #> # abbreviated variable names 1: trend_strength, 2: seasonal_strength_year, #> # 3: seasonal_peak_year, 4: seasonal_trough_year, 5: spikiness, 6: linearity, #> # 7: curvature
In the example above, the components from the STL decomposition has been used to summarise the strength of trend and seasonality components of each series.
These features are particularly useful to show on a plot.
library(ggplot2) %>% tourism_melb_features ggplot(aes(x = trend_strength, y = seasonal_strength_year, colour = Purpose)) + geom_point() + coord_equal() + lims(x = c(0,1), y = c(0,1))
When analysing just four series, a plot of the features does not look
very exciting. It is worth noting that a lot of the individuality seen
in the previous analysis have been lost when each series is summarised
down to just two values. However recall that the analysis has been
working from a very small subset of the complete
data set. Let’s see how Melbourne compares with the other regions in the
<- tourism %>% tourism_features features(Trips, feat_stl) ggplot(mapping = aes(x = trend_strength, y = seasonal_strength_year, colour = Purpose)) + geom_point(data = tourism_features, alpha = 0.3) + geom_point(data = tourism_melb_features, size = 2) + coord_equal() + facet_wrap(vars(Purpose)) + lims(x = c(0,1), y = c(0,1))
Looks like Melbourne is one of the trendiest places around Australia!
More information about time series analysis using the feasts package can be found in Forecasting: Principles and Practices (3rd Ed.) and in the pkgdown site.