minutos_ano <- 252 * 390 # NYSE: ~390 min/dia x 252 dias uteis
calc_RV <- function(r) sum(r^2, na.rm = TRUE)
calc_BV <- function(r) {
if (length(r) < 2) return(NA_real_)
(pi / 2)^(-1) * sum(abs(r[-1]) * abs(r[-length(r)]), na.rm = TRUE)
}
annualize_vol_from_rv <- function(rv, n_obs) {
if (n_obs <= 0) return(NA_real_)
sqrt(rv / n_obs) * sqrt(minutos_ano)
}
bootstrap_RV <- function(r, R = 1500) {
r <- na.omit(as.numeric(r))
if (length(r) < 5) return(NULL)
b <- boot::boot(data = r, statistic = function(d, i) sum(d[i]^2), R = R)
ci <- boot::boot.ci(b, type = c("perc", "bca"))
list(boot = b, ci = ci)
}
perm_test_RV_diff <- function(r1, r2, R = 2000) {
r1 <- na.omit(as.numeric(r1)); r2 <- na.omit(as.numeric(r2))
obs_diff <- calc_RV(r2) - calc_RV(r1)
pooled <- c(r1, r2); n1 <- length(r1); n2 <- length(r2)
set.seed(123)
perm_diffs <- replicate(R, {
p <- sample(pooled)
calc_RV(p[(n1+1):(n1+n2)]) - calc_RV(p[1:n1])
})
list(obs_diff = obs_diff, p_value = mean(abs(perm_diffs) >= abs(obs_diff)))
}
calc_for_segment <- function(name, idx, ret_vec) {
r <- as.numeric(na.omit(ret_vec[idx]))
rv <- calc_RV(r); bv <- calc_BV(r)
jmp <- max(rv - bv, 0)
list(name = name, n = length(r),
RV = rv, BV = bv, Jump = jmp,
JumpShare = ifelse(rv > 0, jmp / rv, NA_real_),
VolAnn = annualize_vol_from_rv(rv, length(r)),
mean = mean(r), sd = sd(r),
skewness = moments::skewness(r), kurtosis = moments::kurtosis(r),
boot = bootstrap_RV(r, R = 1500))
}
all_rv_results <- vector("list", 3)
for (k in 1:3) {
ev <- events[[k]]
janelas <- list(pre = ev$pre, during = ev$during, post = ev$post)
results <- lapply(names(janelas), function(nm) {
if (is.null(janelas[[nm]])) return(NULL)
calc_for_segment(nm, janelas[[nm]], rets_full)
})
names(results) <- names(janelas)
results <- Filter(Negate(is.null), results)
all_rv_results[[k]] <- results
cat("\n\n===", ev$event_label, "| Crash:", format(ev$crash_date, "%Y-%m-%d"), "===\n")
summary_tbl <- do.call(rbind, lapply(results, function(x) {
data.frame(janela = x$name, n = x$n, RV = x$RV, BV = x$BV,
Jump = x$Jump, JumpShare = x$JumpShare, VolAnn = x$VolAnn,
mean = x$mean, sd = x$sd,
skewness = x$skewness, kurtosis = x$kurtosis)
}))
print(knitr::kable(summary_tbl, digits = 6,
caption = paste0(ev$event_label, " - RV, Bipower, Jump, Vol Anualizada")))
# Diferencas de RV
pairs <- list(c("pre","during"), c("pre","post"))
for (pr in pairs) {
a <- pr[1]; b <- pr[2]
if (!is.null(results[[a]]) && !is.null(results[[b]])) {
diff_rv <- results[[b]]$RV - results[[a]]$RV
pct_rv <- diff_rv / abs(results[[a]]$RV) * 100
cat(sprintf(" RV (%s - %s): %+.6f (%+.2f%%)\n", b, a, diff_rv, pct_rv))
}
}
# Testes permutacionais
for (pr in pairs) {
a <- pr[1]; b <- pr[2]
if (!is.null(janelas[[a]]) && !is.null(janelas[[b]])) {
perm_res <- perm_test_RV_diff(rets_full[janelas[[a]]], rets_full[janelas[[b]]], R = 2000)
cat(sprintf(" Perm test (%s vs %s): p = %.4f\n", b, a, perm_res$p_value))
}
}
# Levene
seg_rets <- unlist(lapply(names(janelas), function(nm) {
if (!is.null(janelas[[nm]])) as.numeric(rets_full[janelas[[nm]]]) else NULL
}))
seg_labs <- unlist(lapply(names(janelas), function(nm) {
if (!is.null(janelas[[nm]])) rep(nm, length(janelas[[nm]])) else NULL
}))
combined <- data.frame(ret = seg_rets,
seg = factor(seg_labs, levels = c("pre","during","post")))
lev <- car::leveneTest(ret ~ seg, data = combined)
cat(sprintf(" Levene: F = %.4f | p = %.6f\n", lev[1,"F value"], lev[1,"Pr(>F)"]))
# KS
for (pr in pairs) {
a <- pr[1]; b <- pr[2]
if (!is.null(janelas[[a]]) && !is.null(janelas[[b]])) {
ks <- ks.test(as.numeric(rets_full[janelas[[a]]]),
as.numeric(rets_full[janelas[[b]]]))
cat(sprintf(" KS (%s vs %s): D = %.4f | p = %.6f\n", b, a, ks$statistic, ks$p.value))
}
}
}
