Forecasting in mlr

Forecasting

The standard objective in forecasting is, at time period $t$, make predictions for $t+h$ periods into the future. Forecasting tasks are most suitable when past patterns in the data will continue on into the future. While R has many wonderful packages for forecasting, there is no package on CRAN that gives users a standardized form for forecasting models, such as mlr and caret for machine learning models. The purpose of this package extension is to give users of mlr the opportunity to safely and productively train, optimize, and deploy forecasting models.

Forecasting Tasks

For our purposes we will simulate some data.

set.seed(1234)
dat = arima.sim(model = list(ar = c(.5,.2), ma = c(.4), order = c(2,0,1)), n = 200)
times = (as.POSIXlt("1992-01-14")) + lubridate::days(1:200)
dat = xts::xts(dat,order.by = times)
colnames(dat) = c("arma_test")
dat.train = dat[1:190,]
dat.test = dat[191:200,]

To specify a forecast task we must pass an xts object containing the data, a target, and the frequency of the data. The frequency of the data can be thought of as the seasonality. For example, a frequency of seven on daily data would be a weekly seasonality. A frequency of fifty-two on weekly data would indicate a yearly seasonality.

library(mlr)

## Loading required package: ParamHelpers

fcregr.task = makeForecastRegrTask(id = "test", data = dat.train, target = "arma_test",
                                     frequency = 7L)

## Warning in makeForecastRegrTask(id = "test", data = dat.train, target =
## "arma_test", : Provided data for task is not a pure data.frame but from
## class xts, hence it will be converted.

fcregr.task

## Task: test
## Type: fcregr
## Target: arma_test
## Observations: 190
## Dates:
##  Start: 1992-01-15 
##  End:   1992-07-22
## Frequency: 7
## Features:
##    numerics     factors     ordered functionals 
##           0           0           0           0 
## Missings: FALSE
## Has weights: FALSE
## Has blocking: FALSE

Like a regression task, this records the type of the learning problem and basic information about the data set. This task also returns the start and end dates of the time series as well as the frequency.

Multivariate Forecasting Tasks

One common problem with forecasting is that it is difficult to use additional explanatory variables or forecast multiple targets that are dependent on one another. If we are at time $t$ and want to forecast 10 periods in the future, we need to know the values of the explanatory variables at time $t+10$, which is often not possible. A new set of models which treats explanatory variables endogenously instead of exogenously allows us to forecast not only our target, but additional explanatory variables. This is done by treating all the variables as targets, making them endogenous to the model. To use these models, we create a multivariate forecasting task.

data("EuStockMarkets")
EuStockMarkets.time = lubridate::date_decimal(as.numeric(time(EuStockMarkets)))
EuStockMarkets  = xts::xts(as.data.frame(EuStockMarkets), order.by = EuStockMarkets.time)
EuStockMarkets.train = EuStockMarkets[1:1850,]
EuStockMarets.test = EuStockMarkets[1851:1860,]
mfcregr.task = makeMultiForecastRegrTask(id = "bigvar", data = EuStockMarkets.train, target = "all", frequency = 7L)

## Warning in makeMultiForecastRegrTask(id = "bigvar", data =
## EuStockMarkets.train, : Provided data for task is not a pure data.frame but
## from class xts, hence it will be converted.

mfcregr.task

## Task: bigvar
## Type: mfcregr
## Target: DAX SMI CAC FTSE
## Observations: 1850
## Dates:
##  Start: 1991-07-01 02:18:27 
##  End:   1998-08-10 19:23:04
## Frequency: 7
## Features:
##    numerics     factors     ordered functionals 
##           0           0           0           0 
## Missings: FALSE
## Has weights: FALSE
## Has blocking: FALSE

This task prints the same information as a univariate forecast task.

If we would like to specify a target variable while still forecasting the rest of the data in an endogenous manner, we change our target to one of the variables.

mfcregr.task = makeMultiForecastRegrTask(id = "bigvar", data = EuStockMarkets.train, target = "DAX", frequency = 7L)

## Warning in makeMultiForecastRegrTask(id = "bigvar", data =
## EuStockMarkets.train, : Provided data for task is not a pure data.frame but
## from class xts, hence it will be converted.

mfcregr.task

## Task: bigvar
## Type: mfcregr
## Target: DAX
## Observations: 1850
## Dates:
##  Start: 1991-07-01 02:18:27 
##  End:   1998-08-10 19:23:04
## Frequency: 7
## Features:
##    numerics     factors     ordered functionals 
##           3           0           0           0 
## Missings: FALSE
## Has weights: FALSE
## Has blocking: FALSE

