Navigate
Courses Tutorials Practice Pricing Tools Upgrade

EDA for time series

Today let's understand how to look at a time series properly, before you ever fit a model to it, using one real dataset.

Here is 37 years of real data: the monthly turnover of cafes, restaurants and takeaway food outlets in Victoria, Australia, from April 1982 to December 2018, in $ million. That's 441 months in a row, published by the Australian Bureau of Statistics, plotted once, in full.

Look at that shape for a second. It climbs for nearly four decades, but it is not a smooth climb: small bumps ride on top of it every single year. Everything in this lesson is about reading that shape properly, one piece at a time, before fitting anything to it.

The plots-first workflow and autoplot()

Before you fit anything to a time series, look at it properly first. That is exploratory data analysis, EDA for short, applied to a time series: a short, repeatable checklist you run before any model touches the data.

For a time series, that checklist covers six things, always in this order:

  1. Shape - what does the series look like, overall?
  2. Season-and-year summaries - how does it move within a year, and across years?
  3. Outliers - which points do not fit the pattern the rest of the series follows?
  4. Calendar effects - do those "outliers" actually repeat every year, at the same point in the calendar?
  5. Structural breaks - did the whole series shift level at some point?
  6. Missing runs - are there gaps in the data itself?

You already saw the first item on the cover: the raw shape of the series. To make that same plot in R, use autoplot(), the plotting function built for a tsibble that skips ggplot2's usual setup. Build the tsibble first.

RInteractive R
# Build the Victoria cafes and restaurants tsibble from the ABS retail data library(tsibble) library(tsibbledata) library(dplyr) options(pillar.sigfig = 6) cafe <- tsibbledata::aus_retail %>% filter(State == "Victoria", Industry == "Cafes, restaurants and takeaway food services") %>% select(Month, Turnover) cafe #> # A tsibble: 441 x 2 [1M] #> Month Turnover #> <mth> <dbl> #> 1 1982 Apr 85.1 #> 2 1982 May 85.1 #> 3 1982 Jun 82.8 #> 4 1982 Jul 82.1 #> 5 1982 Aug 81.8 #> 6 1982 Sep 84.6 #> 7 1982 Oct 91.7 #> 8 1982 Nov 97.7 #> 9 1982 Dec 109.3 #> 10 1983 Jan 94.6 #> # i 431 more rows

  

  





aus_retail is a large keyed tsibble bundled with the tsibbledata package: one series per Australian state and industry. Filtering down to Victoria's cafes, restaurants and takeaway food services leaves one series, 441 months, from April 1982 to December 2018. [1M] in the header says its index moves in fixed 1-month steps.



Now plot it with autoplot().




  

    
RInteractive R

    
# The tsibble-aware plot: no ggplot() setup needed
library(feasts)

autoplot(cafe, Turnover)

    

      
      
    

    

  

  





autoplot() reads the tsibble's index directly, so you only name the column to plot, Turnover, and it draws the line correctly ordered by Month without you writing a single ggplot() call. The plot is the same one from the cover: a level that rises for nearly four decades, an unmistakable trend, with a small wave riding on top of it every year.







Summaries by season and year



Shape tells you the series rises. It does not tell you which months are typically strong and which are typically weak, or whether that pattern holds steady from one year to the next. For that, summarise the series by season and by year.



Pivot the tsibble so each row is a calendar month and each column is a year, for the last three years in the data.




  

    
RInteractive R

    
# Pivot the tsibble to a Year x Month table for 2016-2018
library(tidyr)
library(lubridate)

season_wide <- cafe %>%
  as_tibble() %>%
  mutate(Year = year(Month), Mon = month(Month, label = TRUE, abbr = TRUE)) %>%
  filter(Year %in% 2016:2018) %>%
  select(Year, Mon, Turnover) %>%
  pivot_wider(names_from = Year, values_from = Turnover)

