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
"cvforecast"
. It must have an argumentx
for original univariate time series, an argumentMEAN
for point forecasts andERROR
for forecast errors on validation set. See the results of a call tocvforecast
.- alpha
A numeric vector of significance levels to achieve a desired coverage level \(1-\alpha\).
- symmetric
If
TRUE
, symmetric nonconformity scores (i.e. \(|e_{t+h|t}|\)) are used. IfFALSE
, asymmetric nonconformity scores (i.e. \(e_{t+h|t}\)) are used, and then upper bounds and lower bounds are produced separately.- ncal
Length of the calibration set. If
rolling = FALSE
, it denotes the initial period of calibration sets. Otherwise, it indicates the period of every rolling calibration set.- rolling
If
TRUE
, a rolling window strategy will be adopted to form the calibration set. Otherwise, expanding window strategy will be used.- 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
weighted_quantile
function in the ggdist package.- 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
NULL
, equal weights will be used for sample quantile computation. Currently, only non-data-dependent weights are supported.- kess
If
TRUE
, Kish's effective sample size is used for sample quantile computation.- update
If
TRUE
, the function will be compatible with theupdate
(update.cpforecast) function, allowing for easy updates of conformal prediction.- na.rm
If
TRUE
, corresponding entries in sample values and weights are removed if either isNA
when calculating sample quantile.- ...
Other arguments are passed to
weightfun
.
Value
A list of class c("scp", "cvforecast", "forecast")
with the following components:
- x
The original time series.
- series
The name of the series
x
.- 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 \(h\)th column holds the point forecasts for forecast horizon \(h\). The time index corresponds to the period for which the forecast is produced.
- ERROR
Forecast errors given by \(e_{t+h|t} = y_{t+h}-\hat{y}_{t+h|t}\).
- LOWER
A list containing lower bounds for prediction intervals for each
level
. Each element within the list will be a multivariate time series with the same dimensional characteristics asMEAN
.- UPPER
A list containing upper bounds for prediction intervals for each
level
. Each element within the list will be a multivariate time series with the same dimensional characteristics asMEAN
.- level
The confidence values associated with the prediction intervals.
- call
The matched call.
- model
A list containing detailed information about the
cvforecast
andconformal
models.
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.
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}]\).
Examples
# Simulate time series from an AR(2) model
library(forecast)
series <- arima.sim(n = 1000, 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 = c(80, 95),
forward = TRUE, initial = 1, window = 100)
# Classical conformal prediction with equal weights
scpfc <- scp(fc, symmetric = FALSE, ncal = 100, rolling = TRUE)
print(scpfc)
#> SCP
#>
#> Call:
#> scp(object = fc, symmetric = FALSE, ncal = 100, rolling = TRUE)
#>
#> cp_times = 799 (the forward step included)
#>
#> Forecasts of the forward step:
#> Cross-validation
#>
#> Call:
#> scp(object = fc, symmetric = FALSE, ncal = 100, rolling = TRUE)
#>
#> fit_times = (the forward step included)
#>
#> Forecasts of the forward step:
#> Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
#> 1001 1.2540418 -0.1490929 2.520996 -1.065434 3.817075
#> 1002 0.3784289 -1.2282025 2.186605 -2.834616 3.380624
#> 1003 -0.3086171 -1.8457568 1.494619 -3.414272 2.618403
summary(scpfc)
#> SCP
#>
#> Call:
#> scp(object = fc, symmetric = FALSE, ncal = 100, rolling = TRUE)
#>
#> cp_times = 799 (the forward step included)
#>
#> Forecasts of the forward step:
#> Cross-validation
#>
#> Call:
#> scp(object = fc, symmetric = FALSE, ncal = 100, rolling = TRUE)
#>
#> fit_times = (the forward step included)
#>
#> Forecasts of the forward step:
#> Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
#> 1001 1.2540418 -0.1490929 2.520996 -1.065434 3.817075
#> 1002 0.3784289 -1.2282025 2.186605 -2.834616 3.380624
#> 1003 -0.3086171 -1.8457568 1.494619 -3.414272 2.618403
#>
#> Cross-validation error measures:
#> ME MAE MSE RMSE MPE MAPE MASE RMSSE Winkler_95 MSIS_95
#> CV -0.011 0.958 1.46 1.08 42.611 228.563 0.934 0.837 5.944 5.731
# Classical conformal prediction with exponential weights
expweight <- function(n) {
0.99^{n+1-(1:n)}
}
scpfc_exp <- scp(fc, symmetric = FALSE, ncal = 100, rolling = TRUE,
weightfun = expweight, kess = TRUE)