---
title: "Testing the last observation for instability"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Testing the last observation for instability}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- Render time: ~12 s under rmarkdown::render() with ggplot2 installed,
     most of it in the size study; macOS arm64 (Apple silicon), R 4.6.1,
     one core. The comparison with strucchange, changepoint, KFAS and dlm
     is read from results/end_of_sample.rds, built by
     data-raw/vignette_results/end_of_sample.R. -->

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  dpi = 150,
  out.width = "100%"
)
```

```{r library}
library(proxymix)
```

```{r engines}
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
```

```{r stored-results, include = FALSE}
## The comparison table reads stored simulation results. They must come
## from the same major.minor version of proxymix as this build.
res <- readRDS("results/end_of_sample.rds")
major_minor <- function(v) paste(unlist(package_version(v))[1:2],
                                 collapse = ".")
if (major_minor(res$proxymix_version) !=
    major_minor(as.character(packageVersion("proxymix")))) {
  stop("results/end_of_sample.rds was built under proxymix ",
       res$proxymix_version, ", but this is proxymix ",
       packageVersion("proxymix"), ". Rerun the simulation and ",
       "data-raw/vignette_results/end_of_sample.R.", call. = FALSE)
}

## Small numbers are written as plain decimals rather than in the
## scientific notation that knitr's inline hook would otherwise use.
fixed <- function(v, digits) {
  format(round(v, digits), nsmall = digits, scientific = FALSE)
}
```

## The problem

Suppose you follow a series over time, such as monthly sales, and you have
a model fitted to every value so far. One new value has just arrived. You
want to know whether it is consistent with the model, or whether something
has changed.

Most tests for a change in a series need data on both sides of the change.
The breakpoint test of Chow (1960) and the sup-Wald test of Andrews (1993),
which tries every possible change point, estimate the model before and
after a change and compare the two.
With one new value, or a handful, the model after the change cannot be
estimated. Chow (1960) also proposed a predictive test that works in this
case, but it assumes that the errors follow a normal distribution
(Andrews, 2003).

Andrews (2003) proposed end-of-sample tests for exactly this case.
proxymix implements a version of them for state-space models. These models
describe a series through an unobserved state, such as a level, that
changes over time.

## Package capabilities

- `gmm_filter()` runs the Kalman filter, which tracks the unobserved state
  of a series one time step at a time. Before each new value arrives, the
  filter forecasts it. The difference between the value and its forecast is
  called the innovation.
- `gmm_eos_test()` tests whether the last `m` values of a series are
  consistent with the model. It takes the model in the same form as
  `gmm_filter()`. It returns a p-value, a decision at the 5% level, and the
  numbers used to compute them.

The test statistic is built from the innovations. Each innovation is
divided by its standard deviation, which the filter also reports. If the
model is correct and nothing has changed, the result, $z_t$, follows a
standard normal distribution. The statistic adds up the squares of the last
`m` of them:
$$ \mathrm{EoS}_m = \sum_{t = n-m+1}^{n} z_t^2. $$
A change in the last `m` values makes their forecasts worse and the
statistic larger.

`gmm_eos_test()` offers two ways, or calibrations, of turning the statistic
into a p-value.

- `method = "chisq"`, the default, compares the statistic with a chi-square
  distribution, the distribution of a sum of squared standard normal values.
  Its degrees of freedom, the number of squared values in the sum, are `m`
  times the number of variables observed at each time step. The p-value is exact when the innovations are normal and
  the model's parameters are known.
- `method = "andrews"` compares the statistic with the same statistic
  computed on every earlier block of `m` consecutive values of the series.
  The p-value is the share of blocks at least as large as the final block,
  counting the final block itself. Taking p-values from the data in this
  way is called subsampling. It follows the P-test of Andrews (2003), with
  two differences: the p-value counts the tested block, and the earlier
  blocks are computed from the supplied model rather than from a model
  re-estimated without each block.

The subsampling calibration does not assume normal innovations. It does
make other assumptions. Its p-value is valid only as the series grows
long with `m` fixed, and only if the innovations before the tested block
are stationary, meaning that their distribution does not change over time.
Its size in a short series can differ from nominal in either direction:
the discreteness of the p-value makes it conservative at `m = 1` with the
model's parameters known, but it can reject too often when the blocks are
no longer exchangeable, such as with estimated parameters or the
overlapping blocks of `m > 1`.

## Addressing the problem

```{r seed}
set.seed(20260621)
```

### A model and a stable series

The model is a local-level model. An unobserved level moves as a random
walk, taking a small normal step at each time point. Each observation is
the level plus normal noise. In `gmm_filter()` form, the prior gives the
starting level, `dynamics` gives the steps of the level, and `measurement`
gives the noise of the observations. `gmm()` builds a mixture of normal
distributions, and with `weights = 1` it is a single normal distribution.
A prior variance of 10 means the starting level is known only roughly.

```{r model}
prior <- gmm(weights = 1, means = list(0), covariances = list(matrix(10)))
dynamics <- list(A = matrix(1), Q = matrix(0.04))  # the level's random walk
measurement <- list(C = matrix(1), R = matrix(1))  # the observation noise