season_wide
#> # A tibble: 12 x 4
#>    Mon   `2016` `2017` `2018`
#>    <ord>  <dbl>  <dbl>  <dbl>
#>  1 Jan    790.8  861.4  887.3
#>  2 Feb    732.4  753.6  810.6
#>  3 Mar    793.3  868.4  909.7
#>  4 Apr    805.4  859.7  898.6
#>  5 May    808.9  845.2  866.4
#>  6 Jun    769    828.7  873.1
#>  7 Jul    828.7  857.6  916.1
#>  8 Aug    866.3  874.1  943.8
#>  9 Sep    878.3  872.4  933.8
#> 10 Oct    886.3  917.7  960
#> 11 Nov    879.5  925.7  981.8
#> 12 Dec    961.8 1008.1 1066.2

    

      
      
    

    

  

  





Read down any one column and turnover climbs from a January dip to a December peak, the same shape every year. Read across any one row and the number grows a little every year, the trend you already saw on the plot. December 2018, at $1,066.2 million, is the single highest month in this slice; February 2016, at $732.4 million, is the lowest.



Three years hints at a pattern. Thirty-seven years confirms it. Average each calendar month over the full 1982-2018 history and rank the months from strongest to weakest.




  

    
RInteractive R

    
# The 37-year average for each calendar month, highest to lowest
month_avg <- cafe %>%
  as_tibble() %>%
  mutate(Mon = month(Month, label = TRUE, abbr = TRUE)) %>%
  group_by(Mon) %>%
  summarise(avg_turnover = round(mean(Turnover), 1), .groups = "drop") %>%
  arrange(desc(avg_turnover))

month_avg
#> # A tibble: 12 x 2
#>    Mon   avg_turnover
#>    <ord>        <dbl>
#>  1 Dec          461.4
#>  2 Nov          416
#>  3 Oct          412
#>  4 Mar          405
#>  5 Aug          404.4
#>  6 Sep          400.7
#>  7 Jul          396.8
#>  8 Jan          396.6
#>  9 May          390.2
#> 10 Apr          389
#> 11 Jun          376.9
#> 12 Feb          363.5

    

      
      
    

    

  

  





December averages $461.4 million across all 37 years, the highest of any month by a wide margin. February averages $363.5 million, the lowest. That is the seasonal pattern this series carries every single year: a December peak, driven by Christmas trading, and a February trough, a shorter month with fewer trading days. Here is the three-year slice again, as a report-ready table.









Spotting outliers with the STL remainder



Season-and-year summaries tell you the typical pattern. To find the months that broke it, you need to strip that pattern away and look at what is left over. STL, short for Seasonal-Trend decomposition using Loess, does exactly that: it splits a series into a trend (the slow-moving overall level), a season component (the repeating yearly wave) and a remainder (whatever is left once both are subtracted out).



Fit STL on cafe, with season(window = "periodic") telling it to use one fixed seasonal shape for the whole 37 years, rather than letting that shape drift year to year.




  

    
RInteractive R

    
# Split the series into trend, yearly season and remainder with STL
dcmp <- cafe %>%
  model(STL(Turnover ~ season(window = "periodic"))) %>%
  components()

dcmp %>%
  as_tibble() %>%
  select(Month, Turnover, trend, season_year, remainder) %>%
  mutate(across(trend:remainder, ~ round(.x, 1)))
#> # A tibble: 441 x 5
#>       Month Turnover trend season_year remainder
#>       <mth>    <dbl> <dbl>       <dbl>     <dbl>
#>  1 1982 Apr     85.1  88.8        -5.4       1.8
#>  2 1982 May     85.1  88.8        -6.1       2.5
#>  3 1982 Jun     82.8  88.7       -21.5      15.5
#>  4 1982 Jul     82.1  88.9        -3.3      -3.5
#>  5 1982 Aug     81.8    89         2.6      -9.8
#>  6 1982 Sep     84.6  89.2        -2.9      -1.7
#>  7 1982 Oct     91.7  89.5         6.6      -4.4
#>  8 1982 Nov     97.7  89.6         8.5      -0.4
#>  9 1982 Dec    109.3  89.8        51.8     -32.3
#> 10 1983 Jan     94.6  90.2         0.3       4.1
#> # i 431 more rows

    

      
      
    

    

  

  





