---
title: "Monte Carlo simulation study"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Monte Carlo simulation study}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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")
}
```

This study checks that the rewrite (2.0.0) estimates the ARMA-GARCH/APARCH
models with stable and GEV innovations **at least as accurately as** the 2012
version (do Rego Sousa, master thesis), and adds the consistency and
standard-error calibration evidence the thesis did not report. The reproducible
scripts live in `dev/benchmarks/arma-garch-mc/`; the figures below are built
from precomputed results (the study takes hours and cannot run during a package
build).

## Part A — non-degradation against the 2012 thesis

Every maximum-likelihood simulation table of the thesis appendix (Tables 9.2,
9.7, 9.11–9.17) is reproduced with 2.0.0: **100 series of length n = 2500**,
the same parameter sets. For each parameter we compare the new RMSE to the 2012
RMSE. The key picture: plot the two against each other. A point **on or below
the diagonal** is as accurate as, or better than, 2012.

```{r 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)
```

Most points sit on or below the diagonal. Summarised by how the ratio
`rmse_new / rmse_2012` falls (ties are within Monte Carlo error at R = 100):

```{r 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)
```

The stable family (the package's headline) has **no regressions**; the GEV
family is median-unbiased, with a few higher-variance parameters in the
ARMA(2,2) cells — which the next section shows are a small-sample effect.

<details>
<summary>Full per-parameter table (173 rows)</summary>

```{r 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)
```

</details>

## Consistency

For the GEV ARMA(2,2) cells (where Part A showed the only higher variance), RMSE
is tracked across sample sizes n = 2500, 5000, 10000 (R = 100). Every parameter
(one line) declines with n — the empirical face of consistency.

```{r 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)")
```

## Part B — standard-error calibration (Wald coverage)

On a representative subset (two AR(1)-GARCH per family and one APARCH each),
fits at n = 500, 1000, 2500, 5000 (R = 300) record the 95% Wald coverage. The
shaded band is a reasonable [0.90, 0.98] window around the 0.95 target.
Mean-equation parameters (`mu`, `ar1`, `ar2`) are highlighted.

```{r 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)
```

Coverage is near nominal for the well-identified parameters. Two limitations
stand out, both expected:

* **Stable mean/AR parameters undercover with Wald** — in `9.7:phi1` the
  coverage of `mu` and `ar1` is ~0.35–0.44 at every n, because under
  infinite-variance stable innovations these estimators have a non-Gaussian
  limit, so symmetric Wald intervals are invalid. The parametric bootstrap of
  `gs_bootstrap()` restores calibration (next section).
* **Near-non-stationary GEV cells** (`9.2:theta2`, `theta3`) have mostly
  unavailable Hessian standard errors: their persistence exceeds one under the
  scale-1 GEV standardization (infinite unconditional variance), for which
  Hessian SEs are meaningless. Point estimates stay consistent. (Flagged in
  `TODO.md` for investigation.)

## Bootstrap vs Wald for stable mean/AR parameters

For `9.7:phi1` at n = 2500, the 95% parametric-bootstrap coverage next to the
Wald coverage (R = 80 replications, B = 99 resamples). The bootstrap lifts the
`mu`/`ar1` coverage back toward 0.95.

```{r 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)
```

```{r 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.)*")
```

## Reference

do Rego Sousa, T. (2012). *Modelos combinados AR-GARCH*. Master thesis,
Universidade de Brasília. Appendix, Chapter 9 (MLE Tables 9.2, 9.7, 9.11–9.17).