CUNY DATA 624 Spring 2021 HW1

2.10 Question No. 1

1. Use the help function to explore what the series gold, woolyrnq and gas represent.

ANSWER:

For question one, the author is asking to use exploratory data techniques to get insight and visualize each time series.

Gold

We can confirm that the Gold dataset is a Time Series dataset by applying the class method. The result “ts” confirms that the dataset is a time series.

class(gold)

## [1] "ts"

The function help gives more details about the time series.

help(gold)

## starting httpd help server ... done

Below is the output from the help documentation on the series:

Gold

Daily morning gold prices
Description
Daily morning gold prices in US dollars. 1 January 1985- 31 March 1989.

Usage
gold
Format
Time series data

The “summary” function gives us basic descriptive statistics for the time series including the number of null values (“NA’s”) which is 34.

summary(gold)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   285.0   337.7   403.2   392.5   443.7   593.7      34

The “head” and “tail” functions show us the firt and last six entries as well as the frequency of the time series. In this case, it’s a monthly time series.

head(gold)

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1] 306.25 299.50 303.45 296.75 304.40 298.35

tail(gold)

## Time Series:
## Start = 1103 
## End = 1108 
## Frequency = 1 
## [1]     NA     NA 391.25 383.30 384.00 382.30

The same exploratory analysis is performed on the other time series.

woolyrnq

class(woolyrnq)

## [1] "ts"

help(woolyrnq)

Quarterly production of woollen yarn in Australia
Description
Quarterly production of woollen yarn in Australia: tonnes. Mar 1965 – Sep 1994.

Usage
woolyrnq
Format
Time series data

Source Time Series Data Library. https://pkg.yangzhuoranyang.com/tsdl/

summary(woolyrnq)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3324    4882    5466    5658    6646    7819

head(woolyrnq)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 1965 6172 6709 6633 6660
## 1966 6786 6800

tail(woolyrnq)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 1993      4588 5309 4732
## 1994 4837 6135 6396

We see with this time series that it’s a quarterly time series from 1965 though the 3rd quarter of 1994.

gas

class(gas)

## [1] "ts"

help(gas)

Australian monthly gas production
Description
Australian monthly gas production: 1956–1995.

Usage
gas
Format
Time series data

Source
Australian Bureau of Statistics.

summary(gas)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1646    2675   16788   21415   38629   66600

head(gas)

##       Jan  Feb  Mar  Apr  May  Jun
## 1956 1709 1646 1794 1878 2173 2321

tail(gas)

##        Mar   Apr   May   Jun   Jul   Aug
## 1995 46287 49013 56624 61739 66600 60054

We see that the time series is monthly from 1956 to 1995.

Next, we use autoplot() to plot each of these in separate plots.

GOLD:

autoplot(gold) +
    ggtitle("Daily Morning Gold Prices") +
    xlab("Date Range: 1/1/1985 to 3/31/1989") +
    ylab("Prices in US Dollars $")

WOOLYRNQ:

autoplot(woolyrnq) +
    ggtitle("Quarterly production of woollen yarn in Australia") +
    xlab("Mar 1965 to Sep 1994") +
    ylab("Woollen Yarn in Tonnes.")

GAS:

autoplot(gas) +
    ggtitle(str_to_title("Australian monthly gas production")) +
    xlab(str_to_title("monthly gas production: 1956–1995")) +
    ylab("Units.")

b. What is the frequency of each series? Hint: apply the frequency() function

The frequency of the gold,woolyrnq,and gas time serieses is daily, quarterly, and monthly respectively.

frequency(gold)

## [1] 1

frequency(woolyrnq)

## [1] 4

frequency(gas)

## [1] 12

c. Use which.max() to spot the outlier in the gold series. Which observation was it?

observation = which.max(gold)
max_price <- gold[c(observation)]

print(paste0("The observation at: ",  observation, " shows that the maximum price of gold (outlier) is: ",  max_price  ))

## [1] "The observation at: 770 shows that the maximum price of gold (outlier) is: 593.7"

2.10 Question No. 2

For this question, the author tasks us with loading data from a csv file and converting it into a time series using the function, “ts”.

tute1 <- read.csv("tute1.csv", header=TRUE)
View(tute1)

Here, the frequency argument of the ts function is set to “4” for quarterly analysis.

mytimeseries <- ts(tute1[,-1], start=1981, frequency=4)

When we plot the series with autoplot and set the “facets” argument to “TRUE”, the plot shows different scales on the Y-axis for the Sales, AdBudget, and GDP data.

