README.Rmd

---
output: github_document
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

lgr::get_logger("mlr3")$set_threshold("warn")
lgr::get_logger("bbotk")$set_threshold("warn")
options(datatable.print.class = FALSE, datatable.print.keys = FALSE)
library(data.table)
library(mlr3misc)
```


# mlr3forecast

Extending mlr3 to time series forecasting.

<!-- badges: start -->
[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
[![RCMD Check](https://github.com/mlr-org/mlr3forecast/actions/workflows/rcmdcheck.yaml/badge.svg)](https://github.com/mlr-org/mlr3forecast/actions/workflows/rcmdcheck.yaml)
[![CRAN status](https://img.shields.io/badge/CRAN-not%20published-orange)](https://CRAN.R-project.org/package=mlr3forecast)
[![StackOverflow](https://img.shields.io/badge/stackoverflow-mlr3-orange.svg)](https://stackoverflow.com/questions/tagged/mlr3)
[![Mattermost](https://img.shields.io/badge/chat-mattermost-orange.svg)](https://lmmisld-lmu-stats-slds.srv.mwn.de/mlr_invite/)
<!-- badges: end -->

> This package is in an early stage of development and should be considered experimental.
> If you are interested in experimenting with it, we welcome your feedback!

## Installation

Install the development version from [GitHub](https://github.com/):

```{r, eval = FALSE}
# install.packages("pak")
pak::pak("mlr-org/mlr3forecast")
```

## Usage

mlr3forecast extends the [mlr3](https://mlr-org.com/) ecosystem to time series
forecasting. It introduces a forecasting task, forecasting learners, temporal
resampling strategies, forecasting measures, and feature-engineering pipe
operators, so that forecasters behave like any other mlr3 learner — ready for
tuning, benchmarking, pipelines, and ensembling.

At a glance, mlr3forecast provides:

- **Classical forecasters** wrapping forecast, smooth, prophet, and tscount —
  over 30 in total, from baselines (`fcst.mean`, `fcst.random_walk`) to ARIMA
  (`fcst.auto_arima`), ETS, theta, TBATS, prophet, neural nets (`fcst.nnetar`),
  and count models (`fcst.tscount`). See `mlr_learners$keys("^fcst")` for the
  full list.

- **Machine learning forecasting** that turns any `regr` learner into a
  forecaster via lag features, with both **recursive** (one model applied
  iteratively) and **direct** (one model per horizon) strategies.
- **Forecasting tasks and temporal resamplings** (`fcst.holdout`, `fcst.cv`)
  that respect the order of observations, plus global (longitudinal)
  forecasting across many series.
- **Feature-engineering pipe operators** such as `fcst.lags`, `fcst.rolling`,
  `fcst.fourier`, `fcst.feasts`, and `fcst.tsfeats`.
- **Forecasting measures** including MASE, RMSSE, Pinball, Winkler, coverage,
  and MSIS (see `mlr_measures$keys("^fcst")` for the full list).
- **Full mlr3 integration**: tuning with mlr3tuning, benchmarking, target
  transformations, and ensembling with mlr3pipelines.
- **Example tasks** to get started — `airpassengers`, `electricity`,
  `livestock`, `lynx`, and `usaccdeaths`. List them with
  `as.data.table(mlr_tasks)[task_type == "fcst"]`.

For now the forecasting task and learner are restricted to time series
regression, but may be extended to classification in the future.

Jump to the examples for:

- [Classical forecasters](#classical-forecasters)
- [Machine learning forecasters](#machine-learning-forecasters)
- [Benchmarking, ensembling, and tuning](#benchmarking-ensembling-and-tuning)
- [Global forecasting](#global-forecasting)

### Classical forecasters

Native forecasting learners are provided by packages such as forecast, smooth,
prophet, and tscount.

```{r, message = FALSE, dpi = 300}
library(mlr3forecast)
library(mlr3pipelines)

task = tsk("airpassengers")
task

# or plot the task
autoplot(task)

