# ==========================================================
# Hotelling T^2 chart  Phase I and Phase II (Example 7.1, Qiu)
# 3-dimensional process, 30 observations (Table 7.1)
# ==========================================================

# --- Data (Table 7.1, p = 3) ---
X1 <- c(-0.224, -0.082, -0.436,  0.455,  0.437,  0.669, -0.186,  0.672,
        -0.121, -0.683,  0.063, -1.029, -0.080,  0.943,  0.147,  0.217,
        -0.611,  0.473, -1.130,  1.379,  0.671,  1.467,  1.077,  2.126,
        2.237,  0.362,  2.119,  1.147,  0.667,  1.944)

X2 <- c(-0.464, -0.203,  0.342,  0.715,  0.237, -0.002, -0.480,  0.384,
        -0.812,  0.093,  0.193,  0.503, -0.238,  0.272,  0.655, -0.260,
        -0.048,  0.066, -1.214,  2.443, -0.466,  0.945,  0.949,  3.964,
        0.959,  0.344,  1.569,  1.003,  0.525,  0.935)

X3 <- c(-0.662,  0.682, -0.174,  1.229, -0.377,  0.230, -0.907,  0.903,
        -1.279, -0.436, -0.028,  0.728, -0.218,  0.831, -0.543, -0.005,
        -0.428,  0.890, -0.064,  2.350,  1.120,  0.402, -0.085,  3.775,
        0.175,  0.375,  1.031,  1.417,  1.676, -0.130)

X <- cbind(X1, X2, X3)
M <- nrow(X)       # 30
p <- 3
alpha <- 0.005

# ==========================================================
# Phase I  all 30 observations used for estimation
# T^2_{1,i} = (X_i - Xbar)' S^{-1} (X_i - Xbar)
# ==========================================================

Xbar  <- colMeans(X)
S     <- cov(X)
S_inv <- solve(S)

T2_phase1 <- numeric(M)
for (i in 1:M) {
  d <- X[i, ] - Xbar
  T2_phase1[i] <- as.numeric(t(d) %*% S_inv %*% d)
}

# Control limit (Beta quantile)
cl_phase1 <- (M - 1)^2 / M * qbeta(1 - alpha, p / 2, (M - p - 1) / 2)

# ==========================================================
# Phase II  first 20 observations IC, last 10 monitored
# T^2_{2,n} = (X_n - Xbar_IC)' S_IC^{-1} (X_n - Xbar_IC)
# ==========================================================

n_ic  <- 20
n_ph2 <- M - n_ic   # 10

Xbar_ic  <- colMeans(X[1:n_ic, ])
S_ic     <- cov(X[1:n_ic, ])
S_ic_inv <- solve(S_ic)

T2_phase2 <- numeric(n_ph2)
for (j in 1:n_ph2) {
  idx <- n_ic + j
  d   <- X[idx, ] - Xbar_ic
  T2_phase2[j] <- as.numeric(t(d) %*% S_ic_inv %*% d)
}

# Control limit (F quantile)
cl_phase2 <- p * (n_ic - 1) * (n_ic + 1) / ((n_ic - p) * n_ic) *
  qf(1 - alpha, p, n_ic - p)

# --- OOC points ---
ooc1 <- T2_phase1 > cl_phase1
ooc2 <- T2_phase2 > cl_phase2

cat("Phase I:  first signal at i =",
    ifelse(any(ooc1), min(which(ooc1)), NA), "\n")
cat("Phase II: first signal at n =",
    ifelse(any(ooc2), min(which(ooc2)), NA), "\n")

# --- Colours ---
ic_col  <- "dodgerblue2"
ooc_col <- "firebrick2"

# ==========================================================
# Plot: 2 rows x 4 columns
# Row 1: X_i1 | X_i2 | X_i3 | Phase I T^2   (all 30 obs)
# Row 2: X_n1 | X_n2 | X_n3 | Phase II T^2  (Phase II obs only)
# ==========================================================

pdf("figure_hotelling.pdf", width = 14, height = 6.5)
par(mfrow = c(2, 4), mar = c(4, 4, 2.5, 1))

# --- Row 1: Phase I (i = 1..30) ---
ii <- 1:M

plot(ii, X1, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "i", ylab = "", main = expression(X[i1]))

plot(ii, X2, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "i", ylab = "", main = expression(X[i2]))

plot(ii, X3, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "i", ylab = "", main = expression(X[i3]))

plot(ii, T2_phase1, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "i", ylab = expression(T[list(1,i)]^2),
     main = expression("Phase I  " * T[list(1,i)]^2))
abline(h = cl_phase1, lty = 2)
points(ii[ooc1], T2_phase1[ooc1], pch = 16, col = ooc_col)

# --- Row 2: Phase II (Phase II observations, plotted as n = 1..10) ---
nn     <- 1:n_ph2
X1_ph2 <- X1[(n_ic + 1):M]
X2_ph2 <- X2[(n_ic + 1):M]
X3_ph2 <- X3[(n_ic + 1):M]