Add up trend, season_year and remainder for any row and you get back Turnover exactly. trend is a smooth, slow-moving line. season_year is the same fixed December-up, February-down wave every year. remainder is what is left over: how far off trend + season_year sits from the real number, for that one month.



Most months have a small remainder, a few $million either way. A month whose remainder is unusually large, relative to the rest, is a candidate outlier. "Unusually large" needs a threshold, and standard deviation gives you a simple one: flag any month whose remainder sits more than 2 standard deviations from zero.




  

    
RInteractive R

    
# Flag any month whose remainder sits more than 2 standard deviations from zero
remainder_sd <- sd(dcmp$remainder)
remainder_sd
#> [1] 16.9647

flagged <- dcmp %>% filter(abs(remainder) > 2 * remainder_sd)
nrow(flagged)
#> [1] 29

flagged %>%
  as_tibble() %>%
  select(Month, remainder) %>%
  mutate(remainder = round(remainder, 1))
#> # A tibble: 29 x 2
#>       Month remainder
#>       <mth>     <dbl>
#>  1 1984 Dec     -35.9
#>  2 1985 Dec     -34.2
#>  3 1988 Dec     -36.1
#>  4 1998 Dec     -34.2
#>  5 2003 Jan      36.7
#>  6 2008 Sep     -34.4
#>  7 2009 Jan      43.4
#>  8 2010 Aug      41.4
#>  9 2010 Sep      34.4
#> 10 2010 Dec      54.9
#> # i 19 more rows

    

      
      
    

    

  

  





The remainder's standard deviation across all 441 months is about 17. Twice that is 34, so any month whose remainder is more than 34 $million above or below zero gets flagged. That rule catches 29 of the 441 months, about 1 in 15. Here is the full series again, with those 29 months picked out.




Switch that widget to points and the 29 flagged months stand out clearly against the rest. But a flag is only a candidate. The next question is whether each one is a genuine surprise, or something the calendar explains perfectly well.







Quick check: the table and the flagging rule









Right on both counts. December leads the table in every one of the three years, and 40 is bigger than the 34 threshold, so it gets flagged.
December really is the highest month in the table for all three years, that part is right. But the threshold check is wrong: 2 x 17 = 34, and 25 is below 34, so a remainder of 25 would not be flagged. Only a remainder past 34, like 40, would.









Calendar effects: when a big residual isn't an outlier



29 flagged months is a starting list, not a final answer. Before you call any of them a genuine surprise, check whether it repeats. A remainder that comes back at the same calendar position, year after year, is not one-off at all: it is the fixed seasonal shape season(window = "periodic") assumed failing to fit that particular month exactly. That is a calendar effect, not an outlier.



Count the 29 flagged months by which calendar month they fall on.




  

    
RInteractive R

    
# Count how many flagged months fall on each calendar month
library(lubridate)

flagged_named <- flagged %>%
  as_tibble() %>%
  mutate(Mon = month(Month, label = TRUE, abbr = TRUE))

flagged_named %>% count(Mon, sort = TRUE)
#> # A tibble: 7 x 2
#>   Mon       n
#>   <ord> <int>
#> 1 Dec      11
#> 2 Feb       8
#> 3 Sep       4
#> 4 Jan       3
#> 5 May       1
#> 6 Jun       1
#> 7 Aug       1

    

      
      
    

    

  

  





Two months dominate the list: December, 11 of the 29 flags, and February, 8 of them. Together that is 19 of the 29, two out of every three flagged months, all landing on the same two calendar positions. Look at each one closely, with its year and its sign.




  

    
RInteractive R

    
# Look at every flagged February and December, with its year and sign
flagged_named %>%
  mutate(Year = year(Month),
         sign = ifelse(remainder > 0, "positive", "negative"),
         remainder = round(remainder, 1)) %>%
  filter(Mon %in% c("Feb", "Dec")) %>%
  arrange(Mon, Year) %>%
  select(Month, Year, remainder, sign)
