## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE,
                      fig.width = 7, fig.height = 4.2, dpi = 120,
                      fig.align = "center", collapse = TRUE, comment = "#>")
has_gg <- requireNamespace("ggplot2", quietly = TRUE)
read_mc <- function(f) read.csv(file.path("mc-tables", f), check.names = FALSE)
stable_tables <- c("9.7", "9.11", "9.12", "9.13", "9.14")
if (has_gg) {
  library(ggplot2)
  theme_set(theme_minimal(base_size = 11) +
              theme(panel.grid.minor = element_blank(),
                    legend.position = "top", plot.title = element_text(face = "bold")))
  pal <- c(stable = "#2b6cb0", GEV = "#c05621", `mean/AR` = "#c53030", other = "#718096")
}

## ----partA-scatter, eval = has_gg---------------------------------------------
a <- read_mc("partA_R100_n2500_comparison.csv")
a$family <- ifelse(a$table %in% stable_tables, "stable", "GEV")
cmp <- a[a$comparable & is.finite(a$ratio), ]

ggplot(cmp, aes(rmse_2012, rmse, colour = family)) +
  geom_abline(slope = 1, linetype = 2, colour = "grey40") +
  geom_abline(slope = 1.1, linetype = 3, colour = "grey60") +
  geom_point(alpha = 0.8, size = 2) +
  scale_x_log10() + scale_y_log10() +
  scale_colour_manual(values = pal) +
  annotate("text", x = min(cmp$rmse_2012) * 1.1, y = max(cmp$rmse) * 0.9,
           label = "above line = worse than 2012", hjust = 0, size = 3, colour = "grey40") +
  labs(title = "Per-parameter RMSE: 2.0.0 vs 2012",
       subtitle = "Dashed = equal; dotted = the 1.10 non-degradation gate",
       x = "RMSE 2012 (log scale)", y = "RMSE 2.0.0 (log scale)", colour = NULL)

## ----partA-bars, eval = has_gg, fig.height = 2.6------------------------------
cmp$band <- cut(cmp$ratio, c(0, 0.9, 1.15, 1.3, Inf),
                labels = c("better\n(<0.9)", "tie\n(0.9–1.15)",
                           "mild\n(1.15–1.3)", "worse\n(>1.3)"))
tab <- as.data.frame(table(family = cmp$family, band = cmp$band))
ggplot(tab, aes(band, Freq, fill = family)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(0.9), vjust = -0.3, size = 3) +
  scale_fill_manual(values = pal) +
  labs(title = "How the RMSE ratio falls, by family",
       x = NULL, y = "number of parameters", fill = NULL)

## ----partA-table--------------------------------------------------------------
ft <- a[order(a$table, a$param_set, a$param),
        c("cell", "param", "true", "rmse", "rmse_2012", "ratio", "comparable", "pass")]
names(ft) <- c("Cell", "Param", "True", "RMSE 2.0.0", "RMSE 2012", "Ratio", "Comparable", "Pass")
knitr::kable(ft, digits = c(0, 0, 3, 4, 4, 2, 0, 0), row.names = FALSE)

## ----consistency, eval = has_gg, fig.height = 4.6-----------------------------
ct <- read.csv(file.path("mc-tables", "consistency_trend.csv"), check.names = FALSE)
long <- reshape(ct[c("cell", "param", "rmse_n2500", "rmse_n5000", "rmse_n10000")],
                varying = list(3:5), v.names = "rmse",
                times = c(2500, 5000, 10000), timevar = "n", direction = "long")
ggplot(long, aes(n, rmse, group = interaction(cell, param))) +
  geom_line(alpha = 0.35, colour = "#2b6cb0") +
  geom_point(alpha = 0.5, size = 1, colour = "#2b6cb0") +
  facet_wrap(~cell) +
  scale_x_log10(breaks = c(2500, 5000, 10000)) + scale_y_log10() +
  labs(title = "RMSE falls with sample size (every parameter)",
       x = "n (log scale)", y = "RMSE (log scale)")

## ----coverage, eval = has_gg, fig.height = 4.8--------------------------------
b <- read_mc("partB_coverage.csv")
b$kind <- ifelse(b$param %in% c("mu", "ar1", "ar2"), "mean/AR", "other")
b <- b[is.finite(b$coverage95), ]
ggplot(b, aes(n, coverage95, group = interaction(cell, param), colour = kind)) +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.90, ymax = 0.98,
           fill = "grey85", alpha = 0.5) +
  geom_hline(yintercept = 0.95, linetype = 2, colour = "grey40") +
  geom_line(alpha = 0.6) + geom_point(size = 1, alpha = 0.8) +
  facet_wrap(~cell) +
  scale_x_log10(breaks = c(500, 1000, 2500, 5000)) +
  scale_colour_manual(values = pal) +
  labs(title = "95% Wald coverage by sample size",
       x = "n (log scale)", y = "coverage", colour = NULL)

## ----bootstrap, eval = has_gg && file.exists(file.path("mc-tables", "bootstrap_spotcheck_phi1_n2500.csv")), fig.height = 3.4----
bs <- read.csv(file.path("mc-tables", "bootstrap_spotcheck_phi1_n2500.csv"), check.names = FALSE)
names(bs)[1:3] <- c("param", "Wald", "Bootstrap")
bl <- reshape(bs[c("param", "Wald", "Bootstrap")], varying = list(2:3),
              v.names = "coverage", times = c("Wald", "Bootstrap"),
              timevar = "method", direction = "long")
bl$method <- factor(bl$method, levels = c("Wald", "Bootstrap"))
ggplot(bl, aes(param, coverage, fill = method)) +
  geom_col(position = "dodge") +
  geom_hline(yintercept = 0.95, linetype = 2, colour = "grey40") +
  scale_fill_manual(values = c(Wald = "#c53030", Bootstrap = "#2b6cb0")) +
  labs(title = "9.7:phi1 (n = 2500): Wald vs bootstrap 95% coverage",
       x = NULL, y = "coverage", fill = NULL)

## ----bootstrap-missing, eval = !file.exists(file.path("mc-tables", "bootstrap_spotcheck_phi1_n2500.csv")), results = "asis"----
# cat("*(Bootstrap spot-check figure pending — will appear once the run completes.)*")