n <- 120L
level <- cumsum(c(0, rnorm(n - 1L, 0, sqrt(0.04))))
y_stable <- level + rnorm(n, 0, 1)
```

### The same series with a broken final value

A copy of the series has its final value moved up by five standard
deviations of the observation noise. Both series are tested at `m = 1`
with the subsampling calibration.

```{r break}
y_break <- y_stable
y_break[n] <- y_break[n] + 5

test_stable <- gmm_eos_test(
  prior, dynamics, measurement, y_stable, m = 1L, method = "andrews"
)
test_break <- gmm_eos_test(
  prior, dynamics, measurement, y_break, m = 1L, method = "andrews"
)
```

```{r tests-kable, echo = FALSE}
eos_row <- function(x) {
  c(
    round(x$statistic, 3L), round(x$p_value, 4L), x$reject,
    x$method, x$m
  )
}
knitr::kable(
  data.frame(
    field = c("statistic", "p-value", "reject at 0.05", "calibration",
              "window m"),
    stable = eos_row(test_stable),
    broken = eos_row(test_break)
  ),
  col.names = c("Result", "Stable series", "Broken final value"),
  caption = paste(
    "The end-of-sample test on the same series before and after its",
    "final value is moved up by five standard deviations of the",
    "observation noise."
  )
)
```

### What the filter expected

The figure shows the forecasts behind the statistic. Before each value
arrives, the filter forecasts it from the values before it. The shaded
band spans two forecast standard deviations either side of each forecast.
The statistic divides each innovation by this forecast standard deviation.

```{r filter-path}
filtered <- gmm_filter(
  prior, dynamics, measurement, y_stable, ridge_eps = 0
)
# ridge_eps adds a tiny amount to each covariance for numerical stability;
# 0 keeps the filter's recursion exact
path <- filtered$summary
# forecast of y_t and its standard deviation, from the filter at t - 1
pred_mean <- c(prior@means[[1L]], path$mean_1[-n])
pred_sd <- sqrt(c(prior@covariances[[1L]][1L, 1L], path$sd_1[-n]^2) +
                  dynamics$Q[1L, 1L] + measurement$R[1L, 1L])
outside <- sum(abs(y_stable - pred_mean)[-1L] > 2 * pred_sd[-1L])
```

```{r fig-series, eval = has_ggplot2, echo = has_ggplot2, fig.height = 4, fig.cap = "The simulated series, the unobserved level that generated it, the level estimated by the filter, and a band of two forecast standard deviations either side of each forecast. The band starts at the second time step. The broken final value lies far outside the band.", fig.alt = "A time series of noisy observations in grey with the true and estimated level overlaid, a shaded forecast band around them, and a single isolated point at the right-hand end lying well above the band."}
series_df <- data.frame(
  t = seq_len(n),
  observed = y_stable,
  truth = level,
  filtered = path$mean_1,
  lo = pred_mean - 2 * pred_sd,
  hi = pred_mean + 2 * pred_sd
)
break_df <- data.frame(t = n, y = y_break[n])
ggplot2::ggplot(series_df, ggplot2::aes(t)) +
  ggplot2::geom_ribbon(
    data = series_df[-1L, ], ggplot2::aes(ymin = lo, ymax = hi),
    fill = "#0072B2", alpha = 0.18
  ) +
  ggplot2::geom_point(
    ggplot2::aes(y = observed, colour = "observed (stable)"),
    size = 1.1, alpha = 0.7
  ) +
  ggplot2::geom_line(
    ggplot2::aes(y = truth, colour = "true level"), linewidth = 0.8
  ) +
  ggplot2::geom_line(
    ggplot2::aes(y = filtered, colour = "estimated level"), linewidth = 0.8
  ) +
  ggplot2::geom_point(
    data = break_df, ggplot2::aes(t, y, colour = "broken final value"),
    size = 2.6
  ) +
  ggplot2::scale_colour_manual(
    name = NULL,
    values = c(
      "observed (stable)" = "grey60",
      "true level" = "#009E73",
      "estimated level" = "#0072B2",
      "broken final value" = "#D55E00"
    )
  ) +
  ggplot2::labs(
    x = "time step", y = "observation",
    title = "A break in the last observation"
  ) +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::theme(legend.position = "top")
