Title:	Conformal Prediction Methods for Multistep-Ahead Time Series Forecasting
Version:	0.1.0
Description:	Methods and tools for performing multistep-ahead time series forecasting using conformal prediction methods including classical conformal prediction, adaptive conformal prediction, conformal PID (Proportional-Integral-Derivative) control, and autocorrelated multistep-ahead conformal prediction. The methods were described by Wang and Hyndman (2024) <doi:10.48550/arXiv.2410.13115>.
License:	GPL-3
URL:	https://github.com/xqnwang/conformalForecast, https://xqnwang.github.io/conformalForecast/
BugReports:	https://github.com/xqnwang/conformalForecast/issues
Imports:	forecast, ggdist, rlang, stats, zoo
Suggests:	dplyr, ggplot2, knitr, rmarkdown, testthat (≥ 3.0.0), tibble, tsibble
Config/testthat/edition:	3
Encoding:	UTF-8
RoxygenNote:	7.3.2
VignetteBuilder:	knitr
Depends:	R (≥ 4.1.0)
NeedsCompilation:	no
Packaged:	2025-09-30 08:45:07 UTC; xiaoqianwang
Author:	Xiaoqian Wang [aut, cre, cph], Rob Hyndman [aut]
Maintainer:	Xiaoqian Wang <Xiaoqian.Wang@amss.ac.cn>
Repository:	CRAN
Date/Publication:	2025-10-06 08:10:12 UTC

conformalForecast: Conformal Prediction Methods for Multistep-Ahead Time Series Forecasting

Description

Methods and tools for performing multistep-ahead time series forecasting using conformal prediction methods including classical conformal prediction, adaptive conformal prediction, conformal PID (Proportional-Integral-Derivative) control, and autocorrelated multistep-ahead conformal prediction. The methods were described by Wang and Hyndman (2024) doi:10.48550/arXiv.2410.13115.

Author(s)

Maintainer: Xiaoqian Wang Xiaoqian.Wang@amss.ac.cn (ORCID) [copyright holder]

Authors:

• Rob Hyndman Rob.Hyndman@monash.edu (ORCID)

See Also

Useful links:

• https://github.com/xqnwang/conformalForecast

• https://xqnwang.github.io/conformalForecast/

• Report bugs at https://github.com/xqnwang/conformalForecast/issues

Point estimate accuracy measures

Description

Accuracy measures for point forecast residuals.

Usage

ME(resid, na.rm = TRUE)

MAE(resid, na.rm = TRUE, ...)

MSE(resid, na.rm = TRUE, ...)

RMSE(resid, na.rm = TRUE, ...)

MPE(resid, actual, na.rm = TRUE, ...)

MAPE(resid, actual, na.rm = TRUE, ...)

MASE(
  resid,
  train,
  demean = FALSE,
  na.rm = TRUE,
  period,
  d = period == 1,
  D = period > 1,
  ...
)

RMSSE(
  resid,
  train,
  demean = FALSE,
  na.rm = TRUE,
  period,
  d = period == 1,
  D = period > 1,
  ...
)

point_measures

Arguments

Format

An object of class list of length 8.

Value

For the individual functions (ME, MAE, MSE, RMSE, MPE, MAPE, MASE, RMSSE), returns a single numeric scalar giving the requested accuracy measure.

For the exported object point_measures, returns a named list of functions that can be supplied to higher-level accuracy routines.

Examples

# Toy residuals and data
set.seed(1)
y_train <- rnorm(50)
y_test  <- rnorm(10)
fcast   <- y_test + rnorm(10, sd = 0.2)
resid   <- y_test - fcast

# Basic measures
ME(resid)
MAE(resid)
RMSE(resid)

# Percentage measures require 'actual'
MPE(resid, actual = y_test)
MAPE(resid, actual = y_test)

# Scaled measures require training data (and seasonal period if applicable)
MASE(resid, train = y_train, period = 1)
RMSSE(resid, train = y_train, period = 1)

Interval estimate accuracy measures

Description

Accuracy measures for interval forecasts.

Usage

MSIS(
  lower,
  upper,
  actual,
  train,
  level = 95,
  period,
  d = period == 1,
  D = period > 1,
  na.rm = TRUE,
  ...
)

