Univariate time series forecasting (autoregressive models) and recursive forecasts

A univariate time series forecasting method is one in which the future values of a series are predicted using only information from the series itself—for example, using its historical mean value as a forecast. An advantage of this type of prediction is that, apart from the series being forecasted, there is no need to collect any further information to train the forecasting model.

An autoregressive model is a type of univariate time series forecasting model in which a value of a time series is expressed as a function of some of its past values. That is, an autoregressive model is a regression model in which the independent variables are lagged values (previous values) of the response variable. For example, given a time series with the following historical values: \(t = \{1, 3, 6, 7, 9, 11, 16\}\), suppose that we want to develop an autoregressive model in which a target “is explained” by its first, second, and fourth past values (in this context, a previous value is also called a lag, so lag 1 is the value immediately preceding a given value in the series). Given the series \(t\) and lags (1, 2, and 4), the training set would be:

The utsf package

The utsf package simplifies the use of classical regression models for univariate time series forecasting, employing the autoregressive approach and the recursive prediction strategy explained in the previous section. All the supported models are applied using a uniform interface:

• the create_model() function to build the forecasting model, and
• the forecast() function for using the trained model to predict future values.

Let us see an example in which a regression tree model is used to forecast the future values of a time series:

m <- create_model(AirPassengers, lags = 1:12, method = "rt")
f <- forecast(m, h = 12)

In this example, an autoregressive tree model (method = "rt", Regression Tree) is trained using the historical values of the AirPassengers time series and a forecast for its next 12 values (h = 12) is generated. The create_model() function returns an S3 object of class utsf with information about the trained model. The forecast information is included in the object returned by the forecast() function, as a component named pred of class ts (a time series):

f$pred
#>           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
#> 1961 460.2980 428.4915 467.0304 496.9833 499.9819 554.7891 627.5849 628.0503
#>           Sep      Oct      Nov      Dec
#> 1961 533.2803 482.4221 448.7926 453.6920
library(ggplot2)
autoplot(f)

The training set used to fit the model is built from the historical values of the time series using the autoregressive approach explained in the previous section. The lags parameter of the create_model() function specifies the autoregressive lags. In this example: lags = 1:12, meaning each target is a function of its 12 previous values. Next, we check the first targets (and their associated features) used to train the regression model:

head(m$targets)  # first targets
#> [1] -11.6666667  -0.9166667  13.4166667   6.6666667  -3.8333333  19.8333333
head(m$features) # and its associated features
#>         Lag12     Lag11     Lag10      Lag9       Lag8       Lag7       Lag6
#> 1 -14.6666667 -8.666667  5.333333  2.333333  -5.666667   8.333333  21.333333
#> 2  -8.9166667  5.083333  2.083333 -5.916667   8.083333  21.083333  21.083333
#> 3   4.4166667  1.416667 -6.583333  7.416667  20.416667  20.416667   8.416667
#> 4   0.6666667 -7.333333  6.666667 19.666667  19.666667   7.666667  -9.333333
#> 5  -7.8333333  6.166667 19.166667 19.166667   7.166667  -9.833333 -24.833333
#> 6   5.8333333 18.833333 18.833333  6.833333 -10.166667 -25.166667 -11.166667
#>         Lag5       Lag4       Lag3       Lag2       Lag1
#> 1  21.333333   9.333333  -7.666667 -22.666667  -8.666667
#> 2   9.083333  -7.916667 -22.916667  -8.916667 -11.916667
#> 3  -8.583333 -23.583333  -9.583333 -12.583333  -1.583333
#> 4 -24.333333 -10.333333 -13.333333  -2.333333  12.666667
#> 5 -10.833333 -13.833333  -2.833333  12.166667   6.166667
#> 6 -14.166667  -3.166667  11.833333   5.833333  -4.166667

The curious reader may have noticed that the features and targets are not on the same scale as the original time series. This is because, by default, a transformation is applied to the examples of the training set. This transformation will be explained later.