```

```{r fig-series-skip, eval = !has_ggplot2, echo = FALSE, results = "asis"}
cat("ggplot2 is not installed on this build, so the series figure is",
    "skipped.\n")
```

### How the subsampling p-value is reached

The subsampling calibration ranks the final block among the earlier blocks
of the same series. The next figure shows that ranking for both series.

```{r fig-blocks, eval = has_ggplot2, echo = has_ggplot2, fig.height = 4, fig.cap = "The cumulative distribution of the earlier block statistics in each series, with the final block's statistic as a dashed line. The p-value is $(1 + k) / (n - 2m + 2)$, where $k$ is the number of earlier blocks at or beyond the dashed line. The stable series' final block lies inside the range of earlier blocks. The broken series' final block lies beyond every earlier block.", fig.alt = "Two panels, each showing a rising step function of the cumulative share of earlier block statistics against a logarithmic axis, with a dashed vertical line. In the left panel the line falls inside the curve. In the right panel it lies far beyond the curve's right-hand end."}
blocks_df <- rbind(
  data.frame(
    block = test_stable$in_sample_blocks, series = "stable series"
  ),
  data.frame(
    block = test_break$in_sample_blocks, series = "broken final value"
  )
)
blocks_df$series <- factor(
  blocks_df$series, levels = c("stable series", "broken final value")
)
rule_df <- data.frame(
  series = factor(
    c("stable series", "broken final value"),
    levels = c("stable series", "broken final value")
  ),
  statistic = c(test_stable$statistic, test_break$statistic)
)
ggplot2::ggplot(blocks_df, ggplot2::aes(block)) +
  ggplot2::stat_ecdf(geom = "step", linewidth = 0.8, colour = "#0072B2") +
  ggplot2::geom_vline(
    data = rule_df, ggplot2::aes(xintercept = statistic),
    colour = "#D55E00", linetype = "dashed", linewidth = 0.8
  ) +
  ggplot2::facet_wrap(~ series) +
  ggplot2::scale_x_log10(labels = function(v) {
    format(v, scientific = FALSE, drop0trailing = TRUE, trim = TRUE)
  }) +
  ggplot2::labs(
    x = "block statistic (log scale)",
    y = "cumulative share of earlier blocks",
    title = "The final block against the earlier blocks of the same series"
  ) +
  ggplot2::theme_minimal(base_size = 11)
```

```{r fig-blocks-skip, eval = !has_ggplot2, echo = FALSE, results = "asis"}
cat("ggplot2 is not installed on this build, so the block figure is",
    "skipped.\n")
```

### How often each calibration rejects when the model is known

The size of a test is how often it rejects when nothing has changed. A test
at the 5% level should have a size of 0.05. The chunk below simulates
stable series of 30 values from the same model, tests each one, and records
the share of p-values below 0.05. The model's parameters are known here,
as in the example above.

```{r size-study}
simulate_null <- function(n_obs) {
  lv <- cumsum(c(0, rnorm(n_obs - 1L, 0, sqrt(0.04))))
  lv + rnorm(n_obs, 0, 1)
}
empirical_size <- function(n_obs, method, n_rep) {
  p <- vapply(seq_len(n_rep), function(i1) {
    gmm_eos_test(
      prior, dynamics, measurement, simulate_null(n_obs),
      m = 1L, method = method
    )$p_value
  }, numeric(1L))
  size <- mean(p < 0.05)
  c(size = size, se = sqrt(size * (1 - size) / n_rep))
}

