## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 7,
  fig.height = 4
)

## ----setup, message = FALSE---------------------------------------------------
library(forecastdom)
data(quaedvlieg2021)
str(quaedvlieg2021, max.level = 1)

## ----mean-table---------------------------------------------------------------
H <- ncol(quaedvlieg2021$uspa)

means <- data.frame(
  h    = seq_len(H),
  uspa = colMeans(quaedvlieg2021$uspa),
  aspa = colMeans(quaedvlieg2021$aspa)
)

knitr::kable(
  means, digits = 3, row.names = FALSE,
  col.names = c("$h$", "uSPA dataset", "aSPA dataset"))

## ----mean-plot, fig.height = 3.5, fig.alt = "Mean loss differential by horizon for the uspa and aspa example matrices"----
library(ggplot2)

df <- data.frame(
  horizon = rep(seq_len(H), 2),
  d_bar   = c(colMeans(quaedvlieg2021$uspa),
              colMeans(quaedvlieg2021$aspa)),
  dataset = rep(c("uspa", "aspa"), each = H)
)

ggplot(df, aes(horizon, d_bar, colour = dataset)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey60") +
  geom_line() + geom_point() +
  labs(x = "Horizon h", y = expression(bar(d)[h]),
       title = "Mean loss differential by horizon") +
  theme_minimal()

## ----uspa, results = "hold"---------------------------------------------------
set.seed(1)
uspa_uspa <- uspa_mh_test(quaedvlieg2021$uspa, L = 3, B = 999)

set.seed(1)
uspa_aspa <- uspa_mh_test(quaedvlieg2021$aspa, L = 3, B = 999)

uspa_uspa
uspa_aspa

## ----aspa-uniform, results = "hold"-------------------------------------------
w_unif <- rep(1 / H, H)

set.seed(1)
aspa_uspa <- aspa_mh_test(quaedvlieg2021$uspa, weights = w_unif, L = 3, B = 999)

set.seed(1)
aspa_aspa <- aspa_mh_test(quaedvlieg2021$aspa, weights = w_unif, L = 3, B = 999)

aspa_uspa
aspa_aspa

## ----aspa-down, results = "hold"----------------------------------------------
w_down <- c(rep(0, 4), rep(1, 16)) / 16 # zero weight on h = 1..4

set.seed(1)
aspa_mh_test(quaedvlieg2021$aspa, weights = w_down, L = 3, B = 999)

## ----aspa-up, results = "hold"------------------------------------------------
w_up <- c(rep(4, 4), rep(0, 16)) / 16 # all weight on h = 1..4

set.seed(1)
aspa_mh_test(quaedvlieg2021$aspa, weights = w_up, L = 3, B = 999)

## ----L-sensitivity------------------------------------------------------------
Ls <- c(2, 3, 5, 8, 12)

sens <- do.call(rbind, lapply(Ls, function(L) {

  set.seed(1)
  u <- uspa_mh_test(quaedvlieg2021$uspa, L = L, B = 499)
  
  set.seed(1)
  a <- aspa_mh_test(quaedvlieg2021$uspa, weights = w_unif, L = L, B = 499)
  
  data.frame(L = L,
             uspa_stat = u$statistic, uspa_p = u$pvalue,
             aspa_stat = a$statistic, aspa_p = a$pvalue)

}))

knitr::kable(
  sens, digits = 3, row.names = FALSE,
  col.names = c("$L$",
                "uSPA stat", "uSPA $p$",
                "aSPA stat", "aSPA $p$"))