Let us inspect the training set associated with the example from the previous section:

t <- ts(c(1, 3, 6, 7, 9, 11, 16))
out <- create_model(t, lags = c(1, 2, 4), trend = "none")
cbind(out$features, Target = out$targets)
#>   Lag4 Lag2 Lag1 Target
#> 1    1    6    7      9
#> 2    3    7    9     11
#> 3    6    9   11     16

No transformation has been applied this time (trend = "none").

Prediction intervals

Prediction intervals for a forecast can be computed:

m <- create_model(USAccDeaths, method = "mt")
f <- forecast(m, h = 12, PI = TRUE, level = 90)
f
#>          Point Forecast    Lo 90     Hi 90
#> Jan 1979       8195.911 7654.387  8750.039
#> Feb 1979       7446.227 6885.823  8024.522
#> Mar 1979       7775.062 7072.970  8419.338
#> Apr 1979       8289.710 7645.091  8977.396
#> May 1979       9048.237 8329.607  9716.592
#> Jun 1979       9237.905 8561.585  9916.638
#> Jul 1979      10297.351 9014.763 10837.597
#> Aug 1979       9741.415 8632.913 10399.199
#> Sep 1979       9103.687 8274.895  9827.232
#> Oct 1979       8948.850 8141.232  9703.285
#> Nov 1979       8750.469 8083.375  9465.097
#> Dec 1979       9247.210 7985.251  9815.854
library(ggplot2)
autoplot(f)

Prediction intervals are calculated following the guidelines in Forecasting: Principles and Practice. Random errors are assumed to follow a normal distribution. In this example, a 90% prediction interval (level = 90) was computed.

Supported models

The create_model() and forecast() functions provide a unified interface for applying an autoregressive approach to time series forecasting with various regression models. These models are implemented across several R packages. Currently, our project focuses primarily on regression tree models, supporting the following approaches:

• k-nearest neighbors: In this case no model is trained and the function FNN::knn.reg() is used, as regression function, to recursively predict the future values of the time series.
• Linear models: The model is trained using the function stats::lm() and its associated method stats::predict.lm() is applied recursively for the forecasts, i.e., as regression function.
• Regression trees: The model is trained using the function rpart::rpart() and its associated method rpart::predict.rpart() is used for the forecasts.
• Model trees: The model is trained with the function Cubist::cubist() and its associated method Cubist::predict.cubist() is used for predictions.
• Bagging: The model is trained with the function ipred::bagging() and its associated method ipred::predict.regbagg() is used for forecasting.
• Random forest: The model is trained with the function ranger::ranger() and its associated method ranger::predict.ranger() is used for predictions.
• Extreme gradient boosting: The model is trained with the function xgboost::xgboost() and its associated method xgboost::predict.xgboost() is used for predictions.

The S3 object of class utsf returned by the create_model() function contains a component with the trained autoregressive model:

m <- create_model(fdeaths, lags = 1:12, method = "rt")
m$model
#> n= 60 
#> 
#> node), split, n, deviance, yval
#>       * denotes terminal node
#> 
#> 1) root 60 1967124.00   -6.465278  
#>   2) Lag12< 73 38  212851.90 -124.414500  
#>     4) Lag6>=-66.45833 30   57355.07 -153.847200  
#>       8) Lag12< -170.7083 10   13293.09 -195.991700 *
#>       9) Lag12>=-170.7083 20   17419.67 -132.775000 *
#>     5) Lag6< -66.45833 8   32051.01  -14.041670 *
#>   3) Lag12>=73 22  312482.00  197.265200  
#>     6) Lag5>=-131.7917 7   24738.12  114.500000 *
#>     7) Lag5< -131.7917 15  217416.40  235.888900 *

In this case, the model is obtained by training a regression tree using the rpart::rpart() function on a training set composed of the features m$features and targets m$targets. Once trained, the rpart::predict.rpart() function can be used recursively within forecast() to predict future values of the time series.

Using your own models

