Imports
library(fpp3)## Warning: package 'fpp3' was built under R version 4.1.2## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──## ✓ tibble 3.1.6 ✓ tsibble 1.1.3
## ✓ dplyr 1.0.8 ✓ tsibbledata 0.4.1
## ✓ tidyr 1.2.0 ✓ feasts 0.3.0
## ✓ lubridate 1.8.0 ✓ fable 0.3.2
## ✓ ggplot2 3.3.5 ✓ fabletools 0.3.2## Warning: package 'dplyr' was built under R version 4.1.2## Warning: package 'tidyr' was built under R version 4.1.2## Warning: package 'tsibble' was built under R version 4.1.2## Warning: package 'tsibbledata' was built under R version 4.1.2## Warning: package 'feasts' was built under R version 4.1.2## Warning: package 'fable' was built under R version 4.1.2## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag() masks stats::lag()
## x tsibble::setdiff() masks base::setdiff()
## x tsibble::union() masks base::union()Do the exercises 9.1, 9.2, 9.3, 9.5, 9.6, 9.7, 9.8 in Hyndman. Please submit both the Rpubs link as well as your .rmd file.
Exercise 9.1
Book Image
- Explain the differences among these figures. Do they all indicate that the data are white noise?
The difference between these images has to do with the magnitude of each spike. That being said, they all still look like white noise processes because dspite the magnitude differences, we see that the spikes never cross the 95% confidence interval, or bounds.
- Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they refer to white noise?
The critical values are different in the above images because of how they are calculated. Part of the calculation includes the size of the time series, T, so as the size of a time series increases we tend to see less and less critical values. Another way to think about this is to consider the law of large numbers. As the number of inputs increase, we tend to see a convergence in terms of means of a population. We can think of something similar for outlier detection. As the size of a series increases, we become more and more confident in the placement of our bounds, and less and less surprised by outliers.
Exercise 9.2
head(gafa_stock)amazon <- gafa_stock %>%
filter(Symbol == "AMZN")gg_tsdisplay(amazon, Close, plot_type = "partial")## Warning: Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.
## Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.
From the above, we can see that the time series for Closing Amazon stock price is trending upwards and therefore note stationary. The pACF demonstrates further that the series is not stationary by the placement of the first lag at 1. With this in mind, it would make sense to difference the series once to stabalize.
Exercise 9.3
- Turkish GDP from global_economy
turkey <- global_economy %>%
filter(Country == "Turkey")
# calculate lambda
lambda <- turkey %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
print(lambda)## [1] 0.1572187# unit root test
turkey %>%
features(box_cox(GDP,lambda), unitroot_ndiffs) For the Turkey data, a lambda of .157 and differencing of 1 would produce stationary data.
- Accommodation takings in the state of Tasmania from aus_accommodation
tasmania <- aus_accommodation %>%
filter(State == "Tasmania")
# calculate lambda
lambda <- tasmania %>%
features(Takings, features = guerrero) %>%
pull(lambda_guerrero)
print(lambda)## [1] 0.001819643# unit root test
tasmania %>%
features(box_cox(Takings,lambda), unitroot_ndiffs) For the Tasmania data, a lambda of .0018 and differencing of 1 would produce stationary data.
- Monthly sales from souvenirs
# calculate lambda
lambda <- souvenirs %>%
features(Sales, features = guerrero) %>%
pull(lambda_guerrero)
print(lambda)## [1] 0.002118221# unit root test
souvenirs %>%
features(box_cox(Sales,lambda), unitroot_ndiffs) For the souvenirs data, a lambda of .002 and differencing of 1 would produce stationary data.
Exercise 9.5
For your retail data (from Exercise 8 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
# Get Lambdas for box cox
lambda <- myseries %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
print(lambda)## [1] -0.08866888# KPSS Unit Root Test
myseries %>%
features(box_cox(Turnover, lambda), unitroot_nsdiffs) print(lambda)## [1] -0.08866888Based on the above tests, it would appear that a box cox transformation using lambda = -0.20 and applying differencing of 1 would be appropriate for this time series.
Exercise 9.6
- use R code to generate data:
phi_func <- function(phi) {
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
y[i] <- phi*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
return(sim)
}- Produce a time plot for the series. How does the plot change as you change phi?
Choosing a value of 1 produces a random walk. Choosing a value of -1 produces a strange oscillation around the center.
# plotting with phi = 0.6
phi_func(0.6) %>%
autoplot(y)
# plotting with phi = 0.6
phi_func(1) %>%
autoplot(y)
# plotting with phi = 0.6
phi_func(-1) %>%
autoplot(y)
- Write your own code to generate data from an MA(1) model with θ1=0.6 and σ2=1
ma1_func <- function(theta) {
y <- numeric(100)
e <- rnorm(100)
e[1] <- 0
for(i in 2:100)
y[i] <- theta*e[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
return(sim)
}- Produce a time plot for the series, How does it change when you change theta?
ma1_func(0.6) %>%
autoplot(y)
ma1_func(1) %>%
autoplot(y)
ma1_func(-1) %>%
autoplot(y)
When changing theta, I don’t see much chante in the shape of the graphs. However, I am noticing that the frequency of spikes is reduced as we increase theta.
- generate data from Arima(1,1) model with phi = 0.6, theta = 0.6, o2 = 1
AR1 <- function(phi, theta){
set.seed(42)
y <- numeric(100)
e <- rnorm(100)
e[1] <- 0
for(i in 2:100)
y[i] <- phi*y[i-1] + theta*e[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
return(sim)
}- Generate data from AR(2) nodel with phi = -0.8, theta = 0.3, o2 = 1 (Note that these parameters will give a non-stationary series.)
AR2 <- function(phi1, phi2){
set.seed(42)
y <- numeric(100)
e <- rnorm(100)
for(i in 3:100)
y[i] <- phi1*y[i-1] + phi2*y[i-2] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)
return(sim)
}- graph the two series and compare them
autoplot(AR1(.6, .6), y) +
autolayer(AR2(-.8, .3),y)
Exercise 9.7
Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.
- use ARIMA to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
head(aus_airpassengers)fit <- aus_airpassengers %>%
filter(Year < 2012) %>%
model(ARIMA(Passengers))
report(fit)## Series: Passengers
## Model: ARIMA(0,2,1)
##
## Coefficients:
## ma1
## -0.8756
## s.e. 0.0722
##
## sigma^2 estimated as 4.671: log likelihood=-87.8
## AIC=179.61 AICc=179.93 BIC=182.99fit %>%
forecast(h=10) %>%
autoplot(aus_airpassengers)
fit %>%
gg_tsresiduals()
- write a model in terms of backshift operator
(1−B)^2 * yt = (1−0.87B)et (subscript)
- plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a
fit2 <-aus_airpassengers %>%
filter(Year < 2012) %>%
model(ARIMA(Passengers ~ pdq(0,1,0)))
fit2 %>%
forecast(h=10) %>%
autoplot(aus_airpassengers)
Compared to the part a model, this model seems to under-forecast values.
- plot forecasts from an ARIMA(2,1,2) model with drift and compare these to part a and c
fit3 <-aus_airpassengers %>%
filter(Year < 2012) %>%
model(ARIMA(Passengers ~ pdq(2,1,2)))
fit3 %>%
forecast(h=10) %>%
autoplot(aus_airpassengers)
#removing constant
fit4 <-aus_airpassengers %>%
filter(Year < 2012) %>%
model(ARIMA(Passengers ~ 0 + pdq(2,1,2)))## Warning: 1 error encountered for ARIMA(Passengers ~ 0 + pdq(2, 1, 2))
## [1] non-stationary AR part from CSSThis plot more closelyt resembles the plot from part a, resembling white noise. I was unable to produce a model when removing the constant.
- plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
fit5 <-aus_airpassengers %>%
filter(Year < 2012) %>%
model(ARIMA(Passengers ~ 1 + pdq(0,2,1)))## Warning: Model specification induces a quadratic or higher order polynomial trend.
## This is generally discouraged, consider removing the constant or reducing the number of differences.fit5 %>%
forecast(h=10) %>%
autoplot(aus_airpassengers)
“Warning: Model specification induces a quadratic or higher order polynomial trend. This is generally discouraged, consider removing the constant or reducing the number of differences.”
Exercise 9.8
For the United States GDP series (from global_economy):
- if necessary, find a suitable Box-Cox transformation for the data;
global_economy %>%
filter(Code == "USA") %>%
gg_tsdisplay(GDP, plot_type='partial')
Because we see non-stable trends in data, we should perform boxcox:
lambda <-global_economy %>%
filter(Code == "USA") %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)- fit a suitable ARIMA model to the transformed data using ARIMA();
fit6 <- global_economy %>%
filter(Code == "USA") %>%
model(ARIMA(box_cox(GDP, lambda)))
report(fit6)## Series: GDP
## Model: ARIMA(1,1,0) w/ drift
## Transformation: box_cox(GDP, lambda)
##
## Coefficients:
## ar1 constant
## 0.4586 118.1822
## s.e. 0.1198 9.5047
##
## sigma^2 estimated as 5479: log likelihood=-325.32
## AIC=656.65 AICc=657.1 BIC=662.78- try some other plausible models by experimenting with the orders chosen;
# modeling several
usa_fit <- global_economy %>%
filter(Code == "USA") %>%
model(arima1 = ARIMA(box_cox(GDP,lambda) ~ pdq(1,1,0)),
arima2 = ARIMA(box_cox(GDP,lambda) ~ pdq(1,2,0)),
arima3 = ARIMA(box_cox(GDP,lambda) ~ pdq(2,1,0)),
arima4 = ARIMA(box_cox(GDP,lambda) ~ pdq(2,1,2)),
arima5 = ARIMA(box_cox(GDP,lambda) ~ pdq(1,1,1)))
glance(usa_fit) %>% arrange(AICc) %>% select(.model:BIC)- choose what you think is the best model and check the residual diagnostics;
After choosing a wide range of ARIMA setups, it is clear that model Arima2 (1,2,0) produces the smallest AIC.
usa_fit %>%
select(arima2) %>%
gg_tsresiduals()
Residuals suggest white noise.
- produce forecasts of your fitted model. Do the forecasts look reasonable?
usa_fit %>%
forecast(h=5) %>%
filter(.model=='arima2') %>%
autoplot(global_economy)
Forecase looks plausibleas it matches the trend in previous visualiation
- compare the results with what you would obtain using ETS() (with no transformation).
ets_ <- global_economy %>%
filter(Code == "USA") %>%
model(ETS(GDP))
report(ets_)## Series: GDP
## Model: ETS(M,A,N)
## Smoothing parameters:
## alpha = 0.9990876
## beta = 0.5011949
##
## Initial states:
## l[0] b[0]
## 448093333334 64917355687
##
## sigma^2: 7e-04
##
## AIC AICc BIC
## 3190.787 3191.941 3201.089ets_ %>%
forecast(h=5) %>%
autoplot(global_economy)
The ETS Model report shows significantly higher rates for AIC, BIC, so we can infer that this model does not perform as well.