n_rep <- 250L
size_chisq_30 <- empirical_size(30L, "chisq", n_rep)
size_andrews_30 <- empirical_size(30L, "andrews", n_rep)
```

The subsampling p-value can take only the values $(1 + k) / (n - 2m + 2)$
for $k = 0, 1, 2, \ldots$. Its smallest value is therefore
$1 / (n - 2m + 2)$. If the final block is equally likely to take any rank
among the blocks, its size is the share of those values below 0.05. This
holds at $m = 1$ with known parameters, apart from the first few
innovations, which the rough starting level affects. It does not hold when
the parameters are estimated or when $m > 1$.

```{r grid-size}
grid_size <- function(n_obs, m = 1L, alpha = 0.05) {
  n_block <- n_obs - 2L * m + 1L
  p_grid <- (1L + seq.int(0L, n_block)) / (1L + n_block)
  mean(p_grid < alpha)
}
grid_30 <- grid_size(30L)
grid_120 <- grid_size(120L)
```

```{r size-checks, include = FALSE}
## the prose below says the chi-square size does not differ detectably
## from 0.05 and the subsampling size agrees with its grid value
stopifnot(abs(size_chisq_30[["size"]] - 0.05) < 2 * size_chisq_30[["se"]],
          abs(size_andrews_30[["size"]] - grid_30) <
            2 * size_andrews_30[["se"]])
```

```{r size-kable, echo = FALSE}
old_opts <- options(knitr.kable.NA = "--")
knitr::kable(
  data.frame(
    calibration = c("chi-square", "subsampling", "subsampling",
                    "subsampling"),
    n_obs = c(30L, 30L, 30L, 120L),
    how = c("simulated", "simulated", "from the p-value grid",
            "from the p-value grid"),
    size = round(c(size_chisq_30[["size"]], size_andrews_30[["size"]],
                   grid_30, grid_120), 4L),
    se = round(c(size_chisq_30[["se"]], size_andrews_30[["se"]],
                 NA_real_, NA_real_), 4L)
  ),
  col.names = c("Calibration", "Series length", "Obtained", "Size",
                "Standard error"),
  caption = paste0(
    "Size at a nominal 0.05 with the model's parameters known, for m = 1. ",
    "Simulated sizes use ", n_rep, " stable series each. Sizes from the ",
    "p-value grid are exact when the final block's rank is equally likely ",
    "to be any rank, so they have no standard error."
  )
)
options(old_opts)
```

### Comparison with strucchange, changepoint, KFAS and dlm

```{r compare-facts, include = FALSE}
sim_value <- function(m, scenario, method) {
  s1 <- res$sim_tab$m == m & res$sim_tab$scenario == scenario &
    res$sim_tab$method == method
  res$sim_tab$rate[s1]
}
size_at <- function(method, m) fixed(sim_value(m, "no shift", method), 3)
power_at <- function(method, m) fixed(sim_value(m, "shift", method), 3)
mcse <- sqrt(0.05 * 0.95 / res$n_rep)
ms_of <- function(pkg) fixed(1000 * res$time_secs[[pkg]], 0)

## the prose below rests on these orderings in the stored results
size_1 <- vapply(unique(res$sim_tab$method), function(k) {
  sim_value(1L, "no shift", k)
}, numeric(1L))
stopifnot(identical(names(size_1)[!is.na(size_1) & abs(size_1 - 0.05) <= mcse],
                    "proxymix subsampling"))
stopifnot(abs(sim_value(5L, "no shift", "strucchange sup-F") - 0.05) <
            min(abs(sim_value(5L, "no shift", "proxymix chi-square") - 0.05),
                abs(sim_value(5L, "no shift", "proxymix subsampling") - 0.05)))
stopifnot(sim_value(1L, "shift", "proxymix subsampling") <
            min(sim_value(1L, "shift", "proxymix chi-square"),
                sim_value(1L, "shift", "KFAS forecast interval"),
                sim_value(1L, "shift", "dlm forecast interval")))
stopifnot(names(which.max(res$time_secs)) == "KFAS",
          res$time_secs[["strucchange"]] < 0.001,
          res$time_secs[["changepoint"]] < 0.001)
