##libraries used
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tsdl)
library(devtools)
## Loading required package: usethis
library(ggplot2)

This week, I used Hyndman’s Time Series Data Library package (tsdl()). I chose to use the data set containing monthly traffic fatalities in Ontario between 1960 and 1974.

data <- tsdl[[10]]

First, I’d like to know if I should be focusing on additive or mutliplicative decomposition. So, I made a few visualizations to get myself started.

autoplot(data, main = "Monthly Traffic Fatalities in Ontario 1960 - 1974", 
         xlab = "Year", 
         ylab = "Traffic Fatalities")

When I look at this, there seems to be seasonality, as well as a positive trend. Let’s take a look at the decompositions to see what this breaks down to be.

##multiplicative
data %>% 
  decompose(type = "multiplicative") %>%
  autoplot() + 
  xlab("Year") + 
  ggtitle("Classical Multiplicative Decomposition of Monthly Traffic Fatalities")

##additive
data %>% 
  decompose(type = "additive") %>%
  autoplot() + 
  xlab("Year") + 
  ggtitle("Classical Additive Decomposition of Monthly Traffic Fatalities")

Based on these two figures, I think that this is more likely to be best decomposed using additive methods, since the remainders seems to be smaller with the additive decomposition and the seasonal trend is consistent throughout the time series. But, I also want to check my assumptions quantitatively.

First, I’ll do the metrics the old way:

##metrics
dec1 <- decompose(data, type = "additive")
dec2 <- decompose(data, type = "multiplicative")

error1<-na.omit(dec1$random)  #generate residuals without NA
error2<-na.omit(dec2$random)  #generate residuals without NA
error1<-as.vector(error1)       #set residuals as a vector, not a time series
error2<-as.vector(error2)

datavector<-length(dec2$random)-length(dec2$random[is.na(dec2$random)])

abserror1<-abs(error1) 

abserror2<-abs(error2)
sqerror1<-error1^2
sqerror2<-error2^2
pererror1<-(abserror1/datavector)*100
pererror2<-(abserror2/datavector)*100
temp1=temp2=temp3=temp4=k=0

#Calculate the denominator for MASE, which compares MAE of model with MAE of Naive or SNaive
for (i in seq(13,length(na.omit(dec2$random)), by=12)){
  temp1=abs(dec1$x[i]-dec1$x[i-12])
  temp2=abs(dec2$x[i]-dec2$x[i-12])
  temp3=temp3+temp1
  temp4=temp4+temp2
  k=k+1
}
denom1=temp3/k
denom2=temp4/k

#Calculate Statisics
me<-c(mean(error1),mean(error2))
mad<-c(mean(abserror1), mean(abserror2))
mse<-c(mean(sqerror1),mean(sqerror2))
rmse<-sqrt(mse)
mape<-c(mean(pererror1),mean(pererror2))
mase=c(mad[1]/denom1, mad[2]/denom2)

rbind(me,mad,mse,rmse,mape, mase)
##              [,1]       [,2]
## me    -0.01190476 1.00167540
## mad   11.96435658 1.00167540
## mse  226.53663401 1.01457619
## rmse  15.05113398 1.00726173
## mape   7.12164082 0.59623536
## mase   0.65078090 0.05448444

So, it’s a good thing that I checked myself, because based off of the calculated statistics, this model is better composed using multiplicative decomposition than additive decomposition, as can be seen by the much smaller root mean square error amounts.

Going forward, however, I will be able to just use the accuracy function for my forecasts.

What is an ETS model?

Now that I know I should be using multiplicative decomposition, I can move on to making my ETS models.

ETS stands for error, trend, and space. We use these three classifications to distinguish between models errors. In other words, this helps us to determine whether the errors are additive or multiplicative. ETS models are state space models, or models that consist of measurement equations that either describe the observed data or unobserved components/states change over time. The possible components for each component are:
Error {A,M}
Trend {N,A,Ad}
Seasonal {N,A,M}

The ets() function does not produce forecasts. It estimates the model parameters and returns information about the fitted model.

Here I set up my test and training data.

#training and testing
str(data)
##  Time-Series [1:180] from 1960 to 1975: 61 65 55 56 91 80 135 129 129 130 ...
##  - attr(*, "source")= chr "Abraham & Ledolter (1983)"
##  - attr(*, "description")= chr "Monthly traffic fatalities in Ontario 1960-1974"
##  - attr(*, "subject")= chr "Transport and tourism"
train <- ts(data[1:144], start = c(1960, 01), frequency = 12)
test <- ts(data[145:180], frequency = 12)

Now I can run and test different ETS models.

#ets model 1
ets1 <- ets(train)
ets1
## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = train) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9779 
## 
##   Initial states:
##     l = 89.4499 
##     b = 1.1857 
##     s = 1.1123 1.1094 1.2716 1.2003 1.2928 1.1775
##            1.0352 0.9433 0.7781 0.7083 0.6709 0.7002
## 
##   sigma:  0.1223
## 
##      AIC     AICc      BIC 
## 1512.266 1517.738 1565.723
#forecast
fc1 <- forecast(ets1,12) # forecasting foreward 12 periods (12 months)
autoplot(fc1, 
         main = "Forecast ETS(M,Ad,M)", 
         ylab = "Traffic Fatalities")

#accuracy
acc1 <- accuracy(fc1, test[1:12])
acc1
##                      ME     RMSE      MAE       MPE      MAPE      MASE
## Training set -0.3595755 14.95825 11.68893 -1.768907  9.592341 0.5494795
## Test set     21.3379423 28.47060 22.11585 12.003967 12.878026 1.0396342
##                   ACF1
## Training set 0.1832511
## Test set            NA
#ets model 2
ets2 <- ets(train, model = "MMM")
ets2
## ETS(M,Md,M) 
## 
## Call:
##  ets(y = train, model = "MMM") 
## 
##   Smoothing parameters:
##     alpha = 2e-04 
##     beta  = 1e-04 
##     gamma = 3e-04 
##     phi   = 0.9749 
## 
##   Initial states:
##     l = 89.6707 
##     b = 1.0117 
##     s = 1.1058 1.1257 1.2613 1.1609 1.3111 1.1929
##            1.0463 0.9482 0.786 0.7038 0.66 0.6979
## 
##   sigma:  0.1225
## 
##      AIC     AICc      BIC 
## 1511.899 1517.371 1565.356
#forecast
fc2 <- forecast(ets2,12)
autoplot(fc2, 
         main = "Forecast ETS(M,Md,M", 
         ylab = "Traffic Fatalities")

#accuracy
acc2 <- accuracy(fc2, test[1:12])
acc2
##                       ME     RMSE      MAE       MPE      MAPE      MASE
## Training set -0.06279213 14.89191 11.59640 -1.452766  9.500212 0.5451298
## Test set     21.54070063 28.89636 22.04472 12.178652 12.744970 1.0362904
##                   ACF1
## Training set 0.1804865
## Test set            NA

As you can see, the two models are actually very close to one another. However, by looking at the plots, neither seems to forecast very well. In fact, in both cases, a naive forecast would actually perform better, considering that the MASE value is greater than one when used with the test set. These are also biased forecasts, as the mean errors are not zero.