This week, I’ve decided to create ETS models using the Monthly corticosteroid drug sales in Australia from July 1991 to June 2008 found in the fpp2 package created by Hyndman. H02 accounts for the total monthly scripts for pharmaceutical products falling under ATC code H02, as recorded by the Australian Health Insurance Commission. Measured in millions of scripts.

library(fpp2)
## Loading required package: ggplot2
## Loading required package: forecast
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Loading required package: fma
## Loading required package: expsmooth
library(ggplot2)
library(forecast)
library(tseries)
str(h02)
##  Time-Series [1:204] from 1992 to 2008: 0.43 0.401 0.432 0.493 0.502 ...
head(h02)
##           Jul      Aug      Sep      Oct      Nov      Dec
## 1991 0.429795 0.400906 0.432159 0.492543 0.502369 0.602652
summary(h02)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3362  0.5889  0.7471  0.7682  0.9376  1.2572
autoplot(h02)

ggseasonplot(h02)

ETS MODELS
Model 1= ETS Automatic Selection 
Model1 <- ets(h02, model ="ZZZ")
Model1 ### M, AD AND M were chosen for MODEL 1
## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = h02, model = "ZZZ") 
## 
##   Smoothing parameters:
##     alpha = 0.1953 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9798 
## 
##   Initial states:
##     l = 0.3945 
##     b = 0.0085 
##     s = 0.874 0.8197 0.7644 0.7693 0.6941 1.2838
##            1.326 1.1765 1.1621 1.0955 1.0422 0.9924
## 
##   sigma:  0.0676
## 
##        AIC       AICc        BIC 
## -122.90601 -119.20871  -63.17985
summary(Model1)
## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = h02, model = "ZZZ") 
## 
##   Smoothing parameters:
##     alpha = 0.1953 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9798 
## 
##   Initial states:
##     l = 0.3945 
##     b = 0.0085 
##     s = 0.874 0.8197 0.7644 0.7693 0.6941 1.2838
##            1.326 1.1765 1.1621 1.0955 1.0422 0.9924
## 
##   sigma:  0.0676
## 
##        AIC       AICc        BIC 
## -122.90601 -119.20871  -63.17985 
## 
## Training set error measures:
##                       ME       RMSE        MAE       MPE     MAPE      MASE
## Training set 0.003872706 0.05096923 0.03903608 0.1125149 5.046201 0.6439654
##                    ACF1
## Training set 0.00612504
##Model1 Forecast 
Model1.forecast <- forecast(Model1,h=8)
Model1.forecast
##          Point Forecast     Lo 80     Hi 80     Lo 95    Hi 95
## Jul 2008      0.9523893 0.8698435 1.0349350 0.8261464 1.078632
## Aug 2008      1.0003428 0.9119936 1.0886921 0.8652243 1.135461
## Sep 2008      1.0516291 0.9570486 1.1462095 0.9069807 1.196277
## Oct 2008      1.1157629 1.0136392 1.2178867 0.9595782 1.271948
## Nov 2008      1.1297979 1.0246211 1.2349748 0.9689438 1.290652
## Dec 2008      1.2734870 1.1529713 1.3940027 1.0891742 1.457800
## Jan 2009      1.2331240 1.1145566 1.3516914 1.0517908 1.414457
## Feb 2009      0.6668077 0.6016956 0.7319197 0.5672273 0.766388
plot(forecast(Model1,h=8))

hist(residuals(Model1))

MODEL 2 : ETS(A,N,N) 
Simple Exponential Smoothing with Additive Errors
Model2<- ets(h02,model = "ANN")
Model2
## ETS(A,N,N) 
## 
## Call:
##  ets(y = h02, model = "ANN") 
## 
##   Smoothing parameters:
##     alpha = 0.8822 
## 
##   Initial states:
##     l = 0.4266 
## 
##   sigma:  0.1553
## 
##      AIC     AICc      BIC 
## 329.0028 329.1228 338.9571
summary(Model2)
## ETS(A,N,N) 
## 
## Call:
##  ets(y = h02, model = "ANN") 
## 
##   Smoothing parameters:
##     alpha = 0.8822 
## 
##   Initial states:
##     l = 0.4266 
## 
##   sigma:  0.1553
## 
##      AIC     AICc      BIC 
## 329.0028 329.1228 338.9571 
## 
## Training set error measures:
##                       ME     RMSE        MAE       MPE     MAPE     MASE
## Training set 0.001899398 0.154527 0.09506631 -2.405843 14.22123 1.568278
##                     ACF1
## Training set 0.007057982
##Model2 Forecasts
Model2.forecast <- forecast(Model2,h=8)
Model2.forecast
##          Point Forecast     Lo 80     Hi 80        Lo 95    Hi 95
## Jul 2008      0.7684112 0.5693990 0.9674234  0.464048244 1.072774
## Aug 2008      0.7684112 0.5030178 1.0338046  0.362526982 1.174295
## Sep 2008      0.7684112 0.4501951 1.0866273  0.281741714 1.255081
## Oct 2008      0.7684112 0.4049703 1.1318521  0.212576352 1.324246
## Nov 2008      0.7684112 0.3647813 1.1720411  0.151112595 1.385710
## Dec 2008      0.7684112 0.3282466 1.2085758  0.095237554 1.441585
## Jan 2009      0.7684112 0.2945202 1.2423022  0.043657484 1.493165
## Feb 2009      0.7684112 0.2630396 1.2737828 -0.004487977 1.541310
plot(forecast(Model2,h=8))

