 Chapter 2 Time series graphics 

Chapter 2 Time series graphics

Load common libraries and settings

Load additional libraries

The first thing to do in any data analysis task is to plot the data. Graphs enable many features of the
data to be visualised, including patterns, unusual observations, changes over time, and
relationships between variables. The features that are seen in plots of the data must then be
incorporated, as much as possible, into the forecasting methods to be used. Just as the type of data
determines what forecasting method to use, it also determines what graphs are appropriate. But
before we produce graphs, we need to set up our time series in Python.

2.1 DataFrame objects

A time series consists of a sequence of numerical observations along with information about when
those numbers were recorded. Sometimes a primary observation is accompanied by related data
known as covariates or exogenous variables. This data is commonly stored in a two-dimensional
tabular data structure called a dataframe. In Python, tables are often created as pandas DataFrame
objects, or more recently, as Polars DataFrame objects. If the dataframe takes too much time or
memory to process, distributed dataframe implementations such as Dask, Spark, or Ray can also be
used. We will focus on the pandas DataFrame in this book.

Index and Columns

A pandas DataFrame can be thought of as an array of values arranged along two axes: the index
which identifies rows, and the columns . For example, suppose you have annual observations for
some recent years:
 Year Observation
0 2015 123

1 2016 39

2 2017 78

3 2018 52

4 2019 110

We can store the information as a DataFrame using the [Link]() constructor:

df = [Link]({
"Year": list(range(2015, 2020)),
"Observation": [123, 39, 78, 52, 110],
})

By inspecting the dataframe’s properties, we can see that the index runs from 0 up to (but not
including) 5, and the columns include Year and Observation :

print([Link])
print([Link])

RangeIndex(start=0, stop=5, step=1)

Index(['Year', 'Observation'], dtype='object')

Data types are tracked on a per-column basis:

[Link]

Year int64
Observation int64
dtype: object

Individual columns can be accessed and manipulated as Series objects:

print(type(df["Year"]))
print(df["Year"])

<class '[Link]'>
0 2015
1 2016
2 2017
3 2018
4 2019
Name: Year, dtype: int64
The index can be replaced by one of the columns using the .set_index() method. For example,

year_df = df.set_index("Year")
year_df

Observation
Year
2015 123

2016 39

2017 78

2018 52

2019 110

Timestamps and Periods

Time series observations may be associated with instants in time or spans of time. pandas
supports these concepts through [Link] and [Link] classes, respectively:

print(repr([Link]("2020-01")))
print(repr([Link]("2020-01")))

Timestamp('2020-01-01 00:00:00')
Period('2020-01', 'M')

pandas provides many utilities for converting between text strings and dedicated timestamp or
period representations. Internally, timestamps are stored with nanosecond precision by default, but
for convenience a partial specification is sufficient to refer to the earliest instant matching the
input. For periods, the same logic is used to parse the start of the time span, and the duration is
inferred from the precision provided unless the frequency is indicated via the freq= argument:

print(repr([Link]("2020")))
print(repr([Link]("2020-01-01 12:34")))
print(repr([Link]("2020-01-01")))
print(repr([Link]("2020-01-01", freq="M")))

Timestamp('2020-01-01 00:00:00')
Timestamp('2020-01-01 12:34:00')
Period('2020-01-01', 'D')
Period('2020-01', 'M')
Timestamp sequences can be conveniently constructed using pd.to_datetime() or
pd.date_range() :

ts_few = pd.to_datetime(["2020-01-01", "2020-01-02", "2020-01-03"])

ts_range = pd.date_range("2020-01-01", "2020-01-02", freq="8h")
print(ts_few)
print(ts_range)

DatetimeIndex(['2020-01-01', '2020-01-02',
'2020-01-03'],
dtype='datetime64[ns]', freq=None)
DatetimeIndex(['2020-01-01 00:00:00', '2020-01-01 08:00:00',
'2020-01-01 16:00:00', '2020-01-02 00:00:00'],
dtype='datetime64[ns]', freq='8h')

The .to_period() and .to_timestamp() methods allow conversions between timestamp and
period. The .to_period() method will infer the period duration unless the freq= argument is
specified. For example,

print(ts_range.to_period())
print(ts_range.to_period(freq="D"))
print(ts_range.to_period(freq="W"))
print(ts_range.to_period().to_timestamp())

PeriodIndex(['2020-01-01 00:00',
'2020-01-01 08:00',
'2020-01-01 16:00',
'2020-01-02 00:00'],
dtype='period[8h]')
PeriodIndex(['2020-01-01', '2020-01-01',
'2020-01-01', '2020-01-02'],
dtype='period[D]')
PeriodIndex(['2019-12-30/2020-01-05',
'2019-12-30/2020-01-05',
'2019-12-30/2020-01-05',
'2019-12-30/2020-01-05'],
dtype='period[W-SUN]')
DatetimeIndex(['2020-01-01 00:00:00',
'2020-01-01 08:00:00',
'2020-01-01 16:00:00',
'2020-01-02 00:00:00'],
dtype='datetime64[ns]', freq='8h')

Note that W-SUN indicates weeks starting on Monday and ending on Sunday.

Both timestamps and period start times can be converted to strings with custom formatting using
.strftime() :
 print(ts_few.strftime("%m/%d/%Y"))
print(ts_range.to_period().strftime("%Y %b ~ %H:%M"))

Index(['01/01/2020', '01/02/2020', '01/03/2020'], dtype='object')

Index(['2020 Jan ~ 00:00', '2020 Jan ~ 08:00', '2020 Jan ~ 16:00',
'2020 Jan ~ 00:00'],
dtype='object')

pandas automatically converts between Index and Series types as needed when these sequences
are used as the data or index in a DataFrame . In Series form, date-time related methods and
attributes are nested under the .dt accessor. For example:

df = [Link]({"ts": ts_few})
df = [Link](
period=df["ts"].dt.to_period(),
yr=df["ts"].[Link],
str=df["ts"].[Link]("%A, %B %-d"),
).set_index("ts")
df

period yr str
ts
2020-01-01 2020-01-01 2020 Wednesday, January 1