Making Learners for Forecasting

Several new models have been included from the forecast package and well as rugarch:

1. Exponential smoothing state space model with Box-Cox transformation (bats)
2. Exponential smoothing state space model with Box-Cox transformation, ARMA errors, Trend and Seasonal Fourier components (tbats)
3. Exponential smoothing state space model (ets)
4. Neural Network Autoregressive model (nnetar)
5. Automated Arima (auto.arima)
6. General Autoregressive Conditional Heteroskedasticity models (GARCH)
7. BigVar for multivariate time series

These all operate the same as the other models in mlr, with a very important parameter. Models will either have an h or n.ahead parameter, which is the number of periods you want to forecast into the future. Note that this should be equal to the horizon you set in your growing or fixed window resampling strategy.

We create a GARCH model from the package rugarch using makeLearner() by calling the learner class fcregr.garch. An important parameter is the n.ahead parameter, which is used to specify that we are forecasting 10 periods into the future. We can view all of the possible parameters that can be tuned with getLearnerParamSet("fcregr.garch")

garch.mod = makeLearner("fcregr.garch",
                        model = "sGARCH", garchOrder = c(2,1),
                        n.ahead = 10L, include.mean = FALSE)
garch.mod

## Learner fcregr.garch from package rugarch
## Type: fcregr
## Name: Generalized AutoRegressive Conditional Heteroskedasticity; Short name: garch
## Class: fcregr.garch
## Properties: numerics,quantile
## Predict-Type: response
## Hyperparameters: model=sGARCH,garchOrder=2,1,n.ahead=10,include.mean=FALSE

Once we set up our task and model train(), predict(), and performance() can be called to build and evaluate the model.

garch.train = train(garch.mod, fcregr.task)
garch.pred = predict(garch.train, newdata = dat.test)

## Warning in predict.WrappedModel(garch.train, newdata = dat.test): Provided
## data for prediction is not a pure data.frame but from class xts, hence it
## will be converted.

performance(garch.pred, measure = mase, task = fcregr.task, model = garch.train)

##      mase 
## 0.9329836

This standard evaluation method is user friendly and in addition, we can now use mlr’s built in resampling and tuning methods to tune a garch model for the data.

Resampling

mlr now has two new cross validation resampling strategies, GrowingCV and FixedCV. They are both rolling forecasting origin techniques established in Hyndman and Athanasopoulos (2013) and first widely available for machine learning in R by the caret package’s createTimeSlices() function. We specify:

1. horizon - the number of periods to forecast
2. initial.window - The proportion of data that will be used in the initial window
3. size - The number of rows in the training set
4. skip - the proportion of windows to skip over, which can be used to save time

resamp.desc = makeResampleDesc("GrowingCV", horizon = 10L,
                               initial.window = .90)
resamp.desc

## Window description:
##  Growing: 90.00 % in initial window, horizon of 10, and skipping 9 windows.
## Predict: test
## Stratification: FALSE

Note that we should not specify stratification, as it does not really make sense in the context of time series to stratify our data (unless a future implementation can somehow use this for panel data). The wonderful graphic posted below comes from the caret package’s website and gives an intuitive idea of the sliding windows for both the growing and fixed options.

The top two graphs correspond to fixed window while the bottom two graphs are for growing window

Using the resampling strategy we can do a windowing cross-validation with our previous GARCH model to evaluate the general performance of the created GARCH model.

garch.resample = resample(learner = garch.mod, task = fcregr.task,
                          resampling = resamp.desc, measures = mase)

## Resampling: Growing

## Measures:             mase

## [Resample] iter 1:    0.7426964

## 

## Aggregated Result: mase.test.mean=0.7426964

## 

garch.resample

## Resample Result
## Task: test
## Learner: fcregr.garch
## Aggr perf: mase.test.mean=0.7426964
## Runtime: 1.08794

Tuning

The forecasting features fully integrate into mlr, allowing us to also make a parameter set to tune over. Here we tune a GARCH model with F1-racing used to tune our parameters. In addition, we can tune the model in parallel using parallelMap.

par_set = makeParamSet(
  makeDiscreteParam(id = "model", values = c("sGARCH", "fGARCH")),
  makeDiscreteParam(id = "submodel", values = c("GARCH", "NGARCH"),
                    requires = quote(model == "fGARCH")),
  makeIntegerVectorParam(id = "garchOrder", len = 2L, lower = 1, upper = 3))

#Specify tune by Irace
ctrl = makeTuneControlIrace(maxExperiments = 180L)

