Objective

Complete and submit exercises 2.1, 2.2, 2.3, 2.4, 2.5 and 2.8 for Homework 1 from the Hyndman online Forecasting book. We have to submit both our Rpubs link as well as attach the .pdf file with our code.

Overview

This homework focuses on time series analytics. A time series is a sequence of historical numerical observations where the past carries information about the present that relationship between observations across time is what distinguishes this data type from others. The goal of the exercises below is to identify and understand that relationship visually.However, before plotting anything, though, there is a prior step: understanding what the data actually measures, what its units are, and how frequently it was observed. Those facts determine which plots are meaningful and what patterns are even possible to see. The exercises follow that order, inspect first, visualize second.

Pre-Requisite

Following two datasets need to be downloaded into the working directory:

  1. tute1.csv for Exercise 2.3 and

  2. tourism.xlsx for Exercise 2.5.

In addition to the datasets four packages are required as follows:

  1. fpp3 is a meta-package that loads the core forecasting functions and most of the datasets used in these exercises in one call.
  2. readxl handles the Excel file in Exercise 2.5.
  3. USgas contains the natural gas consumption data for Exercise 2.4.
  4. readr provides read_csv() for Exercise 2.3.
if (!file.exists("tute1.csv")) {
  download.file("https://bit.ly/fpptute1", "tute1.csv")
}
tute1 <- readr::read_csv("tute1.csv")
head(tute1)
## # A tibble: 6 × 4
##   Quarter    Sales AdBudget   GDP
##   <date>     <dbl>    <dbl> <dbl>
## 1 1981-03-01 1020.     659.  252.
## 2 1981-06-01  889.     589   291.
## 3 1981-09-01  795      512.  291.
## 4 1981-12-01 1004.     614.  292.
## 5 1982-03-01 1058.     647.  279.
## 6 1982-06-01  944.     602   254
if (!file.exists("tourism.xlsx")) {
  download.file("https://bit.ly/fpptourism", "tourism.xlsx", mode = "wb")
}
tourism_xl <- readxl::read_excel("tourism.xlsx")
head(tourism_xl)
## # A tibble: 6 × 5
##   Quarter    Region   State           Purpose  Trips
##   <chr>      <chr>    <chr>           <chr>    <dbl>
## 1 1998-01-01 Adelaide South Australia Business  135.
## 2 1998-04-01 Adelaide South Australia Business  110.
## 3 1998-07-01 Adelaide South Australia Business  166.
## 4 1998-10-01 Adelaide South Australia Business  127.
## 5 1999-01-01 Adelaide South Australia Business  137.
## 6 1999-04-01 Adelaide South Australia Business  200.

Exercise 2.1

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

Approach.

The first thing I want to know about any time series is how often the data was measured. That one fact decides everything that comes later. So my plan is as follows.

  • First I read the helppage for each dataset so I know what the numbers actually mean.
  • Second I use R to find out the time interval instead of guessing by eye, because interval() gives me the real answer.
  • Third and last I plot, so I already know what I am looking at.

These are the help calls the question asks for. I set this chunk to eval = FALSE on purpose. The ? command opens the Help pane inside RStudio, so it shows me the documentation on screen while I work, but it does not return anything that can be printed into a knitted document. Running it here would add nothing to the PDF.

# Help pages for context on each dataset.
# These open in the RStudio Help pane when run interactively.
?aus_production
?pelt
?gafa_stock
?vic_elec

So that the documentation actually appears in this report, I pull the help file in as text and print it. This is the same content the Help pane shows, just captured in a way the document can display.

# Small helper: print the first n lines of a dataset's help page as plain text.
show_help <- function(topic, package, n = 22) {
  db  <- utils::help(topic, package = (package))
  txt <- utils::capture.output(
    tools::Rd2txt(utils:::.getHelpFile(db),
                  options = list(underline_titles = FALSE))
  )
  cat(head(txt, n), sep = "\n")
}