2020-01-02 2020-01-02 2020 Thursday, January 2

2020-01-03 2020-01-03 2020 Friday, January 3

In this book, we will usually store observation times as timestamps. In some cases we will use
integer types such as the year or number of days from some reference date. Usually timing data will
be stored in a column, often with the name ds as in “date stamp”, but occasionally we will use
set_index() to access certain pandas functions.

Key variables
A DataFrame allows multiple time series to be stored in a single object. Suppose you are interested
in a dataset containing the fastest running times for women’s and men’s track races at the
Olympics, from 100m to 10000m:

olympic_running = pd.read_csv("data/olympic_running_unparsed.csv")
olympic_running.head(8)
 Year Length Sex Time
0 1896 100 men 12.0

1 1900 100 men 11.0

2 1904 100 men 11.0

3 1908 100 men 10.8

4 1912 100 men 10.8

5 1916 100 men NaN

6 1920 100 men 10.8

7 1924 100 men 10.6

This DataFrame contains 312 rows and 4 columns. The data is recorded every four years, and there
are 14 separate time series in this DataFrame . A preview of the first 8 observations is also shown, in
which we can see a missing value occurs in 1916. This is because the Olympics were not held during
World War I.

The 14 time series are uniquely identified by the Length and Sex key variables. The .unique()
method can be used to show the categories of a variable as an array:

print(olympic_running["Sex"].unique())

['men' 'women']

The .drop_duplicates() method can be used to show distinct values in a Series or distinct
combinations of values in multiple Series :

print(olympic_running[["Sex", "Length"]].drop_duplicates())

Sex Length
0 men 100
31 women 100
54 men 200
84 women 200
102 men 400
133 women 400
147 men 800
178 women 800
201 men 1500
232 women 1500
244 men 5000
271 women 5000
277 men 10000
304 women 10000

Working with time series dataframes

We can manipulate DataFrame objects using pandas methods such as .assign() , .rename() ,
.drop() , .groupby() , and .agg() , along with slicing operations such as [] and .loc[] . To

illustrate these, we will use the PBS dataset, containing sales data on pharmaceutical products in
Australia.

pbs = (
pd.read_csv("data/PBS_unparsed.csv", parse_dates=["Month"])
[["Month", "Concession", "Type", "ATC1", "ATC2", "Scripts", "Cost"]]
)
pbs

Month Concession Type ATC1 ATC2 Scripts Cost

0 1991-07-01 Concessional Co-payments A A01 18228 67877

1 1991-08-01 Concessional Co-payments A A01 15327 57011

2 1991-09-01 Concessional Co-payments A A01 14775 55020

3 1991-10-01 Concessional Co-payments A A01 15380 57222

... ... ... ... ... ... ...

67592 2008-03-01 General Safety net Z Z 15 276

67593 2008-04-01 General Safety net Z Z 11 165

67594 2008-05-01 General Safety net Z Z 21 278

67595 2008-06-01 General Safety net Z Z 57 491

This contains monthly data on Medicare Australia prescription data from July 1991 to June 2008.
These are classified according to various concession types, and Anatomical Therapeutic Chemical
(ATC) indexes. For this example, we are interested in the Cost time series (total cost of scripts in
Australian dollars).

We have already used [] to select a subset of the available columns. We can use .loc[] to extract
rows pertaining to A10 scripts.

a10 = [Link][pbs["ATC2"] == "A10"]

a10
 Month Concession Type ATC1 ATC2 Scripts Cost
1524 1991-07-01 Concessional Co-payments A A10 89733 2092878

1525 1991-08-01 Concessional Co-payments A A10 77101 1795733

1526 1991-09-01 Concessional Co-payments A A10 76255 1777231

1527 1991-10-01 Concessional Co-payments A A10 78681 1848507

... ... ... ... ... ... ...

52340 2008-03-01 General Safety net A A10 1119 51773

52341 2008-04-01 General Safety net A A10 721 36289

52342 2008-05-01 General Safety net A A10 1947 101233

52343 2008-06-01 General Safety net A A10 4331 193179

Next we can simplify the resulting object by using .drop() to remove columns not needed in
subsequent analysis.

a10 = (
[Link][pbs["ATC2"] == "A10"]
.drop(columns=["ATC1", "ATC2"])
)
a10

Month Concession Type Scripts Cost

1524 1991-07-01 Concessional Co-payments 89733 2092878

1525 1991-08-01 Concessional Co-payments 77101 1795733

1526 1991-09-01 Concessional Co-payments 76255 1777231

1527 1991-10-01 Concessional Co-payments 78681 1848507

... ... ... ... ...

52340 2008-03-01 General Safety net 1119 51773

52341 2008-04-01 General Safety net 721 36289

52342 2008-05-01 General Safety net 1947 101233

52343 2008-06-01 General Safety net 4331 193179

Another useful method is .agg() which allows us to aggregate data either across all rows or, when
used with .groupby() , within groups based on one or more columns in the DataFrame . For
example, we may wish to compute total cost per month regardless of the Concession or Type keys.

total_cost_df = (
[Link][pbs["ATC2"] == "A10"]
.drop(columns=["ATC1", "ATC2"])
.groupby("Month", as_index=False)
.agg({"Cost": "sum"})
.rename(columns={"Cost": "TotalC"})
)
total_cost_df

Month TotalC
0 1991-07-01 3526591

1 1991-08-01 3180891

2 1991-09-01 3252221

3 1991-10-01 3611003

... ...

200 2008-03-01 18264945

201 2008-04-01 23107677

202 2008-05-01 22912510

203 2008-06-01 19431740

The new column, which we named TotalC using the .rename() method, represents the sum of all
Cost values for each month.

We can create or update columns using the assign() method. Here we remove the rename
operation and change the units of Cost from dollars to millions of dollars:

total_cost_df = (
[Link][pbs["ATC2"] == "A10"]
.drop(columns=["ATC1", "ATC2"])
.groupby("Month", as_index=False)
.agg({"Cost": "sum"})
.assign(Cost=lambda x: (x["Cost"] / 1e6).round(2))
)
total_cost_df

Month Cost
0 1991-07-01 3.5
 Month Cost
