The “forecastHybrid” package provides functions to build composite models using multiple individual component models from the “forecast” package. These hybridModel objects can then be manipulated with many of the familiar functions from the “forecast” and “stats” packages including forecast(), plot(), accuracy(), residuals(), and fitted().

Installation

The stable release of the package is hosted on CRAN and can be installed as usual.

install.packages("forecastHybrid")

The latest development version can be installed using the “devtools” package.

devtools::install_github("ellisp/forecastHybrid/pkg")

Version updates to CRAN will be published frequently after new features are implemented, so the development version is not recommended unless you plan to modify the code.

Basic usage

First load the package.

library(forecastHybrid)

Quick start

If you don’t have time to read the whole guide and want to get started immediatly with sane default settings to forecast the USAccDeaths timeseries, run the following:

quickModel <- hybridModel(USAccDeaths)
## Fitting the auto.arima model
## Fitting the ets model
## Fitting the thetam model
## Fitting the nnetar model
## Fitting the stlm model
## Fitting the tbats model
forecast(quickModel)
##          Point Forecast    Lo 80     Hi 80    Lo 95     Hi 95
## Jan 1979       8349.316 7924.712  8962.828 7689.475  9234.317
## Feb 1979       7530.656 6822.403  8087.708 6495.787  8400.442
## Mar 1979       8246.344 7247.559  8886.679 6972.738  9146.115
## Apr 1979       8531.061 7587.667  9194.629 7319.658  9500.477
## May 1979       9331.277 8152.325 10112.349 7861.785 10442.376
## Jun 1979       9772.018 8530.086 10525.745 8167.672 10878.297
## Jul 1979      10680.651 9173.315 11613.448 8793.448 11987.171
## Aug 1979       9995.591 9040.758 10830.331 8794.328 11224.086
## Sep 1979       8998.457 8281.314  9944.791 7892.801 10357.609
## Oct 1979       9246.975 8249.186 10198.510 7994.536 10629.548
## Nov 1979       8774.772 8023.765  9732.246 7589.183 10180.765
## Dec 1979       9095.710 8255.912 10255.628 7816.759 10720.971
## Jan 1980       8365.816 7485.300  9496.297 7033.420 10011.749
## Feb 1980       7576.833 6687.306  8749.398 6141.503  9295.201
## Mar 1980       8266.173 7112.862  9586.522 6813.434 10161.075
## Apr 1980       8546.281 7394.330  9940.466 7064.387 10542.396
## May 1980       9333.430 7817.080 10861.976 7408.566 11490.092
## Jun 1980       9755.489 8138.164 11280.306 7766.472 11933.559
## Jul 1980      10648.213 8608.503 12373.743 8119.529 13051.201
## Aug 1980       9966.667 8647.591 11596.927 8248.308 12297.755
## Sep 1980       8986.324 7949.455 10718.099 7261.423 11441.542
## Oct 1980       9239.293 7942.495 10978.824 7417.369 11724.196
## Nov 1980       8773.647 7623.215 10519.779 6856.540 11286.453
## Dec 1980       9092.796 7951.715 11050.531 7288.262 11837.932
plot(forecast(quickModel), main = "Forecast from auto.arima, ets, thetam, nnetar, stlm, and tbats model")

Fitting a model

The workhorse function of the package is hybridModel(), a function that combines several component models from the “forecast” package. At a minimum, the user must supply a ts or numeric vector for y. In this case, the ensemble will include all six component models: auto.arima(), ets(), thetam(), nnetar(), stlm(), and tbats(). To instead use only a subset of these models, pass a character string to the models argument with the first letter of each model to include. For example, to build an ensemble model on a simulated dataset with auto.arima(), ets(), and tbats() components, run

# Build a hybrid forecast on a simulated dataset using auto.arima, ets, and tbats models.
# Each model is given equal weight 
set.seed(12345)
series <- ts(rnorm(18), f = 2)
hm1 <- hybridModel(y = series, models = "aet", weights = "equal")
## Fitting the auto.arima model
## Fitting the ets model
## Fitting the tbats model

The individual component models are stored inside the hybridModel objects and can viewed in their respective slots, and all the regular methods from the “forecast” package could be applied to these individual component models.

# View the individual models 
hm1$auto.arima
## Series: structure(c(0.585528817843856, 0.709466017509524, -0.109303314681054,  -0.453497173462763, 0.605887455840394, -1.81795596770373, 0.630098551068391,  -0.276184105225216, -0.284159743943371, -0.919322002474128, -0.116247806352002,  1.81731204370422, 0.370627864257954, 0.520216457554957, -0.750531994502331,  0.816899839520583, -0.886357521243213, -0.331577589942552), .Tsp = c(1,  9.5, 2), class = "ts") 
## ARIMA(0,0,0) with zero mean 
## 
## sigma^2 estimated as 0.6659:  log likelihood=-21.88
## AIC=45.76   AICc=46.01   BIC=46.65
# See forecasts from the auto.arima model
plot(forecast(hm1$auto.arima))

Model diagnostics

The hybridModel() function produces an S3 object of class forecastHybrid.

class(hm1) 
## [1] "hybridModel"
is.hybridModel(hm1)
## [1] TRUE

The print() and summary() methods print information about the ensemble model including the weights assigned to each individual component model.

print(hm1) 
## Hybrid forecast model comprised of the following models: auto.arima, ets, tbats
## ############
## auto.arima with weight 0.333 
## ############
## ets with weight 0.333 
## ############
## tbats with weight 0.333
summary(hm1)
##            Length Class  Mode     
## auto.arima 18     ARIMA  list     
## ets        19     ets    list     
## tbats      21     bats   list     
## weights     3     -none- numeric  
## frequency   1     -none- numeric  
## x          18     ts     numeric  
## xreg        1     -none- list     
## models      3     -none- character
## fitted     18     ts     numeric  
## residuals  18     ts     numeric

Two types of plots can be created for the created ensemble model: either a plot showing the actual and fitted value of each component model on the data or individual plots of the component models as created by their regular S3 plot() methods. Note that a plot() method does not exist in the “forecast” package for objects generated with stlm(), so this component model will be ignored when type = "models", but the other component models will be plotted regardless.

plot(quickModel, type = "fit")

plot(quickModel, type = "models")