# train a forecast learner
learner = lrn("fcst.auto_arima")$train(task)
prediction = learner$predict(task, 140:144)
prediction
prediction$score(msr("regr.rmse"))
```

To forecast beyond the observed data, `generate_newdata()` builds the future
rows (with missing targets) and `predict_newdata()` fills them in:

```{r, message = FALSE}
# generate new data to forecast unseen data
newdata = generate_newdata(task, 12L)
head(newdata)
prediction = learner$predict_newdata(newdata, task)
prediction
```

The `forecast()` helper combines these two steps, generating the future rows and
predicting them in a single call:

```{r, message = FALSE}
forecast(learner, task, 12L)
```

Target transformations can be applied by wrapping the learner in
`ppl("targettrafo")`:

```{r, message = FALSE}
# add a target log transformation
learner = as_learner(ppl(
  "targettrafo",
  graph = lrn("fcst.auto_arima"),
  targetmutate.trafo = function(x) log(x),
  targetmutate.inverter = function(x) list(response = exp(x$response))
))
prediction = learner$train(task)$predict(task, 140:144)
prediction$score(msr("regr.rmse"))
```

Which predict types a learner supports (e.g. `"quantiles"`, `"se"`) is listed
in its `predict_types`:

```{r, message = FALSE}
lrn("fcst.auto_arima")$predict_types
```

Classical forecasters can then return a predictive distribution as quantiles,
scored with probabilistic measures such as the Pinball loss:

```{r, message = FALSE}
# works with quantile response
learner = lrn(
  "fcst.auto_arima",
  predict_type = "quantiles",
  quantiles = c(0.1, 0.15, 0.5, 0.85, 0.9),
  quantile_response = 0.5
)$train(task, 1:132)
learner$predict_newdata(newdata, task)
learner$predict(task, 133:144)$score(msr("fcst.pinball"))
```

Forecasting resamplings respect the temporal order of the observations:

```{r, message = FALSE}
# resampling, scored with a forecasting measure (MASE)
learner = lrn("fcst.auto_arima")
resampling = rsmp("fcst.holdout", ratio = 0.7)
rr = resample(task, learner, resampling)
rr$aggregate(msrs(c("regr.rmse", "fcst.mase")))
```

### Machine learning forecasters

Any regression learner can be turned into a forecaster with `recursive_forecaster()`,
which adds lag features and forecasts recursively:

```{r, message = FALSE}
library(mlr3learners)

task = tsk("airpassengers")
learner = lrn("regr.ranger")
flrn = recursive_forecaster(learner, lags = 1:12)$train(task)
newdata = generate_newdata(task, 12L)
prediction = flrn$predict_newdata(newdata, task)
prediction
prediction = flrn$predict(task, 140:144)
prediction
prediction$score(msr("regr.rmse"))

flrn = recursive_forecaster(learner, lags = 1:12)
resampling = rsmp("fcst.holdout", ratio = 0.9)
rr = resample(task, flrn, resampling)
rr$aggregate(msr("regr.rmse"))

resampling = rsmp("fcst.cv")
rr = resample(task, flrn, resampling)
rr$aggregate(msr("regr.rmse"))
```

#### Direct forecasting

`recursive_forecaster()` builds a recursive forecaster (one model, applied
iteratively). Use `direct_forecaster()` with `horizons` to train one model per
horizon instead — predictions then come straight from each horizon's model,
with no error accumulation:

```{r, message = FALSE}
task = tsk("airpassengers")
flrn = direct_forecaster(
  lrn("regr.ranger"),
  lags = 1:12,
  horizons = 12
)$train(task, 1:132)
flrn$predict(task, 133:144)$score(msr("regr.rmse"))
```

#### Feature engineering

Lag features can be combined with other transformations using mlr3pipelines:

```{r}
library(mlr3pipelines)

task = tsk("airpassengers")
task$set_col_roles("month", add = "feature")
graph = po("fcst.lags", lags = 1:12) %>>%
  po(
    "datefeatures",
    param_vals = list(
      week_of_year = FALSE,
      day_of_year = FALSE,
      day_of_month = FALSE,
      day_of_week = FALSE
    )
  ) %>>%
  lrn("regr.ranger")
flrn = recursive_forecaster(graph)$train(task)
prediction = flrn$predict(task, 142:144)
prediction$score(msr("regr.rmse"))
```

Use `selector_fcst_lags()` to apply transformations only to the lag features,
e.g. log-transforming lags while leaving date features untouched:

```{r}
task = tsk("airpassengers")
task$set_col_roles("month", add = "feature")
graph = po("fcst.lags", lags = 1:12) %>>%
  po("colapply", applicator = log, affect_columns = selector_fcst_lags()) %>>%
  po(
    "datefeatures",
    param_vals = list(
      week_of_year = FALSE,
      day_of_year = FALSE,
      day_of_month = FALSE,
      day_of_week = FALSE
    )
  ) %>>%
  lrn("regr.ranger")