show_help("aus_production", "tsibbledata")
## Quarterly production of selected commodities in Australia.
## 
## Description:
## 
##      Quarterly estimates of selected indicators of manufacturing
##      production in Australia.
## 
## Format:
## 
##      Time series of class 'tsibble'.
## 
## Details:
## 
##      'aus_production' is a half-hourly 'tsibble' with six values:
## 
##        Beer:         Beer production in megalitres.                     
##        Tobacco:      Tobacco and cigarette production in tonnes.        
##        Bricks:       Clay brick production in millions of bricks.       
##        Cement:       Portland cement production in thousands of tonnes. 
##        Electricity:  Electricity production in gigawatt hours.          
##        Gas:          Gas production in petajoules.                      
## 
show_help("pelt", "tsibbledata")
## Pelt trading records
## 
## Description:
## 
##      Hudson Bay Company trading records for Snowshoe Hare and Canadian
##      Lynx furs from 1845 to 1935. This data contains trade records for
##      all areas of the company.
## 
## Format:
## 
##      Time series of class 'tsibble'
## 
## Details:
## 
##      'pelt' is an annual 'tsibble' with two values:
## 
##        Hare:  The number of Snowshoe Hare pelts traded. 
##        Lynx:  The number of Canadian Lynx pelts traded. 
##       
## Source:
## 
##      Hudson Bay Company
show_help("gafa_stock", "tsibbledata")
## GAFA stock prices
## 
## Description:
## 
##      Historical stock prices from 2014-2018 for Google, Amazon,
##      Facebook and Apple. All prices are in $USD.
## 
## Format:
## 
##      Time series of class 'tsibble'
## 
## Details:
## 
##      'gafa_stock' is a 'tsibble' containing data on irregular trading
##      days:
## 
##        Open:       The opening price for the stock.          
##        High:       The stock's highest trading price.        
##        Low:        The stock's lowest trading price.         
##        Close:      The closing price for the stock.          
##        Adj_Close:  The adjusted closing price for the stock. 
##        Volume:     The amount of stock traded.
show_help("vic_elec", "tsibbledata")
## Half-hourly electricity demand for Victoria, Australia
## 
## Description:
## 
##      'vic_elec' is a half-hourly 'tsibble' with three values:
## 
##        Demand:       Total electricity demand in MWh.               
##        Temperature:  Temperature of Melbourne (BOM site 086071).    
##        Holiday:      Indicator for if that day is a public holiday. 
##       
## Format:
## 
##      Time series of class 'tsibble'.
## 
## Details:
## 
##      This data is for operational demand, which is the demand met by
##      local scheduled generating units, semi-scheduled generating units,
##      and non-scheduled intermittent generating units of aggregate
##      capacity larger than 30 MWh, and by generation imports to the
##      region. The operational demand excludes the demand met by
##      non-scheduled non-intermittent generating units, non-scheduled

