library(fpp3)   

#_________________________ Particulatar functions I woould need to read data________________
library(readr)   # read_csv()  for exercise 2.3
library(readxl)  # read_excel()  for exercise 2.5
library(USgas)   # us_total     for exercise 2.4

Exercise 2.1

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

a) Use ? (or help()) to find out about the data in each series

# eval=FALSE because  opens the help tab on my IDE and does nothing in a knit.
# Thie is the documentation of the code however I will summarise below this chunk what each one provides
?aus_production
?pelt #Anual is derrived from example dataset with the features being year, Hare and Lynk
?gafa_stock #Daily
?vic_elec

From the help pages:

  • aus_production – Quarterly production of selected commodities in Australia (beer, tobacco, bricks, cement, electricity, gas). Bricks is clay brick production in millions of bricks.
  • pelt – Annual pelt trading records of the Hudson Bay Company, 1845–1935. Lynx is the number of Canadian lynx pelts traded.
  • gafa_stock – Daily historical stock prices for Google, Amazon, Facebook and Apple, 2014–2018, from Yahoo Finance. Close is the closing price (USD) for each trading day.
  • vic_elec – Half-hourly electricity demand for Victoria, Australia, 2012–2014. Demand is total electricity demand in MWh.

b) What is the time interval of each series?

# interval() reads the spacing of the index off each tsibble but it may struggle with the daily because weekends is excluded in this case

interval(aus_production)
## <interval[1]>
## [1] 1Q
interval(pelt)
## <interval[1]>
## [1] 1Y
interval(gafa_stock)
## <interval[1]>
## [1] !
interval(vic_elec)
## <interval[1]>
## [1] 30m

Bricks–> (aus_production) –> Quarterly Lynx–> (pelt)–> Annual Close–> (gafa_stock)–> Daily however weekends are excluded Demand–> (vic_elec)–>Half-hourly

c) Use autoplot() to produce a time plot of each series

# autoplot() sees draw a timeplot and since with timeseries the independent variable is always use the index and y in this case is Bricks
aus_production |> autoplot(Bricks)

Overall I would say an upward trend here, seasonality is in play here. There’s no cyclic evidence here.

#again time is the index and it's the x while Lynx is the y
pelt |> autoplot(Lynx)

Definitely cyclical, you see overtime the up and down consistently

# Since the tsibble key is Symbol, autoplot in this case draws a colored line for each stock and I didn't even have to group_by
gafa_stock |> autoplot(Close)

Overall all of them have upward trends, it’s hard to see if there’s any seasonality here.

d) For the last plot, modify the axis labels and title

# autoplot() returns a ggplot object and ggplot uses + as an operator so that's why were not using |> as it
# passes the object as a first argument and I'm using labs() so I can pass in the requested parameters in one call instead 
# of adding of spamming the + operator and adding them individually.
vic_elec |>
  autoplot(Demand) +
  labs(
    title    = "Half-hourly Electricity Demand: Victoria, Australia",
    x        = "Time (30-minute intervals)",
    y        = "Demand (MWh)"
  )

Exercise 2.2

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

gafa_stock |>
  group_by(Symbol) |>              # group for dplyr operations
  filter(Close == max(Close))    # so max(Close) will give us each group is max and we want to filter our rows by the column Close where it matches the roles with the max close price from each group
## # A tsibble: 4 x 8 [!]
## # Key:       Symbol [4]
## # Groups:    Symbol [4]
##   Symbol Date        Open  High   Low Close Adj_Close   Volume
##   <chr>  <date>     <dbl> <dbl> <dbl> <dbl>     <dbl>    <dbl>
## 1 AAPL   2018-10-03  230.  233.  230.  232.      230. 28654800
## 2 AMZN   2018-09-04 2026. 2050. 2013  2040.     2040.  5721100
## 3 FB     2018-07-25  216.  219.  214.  218.      218. 58954200
## 4 GOOG   2018-07-26 1251  1270. 1249. 1268.     1268.  2405600

Each stock’s closing price peaked in the second half of 2018

Exercise 2.3

Download tute1.csv from the book website, open it in Excel, and review its contents. The file has four columns: Quarter, plus three quarterly series — Sales, AdBudget and GDP — for a small company over 1981–2005, all inflation-adjusted.

a) Read the data into R

# Reading straight from the url so it can be reprodusable and it's not limited to my local machine
tute1 <- read_csv("https://otexts.com/fpp3/extrafiles/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

b) Convert the data to a time series

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

c) Construct time series plots of each of the three series

mytimeseries |>
  pivot_longer(-Quarter) |>
  ggplot(aes(x = Quarter, y = value, colour = name)) +
  geom_line() +
  facet_grid(name ~ ., scales = "free_y") +
  labs(title = "tute1: Sales, AdBudget and GDP", y = NULL)

d) Check what happens when you don’t include facet_grid()

mytimeseries |>
  pivot_longer(-Quarter) |>
  ggplot(aes(x = Quarter, y = value, colour = name)) +
  geom_line() +
  labs(title = "tute1 without facet_grid()", y = NULL)

