inst/doc/maxbootR-intro.R

In maxbootR: Efficient Bootstrap Methods for Block Maxima

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)
library(maxbootR)
library(ggplot2)
library(tidyr)
library(dplyr)
library(lubridate)
set.seed(42)

## -----------------------------------------------------------------------------
head(logret_data)
tail(logret_data)
help("logret_data")

## ----fig.alt= "Raw Data Plot"-------------------------------------------------
logret_data %>% 
  ggplot(aes(x = day, y = neg_log_ret)) +
  geom_line(color = "steelblue")

length(logret_data$day) / 30  # approx. number of years
sum(is.na(logret_data$neg_log_ret))  # number of missing values

## ----fig.alt = "Block Maxima of Negative Log-Returns"-------------------------
bsize <- 250
bm_db <- blockmax(logret_data$neg_log_ret, block_size = bsize, type = "db")
bm_sb <- blockmax(logret_data$neg_log_ret, block_size = bsize, type = "sb")

# Time vector per block type
day_db <- logret_data$day[seq(1, length(bm_db) * bsize, by = bsize)]
day_sb <- logret_data$day[1:length(bm_sb)]

# Combine into tidy tibble
df_db <- tibble(day = day_db, value = bm_db, method = "Disjoint Blocks")
df_sb <- tibble(day = day_sb, value = bm_sb, method = "Sliding Blocks")
df_all <- bind_rows(df_db, df_sb)

# Plot
ggplot(df_all, aes(x = day, y = value)) +
  geom_line(color = "steelblue") +
  facet_wrap(~ method, nrow = 1) +
  labs(title = "Block Maxima of Negative Log-Returns",
       x = "Date", y = "Block Maximum")

## -----------------------------------------------------------------------------
bst.bm_db_gev <- maxbootr(
  xx = logret_data$neg_log_ret, est = "gev", block_size = 250, 
  B = 1000, type = "db"
) 
summary(bst.bm_db_gev[, 3])

bst.bm_sb_gev <- maxbootr(
  xx = logret_data$neg_log_ret, est = "gev", block_size = 250, 
  B = 1000, type = "sb"
)
summary(bst.bm_sb_gev[, 3])

## -----------------------------------------------------------------------------
bst.bm_db_q <- maxbootr(
  xx = logret_data$neg_log_ret, est = "quantile", block_size = 250, 
  B = 1000, type = "db", p = 0.99
) 
summary(bst.bm_db_q)

bst.bm_sb_q <- maxbootr(
  xx = logret_data$neg_log_ret, est = "quantile", block_size = 250, 
  B = 1000, type = "sb", p = 0.99
)
summary(bst.bm_sb_q)

## ----fig.alt= "Bootstrap Estimates of Extreme Quantile"-----------------------
# Trim upper 2% of bootstrap replicates
bst.bm_db_q_trimmed <- bst.bm_db_q[bst.bm_db_q < quantile(bst.bm_db_q, 0.98)]
bst.bm_sb_q_trimmed <- bst.bm_sb_q[bst.bm_sb_q < quantile(bst.bm_sb_q, 0.98)]

# Combine for plotting
df_q <- tibble(
  value = c(bst.bm_db_q_trimmed, bst.bm_sb_q_trimmed),
  method = c(rep("Disjoint Blocks", length(bst.bm_db_q_trimmed)),
             rep("Sliding Blocks", length(bst.bm_sb_q_trimmed)))
)

# Histogram plot
ggplot(df_q, aes(x = value)) +
  geom_histogram(fill = "steelblue", color = "white", bins = 30) +
  facet_wrap(~ method, nrow = 1) +
  labs(
    title = "Bootstrap Estimates of Extreme Quantile",
    x = "Estimated Quantile",
    y = "Count"
  )

## -----------------------------------------------------------------------------
# Variance ratio
var(bst.bm_sb_q_trimmed) / var(bst.bm_db_q_trimmed)

## ----fig.alt="Daily Negative Log-Returns with Extreme Quantile"---------------
q99 <- quantile(bst.bm_sb_q_trimmed, 0.5)