```

In practice the model's parameters are estimated from the same series that
is tested. A simulation compared the two calibrations with six other
checks in this setting. Each of `r res$n_rep` series had `r res$n` values
from an autoregressive model, in which each value is `r res$phi` times the
previous one plus standard normal noise. Every check was applied at
$m$ = 1, 3 and 5, where it is defined. proxymix's two calibrations, KFAS, dlm and the t-test
were fitted to all but the last $m$ values and then tested against those
values; the `strucchange` and `cpt.mean()` checks below instead ran on the
whole series. The power is how often a check rejected when the last $m$
values were raised by `r res$delta` noise standard deviations. The size is
how often it rejected when they were not.

The other checks were the sup-F test and the OLS-CUSUM test of
`strucchange` (Zeileis et al., 2002), `cpt.mean()` of `changepoint`
(Killick and Eckley, 2014), and a t-test on the forecast errors of the last
$m$ values. The sup-F test looks for a change near the end of the series
in the regression of each value on the previous one. The OLS-CUSUM test looks for drift, over the whole series, in the
cumulative sum of the errors of that regression. `cpt.mean()` searches for
changes in the mean, and it counts as rejecting when it places a change
among the last $m$ values. `KFAS` (Helske, 2017) and `dlm` (Petris, 2010)
each fitted their own autoregressive model and checked every tested value
against its 95% forecast interval. Each interval had a level of $1 - 0.05 / m$,
which splits the 5% error rate evenly over the $m$ values (a Bonferroni
correction). The check then has a size of at most 0.05 when the parameters
are known.

```{r compare-table, echo = FALSE}
method_order <- c("proxymix chi-square", "proxymix subsampling",
                  "strucchange sup-F", "strucchange OLS-CUSUM",
                  "changepoint cpt.mean", "KFAS forecast interval",
                  "dlm forecast interval", "t-test on forecast residuals")
method_label <- c("proxymix, chi-square (default)", "proxymix, subsampling",
                  "strucchange, sup-F", "strucchange, OLS-CUSUM",
                  "changepoint, cpt.mean()", "KFAS, forecast interval",
                  "dlm, forecast interval", "t-test on forecast errors")
rate_col <- function(scenario, m) {
  v <- vapply(method_order, function(k) sim_value(m, scenario, k),
              numeric(1L))
  ifelse(is.na(v), "--", format(round(v, 3L), nsmall = 3L))
}
m_shown <- c(1L, 5L)
cmp_tbl <- data.frame(method = method_label, stringsAsFactors = FALSE)
for (m1 in m_shown) {
  cmp_tbl[[paste0("size_", m1)]] <- rate_col("no shift", m1)
  cmp_tbl[[paste0("power_", m1)]] <- rate_col("shift", m1)
}
knitr::kable(
  cmp_tbl, row.names = FALSE,
  align = c("l", rep("r", 2L * length(m_shown))),
  col.names = c("Check", as.vector(rbind(paste0("Size, m&nbsp;=&nbsp;", m_shown),
                                         paste0("Power, m&nbsp;=&nbsp;", m_shown)))),
  caption = paste0(
    "Rejection rates at a nominal 0.05 over ", res$n_rep, " simulated ",
    "series of ", res$n, " values, with parameters estimated from each ",
    "series. Size is the rate without a change and should be near 0.05. ",
    "Power is the rate when the last m values were raised by ", res$delta,
    " noise standard deviations. A size near 0.05 has a simulation ",
    "standard error of about ", fixed(mcse, 3), ". A dash marks a check ",
    "that cannot be computed at that m. The extended version of this ",
    "article also reports m = 3."
  )
)
```

At a single final value, only the subsampling calibration had a size
close to 0.05, at `r size_at("proxymix subsampling", 1)`. The default
chi-square calibration had a size of `r size_at("proxymix chi-square", 1)`,
close to that of `dlm` (`r size_at("dlm forecast interval", 1)`), and
that of `KFAS` was `r size_at("KFAS forecast interval", 1)`. All three treat the
estimated parameters as if they were known. Of these four forecast checks,
the subsampling calibration had the lowest power,
`r power_at("proxymix subsampling", 1)` against
`r power_at("proxymix chi-square", 1)` for the chi-square calibration,
`r power_at("KFAS forecast interval", 1)` for `KFAS` and
`r power_at("dlm forecast interval", 1)` for `dlm`. At $m = 5$ neither
calibration held its size (`r size_at("proxymix chi-square", 5)` and
`r size_at("proxymix subsampling", 5)`), and the sup-F test came closer, at
`r size_at("strucchange sup-F", 5)`. Its power was
`r power_at("strucchange sup-F", 5)`, against
`r power_at("proxymix chi-square", 5)` for the chi-square calibration. At $m = 3$, a column the table
omits, the sup-F test had a size of `r size_at("strucchange sup-F", 3)`,
far above 0.05. The
OLS-CUSUM test and `cpt.mean()` rejected fewer than 5% of stable series.
The OLS-CUSUM test rarely detected a change confined to the last few
values, with a power of `r power_at("strucchange OLS-CUSUM", 1)` at $m = 1$
and `r power_at("strucchange OLS-CUSUM", 5)` at $m = 5$. `cpt.mean()`
detected it in `r power_at("changepoint cpt.mean", 1)` of series at
$m = 1$ and `r power_at("changepoint cpt.mean", 5)` at $m = 5$.

`KFAS` was the slowest of the five packages. For one series at $m$ =
`r res$time_m`, it took `r ms_of("KFAS")` ms and `dlm` took
`r ms_of("dlm")` ms, each including its own fit. The two proxymix
calibrations together took `r ms_of("proxymix")` ms, including the `arima()`
fit that supplies the parameters. `strucchange` and `changepoint` each took
less than 1 ms (median of five runs of `r res$time_n_call` calls each, on one
computer).

The code below applies every check to one simulated series with `m = 3`.
It is the simulation code for a single series. It needs
`strucchange`, `changepoint`, `KFAS` and `dlm`, all on CRAN, and it is not
run when this vignette is built.

```{r compare-code, eval = FALSE}
library(proxymix)
library(strucchange)
library(changepoint)
library(KFAS)
library(dlm)