An appealing feature of the utsf package is its extensibility. Beyond the natively supported models, you can leverage the package’s infrastructure to perform autoregressive time series forecasting with custom regression models. This allows user-defined models to benefit from built-in functionality, such as training set construction, recursive forecasting, prediction intervals, preprocessing, accuracy evaluation, and hyperparameter tuning.

To use a custom regression model, pass a function capable of training your model to the method parameter of the create_model() function. We refer to this function as the modeling function, and it must return an object containing the trained regression model. The modeling function must accept three input parameters:

• X: a data frame containing the features of the training examples. This data frame is constructed from the time series according to the autoregressive lags, as explained in a previous section. It is identical to the features component of the object returned by create_model().

• y: a vector containing the targets of the training examples, constructed as described in a previous section. It is identical to the targets component of the object returned by create_model().

• ...: additional arguments passed to fine-tune model behavior. This parameter is not used in this section and will be detailed in the next.

Furthermore, if the modeling function returns an object of class model_class, an S3 method with the signature predict.model_class(object, new_value) must be implemented. This method uses the model to predict a single new value (acting as the prediction function associated with the model). Its parameters are described below:

• object: the object of class model_class returned by your modeling function.
• new_value: a one-row data frame with the same structure as X, containing the features of the instance to be predicted.

In summary, to use a custom regression model, you must provide:

• A modeling function that trains the model given a training set.
• A prediction function that predicts a single target value given a trained model and a new feature vector.

The following examples demonstrate how to use this functionality to apply custom models to autoregressive time series forecasting, leveraging the capabilities of the utsf package.

# Modeling function: in this case (k-NN) just stores the training set
my_knn_model <- function(X, y, ...) { structure(list(X = X, y = y), class = "my_knn") }

# Regression function
predict.my_knn <- function(object, new_value) {
  FNN::knn.reg(train = object$X, test = new_value, y = object$y)$pred
}

m <- create_model(AirPassengers, lags = 1:12, method = my_knn_model)
f <- forecast(m, h = 12)
f$pred
#>           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
#> 1961 455.9167 434.3264 480.7703 490.1678 506.1262 568.0534 640.6689 640.8636
#>           Sep      Oct      Nov      Dec
#> 1961 549.5467 495.4255 441.6554 476.7934
autoplot(f)

# Modeling function
my_knn_model2 <- function(X, y, ...) {
  structure(list(X = X, y = y), class = "my_knn2")
}

# Regression function
predict.my_knn2 <- function(object, new_value) {
  k <- 3 # number of nearest neighbors
  distances <- sapply(1:nrow(object$X), function(i) sum((object$X[i, ] - new_value)^2))
  k_nearest <- order(distances)[1:k]
  mean(object$y[k_nearest])
}

m2 <- create_model(AirPassengers, lags = 1:12, method = my_knn_model2)
forecast(m2, h = 12)$pred
#>           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
#> 1961 455.9167 434.3264 480.7703 490.1678 506.1262 568.0534 640.6689 640.8636
#>           Sep      Oct      Nov      Dec
#> 1961 549.5467 495.4255 441.6554 476.7934

library(randomForest)
#> randomForest 4.7-1.2
#> Type rfNews() to see new features/changes/bug fixes.
#> 
#> Adjuntando el paquete: 'randomForest'
#> The following object is masked from 'package:ggplot2':
#> 
#>     margin
my_model <- function(X, y, ...) { randomForest(x = X, y = y) }
m <- create_model(USAccDeaths, lags = 1:12, method = my_model)
f <- forecast(m, h = 12)
library(ggplot2)
autoplot(f)

print(m$model)
#> 
#> Call:
#>  randomForest(x = X, y = y) 
#>                Type of random forest: regression
#>                      Number of trees: 500
#> No. of variables tried at each split: 4
#> 
#>           Mean of squared residuals: 190450.2
#>                     % Var explained: 77.39