#> # A tibble: 19 x 4
#>       Month  Year remainder sign
#>       <mth> <dbl>     <dbl> <chr>
#>  1 2011 Feb  2011     -62.9 negative
#>  2 2012 Feb  2012     -39.8 negative
#>  3 2013 Feb  2013     -35.1 negative
#>  4 2014 Feb  2014       -46 negative
#>  5 2015 Feb  2015     -40.8 negative
#>  6 2016 Feb  2016     -40.3 negative
#>  7 2017 Feb  2017     -72.6 negative
#>  8 2018 Feb  2018     -53.8 negative
#>  9 1984 Dec  1984     -35.9 negative
#> 10 1985 Dec  1985     -34.2 negative
#> 11 1988 Dec  1988     -36.1 negative
#> 12 1998 Dec  1998     -34.2 negative
#> 13 2010 Dec  2010      54.9 positive
#> 14 2011 Dec  2011      36.7 positive
#> 15 2013 Dec  2013      41.7 positive
#> 16 2015 Dec  2015        45 positive
#> 17 2016 Dec  2016      51.4 positive
#> 18 2017 Dec  2017      63.9 positive
#> 19 2018 Dec  2018      41.7 positive

    

      
      
    

    

  

  





Every one of the 8 flagged Februaries is negative, and every one of them falls in 2011 or later. February is a short month with fewer trading days than a 30 or 31-day month, so its total turnover runs low even when daily spending has not changed at all. As the whole series grew larger in later decades, that same February shortfall turned into a bigger dollar gap, big enough to clear the 34 threshold from 2011 on.



December tells a two-part story. The 4 flagged Decembers from 1984 to 1998 are all negative: the fixed seasonal shape from season(window = "periodic") was set higher than the smaller 1980s and 90s economy actually delivered. The 7 flagged Decembers from 2010 on are all positive: Christmas trading grew to outgrow that same fixed shape. Same calendar month, opposite sign, two different eras.



A flag that recurs at the same calendar position most years, in a direction the calendar itself explains, is a calendar effect. That leaves 10 of the original 29 flags, the ones scattered across September, January, May, June and August, as points still worth a closer look for a genuine, one-off surprise.







Structural breaks: a level shift the trend doesn't explain



Calendar effects repeat every year. A structural break is different: a point where the series shifts to a new level and stays there, something the smooth trend component was not built to catch in one sharp move.



Compute the mean turnover for each of the 37 years, then the percentage change from one year to the next.




  

    
RInteractive R

    
# Year-over-year percent change in the yearly mean turnover
yearly <- cafe %>%
  as_tibble() %>%
  mutate(Year = year(Month)) %>%
  group_by(Year) %>%
  summarise(mean_turnover = round(mean(Turnover), 1), .groups = "drop") %>%
  arrange(Year) %>%
  mutate(pct_change = round((mean_turnover / lag(mean_turnover) - 1) * 100, 1))

yearly
#> # A tibble: 37 x 3
#>     Year mean_turnover pct_change
#>    <dbl>         <dbl>      <dbl>
#>  1  1982          88.9       NA
#>  2  1983          91.3        2.7
#>  3  1984          96.4        5.6
#>  4  1985         103.2        7.1
#>  5  1986         124.4       20.5
#>  6  1987         147.7       18.7
#>  7  1988         160.4        8.6
#>  8  1989         187.3       16.8
#>  9  1990         189.7        1.3
#> 10  1991         183.2       -3.4
#> # i 27 more rows

yearly %>% arrange(desc(pct_change)) %>% head(3)
#> # A tibble: 3 x 3
#>    Year mean_turnover pct_change
#>   <dbl>         <dbl>      <dbl>
#> 1  1999         313.5       26.6
#> 2  1986         124.4       20.5
#> 3  1987         147.7       18.7

    

      
      
    

    

  

  