When the facets argument is set to FALSE, there’s a common y-axis for all three datasets.

autoplot(mytimeseries, facets=TRUE)

autoplot(mytimeseries, facets=FALSE)

2.10 Question No. 3

For this question, the author tasks us to read in an excel file containing Australian retail data.

retaildata <- readxl::read_excel("retail.xlsx", skip=1)

liquor_sales <- ts(retaildata[,"A3349414R"],
  frequency=12, start=c(1982,4))

The autoplot shows monthly liquor sales from April 1982 to December 2013. The plot shows an increasing trend in sales along with a strong seasonality component which grows in strength.

autoplot(liquor_sales) +
    ggtitle(str_to_title("Victoria, Australia  monthly liquor sales")) +
    xlab("April 1982 to December 2013") +
    ylab("Units.")

The ggseasonplot and ggsubseriesplot show the seasonality of liquor sales. The earlier years starting from 1982 are at the bottom of the plot with the later years at the top. The graph shows the yearly increase in sales and towards the last two months of the year, sales tend to spike.

ggseasonplot(liquor_sales)

ggsubseriesplot(liquor_sales)

Setting the polar attribute to TRUE for the ggseasonplot further deomonstrates the seasonal nature of liquor sales with a spike in December sales

ggseasonplot(liquor_sales, year.labels=FALSE, continuous=TRUE, polar = TRUE)

Over the next two plots, the gglagplot and the ggACF plot, we see how observations are plotted against earlier observations or lags against itself. Lag 12 supports the finding that sales spike in the 12 month which is shown in the ggacf plot.

gglagplot(liquor_sales)

ggAcf(liquor_sales)

2.10 Question No. 6

Sales of one-family houses

help('hsales')

Sales of one-family houses
Description
Monthly sales of new one-family houses sold in the USA since 1973.

Usage
hsales
Format
Time series data

Source
Makridakis, Wheelwright and Hyndman (1998) Forecasting: methods and applications, John Wiley & Sons: New York. Chapter 3.

head(hsales)

##      Jan Feb Mar Apr May Jun
## 1973  55  60  68  63  65  61

tail(hsales)

##      Jun Jul Aug Sep Oct Nov
## 1995  64  64  63  55  54  44