library(nnet)
my_model <- function(X, y, ...) {
  nnet(x = X, y = y, size = 5, linout = TRUE, trace = FALSE)
}
m <- create_model(USAccDeaths, lags = 1:12, method = my_model)
f <- forecast(m, h = 12)
library(ggplot2)
autoplot(f)

Setting the parameters of the regression models

Typically, a regression model can be tuned using hyperparameters. By default, the models supported by the package use specific default values defined by their underlying training functions (listed in a previous section). However, users can customize these hyperparameters to adjust model behavior. Below is an example:

# A bagging model set with default parameters
m <- create_model(USAccDeaths, lags = 1:12, method = "bagging")
length(m$model$mtrees) # number of regression trees (25 by default)
#> [1] 25
# A bagging model set with 3 regression trees
m2 <- create_model(USAccDeaths,
                   lags = 1:12,
                   method = "bagging",
                   nbagg = 3
)
length(m2$model$mtrees) # number of regression trees
#> [1] 3

In the previous example, two bagging models (using regression trees) are trained with create_model(). The first model uses 25 trees (the default for ipred::ipredbagg()) while the second sets the number of trees to 3. To customize hyperparameters, consult the parameter documentation of the training function used internally by utsf. In this case, ipred::ipredbagg() uses the nbagg parameter to specify the number of trees in the ensemble.

# Function to train the model
my_knn_model <- function(X, y, ...) {
  param <- list(...)
  k <- if ("k" %in% names(param)) param$k else 3
  structure(list(X = X, y = y, k = k), class = "my_knn")
}

# Regression function for object of class my_knn
predict.my_knn <- function(object, new_value) {
  FNN::knn.reg(train = object$X, test = new_value,
               y = object$y, k = object$k)$pred
}

# The model is trained with default parameters (k = 3)
m <- create_model(AirPassengers, lags = 1:12,  method = my_knn_model)
print(m$model$k)
#> [1] 3
# The model is trained with k = 5
m2 <- create_model(AirPassengers, method = my_knn_model, k = 5)
print(m2$model$k)
#> [1] 5

library(randomForest)
my_model <- function(X, y, ...) {
    param <- list(...)
    args <- list(x = X, y = y, ntree = 200) # default parameters
    args <- args[!(names(args) %in% names(param))]
    args <- c(args, param)
    do.call(randomForest::randomForest, args = args)
}
# The random forest is built with our default parameters
m <- create_model(USAccDeaths, method = my_model)
print(m$model$ntree)
#> [1] 200
print(m$model$mtry)
#> [1] 4
# The random forest is built with ntree = 400 and mtry = 6
m <- create_model(USAccDeaths, method = my_model, ntree = 400, mtry = 6)
print(m$model$ntree)
#> [1] 400
print(m$model$mtry)
#> [1] 6

Estimating forecast accuracy

This section demonstrates how to estimate the forecast accuracy of a regression model when predicting a time series with utsf. Below is an example:

m <- create_model(UKgas, lags = 1:4, method = "knn")
r <- efa(m, h = 4, type = "normal", size = 8)
r$per_horizon
#>       Horizon 1 Horizon 2 Horizon 3 Horizon 4
#> MAE   40.956667 37.205000 42.730417 51.424167
#> MAPE   4.118140  5.030128  7.468098  7.434864
#> sMAPE  4.261377  5.208551  7.814417  7.732037
#> RMSE  58.920772 47.660996 49.733955 55.023038
r$global
#>       MAE      MAPE     sMAPE      RMSE 
#> 43.079063  6.012808  6.254095 52.834690

To assess model forecast accuracy, use the efa() function on the fitted model. In this example, the forecasting method is a \(k\)-nearest neighbors algorithm (method = "knn") with autoregressive lags 1 to 4, applied to the UKgas time series. The function estimates forecast accuracy across a 4-step-ahead horizon (h = 4) using multiple evaluation metrics.

The efa() function returns a list with two components:

• per_horizon: Contains accuracy estimates for each individual forecasting step.
• global: Contains the overall accuracy estimates, calculated as the row means of the per_horizon matrix.