What the help pages tell me.

  • aus_productionis quarterly manufacturing output for Australia, andBricksis clay brick production in millions of bricks. -peltis the Hudson Bay Company trading records, whereLynx` is the number of Canadian lynx pelts traded.

  • gafa_stock holds daily share prices for Google, Amazon, Facebook and Apple, and Close is the closing price.

  • vic_elec is half-hourly electricity demand for Victoria, where Demand is in megawatt hours.

Knowing the units matters. Without the help page I would not know that Bricks is counted in millions or that Demand is in MWh, and I would not be able to label the axes properly later on.

# Ask R for the time interval instead of guessing from the plot.
interval(aus_production)   # Bricks: quarterly
## <interval[1]>
## [1] 1Q
interval(pelt)             # Lynx: annual
## <interval[1]>
## [1] 1Y
interval(gafa_stock)       # Close: irregular, shown as "!" (trading days only)
## <interval[1]>
## [1] !
interval(vic_elec)         # Demand: half-hourly (30 minutes)
## <interval[1]>
## [1] 30m

What the intervals tell me.`

  • aus_productionis quarterly (1Q`)
  • pelt is annual (1Y)
  • vic_elec is half-hourly (30m)
  • gafa_stockcomes back as!`. The exclamation mark is how tsibble reports an irregular interval
aus_production |> autoplot(Bricks)

pelt |> autoplot(Lynx)

gafa_stock |> autoplot(Close)

vic_elec |> autoplot(Demand)

What I see in the plots.

  • Bricks rises for a long time and then falls off after the 1980s, with a repeating up and down inside each year.
  • Lynx has big waves that repeat roughly every ten years, but the waves are not tied to the calendar, so that is a cycle and not seasonality.
  • Close shows four separate stock lines that mostly drift upward. Demand is so dense that it looks like a solid block, which is normal when you plot years of half-hourly data all at once.

Now the last plot with better labels:

head(vic_elec)
## # A tsibble: 6 x 5 [30m] <Australia/Melbourne>
##   Time                Demand Temperature Date       Holiday
##   <dttm>               <dbl>       <dbl> <date>     <lgl>  
## 1 2012-01-01 00:00:00  4383.        21.4 2012-01-01 TRUE   
## 2 2012-01-01 00:30:00  4263.        21.0 2012-01-01 TRUE   
## 3 2012-01-01 01:00:00  4049.        20.7 2012-01-01 TRUE   
## 4 2012-01-01 01:30:00  3878.        20.6 2012-01-01 TRUE   
## 5 2012-01-01 02:00:00  4036.        20.4 2012-01-01 TRUE   
## 6 2012-01-01 02:30:00  3866.        20.2 2012-01-01 TRUE
vic_elec |>
  autoplot(Demand) +
  labs(
    title = "Half-hourly electricity demand: Victoria, Australia",
    x = "Time",
    y = "Demand (MWh)"
  )

Why I labeled it this way.

The default axis labels just repeat the column names, which means nothing to a reader who has not seen the dataset. Saying “Demand (MWh)” tells them the unit, and the title tells them the place. A plot should be readable on its own without the code next to it.

Exercise 2.2

Use filter() to find what days corresponded to the peak closing price for each of the four stocks in gafa_stock.

My Approach

The key words here are “for each of the four stocks.” If I just filtered on the overall maximum, I would only get one row, and it would be whichever company had the highest price in dollars. That is not what the question asks. So I group by Symbol first. Once the data is grouped, max(Close) is computed inside each company separately, and the filter keeps the top day for each one. Then I ungroup so the result behaves like a normal table again.

gafa_stock |>
  group_by(Symbol) |>
  filter(Close == max(Close)) |>
  ungroup() |>
  select(Symbol, Date, Close) |>
  arrange(desc(Close))
## # A tsibble: 4 x 3 [!]
## # Key:       Symbol [4]
##   Symbol Date       Close
##   <chr>  <date>     <dbl>
## 1 AMZN   2018-09-04 2040.
## 2 GOOG   2018-07-26 1268.
## 3 AAPL   2018-10-03  232.
## 4 FB     2018-07-25  218.

What I found.

  • Each of the four stocks peaks in 2018,
  • Amazon and Google have much higher share prices than Apple and Facebook.
  • That price gap is a good reminder that you cannot compare these four lines directly. A one dollar move means something very different for a 2000 dollar stock than for a 200 dollar stock.

Exercise 2.3

Download the file tute1.csv 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.

Approach

There are three steps here as follows. Reading the CSV gives me a plain table, and R has no idea the Quarter column is time. So I run yearquarter() to turn it into a real quarterly value, then as_tsibble() to tell R that this column is the time index. After that, pivot_longer() stacks the three measures into one long column so ggplot can draw them with a single geom_line() call instead of three.

tute1 <- readr::read_csv("tute1.csv")
head(tute1)
## # A tibble: 6 × 4
##   Quarter    Sales AdBudget   GDP
##   <date>     <dbl>    <dbl> <dbl>
## 1 1981-03-01 1020.     659.  252.
## 2 1981-06-01  889.     589   291.
## 3 1981-09-01  795      512.  291.
## 4 1981-12-01 1004.     614.  292.
## 5 1982-03-01 1058.     647.  279.
## 6 1982-06-01  944.     602   254
mytimeseries <- tute1 |>
  mutate(Quarter = yearquarter(Quarter)) |>
  as_tsibble(index = Quarter)

mytimeseries
## # A tsibble: 100 x 4 [1Q]
##    Quarter Sales AdBudget   GDP
##      <qtr> <dbl>    <dbl> <dbl>
##  1 1981 Q1 1020.     659.  252.
##  2 1981 Q2  889.     589   291.
##  3 1981 Q3  795      512.  291.
##  4 1981 Q4 1004.     614.  292.
##  5 1982 Q1 1058.     647.  279.
##  6 1982 Q2  944.     602   254 
##  7 1982 Q3  778.     531.  296.
##  8 1982 Q4  932.     608.  272.
##  9 1983 Q1  996.     638.  260.
## 10 1983 Q2  908.     582.  280.
## # ℹ 90 more rows
mytimeseries |>
  pivot_longer(-Quarter) |>
  ggplot(aes(x = Quarter, y = value, colour = name)) +
  geom_line() +
  facet_grid(name ~ ., scales = "free_y")

# The same plot without facet_grid()
mytimeseries |>
  pivot_longer(-Quarter) |>
  ggplot(aes(x = Quarter, y = value, colour = name)) +
  geom_line()

What changes when I drop facet_grid().

With facets, each series gets its own panel and its own y axis because of scales = "free_y". That lets me see the shape of each series clearly. Without facets, all three lines are squeezed onto one shared y axis. GDP sits on a much smaller scale than Sales and AdBudget, so its line gets flattened and its movement is hard to see.

The lesson I take from this is that a shared axis is only fair when the series are on similar scales. When they are not, separate panels are the honest way to show them.

Exercise 2.4

The USgas package contains data on the demand for natural gas in the US.

My thought process.**

In this data the index is year because that is the time column. The key is state because that is what tells one series apart from another. Without setting the key, R would see several rows sharing the same year and would complain that the rows are not unique. After that, filtering down to the six New England states it is straightforward.

us_total_ts <- us_total |>
  as_tsibble(index = year, key = state)

new_england <- c("Maine", "Vermont", "New Hampshire",
                 "Massachusetts", "Connecticut", "Rhode Island")

us_total_ts |>
  filter(state %in% new_england) |>
  autoplot(y) +
  labs(
    title = "Annual natural gas consumption: New England states",
    x = "Year",
    y = "Consumption (million cubic feet)"
  )

What I see

Massachusetts uses by far the most gas, with Connecticut second, which makes sense because they have much bigger populations than the rest. Vermont uses the least. Most of the states trend upward over the period.

Maine is the one that really stands out. It sits near 6,000 through the late 1990s and then jumps to about 96,000 by 2001, which is roughly fifteen times higher in two years. A change that large in that short a time cannot be homes suddenly using more gas. It almost certainly reflects new gas fired power plants coming online, which is a change in how electricity was generated rather than a change in household demand.

Exercise 2.5

Download tourism.xlsx from the book website and read it into R using readxl::read_excel().

My thought process

The goal is to match the built in tourism object exactly, so I first ask what makes each row unique. A row is one region, one state, one purpose, in one quarter. That means the key needs all three of Region, State, and Purpose together. The Quarter column comes in from Excel as text, so I have to convert it with yearquarter() before it can be used as an index.

For the average part, the words “on average” are what matter. I want the mean across all quarters for each region and purpose pair, not the single biggest quarter.

tourism_xl <- readxl::read_excel("tourism.xlsx")

my_tourism <- tourism_xl |>
  mutate(Quarter = yearquarter(Quarter)) |>
  as_tsibble(index = Quarter, key = c(Region, State, Purpose))

my_tourism
## # A tsibble: 24,320 x 5 [1Q]
## # Key:       Region, State, Purpose [304]
##    Quarter Region   State           Purpose  Trips
##      <qtr> <chr>    <chr>           <chr>    <dbl>
##  1 1998 Q1 Adelaide South Australia Business  135.
##  2 1998 Q2 Adelaide South Australia Business  110.
##  3 1998 Q3 Adelaide South Australia Business  166.
##  4 1998 Q4 Adelaide South Australia Business  127.
##  5 1999 Q1 Adelaide South Australia Business  137.
##  6 1999 Q2 Adelaide South Australia Business  200.
##  7 1999 Q3 Adelaide South Australia Business  169.
##  8 1999 Q4 Adelaide South Australia Business  134.
##  9 2000 Q1 Adelaide South Australia Business  154.
## 10 2000 Q2 Adelaide South Australia Business  169.
## # ℹ 24,310 more rows
# Check it matches the built in tourism tsibble.
identical(dim(my_tourism), dim(tourism))
## [1] TRUE
all.equal(as_tibble(my_tourism), as_tibble(tourism))
## [1] TRUE
# Region and Purpose combination with the highest average overnight trips.
my_tourism |>
  as_tibble() |>
  group_by(Region, Purpose) |>
  summarise(avg_trips = mean(Trips), .groups = "drop") |>
  slice_max(avg_trips, n = 5)
## # A tibble: 5 × 3
##   Region          Purpose  avg_trips
##   <chr>           <chr>        <dbl>
## 1 Sydney          Visiting      747.
## 2 Melbourne       Visiting      619.
## 3 Sydney          Business      602.
## 4 North Coast NSW Holiday       588.
## 5 Sydney          Holiday       550.

What I found

Sydney with the purpose “Visiting” has the highest average number of overnight trips. That fits what I would expect. Sydney is the largest city in Australia, so it has the most friends and relatives for people to travel and visit.

# Total trips by State, with Region and Purpose combined together.
state_trips <- my_tourism |>
  group_by(State) |>
  summarise(Trips = sum(Trips))

state_trips
## # A tsibble: 640 x 3 [1Q]
## # Key:       State [8]
##    State Quarter Trips
##    <chr>   <qtr> <dbl>
##  1 ACT   1998 Q1  551.
##  2 ACT   1998 Q2  416.
##  3 ACT   1998 Q3  436.
##  4 ACT   1998 Q4  450.
##  5 ACT   1999 Q1  379.
##  6 ACT   1999 Q2  558.
##  7 ACT   1999 Q3  449.
##  8 ACT   1999 Q4  595.
##  9 ACT   2000 Q1  600.
## 10 ACT   2000 Q2  557.
## # ℹ 630 more rows
state_trips |>
  autoplot(Trips) +
  labs(
    title = "Total overnight trips by State",
    x = "Quarter",
    y = "Trips (thousands)"
  )

Why this works

Because state_trips is still a tsibble, group_by() and summarise() add up the trips within each quarter instead of collapsing time away. The result keeps Quarter as the index and uses State as the only key, which is exactly what the question asked for.

Exercise 2.8

Use the following graphics functions: autoplot(), gg_season(), gg_subseries(), gg_lag(), 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.

My thought process

Five series is a lot, so instead of running all five plots on all five series I picked the plots that actually say something for each one. Annual data like Hare cannot have a seasonal plot at all, because there is only one observation per year. Recognizing that up front saved me from writing code that would just throw an error.

A note on which plots I skipped and why. gg_season() and gg_subseries() need more than one observation per year, so neither applies to Hare, which is annual. gg_subseries() on us_gasoline would draw roughly 52 panels, one per week of the year, which is technically valid and completely unreadable, so I used the seasonal plot alone for that one. I ran gg_lag() only on Hare below, because it is the series where the lag structure says something the ACF does not already say more clearly.

Total Private Employed (us_employment)

us_private <- us_employment |>
  filter(Title == "Total Private")

us_private |> autoplot(Employed) +
  labs(title = "US Total Private Employment", y = "Employed (thousands)")

us_private |> gg_season(Employed) +
  labs(title = "Seasonal plot: US Total Private Employment")

us_private |> ACF(Employed) |> autoplot() +
  labs(title = "ACF: US Total Private Employment")

Strong upward trend across the whole period. There is mild seasonality with a dip each January after holiday hiring ends. The clearly unusual years are the recessions, and the drop around 2008 and 2009 is the deepest one in the series. The ACF decays very slowly, which is what a strongly trended series always looks like.

Bricks (aus_production)

bricks <- aus_production |>
  filter(!is.na(Bricks))

bricks |> autoplot(Bricks) + labs(title = "Australian brick production")

bricks |> gg_season(Bricks) + labs(title = "Seasonal plot: Bricks")

bricks |> gg_subseries(Bricks) + labs(title = "Subseries plot: Bricks")

bricks |> ACF(Bricks) |> autoplot() + labs(title = "ACF: Bricks")

Bricks rise until about 1980 and then fall, so the trend changes direction partway through. That is a good example of why you cannot describe a whole series with one word. Seasonally, Q1 is the weakest quarter and Q3 is the strongest, which fits the building season. The unusual years are the sharp drops in the early 1980s and again in the early 1990s, both recession periods for Australian construction. I filtered out the missing values first, since Bricks is not recorded for the most recent years in the dataset.

Hare (pelt)

pelt |> autoplot(Hare) + labs(title = "Snowshoe hare pelts traded")

pelt |> ACF(Hare) |> autoplot() + labs(title = "ACF: Hare")

This series is annual, so seasonal and subseries plots do not apply. There is no real trend. What dominates is a strong cycle of roughly nine to ten years, which is the well known predator and prey cycle with the lynx. The ACF shows this clearly by rising and falling in a wave and going negative, instead of decaying steadily the way a trended series does. This is the cleanest example in the whole homework of a cycle that is not seasonality, because the peaks do not land on a fixed calendar period.

H02 Cost (PBS)

h02 <- PBS |>
  filter(ATC2 == "H02") |>
  summarise(Cost = sum(Cost))

h02 |> autoplot(Cost) + labs(title = "H02 drug cost", y = "Cost ($)")

h02 |> gg_season(Cost) + labs(title = "Seasonal plot: H02 Cost")

h02 |> gg_subseries(Cost) + labs(title = "Subseries plot: H02 Cost")

h02 |> ACF(Cost) |> autoplot() + labs(title = "ACF: H02 Cost")

I had to sum across the Concession and Type categories first, because PBS splits H02 into several series and the question asks about H02 overall. There is an upward trend and a very strong yearly pattern. Cost is lowest in February and climbs through the year. That pattern comes from the way the Australian subsidy scheme resets at the start of each calendar year, so it is a rule made by people rather than a weather effect. Seasonality does not have to come from nature.

Barrels (us_gasoline)

us_gasoline |> autoplot(Barrels) +
  labs(title = "US finished motor gasoline product supplied",
       y = "Million barrels per day")

us_gasoline |> gg_season(Barrels) + labs(title = "Seasonal plot: US gasoline")

us_gasoline |> ACF(Barrels) |> autoplot() + labs(title = "ACF: US gasoline")

This is weekly data, so the seasonal period is about 52 rather than 12, and the plots look noisier than the monthly ones. The trend rises until around 2007, flattens, dips during the 2008 recession, and then recovers. Within each year, demand is lowest in winter and highest in summer, which matches the American summer driving season. The unusual years are 2008 and 2009, when the recession and high fuel prices pushed demand down.