# one series of 100 values whose last m values are raised by 3
set.seed(1L)
n <- 100L
m <- 3L
y <- as.numeric(arima.sim(list(ar = 0.6), n = n))
y[(n - m + 1L):n] <- y[(n - m + 1L):n] + 3
y_fit <- y[seq_len(n - m)]
last <- (n - m + 1L):n

# proxymix, with the autoregressive model fitted by arima()
ar <- arima(y_fit, order = c(1L, 0L, 0L))
a1 <- coef(ar)[["ar1"]]
mu <- coef(ar)[["intercept"]]
s2 <- ar$sigma2
prior <- gmm(weights = 1, means = list(mu),
             covariances = list(matrix(s2 / (1 - a1^2))))
dynamics <- list(A = matrix(a1), b = mu * (1 - a1), Q = matrix(s2))
measurement <- list(C = matrix(1), R = matrix(0))
p_chisq <- gmm_eos_test(prior, dynamics, measurement, y, m = m,
                        method = "chisq")$p_value
p_andrews <- gmm_eos_test(prior, dynamics, measurement, y, m = m,
                          method = "andrews")$p_value
resid <- y[last] - predict(ar, n.ahead = m)$pred
p_t <- if (m > 1L) t.test(resid)$p.value else NA_real_

# strucchange needs at least three observations after a break in an AR(1)
reg <- data.frame(y = y[-1L], y_lag = y[-n])
n_reg <- nrow(reg)
p_supf <- if (m >= 3L) {
  fs <- Fstats(y ~ y_lag, data = reg, from = n_reg - m, to = n_reg - 3L)
  unname(sctest(fs)$p.value)
} else NA_real_
p_cusum <- unname(sctest(efp(y ~ y_lag, data = reg,
                             type = "OLS-CUSUM"))$p.value)

cp <- cpts(cpt.mean(y / sd(y_fit)))
p_cpt <- if (any(cp >= n - m)) 0 else 1

