The main objective in this exercise is to explore the following four time series: Bricks from aus_production, Lynx from pelt , Close from gafa_stock, and Demand from vic_elec.
Loading the fpp3 library
I will install the package in my console so that it will not have to install it anytime that i am rendering the file.
library(fpp3)
Warning: package 'fpp3' was built under R version 4.5.3
aus_production is a time series provided by the Australian Bureau of statistics that contains the production indicators of six selected products in Australia such as Beer, Tobacco, Bricks, Cement, electricity and Gas.
The aus_production time interval is Quarterly meaning the data an estimate of the production each quarter.
Case 2: Let’s find out about the data in the pelt time series.
The “gafa_stock” time series is composed of historical stock prices in US dollars of four selected companies including Google, facebook, Apple and Amazon from 2014 to 2018.
Each stock price is recorded daily for the majority, but there are some missing days in the dataset, the time interval for the “gafa_stock” time series will be considered daily irregular.
Case 4 Let’s find out about the data in the vic_elec time series.
The vic_elec time series covers the total electricity demand in Megawatt in Victoria,Australia alongside of the temperature of Melbourne.
The data were recorded every 30 minutes of every single day of 2012. For that reason, the time interval of the Vic_elec series is half-hourly.
Let’s produce a time plot for each series.
To accomplish, we will use the autoplot() function which is very suitable to create time series plot.s
a) Let’s create a time plot for the aus_production series
# Let's create a time plot for the aus_production series with the bricks as variableBricks_plot <- aus_production %>%autoplot(Bricks)+labs(y="number of bricks in Million", title="Quarterly production of Bricks in Australia")Bricks_plot
Warning: Removed 20 rows containing missing values or values outside the scale range
(`geom_line()`).
# Let's create a time plot for the pelt series with Lynx as our variableLynx_plot <- pelt |>autoplot(Lynx) Lynx_plot
# Let's create a time plot for the gafa_stock series with Close as our variableClose_plot <- gafa_stock|>autoplot(Close)Close_plot
# Let's create a time plot for the vic_elec series with Demand as our variableDemand_plot <- vic_elec %>%autoplot(Demand) +labs(y="Electricity Demand in MWh", title="Half-hourly electricity demand for Victoria")Demand_plot
Interpretation
a) Bricks production in the aus_production series
We notice that the production of bricks have been generally increasing throughout the years with the highest peak in the early 1980. We also observe biggest decline in the first quarter of 1983.
b) Lynx furs traded in the pelt series
We observe a seasonal pattern in the the number of canadian Lynx furs traded throughout the years which is understandable because there are most likely more furs traded during cold season than summer.
c) Closing price in the gafa_stock series
We observe a long term significant increase in closing prices of Amazon and Google stocks. Facebook is generally steady with a slightly increase over time and a little decrease by the end of the year 2018.
d) Electricity Demand in the Vic_elec serie
We observe a half-hourly seasonal electricity demands in Victoria,Australia with the highest demand reaching over 8000 megawatt in the early 2014.
Exercise 2.2
Let’s find what days correspond to the peak closing price for each of the four stocks in gafa_stock.
# Let's find what day(s) correspond to the peak closing price for APPLEAAPL_Max_Closing_Price <- gafa_stock |>filter(Symbol=="AAPL") |># we filter only the AAPL stockselect(Symbol,Date,Close) |># we select only the column of interestfilter(Close==max(Close)) # we filter the day with the highest closing priceAAPL_Max_Closing_Price
# A tsibble: 1 x 3 [!]
# Key: Symbol [1]
Symbol Date Close
<chr> <date> <dbl>
1 AAPL 2018-10-03 232.
# Let's find what day(s) correspond to the peak closing price for AMAZONAMZN_Max_Closing_Price <- gafa_stock |>filter(Symbol=="AMZN") |># we filter only the AMZN stockselect(Symbol,Date,Close) |># we select only the column of interestfilter(Close==max(Close)) # we filter the day with the highest closing priceAMZN_Max_Closing_Price
# A tsibble: 1 x 3 [!]
# Key: Symbol [1]
Symbol Date Close
<chr> <date> <dbl>
1 AMZN 2018-09-04 2040.
# Let's find what day(s) correspond to the peak closing price for FACEBOOKFacebook_Max_Closing_Price <- gafa_stock |>filter(Symbol=="FB") |># we filter only the FACEBOOK stockselect(Symbol,Date,Close) |># we select only the column of interestfilter(Close==max(Close)) # we filter the day with the highest closing priceFacebook_Max_Closing_Price
# A tsibble: 1 x 3 [!]
# Key: Symbol [1]
Symbol Date Close
<chr> <date> <dbl>
1 FB 2018-07-25 218.
# Let's find what day(s) correspond to the peak closing price for GOOGLEGoogle_Max_Closing_Price <- gafa_stock |>filter(Symbol=="GOOG") |># we filter only the GOOGLE stockselect(Symbol,Date,Close) |># we select only the column of interestfilter(Close==max(Close)) # we filter the day with the highest closing priceGoogle_Max_Closing_Price
# A tsibble: 1 x 3 [!]
# Key: Symbol [1]
Symbol Date Close
<chr> <date> <dbl>
1 GOOG 2018-07-26 1268.
Exercise 2.3
a) Let’s read the tute1 csv file in R.
In order to read the file,i loaded the file into a public GitHub so that everybody could have access to the file.
The tute1 file is a table of 4 columns labelled Quarter,Sales, AdBudget and GDP where AdBudget is the advertising budget and GDP is the gross domestic product. It represents the quarterly sales for a small company over the period 1981-2005. Since, the sales are recorded quarterly, the index or time interval for the time series that we will create in quarter. For a better visibility, we will convert the month into quarter. For instance,1981-03-01 becomes 1981Quarter1(1981Q1), 1981-06-01 becomes 1981Quarter2(1981Q2),etc…
# A tsibble: 10 x 3 [1Y]
# Key: state [1]
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
c) Let’s Plot the annual natural gas consumption by state for the New England area
us_total_timeseries %>%filter(state=="Maine"| state=="Vermont"| state=="New Hampshire"| state=="Massachusetts"| state=="Connecticut"| state=="Rhode Island") %>%autoplot(y) +labs(y="Quantity of Gas ", title="Annual natural Gas consumption by state for the New England area")
Exercise 2.5:
a) Let’s read the tourism excel file into R
# let's load the readxl library to library(readxl)
Warning: package 'readxl' was built under R version 4.5.3
# Let's read the tourism excel filetourism <-read_excel("C:/Users/herma/OneDrive/Desktop/CUNY SPS GRAD SCHOOL/FALL 2026 CLASSES/DATA 624/HOMEWORKS/Homework 1/tourism.xlsx")head(tourism)
# 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.
tail(tourism)
# A tibble: 6 × 5
Quarter Region State Purpose Trips
<chr> <chr> <chr> <chr> <dbl>
1 2016-07-01 Yorke Peninsula South Australia Visiting 22.8
2 2016-10-01 Yorke Peninsula South Australia Visiting 33.7
3 2017-01-01 Yorke Peninsula South Australia Visiting 46.2
4 2017-04-01 Yorke Peninsula South Australia Visiting 50.6
5 2017-07-01 Yorke Peninsula South Australia Visiting 27.8
6 2017-10-01 Yorke Peninsula South Australia Visiting 46.3
b) Let’s create a tsibble which is identical to the tourism tsibble from the tsibble package.
To reach our goal, we will need to construct a tsibble of five column with a time interval or index equivalent to a Quarter and 3 keys variables including region,State, Purpose and our one measured variable here which is the number trips.
# 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
c) Find what combination of Region and Purpose had the maximum number of overnight trips on average.
# A tsibble: 76 x 4 [1Q]
# Key: Region, Purpose [76]
# Groups: Region [76]
Region Purpose Quarter average_trips
<chr> <chr> <qtr> <dbl>
1 Adelaide Visiting 2017 Q1 270.
2 Adelaide Hills Visiting 2002 Q4 81.1
3 Alice Springs Holiday 1998 Q3 76.5
4 Australia's Coral Coast Holiday 2014 Q3 198.
5 Australia's Golden Outback Business 2017 Q3 174.
6 Australia's North West Business 2016 Q3 297.
7 Australia's South West Holiday 2016 Q1 612.
8 Ballarat Visiting 2004 Q1 103.
9 Barkly Holiday 1998 Q3 37.9
10 Barossa Holiday 2006 Q1 51.0
# ℹ 66 more rows
Create a new tsibble which combines the Purposes and Regions, and just has total trips by State
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.
To knightly solve this exercise, it will be better to do it case by case.
a) Let’s use the autoplot(), gg_season(), gg_subseries(), gg_lag(), ACF() graphics functions to explore the “Total Private”Employed form us_employment time series features.
# # Let's visualize the content of the us_employment time serieshead(us_employment)
# A tsibble: 6 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
tail(us_employment)
# A tsibble: 6 x 4 [1M]
# Key: Series_ID [1]
Month Series_ID Title Employed
<mth> <chr> <chr> <dbl>
1 2019 Apr TEMPHELPN All Employees, Temporary Help Services 2969.
2 2019 May TEMPHELPN All Employees, Temporary Help Services 3010.
3 2019 Jun TEMPHELPN All Employees, Temporary Help Services 3026.
4 2019 Jul TEMPHELPN All Employees, Temporary Help Services 2966.
5 2019 Aug TEMPHELPN All Employees, Temporary Help Services 3040
6 2019 Sep TEMPHELPN All Employees, Temporary Help Services 3109.
# Let's create a new tsibble that contains only Total private employment title.Total_Private_empl <- us_employment %>%select(-Series_ID) %>%# we select only the column of interestfilter(Title=="Total Private")Total_Private_empl
# A tsibble: 969 x 4 [1M]
# Key: Series_ID [1]
Month Title Employed Series_ID
<mth> <chr> <dbl> <chr>
1 1939 Jan Total Private 25338 CEU0500000001
2 1939 Feb Total Private 25447 CEU0500000001
3 1939 Mar Total Private 25833 CEU0500000001
4 1939 Apr Total Private 25801 CEU0500000001
5 1939 May Total Private 26113 CEU0500000001
6 1939 Jun Total Private 26485 CEU0500000001
7 1939 Jul Total Private 26481 CEU0500000001
8 1939 Aug Total Private 26848 CEU0500000001
9 1939 Sep Total Private 27468 CEU0500000001
10 1939 Oct Total Private 27830 CEU0500000001
# ℹ 959 more rows
# Let's use the graphics functions to explore feature from Total_Private_empl tsibble.autoplot(Total_Private_empl)
Plot variable not specified, automatically selected `.vars = Employed`
gg_season(Total_Private_empl)
Plot variable not specified, automatically selected `y = Employed`
gg_subseries(Total_Private_empl)
Plot variable not specified, automatically selected `y = Employed`
gg_lag(Total_Private_empl)
Plot variable not specified, automatically selected `y = Employed`
ACF(Total_Private_empl|>select(Employed))|>autoplot() # The autocorrelation function needs a numerical vector as an entry
Response variable not specified, automatically selected `var = Employed`
Interpretation
Overall,we observe a strong positive trend in the Total private employed from the US employment throughout the years with a sharp decrease around 2008 and 2009 which we believe was the consequence of the great recession. Moreover, the gg_seasonal function shows us a graphic where each year graphs overlap with a clear visible pattern the years. Consequently, we believe that there is a clear annual seasonality.
b) Let’s use the autoplot(), gg_season(), gg_subseries(), gg_lag(), ACF() graphics functions to explore the Bricks from aus_production time series features.
# Let's visualize the content of the aus_production time series.head(aus_production)
# A tsibble: 6 x 7 [1Q]
Quarter Beer Tobacco Bricks Cement Electricity Gas
<qtr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2009 Q1 415 NA NA 1963 58368 196
2 2009 Q2 398 NA NA 2160 57471 238
3 2009 Q3 419 NA NA 2325 58394 252
4 2009 Q4 488 NA NA 2273 57336 210
5 2010 Q1 414 NA NA 1904 58309 205
6 2010 Q2 374 NA NA 2401 58041 236
# Let's create a new tsibble that contains only the number of Bricks produced.Bricks_tsibble <- aus_production %>%select(Quarter,Bricks) # We select only the columns of interestBricks_tsibble
ACF(Bricks_tsibble|>select(Bricks))|>autoplot() # The autocorrelation function needs a numerical vector as an entry
Response variable not specified, automatically selected `var = Bricks`
Interpretation
We observe some similarities in the quarterly production of bricks with a sharp drop around 1983 and 1984. It is important to point out that the third quarter has the highest bricks production in average which we believe could due to a dryer weather that facilitate the productions of bricks. The production of bricks from the aus_production time series can be considered quarterly seasonal.
C) Let’s use the autoplot(), gg_season(), gg_subseries(), gg_lag(), ACF() graphics functions to explore the Hare from pelt time series features.
The pelt time series represents the annual trading records from the Hudson Bay Company for Snowshoe Hare from 1845 to 1935 with no seasonal structure, we will not be able to use both the gg_season() and gg_subseries() because those functions need at least a seasonal observation in the dataset.
# Let's visualize the content of the pelt time series.head(pelt)
# A tsibble: 6 x 3 [1Y]
Year Hare Lynx
<dbl> <dbl> <dbl>
1 1845 19580 30090
2 1846 19600 45150
3 1847 19610 49150
4 1848 11990 39520
5 1849 28040 21230
6 1850 58000 8420
tail(pelt)
# A tsibble: 6 x 3 [1Y]
Year Hare Lynx
<dbl> <dbl> <dbl>
1 1930 4230 6980
2 1931 19520 8310
3 1932 82110 16010
4 1933 89760 24820
5 1934 81660 29700
6 1935 15760 35400
# Let's create a new tsibble that contains only the annual Hare trading record.Hare_tsibble <- pelt %>%select(Year,Hare) # We select only the columns of interestHare_tsibble
# A tsibble: 91 x 2 [1Y]
Year Hare
<dbl> <dbl>
1 1845 19580
2 1846 19600
3 1847 19610
4 1848 11990
5 1849 28040
6 1850 58000
7 1851 74600
8 1852 75090
9 1853 88480
10 1854 61280
# ℹ 81 more rows
# Let's use the graphics functions to explore the annual Hare trading record tsibble feature.autoplot(Hare_tsibble)
Plot variable not specified, automatically selected `.vars = Hare`
gg_lag(Hare_tsibble)
Plot variable not specified, automatically selected `y = Hare`
ACF(Hare_tsibble|>select(Hare))|>autoplot() # The autocorrelation function needs a numerical vector as an entry
Response variable not specified, automatically selected `var = Hare`
Interpretation
There is not a clear visible repeated pattern in the snowshoe Hare trading record throughout the years. Therefore we can conclude that there is no annual seasonal cycle in the snowshoe Hare trading.
d) Let’s use the autoplot(), gg_season(), gg_subseries(), gg_lag(), ACF() graphics functions to explore the “H02” Cost from PBS time series features.
# Let's visualize the content of the PBS time series.head(PBS)
# A tsibble: 6 x 9 [1M]
# Key: Concession, Type, ATC1, ATC2 [1]
Month Concession Type ATC1 ATC1_desc ATC2 ATC2_desc Scripts Cost
<mth> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
1 1991 Jul Concessional Co-paymen… A Alimenta… A01 STOMATOL… 18228 67877
2 1991 Aug Concessional Co-paymen… A Alimenta… A01 STOMATOL… 15327 57011
3 1991 Sep Concessional Co-paymen… A Alimenta… A01 STOMATOL… 14775 55020
4 1991 Oct Concessional Co-paymen… A Alimenta… A01 STOMATOL… 15380 57222
5 1991 Nov Concessional Co-paymen… A Alimenta… A01 STOMATOL… 14371 52120
6 1991 Dec Concessional Co-paymen… A Alimenta… A01 STOMATOL… 15028 54299
tail(PBS)
# A tsibble: 6 x 9 [1M]
# Key: Concession, Type, ATC1, ATC2 [1]
Month Concession Type ATC1 ATC1_desc ATC2 ATC2_desc Scripts Cost
<mth> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
1 2008 Jan General Safety net Z <NA> Z Z 797 9604
2 2008 Feb General Safety net Z <NA> Z Z 135 1591
3 2008 Mar General Safety net Z <NA> Z Z 15 276
4 2008 Apr General Safety net Z <NA> Z Z 11 165
5 2008 May General Safety net Z <NA> Z Z 21 278
6 2008 Jun General Safety net Z <NA> Z Z 57 491
# Let's list the unique values in ATC2 column.unique(PBS$ATC2)
# Let's create a new tsibble that contain only the H02 drug rowsH02_tsibble <- PBS %>%filter(ATC2=="H02") # we pull out only H02 rowsH02_tsibble
# A tsibble: 816 x 9 [1M]
# Key: Concession, Type, ATC1, ATC2 [4]
Month Concession Type ATC1 ATC1_desc ATC2 ATC2_desc Scripts Cost
<mth> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
1 1991 Jul Concessional Co-paym… H Systemic… H02 CORTICOS… 63261 317384
2 1991 Aug Concessional Co-paym… H Systemic… H02 CORTICOS… 53528 269891
3 1991 Sep Concessional Co-paym… H Systemic… H02 CORTICOS… 52822 269703
4 1991 Oct Concessional Co-paym… H Systemic… H02 CORTICOS… 54016 280418
5 1991 Nov Concessional Co-paym… H Systemic… H02 CORTICOS… 49281 268070
6 1991 Dec Concessional Co-paym… H Systemic… H02 CORTICOS… 51798 277139
7 1992 Jan Concessional Co-paym… H Systemic… H02 CORTICOS… 42436 221772
8 1992 Feb Concessional Co-paym… H Systemic… H02 CORTICOS… 52913 272345
9 1992 Mar Concessional Co-paym… H Systemic… H02 CORTICOS… 62908 325700
10 1992 Apr Concessional Co-paym… H Systemic… H02 CORTICOS… 68499 349271
# ℹ 806 more rows
# Let's use the graphics functions to explore the “H02” Cost from PBS time series features.autoplot(H02_tsibble)
Plot variable not specified, automatically selected `.vars = Scripts`
gg_season(H02_tsibble)
Plot variable not specified, automatically selected `y = Scripts`
gg_subseries(H02_tsibble)
Plot variable not specified, automatically selected `y = Scripts`
ACF(H02_tsibble %>%select(Cost))|>autoplot()
Response variable not specified, automatically selected `var = Cost`
Interpretation
The graphics show us an overall repeated pattern of the Corticosteroids (H02)drug prices monthly over the years where each year graphs overlap with a clear underlying pattern visible monthly. We can safely believe that there is a monthly seasonality of the H02 cost throughout the years.
e) Let’s use the autoplot(), gg_season(), gg_subseries(), gg_lag(), ACF() graphics functions to explore the Barrels from us_gasoline. time series features.
# Let's visualize the content of the us_gasoline time series.head(us_gasoline)
# Let's use the graphics functions to explore the Barrels from us_gasoline feature .autoplot(Barrels_tsibble)
Plot variable not specified, automatically selected `.vars = Barrels`
gg_season(Barrels_tsibble)
Plot variable not specified, automatically selected `y = Barrels`
gg_subseries(Barrels_tsibble)
Plot variable not specified, automatically selected `y = Barrels`
gg_lag(Barrels_tsibble)
Plot variable not specified, automatically selected `y = Barrels`
ACF(Barrels_tsibble|>select(Barrels))|>autoplot() # The autocorrelation function needs a numerical vector as an entry
Response variable not specified, automatically selected `var = Barrels`
Interpretation
Looking at our first graph, we can clearly see a steady increases from 1990 to 1999,from 1999 to 2009 and from 2009 upward with a slight difference between 2007 and 2009 which can also be tight to the great recession that strengthen economics conditions. Therefore, we can conclude of an overall confirmation of a weekly seasonality.
