# Reproduce Example 1.13 (Workbook 1) in its original random-draw order.
sv_record <- function(seed = 542, n = 400L, delta = 0.02, mu = -1,
                      phi = 0.95, sigma_eta = 0.3) {
  stopifnot(n >= 1L, n == as.integer(n), abs(phi) < 1, sigma_eta > 0)
  set.seed(seed)
  h <- numeric(n)
  h[1L] <- rnorm(1L, mean = mu, sd = sigma_eta / sqrt(1 - phi^2))
  eta <- rnorm(n, sd = sigma_eta)
  z <- rnorm(n)
  for (t in seq_len(n)[-1L]) {
    h[t] <- mu + phi * (h[t - 1L] - mu) + eta[t]
  }
  delta + exp(h / 2) * z
}

# Integrate over a uniform grid of log variances. Column t of next_weight
# gives weights for h[t + 1] after updating with returns x[1], ..., x[t].
sv_filter <- function(x, delta = 0.02, mu = -1, phi = 0.95, sigma_eta = 0.3,
                      nodes = 241L, span = 7) {
  stopifnot(is.numeric(x), length(x) > 0L, all(is.finite(x)), abs(phi) < 1,
            sigma_eta > 0, nodes >= 3L, nodes == as.integer(nodes), span > 0)
  state_sd <- sigma_eta / sqrt(1 - phi^2)
  h <- mu + state_sd * seq(-span, span, length.out = nodes)
  observation_sd <- exp(h / 2)
  predicted <- dnorm(h, mean = mu, sd = state_sd)
  predicted <- predicted / sum(predicted)
  # Each column is the discrete transition from one current grid point.
  step <- outer(h, mu + phi * (h - mu),
                function(next_h, center) dnorm(next_h, center, sigma_eta))
  step <- sweep(step, 2L, colSums(step), "/")
  next_weight <- matrix(0, nrow = nodes, ncol = length(x))
  log_score <- numeric(length(x))
  for (t in seq_along(x)) {
    log_joint <- log(predicted) + dnorm(x[t], delta, observation_sd, log = TRUE)
    peak <- max(log_joint)
    unnormalized <- exp(log_joint - peak)
    log_score[t] <- peak + log(sum(unnormalized))
    filtered <- unnormalized / sum(unnormalized)
    predicted <- drop(step %*% filtered)
    next_weight[, t] <- predicted
  }
  list(h = h, delta = delta, next_weight = next_weight, log_score = log_score,
       transition = step)
}

# Probability of an absolute deviation greater than threshold on the next day.
sv_tail <- function(fit, threshold = 1) {
  stopifnot(length(threshold) == 1L, is.finite(threshold), threshold >= 0)
  state_tail <- 2 * pnorm(-threshold / exp(fit$h / 2))
  drop(crossprod(state_tail, fit$next_weight))
}

sv_history_plot <- function(x, origins, delta = 0.02, threshold = 1) {
  stopifnot(length(origins) == 2L, min(origins) >= 10L, max(origins) <= length(x),
            length(threshold) == 1L, is.finite(threshold), threshold > 0)
  old <- par(mfrow = c(1, 2), mar = c(3.7, 3.8, 2.2, 0.8),
             mgp = c(2.3, 0.65, 0), las = 1, bty = "l")
  on.exit(par(old))
  colors <- c("#1C5276", "#A75226")
  days <- lapply(origins, function(t) seq(t - 9L, t))
  band <- delta + c(-threshold, threshold)
  y <- range(x[unlist(days)], band)
  for (i in seq_along(days)) {
    tt <- days[[i]]
    plot(tt, x[tt], type = "n", ylim = y, xlab = "Day", ylab = "Return",
         xaxt = "n", main = paste("Through day", origins[i]), cex.main = 1)
    bounds <- par("usr")
    rect(bounds[1], band[1], bounds[2], band[2], col = "#E9EFF2", border = NA)
    axis(1, at = tt[c(1, 4, 7, 10)])
    abline(h = delta, col = "#A7AFB4", lty = 2)
    lines(tt, x[tt], type = "o", pch = 16, lwd = 1.65, cex = 0.7,
          col = colors[i])
  }
}

# Plot forecasts from two separately supplied histories, using only the final
# state prediction of each filter. Both fits refer to the same model.
sv_forecast_plot <- function(fits, threshold = 1) {
  stopifnot(length(fits) == 2L, threshold > 0,
            identical(fits[[1]]$h, fits[[2]]$h),
            identical(fits[[1]]$delta, fits[[2]]$delta))
  origins <- vapply(fits, function(fit) ncol(fit$next_weight), integer(1))
  h <- fits[[1]]$h
  delta <- fits[[1]]$delta
  weights <- vapply(fits, function(fit) fit$next_weight[, ncol(fit$next_weight)],
                    numeric(length(h)))
  old <- par(mfrow = c(1, 2), mar = c(3.7, 3.8, 2.2, 0.8),
             mgp = c(2.3, 0.65, 0), las = 1, bty = "l")
  on.exit(par(old))
  colors <- c("#1C5276", "#A75226")
  scale <- exp(h / 2)
  sd_next <- sqrt(drop(crossprod(scale^2, weights)))
  half_width <- max(3 * sd_next, 1.5 * threshold)
  xx <- delta + sort(unique(c(seq(-half_width, half_width, length.out = 401L),
                              -threshold, threshold)))
  component <- outer(xx, scale, function(value, sd) dnorm(value, delta, sd))
  density <- component %*% weights
  ylim <- c(0, 1.04 * max(density))
  tails <- list(which(xx <= delta - threshold), which(xx >= delta + threshold))
  for (i in seq_along(fits)) {
    plot(xx, density[, i], type = "n", ylim = ylim, yaxs = "i",
         xlab = "Next return", ylab = "Predictive density",
         main = paste("Forecast for day", origins[i] + 1L), cex.main = 1)
    for (part in tails) {
      polygon(c(xx[part[1L]], xx[part], xx[tail(part, 1)]),
              c(0, density[part, i], 0), border = NA, col = adjustcolor(colors[i], 0.47))
    }
    abline(v = delta + c(-1, 1) * threshold, col = "#939BA1", lty = 3)
    lines(xx, density[, i], col = colors[i], lwd = 2.1)
  }
}