1 1991-08-01 3.2

2 1991-09-01 3.2

3 1991-10-01 3.6

... ...

200 2008-03-01 18.3

201 2008-04-01 23.1

202 2008-05-01 22.9

203 2008-06-01 19.4

In this operation, we also demonstrate the use of an unnamed lambda function. These are often
used with .assign() , .loc[] , and elsewhere in method chains to access column data even when
the input dataframe may not be bound to a variable name.

Working with csv files

Real world data is found in diverse storage formats such as spreadsheets, databases, and various
text-based and binary file types. Almost all the data in this book is stored in csv files. The first step
in working with time series data is to read the data into a DataFrame and identify the time column
and any key variables.

For example, suppose we have the following quarterly data stored in a csv file (only the first 10 rows
are shown). This data set provides information on the size of the prison population in Australia,
disaggregated by state, gender, legal status and indigenous status. (Here, ATSI stands for
Aboriginal or Torres Strait Islander.)

Date State Gender Legal Indigenous Count

0 2005-03-01 ACT Female Remanded ATSI 0

1 2005-03-01 ACT Female Remanded Non-ATSI 2

2 2005-03-01 ACT Female Sentenced ATSI 0

3 2005-03-01 ACT Female Sentenced Non-ATSI 5

4 2005-03-01 ACT Male Remanded ATSI 7

5 2005-03-01 ACT Male Remanded Non-ATSI 58

6 2005-03-01 ACT Male Sentenced ATSI 5

 Date State Gender Legal Indigenous Count
7 2005-03-01 ACT Male Sentenced Non-ATSI 101

8 2005-03-01 NSW Female Remanded ATSI 51

9 2005-03-01 NSW Female Remanded Non-ATSI 131

We can read (or load) the file into a DataFrame by calling the pd.read_csv() function. By setting
parse_dates=["Date"] , we ensure the Date column gives a timestamp associated with each
observation. In this case the data is quarterly, and the Date column should be interpreted as the
first date of the last month of the quarter.

prison = (
pd.read_csv("data/prison_population.csv", parse_dates=["Date"])
.rename(columns={"Date": "Quarter"})
.sort_values(by=["State", "Gender", "Legal", "Indigenous"])
)
prison

Quarter State Gender Legal Indigenous Count

0 2005-03-01 ACT Female Remanded ATSI 0

64 2005-06-01 ACT Female Remanded ATSI 1

128 2005-09-01 ACT Female Remanded ATSI 0

192 2005-12-01 ACT Female Remanded ATSI 0

... ... ... ... ... ...

2879 2016-03-01 WA Male Sentenced Non-ATSI 2488

2943 2016-06-01 WA Male Sentenced Non-ATSI 2539

3007 2016-09-01 WA Male Sentenced Non-ATSI 2608

3071 2016-12-01 WA Male Sentenced Non-ATSI 2625

This dataframe contains 64 unique time series corresponding to the unique combinations of the 8
states, 2 genders, 2 legal statuses and 2 indigenous statuses. Each of these series is 48 observations
in length, from 2005 Q1 to 2016 Q4.

For a time series dataframe to be valid for forecasting or other analysis, each distinct time series
must be indicated by a unique combination of identifying key variables; in this case, these are the
State , Gender , Legal and Indigenous columns.
We can also write dataframes into csv files using the .to_csv() method. For example, here we save
the total_cost_df dataframe from the previous section.

total_cost_df.to_csv("data/total_cost_df.csv", index=False)

Make sure to specify the correct path to save the file in your own working directory if you run this
code. If unsure, the current working directory can be checked using [Link]() .

The seasonal period

Seasonalities are particularly important for time series analysis. Some graphics and some models
rely on the repeated seasonal period of the data. The seasonal period is the number of observations
before a seasonal pattern repeats. Some algorithms can detect this automatically using a time
index or column.

Some common periods for different time intervals are shown in the table below:

Minute Hour Day Week Year

Data
Quarters 4

Months 12

Weeks 52

Days 7 365.25

Hours 24 168 8766

Minutes 60 1440 10080 525960

Seconds 60 3600 86400 604800 31557600

For quarterly, monthly and weekly data, there is only one seasonal period — the number of
observations within each year. Let it be noted that actually there are not 52 weeks in a year, but
365.25/7 = 52.18 on average, allowing for a leap year every fourth year. Approximating seasonal
periods to integers can be useful as many seasonal terms in models only support integer seasonal
periods.

If the data is observed more than once per week, then there is often more than one seasonal
pattern in the data. For example, data with daily observations might have weekly (period = 7) or
annual (period = 365.25) seasonal patterns. Similarly, data that are observed every minute might
have hourly (period = 60), daily (period = 24 × 60 = 1440), weekly (period = 24 × 60 × 7 =
10080) and annual seasonality (period = 24 × 60 × 365.25 = 525960).
More complicated (and unusual) seasonal patterns can be analysed using specialised models from
the statsforecast package, such as MSTL , MFLES , and TBATS .

2.2 Time plots

For time series data, the obvious graph to start with is a time plot. That is, the observations are
plotted against the time of observation, with consecutive observations joined by straight lines.
Figure 2.1 shows the weekly economy passenger load on Ansett airlines between Australia’s two
largest cities (Melbourne and Sydney).

melsyd_economy = (
pd.read_csv("data/[Link]", parse_dates=["ds"])
.loc[lambda x: (x["Airports"] == "MEL-SYD")
& (x["Class"] == "Economy")]
.rename(columns={"Airports": "unique_id"})
.assign(y=lambda x: x["y"] / 1000)
)
plot_series(
df=melsyd_economy, id_col="unique_id", time_col="ds", target_col="y",
xlabel="Week [1W]", ylabel="Passengers ('000)",
title="Ansett airlines economy class: Melbourne-Sydney")
 Figure 2.1: Weekly economy passenger load on Ansett Airlines.

For plotting series, we will often use the plot_series() function from the utilsforecast library.
The function needs these parameters:

df : Input dataframe (aka pandas DataFrame)

id_col : Series identifier column (aka unique id)

time_col : Timestamp (aka index)

target_col : Target variable column (aka measurement)