ggplot(logret_data, aes(x = day, y = neg_log_ret)) +
  geom_line(color = "steelblue") +
  geom_hline(yintercept = q99, color = "red", linetype = "dashed") +
  labs(
    title = "Daily Negative Log-Returns with Extreme Quantile",
    x = "Date",
    y = "Negative Log-Return"
  )

## -----------------------------------------------------------------------------
head(temp_data)
tail(temp_data)
help("temp_data")

## ----fig.alt="3 Years of Daily Temperature"-----------------------------------
temp_data %>% 
  filter(lubridate::year(day) %in% c(1900, 1901, 1902)) %>% 
  ggplot(aes(x = day, y = temp)) +
  geom_line(color = "steelblue")

## -----------------------------------------------------------------------------
temp_data_cl <- temp_data %>% 
  filter(lubridate::month(day) %in% c(6, 7, 8))

## ----fig.alt="Block Maxima of Summer Temperatures"----------------------------
bsize <- 92
bm_db_temp <- blockmax(temp_data_cl$temp, block_size = bsize, type = "db")
bm_sb_temp <- blockmax(temp_data_cl$temp, block_size = bsize, type = "sb")

# Create time vectors for plotting
day_db_temp <- temp_data_cl$day[seq(1, length(bm_db_temp) * bsize, by = bsize)]
day_sb_temp <- temp_data_cl$day[1:length(bm_sb_temp)]

# Create tidy tibble for plotting
df_db_temp <- tibble(day = day_db_temp, value = bm_db_temp, method = "Disjoint Blocks")
df_sb_temp <- tibble(day = day_sb_temp, value = bm_sb_temp, method = "Sliding Blocks")
df_all_temp <- bind_rows(df_db_temp, df_sb_temp)

# Plot block maxima
ggplot(df_all_temp, aes(x = day, y = value)) +
  geom_line(color = "steelblue") +
  facet_wrap(~ method, nrow = 1) +
  labs(title = "Block Maxima of Summer Temperatures",
       x = "Date", y = "Block Maximum")

## -----------------------------------------------------------------------------
bst.bm_db_temp_q <- maxbootr(
  xx = temp_data_cl$temp, est = "rl", block_size = bsize, 
  B = 1000, type = "db", annuity = 100
)
summary(bst.bm_db_temp_q)

bst.bm_sb_temp_q <- maxbootr(
  xx = temp_data_cl$temp, est = "rl", block_size = bsize, 
  B = 1000, type = "sb", annuity = 100
)
summary(bst.bm_sb_temp_q)

## ----fig.alt="Bootstrap Estimates of 100-Year Return Level"-------------------
# Combine for plotting
df_q_temp <- tibble(
  value = c(bst.bm_db_temp_q, bst.bm_sb_temp_q),
  method = c(rep("Disjoint Blocks", length(bst.bm_db_temp_q)),
             rep("Sliding Blocks", length(bst.bm_sb_temp_q)))
)

# Histogram plot
ggplot(df_q_temp, aes(x = value)) +
  geom_histogram(fill = "steelblue", color = "white", bins = 30) +
  facet_wrap(~ method, nrow = 1) +
  labs(
    title = "Bootstrap Estimates of 100-Year Return Level",
    x = "Estimated Return Level",
    y = "Count"
  )

## -----------------------------------------------------------------------------
# Compute and display variance ratio
var(bst.bm_sb_temp_q) / var(bst.bm_db_temp_q)

## ----fig.alt="All Temperatures with Estimated 100-Year Return Level"----------
rl <- quantile(bst.bm_sb_temp_q, 0.5)

ggplot(temp_data, aes(x = day, y = temp)) +
  geom_line(color = "steelblue") +
  geom_hline(yintercept = rl, color = "red", linetype = "dashed") +
  labs(
    title = "All Temperatures with Estimated 100-Year Return Level",
    x = "Date",
    y = "Daily Max Temperature"
  )