DATA 624 Homework 2

Exercise 3.1

## usnetelec
# load data
data(usnetelec)

# find lambda --> 0.5167714
lambda1 <- BoxCox.lambda(usnetelec)

# plot
autoplot(BoxCox(usnetelec,lambda1)) +
  labs(title="usnetelec Transformed Series") + 
  theme(
    plot.title= element_text(hjust = 0.5, face = 'bold')
  )

## usgdp
# load data
data(usgdp)

# find lambda --> 0.366352
lambda2 <- BoxCox.lambda(usgdp)

# plot
autoplot(BoxCox(usgdp,lambda2)) +
  labs(title="usgdp Transformed Series") +
  theme(
    plot.title= element_text(hjust = 0.5, face = 'bold')
  )

## mcopper
# load data
data(mcopper)

# find lambda --> 0.1919047
lambda3 <- BoxCox.lambda(mcopper)

# plot
autoplot(BoxCox(mcopper,lambda3)) +
  labs(title="mcopper Transformed Series") +
  theme(
    plot.title= element_text(hjust = 0.5, face = 'bold')
  )

## enplanements
# load data
data(enplanements)

# find lambda --> -0.2269461
lambda4 <- BoxCox.lambda(enplanements)

# plot
autoplot(BoxCox(enplanements,lambda4)) +
  labs(title="enplanements Transformed Series") +
  theme(
    plot.title= element_text(hjust = 0.5, face = 'bold')
  )

Exercise 3.2

# get data
data(cangas)