recessions = read.table(textConnection(
"Peak, Trough
1973-11-01, 1975-03-01
1980-01-01, 1980-07-01
1981-07-01, 1982-11-01
1990-07-01, 1991-03-01"), sep=',',
colClasses=c('Date', 'Date'), header=TRUE)

recessions

##         Peak     Trough
## 1 1973-11-01 1975-03-01
## 2 1980-01-01 1980-07-01
## 3 1981-07-01 1982-11-01
## 4 1990-07-01 1991-03-01

autoplot(hsales) +
    ggtitle("USA Sales of One-Family Houses") +
    xlab(str_to_title("monthly sales from 1973 to 1996")) +
    ylab("No. of House Sales")

ggseasonplot(hsales)

ggsubseriesplot(hsales)

gglagplot(hsales, lags = 24 )

ggAcf(hsales)

Analysis of Sales of one-family house

The above plots show that there is no trend for the sales of one family houses from 1973 to 1995. However, there are cyclic and seasonal trends. We can see the downward trends for sales matching downturns in the economy. From the seasonal plots, we can clearly see that sales are strongest in the months of March, April, and May. Finally from the lag plots, we can see that the most recent prior period (lag=1) is strongly correlated to the current period. Lags that are further out, esp. lags at 18 to 21, are not as well correlated. This shows that housing sales have recency bias.

Accidential Deaths in USA

help('usdeaths')

Accidental deaths in USA
Description
Monthly accidental deaths in USA.

Usage
usdeaths
Format
Time series data

Source
Makridakis, Wheelwright and Hyndman (1998) Forecasting: methods and applications, John Wiley & Sons: New York. Exercises 2.3 and 2.4.

head(usdeaths)

##        Jan   Feb   Mar   Apr   May   Jun
## 1973  9007  8106  8928  9137 10017 10826

tail(usdeaths)

##        Jul   Aug   Sep   Oct   Nov   Dec
## 1978 10484  9827  9110  9070  8633  9240

frequency(usdeaths)

## [1] 12

autoplot(usdeaths) +
    ggtitle("Accidental deaths in the USA") +
    xlab("Monthly 1973 t0 1978") +
    ylab("No. of Accidental Deaths")

ggseasonplot(usdeaths)

ggsubseriesplot(usdeaths)

gglagplot(usdeaths, lags = 12)

ggAcf(usdeaths)

Analysis of Accidental deaths in USA

The strongest signal in the data is a seasonal trend for accidental deaths. Intuition tells us that since these are the summer months for most of the USA, more people are on the roads driving which is probably the leading cause of accidental deaths. There does not appear to be any cyclical or other trends.

Quarterly clay brick production

help(bricksq)

Quarterly clay brick production
Description
Australian quarterly clay brick production: 1956–1994.

Usage
bricksq
Format
Time series data

Source
Makridakis, Wheelwright and Hyndman (1998) Forecasting: methods and applications, John Wiley & Sons: New York. Chapter 1 and Exercise 2.3.

head(bricksq)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 1956  189  204  208  197
## 1957  187  214

tail(bricksq)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 1993       462  476  443
## 1994  421  472  494

frequency(bricksq)

## [1] 4

autoplot(bricksq) +
    ggtitle(str_to_title("Quarterly clay brick production"))+
    xlab("Quarterly 1956 to 1994") +
    ylab("Units")

ggseasonplot(bricksq)

ggsubseriesplot(bricksq)

gglagplot(bricksq, lags = 21)

ggAcf(bricksq)

Analysis of Quarterly clay brick production

There is a long term growth trend for clay brick production from 1956 to 1994 along with a slight seasonal trend which shows higher production in quarters 2 and 3. The most interesting item is that there is a strong auto correlation with prior periods of time. A 20 period lag is still correlated above .5.

Annual average sunspot area (1875-2015)

help(sunspotarea)

Annual average sunspot area (1875-2015)
Description
Annual averages of the daily sunspot areas (in units of millionths of a hemisphere) for the full sun. Sunspots are magnetic regions that appear as dark spots on the surface of the sun. The Royal Greenwich Observatory compiled daily sunspot observations from May 1874 to 1976. Later data are from the US Air Force and the US National Oceanic and Atmospheric Administration. The data have been calibrated to be consistent across the whole history of observations.

Format
Annual time series of class ts.

Source
NASA

head(sunspotarea)

## Time Series:
## Start = 1875 
## End = 1880 
## Frequency = 1 
## [1] 213.13333 109.28333  92.85833  22.21667  36.33333 446.75000

class(sunspotarea)

## [1] "ts"

tail(sunspotarea)

## Time Series:
## Start = 2010 
## End = 2015 
## Frequency = 1 
## [1]  214.2917  749.5667  796.8833  860.7750 1252.1500  618.8083

frequency(sunspotarea)

## [1] 1

autoplot(sunspotarea) +
    ggtitle(str_to_title("Annual average sunspot area"))+
    xlab("1875 to 2010") +
    ylab("Units")

#ggseasonplot(sunspotarea)

#ggsubseriesplot(sunspotarea)

gglagplot(sunspotarea)

ggAcf(sunspotarea)

Analysis of Annual average sunspot area (1875-2015)

Since the data is collected on an annual basis, there is no seasonal component to it. There is no trend, but there is a cycle in the data. It seems as though in odd numbered years, there is more sunspot activity. There also seems to be negative correlations between lags.

US finished motor gasoline product supplied

help(gasoline)

US finished motor gasoline product supplied.
Description
Weekly data beginning 2 February 1991, ending 20 January 2017. Units are “million barrels per day”.

Format
Time series object of class ts.

Source
US Energy Information Administration.

head(gasoline)

## Time Series:
## Start = 1991.1 
## End = 1991.19582477755 
## Frequency = 52.1785714285714 
## [1] 6.621 6.433 6.582 7.224 6.875 6.947

tail(gasoline)

## Time Series:
## Start = 2016.95352498289 
## End = 2017.04934976044 
## Frequency = 52.1785714285714 
## [1] 9.269 9.278 8.465 8.470 8.069 8.039

frequency(gasoline)

## [1] 52.17857

autoplot(gasoline) +
    ggtitle(str_to_title("US finished motor gasoline product supplied."))+
    xlab("Weekly 1991 to 2017") +
    ylab("Units")

ggseasonplot(gasoline)

#ggsubseriesplot(gasoline)

gglagplot(gasoline)

ggAcf(gasoline, lag.max = 52)

Analysis of US finished motor gasoline product supplied

There is an increasing trend in the US finished motor gasoline product supplied from 1991 to 2017. The seasonality trend shows an uptick in the middle of a year during the summer months and lower in the winter months. Perhaps the most important item is the very strong auto correlation between the lags.