By arranging a dataframe with column names matching the defaults anticipated by utilsforecast
and related packages, we can simplify the function call. For example,

plot_series(melsyd_economy,
xlabel="Week [1W]", ylabel="Passengers ('000)",
title="Ansett airlines economy class: Melbourne-Sydney")

For further details, see the utilsforecast documentation.

This plot reveals some interesting features.

There was a period in 1989 when no passengers were carried — this was due to an industrial
dispute.
 There was a period of reduced load in 1992. This was due to a trial in which some economy
class seats were replaced by business class seats.
A large increase in passenger load occurred in the second half of 1991.
There are some large dips in load around the start of each year. These are due to holiday effects.
There is a long-term fluctuation in the level of the series which increases during 1987,
decreases in 1989, and increases again through 1990 and 1991.

Any model will need to take all these features into account in order to effectively forecast the
passenger load into the future.

A simpler time series is shown in Figure 2.2, using the total_cost_df data saved earlier. In this
case, given that we aggregated the data before, we have just one unique time series, and we don’t
have a unique_id column so we need to create one.

plot_series(total_cost_df.assign(unique_id="total_cost"),
time_col="Month", target_col="Cost",
xlabel="Month [1M]", ylabel="$ (millions)",
title="Australian antidiabetic drug sales")

Figure 2.2: Monthly sales of antidiabetic drugs in Australia.

Here, there is a clear and increasing trend. There is also a strong seasonal pattern that increases in
size as the level of the series increases. The sudden drop at the start of each year is caused by a
government subsidisation scheme that makes it cost-effective for patients to stockpile drugs at the
end of the calendar year. Any forecasts of this series would need to capture the seasonal pattern,
and the fact that the trend is changing slowly.

2.3 Time series patterns

In describing these time series, we often use terms like “trend” and “seasonal” which require
further clarification.

Trend
A trend exists when there is a long-term increase or decrease in the data. It does not have to
be linear. Sometimes we will refer to a trend as “changing direction”, when it might go from
an increasing trend to a decreasing trend. There is a trend in the antidiabetic drug sales data
shown in Figure 2.2.

Seasonal
A seasonal pattern occurs when a time series is affected by seasonal factors such as the time
of the year, the day of the week or the hour of the day. Seasonality is always of a fixed and
known period. The monthly sales of antidiabetic drugs (Figure 2.2) shows seasonality which
is induced partly by the change in the cost of the drugs at the end of the calendar year. (Note
that one series can have more than one seasonal pattern.)

Cyclic
A cycle occurs when the data exhibit rises and falls that are not of a fixed frequency. These
fluctuations are usually due to economic conditions, and are often related to the “business
cycle”. The duration of these fluctuations is usually at least 2 years.

Many people confuse cyclic behaviour with seasonal behaviour, but they are really quite different. If
the fluctuations are not of a fixed frequency then they are cyclic; if the frequency is unchanging
and associated with some aspect of the calendar, then the pattern is seasonal. In general, the
average length of cycles is longer than the length of a seasonal pattern, and the magnitudes of
cycles tend to be more variable than the magnitudes of seasonal patterns.

Many time series include trend, cycles and seasonality. When choosing a forecasting method, we
will first need to identify the time series patterns in the data, and then choose a method that is able
to capture the patterns properly.
The examples in Figure 2.3 show different combinations of these components.

Figure 2.3: Four examples of time series showing different patterns.

1. The monthly housing sales (top left) show strong seasonality within each year, as well as some
strong cyclic behaviour with a period of about 6–10 years. There is no apparent trend in the
data over this period.
2. The US treasury bill contracts (top right) show results from the Chicago market for 100
consecutive trading days in 1981. Here there is no seasonality, but an obvious downward trend.
Possibly, if we had a much longer series, we would see that this downward trend is actually part
of a long cycle, but when viewed over only 100 days it appears to be a trend.
3. The Australian quarterly electricity production (bottom left) shows a strong increasing trend,
with strong seasonality. There is no evidence of any cyclic behaviour here.
4. The daily change in the Google closing stock price (bottom right) has no trend, seasonality or
cyclic behaviour. There are random fluctuations which do not appear to be predictable, and no
strong patterns that would help with developing a forecasting model.

2.4 Seasonal plots

A seasonal plot is similar to a time plot except that the data are plotted against the individual
“seasons” in which the data were observed. An example is given in Figure 2.4 showing the
antidiabetic drug sales.

df = total_cost_df.assign(
Month_name=total_cost_df["Month"].[Link]("%b"),
 Year=total_cost_df["Month"].[Link],
Month_num=total_cost_df["Month"].[Link],
)
unique_years = df["Year"].unique()
year_palette = sns.color_palette("husl", n_colors=len(unique_years))
fig, ax = [Link]()
[Link](data=df, x="Month_num", y="Cost",
hue="Year", palette=year_palette, legend=False, ax=ax)

[Link](
title="Seasonal Plot: Antidiabetic Drug Sales",
xlabel="Month",
ylabel="$ (millions)",
xticks=range(1, 13),
xticklabels=[
"Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec",
]
)
min_year = unique_years.min()
for year, subset in [Link]("Year"):
x = subset["Month_num"].iloc[-1] + .1
y = subset["Cost"].iloc[-1]
color = year_palette[year - min_year]
[Link](x, y, str(year),
ha="left", va="center", fontsize=9, weight="bold", color=color)
[Link]()
 Figure 2.4: Seasonal plot of monthly antidiabetic drug sales in Australia.

This is the same data shown earlier, but now the data from each year is overlapped. A seasonal plot
shows the underlying seasonal pattern more clearly and helps identify years when the pattern
changes.

There is a large jump in sales in January each year. These are probably sales in late December as
customers stockpile before the end of the calendar year, but the sales are not registered with the
government until a week or two later. The graph also shows an unusually small number of sales in
March 2008 (most other years show an increase between February and March). The small number
of sales in June 2008 is probably due to incomplete counting of sales when the data were collected.

Multiple seasonal periods