winkler_score(lower, upper, actual, level = 95, na.rm = TRUE, ...)

interval_measures

Arguments

Format

An object of class list of length 2.

Value

For winkler_score and MSIS, returns a single numeric scalar giving the average interval score (Winkler or mean scaled interval score).

For the exported object interval_measures, returns a named list of functions that can be supplied to higher-level accuracy routines.

Examples

set.seed(1)
actual <- rnorm(10)
lower  <- actual - runif(10, 0.5, 1)
upper  <- actual + runif(10, 0.5, 1)
train  <- rnorm(50)

# Winkler score at 95%
winkler_score(lower, upper, actual, level = 95)

# Mean scaled interval score (needs training data and period)
MSIS(lower, upper, actual, train, level = 95, period = 1)

Accuracy measures for a cross-validation model and a conformal prediction model

Description

Return range of summary measures of the out-of-sample forecast accuracy. If x is given, the function also measures test set forecast accuracy. If x is not given, the function only produces accuracy measures on validation set.

Usage

## Default S3 method:
accuracy(
  object,
  x,
  CV = TRUE,
  period = NULL,
  measures = interval_measures,
  byhorizon = FALSE,
  ...
)

Arguments

Details

The measures calculated are:

• ME: Mean Error

• MAE: Mean Absolute Error

• MSE: Mean Squared Error

• RMSE: Root Mean Squared Error

• MPE: Mean Percentage Error

• MAPE: Mean Absolute Percentage Error

• MASE: Mean Absolute Scaled Error

• RMSSE: Root Mean Squared Scaled Error

• winkler_score: Winkler Score

• MSIS: Mean Scaled Interval Score

Value

A matrix giving mean out-of-sample forecast accuracy measures.

See Also

point_measures, interval_measures

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting with a rolling window
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# Out-of-sample forecast accuracy on validation set
accuracy(fc, measures = point_measures, byhorizon = TRUE)
accuracy(fc, measures = interval_measures, level = 95, byhorizon = TRUE)

# Out-of-sample forecast accuracy on test set
accuracy(fc, x = c(1, 0.5, 0), measures = interval_measures,
         level = 95, byhorizon = TRUE)

Autocorrelated multistep-ahead conformal prediction method

Description

Compute prediction intervals and other information by applying the Autocorrelated Multistep-ahead Conformal Prediction (AcMCP) method. The method can only deal with asymmetric nonconformity scores, i.e., forecast errors.

Usage

acmcp(
  object,
  alpha = 1 - 0.01 * object$level,
  ncal = 10,
  rolling = FALSE,
  integrate = TRUE,
  scorecast = TRUE,
  lr = 0.1,
  Tg = NULL,
  delta = NULL,
  Csat = 2/pi * (ceiling(log(Tg) * delta) - 1/log(Tg)),
  KI = max(abs(object$errors), na.rm = TRUE),
  ...
)

Arguments

Details

Similar to the PID method, the AcMCP method also integrates three modules (P, I, and D) to form the final iteration. However, instead of performing conformal prediction for each individual forecast horizon h separately, AcMCP employs a combination of an MA(h-1) model and a linear regression model of e_{t+h|t} on e_{t+h-1|t},\dots,e_{t+1|t} as the scorecaster. This allows the AcMCP method to capture the relationship between the h-step ahead forecast error and past errors.

Value

A list of class c("acmcp", "cpforecast", "forecast") with the following components:

If mean is included in the object, the components mean, lower, and upper will also be returned, showing the information about the test set forecasts generated using all available observations.

References

Wang, X., and Hyndman, R. J. (2024). "Online conformal inference for multi-step time series forecasting", arXiv preprint arXiv:2410.13115.

See Also

pid

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# AcMCP setup
Tg <- 200; delta <- 0.01
Csat <- 2 / pi * (ceiling(log(Tg) * delta) - 1 / log(Tg))
KI <- 2
lr <- 0.1

# AcMCP with integrator and scorecaster
acmcpfc <- acmcp(fc, ncal = 50, rolling = TRUE,
             integrate = TRUE, scorecast = TRUE,
             lr = lr, KI = KI, Csat = Csat)
print(acmcpfc)
summary(acmcpfc)

Adaptive conformal prediction method

Description