library(parallelMap)
parallelStartSocket(6)
configureMlr(on.learner.error = "warn")
res = tuneParams(garch.mod, task = fcregr.task,
                 resampling = resamp.desc, par.set = par_set,
                 control = ctrl, measures = mase)
parallelStop()
res

## Tune result:
## Op. pars: model=fGARCH; submodel=GARCH; garchOrder=1,2
## mase.test.mean=1.8285587

Once we have tuned the model we get back the hyperparameters of the best model through setHyperPars() and then train the model like any other.

garch.hyp  = setHyperPars(makeLearner("fcregr.garch", n.ahead = 10L),
                          par.vals = res$x)
garch.best = train(garch.hyp, fcregr.task)
garch.pred = predict(garch.best, newdata = dat.test)

## Warning in predict.WrappedModel(garch.best, newdata = dat.test): Provided
## data for prediction is not a pure data.frame but from class xts, hence it
## will be converted.

performance(garch.pred, measures = mase, task = fcregr.task, model = garch.train)

##      mase 
## 0.9595138

Updating Models

A new function updateModel() has been implemented that updates the model given new data. This function is currently only implemented for ets, Arima, auto.arima, bats, tbats, and nnetar.

fcregr.task = makeForecastRegrTask(id = "test", data = dat.train,
                                   target = "arma_test", frequency = 7L)

## Warning in makeForecastRegrTask(id = "test", data = dat.train, target =
## "arma_test", : Provided data for task is not a pure data.frame but from
## class xts, hence it will be converted.

arm = makeLearner("fcregr.Arima", order = c(2L,0L,1L),
                  h = 10L, include.mean = FALSE)
arm

## Learner fcregr.Arima from package forecast
## Type: fcregr
## Name: AutoRegressive Integrated Moving Average; Short name: Arima
## Class: fcregr.Arima
## Properties: numerics,quantile
## Predict-Type: response
## Hyperparameters: order=2,0,1,h=10,include.mean=FALSE

armMod = train(arm, fcregr.task)
updateArmMod = updateModel(armMod, fcregr.task, newdata = dat.test)

## Warning in updateModel(armMod, fcregr.task, newdata = dat.test): Provided
## data for update is not a pure data.frame but from class function, hence it
## will be converted.

## Warning in max(j, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

updateArmMod

## Model for learner.id=fcregr.Arima; learner.class=fcregr.Arima
## Trained on: task.id = test; obs = 10; features = 0
## Hyperparameters: order=2,0,1,h=10,include.mean=FALSE

This works by making a call to updateLearner.fcregr.Arima() and updating the model and task data with newdata. predict() now forcasts the next 10 periods.

predict(updateArmMod, task = fcregr.task)

## Prediction: 10 observations
## predict.type: response
## threshold: 
## time: 0.00
##             response
## 1992-07-23 1.3631309
## 1992-07-24 1.0226897
## 1992-07-25 0.7818986
## 1992-07-26 0.5996956
## 1992-07-27 0.4601914
## 1992-07-28 0.3531699
## ... (#rows: 10, #cols: 1)

Pre-processing

library(xts)
library(lubridate)
date.col = as.POSIXct(index(dat.train))
dat.train.reg = as.data.frame(dat.train)
regr.task = makeRegrTask(data = dat.train.reg, target = "arma_test")
regr.task.lag = createLagDiffFeatures(regr.task,lag = 2L:4L, difference = 1L, 
                                          seasonal.lag = 1L:2L, frequency = 5, date.col = date.col)
regr.task.lag

## Supervised task: dat.train.reg
## Type: regr
## Target: arma_test
## Observations: 180
## Features:
##    numerics     factors     ordered functionals 
##           7           0           0           0 
## Missings: FALSE
## Has weights: FALSE
## Has blocking: FALSE
## Is spatial: FALSE

regrGbm <- makeLearner("regr.gbm", par.vals = list(n.trees = 100))
gbm.mod = train(regrGbm, regr.task.lag)

gbm.fore = forecast(gbm.mod, h = 10, newdata = dat.test)

## Warning in forecast.WrappedModel(gbm.mod, h = 10, newdata = dat.test):
## Provided data for prediction is not a pure data.frame but from class xts,
## hence it will be converted.

performance(gbm.fore, measures = mase, task = regr.task.lag, model = gbm.mod)

##     mase 
## 2.066958

Classification Forecasting Example

For developing trading strategies, we normally have a discrete set of choices such as to buy, sell, or hold onto a stock. Using forecasting in mlr we can now train classification models that forecast these choices. He we try to predict whether a ‘stock’ will go up 5%, down 5%, or neither. Below we simulate data from an ARIMA process, take the percentage change, and descritize the data into either a ‘buy’, ‘sell’, or ‘hold’.