Where the data has more than one seasonal pattern, we can group the data according to different
seasonalities and plot. The vic_elec data contains half-hourly ( freq="30min" ) electricity demand
for the state of Victoria, Australia. We can plot the daily pattern, weekly pattern or yearly pattern by
grouping according to date, week and year as shown in Figures 2.5–2.7.

vic_elec_df = pd.read_csv("data/vic_elec.csv", parse_dates=["ds"])

vic_elec_demand = vic_elec_df[vic_elec_df["unique_id"] == "Demand"]
df = vic_elec_demand.assign(
hour_minute=lambda x: x["ds"].[Link]("%H:%M"),
 day=lambda x: x["ds"].[Link],
)

fig, ax = [Link]()
[Link](data=df, x="hour_minute", y="y",
hue="day", palette="husl", legend=False, ax=ax)
[Link](
title="Electricity Demand: Victoria",
xlabel="Time",
ylabel="MWh",
)
unique_ticks = df["hour_minute"].unique()
ticks = range(0, len(unique_ticks), 2)
ticklabels = unique_ticks[::2]
ax.set_xticks(ticks, labels=ticklabels, rotation=45)
[Link]()

Figure 2.5: Seasonal plot showing daily seasonal patterns for Victorian electricity demand.

df = (
vic_elec_demand
.loc[lambda x: x["ds"].between("2012-01-02", "2014-12-28 23:59")]
.assign(
day_of_week=lambda x: x["ds"].dt.day_name(),
week=lambda x: x["ds"].dt.to_period("W").dt.start_time,
)
)
unique_weeks = df["week"].unique()
palette = sns.color_palette("husl", n_colors=len(unique_weeks))
color_map = dict(zip(unique_weeks, palette))

fig, ax = [Link]()
for week, df_week in [Link]("week"):
df_week.plot(x="day_of_week", y="y",
ax=ax, color=color_map[week])
ax.get_legend().remove()
[Link](
title="Electricity Demand: Victoria",
xlabel="Time",
ylabel="MWh",
)
[Link]()

Figure 2.6: Seasonal plot showing weekly seasonal patterns for Victorian electricity demand.

df = vic_elec_demand.assign(
day_of_year=lambda x: x["ds"].[Link]("%m-%d"),
year=lambda x: x["ds"].[Link],
)

palette = sns.color_palette("husl", n_colors=df["year"].nunique())

fig, ax = [Link]()
 for i_year, (year, df_year) in enumerate([Link]("year")):
df_year.plot(x="day_of_year", y="y", ax=ax,
label=str(year), color=palette[i_year])
[Link](
title="Electricity Demand: Victoria",
ylabel="MWh",
xlabel="Time",
)
[Link]()

Figure 2.7: Seasonal plot showing yearly seasonal patterns for Victorian electricity demand.

2.5 Seasonal subseries plots

An alternative plot that emphasises the seasonal patterns is where the data for each season are
collected together in separate mini time plots.

df = total_cost_df.assign(
year=total_cost_df["Month"].[Link],
month_name=total_cost_df["Month"].dt.month_name(),
month_idx=total_cost_df["Month"].[Link],
)
 fig, axes = [Link](1, 12, figsize=(9, 3), sharey=True)
for ax, ((_, month_name), month_df) in zip(
axes, [Link](["month_idx", "month_name"])):
mean_cost = month_df["Cost"].mean()
[Link](month_df["year"], month_df["Cost"], color="k")
[Link](mean_cost, color="b", linewidth=1)
[Link](title=month_name[:3], xlabel="")
ax.tick_params(axis="x", rotation=90)

[Link]("Australian Antidiabetic Drug Sales")

[Link]("Month")
[Link]("$(millions)")
[Link]()

Figure 2.8: Seasonal subseries plot of monthly antidiabetic drug sales in Australia.

The blue horizontal lines indicate the means for each month. This form of plot enables the
underlying seasonal pattern to be seen clearly, and also shows the changes in seasonality over time.
It is especially useful in identifying changes within particular seasons. In this example, the plot is
not particularly revealing; but in some cases, this is the most useful way of viewing seasonal
changes over time.

Example: Australian holiday tourism

Australian quarterly vacation data provides an interesting example of how these plots can reveal
information. First we need to extract the relevant data from the tourism dataframe. All the usual
pandas wrangling functions apply. To get the total visitor nights spent on Holiday by State for each
quarter (i.e., ignoring Regions) we can use the following code.

tourism = pd.read_csv("data/[Link]", parse_dates=["ds"])

trips = (
[Link][lambda x: x["Purpose"] == "Holiday"]
.groupby(["State", "ds"], as_index=False)
 .agg({"y": "sum"})
)
trips

State ds y
0 ACT 1998-01-01 196

1 ACT 1998-04-01 127

2 ACT 1998-07-01 111

3 ACT 1998-10-01 170

... ... ...

636 Western Australia 2017-01-01 1134

637 Western Australia 2017-04-01 998

638 Western Australia 2017-07-01 880

639 Western Australia 2017-10-01 1026

Time plots of each series show that there is strong seasonality for most states, but that the seasonal
peaks do not coincide.

fig, ax = [Link]()
palette = sns.color_palette("husl", trips["State"].nunique())
[Link](data=trips, x="ds", y="y", hue="State", palette=palette)
[Link](
title="Australian domestic holidays",
ylabel="Overnight trips ('000)",
xlabel="Quarter [1Q]",
)
[Link](loc="center left", bbox_to_anchor=(1.02, 0.5),
frameon=False, borderaxespad=0, title="State")
[Link]()
 Figure 2.9: Time plots of Australian domestic holidays by state.

To see the timing of the seasonal peaks in each state, we can use a season plot. Figure 2.10 makes it
clear that the southern states of Australia (Tasmania, Victoria and South Australia) have the
strongest tourism in Q1 (their summer), while the northern states (Queensland and the Northern
Territory) have the strongest tourism in Q3 (their dry season).

