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 |
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:
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
resid |
A numeric vector of residuals from either the validation or test data. |
na.rm |
If |
... |
Additional arguments for each measure. |
actual |
A numeric vector of responses matching the forecasts (for percentage measures). |
train |
A numeric vector of responses used to train the model (for scaled measures). |
demean |
Should the response be demeaned (for MASE and RMSSE)? |
period |
The seasonal period of the data. |
d |
Should the response model include a first difference? |
D |
Should the response model include a seasonal difference? |
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
lower |
A numeric vector of lower bounds of interval forecasts. |
upper |
A numeric vector of upper bounds of interval forecasts. |
actual |
A numeric vector of realised values. |
train |
A numeric vector of responses used to train the model (for scaled scores). |
level |
The nominal level of the forecast interval (e.g., 95 or 0.95). |
period |
The seasonal period of the data. |
d |
Should the response model include a first difference? |
D |
Should the response model include a seasonal difference? |
na.rm |
If |
... |
Additional arguments for each measure. |
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
object |
An object of class |
x |
An optional numerical vector containing actual values of the same
length as |
CV |
If |
period |
The seasonal period of the data. |
measures |
A list of accuracy measure functions to compute (such as point_measures or interval_measures). |
byhorizon |
If |
... |
Additional arguments depending on the specific measure. |
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
object |
An object of class |
alpha |
A numeric vector of significance levels to achieve a desired
coverage level |
ncal |
Length of the burn-in period for training the scorecaster.
If |
rolling |
If |
integrate |
If |
scorecast |
If |
lr |
Initial learning rate used for quantile tracking. |
Tg |
The time that is set to achieve the target absolute coverage guarantee before this. |
delta |
The target absolute coverage guarantee is set to |
Csat |
A positive constant ensuring that by time |
KI |
A positive constant to place the integrator on the same scale as the scores. |
... |
Other arguments are passed to the function. |
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:
x |
The original time series. |
series |
The name of the series |
method |
A character string "acmcp". |
cp_times |
The number of times the conformal prediction is performed in cross-validation. |
MEAN |
Point forecasts as a multivariate time series, where the |
ERROR |
Forecast errors given by
|
LOWER |
A list containing lower bounds for prediction intervals for
each |
UPPER |
A list containing upper bounds for prediction intervals for
each |
level |
The confidence values associated with the prediction intervals. |
call |
The matched call. |
model |
A list containing information abouth the conformal prediction model. |
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
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
object |
An object of class |
alpha |
A numeric vector of significance levels to achieve a desired
coverage level |
gamma |
The step size parameter |
symmetric |
If |
ncal |
Length of the calibration set. If |
rolling |
If |
quantiletype |
An integer between 1 and 9 determining the type of
quantile estimator to be used. Types 1 to 3 are for discontinuous quantiles,
types 4 to 9 are for continuous quantiles. See the
|
update |
If |
na.rm |
If |
... |
Other arguments are passed to the
|
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:
x |
The original time series. |
series |
The name of the series |
method |
A character string "acp". |
cp_times |
The number of times the conformal prediction is performed in cross-validation. |
MEAN |
Point forecasts as a multivariate time series, where the |
ERROR |
Forecast errors given by
|
LOWER |
A list containing lower bounds for prediction intervals for
each |
UPPER |
A list containing upper bounds for prediction intervals for
each |
level |
The confidence values associated with the prediction intervals. |
call |
The matched call. |
model |
A list containing information abouth the conformal prediction model. |
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
object |
An object of class |
... |
Additional arguments to be passed to the selected conformal method. |
method |
A character string specifying the conformal method to be applied.
Possible options include |
Value
An object whose class depends on the method invoked.
See Also
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
object |
An object of class |
... |
Additional inputs if |
level |
Target confidence level for prediction intervals. |
window |
If not |
na.rm |
A logical indicating whether |
Value
A list of class "coverage"
with the following components:
mean |
Mean coverage across the validation set. |
ifinn |
A indicator matrix as a multivariate time series, where the |
rollmean |
If |
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
y |
Univariate time series. |
forecastfun |
Function to return an object of class |
h |
Forecast horizon. |
level |
Confidence level for prediction intervals.
If |
forward |
If |
xreg |
Exogenous predictor variables passed to |
initial |
Initial period of the time series where no cross-validation forecasting is performed. |
window |
Length of the rolling window. If |
... |
Other arguments are passed to |
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:
x |
The original time series. |
series |
The name of the series |
xreg |
Exogenous predictor variables used in the model, if applicable. |
method |
A character string "cvforecast". |
fit_times |
The number of times the model is fitted in cross-validation. |
MEAN |
Point forecasts as a multivariate time series, where the |
ERROR |
Forecast errors given by
|
LOWER |
A list containing lower bounds for prediction intervals for
each |
UPPER |
A list containing upper bounds for prediction intervals for
each |
level |
The confidence values associated with the prediction intervals. |
call |
The matched call. |
forward |
Whether |
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
x |
A matrix or multivariate time series. |
lag |
A vector of lags (positive values) or leads (negative values) with
a length equal to the number of columns of |
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
object |
An object of class |
alpha |
A numeric vector of significance levels to achieve a desired
coverage level |
symmetric |
If |
ncal |
Length of the burn-in period for training the scorecaster.
If |
rolling |
If |
integrate |
If |
scorecast |
If |
scorecastfun |
A scorecaster function to return an object of class
|
lr |
Initial learning rate used for quantile tracking. |
Tg |
The time that is set to achieve the target absolute coverage guarantee before this. |
delta |
The target absolute coverage guarantee is set to |
Csat |
A positive constant ensuring that by time |
KI |
A positive constant to place the integrator on the same scale as the scores. |
... |
Other arguments are passed to the |
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)
, whereq_{1+h|1}
is set to 0 without a loss of generality,\mathrm{err}_{t|t-h}=1
ifs_{t|t-h}>q_{t|t-h}
, and\mathrm{err}_{t|t-h}=0
ifs_{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 functionr_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
forx \notin[-\pi / 2, \pi / 2]
, andC_{\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 timet
.
Value
A list of class c("pid", "cpforecast", "forecast")
with the following components:
x |
The original time series. |
series |
The name of the series |
method |
A character string "pid". |
cp_times |
The number of times the conformal prediction is performed in cross-validation. |
MEAN |
Point forecasts as a multivariate time series, where the |
ERROR |
Forecast errors given by
|
LOWER |
A list containing lower bounds for prediction intervals for
each |
UPPER |
A list containing upper bounds for prediction intervals for
each |
level |
The confidence values associated with the prediction intervals. |
call |
The matched call. |
model |
A list containing information abouth the conformal prediction model. |
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
object |
An object of class |
alpha |
A numeric vector of significance levels to achieve a desired
coverage level |
symmetric |
If |
ncal |
Length of the calibration set. If |
rolling |
If |
quantiletype |
An integer between 1 and 9 determining the type of
quantile estimator to be used. Types 1 to 3 are for discontinuous quantiles,
types 4 to 9 are for continuous quantiles. See the
|
weightfun |
Function to return a vector of weights used for sample quantile
computation. Its first argument must be an integer indicating the number of
observations for which weights are generated. If |
kess |
If |
update |
If |
na.rm |
If |
... |
Other arguments are passed to |
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:
x |
The original time series. |
series |
The name of the series |
xreg |
Exogenous predictor variables used, if applicable. |
method |
A character string "scp". |
cp_times |
The number of times the conformal prediction is performed in cross-validation. |
MEAN |
Point forecasts as a multivariate time series, where the |
ERROR |
Forecast errors given by
|
LOWER |
A list containing lower bounds for prediction intervals for
each |
UPPER |
A list containing upper bounds for prediction intervals for
each |
level |
The confidence values associated with the prediction intervals. |
call |
The matched call. |
model |
A list containing detailed information about the |
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
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
object |
An object of class |
new_data |
A vector of newly available data. |
forecastfun |
Function to return an object of class |
new_xreg |
Newly available exogenous predictor variables passed to
|
... |
Other arguments are passed to |
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
object |
An object of class |
... |
Additional inputs if |
level |
Target confidence level for prediction intervals. |
includemedian |
If |
window |
If not |
na.rm |
A logical indicating whether |
Value
A list of class "width"
with the following components:
width |
Forecast interval width as a multivariate time series, where the |
mean |
Mean interval width across the validation set. |
rollmean |
If |
median |
Median interval width across the validation set. |
rollmedian |
If |
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)