library(xts)
library(lubridate)
set.seed(1234)
# Generate an ARIMA based data set
dat = arima.sim(model = list(ar = c(.5,.2,.1), ma = c(.5,.3), order = c(3,0,2)), n = 500)
dat = dat/lag(dat,-1) - 1 
jump = data.frame(jump = ifelse(diff(dat) > .5, "Buy",ifelse(diff(dat) < -.5, "Sell","Hold")))
times = (as.POSIXct("1992-01-14")) + lubridate::days(1:498)
rownames(jump) = times
jump.train = jump[1:488,,drop = FALSE]
jump.test  = jump[489:498,,drop = FALSE]

# Make the classif task
classif.task = makeClassifTask(data = jump.train,target = "jump")
classif.task.lag = createLagDiffFeatures(classif.task, lag = 1L:15L, 
                                         date.col = times[1:488])

classif.learn = makeLearner("classif.boosting")
classif.train = train(classif.learn,classif.task.lag)
classif.fore = forecast(classif.train, h = 10, newdata = jump.test)
performance(classif.fore)

## mmce 
##  0.7

So our model does pretty bad, but that is because it was not tuned. Given the structure in mlr, parallel computing, and time, it would be possible to use this type of model to make forecast classifiers.

Ensembles of Forecasts

It’s known that ensembles of forecasts tend to outperform standard forecasting techniques. Here we use mlr’s stacked modeling functionality to ensemble multiple forecast techniques.

# Create some data
library(xts)
library(lubridate)
dat = arima.sim(model = list(ar = c(.5,.2), ma = c(.4), order = c(2,0,1)), n = 500)
times = (as.POSIXlt("1992-01-14")) + lubridate::days(1:500)
dat = xts(dat,order.by = times)
colnames(dat) = c("arma_test")

# seperate into test and train
dat.train = dat[1:490,]
dat.test  = dat[491:500,]

fcregr.task = makeForecastRegrTask(id = "test", data = dat.train, target = "arma_test",
                                     frequency = 7L)

## Warning in makeForecastRegrTask(id = "test", data = dat.train, target =
## "arma_test", : Provided data for task is not a pure data.frame but from
## class xts, hence it will be converted.

resamp.sub = makeResampleDesc("GrowingCV",
                          horizon = 10L,
                          initial.window = .90
                          )
                          
resamp.super = makeResampleDesc("CV", iters = 5)

lrns = makeLearners(c("fcregr.tbats","fcregr.garch"))

                            
stack.forecast = makeStackedLearner(base.learners = lrns,
                       predict.type = "response",
                       method = "average")

# Simple param set for tuning sub learners
ps = makeParamSet(
  makeDiscreteParam("fcregr.tbats.h", values = 10),
  makeDiscreteParam("fcregr.garch.n.ahead", values = 10)
)

## tuning
fore.tune = tuneParams(stack.forecast, fcregr.task, resampling = resamp.sub,
                   par.set = ps, control = makeTuneControlGrid(),
                   measures = mase, show.info = FALSE)
fore.tune

## Tune result:
## Op. pars: fcregr.tbats.h=10; fcregr.garch.n.ahead=10
## mase.test.mean=0.6834574

# get hyper params
stack.forecast.tune  = setHyperPars2(stack.forecast,fore.tune$x)
stack.forecast.tune

## Learner stack from package forecast,rugarch
## Type: fcregr
## Name: ; Short name: 
## Class: StackedLearner
## Properties: numerics,quantile
## Predict-Type: response
## Hyperparameters: fcregr.tbats.h=10,fcregr.garch.n.ahead=10

# Train the final best models and predict
stack.forecast.mod = train(stack.forecast.tune,fcregr.task)
stack.forecast.pred = predict(stack.forecast.mod, newdata = dat.test)

## Warning in predict.WrappedModel(stack.forecast.mod, newdata = dat.test):
## Provided data for prediction is not a pure data.frame but from class xts,
## hence it will be converted.

stack.forecast.pred

## Prediction: 10 observations
## predict.type: response
## threshold: 
## time: 0.00
##                truth  response
## 1993-05-19 0.5473146 0.3755373
## 1993-05-20 0.3629478 0.2787431
## 1993-05-21 0.3217870 0.2242204
## 1993-05-22 1.6756571 0.1845821
## 1993-05-23 1.7294039 0.1557456
## 1993-05-24 0.8262623 0.1347527
## ... (#rows: 10, #cols: 2)

Please note that this an experimental branch of mlr. Please direct any issues to mlr’s issue tracker on github