| Title: | Time Series Forecasting with uncertainty quantification |
|---|---|
| Description: | Univariate and multivariate time series forecasting with uncertainty quantification. |
| Authors: | T. Moudiki |
| Maintainer: | T. Moudiki <[email protected]> |
| License: | BSD_3_clause + file LICENSE |
| Version: | 0.37.2 |
| Built: | 2026-05-23 10:40:18 UTC |
| Source: | https://github.com/Techtonique/ahead |
ANY MODEL+GARCH(1, 1) forecasting
agnosticgarchf( y, h = 5, level = 95, FUN = forecast::auto.arima, seed = 123, ... )agnosticgarchf( y, h = 5, level = 95, FUN = forecast::auto.arima, seed = 123, ... )
y |
a univariate time series |
h |
number of periods for forecasting |
level |
confidence level for prediction intervals |
FUN |
forecasting function for the main model;
default |
seed |
reproducibility seed |
... |
Additional arguments passed to ... |
An object of class "forecast"; a list containing the following elements:
model |
A list containing information about the fitted model |
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
x |
The original time series |
sims |
Simulations of ANYMODEL+GARCH(1, 1) |
T. Moudiki
y <- datasets::EuStockMarkets[ , "DAX"] # require(forecast) # z <- ahead::agnosticgarchf(y=y, h=20) # plot(z)y <- datasets::EuStockMarkets[ , "DAX"] # require(forecast) # z <- ahead::agnosticgarchf(y=y, h=20) # plot(z)
ARMA(1, 1)-GARCH(1, 1) forecasting (with simulation)
armagarchf( y, h = 5, level = 95, B = NULL, cl = 1L, dist = c("student", "gaussian"), seed = 123, ... )armagarchf( y, h = 5, level = 95, B = NULL, cl = 1L, dist = c("student", "gaussian"), seed = 123, ... )
y |
a univariate time series |
h |
number of periods for forecasting |
level |
confidence level for prediction intervals |
B |
number of simulations for |
cl |
an integer; the number of clusters for parallel execution |
dist |
distribution of innovations ("student" or "gaussian") |
seed |
reproducibility seed |
An object of class "forecast"; a list containing the following elements:
model |
A list containing information about the fitted model |
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
x |
The original time series |
sims |
Simulations of ARMA(1, 1)-GARCH(1, 1) |
T. Moudiki
y <- datasets::EuStockMarkets[ , "DAX"] # require(forecast) # z <- ahead::armagarchf(y=y, h=200) # plot(z)y <- datasets::EuStockMarkets[ , "DAX"] # require(forecast) # z <- ahead::armagarchf(y=y, h=200) # plot(z)
Basic forecasting functions for multivariate time series
basicf( y, h = 5, level = 95, method = c("mean", "median", "rw"), type_pi = c("gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap"), block_length = NULL, seed = 1, B = 100, show_progress = TRUE )basicf( y, h = 5, level = 95, method = c("mean", "median", "rw"), type_pi = c("gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap"), block_length = NULL, seed = 1, B = 100, show_progress = TRUE )
y |
A multivariate time series of class |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
method |
forecasting method, either "mean", "median", or random walk ("rw") |
type_pi |
type of prediction interval currently, "gaussian", "bootstrap", "blockbootstrap" or "movingblockbootstrap" |
block_length |
length of block for (circular) "blockbootstrap" or "movingblockbootstrap" |
seed |
reproducibility seed for |
B |
Number of bootstrap replications for |
show_progress |
A boolean; show progress bar for bootstrapping? Default is TRUE. |
An object of class "mtsforecast"; a list containing the following elements:
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
sims |
Model simulations for bootstrapping (basic, or block) |
x |
The original time series |
residuals |
Residuals from the fitted model |
coefficients |
Regression coefficients for |
require(fpp2) res <- ahead::basicf(fpp2::insurance, h=10) par(mfrow=c(1, 2)) plot(res, "TV.advert") plot(res, "Quotes") res <- ahead::basicf(fpp2::insurance, method="rw", h=10) par(mfrow=c(1, 2)) plot(res, "TV.advert") plot(res, "Quotes") # block bootstrap res3 <- ahead::basicf(fpp2::insurance, h=10, type_pi = "bootstrap", B=10) res5 <- ahead::basicf(fpp2::insurance, h=10, type_pi = "blockbootstrap", B=10, block_length = 4) print(res3$sims[[2]]) print(res5$sims[[2]]) par(mfrow=c(2, 2)) plot(res3, "Quotes") plot(res3, "TV.advert") plot(res5, "Quotes") plot(res5, "TV.advert") # moving block bootstrap res6 <- ahead::basicf(fpp2::insurance, h=10, type_pi = "movingblockbootstrap", B=10, block_length = 4, method = "rw") par(mfrow=c(1, 2)) plot(res6, "Quotes") plot(res6, "TV.advert")require(fpp2) res <- ahead::basicf(fpp2::insurance, h=10) par(mfrow=c(1, 2)) plot(res, "TV.advert") plot(res, "Quotes") res <- ahead::basicf(fpp2::insurance, method="rw", h=10) par(mfrow=c(1, 2)) plot(res, "TV.advert") plot(res, "Quotes") # block bootstrap res3 <- ahead::basicf(fpp2::insurance, h=10, type_pi = "bootstrap", B=10) res5 <- ahead::basicf(fpp2::insurance, h=10, type_pi = "blockbootstrap", B=10, block_length = 4) print(res3$sims[[2]]) print(res5$sims[[2]]) par(mfrow=c(2, 2)) plot(res3, "Quotes") plot(res3, "TV.advert") plot(res5, "Quotes") plot(res5, "TV.advert") # moving block bootstrap res6 <- ahead::basicf(fpp2::insurance, h=10, type_pi = "movingblockbootstrap", B=10, block_length = 4, method = "rw") par(mfrow=c(1, 2)) plot(res6, "Quotes") plot(res6, "TV.advert")
Computes forecast combination weights using GLMNET Regression (OLS) regression.
comb_GLMNET(x, custom_error = NULL)comb_GLMNET(x, custom_error = NULL)
x |
An object of class 'foreccomb'. Contains training set (actual values + matrix of model forecasts) and optionally a test set. |
The function integrates the GLMNET Regression forecast combination implementation of the ForecastCombinations package into ForecastComb.
The results are stored in an object of class 'ForecastComb::foreccomb_res', for which separate plot and summary functions are provided.
Returns an object of class ForecastComb::foreccomb_res with the following components:
Method |
Returns the best-fit forecast combination method. |
Models |
Returns the individual input models that were used for the forecast combinations. |
Weights |
Returns the combination weights obtained by applying the combination method to the training set. |
Intercept |
Returns the intercept of the linear regression. |
Fitted |
Returns the fitted values of the combination method for the training set. |
Accuracy_Train |
Returns range of summary measures of the forecast accuracy for the training set. |
Forecasts_Test |
Returns forecasts produced by the combination method for the test set. Only returned if input included a forecast matrix for the test set. |
Accuracy_Test |
Returns range of summary measures of the forecast accuracy for the test set. Only returned if input included a forecast matrix and a vector of actual values for the test set. |
Input_Data |
Returns the data forwarded to the method. |
Forecast_comb,
foreccomb,
plot.ForecastComb::foreccomb_res,
summary.ForecastComb::foreccomb_res,
accuracy
library(ForecastComb) data(electricity) print(head(electricity)) forecasting_methods <- colnames(electricity)[1:5] train_obs <- electricity[1:84, "Actual"] train_pred <- electricity[1:84, forecasting_methods] test_obs <- electricity[85:123, "Actual"] test_pred <- electricity[85:123, forecasting_methods] data <- ForecastComb::foreccomb(train_obs, train_pred, test_obs, test_pred) # obj <- ahead::comb_GLMNET(data))library(ForecastComb) data(electricity) print(head(electricity)) forecasting_methods <- colnames(electricity)[1:5] train_obs <- electricity[1:84, "Actual"] train_pred <- electricity[1:84, forecasting_methods] test_obs <- electricity[85:123, "Actual"] test_pred <- electricity[85:123, forecasting_methods] data <- ForecastComb::foreccomb(train_obs, train_pred, test_obs, test_pred) # obj <- ahead::comb_GLMNET(data))
Computes forecast combination weights using ordinary least squares (OLS) regression.
comb_OLS(x, custom_error = NULL)comb_OLS(x, custom_error = NULL)
x |
An object of class 'foreccomb'. Contains training set (actual values + matrix of model forecasts) and optionally a test set. |
The function integrates the ordinary least squares (OLS) forecast combination implementation of the ForecastCombinations package into ForecastComb.
The OLS combination method (Granger and Ramanathan (1984)) uses ordinary least squares to
estimate the weights, , as well as an intercept, , for the combination of
the forecasts.
Suppose that there are not perfectly collinear predictors ,
then the forecast combination for one data point can be represented as:
An appealing feature of the method is its bias correction through the intercept – even if one or more of the individual
predictors are biased, the resulting combined forecast is unbiased. A disadvantage of the method is that it places no
restriction on the combination weights (i.e., they do not add up to 1 and can be negative), which can make interpretation
hard. Another issue, documented in Nowotarski et al. (2014), is the method's unstable behavior
when predictors are highly correlated (which is the norm in forecast combination): Minor fluctuations in the sample
can cause major shifts of the coefficient vector (‘bouncing betas’) – often causing poor out-of-sample performance.
This issue is addressed by the comb_LAD method that is more robust to outliers.
The results are stored in an object of class 'ForecastComb::foreccomb_res', for which separate plot and summary functions are provided.
Returns an object of class ForecastComb::foreccomb_res with the following components:
Method |
Returns the best-fit forecast combination method. |
Models |
Returns the individual input models that were used for the forecast combinations. |
Weights |
Returns the combination weights obtained by applying the combination method to the training set. |
Intercept |
Returns the intercept of the linear regression. |
Fitted |
Returns the fitted values of the combination method for the training set. |
Accuracy_Train |
Returns range of summary measures of the forecast accuracy for the training set. |
Forecasts_Test |
Returns forecasts produced by the combination method for the test set. Only returned if input included a forecast matrix for the test set. |
Accuracy_Test |
Returns range of summary measures of the forecast accuracy for the test set. Only returned if input included a forecast matrix and a vector of actual values for the test set. |
Input_Data |
Returns the data forwarded to the method. |
Forecast_comb,
foreccomb,
plot.ForecastComb::foreccomb_res,
summary.ForecastComb::foreccomb_res,
accuracy
obs <- rnorm(100) preds <- matrix(rnorm(1000, 1), 100, 10) train_o<-obs[1:80] train_p<-preds[1:80,] test_o<-obs[81:100] test_p<-preds[81:100,] data<-ForecastComb::foreccomb(train_o, train_p, test_o, test_p) ahead::comb_OLS(data)obs <- rnorm(100) preds <- matrix(rnorm(1000, 1), 100, 10) train_o<-obs[1:80] train_p<-preds[1:80,] test_o<-obs[81:100] test_p<-preds[81:100,] data<-ForecastComb::foreccomb(train_o, train_p, test_o, test_p) ahead::comb_OLS(data)
Computes forecast combination weights using Ridge Regression (OLS) regression.
comb_Ridge(x, custom_error = NULL)comb_Ridge(x, custom_error = NULL)
x |
An object of class 'foreccomb'. Contains training set (actual values + matrix of model forecasts) and optionally a test set. |
The function integrates the Ridge Regression forecast combination implementation of the ForecastCombinations package into ForecastComb.
The results are stored in an object of class 'ForecastComb::foreccomb_res', for which separate plot and summary functions are provided.
Returns an object of class ForecastComb::foreccomb_res with the following components:
Method |
Returns the best-fit forecast combination method. |
Models |
Returns the individual input models that were used for the forecast combinations. |
Weights |
Returns the combination weights obtained by applying the combination method to the training set. |
Intercept |
Returns the intercept of the linear regression. |
Fitted |
Returns the fitted values of the combination method for the training set. |
Accuracy_Train |
Returns range of summary measures of the forecast accuracy for the training set. |
Forecasts_Test |
Returns forecasts produced by the combination method for the test set. Only returned if input included a forecast matrix for the test set. |
Accuracy_Test |
Returns range of summary measures of the forecast accuracy for the test set. Only returned if input included a forecast matrix and a vector of actual values for the test set. |
Input_Data |
Returns the data forwarded to the method. |
Forecast_comb,
foreccomb,
plot.ForecastComb::foreccomb_res,
summary.ForecastComb::foreccomb_res,
accuracy
obs <- rnorm(100) preds <- matrix(rnorm(1000, 1), 100, 10) train_o<-obs[1:80] train_p<-preds[1:80,] test_o<-obs[81:100] test_p<-preds[81:100,] data<-ForecastComb::foreccomb(train_o, train_p, test_o, test_p) ahead::comb_Ridge(data)obs <- rnorm(100) preds <- matrix(rnorm(1000, 1), 100, 10) train_o<-obs[1:80] train_p<-preds[1:80,] test_o<-obs[81:100] test_p<-preds[81:100,] data<-ForecastComb::foreccomb(train_o, train_p, test_o, test_p) ahead::comb_Ridge(data)
Compute global attention weights and context vectors for time series
computeattention( series, attention_type = c("dot_product", "scaled_dot_product", "cosine", "exponential", "gaussian", "linear", "value_based", "hybrid", "parametric"), window_size = 3, decay_factor = 5, temperature = 1, sigma = 1, sensitivity = 1, alpha = 0.5, beta = 0.5 )computeattention( series, attention_type = c("dot_product", "scaled_dot_product", "cosine", "exponential", "gaussian", "linear", "value_based", "hybrid", "parametric"), window_size = 3, decay_factor = 5, temperature = 1, sigma = 1, sensitivity = 1, alpha = 0.5, beta = 0.5 )
series |
Numeric vector containing the time series of length n |
attention_type |
String specifying the type of attention mechanism to use. Options are: "cosine", "exponential", "dot_product", "scaled_dot_product", "gaussian", "linear", "value_based", "hybrid", "parametric". Default is "cosine". |
window_size |
Integer parameter for window size (applicable for "cosine" attention). |
decay_factor |
Double for decay factor (applicable for "exponential" attention). |
temperature |
Double for temperature (applicable for "scaled_dot_product" attention). |
sigma |
Double for sigma (applicable for "gaussian" attention). |
sensitivity |
Double for sensitivity (applicable for "value_based" or "hybrid" attention). |
alpha |
Double for alpha (applicable for "parametric" attention). |
beta |
Double for beta (applicable for "parametric" attention). |
List containing:
attention_weights |
n × n matrix where entry (i,j) represents the attention weight of time j on time i. Only entries j <= i are non-zero (causal attention). |
context_vectors |
Vector of length n where each entry i is the weighted sum of all values up to time i, using the attention weights. |
# For a series of length 5 using "cosine" attention series <- c(1, 2, 3, 4, 5) result <- computeattention(series, attention_type = "cosine", window_size = 3) # attention_weights will be 5x5 matrix # context_vectors will be length 5 dim(result$attention_weights) # [1] 5 5 length(result$context_vectors) # [1] 5# For a series of length 5 using "cosine" attention series <- c(1, 2, 3, 4, 5) result <- computeattention(series, attention_type = "cosine", window_size = 3) # attention_weights will be 5x5 matrix # context_vectors will be length 5 dim(result$attention_weights) # [1] 5 5 length(result$context_vectors) # [1] 5
Combines a mean model and a volatility model to produce probabilistic forecasts. Uncertainty is propagated from two sources: simulated innovation draws and simulated volatility paths derived from the variance model's residuals via bootstrap.
condvolf( y, h = 10L, mean_model = forecast::auto.arima, sigma_model = forecast::auto.arima, innovation = c("empirical", "gaussian", "student"), level = 95, B = 1000L, bootstrap_vol = TRUE, ... )condvolf( y, h = 10L, mean_model = forecast::auto.arima, sigma_model = forecast::auto.arima, innovation = c("empirical", "gaussian", "student"), level = 95, B = 1000L, bootstrap_vol = TRUE, ... )
y |
A numeric vector or time series of class |
h |
Forecasting horizon. |
mean_model |
Function to fit the mean model (default: |
sigma_model |
Function to fit the variance model on squared residuals
(default: |
innovation |
Distribution for standardized innovations. One of
|
level |
Confidence level for prediction intervals (default: 95). |
B |
Number of simulated paths (default: 2000). |
bootstrap_vol |
Whether to bootstrap volatility residuals (default: TRUE). If FALSE, uses parametric normal approximation for volatility. |
... |
Additional arguments passed to |
A forecast object with components:
xOriginal time series.
meanPoint forecast (mean of simulations).
lowerLower prediction interval bound.
upperUpper prediction interval bound.
levelConfidence level.
simsMatrix of simulated forecast paths (h x B).
modelList containing mean and sigma model fits.
residualsResiduals from the mean model.
standardized_residualsScaled standardized residuals used for innovation fitting.
methodCharacter string describing the method used.
library(forecast) # Basic usage with Google stock data y <- fpp2::goog200 # Gaussian innovations (fastest, often sufficient) fc1 <- condvolf(y, h = 20, innovation = "gaussian") plot(fc1) # Student-t innovations (heavy-tailed) fc2 <- condvolf(y, h = 20, innovation = "student") plot(fc2) # Empirical innovations (non-parametric) fc3 <- condvolf(y, h = 20, innovation = "empirical") plot(fc3) # Using different mean and volatility models fc4 <- condvolf(y, h = 20, mean_model = forecast::thetaf, sigma_model = forecast::ets, innovation = "gaussian") # Compare prediction intervals par(mfrow = c(2, 2)) plot(fc1, main = "Gaussian innovations") plot(fc2, main = "Student-t innovations") plot(fc3, main = "Empirical innovations") plot(fc4, main = "Different mean and volatility models") par(mfrow = c(1, 1))library(forecast) # Basic usage with Google stock data y <- fpp2::goog200 # Gaussian innovations (fastest, often sufficient) fc1 <- condvolf(y, h = 20, innovation = "gaussian") plot(fc1) # Student-t innovations (heavy-tailed) fc2 <- condvolf(y, h = 20, innovation = "student") plot(fc2) # Empirical innovations (non-parametric) fc3 <- condvolf(y, h = 20, innovation = "empirical") plot(fc3) # Using different mean and volatility models fc4 <- condvolf(y, h = 20, mean_model = forecast::thetaf, sigma_model = forecast::ets, innovation = "gaussian") # Compare prediction intervals par(mfrow = c(2, 2)) plot(fc1, main = "Gaussian innovations") plot(fc2, main = "Student-t innovations") plot(fc3, main = "Empirical innovations") plot(fc4, main = "Different mean and volatility models") par(mfrow = c(1, 1))
This function allows to conformalize any forecasting function.
conformalize( FUN, y, h, level = 95, method = c("block-bootstrap", "surrogate", "kde", "bootstrap", "fitdistr", "meboot"), nsim = 100L, block_size = 5, seed = 123L, B = NULL, ... )conformalize( FUN, y, h, level = 95, method = c("block-bootstrap", "surrogate", "kde", "bootstrap", "fitdistr", "meboot"), nsim = 100L, block_size = 5, seed = 123L, B = NULL, ... )
FUN |
A forecasting function. |
y |
A time series ( |
h |
Forecasting horizon. |
level |
Confidence level. |
method |
Method to be used for conformalization (simulation of calibrated residuals). |
nsim |
Number of simulations. |
block_size |
Block size for block-bootstrap. |
seed |
Seed for reproducibility. |
B |
Alias for |
... |
Additional arguments to be passed to the forecasting function. |
An object of class forecast.
y <- fdeaths h <- 25L obj <- conformalize(FUN=forecast::ets, y, h); plot(obj) obj <- conformalize(FUN=HoltWinters, y=y, h=h, seasonal = "mult"); plot(obj)y <- fdeaths h <- 25L obj <- conformalize(FUN=forecast::ets, y, h); plot(obj) obj <- conformalize(FUN=HoltWinters, y=y, h=h, seasonal = "mult"); plot(obj)
This function performs ridge regression forecasting using attention-based context vectors as external regressors. Context vectors are computed from the training data using various attention mechanisms and then used to enhance the forecast model.
contextridge2f( y, h = 5L, attention_type = "exponential", window_size = 3, decay_factor = 5, temperature = 1, sigma = 1, sensitivity = 1, alpha = 0.5, beta = 0.5, ... )contextridge2f( y, h = 5L, attention_type = "exponential", window_size = 3, decay_factor = 5, temperature = 1, sigma = 1, sensitivity = 1, alpha = 0.5, beta = 0.5, ... )
y |
A multivariate time series object. |
h |
Integer. Forecast horizon. Defaults to 5. |
attention_type |
String specifying the type of attention mechanism. Options are: "cosine", "exponential", "dot_product", "scaled_dot_product", "gaussian", "linear", "value_based", "hybrid", "parametric". Default is "exponential". |
window_size |
Integer parameter for window size (applicable for "cosine" attention). Defaults to 3. |
decay_factor |
Double for decay factor (applicable for "exponential" and "hybrid" attention). Defaults to 5.0. |
temperature |
Double for temperature (applicable for "scaled_dot_product" attention). Defaults to 1.0. |
sigma |
Double for sigma (applicable for "gaussian" attention). Defaults to 1.0. |
sensitivity |
Double for sensitivity (applicable for "value_based" or "hybrid" attention). Defaults to 1.0. |
alpha |
Double for alpha (applicable for "parametric" attention). Defaults to 0.5. |
beta |
Double for beta (applicable for "parametric" attention). Defaults to 0.5. |
... |
Additional arguments passed to |
This approach allows the model to leverage temporal dependencies captured by attention mechanisms, potentially improving forecast accuracy by incorporating weighted historical information.
An object returned by ahead::ridge2f for the forecast,
typically a list including mean and prediction intervals.
ahead::computeattention
plot(contextridge2f(AirPassengers, h = 15, lags = 15, attention_type = "exponential")) plot(contextridge2f(fdeaths, h = 20, lags = 15, attention_type = "exponential"))plot(contextridge2f(AirPassengers, h = 15, lags = 15, attention_type = "exponential")) plot(contextridge2f(fdeaths, h = 20, lags = 15, attention_type = "exponential"))
Create simple 2-level hierarchical time series
create_2level_hts(bts)create_2level_hts(bts)
bts |
Bottom-level time series matrix (T x m) |
List with total and bottom series
Create trend and seasonality features for univariate time series
createtrendseason(y)createtrendseason(y)
y |
a univariate time series |
a vector or matrix of features
y <- ts(rnorm(100), start = 1, frequency = 12) createtrendseason(y) createtrendseason(USAccDeaths)y <- ts(rnorm(100), start = 1, frequency = 12) createtrendseason(y) createtrendseason(USAccDeaths)
Returns forecasts and prediction intervals for a context-aware theta method forecast.
ctxthetaf( y, h = if (frequency(y) > 1) 2 * frequency(y) else 10, level = 95, theta = 0.5, fit_func = lm, predict_func = predict, type_pi = c("gaussian", "block-bootstrap", "surrogate", "kde", "bootstrap", "fitdistr", "meboot"), nsim = 100L, block_size = 5, seed = 123L, B = NULL, x = y, ... )ctxthetaf( y, h = if (frequency(y) > 1) 2 * frequency(y) else 10, level = 95, theta = 0.5, fit_func = lm, predict_func = predict, type_pi = c("gaussian", "block-bootstrap", "surrogate", "kde", "bootstrap", "fitdistr", "meboot"), nsim = 100L, block_size = 5, seed = 123L, B = NULL, x = y, ... )
y |
Time series (ts object or numeric vector) |
h |
Forecast horizon (default: 2*frequency for seasonal, 10 for non-seasonal) |
level |
Confidence level(s) for prediction intervals (default: 95) |
theta |
Theta coefficient for drift weighting (default: 0.5, giving theta line = 2)
|
fit_func |
Model fitting function (default: lm) |
predict_func |
Prediction function (default: predict) |
type_pi |
Default is "gaussian", otherwise, uses sequential split conformal prediction as in Moudiki, T. (2025) |
nsim |
Number of simulations when |
block_size |
Size of block when |
seed |
Reproducibility seed |
B |
Alias for |
x |
Deprecated, use y instead |
... |
Additional arguments passed to fit_func |
The classical theta method of Assimakopoulos and Nikolopoulos (2000) is equivalent to simple exponential smoothing with drift, using theta = 2. This implementation extends the method by:
Allowing flexible specification of theta parameter (default 0.5, which gives theta = 2)
Using machine learning models to capture non-linear drift patterns
Computing time-varying slopes via numerical differentiation
An object of class forecast
T. Moudiki (based on RJH)
Hyndman, R.J., and Billah, B. (2003) Unmasking the Theta method. International J. Forecasting, 19, 287-290.
Moudiki, T. (2025). Conformal Predictive Simulations for Univariate Time Series. Proceedings of Machine Learning Research, 266, 1-2.
forecast::thetaf(), forecast::ses()
# Classical theta-like behavior (theta = 0.5) fit1 <- ctxthetaf(AirPassengers) plot(fit1) # No drift (theta = 0) fit2 <- ctxthetaf(AirPassengers, theta = 0) plot(fit2) # Amplified drift (theta = 1) fit3 <- ctxthetaf(AirPassengers, theta = 1) plot(fit3) # With Random Forest fit4 <- ctxthetaf(AirPassengers, fit_func = randomForest::randomForest) plot(fit4)# Classical theta-like behavior (theta = 0.5) fit1 <- ctxthetaf(AirPassengers) plot(fit1) # No drift (theta = 0) fit2 <- ctxthetaf(AirPassengers, theta = 0) plot(fit2) # Amplified drift (theta = 1) fit3 <- ctxthetaf(AirPassengers, theta = 1) plot(fit3) # With Random Forest fit4 <- ctxthetaf(AirPassengers, fit_func = randomForest::randomForest) plot(fit4)
Direct sampling
direct_sampling( data = NULL, n = 100L, method = c("kde", "surrogate", "bootstrap"), kde = NULL, seed = NULL, ... )direct_sampling( data = NULL, n = 100L, method = c("kde", "surrogate", "bootstrap"), kde = NULL, seed = NULL, ... )
data |
A numeric vector or matrix. |
n |
The number of samples to draw. |
method |
The method to use for sampling. |
kde |
The kernel density estimate to use for sampling. |
seed |
The seed to use for sampling. |
... |
Additional arguments to pass to the density function. |
Adapted from forecast::nnetar, with alternative fitting functions (see examples)
dynrmf( y, h = 5, level = 95, fit_func = ahead::ridge, predict_func = predict, fit_params = NULL, type_pi = c("gaussian", "E", "A", "T"), xreg_fit = NULL, xreg_predict = NULL, ... )dynrmf( y, h = 5, level = 95, fit_func = ahead::ridge, predict_func = predict, fit_params = NULL, type_pi = c("gaussian", "E", "A", "T"), xreg_fit = NULL, xreg_predict = NULL, ... )
y |
A numeric vector or time series of class |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
fit_func |
Fitting function (Statistical/ML model). Default is Ridge regression. |
predict_func |
Prediction function (Statistical/ML model) |
fit_params |
a list of additional parameters for the fitting function |
type_pi |
Type of prediction interval (currently "gaussian", ETS: "E", Arima: "A" or Theta: "T") |
xreg_fit |
Optionally, a vector or matrix of external regressors, which must have the same number of rows as y. Must be numeric. |
xreg_predict |
Future values of external regressor variables. |
... |
additional parameters |
a list; an object of class forecast.
The function summary is used to obtain and print a summary of the
results.
The generic accessor functions fitted.values and residuals
extract useful features.
T. Moudiki
Hyndman, R. J., & Athanasopoulos, G. (2018). Forecasting: principles and practice. OTexts.
Hyndman R, Athanasopoulos G, Bergmeir C, Caceres G, Chhay L,
O'Hara-Wild M, Petropoulos F, Razbash S, Wang E, Yasmeen F (2021).
forecast: Forecasting functions for time series and linear models. R
package version 8.14, <URL: https://pkg.robjhyndman.com/forecast/>.
# Example 0: with Ridge regression par(mfrow=c(3, 2)) plot(dynrmf(USAccDeaths, h=20, level=95)) plot(dynrmf(AirPassengers, h=20, level=95)) plot(dynrmf(lynx, h=20, level=95)) plot(dynrmf(WWWusage, h=20, level=95)) plot(dynrmf(Nile, h=20, level=95)) plot(dynrmf(fdeaths, h=20, level=95)) # Example 1: with Random Forest ## Not run: require(randomForest) par(mfrow=c(3, 2)) plot(dynrmf(USAccDeaths, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(AirPassengers, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(lynx, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(WWWusage, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(Nile, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) ## End(Not run) # Example 2: with SVM ## Not run: require(e1071) par(mfrow=c(2, 2)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "linear"), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "polynomial"), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "radial"), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "sigmoid"), predict_func = predict)) ## End(Not run)# Example 0: with Ridge regression par(mfrow=c(3, 2)) plot(dynrmf(USAccDeaths, h=20, level=95)) plot(dynrmf(AirPassengers, h=20, level=95)) plot(dynrmf(lynx, h=20, level=95)) plot(dynrmf(WWWusage, h=20, level=95)) plot(dynrmf(Nile, h=20, level=95)) plot(dynrmf(fdeaths, h=20, level=95)) # Example 1: with Random Forest ## Not run: require(randomForest) par(mfrow=c(3, 2)) plot(dynrmf(USAccDeaths, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(AirPassengers, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(lynx, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(WWWusage, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(Nile, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = randomForest::randomForest, fit_params = list(ntree = 50), predict_func = predict)) ## End(Not run) # Example 2: with SVM ## Not run: require(e1071) par(mfrow=c(2, 2)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "linear"), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "polynomial"), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "radial"), predict_func = predict)) plot(dynrmf(fdeaths, h=20, level=95, fit_func = e1071::svm, fit_params = list(kernel = "sigmoid"), predict_func = predict)) ## End(Not run)
This function computes the sensitivity of a forecasted time series (y)
to small perturbations in an external regressor (xreg) using the dynrmf function.
dynrmf_sensi(y, xreg, h = NULL, level = 99, zero = 1e-04, ...)dynrmf_sensi(y, xreg, h = NULL, level = 99, zero = 1e-04, ...)
y |
A univariate time series object to forecast. |
xreg |
A time series of external regressors (same length as |
h |
Forecast horizon (length of test set). If NULL, uses length of xreg. |
level |
Confidence level for |
zero |
Small positive constant to avoid division by zero (default 1e-4). |
... |
Additional parameters to be passed to |
A list containing:
Sensitivity effects on the mean forecast.
Forecast without xreg.
Forecast with xreg.
Forecast with positive perturbation on xreg.
Forecast with negative perturbation on xreg.
Forecast horizon used.
Test set of xreg used in forecasting.
Perturbation applied to xreg.
List with training and testing sets of y.
sensitivity_results_auto <- ahead::dynrmf_sensi( y = fpp2::uschange[, "Consumption"], xreg = fpp2::uschange[, "Income"], h = 10 ) plot1 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_auto, title = "Sensitivity of Consumption to Income", y_label = "Effect (DeltaConsumption / DeltaIncome)") print(plot1)sensitivity_results_auto <- ahead::dynrmf_sensi( y = fpp2::uschange[, "Consumption"], xreg = fpp2::uschange[, "Income"], h = 10 ) plot1 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_auto, title = "Sensitivity of Consumption to Income", y_label = "Effect (DeltaConsumption / DeltaIncome)") print(plot1)
Compute exact Shapley values for an ahead::dynrmf model
dynrmf_shap( y, xreg_fit, xreg_predict, fit_func = ahead::ridge, verbose = TRUE, tol = 1e-10 )dynrmf_shap( y, xreg_fit, xreg_predict, fit_func = ahead::ridge, verbose = TRUE, tol = 1e-10 )
y |
Training time series (passed to ahead::dynrmf as y) |
xreg_fit |
Regressor matrix used during training |
xreg_predict |
Regressor matrix for the forecast horizon (h x p); h and feature names are derived from this matrix. |
fit_func |
Fitting function for ahead::dynrmf (default: ahead::ridge) |
verbose |
Print progress per horizon and efficiency check (default: TRUE) |
tol |
Tolerance for the efficiency-axiom check (default: 1e-10) |
A list of class "dynrmf_shap" with: $S : h x p matrix of Shapley values (per variable, per horizon) $X : test regressor matrix $baseline : scalar baseline (mean of empty-coalition forecasts) $baselines : per-horizon baseline forecasts $predictions : per-horizon full-coalition forecasts $feature_names $h
Combined ets, arima, and theta (eat) forecasting (uses forecast::ets,
forecast::auto.arima, forecast::thetaf)
eatf( y, h = 5, level = 95, method = c("EAT", "E", "A", "T"), weights = rep(1/3, 3), type_pi = c("gaussian", "E", "A", "T"), ... )eatf( y, h = 5, level = 95, method = c("EAT", "E", "A", "T"), weights = rep(1/3, 3), type_pi = c("gaussian", "E", "A", "T"), ... )
y |
a univariate time series |
h |
number of periods for forecasting |
level |
confidence level for prediction intervals |
method |
forecasting method: "E" for |
weights |
weights for each method, in method |
type_pi |
type of prediction interval: currently ETS: "E", Auto.Arima: "A" or Theta: "T" |
... |
additional parameters to be passed to |
ensemble forecasts obtained from forecast::ets,
forecast::auto.arima and forecast::theta (with weights)
An object of class "forecast"; a list containing the following elements:
model |
A list containing information about the fitted model |
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
x |
The original time series |
residuals |
Residuals from the fitted model |
T. Moudiki
Hyndman R, Athanasopoulos G, Bergmeir C, Caceres G, Chhay L,
O'Hara-Wild M, Petropoulos F, Razbash S, Wang E, Yasmeen F (2021).
forecast: Forecasting functions for time series and linear models. R
package version 8.14, <URL: https://pkg.robjhyndman.com/forecast/>.
Hyndman RJ, Khandakar Y (2008). 'Automatic time series forecasting: the forecast package for R.' Journal of Statistical Software, 26 (3), 1-22. <URL: https://www.jstatsoft.org/article/view/v027i03>.
Assimakopoulos, V. and Nikolopoulos, K. (2000). The theta model: a decomposition approach to forecasting. International Journal of Forecasting 16, 521-530.
Hyndman, R.J., and Billah, B. (2003) Unmasking the Theta method. International J. Forecasting, 19, 287-290.
require(forecast) ## Not run: print(ahead::eatf(WWWusage, method = "EAT", weights = c(0.5, 0, 0.5))) print(ahead::eatf(WWWusage, method = "EAT")) print(ahead::eatf(WWWusage, method = "E")) print(ahead::eatf(WWWusage, method = "A")) print(ahead::eatf(WWWusage, method = "T")) obj <- ahead::eatf(WWWusage, method = "EAT", weights = c(0, 0.5, 0.5), h=10, type_pi = "T") plot(obj) obj <- ahead::eatf(WWWusage, method = "EAT", weights = c(0, 0.5, 0.5), h=10, type_pi="A") plot(obj) ## End(Not run) par(mfrow=c(3, 2)) plot(ahead::eatf(USAccDeaths, h=10, level=95)) plot(ahead::eatf(AirPassengers, h=10, level=95, type_pi = "T")) plot(ahead::eatf(lynx, h=10, level=95, type_pi = "A")) plot(ahead::eatf(WWWusage, h=10, level=95, type_pi = "E")) plot(ahead::eatf(Nile, h=10, level=95)) plot(ahead::eatf(fdeaths, h=10, level=95))require(forecast) ## Not run: print(ahead::eatf(WWWusage, method = "EAT", weights = c(0.5, 0, 0.5))) print(ahead::eatf(WWWusage, method = "EAT")) print(ahead::eatf(WWWusage, method = "E")) print(ahead::eatf(WWWusage, method = "A")) print(ahead::eatf(WWWusage, method = "T")) obj <- ahead::eatf(WWWusage, method = "EAT", weights = c(0, 0.5, 0.5), h=10, type_pi = "T") plot(obj) obj <- ahead::eatf(WWWusage, method = "EAT", weights = c(0, 0.5, 0.5), h=10, type_pi="A") plot(obj) ## End(Not run) par(mfrow=c(3, 2)) plot(ahead::eatf(USAccDeaths, h=10, level=95)) plot(ahead::eatf(AirPassengers, h=10, level=95, type_pi = "T")) plot(ahead::eatf(lynx, h=10, level=95, type_pi = "A")) plot(ahead::eatf(WWWusage, h=10, level=95, type_pi = "E")) plot(ahead::eatf(Nile, h=10, level=95)) plot(ahead::eatf(fdeaths, h=10, level=95))
dynrmf)Fit univariate time series using caret ML model (for use with dynrmf)
fit_func( x, y, method = "ranger", initial_window = 10L, horizon = 10L, fixed_window = FALSE, tune_length = 5, summary_function = NULL, verbose = TRUE )fit_func( x, y, method = "ranger", initial_window = 10L, horizon = 10L, fixed_window = FALSE, tune_length = 5, summary_function = NULL, verbose = TRUE )
x |
A matrix of predictors |
y |
A vector of responses |
method |
The caret method to use for fitting the model |
initial_window |
The initial window size |
horizon |
The forecast horizon |
fixed_window |
Whether to use a fixed window size |
tune_length |
Length of the tuning grid |
verbose |
Whether to print the model summary |
A model object
Fit and forecast for benchmarking purposes
fitforecast( y, h = NULL, pct_train = 0.9, pct_calibration = 0.5, method = c("thetaf", "arima", "ets", "te", "tbats", "tslm", "dynrmf", "ridge2f", "naive", "snaive", "arimagarch"), level = 95, B = 1000L, seed = 17223L, graph = TRUE, conformalize = FALSE, type_calibration = c("splitconformal", "cv1", "loocv"), gap = 3L, agg = c("mean", "median"), vol = c("constant", "garch"), type_sim = c("kde", "surrogate", "bootstrap", "meboot"), ... )fitforecast( y, h = NULL, pct_train = 0.9, pct_calibration = 0.5, method = c("thetaf", "arima", "ets", "te", "tbats", "tslm", "dynrmf", "ridge2f", "naive", "snaive", "arimagarch"), level = 95, B = 1000L, seed = 17223L, graph = TRUE, conformalize = FALSE, type_calibration = c("splitconformal", "cv1", "loocv"), gap = 3L, agg = c("mean", "median"), vol = c("constant", "garch"), type_sim = c("kde", "surrogate", "bootstrap", "meboot"), ... )
y |
A univariate time series of class |
h |
Forecasting horizon (default is |
pct_train |
Percentage of data in the training set, when |
pct_calibration |
Percentage of data in the calibration set for conformal prediction |
method |
For now "thetaf" (default), "arima", "ets", "tbats", "tslm", "dynrmf" (from ahead), "ridge2f" (from ahead), "naive", "snaive" |
level |
Confidence level for prediction intervals in %, default is 95 |
B |
Number of bootstrap replications or number of simulations (yes, 'B' is unfortunate) |
seed |
Reproducibility seed |
graph |
Plot or not? |
conformalize |
Calibrate or not? |
type_calibration |
"splitconformal" (default conformal method), "cv1" (do not use), "loocv" (do not use) |
gap |
length of training set for loocv conformal (do not use) |
agg |
"mean" or "median" (aggregation method) |
vol |
"constant" or "garch" (type of volatility modeling for calibrated residuals) |
type_sim |
"kde", "surrogate", "bootstrap" (type of simulation for calibrated residuals) |
... |
additional parameters |
an object of class 'forecast' with additional information
par(mfrow=c(2, 2)) obj1 <- ahead::fitforecast(AirPassengers) obj2 <- ahead::fitforecast(AirPassengers, conformalize = TRUE) plot(AirPassengers) plot(obj1) obj2$plot() obj2$plot("simulations")par(mfrow=c(2, 2)) obj1 <- ahead::fitforecast(AirPassengers) obj2 <- ahead::fitforecast(AirPassengers, conformalize = TRUE) plot(AirPassengers) plot(obj1) obj2$plot() obj2$plot("simulations")
Generate synthetic time series via model-based residual bootstrap
generate_synthetic_ts( y, model_type = "auto", model_func = NULL, formula = NULL, order = NULL, seasonal = list(order = c(0, 0, 0), period = NA), length_out = NULL, n_sim = 1000, seed = NULL, ... )generate_synthetic_ts( y, model_type = "auto", model_func = NULL, formula = NULL, order = NULL, seasonal = list(order = c(0, 0, 0), period = NA), length_out = NULL, n_sim = 1000, seed = NULL, ... )
y |
Time series data (ts or mts object) |
model_type |
One of "auto", "auto.arima", "tslm", "ets", "arima", "custom" |
model_func |
Custom model function (required if model_type = "custom") |
formula |
Formula for tslm (default: trend + season) |
order |
ARIMA order (p,d,q) for model_type = "arima" |
seasonal |
List with order and period for seasonal ARIMA |
length_out |
Length of synthetic series to generate (default = length of y) |
n_sim |
Number of synthetic series to generate |
seed |
Random seed for reproducibility (optional) |
... |
Extra arguments passed to the modelling function |
A list with synthetic series and model information
This function allows to call any function "of class forecast" in a unified way.
genericforecast(FUN, y, h, level = 95, ...)genericforecast(FUN, y, h, level = 95, ...)
FUN |
A forecasting function. |
y |
A time series ( |
h |
Forecasting horizon. |
level |
The confidence level. |
... |
Additional arguments to be passed to the forecasting function. |
An object of class forecast.
y <- fdeaths h <- 25L plot(ahead::genericforecast(FUN=forecast::thetaf, y, h)) plot(ahead::genericforecast(FUN=ahead::dynrmf, y, h)) plot(ahead::genericforecast(FUN=forecast::tbats, y=y, h=h, use.box.cox = TRUE, use.trend=FALSE)) plot(ahead::genericforecast(FUN=forecast::auto.arima, y, h))y <- fdeaths h <- 25L plot(ahead::genericforecast(FUN=forecast::thetaf, y, h)) plot(ahead::genericforecast(FUN=ahead::dynrmf, y, h)) plot(ahead::genericforecast(FUN=forecast::tbats, y=y, h=h, use.box.cox = TRUE, use.trend=FALSE)) plot(ahead::genericforecast(FUN=forecast::auto.arima, y, h))
This function calculates various error metrics for a given prediction object.
geterror(obj, actual, level = 95)geterror(obj, actual, level = 95)
obj |
A prediction object containing mean predictions and intervals. |
actual |
A numeric vector of actual values. |
level |
The confidence level for the prediction intervals. |
A numeric vector containing the error metrics.
obj <- list(mean = 10, lower = 8, upper = 12) actual <- 11 geterror(obj, actual)obj <- list(mean = 10, lower = 8, upper = 12) actual <- 11 geterror(obj, actual)
Calculate returns or log-returns for multivariate time series
getreturns(x, type = c("basic", "log"))getreturns(x, type = c("basic", "log"))
x |
Multivariate time series |
type |
Type of return: basic return ("basic") or log-return ("log") |
The returns
returns <- getreturns(EuStockMarkets) log_returns <- getreturns(EuStockMarkets, type = "log") par(mfrow=c(1, 3)) matplot(EuStockMarkets, type = 'l', main = "Closing Prices of \n European stocks (1991-1998)", xlab = "time") matplot(returns, type = 'l', main = "Returns", xlab = "time") matplot(log_returns, type = 'l', main = "Log-returns", xlab = "time")returns <- getreturns(EuStockMarkets) log_returns <- getreturns(EuStockMarkets, type = "log") par(mfrow=c(1, 3)) matplot(EuStockMarkets, type = 'l', main = "Closing Prices of \n European stocks (1991-1998)", xlab = "time") matplot(returns, type = 'l', main = "Returns", xlab = "time") matplot(log_returns, type = 'l', main = "Log-returns", xlab = "time")
Obtain simulations (when relevant) from a selected time series
getsimulations(obj, selected_series, transpose = FALSE)getsimulations(obj, selected_series, transpose = FALSE)
obj |
result from ridge2f (multivariate time series forecast with simulations) |
selected_series |
name of the time series selected |
transpose |
return a transposed time series |
require(fpp2) obj <- ahead::ridge2f(fpp2::insurance, h = 7, type_pi = "bootstrap", B = 5) print(getsimulations(obj, selected_series = "TV.advert")) print(getsimulations(obj, selected_series = "Quotes")) print(getsimulations(obj, selected_series = "TV.advert", transpose = TRUE)) print(getsimulations(obj, selected_series = "Quotes", transpose = TRUE))require(fpp2) obj <- ahead::ridge2f(fpp2::insurance, h = 7, type_pi = "bootstrap", B = 5) print(getsimulations(obj, selected_series = "TV.advert")) print(getsimulations(obj, selected_series = "Quotes")) print(getsimulations(obj, selected_series = "TV.advert", transpose = TRUE)) print(getsimulations(obj, selected_series = "Quotes", transpose = TRUE))
This function implements the Theta method using a Generalized Linear Model (GLM)
glmthetaf( y, h = ifelse(frequency(y) > 1, 2 * frequency(y), 10), level = 95L, fit_func = stats::glm, predict_func = predict, fan = FALSE, x = y, type_pi = c("gaussian", "conformal-split", "conformal-surrogate", "conformal-bootstrap", "conformal-block-bootstrap", "conformal-kde", "conformal-fitdistr"), attention = TRUE, attention_type = c("dot_product", "scaled_dot_product", "cosine", "exponential", "gaussian", "linear", "value_based", "hybrid", "parametric"), attention_method = c("heuristic", "historical"), scale_ctxt = 1, historical_lookback = 50L, max_adjustment = 0.5, B = 250L, nsim = B, ... )glmthetaf( y, h = ifelse(frequency(y) > 1, 2 * frequency(y), 10), level = 95L, fit_func = stats::glm, predict_func = predict, fan = FALSE, x = y, type_pi = c("gaussian", "conformal-split", "conformal-surrogate", "conformal-bootstrap", "conformal-block-bootstrap", "conformal-kde", "conformal-fitdistr"), attention = TRUE, attention_type = c("dot_product", "scaled_dot_product", "cosine", "exponential", "gaussian", "linear", "value_based", "hybrid", "parametric"), attention_method = c("heuristic", "historical"), scale_ctxt = 1, historical_lookback = 50L, max_adjustment = 0.5, B = 250L, nsim = B, ... )
y |
The time series data |
h |
The number of periods to forecast |
level |
The confidence level for the forecast intervals |
fit_func |
The function to use for fitting the GLM (and other models) |
predict_func |
The function to use for predict the other models (not GLMs) |
fan |
Logical flag for fan plot |
x |
The time series data |
attention |
Logical flag for using attention mechanism |
attention_type |
The type of attention mechanism to use |
attention_method |
The method for computing attention-based adjustments: "heuristic" (original) or "historical" (new evidence-based) |
scale_ctxt |
Scaling coefficient for context vector |
historical_lookback |
Number of periods to consider for historical matching |
max_adjustment |
Maximum allowed adjustment percentage (0-1) |
B |
Number of bootstrap replications or number of simulations (alias: nsim) |
nsim |
Alias for B |
... |
Additional arguments to pass to the fit_func |
A forecast object
Loess forecasting
loessf( y, h = 5, level = 95, span = 0.75, degree = 2, type_pi = c("bootstrap", "blockbootstrap", "movingblockbootstrap"), b = NULL, B = 250, type_aggregation = c("mean", "median"), seed = 123 )loessf( y, h = 5, level = 95, span = 0.75, degree = 2, type_pi = c("bootstrap", "blockbootstrap", "movingblockbootstrap"), b = NULL, B = 250, type_aggregation = c("mean", "median"), seed = 123 )
y |
A numeric vector or time series of class |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
span |
the parameter which controls the degree of smoothing |
degree |
the degree of the polynomials to be used, normally 1 or 2. (Degree 0 is also allowed, but see |
type_pi |
Type of prediction interval currently (independent) "bootstrap", (circular) "blockbootstrap", or "movingblockbootstrap" |
b |
block length for circular block bootstrap |
B |
number of bootstrap replications |
type_aggregation |
Type of aggregation, ONLY for bootstrapping; either "mean" or "median" |
seed |
reproducibility seed |
An object of class "forecast"; a list containing the following elements:
model |
A list containing information about the fitted model |
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
x |
The original time series |
residuals |
Residuals from the fitted model |
sims |
Model simulations for bootstrapping |
par(mfrow = c(3, 1)) plot(loessf(Nile, h=20, level=95, B=10)) plot(loessf(Nile, h=20, level=95, B=10, type_pi = "blockbootstrap")) plot(loessf(Nile, h=20, level=95, B=10, type_pi = "movingblockbootstrap"))par(mfrow = c(3, 1)) plot(loessf(Nile, h=20, level=95, B=10)) plot(loessf(Nile, h=20, level=95, B=10, type_pi = "blockbootstrap")) plot(loessf(Nile, h=20, level=95, B=10, type_pi = "movingblockbootstrap"))
LOOCV for Random Vector functional link network model with 2 regularization parameters
loocvridge2f( y, xreg = NULL, h = 5, level = 95, lags = 1, nb_hidden = 5, nodes_sim = c("sobol", "halton", "unif"), activ = c("relu", "sigmoid", "tanh", "leakyrelu", "elu", "linear"), a = 0.01, lambda_1 = 0.1, lambda_2 = 0.1, dropout = 0, type_forecast = c("recursive", "direct"), type_pi = c("gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap", "rvinecopula", "splitconformal"), block_length = NULL, margins = c("gaussian", "empirical", "student"), seed = 1, B = 100L, type_aggregation = c("mean", "median"), centers = NULL, type_clustering = c("kmeans", "hclust"), ym = NULL, cl = 1L, show_progress = TRUE, ... )loocvridge2f( y, xreg = NULL, h = 5, level = 95, lags = 1, nb_hidden = 5, nodes_sim = c("sobol", "halton", "unif"), activ = c("relu", "sigmoid", "tanh", "leakyrelu", "elu", "linear"), a = 0.01, lambda_1 = 0.1, lambda_2 = 0.1, dropout = 0, type_forecast = c("recursive", "direct"), type_pi = c("gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap", "rvinecopula", "splitconformal"), block_length = NULL, margins = c("gaussian", "empirical", "student"), seed = 1, B = 100L, type_aggregation = c("mean", "median"), centers = NULL, type_clustering = c("kmeans", "hclust"), ym = NULL, cl = 1L, show_progress = TRUE, ... )
y |
A multivariate time series of class |
xreg |
External regressors. A data.frame (preferred) or a |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
lags |
Number of lags |
|
Number of nodes in hidden layer |
|
nodes_sim |
Type of simulation for nodes in the hidden layer |
activ |
Activation function |
a |
Hyperparameter for activation function "leakyrelu", "elu" |
lambda_1 |
Regularization parameter for original predictors |
lambda_2 |
Regularization parameter for transformed predictors |
dropout |
dropout regularization parameter (dropping nodes in hidden layer) |
type_forecast |
Recursive or direct forecast |
type_pi |
Type of prediction interval currently "gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap", "splitconformal" (very experimental right now), "rvinecopula" (with Gaussian margins for now, Student-t coming soon) |
block_length |
Length of block for circular or moving block bootstrap |
margins |
Distribution of margins: "gaussian", "empirical", "student" (postponed or
never) for |
seed |
Reproducibility seed for random stuff |
B |
Number of bootstrap replications or number of simulations (yes, 'B' is unfortunate) |
type_aggregation |
Type of aggregation, ONLY for bootstrapping; either "mean" or "median" |
centers |
Number of clusters for |
type_clustering |
"kmeans" (K-Means clustering) or "hclust" (Hierarchical clustering) |
ym |
Univariate time series ( |
cl |
An integer; the number of clusters for parallel execution, for bootstrap |
show_progress |
A boolean; show progress bar for bootstrapping? Default is TRUE. |
... |
An object of class "mtsforecast"; a list containing the following elements:
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
sims |
Model simulations for bootstrapping (basic, or block) |
x |
The original time series |
residuals |
Residuals from the fitted model |
coefficients |
Regression coefficients for |
T. Moudiki
Moudiki, T., Planchet, F., & Cousin, A. (2018).
Multiple time series forecasting using quasi-randomized
functional link neural networks. Risks, 6(1), 22.
require(fpp2) print(ahead::loocvridge2f(fpp2::insurance)) #foo <- function(xx) ahead::loocvridge2f(fpp2::insurance, lambda_1=10^xx[1], lambda_2=10^xx[2]) #(opt <- stats::nlminb(objective=foo, lower=c(-10,-10), upper=c(10,10), start=c(0, 0))) #print(ahead::loocvridge2f(fpp2::insurance, lambda_1=10^opt$par[1], lambda_2=10^opt$par[2]))require(fpp2) print(ahead::loocvridge2f(fpp2::insurance)) #foo <- function(xx) ahead::loocvridge2f(fpp2::insurance, lambda_1=10^xx[1], lambda_2=10^xx[2]) #(opt <- stats::nlminb(objective=foo, lower=c(-10,-10), upper=c(10,10), start=c(0, 0))) #print(ahead::loocvridge2f(fpp2::insurance, lambda_1=10^opt$par[1], lambda_2=10^opt$par[2]))
Generates bootstrap replicates of a time series using the maximum entropy bootstrap algorithm with Rcpp implementation for improved performance. This method is particularly useful for non-stationary time series and preserves the dependence structure of the original data.
meboot( x, reps = 999, trim = list(trim = 0.1, xmin = NULL, xmax = NULL), reachbnd = TRUE, expand.sd = TRUE, force.clt = TRUE, scl.adjustment = FALSE, sym = FALSE, colsubj, coldata, coltimes, ... )meboot( x, reps = 999, trim = list(trim = 0.1, xmin = NULL, xmax = NULL), reachbnd = TRUE, expand.sd = TRUE, force.clt = TRUE, scl.adjustment = FALSE, sym = FALSE, colsubj, coldata, coltimes, ... )
x |
A numeric vector or time series object to be bootstrapped. |
reps |
Number of bootstrap replicates to generate (default: 999). |
trim |
Controls tail behavior. Can be a single numeric value specifying the trim proportion (e.g., 0.10 for 10% trimming), or a list with components:
|
reachbnd |
Logical indicating whether to allow generated values to reach the boundaries xmin and xmax (default: TRUE). |
expand.sd |
Logical indicating whether to expand the standard deviation of the ensemble (default: TRUE). |
force.clt |
Logical indicating whether to force the central limit theorem compliance by centering each replicate (default: TRUE). |
scl.adjustment |
Logical indicating whether to adjust the scale of the ensemble to match the original data's variance (default: FALSE). |
sym |
Logical indicating whether to force symmetry in the maximum entropy density (default: FALSE). |
colsubj |
Deprecated parameter from original meboot (included for compatibility). |
coldata |
Deprecated parameter from original meboot (included for compatibility). |
coltimes |
Deprecated parameter from original meboot (included for compatibility). |
... |
Additional arguments passed to expansion functions. |
The maximum entropy bootstrap algorithm generates replicates that:
Preserve the dependence structure of the original time series
Can handle non-stationary time series
Satisfy the ergodic theorem and central limit theorem
Maintain the mean and autocorrelation structure
The Rcpp implementation provides significant performance improvements over the original R implementation, especially for large datasets and many replications.
A list with the following components:
x |
Original time series data |
ensemble |
Matrix of bootstrap replicates (n x reps) |
xx |
Sorted original data |
z |
Intermediate points between sorted values |
dv |
Absolute differences between consecutive observations |
dvtrim |
Trimmed mean of differences |
xmin |
Lower bound used for generation |
xmax |
Upper bound used for generation |
desintxb |
Interval means satisfying mean-preserving constraint |
ordxx |
Ordering index of original data |
kappa |
Scale adjustment factor (if scl.adjustment = TRUE) |
Vinod, H. D., & Lopez-de-Lacalle, J. (2009). Maximum entropy bootstrap for time series: The meboot R package. Journal of Statistical Software, 29(5), 1-19.
# Basic usage with a time series set.seed(123) x <- ts(rnorm(100), start = c(2000, 1), frequency = 12) result <- meboot(x, reps = 1000) # Plot first few replicates matplot(result$ensemble[, 1:5], type = "l", lty = 1) lines(result$x, col = "black", lwd = 2) # With custom bounds result_bounded <- meboot(x, reps = 100, trim = list(trim = 0.1, xmin = -3, xmax = 3)) # With scale adjustment result_scaled <- meboot(x, reps = 100, scl.adjustment = TRUE)# Basic usage with a time series set.seed(123) x <- ts(rnorm(100), start = c(2000, 1), frequency = 12) result <- meboot(x, reps = 1000) # Plot first few replicates matplot(result$ensemble[, 1:5], type = "l", lty = 1) lines(result$x, col = "black", lwd = 2) # With custom bounds result_bounded <- meboot(x, reps = 100, trim = list(trim = 0.1, xmin = -3, xmax = 3)) # With scale adjustment result_scaled <- meboot(x, reps = 100, scl.adjustment = TRUE)
Forecasting using Machine Leaning models
ml_forecast( y, h, level = 95, lags = 1, fit_func = ahead::ridge, predict_func = predict, xreg = NULL, coeffs = NULL, ... )ml_forecast( y, h, level = 95, lags = 1, fit_func = ahead::ridge, predict_func = predict, xreg = NULL, coeffs = NULL, ... )
y |
A numeric vector or time series of class |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
lags |
Number of lags of the input time series considered in the regression |
fit_func |
Fitting function (Statistical/ML model). Default is Ridge regression. |
predict_func |
Prediction function (Statistical/ML model) |
xreg |
External regressor variable |
coeffs |
Coefficients of the fitted model. If provided, a linear combination with the coefficients is used to compute the prediction. |
... |
additional parameters passed to the fitting function |
An object of class 'forecast'
## Not run: plot(ahead::ml_forecast(AirPassengers, h=20L)) plot(ahead::ml_forecast(AirPassengers, h=25L, lags=20L, fit_func=glmnet::cv.glmnet)) res <- ahead::ml_forecast(USAccDeaths, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(USAccDeaths, fit_func = glmnet::cv.glmnet, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(USAccDeaths, fit_func = e1071::svm, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(mdeaths, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(fdeaths, fit_func = glmnet::cv.glmnet, h=15L, lags=25L) plot(res) res <- ahead::ml_forecast(fdeaths, fit_func = randomForest::randomForest, h=15L, lags=15L) plot(res) ## End(Not run)## Not run: plot(ahead::ml_forecast(AirPassengers, h=20L)) plot(ahead::ml_forecast(AirPassengers, h=25L, lags=20L, fit_func=glmnet::cv.glmnet)) res <- ahead::ml_forecast(USAccDeaths, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(USAccDeaths, fit_func = glmnet::cv.glmnet, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(USAccDeaths, fit_func = e1071::svm, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(mdeaths, h=15L, lags=15L) plot(res) res <- ahead::ml_forecast(fdeaths, fit_func = glmnet::cv.glmnet, h=15L, lags=25L) plot(res) res <- ahead::ml_forecast(fdeaths, fit_func = randomForest::randomForest, h=15L, lags=15L) plot(res) ## End(Not run)
Conformalized Forecasting using Machine Learning (and statistical) models with ARCH effects
mlarchf( y, h = 10L, mean_model = forecast::auto.arima, model_residuals = forecast::meanf, fit_func = ahead::ridge, predict_func = predict, conformal = FALSE, lags_vol = 10L, type_pi = c("surrogate", "bootstrap", "kde"), type_sim_conformalize = c("surrogate", "block-bootstrap", "bootstrap", "kde", "fitdistr"), ml_method = NULL, level = 95, B = 250L, ml = TRUE, stat_model = NULL, n_clusters = 0, clustering_dist = "euclidean", clustering_method = "kmeans", ... )mlarchf( y, h = 10L, mean_model = forecast::auto.arima, model_residuals = forecast::meanf, fit_func = ahead::ridge, predict_func = predict, conformal = FALSE, lags_vol = 10L, type_pi = c("surrogate", "bootstrap", "kde"), type_sim_conformalize = c("surrogate", "block-bootstrap", "bootstrap", "kde", "fitdistr"), ml_method = NULL, level = 95, B = 250L, ml = TRUE, stat_model = NULL, n_clusters = 0, clustering_dist = "euclidean", clustering_method = "kmeans", ... )
y |
A numeric vector or time series of class |
h |
Forecasting horizon |
mean_model |
Function to fit the mean model (default: |
model_residuals |
Function to model the residuals (default: |
fit_func |
Fitting function for the variance model (default: |
predict_func |
Prediction function for the variance model (default: |
conformal |
A boolean; either use conformal prediction or not |
lags_vol |
Number of lags for squared residuals regression |
type_pi |
Type of prediction interval ("kde", "surrogate", or "bootstrap") for volatility modeling |
type_sim_conformalize |
Type of simulation for conformalization of standardized residuals ("block-bootstrap", "surrogate", "kde", "bootstrap", or "fitdistr") |
ml_method |
caret package Machine learning method to use, if |
level |
Confidence level for prediction intervals |
B |
Number of bootstrap replications or simulations |
ml |
If |
stat_model |
A statistical model, e.g |
n_clusters |
Number of clusters for residuals (default is 0) for |
clustering_dist |
Only if |
clustering_method |
Only if |
... |
Additional parameters to be passed to |
A forecast object containing predictions and prediction intervals
y <- fpp2::goog200 par(mfrow=c(2, 2)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=TRUE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=FALSE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=TRUE, stat_model=forecast::thetaf, ml = FALSE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=FALSE, stat_model=forecast::thetaf, ml = FALSE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) par(mfrow=c(1, 1))y <- fpp2::goog200 par(mfrow=c(2, 2)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=TRUE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=FALSE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=TRUE, stat_model=forecast::thetaf, ml = FALSE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) (obj_ridge <- ahead::mlarchf(y, h=20L, B=500L, conformal=FALSE, stat_model=forecast::thetaf, ml = FALSE)) plot(obj_ridge) print(mean(obj_ridge$upper - obj_ridge$lower)) par(mfrow=c(1, 1))
Conformalized Forecasting using Machine Leaning models
mlf( y, h = 5, level = 95, lags = 15L, fit_func = ahead::ridge, predict_func = predict, stack = FALSE, stacking_models = NULL, coeffs = NULL, type_pi = c("surrogate", "bootstrap", "kde"), B = 250L, agg = c("mean", "median"), seed = 123, show_progress = TRUE, ... )mlf( y, h = 5, level = 95, lags = 15L, fit_func = ahead::ridge, predict_func = predict, stack = FALSE, stacking_models = NULL, coeffs = NULL, type_pi = c("surrogate", "bootstrap", "kde"), B = 250L, agg = c("mean", "median"), seed = 123, show_progress = TRUE, ... )
y |
A numeric vector or time series of class |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
lags |
Number of lags of the input time series considered in the regression |
fit_func |
Fitting function (Statistical/ML model). Default is Ridge regression. |
predict_func |
Prediction function (Statistical/ML model) |
stack |
Boolean, use stacking regression or not |
stacking_models |
A list of |
coeffs |
Coefficients of the fitted model. If provided, a linear combination with the coefficients is used to compute the prediction. |
type_pi |
Type of prediction interval |
B |
Number of bootstrap replications or number of simulations |
agg |
"mean" or "median" (aggregation method) |
show_progress |
show progress bar for stacking, if |
... |
additional parameters passed to the fitting function |
An object of class 'forecast'
This function plots the computed sensitivity effects over the forecast horizon.
plot_dynrmf_sensitivity( sensitivity_results, title = "First-Order Sensitivity Effect", y_label = "Effect (DeltaY / DeltaX)", x_label = "Forecast Horizon" )plot_dynrmf_sensitivity( sensitivity_results, title = "First-Order Sensitivity Effect", y_label = "Effect (DeltaY / DeltaX)", x_label = "Forecast Horizon" )
sensitivity_results |
Output list from |
title |
Character. Plot title. |
y_label |
Character. Label for y-axis. |
x_label |
Character. Label for x-axis (default "Forecast Horizon"). |
ggplot2 object displaying the sensitivity effects.
## Not run: plot1 <- plot_dynrmf_sensitivity(res, title = "Sensitivity of Consumption to Income", y_label = "Effect (DeltaConsumption / DeltaIncome)") print(plot1) ## End(Not run)## Not run: plot1 <- plot_dynrmf_sensitivity(res, title = "Sensitivity of Consumption to Income", y_label = "Effect (DeltaConsumption / DeltaIncome)") print(plot1) ## End(Not run)
Aggregates SHAP values across horizons and draws a horizontal waterfall.
plot_dynrmf_shap_waterfall( shap, title = "SHAP Waterfall", agg_fun = sum, colors = list(pos = "#4575b4", neg = "#d73027", base = "#878787") )plot_dynrmf_shap_waterfall( shap, title = "SHAP Waterfall", agg_fun = sum, colors = list(pos = "#4575b4", neg = "#d73027", base = "#878787") )
shap |
A "dynrmf_shap" object from dynrmf_shap() |
title |
Plot title string |
agg_fun |
How to aggregate across horizons: sum (default) or mean |
colors |
Named list with elements "pos", "neg", "base" |
A ggplot2 object
Plot forecast results with simulation intervals
plot_hts_forecast(result, hts, series_type = "total", series_idx = 1)plot_hts_forecast(result, hts, series_type = "total", series_idx = 1)
Plot simulation paths
plot_simulations( result, series_type = "total", series_idx = 1, n_paths = 50, main = NULL )plot_simulations( result, series_type = "total", series_idx = 1, n_paths = 50, main = NULL )
Produces plots for the results of a forecast combination method. Either
an actual vs. fitted plot (which = 1) or a barplot of the combination weights
(which = 2).
## S3 method for class 'foreccomb_res' plot(x, which = 1, ...)## S3 method for class 'foreccomb_res' plot(x, which = 1, ...)
x |
An object of class 'foreccomb_res'. |
which |
Type of plot: 1 = actual vs. fitted, 2 = combination weights. |
... |
Other arguments passing to |
A plot for the foreccomb_res class.
adapted from Christoph E. Weiss and Gernot R. Roetzer (ForecastComb)
foreccomb,
summary.foreccomb_res
Plot multivariate time series forecast or residuals
## S3 method for class 'mtsforecast' plot(x, selected_series, type = c("pi", "dist", "sims"), level = 95, ...)## S3 method for class 'mtsforecast' plot(x, selected_series, type = c("pi", "dist", "sims"), level = 95, ...)
x |
result from |
selected_series |
name of the time series selected for plotting |
type |
"pi": basic prediction intervals; "dist": a distribution of predictions; "sims": the simulations |
level |
confidence levels for prediction intervals |
... |
additional parameters to be passed to |
require(fpp2) fit_obj_VAR <- ahead::varf(fpp2::insurance, lags = 2, h = 10, level = 95) fit_obj_ridge2 <- ahead::ridge2f(fpp2::insurance, lags = 2, h = 10, level = 95) par(mfrow=c(2, 2)) plot(fit_obj_VAR, "Quotes") plot(fit_obj_VAR, "TV.advert") plot(fit_obj_ridge2, "Quotes") plot(fit_obj_ridge2, "TV.advert") obj <- ahead::ridge2f(fpp2::insurance, h = 10, type_pi = "blockbootstrap", block_length=5, B = 10) par(mfrow=c(1, 2)) plot(obj, selected_series = "Quotes", type = "sims", main = "Predictive simulation for Quotes") plot(obj, selected_series = "TV.advert", type = "sims", main = "Predictive simulation for TV.advert") par(mfrow=c(1, 2)) plot(obj, selected_series = "Quotes", type = "dist", main = "Predictive simulation for Quotes") plot(obj, selected_series = "TV.advert", type = "dist", main = "Predictive simulation for TV.advert")require(fpp2) fit_obj_VAR <- ahead::varf(fpp2::insurance, lags = 2, h = 10, level = 95) fit_obj_ridge2 <- ahead::ridge2f(fpp2::insurance, lags = 2, h = 10, level = 95) par(mfrow=c(2, 2)) plot(fit_obj_VAR, "Quotes") plot(fit_obj_VAR, "TV.advert") plot(fit_obj_ridge2, "Quotes") plot(fit_obj_ridge2, "TV.advert") obj <- ahead::ridge2f(fpp2::insurance, h = 10, type_pi = "blockbootstrap", block_length=5, B = 10) par(mfrow=c(1, 2)) plot(obj, selected_series = "Quotes", type = "sims", main = "Predictive simulation for Quotes") plot(obj, selected_series = "TV.advert", type = "sims", main = "Predictive simulation for TV.advert") par(mfrow=c(1, 2)) plot(obj, selected_series = "Quotes", type = "dist", main = "Predictive simulation for Quotes") plot(obj, selected_series = "TV.advert", type = "dist", main = "Predictive simulation for TV.advert")
Plot method for synthetic_ts objects
## S3 method for class 'synthetic_ts' plot(x, type = "series", which_sim = 1, series_index = 1, ...)## S3 method for class 'synthetic_ts' plot(x, type = "series", which_sim = 1, series_index = 1, ...)
x |
synthetic_ts object |
type |
Type of plot: "series", "residuals", "acf", "density", "qq" |
which_sim |
Which simulation to plot (for type = "series") |
series_index |
Which series to plot (for multivariate data) |
... |
Additional arguments passed to plot functions |
dynrmf)Predict univariate time series using caret ML model(for use with dynrmf)
predict_func(obj, newx)predict_func(obj, newx)
x |
A matrix of predictors |
y |
A vector of responses |
method |
The caret method to use for fitting the model |
verbose |
Whether to print the model summary |
A model object
prediction method for class ‘foreccomb_res’. Uses the previously created forecast combination
result to predict the combination for a newly provided prediction dataset.
## S3 method for class 'foreccomb_res' predict(object, newpreds, newobs = NULL, simplify = TRUE, byrow = FALSE, ...)## S3 method for class 'foreccomb_res' predict(object, newpreds, newobs = NULL, simplify = TRUE, byrow = FALSE, ...)
object |
An object of class 'foreccomb'. Contains training set (actual values + matrix of model forecasts) and optionally a test set. |
newpreds |
A matrix or multivariate time series; contains individual model forecasts if a test set is used (optional). Does not
require specification of |
newobs |
A vector or univariate time series; contains ‘actual values’ if a test set is used (optional). |
simplify |
logical. The default ( |
byrow |
logical. The default ( |
... |
potential further arguments (require by generic) |
Adapted from Christoph E. Weiss and Gernot R. Roetzer (ForecastComb)
Print method for summary.synthetic_ts objects
## S3 method for class 'summary.synthetic_ts' print(x, ...)## S3 method for class 'summary.synthetic_ts' print(x, ...)
x |
summary.synthetic_ts object |
... |
Additional arguments |
Simulate from parametric distribution
rfitdistr(x, n = length(x), p = 1)rfitdistr(x, n = length(x), p = 1)
Generates samples from a Gaussian kernel density estimate of a numeric vector.
rgaussiandens(x, n = length(x), p = 1, seed = 32124)rgaussiandens(x, n = length(x), p = 1, seed = 32124)
x |
Numeric vector to estimate density from. |
n |
Integer. Number of samples per replicate (default: length of x). |
p |
Integer. Number of replicates (default: 1). |
seed |
Integer. Random seed for reproducibility (default: 123). |
method |
Character. Sampling method: "antithetic" or "traditional". |
A vector or matrix of samples from the estimated density.
x <- rnorm(10) rgaussiandens(x, n = 10, p = 3)x <- rnorm(10) rgaussiandens(x, n = 10, p = 3)
Random Vector functional link network model with 2 regularization parameters
ridge2f( y, h = 5, level = 95, xreg = NULL, lags = 1, nb_hidden = 5, nodes_sim = c("sobol", "halton", "unif"), activ = c("relu", "sigmoid", "tanh", "leakyrelu", "elu", "linear"), a = 0.01, lambda_1 = 0.1, lambda_2 = 0.1, dropout = 0, type_forecast = c("recursive", "direct"), type_pi = c("gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap", "rvinecopula", "conformal-split", "conformal-bootstrap", "conformal-block-bootstrap"), block_length = NULL, margins = c("gaussian", "empirical"), seed = 1, B = 100L, type_aggregation = c("mean", "median"), centers = NULL, type_clustering = c("kmeans", "hclust"), ym = NULL, cl = 1L, show_progress = TRUE, ... )ridge2f( y, h = 5, level = 95, xreg = NULL, lags = 1, nb_hidden = 5, nodes_sim = c("sobol", "halton", "unif"), activ = c("relu", "sigmoid", "tanh", "leakyrelu", "elu", "linear"), a = 0.01, lambda_1 = 0.1, lambda_2 = 0.1, dropout = 0, type_forecast = c("recursive", "direct"), type_pi = c("gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap", "rvinecopula", "conformal-split", "conformal-bootstrap", "conformal-block-bootstrap"), block_length = NULL, margins = c("gaussian", "empirical"), seed = 1, B = 100L, type_aggregation = c("mean", "median"), centers = NULL, type_clustering = c("kmeans", "hclust"), ym = NULL, cl = 1L, show_progress = TRUE, ... )
y |
A univariate of multivariate time series of class |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
xreg |
External regressors. A data.frame (preferred) or a |
lags |
Number of lags |
|
Number of nodes in hidden layer |
|
nodes_sim |
Type of simulation for nodes in the hidden layer |
activ |
Activation function |
a |
Hyperparameter for activation function "leakyrelu", "elu" |
lambda_1 |
Regularization parameter for original predictors |
lambda_2 |
Regularization parameter for transformed predictors |
dropout |
dropout regularization parameter (dropping nodes in hidden layer) |
type_forecast |
Recursive or direct forecast |
type_pi |
Type of prediction interval currently "gaussian", "bootstrap", "blockbootstrap", "movingblockbootstrap", "rvinecopula" (with Gaussian and empirical margins for now) "conformal-split", "conformal-bootstrap", "conformal-block-bootstrap" |
block_length |
Length of block for circular or moving block bootstrap |
margins |
Distribution of margins: "gaussian", "empirical" for |
seed |
Reproducibility seed for random stuff |
B |
Number of bootstrap replications or number of simulations (yes, 'B' is unfortunate) |
type_aggregation |
Type of aggregation, ONLY for bootstrapping; either "mean" or "median" |
centers |
Number of clusters for |
type_clustering |
"kmeans" (K-Means clustering) or "hclust" (Hierarchical clustering) |
ym |
Univariate time series ( |
cl |
An integer; the number of clusters for parallel execution, for bootstrap |
show_progress |
A boolean; show progress bar for bootstrapping? Default is TRUE. |
... |
An object of class "mtsforecast"; a list containing the following elements:
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
sims |
Model simulations for bootstrapping (basic, or block) |
x |
The original time series |
residuals |
Residuals from the fitted model |
coefficients |
Regression coefficients for |
T. Moudiki
Moudiki, T., Planchet, F., & Cousin, A. (2018).
Multiple time series forecasting using quasi-randomized
functional link neural networks. Risks, 6(1), 22.
require(fpp2) print(ahead::ridge2f(fpp2::insurance)$mean) res <- ahead::ridge2f(fpp2::insurance, h=10, lags=2) par(mfrow=c(1, 2)) plot(res, "Quotes") plot(res, "TV.advert") # include a trend (just for the example) xreg <- as.numeric(time(fpp2::insurance)) res2 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=2) par(mfrow=c(1, 2)) plot(res2, "Quotes") plot(res2, "TV.advert") # block bootstrap xreg <- as.numeric(time(fpp2::insurance)) res3 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=1L, type_pi = "bootstrap", B=10) res5 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=1L, type_pi = "blockbootstrap", B=10, block_length = 4) print(res3$sims[[2]]) print(res5$sims[[2]]) par(mfrow=c(2, 2)) plot(res3, "Quotes") plot(res3, "TV.advert") plot(res5, "Quotes") plot(res5, "TV.advert") # moving block bootstrap xreg <- as.numeric(time(fpp2::insurance)) res6 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=1L, type_pi = "movingblockbootstrap", B=10, block_length = 4) print(res6$sims[[2]]) par(mfrow=c(1, 2)) plot(res6, "Quotes") plot(res6, "TV.advert")require(fpp2) print(ahead::ridge2f(fpp2::insurance)$mean) res <- ahead::ridge2f(fpp2::insurance, h=10, lags=2) par(mfrow=c(1, 2)) plot(res, "Quotes") plot(res, "TV.advert") # include a trend (just for the example) xreg <- as.numeric(time(fpp2::insurance)) res2 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=2) par(mfrow=c(1, 2)) plot(res2, "Quotes") plot(res2, "TV.advert") # block bootstrap xreg <- as.numeric(time(fpp2::insurance)) res3 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=1L, type_pi = "bootstrap", B=10) res5 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=1L, type_pi = "blockbootstrap", B=10, block_length = 4) print(res3$sims[[2]]) print(res5$sims[[2]]) par(mfrow=c(2, 2)) plot(res3, "Quotes") plot(res3, "TV.advert") plot(res5, "Quotes") plot(res5, "TV.advert") # moving block bootstrap xreg <- as.numeric(time(fpp2::insurance)) res6 <- ahead::ridge2f(fpp2::insurance, xreg=xreg, h=10, lags=1L, type_pi = "movingblockbootstrap", B=10, block_length = 4) print(res6$sims[[2]]) par(mfrow=c(1, 2)) plot(res6, "Quotes") plot(res6, "TV.advert")
Simulate multivariate data
rmultivariate( data, method = c("bootstrap", "block-bootstrap"), n = 100L, block_size = 5 )rmultivariate( data, method = c("bootstrap", "block-bootstrap"), n = 100L, block_size = 5 )
data |
A numeric vector or matrix. |
method |
The method to use for sampling. |
n |
The number of samples to draw. |
block_size |
The size of the blocks to use for the block bootstrap. |
... |
Additional arguments to pass to the density function. |
Simulate using surrogate data
rsurrogate(x, n = length(x), p = 1, seed = 123)rsurrogate(x, n = length(x), p = 1, seed = 123)
Sequential split conformal prediction for hierarchical forecasting Returns simulations at both total and bottom levels
sequential_conformal_hts( hts, h = 12, split_ratio = 0.6, n_sim = 1000, alpha = 0.1 )sequential_conformal_hts( hts, h = 12, split_ratio = 0.6, n_sim = 1000, alpha = 0.1 )
hts |
2-level HTS object |
h |
Forecast horizon |
split_ratio |
Training/calibration split ratio |
n_sim |
Number of simulation samples |
alpha |
Significance level for prediction intervals |
Forecast with prediction intervals and simulations
Simple ETS forecast function
simple_forecast(y, h)simple_forecast(y, h)
y |
Time series vector |
h |
Forecast horizon |
Forecast vector
This function allows to obtain predictive simulations from any forecasting function.
simulator( FUN, y, h, level = 95, method = c("block-bootstrap", "surrogate", "kde", "bootstrap", "fitdistr", "meboot"), nsim = 100L, block_size = 5, seed = 123L, B = NULL, ... )simulator( FUN, y, h, level = 95, method = c("block-bootstrap", "surrogate", "kde", "bootstrap", "fitdistr", "meboot"), nsim = 100L, block_size = 5, seed = 123L, B = NULL, ... )
FUN |
A forecasting function. |
y |
A time series ( |
h |
Forecasting horizon. |
level |
Confidence level. |
method |
Method to be used for the simulation of in-sample residuals. |
nsim |
Number of simulations. |
block_size |
Block size for block-bootstrap. |
seed |
Seed for reproducibility. |
B |
Alias for |
... |
Additional arguments to be passed to the forecasting function. |
An object of class forecast, with a slot sims
y <- fdeaths h <- 25L obj <- simulator(FUN=forecast::thetaf, y, h); plot(obj) obj <- simulator(FUN=forecast::auto.arima, y, h); plot(obj) obj <- simulator(FUN=forecast::ets, y, h); plot(obj)y <- fdeaths h <- 25L obj <- simulator(FUN=forecast::thetaf, y, h); plot(obj) obj <- simulator(FUN=forecast::auto.arima, y, h); plot(obj) obj <- simulator(FUN=forecast::ets, y, h); plot(obj)
Partition a time series object
splitts(y, split_prob = 0.5, return_indices = FALSE, ...)splitts(y, split_prob = 0.5, return_indices = FALSE, ...)
y |
A time series object |
split_prob |
Splitting ratio |
return_indices |
if TRUE, returns series' indices, otherwise, time series objects |
misc::splitts(ts(1:10))misc::splitts(ts(1:10))
This function performs a two-stage stacked forecast using ridge regression. In the first stage, it generates base forecasts from the training portion of the time series. In the second stage, these forecasts are used as external regressors to produce the final forecast on the testing set.
stackridge2f(y, h = 5L, split_fraction = 0.5, ...)stackridge2f(y, h = 5L, split_fraction = 0.5, ...)
y |
A multivariate time series object. |
h |
Integer. Forecast horizon for the second-stage (final) forecast. Defaults to 5. |
split_fraction |
Numeric between 0 and 1. Fraction of the time series used for training. The rest is used for testing and generating stacking features. Defaults to 0.5. |
... |
Additional arguments passed to |
The function works as follows:
Split the time series into training and testing sets according to
split_fraction.
Generate base forecasts from the training set using ahead::ridge2f
Use these base forecasts as external regressors (xreg) to predict
the testing set with a second ahead::ridge2f model.
This approach allows stacking of forecasts to potentially improve accuracy by leveraging the predictions of multiple first-stage models.
For multivariate base learners (e.g., tslm with multiple dependent
variables), the function automatically extracts forecasts from each series
and combines them into the feature matrix.
An object returned by ahead::ridge2f for the stacked forecast,
typically a list including mean and prediction intervals.
## Not run: # Univariate example with default ridge2f base learner result1 <- stackridge2f(fpp2::insurance, h = 10, split_fraction = 0.5) ## End(Not run)## Not run: # Univariate example with default ridge2f base learner result1 <- stackridge2f(fpp2::insurance, h = 10, split_fraction = 0.5) ## End(Not run)
summary method for class ‘foreccomb_res’. Includes information about combination method,
combination weights assigned to the individual forecast models, as well as an accuracy evaluation of the combined
forecast.
## S3 method for class 'foreccomb_res' summary(object, ...) ## S3 method for class 'foreccomb_res_summary' print(x, ...)## S3 method for class 'foreccomb_res' summary(object, ...) ## S3 method for class 'foreccomb_res_summary' print(x, ...)
object |
An object of class 'foreccomb'. Contains training set (actual values + matrix of model forecasts) and optionally a test set. |
... |
potential further arguments (require by generic) |
x |
An object of class 'foreccomb'. Contains training set (actual values + matrix of model forecasts) and optionally a test set. |
Christoph E. Weiss and Gernot R. Roetzer
foreccomb,
plot.foreccomb_res,
Summary method for synthetic_ts objects
## S3 method for class 'synthetic_ts' summary(object, ...)## S3 method for class 'synthetic_ts' summary(object, ...)
object |
synthetic_ts object |
... |
Additional arguments |
Top-down forecast using historical proportions
topdown_forecast(hts, h)topdown_forecast(hts, h)
hts |
2-level HTS object |
h |
Forecast horizon |
List with total and bottom forecasts
Vector Autoregressive model adapted from vars::VAR (only for benchmarking)
varf( y, h = 5, level = 95, lags = 1, type_VAR = c("const", "trend", "both", "none"), ... )varf( y, h = 5, level = 95, lags = 1, type_VAR = c("const", "trend", "both", "none"), ... )
y |
A multivariate time series of class |
h |
Forecasting horizon |
level |
Confidence level for prediction intervals |
lags |
Number of lags |
type_VAR |
Type of deterministic regressors to include. |
... |
Additional parameters to be passed to vars::VAR. |
An object of class "mtsforecast"; a list containing the following elements:
method |
The name of the forecasting method as a character string |
mean |
Point forecasts for the time series |
lower |
Lower bound for prediction interval |
upper |
Upper bound for prediction interval |
x |
The original time series |
residuals |
Residuals from the fitted model |
T. Moudiki
Bernhard Pfaff (2008). VAR, SVAR and SVEC Models: Implementation
Within R Package vars. Journal of Statistical Software 27(4). URL
http://www.jstatsoft.org/v27/i04/.
Pfaff, B. (2008) Analysis of Integrated and Cointegrated Time Series with R. Second Edition. Springer, New York. ISBN 0-387-27960-1
require(fpp2) print(varf(fpp2::insurance, lags=2, h=10))require(fpp2) print(varf(fpp2::insurance, lags=2, h=10))