plot(nn, X1_ph2, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n1]))

plot(nn, X2_ph2, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n2]))

plot(nn, X3_ph2, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n3]))

plot(nn, T2_phase2, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = expression(T[list(2,n)]^2),
     main = expression("Phase II  " * T[list(2,n)]^2))
abline(h = cl_phase2, lty = 2)
points(nn[ooc2], T2_phase2[ooc2], pch = 16, col = ooc_col)

dev.off()
cat("Saved figure: figure_hotelling.pdf\n")

# ==========================================================
# Crosier's MCUSUM chart (Crosier, 1988)
# Simulated bivariate process with mean shift at n = 21
# ==========================================================

set.seed(42)

# --- Simulation setup ---
p       <- 2
n_total <- 40
tau     <- 21               # shift starts here

mu0 <- c(0, 0)
mu1 <- c(1, 1)

Sigma0 <- matrix(c(
  1.0, 0.6,
  0.6, 1.0
), nrow = 2, byrow = TRUE)

Sigma0_inv <- solve(Sigma0)

# --- Generate bivariate normal data  X_n ~ N_2(mu_n, Sigma0) ---
L <- t(chol(Sigma0))        # lower triangular,  Sigma0 = L L'

X <- matrix(0, nrow = n_total, ncol = p)
for (i in 1:n_total) {
  mu_i   <- if (i < tau) mu0 else mu1
  X[i, ] <- mu_i + L %*% rnorm(p)
}

# ==========================================================
# Crosier's MCUSUM recursion
#
#   Y_n = [ (U_{n-1} + X_n - mu0)' Sigma0^{-1} (.) ]^{1/2}
#   U_n = 0                                    if Y_n <= k
#   U_n = (U_{n-1} + X_n - mu0) * (1 - k/Y_n)  otherwise
#   C_n = (U_n' Sigma0^{-1} U_n)^{1/2}
#   Signal when C_n > h
# ==========================================================

k <- 0.5
h <- 5.5                    # p = 2, k = 0.5, ARL0 ~ 200

U      <- matrix(0, nrow = n_total, ncol = p)
C      <- numeric(n_total)
U_prev <- rep(0, p)

for (i in 1:n_total) {
  v <- U_prev + X[i, ] - mu0
  Y <- sqrt(as.numeric(t(v) %*% Sigma0_inv %*% v))
  
  if (Y <= k) {
    U[i, ] <- rep(0, p)
  } else {
    U[i, ] <- v * (1 - k / Y)
  }
  
  C[i]   <- sqrt(as.numeric(t(U[i, ]) %*% Sigma0_inv %*% U[i, ]))
  U_prev <- U[i, ]
}

# --- First signal ---
ooc <- C > h
cat("Signal at n =",
    ifelse(any(ooc), min(which(ooc)), NA), "\n")
cat("Mahalanobis length of shift =",
    round(sqrt(as.numeric(
      t(mu1 - mu0) %*% Sigma0_inv %*% (mu1 - mu0))), 3), "\n")

# --- Colours ---
ic_col  <- "dodgerblue2"
ooc_col <- "firebrick2"

# ==========================================================
# Plot layout
#   Top row:    [ X_1 panel ] [ X_2 panel ]
#   Bottom row: [    MCUSUM chart (full width)    ]
# ==========================================================

pdf("figure_crosier.pdf", width = 10.5, height = 6.5)
layout(matrix(c(1, 2,
                3, 3), nrow = 2, byrow = TRUE),
       heights = c(1, 1.2))
par(mar = c(4, 4, 2.5, 1))

nn <- 1:n_total

# --- Top-left: X_1 time series ---
plot(nn, X[, 1], type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n1]))
abline(v = tau - 0.5, lty = 3, col = "grey50")

# --- Top-right: X_2 time series ---
plot(nn, X[, 2], type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n2]))
abline(v = tau - 0.5, lty = 3, col = "grey50")

# --- Bottom: Crosier's MCUSUM C_n with control limit h ---
plot(nn, C, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = expression(C[n]),
     main = expression("Crosier's MCUSUM  " * C[n]))
abline(h = h,         lty = 2)
abline(v = tau - 0.5, lty = 3, col = "grey50")
points(nn[ooc], C[ooc], pch = 16, col = ooc_col)

dev.off()
cat("Saved figure: figure_crosier.pdf\n")



# ==========================================================
# Multivariate CPD chart (Zamba & Hawkins, 2006)
# 3-dimensional process, IC parameters unknown
# Qiu (2014), Example 7.8 / Figure 7.11
# ==========================================================

set.seed(42)

# --- IC parameters (unknown to the chart, used only to generate data) ---
p   <- 3
mu0 <- c(0, 0, 0)

Sigma0 <- matrix(c(
  1.0, 0.8, 0.5,
  0.8, 1.0, 0.8,
  0.5, 0.8, 1.0
), nrow = 3, byrow = TRUE)

