Probabilistic Time Series Cross-Validation with R package crossvalidation

R-bloggers 2026-05-16

[This article was first published on T. Moudiki's Webpage - R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

A previous post introduced the crossvalidation package for R. This time, the focus is on probabilistic forecasting — evaluating not just how accurate point forecasts are, but how well-calibrated prediction intervals are, using empirical coverage rates and Winkler scores – and crossvalidation.

install.packages("remotes")install.packages("forecast")remotes::install_github("Techtonique/crossvalidation")library(crossvalidation)

Example 1

require(forecast)data("AirPassengers")eval_metric <- function(predicted, observed){  error <- observed - predicted$mean  me <- mean(error)  rmse <- sqrt(mean(error^2))  mae <- mean(abs(error))  # ----- 80% interval -----  lower80 <- predicted$lower[, 1]  upper80 <- predicted$upper[, 1]  coverage80 <- mean(    observed >= lower80 & observed <= upper80  )  alpha80 <- 0.20  winkler80 <- ifelse(    observed < lower80,    (upper80 - lower80) + (2 / alpha80) * (lower80 - observed),    ifelse(      observed > upper80,      (upper80 - lower80) + (2 / alpha80) * (observed - upper80),      (upper80 - lower80)    )  )  # ----- 95% interval -----  lower95 <- predicted$lower[, 2]  upper95 <- predicted$upper[, 2]  coverage95 <- mean(    observed >= lower95 & observed <= upper95  )  alpha95 <- 0.05  winkler95 <- ifelse(    observed < lower95,    (upper95 - lower95) + (2 / alpha95) * (lower95 - observed),    ifelse(      observed > upper95,      (upper95 - lower95) + (2 / alpha95) * (observed - upper95),      (upper95 - lower95)    )  )  c(    ME = me,    RMSE = rmse,    MAE = mae,    Coverage80 = coverage80,    Winkler80 = mean(winkler80),    Coverage95 = coverage95,    Winkler95 = mean(winkler95)  )}(res <- crossval_ts(y=AirPassengers, initial_window = 10,horizon = 3, fcast_func = forecast::thetaf, eval_metric = eval_metric))print(colMeans(res))Loading required package: forecast  |======================================================================| 100%

A matrix: 132 × 7 of type dblMERMSEMAECoverage80Winkler80Coverage95Winkler95result.1-28.79466029.30028728.7946600.0000000153.109920.3333333207.58384result.2 16.19852616.89430216.1985261.0000000 45.017951.0000000 68.84902result.3 11.20149415.99335912.5782761.0000000 45.059961.0000000 68.91326result.4 21.43012522.48389521.4301250.6666667 63.012071.0000000 68.84778result.5 10.05576511.52774610.0557651.0000000 45.999671.0000000 70.35043result.6 -2.64082210.676714 9.9994661.0000000 46.569071.0000000 71.22125result.7 14.29643423.70913220.5311350.6666667 75.041861.0000000 67.58381result.8 38.24749739.52999838.2474970.0000000198.749900.3333333212.44029result.9 23.04315923.94763023.0431590.3333333 93.834631.0000000 64.19366result.10-21.68906727.90756021.6890670.6666667 90.233771.0000000 84.12361result.11-41.78215746.66419941.7821570.3333333222.063100.3333333345.16553result.12-34.93483136.51208134.9348310.3333333162.380920.6666667212.58117result.13 -4.00270012.728771 9.9991001.0000000 59.644751.0000000 91.21878result.14 30.34958230.58876130.3495820.6666667 72.143551.0000000 99.76932result.15 21.19234925.80671221.1923490.6666667 71.390941.0000000101.02401result.16 23.19314325.91487523.1931430.6666667 91.706601.0000000 76.57925result.17 30.08154230.67996030.0815420.3333333111.586891.0000000 75.78459result.18 -6.530509 9.111376 6.9990591.0000000 69.517041.0000000106.31714result.19 19.90758623.01076219.9075861.0000000 67.035061.0000000102.52128result.20 17.63108919.82935517.6310891.0000000 67.975731.0000000103.95991result.21 11.73802214.71818512.2298461.0000000 61.616171.0000000 94.23380result.22-21.78749028.48950921.7874900.6666667 93.309201.0000000 88.70090result.23-43.55757147.52724443.5575710.3333333206.773680.6666667216.50078result.24-34.47355835.51415534.4735580.3333333146.632880.6666667173.17046result.25 -4.69936010.498595 7.2012241.0000000 60.075501.0000000 91.87755result.26 25.97413826.58127225.9741381.0000000 63.019421.0000000 96.37989result.27 16.90510919.47460016.9051091.0000000 58.044721.0000000 88.77173result.28 15.21876016.35291715.2187601.0000000 55.277211.0000000 84.53920result.29 7.625241 8.933828 7.6252411.0000000 55.277181.0000000 84.53916result.30 2.26197017.59532615.6662121.0000000 57.132921.0000000 87.37725⋮⋮⋮⋮⋮⋮⋮⋮result.103 95.047754111.26440 95.047750.3333333 485.70960.6666667 594.3549result.104 121.335201125.76554121.335200.0000000 646.57500.3333333 772.4818result.105 27.661546 53.66952 52.335670.6666667 149.46691.0000000 226.7499result.106 -82.928463106.53675 87.398380.3333333 439.04760.6666667 391.0034result.107-168.429957174.86402168.429960.00000001125.85340.00000002680.3671result.108 -86.047368 89.34969 86.047370.6666667 241.50861.0000000 281.3325result.109 -35.392983 38.64620 35.392981.0000000 192.33141.0000000 294.1455result.110 32.273683 33.69167 32.273681.0000000 199.99781.0000000 305.8702result.111 35.911969 45.52857 35.911971.0000000 195.20691.0000000 298.5432result.112 28.584481 41.79144 38.166541.0000000 196.54091.0000000 300.5833result.113 78.144295 79.31310 78.144301.0000000 196.93431.0000000 301.1850result.114 37.152546 52.61404 39.210441.0000000 192.54871.0000000 294.4778result.115 95.078342110.88602 95.078340.6666667 366.36761.0000000 274.9151result.116 109.166178116.17612109.166180.3333333 406.73971.0000000 277.4405result.117 41.289554 62.02085 57.334900.3333333 215.15771.0000000 222.4127result.118 -92.399494116.61777 92.824070.3333333 466.72850.6666667 445.3571result.119-175.618445183.27955175.618450.00000001143.54790.00000002574.2409result.120 -94.580461 97.36039 94.580460.6666667 277.78471.0000000 293.2590result.121 -27.751828 32.93559 27.751831.0000000 202.13741.0000000 309.1425result.122 36.177008 38.16646 36.177011.0000000 208.63521.0000000 319.0800result.123 5.992278 14.16185 13.997431.0000000 200.00981.0000000 305.8885result.124 12.637863 33.65269 27.980301.0000000 200.18281.0000000 306.1532result.125 71.834372 76.95073 71.834371.0000000 200.57531.0000000 306.7534result.126 85.518711 93.75094 85.518710.6666667 252.56381.0000000 295.0496result.127 94.429064115.52397 94.429060.6666667 407.36360.6666667 417.2566result.128 173.325805177.66652173.325800.00000001129.61410.00000002547.8618result.129 33.890665 63.84191 61.668610.6666667 242.69011.0000000 230.3885result.130-119.059067137.73685119.059070.3333333 619.41660.3333333 668.9786result.131-180.821172190.45241180.821170.00000001152.49490.00000002469.3936result.132-103.156396108.61881103.156400.6666667 330.04001.0000000 302.1675

         ME        RMSE         MAE  Coverage80   Winkler80  Coverage95   2.6570822  51.4271704  46.5118747   0.6590909 218.4527816   0.8459596   Winkler95 312.1383104

Example 2

eval_metric <- function(predicted, observed){  error <- observed - predicted$mean  me <- mean(error)  rmse <- sqrt(mean(error^2))  mae <- mean(abs(error))  # Only one interval returned  lower <- predicted$lower  upper <- predicted$upper  coverage <- mean(    observed >= lower & observed <= upper  )  alpha <- 0.05  winkler <- ifelse(    observed < lower,    (upper - lower) + (2 / alpha) * (lower - observed),    ifelse(      observed > upper,      (upper - lower) + (2 / alpha) * (observed - upper),      (upper - lower)    )  )  c(    ME = me,    RMSE = rmse,    MAE = mae,    Coverage95 = coverage,    Winkler95 = mean(winkler)  )}fcast_func <- function(y, h, ...){  forecast::thetaf(    y,    h = h,    level = 95  )}res <- crossval_ts(  y = AirPassengers,  initial_window = 10,  horizon = 3,  fcast_func = fcast_func,  eval_metric = eval_metric)print(colMeans(res))  |======================================================================| 100%         ME        RMSE         MAE  Coverage95   Winkler95   2.6570822  51.4271704  46.5118747   0.8459596 312.1383104 boxplot(res[, "Coverage95"])

image-title-here

To leave a comment for the author, please follow the link and comment on their blog: T. Moudiki's Webpage - R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: Probabilistic Time Series Cross-Validation with R package crossvalidation