Currently, the following forecast accuracy measures are computed:

• MAE: mean absolute error
• MAPE: mean absolute percentage error
• sMAPE: symmetric MAPE
• RMSE: root mean squared error

Next, we describe how the forecasting accuracy measures are computed for a forecasting horizon \(h\) (\(y_t\) and \(\hat{y}_t\) are the actual future value and its forecast for horizon \(t\) respectively):

\[ MAE = \frac{1}{h}\sum_{t=1}^{h} |y_t-\hat{y}_t| \]

\[ MAPE = \frac{1}{h}\sum_{t=1}^{h} 100\frac{|y_t-\hat{y}_t|}{y_t} \] \[ sMAPE = \frac{1}{h}\sum_{t=1}^{h} 200\frac{\left|y_{t}-\hat{y}_{t}\right|}{|y_t|+|\hat{y}_t|} \]

\[ RMSE = \sqrt{\frac{1}{h}\sum_{t=1}^{h} (y_t-\hat{y}_t)^2} \]

Forecast accuracy is estimated using the well-known rolling origin evaluation strategy. This technique generates multiple training and test sets from the time series, as illustrated below for a 4-step-ahead forecasting horizon (h = 4):

In this setup, the last nine observations of the time series are reserved to construct the test sets. Six models are fitted using the expanding training sets, and their forecasts are evaluated against the corresponding test values. Consequently, the forecast error for each step of the horizon is averaged across six evaluation folds.

The number of trailing observations used to form the test sets can be controlled via the size or prop parameters:

m <- create_model(UKgas, lags = 1:4, method = "knn")
# Use the last 9 observations of the series to build test sets
r1 <- efa(m, h = 4, type = "normal", size = 9)
# Use the last 20% observations of the series to build test sets
r2 <- efa(m, h = 4, type = "normal", prop = 0.2)

If the series is short, the training sets may be too small to fit a meaningful model. In such cases, an alternative rolling origin scheme can be applied by setting the type parameter to "minimum":

m <- create_model(UKgas, lags = 1:4, method = "knn")
r <- efa(m, h = 4, type = "minimum")
r$per_horizon
#>       Horizon 1 Horizon 2 Horizon 3 Horizon 4
#> MAE   47.893750 44.934722 44.918229 26.961198
#> MAPE   5.745458  6.757300  8.250563  3.444200
#> sMAPE  5.761843  6.791762  8.329877  3.385892
#> RMSE  56.558584 52.858868 46.912818 26.961198
r$global
#>       MAE      MAPE     sMAPE      RMSE 
#> 41.176975  6.049380  6.067343 45.822867

When this option is selected, the last h observations (4 in this example) are used to construct the test sets, as illustrated below for h = 4:

In this case, the estimated forecast error for horizon 1 is based on 4 errors, for horizon 2 on 3 errors, for horizon 3 on 2 errors, and for horizon 4 on 1 error.

The efa() function also returns the test sets (and their associated predictions) used during accuracy assessment:

timeS <- ts(1:25)
m <- create_model(timeS, lags = 1:3, method = "mt")
r <- efa(m, h = 5, size = 7)
r$test_sets
#>      h=1 h=2 h=3 h=4 h=5
#> [1,]  19  20  21  22  23
#> [2,]  20  21  22  23  24
#> [3,]  21  22  23  24  25
r$predictions
#>      h=1 h=2 h=3 h=4 h=5
#> [1,]  19  20  21  22  23
#> [2,]  20  21  22  23  24
#> [3,]  21  22  23  24  25

Each row in these tables represents a test set or a set of predictions. Below are the test sets generated when using the "minimum" evaluation scheme:

timeS <- ts(1:25)
m <- create_model(timeS, lags = 1:3, method = "mt")
r <- efa(m, h = 3, type = "minimum")
r$test_sets
#>      [,1] [,2] [,3]
#> [1,]   25   NA   NA
#> [2,]   24   25   NA
#> [3,]   23   24   25
r$predictions
#>      [,1] [,2] [,3]
#> [1,]   25   NA   NA
#> [2,]   24   25   NA
#> [3,]   23   24   25

Parameter tuning

Another key feature of utsf is parameter tuning. The tune_grid() function estimates the forecast accuracy of a model across various parameter combinations. The optimal configuration (the one achieving the best forecast accuracy) is then used to train a final model on the full historical series and generate future forecasts. Below is an example:

m <- create_model(UKgas, lags = 1:4, method = "knn")
r <- tune_grid(m, h = 4, tuneGrid = expand.grid(k = 1:7), type = "normal", size = 8)
# To check the estimated forecast accuracy with the different configurations
r$tuneGrid
#>   k      MAE     MAPE    sMAPE     RMSE
#> 1 1 27.50551 4.195273 4.344649 38.72080
#> 2 2 41.67751 5.782108 5.995799 50.67769
#> 3 3 43.07906 6.012808 6.254095 52.83469
#> 4 4 43.58916 6.206692 6.382526 55.14836
#> 5 5 46.03185 6.332209 6.503322 58.41501
#> 6 6 47.64410 6.441446 6.583134 62.46201
#> 7 7 48.95917 6.841679 6.886563 63.75997

In this example, the tuneGrid parameter specifies (as a data frame) the set of hyperparameters to evaluate. The forecast accuracy for each parameter combination is estimated as described in the previous section, using the most recent observations of the time series as a validation set. The size parameter controls the length of this test set.

The tuneGrid component of the list returned by tune_grid() contains these estimation results. Here, the \(k\)-nearest neighbors model with \(k = 1\) achieves the best performance across all forecast accuracy metrics, which is driven by the increasing variance of the series. The optimal parameter combination according to RMSE is then used to generate the final time series forecast:

r$best
#> $k
#> [1] 1
r$forecast
#>           Qtr1      Qtr2      Qtr3      Qtr4
#> 1987 1217.9250  661.4063  388.1828  817.3785

Below, we plot the \(k\) values against their estimated forecast accuracy in terms of RMSE:

plot(r$tuneGrid$k, r$tuneGrid$RMSE, type = "o", pch = 19, xlab = "k (number of nearest neighbors)", ylab = "RMSE", main = "Estimated accuracy")

The following example illustrates tuning a model with multiple hyperparameters. Here, we evaluate a random forest model using three tuning parameters:

m <- create_model(UKDriverDeaths, lags = 1:12, method = "rf")
r <- tune_grid(m, 
               h = 12, 
               tuneGrid = expand.grid(num.trees = c(200, 500), replace = c(TRUE, FALSE), mtry = c(4, 8)),
               type = "normal", 
               size = 12
)
r$tuneGrid
#>   num.trees replace mtry      MAE     MAPE    sMAPE     RMSE
#> 1       200    TRUE    4 103.4673 7.003647 7.386769 103.4673
#> 2       500    TRUE    4 105.3515 7.058263 7.479477 105.3515
#> 3       200   FALSE    4 110.0918 7.297095 7.771264 110.0918
#> 4       500   FALSE    4 107.0101 7.159315 7.578507 107.0101
#> 5       200    TRUE    8 105.4484 7.173581 7.590381 105.4484
#> 6       500    TRUE    8 107.4666 7.316988 7.739419 107.4666
#> 7       200   FALSE    8 107.9476 7.327047 7.740572 107.9476
#> 8       500   FALSE    8 106.6378 7.219095 7.632296 106.6378

Preprocessings and transformations

Model forecast accuracy can often be improved by applying transformations or preprocessing techniques to the target time series. Currently, utsf focuses on preprocessing methods designed for trended time series. This is because most models currently integrated into the package (such as regression trees and \(k\)-nearest neighbors) predict averages of the target values in the training set (i.e., historical observations). When a series exhibits a trend, these historical averages typically fall outside the range of future values.