# plot data to see behavior
autoplot(cangas) +
  labs(
    title = "cangas"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

A Box-Cox transformation is not helpful for the cangas data because its variance is already stable and the main feature of the series is a long-term upward trend. This wouldn’t be addressed by a variance-stabilizing transformation like the Box-Cox transformation.

Exercise 3.3

# read excel file
retail <- readxl::read_excel("/Users/kristinlussi/Documents/MSDS/DATA 624/retail.xlsx", skip = 1)

# convert data frame to time series object, select just A3349398A time series
myts <- ts(retail[,"A3349398A"], frequency = 12, start = c(1982, 4))

# find lambda --> lambda is 0.1231563
lambda_retail <- BoxCox.lambda(myts)

# plot
autoplot(BoxCox(myts, lambda_retail)) +
  labs(
    title = "Retail Transformed Series"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

A lambda of 0.1231563 would be used for this Box-Cox transformation.

Exercise 3.4

dole

# get data
data(dole)

# plot the data
autoplot(dole) +
  labs(
    title = "Unemployment Benefits in Australia\n(No Transformation)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

# find lambda 
lambda_dole <- BoxCox.lambda(dole)
lambda_dole

## [1] 0.3290922

Visually, the plot shows a strong upward trend with little visual variance. However, the BoxCox.lambda value is 0.3290922, which suggests that a Box-Cox transformation could stabilize the variance.

autoplot(BoxCox(dole, lambda_dole)) +
  labs(
    title = "Unemployment Benefits in\nAustralia (Box-Cox Transformation)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

Comparing vs. the original plot, we can see that the transformation stabilizes the variance slightly, but there is still a clear upward trend.

usdeaths

# get data
data(usdeaths)

# plot data
autoplot(usdeaths) +
  labs(
    title ="Accidental Deaths in USA (Monthly)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

Visually, the variance in the above plot is already pretty stable. Since the data is monthly accidental deaths in the USA and not every month has the same number of days, it may be better to take a look at the average daily death count in the USA.

dframe <- cbind(Monthly = usdeaths,
                DailyAverage = usdeaths/monthdays(usdeaths))

autoplot(dframe, facet=TRUE) +
  xlab("Years") + ylab("Deaths") +
  labs( 
    title= "Accidental Deaths in USA"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

The bottom plot removes the calendar driven variance so the seasonal pattern represents the genuine changes in accident deaths.

bricksq

# get data
data(bricksq)

# plot data
autoplot(bricksq) +
  labs(
    title = "Quarterly Clay Brick Production"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

Visually, there is an upward trend with possible seasonality. Let’s see if a Box-Cox transformation would stabilize any variance.

# find lambda
lambda_bricksq <- BoxCox.lambda(bricksq)

lambda_bricksq # lambda is 0.2548929 which suggests that a Box-Cox transformation may stabilize the variance.

## [1] 0.2548929

# plot box-cox transformation
autoplot(BoxCox(bricksq, lambda_bricksq)) +
  labs(
    title = "Quarterly Clay Brick Production\n(Box-Cox Transformation)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

The transformation slightly reduces the variance. There is still a clear upward trend and possible seasonal pattern.

Exercise 3.5

# code from textbook
beer <- window(ausbeer, start=1992)
fc <- snaive(beer)
autoplot(fc)

res <- residuals(fc)
autoplot(res)

# check residuals
checkresiduals(res)

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 32.269, df = 8, p-value = 8.336e-05
## 
## Model df: 0.   Total lags used: 8

The p-value is very small which means that there is significant autocorrelation, so we can conclude that the residuals are not white noise. In the histogram of the residuals, we can see that there is a longer right tail, which indicates that the residuals are not normally distributed.

Exercise 3.7

Are the following statements true or false?

a. Good forecast methods should have normally distributed residuals. This is mostly true. It is not necessary for the residuals ot be normally distributed, but it is useful. It helps make the calculation of prediction intervals easier.

b. A model with small residuals will give good forecasts. False. A model with small residuals may indicate overfitting.

c. The best measure of forecast accuracy is MAPE. False. Since this measure is based on percentage error, it runs the risk of being infinite or undefined.

d. If your model doesn’t forecast well, you should make it more complicated. False. Making the model more complicated could risk it overfitting.

e. Always choose the model with the best forecast accuracy as measured on the test set. True. The model with the best test-set accuracy should be chosen.

Exercise 3.8

# a. split retail data into train & testing
myts.train <- window(myts, end=c(2010,12))
myts.test <- window(myts, start=2011)

# b. check to see if data has been split appropriately
autoplot(myts) +
  autolayer(myts.train, series="Training") +
  autolayer(myts.test, series="Test")

# c. calculate forecasts
fc <- snaive(myts.train)

# d. calculate accuracy
accuracy(fc,myts.test)

##                     ME      RMSE       MAE      MPE     MAPE     MASE      ACF1
## Training set  73.94114  88.31208  75.13514 6.068915 6.134838 1.000000 0.6312891
## Test set     115.00000 127.92727 115.00000 4.459712 4.459712 1.530576 0.2653013
##              Theil's U
## Training set        NA
## Test set     0.7267171

# e. check residuals
checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 671.41, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The p-value is extremely small, which indicates that there is autocorrelation and that the residuals are not white noise. The histogram of the residuals is slightly normal, but there is a longer right tail. The ACF plot shows recurring spikes which indicates seasonality. This means that the forecast has not captured all the structure in the data.

How sensitive are the accuracy measures to the training/test split?

I will redo the forecast 3 more times with different train/test splits to assess the accuracy measure sensitivity.

## Forecast 1
# split retail data into train & testing
myts1.train <- window(myts, end=c(2012,12))
myts1.test <- window(myts, start=2013)

# calculate forecasts
fc1 <- snaive(myts1.train)

# calculate accuracy
accuracy(fc1,myts1.test)

##                     ME      RMSE       MAE      MPE     MAPE     MASE
## Training set  73.97955  87.99666  75.21541 5.856189 5.922392 1.000000
## Test set     101.06667 108.76875 101.06667 3.637645 3.637645 1.343696
##                    ACF1 Theil's U
## Training set  0.5956535        NA
## Test set     -0.5562806 0.4962007

## Forecast 2
# split retail data into train & testing
myts2.train <- window(myts, end=c(2009,12))
myts2.test <- window(myts, start=2010)

# calculate forecasts
fc2 <- snaive(myts2.train)

# calculate accuracy
accuracy(fc2,myts2.test)

##                     ME      RMSE      MAE      MPE     MAPE     MASE      ACF1
## Training set  74.29751  88.86355  75.4676 6.197542 6.263306 1.000000 0.6392016
## Test set     104.89583 120.93167 106.6792 4.180821 4.249077 1.413576 0.4949489
##              Theil's U
## Training set        NA
## Test set     0.6631145

## Forecast 3
# split retail data into train & testing
myts3.train <- window(myts, end=c(2005,12))
myts3.test <- window(myts, start=2006)

# calculate forecasts
fc3 <- snaive(myts3.train)

# calculate accuracy
accuracy(fc3,myts3.test)

##                     ME      RMSE       MAE      MPE     MAPE    MASE      ACF1
## Training set  67.23443  79.20459  67.80073 6.396519 6.436720 1.00000 0.5641847
## Test set     209.02083 224.63310 209.02083 9.530367 9.530367 3.08287 0.8273582
##              Theil's U
## Training set        NA
## Test set      1.554264

The accuracy measures for each forecast vary, showing that the accuracy measures are sensitive to the train/test split.

Exercise 3.9

data(visnights)

# Create 3 training sets
train1 <- window(visnights[, "QLDMetro"], end = c(2015, 4)) # omit last 1 year
train2 <- window(visnights[,"QLDMetro"], end = c(2014, 4)) # omit last 2 years
train3 <- window(visnights[, "QLDMetro"], end = c(2013, 4)) # omit last 3 years

# get the corresponding test sets

test1 <- window(visnights[,"QLDMetro"], start = c(2016,1))
test2 <- window(visnights[,"QLDMetro"], start = c(2015,1), end = c(2015, 4))
test3 <- window(visnights[,"QLDMetro"], start = c(2014,1), end = c(2014,4))

# Computer 1 year forecasts for each training set
fc1 <- snaive(train1, h = 4) 
fc2 <- snaive(train2, h =4)
fc3 <- snaive(train3, h=4)

# Compute the accuracy for each forecast
accuracy(fc1,test1)['Test set','MAPE']

## [1] 6.159821

accuracy(fc2,test2)['Test set','MAPE']

## [1] 3.057463

accuracy(fc3,test3)['Test set','MAPE']

## [1] 8.476977

MAPE varies with the training/test split because each test year has different conditions. Training set 2 has the lowest MAPE, which omitted 2 years, suggesting that sometimes excluding very recent data can improve accuracy if those years contained atypical behavior. Training set 3 had 3 years omitted and yielded the highest MAPE. This may indicate that the history does not reflect current conditions. In training set 1, only the most recent year was omitted and yielded the second lowest MAPE. This indicates that the most recent year may not always be the most representative. Overall, this shows that forecast accuracy not only depends on the amount of history used, but also whether the history reflects the conditions of the forecast period.