flrn = recursive_forecaster(graph)$train(task)
prediction = flrn$predict(task, 142:144)
prediction$score(msr("regr.rmse"))
```

#### Target transformations

Target transformations can be applied by wrapping the forecast learner in
`ppl("targettrafo")`. The lags are created from the transformed target and
predictions are automatically inverted back to the original scale:

```{r}
task = tsk("airpassengers")
graph = po("fcst.lags", lags = 1:12) %>>% lrn("regr.ranger")
pipeline = ppl(
  "targettrafo",
  graph = recursive_forecaster(graph),
  targetmutate.trafo = function(x) log(x),
  targetmutate.inverter = function(x) list(response = exp(x$response))
)
learner = as_learner(pipeline)$train(task)
prediction = learner$predict(task, 142:144)
prediction$score(msr("regr.rmse"))
```

Ready-made `po("fcst.targetboxcox")` and `po("fcst.targetdiff")` pipeops are
also available for Box-Cox transformation and differencing.

#### Exogenous covariates

Forecasting tasks can include exogenous covariates. Here electricity demand is
forecast from its own lags, calendar features, and external regressors
(temperature, holiday) supplied for the forecast horizon:

```{r, message = FALSE}
library(mlr3learners)
library(mlr3pipelines)

task = tsk("electricity")
task$set_col_roles("date", add = "feature")
graph = po("fcst.lags", lags = 1:3) %>>%
  po("datefeatures", param_vals = list(year = FALSE)) %>>%
  lrn("regr.ranger")
flrn = recursive_forecaster(graph)$train(task)

max_date = task$data()[.N, date]
newdata = data.table(
  date = max_date + 1:14,
  demand = rep(NA_real_, 14L),
  temperature = 26,
  holiday = c(TRUE, rep(FALSE, 13L))
)
prediction = flrn$predict_newdata(newdata, task)
prediction
```

### Benchmarking, ensembling, and tuning

#### Comparing classical and ML forecasters

ML forecasters declare `task_type = "fcst"`, so they can be benchmarked
side-by-side with classical learners on the same task in a single
`benchmark()` call:

```{r, message = FALSE}
task = tsk("airpassengers")
resampling = rsmp("fcst.holdout", ratio = 0.9)$instantiate(task)
n_test = length(resampling$test_set(1L))

learners = list(
  lrn("fcst.arima", id = "arima"),
  recursive_forecaster(lrn("regr.ranger"), lags = 1:12, id = "ranger_recursive"),
  direct_forecaster(
    lrn("regr.ranger"),
    lags = 1:12,
    horizons = n_test,
    id = "ranger_direct"
  )
)
design = benchmark_grid(task, learners, resampling)
bmr = benchmark(design)
bmr$aggregate(msr("regr.rmse"))[, .(learner_id, regr.rmse)]
```

#### Ensemble forecasting

Forecast learners produce regression predictions under the hood, so the standard
mlr3pipelines ensemble pattern works directly: branch to several forecasters with
`gunion()` and average their forecasts with `po("regravg")`. This mirrors the idea
behind the forecastHybrid package, but with any mix of classical or ML learners.

```{r, message = FALSE}
task = tsk("airpassengers")

graph = gunion(list(
  po("learner", lrn("fcst.auto_arima"), id = "arima"),
  po("learner", lrn("fcst.ets"), id = "ets"),
  po("learner", lrn("fcst.theta"), id = "theta")
)) %>>%
  po("regravg")
flrn = as_learner(graph)$train(task)
forecast(flrn, task, 12L)
flrn$predict(task, 140:144)$score(msr("regr.rmse"))

# weight the members instead of averaging equally
graph$param_set$set_values(regravg.weights = c(0.5, 0.3, 0.2))
flrn = as_learner(graph)$train(task)
flrn$predict(task, 140:144)$score(msr("regr.rmse"))
```

#### Tuning a forecaster

Forecast learners are regular mlr3 learners, so they plug into the standard
[mlr3tuning](https://mlr3tuning.mlr-org.com/) machinery. Mark hyperparameters
with `to_tune()` and wrap the learner in an `auto_tuner()`, using a forecasting
resampling such as `fcst.holdout` or `fcst.cv` to respect the temporal order:

```{r, message = FALSE}
library(mlr3tuning)

task = tsk("airpassengers")