The following example demonstrates forecasting a trended series (the annual revenue passenger miles flown by commercial airlines in the United States) using a random forest model:

m <- create_model(airmiles, lags = 1:4, method = "rf", trend = "none")
f <- forecast(m, h = 4)
autoplot(f)

In this case, no preprocessing is applied to handle the trend (trend = "none"). As shown below, the forecast fails to capture the upward trajectory of the series because regression tree models can only predict average values of the training targets.

The following subsections explain how to handle trended series using three different transformations:

m <- create_model(airmiles, lags = 1:4, method = "rf", trend = "differences", nfd = 1)
f <- forecast(m, h = 4)
autoplot(f)

# The order of first differences is estimated using the ndiffs function
m <- create_model(airmiles, lags = 1:4, method = "rf", trend = "differences", nfd = -1)
m$differences
#> Order of first differences (selected automatically): 2

timeS <- ts(c(1, 3, 7, 9, 10, 12))
m <- create_model(timeS, lags = 1:2, trend = "none")
cbind(m$features, Targets = m$targets)
#>   Lag2 Lag1 Targets
#> 1    1    3       7
#> 2    3    7       9
#> 3    7    9      10
#> 4    9   10      12

timeS <- ts(c(1, 3, 7, 9, 10, 12))
m <- create_model(timeS, lags = 1:2, trend = "additive", transform_features = FALSE)
cbind(m$features, Targets = m$targets)
#>   Lag2 Lag1 Targets
#> 1    1    3     5.0
#> 2    3    7     4.0
#> 3    7    9     2.0
#> 4    9   10     2.5

timeS <- ts(c(1, 3, 7, 9, 10, 12))
m <- create_model(timeS, lags = 1:2, trend = "additive", transform_features = TRUE)
cbind(m$features, Targets = m$targets)
#>   Lag2 Lag1 Targets
#> 1 -1.0  1.0     5.0
#> 2 -2.0  2.0     4.0
#> 3 -1.0  1.0     2.0
#> 4 -0.5  0.5     2.5

m <- create_model(airmiles, lags = 1:4, method = "rf")
f <- forecast(m, h = 4)
autoplot(f)

t <- ts(10 * 1.05^(1:20))
m_m <- create_model(t, lags = 1:3, method = "rf", trend = "multiplicative")
f_m <- forecast(m_m, h = 4)
m_a <- create_model(t, lags = 1:3, method = "rf", trend = "additive")
f_a <- forecast(m_a, h = 4)
library(vctsfr)
plot_predictions(t, predictions = list(Multiplicative = f_m$pred, Additive = f_a$pred))

m_m <- create_model(AirPassengers, lags = 1:12, method = "knn", trend = "multiplicative")
f_m <- forecast(m_m, h = 36)
m_a <- create_model(AirPassengers, lags = 1:12, method = "knn", trend = "additive")
f_a <- forecast(m_a, h = 36)
library(vctsfr)
plot_predictions(AirPassengers, predictions = list(Multiplicative = f_m$pred, Additive = f_a$pred))

Default parameters

The create_model() function requires only one mandatory parameter: the time series used to train the model. Below, we describe the default settings for the remaining parameters:

• lags: An integer vector defining the autoregressive lags. If frequency(ts) == f (where ts is the input time series and \(f > 1\)), the default lags are 1:f (e.g., 1:4 for quarterly data and 1:12 for monthly data), capturing potential seasonal patterns. If frequency(ts) == 1:
  • Significant lags are automatically selected based on the partial autocorrelation function (stats::pacf()).
  • If no lag is statistically significant, lags 1:5 are selected.
  • If only a single lag is significant and the additive or multiplicative transformation is applied to both targets and features, lags 1:5 are used instead. This prevents applying feature-level transformations when only one lag is present, as scaling a single feature by its own mean is uninformative.