df = [Link](
Quarter="Q" + trips["ds"].[Link]("string"),
Year=trips["ds"].[Link],
)
num_states, num_years = df[["State", "Year"]].nunique()
palette = sns.color_palette("husl", num_years)
fig, axes = [Link](num_states, sharex=True, figsize=(8, 10))
for ax, (state, state_df) in zip(axes, [Link]("State")):
[Link](data=state_df, x="Quarter", y="y",
hue="Year", palette=palette, ax=ax)
ax.get_legend().remove()
[Link](ylabel="")
[Link](1.02, 0.5, state, va="center", ha="right", rotation=270,
size="medium", transform=[Link])
handles, labels = ax.get_legend_handles_labels()
for h in handles:
[Link](linewidth=4)
[Link](handles, labels, title="Year", loc="center left",
 bbox_to_anchor=(1.05, 0.5), frameon=False, borderaxespad=0)
[Link]("Australian domestic holidays")
[Link]("Overnight trips ('000)")
[Link]()

Figure 2.10: Seasonal plot of Australian domestic holidays by state.

fig, axes = [Link](num_states, 4,

sharex=True, sharey="row", figsize=(8, 11))
for ax, ((state, quarter), sq_df) in zip(
[Link], [Link](["State", "Quarter"])):
[Link](sq_df["Year"], sq_df["y"], color="k")
[Link](sq_df["y"].mean(), color="b", linewidth=1)
ax.tick_params(axis="x", rotation=90)
xticks = sq_df["Year"].loc[lambda x: (x % 5) == 0]
if ax in axes[0]:
[Link](title=quarter, xticks=xticks)
if ax in axes[:, -1]:
[Link](1.02, 0.5, state, va="center", ha="left", rotation=270,
size="medium", transform=[Link])
[Link]("Australian domestic holidays")
[Link]("Quarter")
[Link]("Overnight trips ('000)")
[Link]()
 Figure 2.11: Seasonal subseries plot of Australian domestic holidays by state.

The corresponding subseries plots are shown in Figure 2.11. This figure makes it evident that
Western Australian tourism has jumped markedly in recent years, while Victorian tourism has
increased in Q1 and Q4 but not in the middle of the year.

2.6 Scatterplots

The graphs discussed so far are useful for visualising individual time series. It is also useful to
explore relationships between time series.

Figures 2.12 and 2.13 show two time series: half-hourly electricity demand (in Gigawatts) and
temperature (in degrees Celsius), for 2014 in Victoria, Australia. The temperatures are for
Melbourne, the largest city in Victoria, while the demand values are for the entire state.

plot_series(vic_elec_df, ids=["Demand"],
max_insample_length=2 * 24 * 365,
xlabel="Time [30m]", ylabel="GW",
title="Half-hourly electricity demand: Victoria")
 Figure 2.12: Half-hourly electricity demand in Victoria, Australia, for 2014.

plot_series(vic_elec_df, ids=["Temperature"],
max_insample_length=2 * 24 * 365,
xlabel="Time [30m]", ylabel="Degrees Celsius",
title="Half-hourly temperature: Melbourne, Australia")

Figure 2.13: Half-hourly temperature in Melbourne, Australia, for 2014.

We can study the relationship between demand and temperature by plotting one series against the
other.

elec_2014 = (
vic_elec_df.loc[lambda x: x["ds"].[Link] == 2014]
.pivot(index="ds", columns="unique_id", values="y")
.reset_index()
)
fig, ax = [Link]()
[Link](data=elec_2014, x="Temperature", y="Demand",
ax=ax, linewidth=0, s=10, legend=False)
[Link](
title="Scatter Plot of Demand vs. Temperature",
xlabel="Temperature (degrees Celsius)",
ylabel="Electricity demand (GW)",
 )
[Link]()

Figure 2.14: Half-hourly electricity demand plotted against temperature for 2014 in Victoria, Australia.

This scatterplot helps us to visualise the relationship between the variables. It is clear that high
demand occurs when temperatures are high due to the effect of air-conditioning. But there is also a
heating effect, where demand increases for low temperatures.

Correlation
It is common to compute correlation coefficients to measure the strength of the linear relationship
between two variables. The correlation between variables x and y is given by

∑(xt − xˉ)(yt − yˉ)

r= .
​ ​ ​

ˉ)2 ∑(yt − yˉ)2