lag(mean_turnover) shifts the column down by one row, so subtracting it from the current row's value and dividing gives each year's change from the year before. Most years land somewhere between roughly 1% and 20%. One year sits well clear of all the others: 1999, up 26.6% on 1998, the single biggest year-over-year jump anywhere in this 37-year series.



A jump that size, that does not repeat and does not fade back, is a structural break: the series moved to a new level around 1999 and stayed there. At this stage in the EDA-for-time-series workflow, the job is to flag that break for whoever builds the model next, not to explain it. Whatever caused it, a chain expanding, a new competitor closing, a change in how the ABS measured the category, is a separate investigation. What matters here is that any model fit across 1999 has to account for that level shift directly, or the jump ends up folded into ordinary trend and season terms that were never built to explain it.







Missing runs: finding gaps before you model



The last item on the checklist has nothing to do with unusual values. It is about rows that are not there at all. Here is a second, much smaller tsibble to show it clearly: a shop's daily sales for the first 17 days of January 2024, with 3 days missing because the shop closed for a stocktake.




  

    
RInteractive R

    
# 14 days of shop sales out of 17, with a 3-day stocktake gap
library(tsibble)

shop <- tsibble(
  date = as.Date("2024-01-01") + c(0:6, 10:16),
  sales = c(41, 38, 45, 52, 49, 33, 29, 47, 50, 44, 39, 36, 31, 46),
  index = date
)

shop
#> # A tsibble: 14 x 2 [1D]
#>    date       sales
#>    <date>     <dbl>
#>  1 2024-01-01    41
#>  2 2024-01-02    38
#>  3 2024-01-03    45
#>  4 2024-01-04    52
#>  5 2024-01-05    49
#>  6 2024-01-06    33
#>  7 2024-01-07    29
#>  8 2024-01-11    47
#>  9 2024-01-12    50
#> 10 2024-01-13    44
#> 11 2024-01-14    39
#> 12 2024-01-15    36
#> 13 2024-01-16    31
#> 14 2024-01-17    46

    

      
      
    

    

  

  





14 rows, not 17. Nothing in that print shouts "gap" on its own; you would have to notice the date column jumps from January 7 straight to January 11. scan_gaps() and count_gaps() do that check for you.




  

    
RInteractive R

    
# Find the missing dates, then summarise them as one run
scan_gaps(shop)
#> # A tsibble: 3 x 1 [1D]
#>   date
#>   <date>
#> 1 2024-01-08
#> 2 2024-01-09
#> 3 2024-01-10

count_gaps(shop)
#> # A tibble: 1 x 3
#>   .from      .to           .n
#>   <date>     <date>     <int>
#> 1 2024-01-08 2024-01-10     3

    

      
      
    

    

  

  





scan_gaps() lists every missing date, one row at a time: January 8, 9 and 10. count_gaps() groups a run of consecutive missing dates into a single summary row instead: .from and .to mark where the run starts and ends, and .n counts how many dates fall inside it, 3 in this case, the exact 3-day stocktake closure.



Neither function changes shop. If you want the missing dates turned into real rows instead of just reported, fill_gaps(shop) does that: it inserts a row for each of the 3 missing dates, with sales set to NA, taking shop from 14 rows to 17. Report first with scan_gaps() and count_gaps(), repair only if the next step in your workflow actually needs a complete series.







Quick check: outlier, calendar effect, or structural break?



A colleague looks at three flagged points in their own retail series and asks you to sort them. The first is a month that comes in far below its neighbours every single year, always at the same point in the calendar. The second is one year where the yearly mean jumps far more than in any other year, and it never happens again. The third is a flagged month with no repeating pattern before or after it, just one unexplained spike.









Right. A flag that repeats every year at the same calendar position is a calendar effect. A one-off jump in the yearly mean that never repeats is a structural break. A flagged point with no repeating pattern around it is a genuine outlier.
The one that repeats every year at the same calendar position is the calendar effect, the same pattern February and December showed. A single jump in the yearly mean that never happens again is the structural break, the same shape 1999 showed. And a flagged point with nothing repeating around it, before or after, is the genuine outlier.









Your turn: find the gap and flag the point