y_bar <- mean(y_fit)
kfas_update <- function(pars, model) {
  part <- SSMarima(ar = 0.999 * tanh(pars[1L]), Q = exp(pars[2L]))
  model["T", "arima"] <- part$T
  model["R", "arima"] <- part$R
  model["Q", "arima"] <- part$Q
  model["P1", "arima"] <- part$P1
  model
}
kfas_fit <- fitSSM(
  SSModel(I(y_fit - y_bar) ~ -1 + SSMarima(ar = 0.5, Q = 1), H = 0),
  inits = c(0, 0), updatefn = kfas_update, method = "BFGS"
)
kfas_par <- kfas_fit$optim.out$par
kfas_model <- SSModel(
  I(y - y_bar) ~ -1 + SSMarima(ar = 0.999 * tanh(kfas_par[1L]),
                                Q = exp(kfas_par[2L])),
  H = 0
)
kfas_pred <- predict(kfas_model, interval = "prediction", level = 0.95,
                     filtered = TRUE)
kfas_z <- (y[last] - y_bar - kfas_pred[last, "fit"]) /
  ((kfas_pred[last, "upr"] - kfas_pred[last, "lwr"]) / (2 * qnorm(0.975)))
p_kfas <- min(1, m * 2 * pnorm(-max(abs(kfas_z))))

dlm_build <- function(pars) {
  dlmModARMA(ar = 0.999 * tanh(pars[1L]), sigma2 = exp(pars[2L]), dV = 0)
}
dlm_fit <- dlmMLE(y_fit - y_bar, parm = c(0, 0), build = dlm_build)
dlm_filt <- dlmFilter(y - y_bar, dlm_build(dlm_fit$par))
dlm_var <- unlist(dlmSvd2var(dlm_filt$U.R, dlm_filt$D.R))
dlm_z <- (y[last] - y_bar - dlm_filt$f[last]) / sqrt(dlm_var[last])
p_dlm <- min(1, m * 2 * pnorm(-max(abs(dlm_z))))

c("proxymix chi-square" = p_chisq, "proxymix subsampling" = p_andrews,
  "strucchange sup-F" = p_supf, "strucchange OLS-CUSUM" = p_cusum,
  "changepoint cpt.mean" = p_cpt, "KFAS forecast interval" = p_kfas,
  "dlm forecast interval" = p_dlm, "t-test on forecast residuals" = p_t)