Compute prediction intervals and other information by applying the adaptive conformal prediction (ACP) method.

Usage

acp(
  object,
  alpha = 1 - 0.01 * object$level,
  gamma = 0.005,
  symmetric = FALSE,
  ncal = 10,
  rolling = FALSE,
  quantiletype = 1,
  update = FALSE,
  na.rm = TRUE,
  ...
)

Arguments

Details

The ACP method considers the online update:

\alpha_{t+h|t}:=\alpha_{t+h-1|t-1}+\gamma(\alpha-\mathrm{err}_{t|t-h}),

for each individual forecast horizon h, respectively, where \mathrm{err}_{t|t-h}=1 if s_{t|t-h}>q_{t|t-h}, and \mathrm{err}_{t|t-h}=0 if s_{t|t-h} \leq q_{t|t-h}.

Value

A list of class c("acp", "cpforecast", "forecast") with the following components:

If mean is included in the object, the components mean, lower, and upper will also be returned, showing the information about the forecasts generated using all available observations.

References

Gibbs, I., and Candes, E. (2021). "Adaptive conformal inference under distribution shift", Advances in Neural Information Processing Systems, 34, 1660–1672.

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# ACP with asymmetric nonconformity scores and rolling calibration sets
acpfc <- acp(fc, symmetric = FALSE, gamma = 0.005, ncal = 50, rolling = TRUE)
print(acpfc)
summary(acpfc)

Conformal prediction

Description

This function allows you to specify the method used to perform conformal prediction.

Usage

conformal(object, ...)

## S3 method for class 'cvforecast'
conformal(object, method = c("scp", "acp", "pid", "acmcp"), ...)

Arguments

Value

An object whose class depends on the method invoked.

See Also

scp, acp, pid, and acmcp

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# Classical conformal prediction with equal weights
scpfc <- conformal(fc, method = "scp", symmetric = FALSE, ncal = 50, rolling = TRUE)
summary(scpfc)

# ACP with asymmetric nonconformity scores and rolling calibration sets
acpfc <- conformal(fc, method = "acp", symmetric = FALSE, gamma = 0.005,
                   ncal = 50, rolling = TRUE)
summary(acpfc)

Calculate interval forecast coverage

Description

Calculate the mean coverage and the ifinn matrix for prediction intervals on validation set. If window is not NULL, a matrix of the rolling means of interval forecast coverage is also returned.

Usage

coverage(object, ..., level = 95, window = NULL, na.rm = FALSE)

Arguments

Value

A list of class "coverage" with the following components:

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting with a rolling window
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# Mean and rolling mean coverage for interval forecasts on validation set
cov_fc <- coverage(fc, level = 95, window = 50)
str(cov_fc)

Time series cross-validation forecasting

Description

Compute forecasts and other information by applying forecastfun to subsets of the time series y using a rolling forecast origin.

Usage

cvforecast(
  y,
  forecastfun,
  h = 1,
  level = c(80, 95),
  forward = TRUE,
  xreg = NULL,
  initial = 1,
  window = NULL,
  ...
)

Arguments

Details

Let y denote the time series y_1,\dots,y_T and let t_0 denote the initial period.

Suppose forward = TRUE. If window is NULL, forecastfun is applied successively to the subset time series y_{1},\dots,y_t, for t=t_0,\dots,T, generating forecasts \hat{y}_{t+1|t},\dots,\hat{y}_{t+h|t}. If window is not NULL and has a length of t_w, then forecastfun is applied successively to the subset time series y_{t-t_w+1},\dots,y_{t}, for t=\max(t_0, t_w),\dots,T.

If forward is FALSE, the last observation used for training will be y_{T-1}.

Value

A list of class c("cvforecast", "forecast") with components:

If forward is TRUE, the components mean, lower, upper, and model will also be returned, showing the information about the final fitted model and forecasts using all available observations, see e.g. forecast.ets for more details.

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Example with a rolling window
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)
print(fc)
summary(fc)

# Example with exogenous predictors
far2_xreg <- function(x, h, level, xreg, newxreg) {
  Arima(x, order=c(2, 0, 0), xreg = xreg) |>
    forecast(h = h, level = level, xreg = newxreg)
}
fc_xreg <- cvforecast(series, forecastfun = far2_xreg, h = 3, level = 95,
                      forward = TRUE, xreg = matrix(rnorm(406), ncol = 2, nrow = 203),
                      initial = 1, window = 50)