Here is a second small series to practice both rules on: a helpdesk's daily ticket count for the 42 days from March 1 to April 11, 2024, with one day missing.




Right: the missing day is March 19, and the only day that clears 2 x 9.05 = 18.1 is March 30, a count of 95 against a remainder of about 38.2.Two calls: count_gaps(tickets) to find the missing date, then filter(tickets_dcmp, abs(remainder) > 2 * remainder_sd) to flag the oversized day.
Show answer
# Find the missing date, then flag the oversized day
count_gaps(tickets)
#> # A tibble: 1 x 3
#>   .from      .to           .n
#>   <date>     <date>     <int>
#> 1 2024-03-19 2024-03-19     1

filter(tickets_dcmp, abs(remainder) > 2 * remainder_sd)
#> # A tibble: 1 x 3
#>   date       count remainder
#>   <date>     <dbl>     <dbl>
#> 1 2024-03-30    95      38.2




Before your turn starts, here is the remainder column already computed for you, the same way you built it in the outlier step: fill the one missing day with the average of its two neighbours, so STL has a complete series to decompose, then run STL exactly as before.




  

    
RInteractive R

    
# Fill the missing day, decompose the series, and read off its remainder sd
library(feasts)

tickets_filled <- fill_gaps(tickets)
gap_row <- which(is.na(tickets_filled$count))
tickets_filled$count[gap_row] <- mean(tickets_filled$count[c(gap_row - 1, gap_row + 1)])

tickets_dcmp <- tickets_filled %>%
  model(STL(count ~ season(window = "periodic"))) %>%
  components() %>%
  as_tibble() %>%
  select(date, count, remainder)

remainder_sd <- sd(tickets_dcmp$remainder)
remainder_sd
#> [1] 9.047743

    

      
      
    

    

  

  





That's a remainder sd of about 9.05. Now it's your turn: find the missing date in tickets with count_gaps(), and flag the oversized day in tickets_dcmp with the same 2 times sd rule you used earlier, abs(remainder) > 2 * remainder_sd.




  

    
RInteractive R

    
# tickets still holds the original 41-day series, with one date missing.
# tickets_dcmp holds the filled, decomposed version, with a remainder column,
# and remainder_sd already holds its standard deviation.
# 1. Find the missing date with count_gaps(tickets).
# 2. Flag the oversized day: filter(tickets_dcmp, abs(remainder) > 2 * remainder_sd).
# Press Check when you have both.

    

      
      
    

    

  

  





March 19 is the missing day, the same one count_gaps() reported earlier. March 30 is the only day whose remainder clears 18.1: a count of 95 tickets against a typical day nearby, well outside anything the rest of the series produced.







References



Before you go, here is where the ideas in this lesson come from.









Putting the EDA-for-time-series workflow together



Run back through the six-item checklist against what this one series actually showed you. Shape: a 37-year climb with a wave riding on top of it every year. Season-and-year summaries: December is the strongest month every year, averaging $461.4 million, and February the weakest, averaging $363.5 million. Outliers: an STL remainder with a standard deviation of about 17 flagged 29 of the 441 months. Calendar effects: 19 of those 29 flags, 8 Februaries and 11 Decembers, turned out to be the fixed seasonal shape missing a short month or an outgrown Christmas peak, not real surprises. Structural breaks: 1999 stood out as a 26.6% jump in the yearly mean, a level shift to flag for whoever models this series next, not to explain away. Missing runs: scan_gaps() and count_gaps() found a clean 3-day gap in the shop series, without needing to repair anything.



That order matters. Run the checks in this sequence and each one narrows what the next has to explain: season-and-year summaries account for the regular wave, the outlier rule finds what is left over, calendar effects clear out the repeating half of those flags, structural breaks catch the one-off level shifts, and missing runs catch what was never recorded in the first place. Skip straight to modelling without this pass and you risk fitting a model to a calendar quirk, a level shift, or a gap it was never built to handle.



Next, you will look at the seasonal shape itself more closely, with plots built specifically to compare one year's pattern against another.