hist(residuals(Model2))

MODEL 3: ETS(M,N,N)
Simple Exponential Smoothing with Multiplicative Errors 
Model3 <- ets(h02, model="MNN")
Model3
## ETS(M,N,N) 
## 
## Call:
##  ets(y = h02, model = "MNN") 
## 
##   Smoothing parameters:
##     alpha = 0.8899 
## 
##   Initial states:
##     l = 0.4148 
## 
##   sigma:  0.1705
## 
##      AIC     AICc      BIC 
## 241.1337 241.2537 251.0881
summary(Model3)
## ETS(M,N,N) 
## 
## Call:
##  ets(y = h02, model = "MNN") 
## 
##   Smoothing parameters:
##     alpha = 0.8899 
## 
##   Initial states:
##     l = 0.4148 
## 
##   sigma:  0.1705
## 
##      AIC     AICc      BIC 
## 241.1337 241.2537 251.0881 
## 
## Training set error measures:
##                       ME      RMSE        MAE       MPE     MAPE     MASE
## Training set 0.001945429 0.1545334 0.09506534 -2.376818 14.22424 1.568262
##                       ACF1
## Training set -0.0001969028
##Model3 Forecast 
Model3.forecast <- forecast(Model3,h=8)
Model3.forecast
##          Point Forecast     Lo 80     Hi 80      Lo 95    Hi 95
## Jul 2008      0.7680207 0.6001766 0.9358647 0.51132535 1.024716
## Aug 2008      0.7680207 0.5418979 0.9941435 0.42219569 1.113846
## Sep 2008      0.7680207 0.4948537 1.0411877 0.35024783 1.185794
## Oct 2008      0.7680207 0.4539350 1.0821064 0.28766804 1.248373
## Nov 2008      0.7680207 0.4169753 1.1190661 0.23114306 1.304898
## Dec 2008      0.7680207 0.3828167 1.1532246 0.17890207 1.357139
## Jan 2009      0.7680207 0.3507552 1.1852862 0.12986819 1.406173
## Feb 2009      0.7680207 0.3203250 1.2157164 0.08332921 1.452712
plot(forecast(Model3,h=8))

hist(residuals(Model3))

par(mfrow=c(3,2))
plot(forecast(Model1,h=8))
hist(residuals(Model1), col = "blue")

plot(forecast(Model2,h=8))
hist(residuals(Model2), col = "green")

plot(forecast(Model3,h=8))
hist(residuals(Model3), col = "red")

Comparing the three models from the graphs alone, it seems Model1 which was the Automatic selection seems
 normally distributed and has the lowest AIC. Followed by Model3 and then Model2.
Double checking that Model1 is indeed the best model by using accuracy()
Model1.accuracy <- accuracy(Model1.forecast)
Model1.accuracy
##                       ME       RMSE        MAE       MPE     MAPE      MASE
## Training set 0.003872706 0.05096923 0.03903608 0.1125149 5.046201 0.6439654
##                    ACF1
## Training set 0.00612504
Model2.accuracy <- accuracy(Model2.forecast)
Model2.accuracy
##                       ME     RMSE        MAE       MPE     MAPE     MASE
## Training set 0.001899398 0.154527 0.09506631 -2.405843 14.22123 1.568278
##                     ACF1
## Training set 0.007057982
Model3.accuracy <- accuracy(Model3.forecast)
Model3.accuracy
##                       ME      RMSE        MAE       MPE     MAPE     MASE
## Training set 0.001945429 0.1545334 0.09506534 -2.376818 14.22424 1.568262
##                       ACF1
## Training set -0.0001969028

As I had stated previously, Model 1 is still the best model used[comparing each Model’s MSAE].