Without facet_grid() all three series share their y-axis. All there series have different scales with GDP being the largest which is why the ones on the smaller scales gets flattened toward the bottom and their variation is hard to see. facet_grid() with scales = "free_y" gives each series its own axis, hence patterns become more visible.

Exercise 2.4

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

a) Install the USgas package

# eval=FALSE so it  install once in the console and not during every knit.
install.packages("USgas")

b) Create a tsibble from us_total with year as the index and state as the key

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

us_gas
## # A tsibble: 1,266 x 3 [1Y]
## # Key:       state [53]
##     year state        y
##    <int> <chr>    <int>
##  1  1997 Alabama 324158
##  2  1998 Alabama 329134
##  3  1999 Alabama 337270
##  4  2000 Alabama 353614
##  5  2001 Alabama 332693
##  6  2002 Alabama 379343
##  7  2003 Alabama 350345
##  8  2004 Alabama 382367
##  9  2005 Alabama 353156
## 10  2006 Alabama 391093
## # ℹ 1,256 more rows

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).

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

#  us_total is  named y in this dataset .
us_gas |>
  filter(state %in% new_england) |>
  autoplot(y) +
  labs(title = "Annual natural gas consumption in New England states",
       x = "Year",
       y = "Consumption")

Massachusetts and Connecticut consumes the most in that region, while Vermont’s uses very little. Most states show a slight upward trend, with Massachusetts and Connecticut rising quickly after 2005.

Exercise 2.5

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

# read_excel does not read from url unlike read_csv() so it has to be downloaded locally
tourism_xl <- 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.

b) Create a tsibble which is identical to the tourism tsibble from the tsibble package

# Each region, state and purpose combination will be its own time series so all three will be key and also we need to do same yearquarter conversion as we did in 2.3 since they are char an tstibble will look at them as days with ~90 day gaps

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

tourism_ts
## # 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 that it matches the built-in version:

# Make sure that it does matches the tourism tibble
all.equal(tourism_ts, tourism)
## [1] TRUE

c) Find what combination of Region and Purpose had the maximum number of overnight trips on average

tourism_ts |>
  as_tibble() |> #converts into a tibble, so from tstible to tible basically thorwing out the time rule 
  group_by(Region, Purpose) |> # Here it will mark #regions * #purpose  groups --> 76 * 4 = 304 groups
  summarise(avg_trips = mean(Trips), .groups = "drop") |> # calculates the average number of trips for each group and drops the group tag with the ".groups = drop"
  slice_max(avg_trips, n = 5) #takes the top 5 avg_trips, sorts them in descending order and drops every other row 
## # 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.

Visiting friends and relatives in Sydney is the largest region and purpose combination on average

d) Create a new tsibble which combines the Purposes and Regions, and just has total trips by State

# Lets group it by states and add up the trips
tourism_state <- tourism_ts |>
  group_by(State) |>
  summarise(Trips = sum(Trips)) |>
  ungroup()

tourism_state
## # 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
tourism_state |>
  autoplot(Trips) +
  labs(title = "Total overnight trips by State",
       y = " # of Trips ")

Exercise 2.8

Use autoplot(), gg_season(), gg_subseries(), gg_lag() and ACF() to 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.

Each series needs a little preparation first (selecting the right key, dropping missing values, or summing across sub-series), so each section starts with that step and then runs the five plots.

Total Private employment (us_employment)

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

total_private
## # A tsibble: 969 x 4 [1M]
## # Key:       Series_ID [1]
##       Month Series_ID     Title         Employed
##       <mth> <chr>         <chr>            <dbl>
##  1 1939 Jan CEU0500000001 Total Private    25338
##  2 1939 Feb CEU0500000001 Total Private    25447
##  3 1939 Mar CEU0500000001 Total Private    25833
##  4 1939 Apr CEU0500000001 Total Private    25801
##  5 1939 May CEU0500000001 Total Private    26113
##  6 1939 Jun CEU0500000001 Total Private    26485
##  7 1939 Jul CEU0500000001 Total Private    26481
##  8 1939 Aug CEU0500000001 Total Private    26848
##  9 1939 Sep CEU0500000001 Total Private    27468
## 10 1939 Oct CEU0500000001 Total Private    27830
## # ℹ 959 more rows
total_private |> autoplot(Employed) +
  labs(title = "US total private employment", y = "Employed (thousands)")

total_private |> gg_season(Employed, labels = "both") +
  labs(title = "Seasonal plot: total private employment", y = "Employed (thousands)")

total_private |> gg_subseries(Employed) +
  labs(title = "Subseries plot: total private employment", y = "Employed (thousands)")

total_private |> gg_lag(Employed, geom = "point", lags = 1:12) +
  labs(x = "lag(Employed, k)")

total_private |> ACF(Employed, lag_max = 48) |> autoplot() +
  labs(title = "ACF: total private employment")