# tune an ML forecaster
flrn = recursive_forecaster(lrn("regr.ranger"), lags = 1:12)
flrn$param_set$set_values(
  regr.ranger.mtry.ratio = to_tune(0.1, 1),
  regr.ranger.num.trees = to_tune(100, 500)
)
at = auto_tuner(
  tuner = tnr("random_search"),
  learner = flrn,
  resampling = rsmp("fcst.cv"),
  measure = msr("regr.rmse"),
  term_evals = 4
)
at$train(task)
at$tuning_result[, .(regr.ranger.mtry.ratio, regr.ranger.num.trees, regr.rmse)]

# the AutoTuner is itself a learner: predict with the best configuration
at$predict(task, 142:144)$score(msr("regr.rmse"))
```

Classical forecasters tune the same way:

```{r, message = FALSE}
flrn = lrn("fcst.auto_arima")
flrn$param_set$set_values(stationary = to_tune(p_lgl()), seasonal = to_tune(p_lgl()))
at = auto_tuner(
  tuner = tnr("grid_search"),
  learner = flrn,
  resampling = rsmp("fcst.holdout", ratio = 0.8),
  measure = msr("regr.rmse")
)
at$train(task)
at$tuning_result[, .(stationary, seasonal, regr.rmse)]
```

### Global forecasting

In machine learning forecasting the difference between forecasting a single time
series and longitudinal data is often referred to as local and global
forecasting. A global model is trained jointly across many series, identified by
a `key`:

```{r, message = FALSE}
library(mlr3learners)
library(mlr3pipelines)
library(tsibble)

dt = setDT(tsibbledata::aus_livestock)
setnames(dt, tolower)
dt[, month := as.Date(month)]
dt = dt[, .(count = sum(count)), by = .(state, month)]
setorder(dt, state, month)
task = as_task_fcst(dt, id = "aus_livestock", target = "count", order = "month", key = "state", freq = "month")
task$set_col_roles("month", add = "feature")

graph = po("fcst.lags", lags = 1:12) %>>%
  po(
    "datefeatures",
    param_vals = list(
      week_of_year = FALSE,
      day_of_week = FALSE,
      day_of_month = FALSE,
      day_of_year = FALSE
    )
  ) %>>%
  lrn("regr.ranger")

flrn = recursive_forecaster(graph)$train(task)
prediction = flrn$predict(task, 4460:4464)
prediction$score(msr("regr.rmse"))

resampling = rsmp("fcst.holdout", ratio = 0.9)
rr = resample(task, flrn, resampling)
rr$aggregate(msr("regr.rmse"))
```

#### Global vs. local forecasting

A single global model can be compared against fitting one local model per series:

```{r, message = FALSE}
retail = setDT(tsibbledata::aus_retail)
setnames(retail, tolower)
retail[, month := as.Date(month)]
vic = retail[state == "Victoria"]
vic[, let(state = NULL, `series id` = NULL)]
vic[, industry := as.factor(industry)]
vic_train = vic[month < as.Date("2015-01-01")]
vic_test = vic[month >= as.Date("2015-01-01")]

# global forecasting
task_train = as_task_fcst(
  vic_train,
  id = "aus_retail_vic",
  target = "turnover",
  order = "month",
  key = "industry",
  freq = "month"
)
task_test = as_task_fcst(
  vic_test,
  id = "aus_retail_vic",
  target = "turnover",
  order = "month",
  key = "industry",
  freq = "month"
)

learner = lrn("regr.ranger", verbose = FALSE)
flrn = recursive_forecaster(learner, lags = 1:12)$train(task_train)
prediction_global = flrn$predict(task_test)
prediction_global
prediction_global$score(msr("regr.rmse"))

# local forecasting
prediction_local = map(split(vic, by = "industry", drop = TRUE), function(dt) {
  task_train = as_task_fcst(
    dt[month < as.Date("2015-01-01")],
    id = "aus_retail_vic_local",
    target = "turnover",
    order = "month",
    freq = "month"
  )
  task_test = as_task_fcst(
    dt[month >= as.Date("2015-01-01")],
    id = "aus_retail_vic_local",
    target = "turnover",
    order = "month",
    freq = "month"
  )
  flrn = recursive_forecaster(learner, lags = 1:12)$train(task_train)
  prediction = flrn$predict(task_test)
  prediction
})
do.call(c, prediction_local)$score(msr("regr.rmse"))
```