Create lags or leads of a matrix

Description

Find a shifted version of a matrix, adjusting the time base backward (lagged) or forward (leading) by a specified number of observations for each column.

Usage

lagmatrix(x, lag)

Arguments

Value

A matrix with the same class and size as x.

Examples

x <- matrix(rnorm(20), nrow = 5, ncol = 4)

# Create lags of a matrix
lagmatrix(x, c(0, 1, 2, 3))

# Create leads of a matrix
lagmatrix(x, c(0, -1, -2, -3))

Conformal PID control method

Description

Compute prediction intervals and other information by applying the conformal PID (Proportional-Integral-Derivative) control method.

Usage

pid(
  object,
  alpha = 1 - 0.01 * object$level,
  symmetric = FALSE,
  ncal = 10,
  rolling = FALSE,
  integrate = TRUE,
  scorecast = !symmetric,
  scorecastfun = NULL,
  lr = 0.1,
  Tg = NULL,
  delta = NULL,
  Csat = 2/pi * (ceiling(log(Tg) * delta) - 1/log(Tg)),
  KI = max(abs(object$errors), na.rm = TRUE),
  ...
)

Arguments

Details

The PID method combines three modules to make the final iteration:

q_{t+h|t}=\underbrace{q_{t+h-1|t-1} + \eta(\mathrm{err}_{t|t-h}-\alpha)}_{\mathrm{P}}+\underbrace{r_t\left(\sum_{i=1}^t\left(\mathrm{err}_{i|i-h}-\alpha\right)\right)}_{\mathrm{I}}+\underbrace{\hat{s}_{t+h|t}}_{\mathrm{D}}

for each individual forecast horizon h, respectively, where

• Quantile tracking part (P) is q_{t+h-1|t-1} + \eta(\mathrm{err}_{t|t-h}-\alpha), where q_{1+h|1} is set to 0 without a loss of generality, \mathrm{err}_{t|t-h}=1 if s_{t|t-h}>q_{t|t-h}, and \mathrm{err}_{t|t-h}=0 if s_{t|t-h} \leq q_{t|t-h}.

• Error integration part (I) is r_t\left(\sum_{i=1}^t\left(\mathrm{err}_{i|i-h}-\alpha\right)\right). Here we use a nonlinear saturation function r_t(x)=K_{\mathrm{I}} \tan \left(x \log (t) /\left(t C_{\text {sat }}\right)\right), where we set \tan (x)=\operatorname{sign}(x) \cdot \infty for x \notin[-\pi / 2, \pi / 2], and C_{\text {sat }}, K_{\mathrm{I}}>0 are constants that we choose heuristically.

• Scorecasting part (D) is \hat{s}_{t+h|t} is forecast generated by training a scorecaster based on nonconformity scores available at time t.

Value

A list of class c("pid", "cpforecast", "forecast") with the following components:

If mean is included in the object, the components mean, lower, and upper will also be returned, showing the information about the forecasts generated using all available observations.

References

Angelopoulos, A., Candes, E., and Tibshirani, R. J. (2024). "Conformal PID control for time series prediction", Advances in Neural Information Processing Systems, 36, 23047–23074.

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))
# Cross-validation forecasting
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)
# PID setup
Tg <- 200; delta <- 0.01
Csat <- 2 / pi * (ceiling(log(Tg) * delta) - 1 / log(Tg))
KI <- 2
lr <- 0.1
# PID without scorecaster
pidfc_nsf <- pid(fc, symmetric = FALSE, ncal = 50, rolling = TRUE,
                 integrate = TRUE, scorecast = FALSE,
                 lr = lr, KI = KI, Csat = Csat)
print(pidfc_nsf)
summary(pidfc_nsf)
# PID with a Naive model for the scorecaster
naivefun <- function(x, h) {
  naive(x) |> forecast(h = h)
}
pidfc <- pid(fc, symmetric = FALSE, ncal = 50, rolling = TRUE,
             integrate = TRUE, scorecast = TRUE, scorecastfun = naivefun,
             lr = lr, KI = KI, Csat = Csat)

Classical split conformal prediction method

Description