• method: Defaults to the \(k\)-nearest neighbors algorithm.
• trend: Defaults to the additive transformation ("additive"). This method has proven effective for trended time series and generally improves performance across most series types.
• transform_features: Defaults to TRUE, ensuring that both features and targets are transformed when using the additive or multiplicative transformation.

More on using your custom models

In previous sections we have explained how to leverage the package’s features to apply your own models for time series forecasting. In this section we describe how to take advantage of other package’s features. We will employ the random forest model implemented in the randomForest package used previously. First, let us examine the computation of prediction intervals:

library(randomForest)
myModel <- function(X, y, ...) {
  param <- list(...)
  args <- list(x = X, y = y, ntree = 200)
  args <- args[!(names(args) %in% names(param))]
  args <- c(args, param)
  do.call(randomForest::randomForest, args = args)
}
m <- create_model(USAccDeaths, method = myModel)
f <- forecast(m, h = 12, level = 90, PI = TRUE)
f
#>          Point Forecast    Lo 90     Hi 90
#> Jan 1979       8140.180 7616.195  8714.630
#> Feb 1979       7293.334 6860.925  7963.039
#> Mar 1979       8056.349 7475.799  8598.615
#> Apr 1979       8280.013 7752.210  8867.134
#> May 1979       9080.205 8481.096  9648.541
#> Jun 1979       9466.613 8784.123 10066.069
#> Jul 1979      10178.016 9331.077 10673.154
#> Aug 1979       9599.156 8840.990 10107.774
#> Sep 1979       9075.167 8472.870  9703.990
#> Oct 1979       8974.975 8306.085  9549.678
#> Nov 1979       8656.216 8103.380  9277.757
#> Dec 1979       9030.775 8161.252  9449.128

As can be seen, all these features are used in the same way as with the built-in models. For example, let us estimate the forecast accuracy of the previous model for a forecasting horizon of 5:

efa(m, h = 5, type = "normal", size = 8)$per_horizon
#>        Horizon 1  Horizon 2  Horizon 3  Horizon 4  Horizon 5
#> MAE   336.606849 512.539793 587.175618 602.430610 618.444333
#> MAPE    3.433026   5.342922   6.146178   6.509639   6.824317
#> sMAPE   3.505668   5.515537   6.363258   6.753901   7.092678
#> RMSE  371.407984 547.846512 613.880093 639.427652 653.309994

By default, the additive transformation is applied. Let us select the multiplicative transformation to estimate the forecast accuracy under this setting:

m <- create_model(USAccDeaths, lags = 1:12, method = myModel,
                        trend = "multiplicative")
efa(m, h = 5, type = "normal", size = 8)$per_horizon
#>        Horizon 1  Horizon 2  Horizon 3  Horizon 4  Horizon 5
#> MAE   332.535554 488.467860 582.657426 602.744565 590.058133
#> MAPE    3.390765   5.085486   6.089902   6.514085   6.510790
#> sMAPE   3.472280   5.248379   6.307826   6.762925   6.757952
#> RMSE  391.876819 533.569133 617.532382 644.986167 627.151852

Finally, we perform parameter tuning over different values for the mtry and ntree arguments:

m <- create_model(USAccDeaths, method = myModel)
r <- tune_grid(m, h = 6,
               tuneGrid = expand.grid(mtry = 3:6, ntree = c(200, 500)),
               type = "normal"
)
r$tuneGrid
#>   mtry ntree      MAE     MAPE    sMAPE     RMSE
#> 1    3   200 626.2959 6.694872 6.952252 626.2959
#> 2    4   200 611.8487 6.549357 6.796386 611.8487
#> 3    5   200 541.4217 5.811333 6.010729 541.4217
#> 4    6   200 538.5158 5.776617 5.972405 538.5158
#> 5    3   500 616.4211 6.589849 6.838681 616.4211
#> 6    4   500 585.6115 6.269034 6.499482 585.6115
#> 7    5   500 574.3956 6.164767 6.391123 574.3956
#> 8    6   500 572.5947 6.142637 6.366733 572.5947