# --- Generate 50 observations ---
# First 30 IC, then mean shift to mu1 = (1, 0, 0)' at n = 31
n_total    <- 50
shift_time <- 31
mu1        <- c(1, 0, 0)

L <- t(chol(Sigma0))         # lower-triangular, Sigma0 = L L'

X <- matrix(0, nrow = n_total, ncol = p)
for (i in 1:n_total) {
  mu_i   <- if (i < shift_time) mu0 else mu1
  X[i, ] <- mu_i + L %*% rnorm(p)
}

# ==========================================================
# CPD chart  (Eqs. 7.54 - 7.57 in Qiu)
#
# For each new observation n, scan all candidate change points
# r = 1,..,n-1 and take the maximum T^2_r.
#
#   Xbar_{0,r} = mean of first r observations
#   Xbar_{0,n} = running mean at time n
#   W_n        = (n-1) * sample covariance of first n observations
#   c_r        = r * n / (n - r)
#   quad       = (Xbar_{0,r} - Xbar_{0,n})' W_n^{-1} (Xbar_{0,r} - Xbar_{0,n})
#   T^2_r      = (n - 2) * c_r * quad / (1 - c_r * quad)
#
# Signal when T^2_max,n > h_n  (control limit from eq. 7.56)
# ==========================================================

n_start <- 26                 # monitoring begins here (first 25 are IC)
alpha   <- 0.005              # => ARL0 ~ 200

# Control limit approximation (Eq. 7.56) for alpha = 0.005
compute_hn <- function(n, p) {
  exp(2.706 + 0.230 * p - (p + 3) / 50 * log(n - 25))
}

T2_max <- numeric(n_total)
r_hat  <- integer(n_total)
hn     <- numeric(n_total)

for (n in n_start:n_total) {
  
  # Running statistics up to time n
  Xbar_n <- colMeans(X[1:n, ])
  W_n    <- (n - 1) * cov(X[1:n, ])     # sum of squared deviations
  Wn_inv <- solve(W_n)
  
  # Scan all candidate change points r = 1, .., n-1
  best_T2 <- -Inf
  best_r  <- 1
  
  for (r in 1:(n - 1)) {
    Xbar_0r <- colMeans(X[1:r, , drop = FALSE])
    d_r     <- Xbar_0r - Xbar_n
    c_r     <- r * n / (n - r)
    
    quad <- as.numeric(t(d_r) %*% Wn_inv %*% d_r)
    if (1 - c_r * quad > 0) {
      T2_r <- (n - 2) * c_r * quad / (1 - c_r * quad)
    } else {
      T2_r <- Inf
    }
    
    if (T2_r > best_T2) {
      best_T2 <- T2_r
      best_r  <- r
    }
  }
  
  T2_max[n] <- best_T2
  r_hat[n]  <- best_r
  hn[n]     <- compute_hn(n, p)
}

# --- First signal ---
idx          <- n_start:n_total
signals      <- which(T2_max[idx] > hn[idx]) + n_start - 1
first_signal <- ifelse(length(signals) > 0, signals[1], NA)

cat("First signal at n =", first_signal, "\n")
if (!is.na(first_signal)) {
  cat("Estimated change point r_hat =", r_hat[first_signal],
      " (tau_hat =", r_hat[first_signal] + 1, ")\n")
}

# --- Colours ---
ic_col  <- "dodgerblue2"
ooc_col <- "firebrick2"

# ==========================================================
# Plot layout
#   Top row:    [ X_1 ] [ X_2 ] [ X_3 ]
#   Bottom row: [   CPD chart T^2_max  (full width)  ]
# ==========================================================

pdf("figure_cpd.pdf", width = 14, height = 7)
layout(matrix(c(1, 2, 3,
                4, 4, 4), nrow = 2, byrow = TRUE),
       heights = c(1, 1.2))
par(mar = c(4, 4, 2.5, 1))

nn <- 1:n_total

# --- Top row: three components of X_n ---
plot(nn, X[, 1], type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n1]))
abline(v = shift_time - 0.5, lty = 3, col = "grey50")

plot(nn, X[, 2], type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n2]))
abline(v = shift_time - 0.5, lty = 3, col = "grey50")

plot(nn, X[, 3], type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = "", main = expression(X[n3]))
abline(v = shift_time - 0.5, lty = 3, col = "grey50")

# --- Bottom: CPD chart  T^2_max,n  with time-varying limit h_n ---
T2  <- T2_max[idx]
h_n <- hn[idx]
ooc <- T2 > h_n

plot(idx, T2, type = "b", pch = 16, col = ic_col, bty = "l",
     xlab = "n", ylab = expression(T[list(max,n)]^2),
     main = "CPD chart")
lines(idx, h_n, lty = 2)                  # time-varying control limit
abline(v = shift_time - 0.5, lty = 3, col = "grey50")
points(idx[ooc], T2[ooc], pch = 16, col = ooc_col)

dev.off()
cat("Saved figure: figure_cpd.pdf\n")