Compute prediction intervals and other information by applying the classical split conformal prediction (SCP) method.

Usage

scp(
  object,
  alpha = 1 - 0.01 * object$level,
  symmetric = FALSE,
  ncal = 10,
  rolling = FALSE,
  quantiletype = 1,
  weightfun = NULL,
  kess = FALSE,
  update = FALSE,
  na.rm = TRUE,
  ...
)

Arguments

Details

Consider a vector s_{t+h|t} that contains the nonconformity scores for the h-step-ahead forecasts.

If symmetric is TRUE, s_{t+h|t}=|e_{t+h|t}|. When rolling is FALSE, the (1-\alpha)-quantile \hat{q}_{t+h|t} are computed successively on expanding calibration sets s_{1+h|1},\dots,s_{t|t-h}, for t=\mathrm{ncal}+h,\dots,T. Then the prediction intervals will be [\hat{y}_{t+h|t}-\hat{q}_{t+h|t}, \hat{y}_{t+h|t}+\hat{q}_{t+h|t}]. When rolling is TRUE, the calibration sets will be of same length ncal.

If symmetric is FALSE, s_{t+h|t}^{u}=e_{t+h|t} for upper interval bounds and s_{t+h|t}^{l} = -e_{t+h|t} for lower bounds. Instead of computing (1-\alpha)-quantile, (1-\alpha/2)-quantiles for lower bound (\hat{q}_{t+h|t}^{l}) and upper bound (\hat{q}_{t+h|t}^{u}) are calculated based on their nonconformity scores, respectively. Then the prediction intervals will be [\hat{y}_{t+h|t}-\hat{q}_{t+h|t}^{l}, \hat{y}_{t+h|t}+\hat{q}_{t+h|t}^{u}].

Value

A list of class c("scp", "cvforecast", "forecast") with the following components:

If mean is included in the object, the components mean, lower, and upper will also be returned, showing the information about the test set forecasts generated using all available observations.

See Also

weighted_quantile

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# Classical conformal prediction with equal weights
scpfc <- scp(fc, symmetric = FALSE, ncal = 50, rolling = TRUE)
print(scpfc)
summary(scpfc)

# Classical conformal prediction with exponential weights
expweight <- function(n) {
  0.99^{n+1-(1:n)}
}
scpfc_exp <- scp(fc, symmetric = FALSE, ncal = 50, rolling = TRUE,
                 weightfun = expweight, kess = TRUE)

Update and repeform cross-validation forecasting and conformal prediction

Description

Update conformal prediction intervals and other information by applying the cvforecast and conformal functions.

Usage

## S3 method for class 'cpforecast'
update(object, new_data, forecastfun, new_xreg = NULL, ...)

Arguments

Value

A refreshed object of class "cpforecast" with updated fields (e.g., x, MEAN, ERROR, LOWER, UPPER, and any method-specific components), reflecting newly appended data and re-computed cross-validation forecasts and conformal prediction intervals.

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# Classical conformal prediction with equal weights
scpfc <- conformal(fc, method = "scp", symmetric = FALSE, ncal = 50, rolling = TRUE)

# Update conformal prediction using newly available data
scpfc_update <- update(scpfc, forecastfun = far2, new_data = c(1.5, 0.8, 2.3))
print(scpfc_update)
summary(scpfc_update)

Calculate interval forecast width

Description

Calculate the mean width of prediction intervals on the validation set. If window is not NULL, a matrix of the rolling means of interval width is also returned. If includemedian is TRUE, the information of the median interval width will be returned.

Usage

width(
  object,
  ...,
  level = 95,
  includemedian = FALSE,
  window = NULL,
  na.rm = FALSE
)

Arguments

Value

A list of class "width" with the following components:

Examples

# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 200, list(ar = c(0.8, -0.5)), sd = sqrt(1))

# Cross-validation forecasting with a rolling window
far2 <- function(x, h, level) {
  Arima(x, order = c(2, 0, 0)) |>
    forecast(h = h, level)
}
fc <- cvforecast(series, forecastfun = far2, h = 3, level = 95,
                 forward = TRUE, initial = 1, window = 50)

# Mean and rolling mean width for interval forecasts on validation set
wid_fc <- width(fc, level = 95, window = 50)
str(wid_fc)