Much of the data that organisations keep is recorded over time: patients admitted each day, products sold each month, rainfall each week. Such a sequence of measurements taken at regular intervals is a time series, and it calls for a different way of thinking from the surveys of the previous chapters. The observations are not independent: this week’s value is related to last week’s, and the order of the observations carries information. The patterns are also of a particular kind. A series usually combines a slow trend, a seasonal pattern that repeats over a fixed period, and irregular noise, and understanding a series largely means separating these parts. Finally, the future of a series can never be known exactly, so a forecast is honestly a range of plausible values, not a single number.
This chapter develops these ideas and the methods built on them, using the fable family of packages. The example comes from outside the thesis. The doctors at the university counselling service have heard about Elaf’s wellbeing study and ask her to analyse their records: five years of weekly visit counts, with the service overwhelmed before exams and quiet in the summer, and staffing planned by guesswork. They want to know how many visits to expect next year, week by week. This is an internal analysis, done for the service rather than for publication, and a common situation for anyone who knows some data analysis: being asked to analyse someone else’s data.
Most of this book’s methods assume that observations are independent: one student’s answers tell you nothing about another’s. In a time series the opposite is true. The number of visits this week is closely related to the number last week, and the order of the observations is essential. The methods of this chapter are built around that dependence.
Two weeks of daily visits to a small clinic, starting on a Monday, show how a time series is stored:
library(dplyr)
library(ggplot2)
library(tsibble)
library(fable)
library(feasts)
clinic <- tsibble(
day = 1:14,
visits = c(12, 10, 9, 11, 7, 3, 2, 14, 11, 10, 12, 8, 4, 2),
index = day
)
clinic# A tsibble: 14 x 2 [1]
day visits
<int> <dbl>
1 1 12
2 2 10
3 3 9
4 4 11
5 5 7
6 6 3
7 7 2
8 8 14
9 9 11
10 10 10
11 11 12
12 12 8
13 13 4
14 14 2A tsibble (time series tibble) is a data frame that knows which column is time, its index. Here the index is the day number; the output header shows that the observations are one unit apart ([1]). tsibbles check that each time appears only once, and fable’s functions use the index to keep everything in order. There is a clear weekly pattern: busy at the start of the week, quiet at the weekend.
A time series is easiest to understand as a sum of parts. The code below builds an artificial weekly series for three years from three ingredients: a trend that rises slowly, a seasonal pattern that repeats every 52 weeks, and random noise. Figure 16.1 shows each part and their sum.
set.seed(4)
week <- 1:156
parts <- tibble(
week = week,
trend = 30 + 0.05 * week,
seasonal = 12 * sin(2 * pi * week / 52),
noise = rnorm(156, sd = 3)
) |>
mutate(series = trend + seasonal + noise)
parts |>
tidyr::pivot_longer(-week, names_to = "part", values_to = "value") |>
mutate(part = factor(part, levels = c("trend", "seasonal", "noise", "series"))) |>
ggplot(aes(x = week, y = value)) +
geom_line(colour = "#2f6793") +
facet_wrap(~ part, ncol = 1, scales = "free_y") +
labs(x = "Week", y = NULL) +
theme_minimal(base_size = 11)
In a real series, only the bottom panel is observed; the parts must be recovered from it. The trend says where the series is heading, the seasonal pattern says what happens at each point in the year, and the noise is what no pattern explains. Forecasting extends the trend and the seasonal pattern into the future. The noise cannot be forecast, and it is the main reason every forecast comes with a range. The decomposition later in this chapter recovers these parts from the counselling data.
The counselling records are in counselling_visits: the date of the Monday of each week, and the number of visits that week.
counselling_visits |> head() week_start visits
1 2020-09-07 32
2 2020-09-14 28
3 2020-09-21 24
4 2020-09-28 37
5 2020-10-05 32
6 2020-10-12 39range(counselling_visits$week_start)[1] "2020-09-07" "2025-08-25"The data covers five academic years, from September 2020 to August 2025, 260 weeks in all. The function yearweek() turns each date into a week, which becomes the index:
visits <- counselling_visits |>
mutate(week = yearweek(week_start)) |>
as_tsibble(index = week)
visits# A tsibble: 260 x 3 [1W]
week_start visits week
<date> <int> <week>
1 2020-09-07 32 2020 W37
2 2020-09-14 28 2020 W38
3 2020-09-21 24 2020 W39
4 2020-09-28 37 2020 W40
5 2020-10-05 32 2020 W41
6 2020-10-12 39 2020 W42
7 2020-10-19 31 2020 W43
8 2020-10-26 36 2020 W44
9 2020-11-02 37 2020 W45
10 2020-11-09 34 2020 W46
# ℹ 250 more rowsThe header now reads [1W]: one week between observations. The first step with any time series is to plot it, and autoplot() from the fable family draws a time plot:
autoplot(visits, visits) +
labs(x = NULL, y = "Visits per week") +
theme_minimal(base_size = 12)
Figure 16.2 shows the parts of a pattern described above, and one feature the doctors did not mention. There is a seasonal pattern that repeats every academic year: peaks before and during the two exam periods, dips during the breaks, and very few visits in the summer, when the service is reduced. There is a gentle upward trend, each year a little busier than the last, and there is noise, week-to-week variation that follows no pattern. Finally, there is one unusual week, in the autumn of 2023, far busier than the weeks around it.
A seasonal plot puts the years on top of each other, so the seasonal pattern and the differences between years are easier to see. The academic year starts in September, so the weeks are numbered from the start of each academic year:
visits |>
mutate(academic_year = paste0(20 + (row_number() - 1) %/% 52, "/", 21 + (row_number() - 1) %/% 52),
week_of_year = (row_number() - 1) %% 52 + 1) |>
ggplot(aes(x = week_of_year, y = visits, colour = academic_year)) +
geom_line() +
scale_colour_viridis_d(end = 0.9) +
labs(x = "Week of the academic year", y = "Visits per week", colour = "Academic year") +
theme_minimal(base_size = 12)
Every year has the same shape, and later years lie slightly above earlier ones. The period of the seasonal pattern is 52 weeks.
Dependence over time can be seen directly by plotting each week’s visits against the visits one week earlier, as in Figure 16.4.
visits |>
mutate(last_week = lag(visits)) |>
ggplot(aes(x = last_week, y = visits)) +
geom_point(alpha = 0.6, colour = "#2f6793") +
labs(x = "Visits in the previous week", y = "Visits this week") +
theme_minimal(base_size = 12)
The points rise from left to right: knowing last week’s visits says a good deal about this week’s. The autocorrelation at lag \(k\) puts a number on this relationship: it is the correlation between the series and itself \(k\) steps earlier, between each week and the week before (lag 1), two weeks before (lag 2), and so on. The function ACF() calculates it:
visits_acf <- visits |> ACF(visits, lag_max = 60)
visits_acf |> slice(c(1, 2, 26, 52))# A tsibble: 4 x 2 [1W]
lag acf
<cf_lag> <dbl>
1 1W 0.690
2 2W 0.572
3 26W -0.244
4 52W 0.684The autocorrelation is 0.69 at lag 1, as the scatter plot suggested, and 0.68 at lag 52: a week tends to be like the same week a year earlier. At lag 26 it is negative: half a year apart, a busy term week is often matched with a quiet summer week. Figure 16.5 shows all the lags. Strong autocorrelation at the seasonal lag is the signature of seasonality; if a series had no autocorrelation at all, there would be nothing to forecast beyond its average.
autoplot(visits_acf) + theme_minimal(base_size = 12)
Decomposition recovers the parts of a pattern from an observed series: a smooth trend, a repeating seasonal component, and the remainder, what is left over, which corresponds to the noise of Figure 16.1. STL (seasonal and trend decomposition using loess) is a flexible, widely used method. season(period = 52) tells it the length of the seasonal pattern, and robust = TRUE stops unusual weeks from distorting the trend and season:
visits_stl <- visits |>
model(STL(visits ~ season(period = 52), robust = TRUE)) |>
components()autoplot(visits_stl) + theme_minimal(base_size = 11)
In Figure 16.6 the series is the sum of the three components below it, just as in the artificial example. The trend rises from about 28 to about 38 visits a week, levelling off in the last year; the seasonal component repeats the academic year; and the remainder is mostly small noise. The largest remainder is the unusual week:
visits_stl |>
as_tibble() |>
slice_max(abs(remainder), n = 3) |>
select(week, visits, trend, season_52, remainder)# A tibble: 3 × 5
week visits trend season_52 remainder
<week> <int> <dbl> <dbl> <dbl>
1 2023 W43 71 35.2 8.08 27.7
2 2024 W48 32 38.4 12.1 -18.5
3 2024 W20 51 37.1 31.3 -17.5In the week of 23 October 2023, there were 71 visits, about 28 more than trend and season would predict. When asked, the doctors remember at once: it was the university’s mental health awareness week, with posters everywhere encouraging students to seek help. The remainder is where such events show up. It is a real week, not an error, so it is kept, and the robust decomposition makes sure it does not distort the forecasts.
Notice also that the seasonal swings grow slightly as the trend rises: the peaks get higher and the summer lows stay low. When the seasonal pattern grows in proportion to the level of the series, it is usual to model the logarithm of the series, on which the swings become constant. fable handles this automatically: write log(visits) in the model, and the forecasts are transformed back to visits.
Every forecasting method should be compared with simple benchmark methods, just as every predictive model in Chapters 11 to 13 was compared with a baseline. The three most common are the mean method, which forecasts the average of all past observations; the naive method, which forecasts the last observed value; and the seasonal naive method, which forecasts the value from the same season in the last cycle (the same day last week, or the same week last year). For the clinic example, with its weekly cycle of 7 days, they give:
clinic_fc <- clinic |>
model(
mean = MEAN(visits),
naive = NAIVE(visits),
snaive = SNAIVE(visits ~ lag(7))
) |>
forecast(h = 7)
clinic_fc |>
as_tibble() |>
select(.model, day, .mean) |>
tidyr::pivot_wider(names_from = .model, values_from = .mean)# A tibble: 7 × 4
day mean naive snaive
<dbl> <dbl> <dbl> <dbl>
1 15 8.21 2 14
2 16 8.21 2 11
3 17 8.21 2 10
4 18 8.21 2 12
5 19 8.21 2 8
6 20 8.21 2 4
7 21 8.21 2 2The function model() fits several models at once, each named; forecast(h = 7) forecasts 7 steps ahead; and .mean holds the point forecasts. The mean method forecasts 8.2 visits every day; the naive method repeats the last value, a quiet Sunday, for the whole week; and the seasonal naive method repeats last week’s pattern, which for this series is clearly the most sensible.
Two families of models go beyond the benchmarks. Exponential smoothing (ETS) forecasts with weighted averages of past observations, with weights that decrease the further back the observations are, so that recent weeks count most. ETS models can track a changing level, a trend, and a seasonal pattern, each estimated from the data; the name stands for the three components it models, error, trend, and season. ARIMA models describe how each observation depends on earlier observations and on earlier random shocks, and use the autocorrelation of the series to forecast it.
The functions ETS() and ARIMA() in fable choose the details of each model automatically, by a criterion like the BIC of Chapter 14. Both work best with short seasonal periods, such as 4 quarters, 7 days, or 12 months. A 52-week season is long, so two strategies are common for weekly data. The first is to decompose first: remove the seasonal pattern with STL, forecast the seasonally adjusted series (trend and remainder) with ETS, and add the seasonal pattern back, all of which decomposition_model() does in one step. The second uses Fourier terms: the seasonal pattern is described with a few smooth waves (sines and cosines) of different lengths, used as predictors in an ARIMA model; fourier(period = 52, K = 6) uses six pairs of waves.
As in machine learning, forecasts must be judged on data the model has not seen. For time series, the test data must come after the training data, because a forecast can only use the past. The models are trained on the first four academic years and tested on the fifth:
visits_train <- visits |> filter(week_start < as.Date("2024-09-01"))
visits_test <- visits |> filter(week_start >= as.Date("2024-09-01"))
visits_fit <- visits_train |>
model(
mean = MEAN(visits),
naive = NAIVE(visits),
snaive = SNAIVE(visits ~ lag(52)),
stl_ets = decomposition_model(
STL(log(visits) ~ season(period = 52), robust = TRUE),
ETS(season_adjust ~ season("N"))
),
fourier = ARIMA(log(visits) ~ fourier(period = 52, K = 6) + PDQ(0, 0, 0))
)In the STL model, ETS(season_adjust ~ season("N")) forecasts the seasonally adjusted series with no seasonal component of its own (“N” for none), because the seasonal pattern is added back from the decomposition. In the Fourier model, PDQ(0, 0, 0) tells ARIMA not to model the season itself, because the Fourier terms do that.
Each model forecasts the 52 weeks of the test year, and accuracy() compares the forecasts with what actually happened, using the RMSE and MAE of Chapter 13:
visits_fc <- visits_fit |> forecast(h = 52)
visits_fc |>
accuracy(visits) |>
select(.model, RMSE, MAE) |>
arrange(RMSE)# A tibble: 5 × 3
.model RMSE MAE
<chr> <dbl> <dbl>
1 stl_ets 7.40 5.57
2 snaive 9.16 6.96
3 fourier 10.3 7.81
4 mean 18.7 16.3
5 naive 27.5 23.5 The STL and ETS model is best: its forecasts are off by 5.6 visits a week on average (MAE), against 7.0 for the seasonal naive benchmark, which simply repeats the previous year. The mean and naive methods, which ignore the seasonal pattern, are far worse. The Fourier model does worse than the seasonal naive benchmark: the counselling service’s pattern has sharp steps (the summer drop happens from one week to the next), and smooth waves struggle to follow them. Figure 16.7 compares the two best models with the actual visits.
visits_fc |>
filter(.model %in% c("stl_ets", "snaive")) |>
autoplot(visits_test, level = NULL) +
labs(x = NULL, y = "Visits per week", colour = "Model") +
theme_minimal(base_size = 12)
The seasonal naive forecast copies last year’s noise along with its pattern: it even repeats the awareness week of October 2023 as a spike in October 2024. The STL and ETS model smooths the noise away and adds the trend, so its forecasts are both smoother and closer.
With the model chosen, it is refitted on all five years, so the forecasts use the most recent information, and the next 52 weeks are forecast:
final_fit <- visits |>
model(stl_ets = decomposition_model(
STL(log(visits) ~ season(period = 52), robust = TRUE),
ETS(season_adjust ~ season("N"))
))
next_year <- final_fit |> forecast(h = 52)next_year |>
autoplot(visits |> filter(week_start >= as.Date("2023-09-01"))) +
labs(x = NULL, y = "Visits per week") +
theme_minimal(base_size = 12)
A forecast is never a single number. The shaded bands in Figure 16.8 are prediction intervals: the ranges within which the actual visits are expected to fall with 80% and 95% probability. The function hilo() shows them as numbers:
next_year |>
hilo(95) |>
select(week, .mean, `95%`) |>
head(4)# A tsibble: 4 x 3 [1W]
week .mean `95%`
<week> <dbl> <hilo>
1 2025 W36 39.0 [23.60596, 60.76414]95
2 2025 W37 35.2 [21.32399, 54.89014]95
3 2025 W38 47.5 [28.74758, 73.99921]95
4 2025 W39 48.8 [29.58119, 76.14503]95In the first week of the new year, about 39 visits are expected, but anything from about 24 to 61 would not be surprising. The doctors also asked for the year as a whole. Adding up the weekly forecasts gives the expected total, but not its uncertainty, because weeks do not vary independently. The function generate() solves this by simulation: it produces 1,000 possible futures from the model, and the totals of those futures show the range of likely totals:
set.seed(2026)
futures <- final_fit |> generate(h = 52, times = 1000)
yearly_totals <- futures |>
as_tibble() |>
summarise(total = sum(.sim), .by = .rep)
quantile(yearly_totals$total, c(0.025, 0.5, 0.975)) |> round() 2.5% 50% 97.5%
1952 2111 2286 The report to the service therefore says: about 2,111 visits are expected in 2025/26 (95% interval 1,952 to 2,286), compared with 1,954 in 2024/25; the busiest weeks will again be just before and during the two exam periods, when they should plan the most staff; and an awareness campaign can raise demand sharply for a week, so extra capacity should be planned for any such event.
Every forecast in this chapter assumes that the patterns of the last five years continue. A change the data cannot know about, such as a new online booking system, a change in the exam calendar, or a crisis like the COVID-19 pandemic, can make any forecast wrong. Report forecasts with their intervals, say what they assume, and update them as new data arrives. (You may also meet Prophet, a forecasting tool from Meta that was popular for business data; it is no longer actively developed, and fable’s models are a well-supported alternative.)
Time series have their own traps, most of them versions of forgetting that time has an order.
generate() gives intervals for totals.Time series, tsibble, index, time plot, trend, seasonality, seasonal period, seasonal plot, autocorrelation, lag, decomposition, STL, remainder, seasonally adjusted series, benchmark forecast, mean method, naive method, seasonal naive method, exponential smoothing (ETS), ARIMA, Fourier terms, forecast horizon, prediction interval, simulation.
The playground has these and more, with hints and solutions.
clinic after adding a steady increase of one visit a day.season_adjust in visits_stl), and describe what it shows that the original series hides.K = 2 and K = 12, and explain how and why the test accuracy changes.generate() to estimate how many visits to expect in the four weeks before the first exam period of 2025/26, with a 95% interval.