## ----include = FALSE---------------------------------------------------------- NOT_CRAN <- identical(Sys.getenv("NOT_CRAN"), "true") knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = NOT_CRAN, echo = TRUE ) ## ----setup-------------------------------------------------------------------- # library(ExactVaRTest) ## ----message=FALSE, warning=FALSE--------------------------------------------- # library(bench) # library(dplyr) # library(tidyr) # library(purrr) # library(knitr) # # n_vec <- c(50, 100, 250, 500, 750, 1000) # alpha_vec <- c(0.01, 0.025, 0.05) # # grid <- expand.grid(n = n_vec, alpha = alpha_vec, KEEP.OUT.ATTRS = FALSE) # # bench_one <- function(n, alpha, fun) { # bm <- bench::mark( # fun(n = n, alpha = alpha), # iterations = 3, # check = FALSE # ) # tibble( # n = n, # alpha = alpha, # median_s = as.numeric(bm$median) # ) # } # # timings_ind <- pmap_dfr(grid, bench_one, fun = lr_ind_dist) # # timings_cc <- pmap_dfr(grid, bench_one, fun = lr_cc_dist) # # fmt_wide <- function(df) { # df %>% # mutate(alpha = sprintf("alpha = %.3f", alpha)) %>% # pivot_wider(names_from = alpha, values_from = median_s) %>% # arrange(n) # } # # table_ind <- fmt_wide(timings_ind) # table_cc <- fmt_wide(timings_cc) # # kable( # table_ind, # digits = 3, # col.names = c("n", "α = 0.010", "α = 0.025", "α = 0.050"), # caption = "Median run-time (seconds) for lr_ind_dist()" # ) # # kable( # table_cc, # digits = 3, # col.names = c("n", "α = 0.010", "α = 0.025", "α = 0.050"), # caption = "Median run-time (seconds) for lr_cc_dist()" # ) ## ----message=FALSE, warning=FALSE--------------------------------------------- # library(xts) # # alpha <- 0.01 # window <- 250 # # local_file <- "inst/extdata/GSPC_close.rds" # file_path <- if (file.exists(local_file)) local_file else # system.file("extdata", "GSPC_close.rds", package = "ExactVaRTest") # # px <- readRDS(file_path) # ret <- diff(log(px), lag = 1) # # VaR <- rollapply( # ret, window, # function(r) quantile(r, alpha, na.rm = TRUE), # by.column = FALSE, align = "right" # ) # # keep <- complete.cases(ret, VaR) # ret <- ret[keep] # VaR <- coredata(VaR[keep]) # # x <- ifelse(coredata(ret) < VaR, 1L, 0L) # # cat("Series length :", length(x), "\n") # cat("Exception rate:", mean(x), "\n\n") # # bench_res <- bench::mark( # LR_ind = backtest_lr(x, alpha, "ind"), # LR_cc = backtest_lr(x, alpha, "cc"), # iterations = 10, # check = FALSE # ) # # suppressWarnings( # print(bench_res[, c("expression", "median")]) # ) # # res_ind <- backtest_lr(x, alpha, "ind") # res_cc <- backtest_lr(x, alpha, "cc") # # cat("\n--- Independence test ---\n") # print(res_ind) # # cat("\n--- Conditional-coverage test ---\n") # print(res_cc) ## ----message=FALSE, warning=FALSE--------------------------------------------- # n <- 250 # alpha <- 0.01 # # d_ind <- lr_ind_dist(n, alpha) # d_cc <- lr_cc_dist(n, alpha) # d_cc$LR <- d_cc$LR_cc # d_cc$prob <- d_cc$prob_cc # # oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1)) # # # LR_ind vs χ²(1) # curve(pchisq(x, df = 1), 0, 20, lty = 2, ylab = "CDF", # xlab = "LR_ind statistic", main = "LR_ind") # lines(stepfun(d_ind$LR, c(0, cumsum(d_ind$prob))), pch = 19) # legend("bottomright", c("Chi-sq(1)", "Exact"), lty = c(2, 1), bty = "n") # # # LR_cc vs χ²(2) # curve(pchisq(x, df = 2), 0, 30, lty = 2, ylab = "CDF", # xlab = "LR_cc statistic", main = "LR_cc") # lines(stepfun(d_cc$LR, c(0, cumsum(d_cc$prob))), pch = 19) # legend("bottomright", c("Chi-sq(2)", "Exact"), lty = c(2, 1), bty = "n") # # par(oldpar) ## ----message=FALSE, warning=FALSE--------------------------------------------- # n <- 250 # alpha <- 0.01 # set.seed(1) # # # LR_cc # p.chi.cc <- replicate( # 1000, # ExactVaRTest::lr_cc_stat(rbinom(n, 1, alpha), alpha) > qchisq(.95, 2) # ) # p.exact.cc <- replicate( # 1000, # { # x <- rbinom(n, 1, alpha) # ExactVaRTest::backtest_lr(x, alpha, "cc")$pval < 0.05 # } # ) # # # LR_ind # p.chi.ind <- replicate( # 1000, # ExactVaRTest::lr_ind_stat(rbinom(n, 1, alpha), alpha) > qchisq(.95, 1) # ) # p.exact.ind <- replicate( # 1000, # { # x <- rbinom(n, 1, alpha) # ExactVaRTest::backtest_lr(x, alpha, "ind")$pval < 0.05 # } # ) # # data.frame( # Test = c("LR_cc Chi^2", "LR_cc Exact", "LR_ind Chi^2", "LR_ind Exact"), # Size = c(mean(p.chi.cc), mean(p.exact.cc), mean(p.chi.ind), mean(p.exact.ind)) # ) ## ----message=FALSE, warning=FALSE--------------------------------------------- # library(ggplot2) # library(dplyr) # library(tidyr) # # alpha <- 0.01 # win <- 250 # # local_file <- "inst/extdata/GSPC_close.rds" # file_path <- if (file.exists(local_file)) local_file else # system.file("extdata", "GSPC_close.rds", package = "ExactVaRTest") # # px <- readRDS(file_path) # px <- px[index(px) >= "2012-01-01"] # # ret <- diff(log(px)) # # VaR <- rollapply( # ret, win, # function(r) quantile(r, alpha, na.rm = TRUE), # by.column = FALSE, align = "right" # ) # # keep <- complete.cases(ret, VaR) # r <- coredata(ret[keep]) # v <- coredata(VaR[keep]) # exc <- ifelse(r < v, 1L, 0L) # # n_seg <- length(exc) - win + 1 # ind_exact <- ind_chi <- cc_exact <- cc_chi <- numeric(n_seg) # # for (i in seq_len(n_seg)) { # seg <- exc[i:(i + win - 1)] # ind_stat <- ExactVaRTest::lr_ind_stat(seg, alpha) # cc_stat <- ExactVaRTest::lr_cc_stat(seg, alpha) # ind_exact[i] <- ExactVaRTest::pval_lr_ind(ind_stat, win, alpha) # cc_exact[i] <- ExactVaRTest::pval_lr_cc(cc_stat, win, alpha) # ind_chi[i] <- pchisq(ind_stat, df = 1, lower.tail = FALSE) # cc_chi[i] <- pchisq(cc_stat, df = 2, lower.tail = FALSE) # } # # df <- tibble( # idx = seq_len(n_seg), # ind_exact, ind_chi, cc_exact, cc_chi # ) %>% # pivot_longer(cols = -idx, # names_to = c("test", "method"), # names_pattern = "(ind|cc)_(.*)", # values_to = "pval") %>% # mutate(test = ifelse(test == "ind", "LRind", "LRcc"), # method = ifelse(method == "exact", "exact", "chi-sq")) # # ggplot(df, aes(idx, pval, colour = method)) + # geom_step() + # geom_hline(yintercept = 0.05, linetype = 2, colour = "red") + # facet_wrap(~ test, ncol = 1, scales = "free_x") + # scale_colour_manual(values = c("chi-sq" = "red", "exact" = "cyan4")) + # labs(x = "window start index", y = "p-value", # title = "Rolling 250-day p-values (α = 1%)") + # theme_bw() ## ----example------------------------------------------------------------------ # library(ExactVaRTest) # # n_set <- c(250, 500, 750, 1000) # alpha_set <- c(0.005, 0.01, 0.025, 0.05) # gamma_set <- c(0.90, 0.95, 0.99) # # q_lr <- function(d, g) d$LR[which(cumsum(d$prob) >= g)[1L]] # # tbl <- expand.grid(n = n_set, # alpha = alpha_set, # gamma = gamma_set, # KEEP.OUT.ATTRS = FALSE, # stringsAsFactors = FALSE) # # tbl$crit_ind <- mapply(function(n, a, g) # q_lr(lr_ind_dist(n, a, prune_threshold = 1e-15), g), # tbl$n, tbl$alpha, tbl$gamma, SIMPLIFY = TRUE) # # tbl$crit_cc <- mapply(function(n, a, g) { # d <- lr_cc_dist(n, a, prune_threshold = 1e-15) # d$LR <- d$LR_cc # d$prob <- d$prob_cc # q_lr(d, g) # }, tbl$n, tbl$alpha, tbl$gamma, SIMPLIFY = TRUE) # # print(tbl, digits = 6) ## ----message=FALSE, warning=FALSE--------------------------------------------- # library(ExactVaRTest) # # set.seed(123) # # P <- 250 # alpha <- 0.05 # alpha_prime <- 0.10 # # inst_flag <- rbinom(P, 1, alpha_prime) # sys_flag <- if (sum(inst_flag)) rbinom(sum(inst_flag), 1, alpha) else integer(0) # # lr_uc <- lr_uc_stat(sys_flag, alpha) # lr_ind <- lr_ind_stat(sys_flag, alpha) # # mix_tail <- function(lr_obs, P, alpha, alpha_prime, # type = c("uc", "ind"), prune = 1e-15) { # type <- match.arg(type) # w <- dbinom(0:P, P, alpha_prime) # # tail_prob <- function(k) { # if (type == "uc") { # if (!k) return(as.numeric(lr_obs <= 0)) # d <- lr_uc_dist(k, alpha) # sum(d$prob[d$LR >= lr_obs]) # } else { # if (k < 2) return(as.numeric(lr_obs <= 0)) # d <- lr_ind_dist(k, alpha, prune) # sum(d$prob[d$LR >= lr_obs]) # } # } # # sum(vapply(0:P, tail_prob, numeric(1)) * w) # } # # p_uc <- mix_tail(lr_uc, P, alpha, alpha_prime, "uc") # p_ind <- mix_tail(lr_ind, P, alpha, alpha_prime, "ind") # # data.frame( # test = c("UC", "IND"), # stat = c(lr_uc, lr_ind), # p = c(p_uc, p_ind) # ) ## ----session-info, echo=FALSE------------------------------------------------- # sessionInfo()