∑(xt − x
​ ​ ​ ​ ​

The value of r always lies between −1 and 1 with negative values indicating a negative relationship
and positive values indicating a positive relationship. The graphs in Figure 2.15 show examples of
data sets with varying levels of correlation.
 Figure 2.15: Examples of data sets with different levels of correlation.

The correlation coefficient only measures the strength of the linear relationship between two
variables, and can sometimes be misleading. For example, the correlation for the electricity
demand and temperature data shown in Figure 2.14 is 0.28, but the non-linear relationship is
stronger than that.
 Figure 2.16: Each of these plots has a correlation coefficient of 0.82. Data from Anscombe (1973).

The plots in Figure 2.16 all have correlation coefficients of 0.82, but they have very different
relationships. This shows how important it is to look at the plots of the data and not simply rely on
correlation values.

Scatterplot matrices
When there are several potential predictor variables, it is useful to plot each variable against each
other variable. Consider the eight time series shown in Figure 2.17, showing quarterly visitor
numbers across states and territories of Australia.

visitors = [Link](["State", "ds"], as_index=False)["y"].sum()

num_states = visitors["State"].nunique()
fig, axes = [Link](num_states, 1, sharex=True, figsize=(8, 10))
for ax, (state, state_df) in zip(axes, [Link]("State")):
[Link](state_df["ds"], state_df["y"])
[Link](1.02, 0.5, state, va="center", ha="left",
rotation=270, transform=[Link])
[Link]("Australian domestic tourism")
[Link]("Quarter")
[Link]("Overnight trips ('000)")
[Link]()
 Figure 2.17: Quarterly visitor nights for the states and territories of Australia.

To see the relationships between these eight time series, we can plot each time series against the
others. These plots can be arranged in a scatterplot matrix, as shown in Figure 2.18.

df = [Link](index="ds", columns="State", values="y")

def corrfunc(x, y, **kws):

r, pvalue = pearsonr(x, y)
ax = [Link]()
[Link](
f"Corr: \n{r:.3f}{'***' if pvalue < 0.05 else ''}",
xy=(0.5, 0.5), xycoords="axes fraction",
ha="center", va="center", fontsize=12)

g = [Link](df, height=1.6)
g.map_lower([Link])
g.map_upper(corrfunc)
g.map_diag([Link], lw=2)
[Link](xlabel="")
for i, col in enumerate([Link]):
[Link][0, i].set_title(col, size="medium")
[Link]()
 Figure 2.18: A scatterplot matrix of the quarterly visitor nights in the states and territories of Australia.

For each panel, the variable on the vertical axis is given by the variable name in that row, and the
variable on the horizontal axis is given by the variable name in that column. There are many
options available to produce different plots within each panel. In the default version, the
correlations are shown in the upper right half of the plot, while the scatterplots are shown in the
lower half. On the diagonal are shown density plots.

The value of the scatterplot matrix is that it enables a quick view of the relationships between all
pairs of variables. In this example, mostly positive relationships are revealed, with the strongest
relationships being between the neighbouring states located in the south and south-east coast of
Australia, namely, New South Wales, Victoria and South Australia. Some negative relationships are
also revealed between the Northern Territory and other regions. The Northern Territory is located
in the north of Australia, famous for its outback desert landscapes visited mostly in winter. Hence,
the peak visitation in the Northern Territory is in the July (winter) quarter in contrast to January
(summer) quarter for the rest of the regions.

2.7 Lag plots

Figure 2.19 displays scatterplots of quarterly Australian beer production (introduced in Figure 1.1),
where the horizontal axis shows lagged values of the time series. Each graph shows yt plotted
​

against yt−k for different values of k .

beer = (
pd.read_csv("data/aus_production.csv", parse_dates=["ds"])
[["ds", "Beer"]]
.loc[lambda x: x["ds"] >= "2000"]
.assign(Season=lambda x: "Q" + x["ds"].[Link]("string"))
.rename(columns={"Beer": "y"})
)
lims = .95 * beer["y"].min(), 1.05 * beer["y"].max()
cmap = plt.get_cmap("viridis")
colors = {
season: cmap(i / 3)
for i, season in enumerate(beer["Season"].unique())
}
fig, axes = [Link](3, 3, figsize=(8, 6), sharex=True, sharey=True)
for i, ax in enumerate([Link]):
lag = i + 1
df = [Link](lag=beer["y"].shift(lag))
[Link](data=df, x="lag", y="y", hue="Season",
palette=colors, ax=ax)
[Link](lims, lims, color=".5", ls="--", lw=1, zorder=-1)
[Link](title=f"lag {lag}", xlabel="", ylabel="", aspect="equal")
ax.get_legend().remove()
axes[1, -1].legend(loc="center left", title="Season",
bbox_to_anchor=(1.05, 0.5), frameon=False, borderaxespad=0)
[Link]("lag(Beer, k)")
[Link]("Beer")
[Link]()
 Figure 2.19: Lagged scatterplots for quarterly beer production.

Here the colours indicate the quarter of the variable on the vertical axis. The relationship is
strongly positive at lags 4 and 8, reflecting the strong seasonality in the data. The negative
relationship seen for lags 2 and 6 occurs because peaks (in Q4) are plotted against troughs (in Q2).

2.8 Autocorrelation

Just as correlation measures the extent of a linear relationship between two variables,
autocorrelation measures the linear relationship between lagged values of a time series.

There are several autocorrelation coefficients, corresponding to each panel in the lag plot. For
example, r1 measures the relationship between yt and yt−1 , r2 measures the relationship between
​ ​ ​

yt and yt−2 , and so on.

​ ​
The value of rk can be written as
​

T
∑ (yt − yˉ)(yt−k − yˉ)
​ ​ ​ ​

t=k+1
rk =
​ ​
,
T
∑ (yt − yˉ)2
​ ​ ​

t=1

where T is the length of the time series. The autocorrelation coefficients make up the
autocorrelation function or ACF.

The autocorrelation coefficients for the beer production data can be computed using the acf()
function.

acf = [Link](beer["y"], nlags=9, fft=False, bartlett_confint=False)

acf_df = [Link](acf, name="ACF").to_frame().rename_axis("lag")
acf_df[1:]

lag 1 2 3 4 5 6 7 8 9
ACF -0.053 -0.758 -0.026 0.802 -0.077 -0.657 0.001 0.707 -0.089

The values in the ACF column are r1 , … , r9 , corresponding to the nine scatterplots in Figure 2.19.
​ ​

We usually plot the ACF to see how the correlations change with the lag k . The plot is sometimes
known as a correlogram.

fig, ax = [Link](figsize=(8, 3))

plot_acf(beer["y"], lags=16, ax=ax,
zero=False, bartlett_confint=False, auto_ylims=True)
[Link](
title="Autocorrelation Function for Beer",
xlabel="lag [1Q]", ylabel="acf",
)
[Link]()
 Figure 2.20: Autocorrelation function for quarterly beer production.

In Figure 2.20, we see that:

r4 is higher than for the other lags. This is due to the seasonal pattern in the data: the peaks
​

tend to be four quarters apart and the troughs tend to be four quarters apart.
r2 is more negative than for the other lags because troughs tend to be two quarters behind
​

peaks.
The gray shaded regions indicate whether the correlations are significantly different from zero
(as explained in Section 2.9).

Trend and seasonality in ACF plots

When data have a trend, the autocorrelations for small lags tend to be large and positive because
observations nearby in time are also nearby in value. So the ACF of a trended time series tends to
have positive values that slowly decrease as the lags increase.

When data are seasonal, the autocorrelations will be larger for the seasonal lags (at multiples of the
seasonal period) than for other lags.

When data are both trended and seasonal, you see a combination of these effects. The
total_cost_df data plotted in Figure 2.2 shows both trend and seasonality. Its ACF is shown in

Figure 2.21. The slow decrease in the ACF as the lags increase is due to the trend, while the
“scalloped” shape is due to the seasonality.

fig, ax = [Link](figsize=(8, 3))

plot_acf(total_cost_df["Cost"], lags=48, ax=ax,
zero=False, bartlett_confint=False, auto_ylims=True)
[Link](
title="Autocorrelation Function for Antidiabetic Drug Sales",
xlabel="lag [1M]", ylabel="acf", ylim=(-0.3, None),
 )
[Link]()

Figure 2.21: Autocorrelation function for monthly antidiabetic drug sales in Australia.

2.9 White noise

Time series that show no autocorrelation are called white noise. Figure 2.22 gives an example of a
white noise series.

[Link]["[Link]"] = [7, 3.5]

random = [Link](30)
wn = [Link]({
"y": [Link](0, 1, 50),
"ds": [Link](1, 51),
"unique_id": "wn",
})
plot_series(wn, target_col="y",
xlabel="sample [1]", title="White noise")
 Figure 2.22: White noise series.

fig, ax = [Link](figsize=(8, 3))

plot_acf(wn["y"], lags=16, ax=ax,
zero=False, bartlett_confint=False, auto_ylims=True)
[Link](
title="Autocorrelation Function for White Noise",
xlabel="lag [1]", ylabel="acf",
)
[Link]()

Figure 2.23: Autocorrelation function for the white noise series.

For white noise series, we expect each autocorrelation to be close to zero. Of course, they will not be
exactly equal to zero as there is some random variation. For a white noise series, we expect 95% of
the spikes in the ACF to lie within ±1.96/ T where T is the length of the time series. It is
​

common to plot these bounds on a graph of the ACF (the gray shaded regions above). If one or more
large spikes are outside these bounds, or if substantially more than 5% of spikes are outside these
bounds, then the series is probably not white noise.

In this example, T = 50 and so the bounds are at ±1.96/ 50 = 0.28. All but the first of the
​

autocorrelation coefficients lie within these limits, and it is only just beyond the bounds,
confirming that the data are white noise.

2.10 Exercises

1. Explore the following four time series: Bricks from aus_production , Lynx from pelt ,
GOOG_Close from gafa_stock , Demand from vic_elec .

Use .info() to find out about the data in each series.

What is the time interval of each series?
Use plot_series() to produce a time plot of each series.
For the last plot, modify the axis labels and title.

2. Use .loc[] or query() to find what days corresponded to the peak closing price for each of the
four stocks in gafa_stock .

3. Download the file [Link] from the book website, open it in Excel (or some other
spreadsheet application), and review its contents. You should find four columns of information.
Columns B through D each contain a quarterly series, labelled Sales, AdBudget and GDP. Sales
contains the quarterly sales for a small company over the period 1981–2005. AdBudget is the
advertising budget and GDP is the gross domestic product. All series have been adjusted for
inflation.

a. You can read the data into Python with the following script:

tute1 = pd.read_csv("[Link]", parse_dates=["ds"])

[Link]()

b. Construct time series plots of each of the three series

plot_series(tute1)

4. The us_total.csv contains data on the demand for natural gas in the US.
 a. Download us_total.csv from the book website read in the csv file using pd.read_csv() .
b. Create a dataframe from us_total with year as the index.
c. Plot the annual natural gas consumption by state for the New England area (comprising
the states of Maine, Vermont, New Hampshire, Massachusetts, Connecticut and Rhode
Island).

5. a. Download [Link] from the book website and read in it using pd.read_excel() .
b. Create a dataframe using [Link] .
c. Find what combination of Region and Purpose had the maximum number of overnight
trips on average.
d. Create a new dataframe which combines the Purposes and Regions, and just has total trips
by State.
6. The aus_arrivals data set comprises quarterly international arrivals to Australia from Japan,
New Zealand, UK and the US.

Use plot_series() to visualise the data.

Use seaborn or matplotlib to create seasonal and subseries plots to compare the
differences between the arrivals from these four countries and identify any unusual
observations.

7. Monthly Australian retail data is provided in aus_retail . Select one of the time series as
follows (but choose your own seed value):

random = [Link](12345678)
random_series_id = [Link](aus_retail["Series ID"].unique())
myseries = aus_retail.loc[aus_retail["Series ID"] == random_series_id]

Explore your chosen retail time series using the following functions:

plot_series() , seasonal_decompose() , lag_plot() , plot_acf()

Can you spot any seasonality, cyclicity and trend? What do you learn about the series?

8. Use the following graphics functions: plot_series() , seasonal_decompose() , lag_plot() ,

plot_acf() and explore features from the following time series: “Total Private” Employed from

us_employment , Bricks from aus_production , Hare from pelt , “H02” Cost from PBS , and

Barrels from us_gasoline .

Can you spot any seasonality, cyclicity and trend?

What do you learn about the series?
What can you say about the seasonal patterns?
Can you identify any unusual years?
 9. The following time plots and ACF plots correspond to four different time series. Your task is to
match each time plot in the first row with one of the ACF plots in the second row.

10. The aus_livestock data contains the monthly total number of pigs slaughtered in Victoria,
Australia, from Jul 1972 to Dec 2018. Use .loc[] to extract pig slaughters in Victoria between
1990 and 1995. Use plot_series() and plot_acf() for this data. How do they differ from white
noise? If a longer period of data is used, what difference does it make to the ACF?

11. a. Use the following code to compute the daily changes in Google closing stock prices.

dgoog = (
pd.read_csv("data/gafa_stock.csv", parse_dates=["ds"])
.loc[lambda x: (x["unique_id"] == "GOOG_Close")
& (x["ds"] >= "2018")
]
.assign(
trading_day=lambda x: [Link](1, len(x) + 1),
diff=lambda x: x["y"].diff(),
)
 .set_index("trading_day")
)

b. Why was it necessary to re-index the dataframe?

c. Plot these differences and their ACF.

d. Do the changes in the stock prices look like white noise?

2.11 Further reading

Cleveland (1993) is a classic book on the principles of visualisation for data analysis. While it is
more than 30 years old, the ideas are timeless.
Unwin (2015) is a modern introduction to graphical data analysis using R. It does not have
much information on time series graphics, but plenty of excellent general advice on using
graphics for data analysis.

2.12 Used modules and classes

StatsForecast
StatsForecast class - Core forecasting engine

AutoETS() model - For automatic exponential smoothing

UtilsForecast
plot_series() utility - For creating time series visualisations

Published by OTexts using Quarto.  2026