The total private employment series has a strong upward trend, rising from about 25 million in 1939 to about 130 million in 2019. There is a small seasonal pattern. January and February are the lowest months, and June to August are the highest. However, the difference is only a few percent, so the trend is much larger than the seasonality. You can also see cycles; the dips align with recessions like 1975, 1982, 1991, 2001, and 2008. The ACF slowly declines and has no spikes at 12 or 24, indicating the seasonality is too small to appear. The most unusual part is 2008 to 2009, when employment drops sharply and takes a few years to recover.

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

bricks
## # A tsibble: 198 x 2 [1Q]
##    Quarter Bricks
##      <qtr>  <dbl>
##  1 1956 Q1    189
##  2 1956 Q2    204
##  3 1956 Q3    208
##  4 1956 Q4    197
##  5 1957 Q1    187
##  6 1957 Q2    214
##  7 1957 Q3    227
##  8 1957 Q4    222
##  9 1958 Q1    199
## 10 1958 Q2    229
## # ℹ 188 more rows
bricks |> autoplot(Bricks) +
  labs(title = "Australian clay brick production", y = "Millions of bricks")

bricks |> gg_season(Bricks, labels = "both") +
  labs(title = "Seasonal plot - brick production", y = "Millions of bricks")

bricks |> gg_subseries(Bricks) +
  labs(title = "Subseries plot -  brick production", y = "Millions of bricks")

bricks |> gg_lag(Bricks, geom = "point", lags = 1:8) +
  labs(x = "lag(Bricks, k)")

bricks |> ACF(Bricks, lag_max = 24) |> autoplot() +
  labs(title = "ACF -  brick production")

Overall it trends upward however it comes with some challenges. First it trends upward up until the 1980s then it starts trending downwards. It has a clear quarterly seasonality, Q1 always being the lowest quarter with Q2 and Q3 being the highest.The series is cyclic because over the long run it rises and falls sort of like waves mainly climbing but crashing during the recesions and those crashes dont come at consistent intervals so its more of cycle which also described the unusual years, where we see those big drops, those all tend to line up with Australia is recesions.

Hare pelts (pelt)

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

pelt |> gg_lag(Hare, geom = "point", lags = 1:12) +
  labs(x = "lag(Hare, k)")

pelt |> ACF(Hare, lag_max = 30) |> autoplot() +
  labs(title = "ACF: hare pelts")

# There keep failing on anual data so im showing code but not runing it 
pelt |> gg_season(Hare)
pelt |> gg_subseries(Hare)

There’s no long term trend as it stays consistent, there’s no way of checking for seasonality as we only have one observation per year. Ciclicity is the feature that stands out as we see it sort of oscilates. The peaks you see in 1860 are by far the highest which is the unusual years.

H02 drug cost (PBS)

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

h02
## # A tsibble: 204 x 2 [1M]
##       Month  Cost
##       <mth> <dbl>
##  1 1991 Jul 0.430
##  2 1991 Aug 0.401
##  3 1991 Sep 0.432
##  4 1991 Oct 0.493
##  5 1991 Nov 0.502
##  6 1991 Dec 0.603
##  7 1992 Jan 0.660
##  8 1992 Feb 0.336
##  9 1992 Mar 0.351
## 10 1992 Apr 0.380
## # ℹ 194 more rows
h02 |> autoplot(Cost) +
  labs(title = "H02 (corticosteroids) drug cost, Australia", y = "$ millions")

h02 |> gg_season(Cost, labels = "both") +
  labs(title = "Seasonal plot -  H02 cost", y = "$ millions")

h02 |> gg_subseries(Cost) +
  labs(title = "Subseries plot -  H02 cost", y = "$ millions")

h02 |> gg_lag(Cost, geom = "point", lags = 1:12) +
  labs(x = "lag(Cost, k)")

h02 |> ACF(Cost, lag_max = 48) |> autoplot() +
  labs(title = "ACF -  H02 cost")

there’s clearly an upward trend, there’s some seasonality features here. If you check Feb you see it falls sharply and then it climbs throughout the year in the seasonal plot. I don’t see no cyclic behavior and in terms of unusual, to be honest, I don’t see anything unusual.

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

us_gasoline |> gg_season(Barrels, labels = "right") +
  labs(title = "Seasonal plot -  gasoline supplied", y = "Million barrels per day")

us_gasoline |> gg_subseries(Barrels) +
  labs(title = "Subseries plot -  gasoline supplied", y = "Million barrels per day")

us_gasoline |> gg_lag(Barrels, geom = "point", lags = 1:9) +
  labs(x = "lag(Barrels, k)")

us_gasoline |> ACF(Barrels, lag_max = 104) |> autoplot() +
  labs(title = " gasoline supplied")

Overall there is an upward trend and it had some change in directions at one point but went upward afterwards. There’s clearly some seasonal features which is low on around January then highest around the summer. There’s cyclic features here but there gas a global financial crisis in 2007 and it recovered in 2012 and is reflected here. In 2008 we notice something unusual, that the summer peaks were much lower compared to other peaks