```

The [extended version of this
article](https://max578.github.io/proxymix/articles/extended/end_of_sample.html)
gives the full simulation, explains each check's size in detail, and
applies all eight checks to the Nile river flow series, which changed level
after 1898.

## Interpretation

On the stable series the statistic is
`r fixed(test_stable$statistic, 3)`, with a subsampling p-value of
`r fixed(test_stable$p_value, 3)`. The test does not reject, because the
last value lies where the filter expected it. Moving that value up by five
noise standard deviations raises the statistic to
`r fixed(test_break$statistic, 1)` and lowers the p-value to
`r fixed(test_break$p_value, 4)`. The test rejects.

In the first figure the broken value lies
`r fixed(sqrt(test_break$statistic), 1)` forecast standard deviations from
its forecast. The square of that distance is the statistic. This
distance differs from five for two reasons. The forecast
standard deviation is larger than the noise standard deviation, because it
also includes the step of the level and the filter's uncertainty about the
level. The stable value was also not exactly at its forecast. Of the
`r n - 1L` stable values from the second time step on, `r outside`
(`r fixed(100 * outside / (n - 1L), 1)` per cent) fall outside the band.
Under the model, a band of two standard deviations leaves out about
`r fixed(100 * 2 * pnorm(-2), 1)` per cent. In the second figure the broken
value's statistic is larger than every one of the
`r length(test_break$in_sample_blocks)` earlier block statistics. Its p-value is therefore the smallest possible,
$1 / (n - 2m + 2) = `r fixed(1 / (n - 2L + 2L), 4)`$.

That smallest p-value limits the subsampling calibration on short series.
At $n = 30$ and $m = 1$ it is $1 / 30$, so the test rejects at the 5% level
only when the last value is the most extreme of the series. Its size with
known parameters is then `r fixed(grid_30, 4)` rather than 0.05. The
simulated size, `r fixed(size_andrews_30[["size"]], 3)` with a standard
error of `r fixed(size_andrews_30[["se"]], 3)`, agrees with this value. The
chi-square calibration has no smallest p-value. Its simulated size at
$n = 30$ is `r fixed(size_chisq_30[["size"]], 3)` with a standard error of
`r fixed(size_chisq_30[["se"]], 3)`, which is not detectably different from
0.05.

Which calibration to use depends on whether the parameters are known. With
known parameters, the default `method = "chisq"` is exact when the noise is
normal. The parameters are usually estimated from the series being tested. The
comparison above used the default in that setting. At a single final
value its size was `r size_at("proxymix chi-square", 1)` rather than 0.05.
For a single final value with estimated parameters, use
`method = "andrews"`. It was the only check in the comparison whose size stayed close
to 0.05, and it detected the change slightly less often. It needs a series long enough that $1 / (n - 2m + 2)$ is well below
0.05. For windows of 3 or 5 values, neither calibration held its size in
the comparison. At $m = 3$ their sizes were
`r size_at("proxymix chi-square", 3)` and
`r size_at("proxymix subsampling", 3)`.

## Limitations

The size study in this article uses `r n_rep` series per calibration. A
size near 0.05 then has a standard error of about
`r fixed(sqrt(0.05 * 0.95 / n_rep), 3)`, so it cannot detect a departure
from 0.05 of one or two percentage points.

The simulations use normal noise throughout. The subsampling calibration
does not assume normal noise, and the chi-square calibration does. Neither
this article nor its extended version tests how the two calibrations behave when the noise has heavier
tails than the normal. The comparison covers one autoregressive model, one
series length, one size of change and a change confined to the tested
values. It does not cover smaller or gradual changes, longer windows or
other models.

The test applies to a linear state-space model with normal noise, given in
the same form as for `gmm_filter()`: a single-component `gmm` prior, a
`dynamics` list and a `measurement` list. Process or measurement noise
given as a mixture of normal distributions is rejected, because neither
calibration is defined for it. The window `m` must be smaller than the
series length, and the test is designed for small windows such as 1, 2
or 3.

The test checks the last values against the model. It does not check the
model. If the noise variances are estimated from the same short series and
come out too small, the test rejects too often.

## Further reading

*The closed-form operator calculus on a mixture* explains the forecast and
update steps that `gmm_filter()` and this test run, and covers the Kalman
filter in more detail.

*Reading the entropy of a fitted mixture* covers summaries of a fitted
mixture once it is in hand.

*Fitting a proxy to a density you cannot sample* introduces the package's
main task, fitting a mixture of normal distributions to a density.

## References

Andrews, D. W. K. (1993). *Tests for parameter instability and structural
change with unknown change point.* Econometrica 61(4), 821--856.
<https://doi.org/10.2307/2951764>.

Andrews, D. W. K. (2003). *End-of-sample instability tests.* Econometrica
71(6), 1661--1694. <https://doi.org/10.1111/1468-0262.00466>.

Chow, G. C. (1960). *Tests of equality between sets of coefficients in two
linear regressions.* Econometrica 28(3), 591--605.
<https://doi.org/10.2307/1910133>.

Helske, J. (2017). *KFAS: Exponential family state space models in R.*
Journal of Statistical Software 78(10), 1--39.
<https://doi.org/10.18637/jss.v078.i10>.

Killick, R. and Eckley, I. A. (2014). *changepoint: An R package for
changepoint analysis.* Journal of Statistical Software 58(3), 1--19.
<https://doi.org/10.18637/jss.v058.i03>.

Petris, G. (2010). *An R package for dynamic linear models.* Journal of
Statistical Software 36(12), 1--16. <https://doi.org/10.18637/jss.v036.i12>.

Zeileis, A., Leisch, F., Hornik, K. and Kleiber, C. (2002). *strucchange:
An R package for testing for structural change in linear regression
models.* Journal of Statistical Software 7(2), 1--38.
<https://doi.org/10.18637/jss.v007.i02>.

## Reproduce

Running the code with `set.seed(20260621)`, as above, reproduces the
example, the figures and the size study exactly. The comparison is read from stored results of a simulation
run under proxymix `r res$proxymix_version`, strucchange
`r res$versions[["strucchange"]]`, changepoint
`r res$versions[["changepoint"]]`, KFAS `r res$versions[["KFAS"]]` and dlm
`r res$versions[["dlm"]]`, which took about
`r round(res$elapsed_secs / 60)` minutes on one core.

```{r session-info, collapse = FALSE, class.output = "session-info"}
sessionInfo()
```