Worked example: time-series workflow

This example fits a structured time-series kernel, labels interpolation and forecast points, and evaluates the model with rolling-origin backtesting. The code below is the script inst/examples/time-series-workflow.R, installed as system.file("examples", "time-series-workflow.R", package = "gaussianprocesses"); it is run as-is to produce this page.

# Time-series Gaussian-process workflow
#
# This example keeps the model structure explicit:
# - a linear covariance component for nonstationary trend;
# - a periodic component;
# - a local Matérn component.
#
# The forecast helper labels interpolation and extrapolation separately.

time <- seq(0, 24, by = 1)
response <- 0.04 * time +
  sin(2 * pi * time / 12) +
  0.15 * cos(2 * pi * time / 3)

kernel <- time_series_kernel(
  period = 12,
  trend_variance = 0.05,
  periodic_variance = 1,
  periodic_length_scale = 0.8,
  local_variance = 0.2,
  local_length_scale = 2
)

model <- fit_time_series_gp(
  time,
  response,
  kernel = kernel,
  noise_variance = 0.05
)

# Includes both interpolation (inside the training range) and extrapolation.
requested_time <- c(6.5, 18.5, 25, 26, 27)

forecast <- forecast_gp(
  model,
  requested_time,
  interval_level = 0.95
)

data.frame(
  time = requested_time,
  region = forecast$region,
  mean = forecast$mean,
  lower = forecast$prediction_interval[, "lower"],
  upper = forecast$prediction_interval[, "upper"]
)
#>   time        region       mean       lower     upper
#> 1  6.5 interpolation 0.05390508 -0.51737637 0.6251865
#> 2 18.5 interpolation 0.53260020 -0.03881691 1.1040173
#> 3 25.0      forecast 1.44644238  0.51979753 2.3730872
#> 4 26.0      forecast 1.79532710  0.60737685 2.9832774
#> 5 27.0      forecast 2.06062944  0.75285126 3.3684076

backtest <- rolling_origin_gp(
  time,
  response,
  kernel = kernel,
  initial_window = 12,
  horizon = 3,
  step = 3,
  noise_variance = 0.05
)

forecast_metrics(backtest)
#> $n_forecasts
#> [1] 12
#> 
#> $rmse
#> [1] 0.2913891
#> 
#> $mae
#> [1] 0.21833
#> 
#> $standardized_rmse
#> [1] 0.3849505
#> 
#> $empirical_coverage
#> [1] 1
#> 
#> $interval_level
#> [1] 0.95
#> 
#> $mean_interval_width
#> [1] 2.69768
#> 
#> $mean_log_predictive_density
#> [1] -0.60062