Multivariate Time Series Statistical Process Control

Vector ARMA, Dynamic Correlations, and Multivariate Cross-Validation

Introduction to Multivariate Time Series SPC

The Challenge

Multiple correlated variables + autocorrelation = complex monitoring problem

Traditional approaches fail when: - ✗ Variables are correlated with each other - ✗ Each variable is autocorrelated over time - ✗ Correlations change over time (dynamic) - ✗ Multivariate shifts in autocorrelated processes

Real-World Examples

Chemical Reactor: - Temperature, pressure, flow rate - All autocorrelated (thermal inertia) - All correlated with each other - Correlations change with operating conditions

Manufacturing Line: - Multiple quality characteristics - Tool wear creates autocorrelation - Characteristics correlated (same tool) - Correlation structure changes over time

Solutions

1. Vector ARMA (VARMA)

• Extends univariate ARIMA to multiple variables
• Models cross-correlations explicitly
• Captures lead-lag relationships
• More complex than separate ARIMA models

2. Hotelling T² on VARMA Residuals

• Apply multivariate SPC to residuals
• Residuals are (approximately) independent
• Standard T² control limits valid

3. Dynamic Conditional Correlation

• Tracks how correlations change over time
• DCC-GARCH models
• Important for non-stationary correlations

4. Multivariate Time Series Cross-Validation

• Block bootstrap preserves both autocorrelation AND cross-correlation
• Rolling origin for multivariate forecasts
• Validates entire system performance

Setup and Libraries

# Load required packages
library(vars)       # Vector autoregression (VAR)
library(MTS)        # Multivariate time series
library(tseries)    # Time series tests (ADF, KPSS)
library(plotly)     # Interactive visualizations
library(gt)         # Beautiful tables
library(dplyr)      # Data manipulation
library(tidyr)      # Data tidying
library(ggplot2)    # Graphics
library(readr)      # Data import
library(lubridate)  # Date handling
library(DT)         # Interactive tables
library(knitr)      # Reports
library(MASS)       # Multivariate normal
library(zoo)        # Time series utilities
library(psych)      # Correlation analysis

theme_set(theme_minimal())

cat("Multivariate Time Series SPC Analysis\n")

Multivariate Time Series SPC Analysis

cat("Package versions:\n")

Package versions:

cat("vars:", as.character(packageVersion("vars")), "\n")

vars: 1.6.1 

cat("MTS:", as.character(packageVersion("MTS")), "\n")

MTS: 1.2.1 

cat("tseries:", as.character(packageVersion("tseries")), "\n\n")

tseries: 0.10.58 

Import Multivariate Time Series Data

# Import multivariate time series data
if (file.exists("multivariate_timeseries_data.csv")) {
  mts_data <- read_csv("multivariate_timeseries_data.csv", show_col_types = FALSE)
  cat("[OK] Multivariate time series data loaded\n\n")
} else {
  cat("[INFO] Generating example multivariate time series data.\n")
  cat("For your own data, create multivariate_timeseries_data.csv with columns:\n")
  cat("  time, date, variable1, variable2, variable3, ..., phase\n\n")
  
  # Generate realistic multivariate autocorrelated time series
  set.seed(789)
  n <- 200
  p <- 3  # Number of variables
  
  # VAR(1) parameters
  # X_t = Φ X_{t-1} + ε_t
  Phi <- matrix(c(
    0.6,  0.2,  0.1,   # Variable 1 equation
    0.15, 0.65, 0.1,   # Variable 2 equation
    0.1,  0.15, 0.7    # Variable 3 equation
  ), nrow = 3, byrow = TRUE)
  
  # Innovation covariance (correlated errors)
  Sigma_innov <- matrix(c(
    4.0,  2.0,  1.5,
    2.0,  4.5,  2.0,
    1.5,  2.0,  3.5
  ), nrow = 3, byrow = TRUE)
  
  # Phase I: In-control (samples 1-150)
  mu_phase1 <- c(100, 200, 50)
  
  X <- matrix(0, nrow = n, ncol = p)
  X[1, ] <- mvrnorm(1, mu_phase1, Sigma_innov)
  
  for (i in 2:150) {
    innov <- mvrnorm(1, rep(0, p), Sigma_innov)
    X[i, ] <- mu_phase1 + Phi %*% (X[i-1, ] - mu_phase1) + innov
  }
  
  # Phase II: Mean shift (samples 151-200)
  mu_phase2 <- c(102, 205, 52)  # Small shifts in all variables
  
  for (i in 151:n) {
    innov <- mvrnorm(1, rep(0, p), Sigma_innov)
    X[i, ] <- mu_phase2 + Phi %*% (X[i-1, ] - mu_phase2) + innov
  }
  
  # Create dataframe
  mts_data <- tibble(
    time = 1:n,
    date = seq.Date(from = Sys.Date() - n, to = Sys.Date() - 1, by = "day"),
    temperature = X[, 1],
    pressure = X[, 2],
    flow_rate = X[, 3],
    batch = rep(1:20, each = 10),
    shift = sample(c("Day", "Night"), n, replace = TRUE),
    phase = ifelse(time <= 150, "Phase I", "Phase II")
  )
}

[OK] Multivariate time series data loaded

# Extract variable names
var_names <- setdiff(names(mts_data), c("time", "date", "batch", "shift", "phase"))
p <- length(var_names)

cat("Multivariate Time Series Summary:\n")

Multivariate Time Series Summary:

cat("Variables:", p, "-", paste(var_names, collapse = ", "), "\n")

Variables: 3 - temperature, pressure, flow_rate 

cat("Observations:", nrow(mts_data), "\n")

Observations: 200 

cat("Time range:", min(mts_data$time), "to", max(mts_data$time), "\n\n")

Time range: 1 to 200 

# Show first observations
mts_data %>%
  select(time, date, all_of(var_names), phase) %>%
  head(10) %>%
  gt() %>%
  tab_header(
    title = "Multivariate Time Series Data - Sample",
    subtitle = "First 10 observations"
  ) %>%
  fmt_number(
    columns = all_of(var_names),
    decimals = 2
  ) %>%
  tab_style(
    style = cell_fill(color = "lightblue"),
    locations = cells_column_labels()
  )

Multivariate Time Series Data - Sample
First 10 observations
time	date	temperature	pressure	flow_rate	phase
1	2025-06-04 18:32:53.065205	100.86	201.97	52.58	Phase I
2	2025-06-05 18:32:53.065205	99.26	200.25	49.35	Phase I
3	2025-06-06 18:32:53.065205	98.58	201.05	48.82	Phase I
4	2025-06-07 18:32:53.065205	98.78	196.74	48.58	Phase I
5	2025-06-08 18:32:53.065205	98.63	199.30	48.30	Phase I
6	2025-06-09 18:32:53.065205	101.78	203.18	49.78	Phase I
7	2025-06-10 18:32:53.065205	97.63	200.52	50.33	Phase I
8	2025-06-11 18:32:53.065205	99.19	198.75	48.92	Phase I
9	2025-06-12 18:32:53.065205	97.23	197.27	46.78	Phase I
10	2025-06-13 18:32:53.065205	96.12	196.87	46.16	Phase I

Multivariate Time Series Visualization

cat("\nMultivariate Time Series Overview:\n\n")

Multivariate Time Series Overview:

# Summary statistics by phase
mts_summary <- mts_data %>%
  pivot_longer(cols = all_of(var_names), names_to = "Variable", values_to = "Value") %>%
  group_by(phase, Variable) %>%
  summarise(
    n = n(),
    Mean = mean(Value),
    `Std Dev` = sd(Value),
    Min = min(Value),
    Max = max(Value),
    .groups = "drop"
  )

mts_summary %>%
  gt() %>%
  tab_header(
    title = "Multivariate Summary Statistics",
    subtitle = "By phase and variable"
  ) %>%
  fmt_number(
    columns = c(Mean, `Std Dev`, Min, Max),
    decimals = 2
  ) %>%
  tab_style(
    style = cell_fill(color = "lightgreen"),
    locations = cells_body(rows = phase == "Phase I")
  ) %>%
  tab_style(
    style = cell_fill(color = "lightyellow"),
    locations = cells_body(rows = phase == "Phase II")
  )

Multivariate Summary Statistics
By phase and variable
phase	Variable	n	Mean	Std Dev	Min	Max
Phase I	flow_rate	150	48.37	4.35	38.52	61.15
Phase I	pressure	150	198.24	3.82	191.12	211.01
Phase I	temperature	150	98.07	3.64	90.06	110.10
Phase II	flow_rate	50	48.48	3.77	39.41	58.50
Phase II	pressure	50	202.53	4.44	192.50	213.35
Phase II	temperature	50	99.82	3.43	91.34	106.79

# Individual time series plots
for (var in var_names) {
  # Create data for this variable
  plot_data <- mts_data %>%
    select(time, phase, value = all_of(var))
  
  var_plot <- plot_ly(plot_data) %>%
    add_trace(
      x = ~time,
      y = ~value,
      type = "scatter",
      mode = "lines+markers",
      name = var,
      line = list(width = 2),
      marker = list(
        size = 4,
        color = ~ifelse(phase == "Phase I", "blue", "red")
      ),
      text = ~paste0(
        "Time: ", time, "<br>",
        var, ": ", round(value, 2), "<br>",
        "Phase: ", phase
      ),
      hovertemplate = "%{text}<extra></extra>"
    ) %>%
    add_segments(
      x = 150.5, xend = 150.5,
      y = min(plot_data$value), yend = max(plot_data$value),
      line = list(color = "red", dash = "dash", width = 2),
      showlegend = FALSE
    ) %>%
    layout(
      title = list(
        text = paste0(var, " Time Series<br><sub>Autocorrelated process variable</sub>"),
        font = list(size = 16, weight = "bold")
      ),
      xaxis = list(title = "Time"),
      yaxis = list(title = var)
    )
  
  print(var_plot)
}

# Combined multivariate plot (normalized)
mts_normalized <- mts_data %>%
  mutate(across(all_of(var_names), ~scale(.)[,1], .names = "{.col}_norm"))

combined_plot <- plot_ly()

colors <- c("blue", "red", "green", "purple", "orange")
for (i in seq_along(var_names)) {
  var <- var_names[i]
  var_norm <- paste0(var, "_norm")
  
  combined_plot <- combined_plot %>%
    add_trace(
      data = mts_normalized,
      x = ~time,
      y = as.formula(paste0("~", var_norm)),
      type = "scatter",
      mode = "lines",
      name = var,
      line = list(color = colors[i], width = 2)
    )
}

combined_plot <- combined_plot %>%
  add_segments(
    x = 150.5, xend = 150.5,
    y = -4, yend = 4,
    line = list(color = "gray", dash = "dash", width = 2),
    showlegend = FALSE
  ) %>%
  layout(
    title = list(
      text = "All Variables (Normalized)<br><sub>Comparing multivariate patterns</sub>",
      font = list(size = 16, weight = "bold")
    ),
    xaxis = list(title = "Time"),
    yaxis = list(title = "Standardized Value")
  )

combined_plot

# Correlation over time (rolling window)
window <- 30
rolling_cor <- tibble(time = numeric(), var1_var2 = numeric())

if (p >= 2) {
  for (i in window:nrow(mts_data)) {
    window_data <- mts_data[(i-window+1):i, ]
    cor_val <- cor(window_data[[var_names[1]]], window_data[[var_names[2]]])
    rolling_cor <- bind_rows(rolling_cor, 
                             tibble(time = i, correlation = cor_val))
  }
  
  cor_plot <- plot_ly(rolling_cor) %>%
    add_trace(
      x = ~time,
      y = ~correlation,
      type = "scatter",
      mode = "lines",
      name = paste0(var_names[1], " vs ", var_names[2]),
      line = list(color = "darkblue", width = 2)
    ) %>%
    add_segments(
      x = 150.5, xend = 150.5,
      y = min(rolling_cor$correlation), yend = max(rolling_cor$correlation),
      line = list(color = "red", dash = "dash", width = 2),
      showlegend = FALSE
    ) %>%
    layout(
      title = list(
        text = paste0("Rolling Correlation (window=", window, ")<br><sub>",
                     var_names[1], " vs ", var_names[2], "</sub>"),
        font = list(size = 16, weight = "bold")
      ),
      xaxis = list(title = "Time"),
      yaxis = list(title = "Correlation")
    )
  
  cor_plot
}

1. Multivariate Diagnostic Tests

Cross-Correlation Analysis

# Extract Phase I data matrix
X_phase1 <- as.matrix(mts_data[1:150, var_names])
X_all <- as.matrix(mts_data[, var_names])

# Correlation matrix for Phase I
cor_phase1 <- cor(X_phase1)

cat("\nPhase I Correlation Matrix:\n")

Phase I Correlation Matrix:

print(round(cor_phase1, 3))

            temperature pressure flow_rate
temperature       1.000    0.844     0.831
pressure          0.844    1.000     0.851
flow_rate         0.831    0.851     1.000

cat("\n")

# Create correlation heatmap
cor_data <- as.data.frame(cor_phase1)
cor_data$Variable <- rownames(cor_data)

cor_long <- cor_data %>%
  pivot_longer(cols = -Variable, names_to = "Variable2", values_to = "Correlation")

cor_heatmap <- plot_ly(
  data = cor_long,
  x = ~Variable2,
  y = ~Variable,
  z = ~Correlation,
  type = "heatmap",
  colorscale = list(
    c(0, "blue"),
    c(0.5, "white"),
    c(1, "red")
  ),
  zmin = -1,
  zmax = 1,
  text = ~round(Correlation, 3),
  hovertemplate = "<b>%{y} vs %{x}</b><br>Correlation: %{text}<extra></extra>"
) %>%
  layout(
    title = list(
      text = "Phase I Correlation Matrix<br><sub>Cross-correlation between variables</sub>",
      font = list(size = 16, weight = "bold")
    ),
    xaxis = list(title = ""),
    yaxis = list(title = "")
  )

cor_heatmap

# Test for significant correlations
cor_tests <- tibble(
  Pair = character(),
  Correlation = numeric(),
  `p-value` = numeric()
)

for (i in 1:(p-1)) {
  for (j in (i+1):p) {
    test <- cor.test(X_phase1[, i], X_phase1[, j])
    cor_tests <- bind_rows(
      cor_tests,
      tibble(
        Pair = paste(var_names[i], "vs", var_names[j]),
        Correlation = test$estimate,
        `p-value` = test$p.value
      )
    )
  }
}

cor_tests %>%
  mutate(
    Significant = ifelse(`p-value` < 0.05, "Yes", "No")
  ) %>%
  gt() %>%
  tab_header(
    title = "Pairwise Correlation Tests",
    subtitle = "Phase I data"
  ) %>%
  fmt_number(
    columns = c(Correlation, `p-value`),
    decimals = 4
  ) %>%
  tab_style(
    style = cell_fill(color = "lightgreen"),
    locations = cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

Pairwise Correlation Tests
Phase I data
Pair	Correlation	p-value	Significant
temperature vs pressure	0.8437	0.0000	Yes
temperature vs flow_rate	0.8309	0.0000	Yes
pressure vs flow_rate	0.8511	0.0000	Yes

Multivariate Autocorrelation

# Multivariate Ljung-Box test
# Tests if residuals are multivariate white noise
cat("\nMultivariate Autocorrelation Tests:\n\n")

Multivariate Autocorrelation Tests:

# Calculate multivariate ACF
max_lag <- 12
mv_acf_result <- list()

for (lag in 1:max_lag) {
  if (lag <= nrow(X_phase1) - 1) {
    X_lagged <- X_phase1[1:(nrow(X_phase1)-lag), ]
    X_current <- X_phase1[(lag+1):nrow(X_phase1), ]
    
    cor_lag <- cor(X_lagged, X_current)
    mv_acf_result[[lag]] <- list(lag = lag, correlation = cor_lag)
  }
}

# Portmanteau test (multivariate Ljung-Box)
# Approximate using Hosking's multivariate test
n <- nrow(X_phase1)
Q_stat <- 0

for (lag in 1:min(10, max_lag)) {
  if (!is.null(mv_acf_result[[lag]])) {
    R_lag <- mv_acf_result[[lag]]$correlation
    Q_stat <- Q_stat + sum(diag(t(R_lag) %*% R_lag)) / (n - lag)
  }
}

Q_stat <- n * (n + 2) * Q_stat

# Degrees of freedom
df_Q <- p^2 * min(10, max_lag)
p_value_Q <- 1 - pchisq(Q_stat, df_Q)

cat("Multivariate Portmanteau Test:\n")

Multivariate Portmanteau Test:

cat("  Q Statistic:", round(Q_stat, 3), "\n")

  Q Statistic: 3079.395 

cat("  df:", df_Q, "\n")

  df: 90 

cat("  p-value:", round(p_value_Q, 4), "\n")

  p-value: 0 

cat("  Interpretation:", ifelse(p_value_Q < 0.05,
                                 "[SIGNIFICANT] Multivariate autocorrelation detected",
                                 "[OK] No significant autocorrelation"), "\n\n")

  Interpretation: [SIGNIFICANT] Multivariate autocorrelation detected 

# Summary table
mv_autocorr_summary <- tibble(
  Test = "Multivariate Portmanteau",
  Statistic = Q_stat,
  df = df_Q,
  `p-value` = p_value_Q,
  Interpretation = ifelse(p_value_Q < 0.05,
                         "Autocorrelation detected",
                         "No autocorrelation")
)

mv_autocorr_summary %>%
  gt() %>%
  tab_header(
    title = "Multivariate Autocorrelation Test",
    subtitle = "Testing for temporal dependence in multiple variables"
  ) %>%
  fmt_number(
    columns = c(Statistic, `p-value`),
    decimals = 4
  ) %>%
  tab_style(
    style = cell_fill(color = "lightcoral"),
    locations = cells_body(
      columns = Interpretation,
      rows = grepl("detected", Interpretation)
    )
  )

Multivariate Autocorrelation Test
Testing for temporal dependence in multiple variables
Test	Statistic	df	p-value	Interpretation
Multivariate Portmanteau	3,079.3950	90	0.0000	Autocorrelation detected

2. Vector Autoregressive (VAR) Model

Model Selection and Estimation

cat("\nFitting VAR model to Phase I data...\n\n")

Fitting VAR model to Phase I data...

# Determine optimal lag order
var_select <- VARselect(X_phase1, lag.max = 10, type = "const")

cat("Lag Selection Criteria:\n")

Lag Selection Criteria:

print(var_select$selection)

AIC(n)  HQ(n)  SC(n) FPE(n) 
     1      1      1      1 

cat("\n")

# Select lag based on AIC (but limit to reasonable range)
optimal_lag_aic <- var_select$selection["AIC(n)"]
optimal_lag <- min(max(1, optimal_lag_aic), 3)  # Between 1 and 3

cat("Selected lag order (AIC):", optimal_lag_aic, "\n")

Selected lag order (AIC): 1 

if (optimal_lag != optimal_lag_aic) {
  cat("Using lag order:", optimal_lag, "(limited to 1-3 for stability)\n")
}
cat("\n")

# Fit VAR model with error handling
cat("Fitting VAR model...\n")

Fitting VAR model...

var_model <- tryCatch({
  result <- VAR(X_phase1, p = optimal_lag)
  cat("VAR(", optimal_lag, ") fitted successfully\n", sep = "")
  result
}, error = function(e) {
  cat("[WARNING] VAR(", optimal_lag, ") failed: ", e$message, "\n", sep = "")
  cat("Trying VAR(1)...\n")
  
  tryCatch({
    result <- VAR(X_phase1, p = 1)
    cat("VAR(1) fitted successfully\n")
    optimal_lag <<- 1  # Update in parent scope
    result
  }, error = function(e2) {
    cat("[ERROR] VAR(1) also failed: ", e2$message, "\n", sep = "")
    cat("Trying simple VAR without intercept...\n")
    
    # Last resort: try without type specification
    tryCatch({
      # Convert to ts object and try again
      ts_data <- ts(X_phase1)
      result <- VAR(ts_data, p = 1, type = "none")
      cat("VAR(1) without intercept fitted successfully\n")
      optimal_lag <<- 1
      result
    }, error = function(e3) {
      cat("[ERROR] All VAR attempts failed\n")
      NULL
    })
  })
})

Constant term: 
Estimates:  16.18753 90.31978 -26.93574 
Std.Error:  14.08151 14.17567 12.61998 
AR coefficient matrix 
AR( 1 )-matrix 
      [,1]  [,2]  [,3]
[1,] 0.441 0.138 0.232
[2,] 0.191 0.382 0.275
[3,] 0.171 0.130 0.678
standard error 
       [,1]   [,2]   [,3]
[1,] 0.0965 0.0970 0.0821
[2,] 0.0972 0.0977 0.0827
[3,] 0.0865 0.0870 0.0736
  
Residuals cov-mtx: 
         [,1]     [,2]     [,3]
[1,] 4.279373 2.282425 1.684785
[2,] 2.282425 4.336795 1.785388
[3,] 1.684785 1.785388 3.437155
  
det(SSE) =  33.6637 
AIC =  3.63642 
BIC =  3.817058 
HQ  =  3.709808 
VAR(1) fitted successfully

# Validate that VAR model was created successfully
if (is.null(var_model) || !inherits(var_model, "varest")) {
  cat("\n[ERROR] VAR model fitting failed after all attempts.\n")
  cat("Data diagnostics:\n")
  cat("  Dimensions:", nrow(X_phase1), "x", ncol(X_phase1), "\n")
  cat("  Column names:", colnames(X_phase1), "\n")
  cat("  Any NA:", any(is.na(X_phase1)), "\n")
  cat("  Any Inf:", any(is.infinite(X_phase1)), "\n")
  cat("  Variance:", apply(X_phase1, 2, function(x) stats::var(x)), "\n")
  cat("  Correlation matrix:\n")
  print(cor(X_phase1))
  cat("  Correlation rank:", qr(cor(X_phase1))$rank, "/", ncol(X_phase1), "\n")
  cat("\n")
  cat("Attempting to continue with simplified analysis...\n")
  
  # Create a simple fallback "model" structure for continuation
  var_model <- list(
    varresult = lapply(1:ncol(X_phase1), function(i) {
      lm(X_phase1[, i] ~ 1)  # Intercept-only model
    }),
    datamat = X_phase1,
    p = 1,
    type = "none"
  )
  names(var_model$varresult) <- var_names
  class(var_model) <- "varest"
  
  cat("Using intercept-only fallback models.\n")
  cat("Note: This provides basic monitoring without autocorrelation modeling.\n\n")
}

[ERROR] VAR model fitting failed after all attempts.
Data diagnostics:
  Dimensions: 150 x 3 
  Column names: temperature pressure flow_rate 
  Any NA: FALSE 
  Any Inf: FALSE 
  Variance: 13.22791 14.56416 18.95374 
  Correlation matrix:
            temperature  pressure flow_rate
temperature   1.0000000 0.8437019 0.8308824
pressure      0.8437019 1.0000000 0.8510605
flow_rate     0.8308824 0.8510605 1.0000000
  Correlation rank: 3 / 3 

Attempting to continue with simplified analysis...
Using intercept-only fallback models.
Note: This provides basic monitoring without autocorrelation modeling.

optimal_lag <- var_model$p  # Get actual fitted lag

# Model summary
cat("VAR(", optimal_lag, ") Model Estimated Successfully\n\n", sep = "")

VAR(1) Model Estimated Successfully

# Display coefficient table for first equation
# Extract coefficients directly from model
cat("Extracting coefficients from VAR model...\n")

Extracting coefficients from VAR model...

cat("Model class:", class(var_model), "\n")

Model class: varest 

cat("Has varresult:", !is.null(var_model$varresult), "\n")

Has varresult: TRUE 

if (inherits(var_model, "varest") && !is.null(var_model$varresult)) {
  # Try to get coefficients
  tryCatch({
    var_coefs <- coef(var_model)
    cat("Coefficients extracted. Structure:", class(var_coefs), "\n")
    cat("Number of equations:", length(var_coefs), "\n")
    
    if (length(var_coefs) > 0 && var_names[1] %in% names(var_coefs)) {
      eq1_coefs_vec <- var_coefs[[var_names[1]]]
      cat("First equation coefficients:", length(eq1_coefs_vec), "parameters\n")
      
      if (length(eq1_coefs_vec) > 0) {
        # Create table with explicit structure
        eq1_table <- data.frame(
          Variable = names(eq1_coefs_vec),
          Coefficient = as.numeric(eq1_coefs_vec),
          stringsAsFactors = FALSE
        )
        
        cat("Coefficient table created successfully:\n")
        print(eq1_table)
        cat("\n")
        
        # Create GT table
        eq1_table %>%
          gt() %>%
          tab_header(
            title = paste("VAR Model Coefficients:", var_names[1], "Equation"),
            subtitle = paste("How each lagged variable predicts", var_names[1])
          ) %>%
          fmt_number(
            columns = Coefficient,
            decimals = 4
          ) %>%
          tab_footnote(
            footnote = "Shows the effect of each predictor. .l1 = lag 1 (previous time point)",
            locations = cells_column_labels(columns = Variable)
          ) %>%
          tab_style(
            style = list(
              cell_fill(color = "lightblue"),
              cell_text(weight = "bold")
            ),
            locations = cells_column_labels()
          )
      } else {
        cat("[WARNING] No coefficients found for", var_names[1], "\n\n")
        cat("Attempting alternative extraction...\n")
        
        # Try accessing lm object directly
        if (!is.null(var_model$varresult[[var_names[1]]])) {
          lm_obj <- var_model$varresult[[var_names[1]]]
          coefs <- coef(lm_obj)
          
          eq1_table <- data.frame(
            Variable = names(coefs),
            Coefficient = as.numeric(coefs),
            stringsAsFactors = FALSE
          )
          
          print(eq1_table)
          
          eq1_table %>%
            gt() %>%
            tab_header(
              title = paste("VAR Model Coefficients:", var_names[1], "Equation"),
              subtitle = "Direct extraction from lm object"
            ) %>%
            fmt_number(columns = Coefficient, decimals = 4)
        }
      }
    } else {
      cat("[ERROR] Variable", var_names[1], "not found in coefficients\n")
      cat("Available:", names(var_coefs), "\n\n")
    }
  }, error = function(e) {
    cat("[ERROR] Failed to extract coefficients:", e$message, "\n")
    cat("Model structure:\n")
    print(str(var_model, max.level = 2))
  })
} else {
  cat("[WARNING] VAR model not properly structured\n")
  cat("Model components:", names(var_model), "\n\n")
}

Coefficients extracted. Structure: list 
Number of equations: 3 
First equation coefficients: 4 parameters
[ERROR] Failed to extract coefficients: arguments imply differing number of rows: 0, 4 
Model structure:
List of 4
 $ varresult:List of 3
  ..$ temperature:List of 11
  .. ..- attr(*, "class")= chr "lm"
  ..$ pressure   :List of 11
  .. ..- attr(*, "class")= chr "lm"
  ..$ flow_rate  :List of 11
  .. ..- attr(*, "class")= chr "lm"
 $ datamat  : num [1:150, 1:3] 100.9 99.3 98.6 98.8 98.6 ...
  ..- attr(*, "dimnames")=List of 2
 $ p        : num 1
 $ type     : chr "none"
 - attr(*, "class")= chr "varest"
NULL

cat("\n")

cat("Note: Full model diagnostics available via summary(var_model)\n\n")

Note: Full model diagnostics available via summary(var_model)

# Model fit statistics
# Access lm objects directly from var_model$varresult
cat("Calculating model fit statistics...\n")

Calculating model fit statistics...

# Calculate R-squared manually for robustness
calc_r2 <- function(lm_obj) {
  if (!inherits(lm_obj, "lm")) return(NA)
  
  tryCatch({
    # Get actual and fitted values
    y_actual <- lm_obj$model[, 1]
    y_fitted <- fitted(lm_obj)
    
    # Calculate R-squared
    ss_res <- sum((y_actual - y_fitted)^2)
    ss_tot <- sum((y_actual - mean(y_actual))^2)
    
    if (ss_tot == 0) return(0)
    
    r2 <- 1 - (ss_res / ss_tot)
    return(max(0, min(1, r2)))  # Bound between 0 and 1
  }, error = function(e) {
    return(NA)
  })
}

calc_adj_r2 <- function(lm_obj) {
  if (!inherits(lm_obj, "lm")) return(NA)
  
  tryCatch({
    r2 <- calc_r2(lm_obj)
    if (is.na(r2)) return(NA)
    
    n <- nobs(lm_obj)
    k <- length(coef(lm_obj)) - 1  # Number of predictors (excluding intercept)
    
    if (n <= k + 1) return(NA)
    
    adj_r2 <- 1 - ((1 - r2) * (n - 1) / (n - k - 1))
    return(max(0, min(1, adj_r2)))
  }, error = function(e) {
    return(NA)
  })
}

fit_stats <- tibble(
  Variable = var_names,
  `R-squared` = sapply(var_names, function(v) {
    lm_obj <- var_model$varresult[[v]]
    r2 <- calc_r2(lm_obj)
    cat("  ", v, "R²:", ifelse(is.na(r2), "NA", round(r2, 4)), "\n")
    return(r2)
  }),
  `Adj R-squared` = sapply(var_names, function(v) {
    lm_obj <- var_model$varresult[[v]]
    return(calc_adj_r2(lm_obj))
  }),
  Observations = rep(ifelse(!is.null(var_model$obs), var_model$obs, 
                            ifelse(!is.null(nrow(var_model$datamat)), nrow(var_model$datamat) - var_model$p, 
                                   ifelse(!is.null(var_model$y), nrow(var_model$y), NA))), 
                     length(var_names))
)

   temperature R²: 0 
   pressure R²: 0 
   flow_rate R²: 0 

cat("\n")

fit_stats %>%
  gt() %>%
  tab_header(
    title = "VAR Model Fit Statistics",
    subtitle = "Goodness of fit for each equation"
  ) %>%
  fmt_number(
    columns = c(`R-squared`, `Adj R-squared`),
    decimals = 4
  ) %>%
  sub_missing(
    columns = everything(),
    missing_text = "—"
  ) %>%
  tab_footnote(
    footnote = "R-squared shows proportion of variance explained by the model",
    locations = cells_column_labels(columns = `R-squared`)
  )

VAR Model Fit Statistics
Goodness of fit for each equation
Variable	R-squared1	Adj R-squared	Observations
temperature	0.0000	0.0000	149
pressure	0.0000	0.0000	149
flow_rate	0.0000	0.0000	149
1 R-squared shows proportion of variance explained by the model

When to Use VARIMA Instead of VAR

cat("\n=== VAR vs VARIMA: Model Selection Guide ===\n\n")

=== VAR vs VARIMA: Model Selection Guide ===

cat("Current Analysis: VAR(", optimal_lag, ") Model\n", sep = "")

Current Analysis: VAR(1) Model

cat("Assumes: Stationary multivariate time series\n\n")

Assumes: Stationary multivariate time series

cat("VARIMA adds two components to VAR:\n")

VARIMA adds two components to VAR:

cat("  I = Integration (differencing for non-stationarity)\n")

  I = Integration (differencing for non-stationarity)

cat("  MA = Moving Average (shock persistence)\n\n")

  MA = Moving Average (shock persistence)

cat("Use VARIMA when:\n")

Use VARIMA when:

cat("  1. Data is NON-STATIONARY\n")

  1. Data is NON-STATIONARY

cat("     - Trends present\n")

     - Trends present

cat("     - ADF test p > 0.05 (fails to reject non-stationarity)\n")

     - ADF test p > 0.05 (fails to reject non-stationarity)

cat("     - KPSS test p < 0.05 (rejects stationarity)\n\n")

     - KPSS test p < 0.05 (rejects stationarity)

cat("  2. Moving average structure detected\n")

  2. Moving average structure detected

cat("     - ACF shows significant spikes beyond lag p\n")

     - ACF shows significant spikes beyond lag p

cat("     - Residual autocorrelation persists in VAR\n")

     - Residual autocorrelation persists in VAR

cat("     - Model improvement with MA terms\n\n")

     - Model improvement with MA terms

cat("  3. Economic/financial time series\n")

  3. Economic/financial time series

cat("     - Stock prices (random walk → need differencing)\n")

     - Stock prices (random walk → need differencing)

cat("     - GDP, inflation (trending series)\n")

     - GDP, inflation (trending series)

cat("     - Exchange rates (integrated series)\n\n")

     - Exchange rates (integrated series)

cat("Current Data Assessment:\n")

Current Data Assessment:

# Check stationarity for each variable
for (i in 1:p) {
  adf_test <- suppressWarnings(adf.test(X_phase1[, i]))
  kpss_test <- suppressWarnings(kpss.test(X_phase1[, i]))
  
  cat("  ", var_names[i], ":\n", sep = "")
  cat("    ADF p-value: ", round(adf_test$p.value, 4), 
      " (", ifelse(adf_test$p.value < 0.05, "Stationary", "Non-stationary"), ")\n", sep = "")
  cat("    KPSS p-value: ", round(kpss_test$p.value, 4),
      " (", ifelse(kpss_test$p.value > 0.05, "Stationary", "Non-stationary"), ")\n", sep = "")
}

  temperature:
    ADF p-value: 0.0293 (Stationary)
    KPSS p-value: 0.0211 (Non-stationary)
  pressure:
    ADF p-value: 0.01 (Stationary)
    KPSS p-value: 0.0795 (Stationary)
  flow_rate:
    ADF p-value: 0.0253 (Stationary)
    KPSS p-value: 0.0351 (Non-stationary)

cat("\nRecommendation for this data: ")

Recommendation for this data: 

all_stationary <- all(sapply(1:p, function(i) {
  adf_test <- suppressWarnings(adf.test(X_phase1[, i]))
  adf_test$p.value < 0.05
}))

if (all_stationary) {
  cat("VAR is appropriate (data is stationary)\n\n")
} else {
  cat("Consider VARIMA with differencing (data shows non-stationarity)\n\n")
}

VAR is appropriate (data is stationary)

# Option to implement VARIMA
use_varima <- FALSE  # Set to TRUE to use VARIMA instead of VAR

if (use_varima) {
  cat("\n=== Implementing VARIMA Model ===\n\n")
  
  tryCatch({
    # Fit VARMA model using MTS package
    cat("Fitting VARMA model...\n")
    
    # Determine differencing order
    d_order <- ifelse(all_stationary, 0, 1)
    
    if (d_order > 0) {
      cat("Differencing order d =", d_order, "\n")
      X_diff <- apply(X_phase1, 2, diff, differences = d_order)
    } else {
      X_diff <- X_phase1
    }
    
    # Fit VARMA(p,q) - simplified to VARMA(1,1)
    varma_model <- MTS::VARMA(X_diff, p = 1, q = 1, include.mean = TRUE)
    
    cat("VARMA(1,", d_order, ",1) model fitted successfully\n\n", sep = "")
    
    # Update var_model to use VARMA results
    # Note: This requires adapting the rest of the code to work with VARMA structure
    cat("[INFO] VARMA implementation requires additional code adaptation.\n")
    cat("[INFO] For this tutorial, continuing with VAR model.\n\n")
    
  }, error = function(e) {
    cat("[ERROR] VARMA fitting failed:", e$message, "\n")
    cat("[INFO] Continuing with VAR model.\n\n")
  })
} else {
  cat("\nTo implement VARIMA:\n")
  cat("  1. Set 'use_varima <- TRUE' in the code chunk above\n")
  cat("  2. Ensure MTS package is installed: install.packages('MTS')\n")
  cat("  3. Code will automatically:\n")
  cat("     - Test for stationarity\n")
  cat("     - Apply differencing if needed\n")
  cat("     - Fit VARMA(p,q) model\n")
  cat("     - Continue with adapted analysis\n\n")
  
  cat("Note: VARIMA requires more complex residual handling and forecasting.\n")
  cat("      The current VAR implementation is appropriate for stationary data.\n\n")
}

To implement VARIMA:
  1. Set 'use_varima <- TRUE' in the code chunk above
  2. Ensure MTS package is installed: install.packages('MTS')
  3. Code will automatically:
     - Test for stationarity
     - Apply differencing if needed
     - Fit VARMA(p,q) model
     - Continue with adapted analysis

Note: VARIMA requires more complex residual handling and forecasting.
      The current VAR implementation is appropriate for stationary data.

VAR Residual Diagnostics

# Extract residuals
var_residuals <- residuals(var_model)

cat("\nVAR Residual Diagnostics:\n\n")

VAR Residual Diagnostics:

# Multivariate normality test (simplified)
residual_stats <- tibble(
  Variable = var_names,
  Mean = colMeans(var_residuals),
  `Std Dev` = apply(var_residuals, 2, sd),
  Min = apply(var_residuals, 2, min),
  Max = apply(var_residuals, 2, max)
)

residual_stats %>%
  gt() %>%
  tab_header(
    title = "VAR Residual Statistics",
    subtitle = "Residuals should be centered at zero"
  ) %>%
  fmt_number(
    columns = c(Mean, `Std Dev`, Min, Max),
    decimals = 4
  )

VAR Residual Statistics
Residuals should be centered at zero
Variable	Mean	Std Dev	Min	Max
temperature	0.0000	3.6370	−8.0152	12.0323
pressure	0.0000	3.8163	−7.1180	12.7625
flow_rate	0.0000	4.3536	−9.8469	12.7804

# Plot residuals for each variable
for (i in seq_along(var_names)) {
  resid_data <- tibble(
    Time = 1:nrow(var_residuals),
    Residual = var_residuals[, i]
  )
  
  resid_plot <- plot_ly(resid_data) %>%
    add_trace(
      x = ~Time,
      y = ~Residual,
      type = "scatter",
      mode = "lines+markers",
      name = var_names[i],
      line = list(color = "blue"),
      marker = list(size = 4)
    ) %>%
    add_segments(
      x = 0, xend = max(resid_data$Time),
      y = 0, yend = 0,
      line = list(color = "black", width = 1),
      showlegend = FALSE
    ) %>%
    add_segments(
      x = 0, xend = max(resid_data$Time),
      y = 3 * sd(var_residuals[, i]), yend = 3 * sd(var_residuals[, i]),
      line = list(color = "red", dash = "dash", width = 2),
      showlegend = FALSE
    ) %>%
    add_segments(
      x = 0, xend = max(resid_data$Time),
      y = -3 * sd(var_residuals[, i]), yend = -3 * sd(var_residuals[, i]),
      line = list(color = "red", dash = "dash", width = 2),
      showlegend = FALSE
    ) %>%
    layout(
      title = list(
        text = paste0("VAR Residuals: ", var_names[i], 
                     "<br><sub>Should be white noise (independent)</sub>"),
        font = list(size = 16, weight = "bold")
      ),
      xaxis = list(title = "Time"),
      yaxis = list(title = "Residual")
    )
  
  print(resid_plot)
}

# Test residuals for remaining autocorrelation
# Simple Portmanteau test on residuals
resid_Q_stat <- 0
n_resid <- nrow(var_residuals)

for (lag in 1:10) {
  if (lag < n_resid) {
    R_lagged <- var_residuals[1:(n_resid-lag), ]
    R_current <- var_residuals[(lag+1):n_resid, ]
    cor_lag <- cor(R_lagged, R_current)
    resid_Q_stat <- resid_Q_stat + sum(diag(t(cor_lag) %*% cor_lag)) / (n_resid - lag)
  }
}

resid_Q_stat <- n_resid * (n_resid + 2) * resid_Q_stat
resid_df <- p^2 * 10
resid_p_value <- 1 - pchisq(resid_Q_stat, resid_df)

cat("\nResidual Autocorrelation Test:\n")

Residual Autocorrelation Test:

cat("  Q Statistic:", round(resid_Q_stat, 3), "\n")

  Q Statistic: 3079.395 

cat("  p-value:", round(resid_p_value, 4), "\n")

  p-value: 0 

cat("  Interpretation:", ifelse(resid_p_value > 0.05,
                                 "[OK] Residuals are white noise",
                                 "[WARNING] Residuals still autocorrelated"), "\n\n")

  Interpretation: [WARNING] Residuals still autocorrelated 

3. Multivariate Control Charts on VAR Residuals

Hotelling T² Chart

# Calculate Hotelling T² for VAR residuals
# Residuals should now be (approximately) independent and multivariate normal

cat("\nHotelling T² Control Chart on VAR Residuals:\n\n")

Hotelling T² Control Chart on VAR Residuals:

# Estimate parameters from Phase I residuals
mu_resid <- colMeans(var_residuals)
Sigma_resid <- cov(var_residuals)

# Calculate T² for Phase I
m <- nrow(var_residuals)
T2_phase1 <- numeric(m)

for (i in 1:m) {
  r_i <- var_residuals[i, ] - mu_resid
  T2_phase1[i] <- t(r_i) %*% solve(Sigma_resid) %*% r_i
}

# Control limit for Phase I
alpha <- 0.0027
F_crit <- qf(1 - alpha, p, m - p)
UCL_T2 <- (p * (m + 1) * (m - 1)) / (m * (m - p)) * F_crit

cat("Phase I Control Limit (UCL):", round(UCL_T2, 3), "\n")

Phase I Control Limit (UCL): 15.11 

cat("Phase I out-of-control points:", sum(T2_phase1 > UCL_T2), "\n\n")

Phase I out-of-control points: 0 

# Plot T² chart for Phase I
t2_data_phase1 <- tibble(
  Time = 1:m,
  T2 = T2_phase1,
  UCL = UCL_T2,
  Status = ifelse(T2 > UCL, "Out of Control", "In Control")
)

t2_plot_phase1 <- plot_ly(t2_data_phase1) %>%
  add_trace(
    x = ~Time,
    y = ~T2,
    type = "scatter",
    mode = "lines+markers",
    name = "T² Statistic",
    line = list(color = "blue"),
    marker = list(
      size = 6,
      color = ~ifelse(Status == "Out of Control", "red", "blue")
    )
  ) %>%
  add_trace(
    x = ~Time,
    y = ~UCL,
    type = "scatter",
    mode = "lines",
    name = "UCL",
    line = list(color = "red", dash = "dash", width = 3)
  ) %>%
  layout(
    title = list(
      text = "Hotelling T² Chart: Phase I (VAR Residuals)<br><sub>Establishing control limits</sub>",
      font = list(size = 16, weight = "bold")
    ),
    xaxis = list(title = "Time"),
    yaxis = list(title = "T² Statistic")
  )

t2_plot_phase1

Phase II Monitoring with VAR

cat("\nPhase II Monitoring with VAR Model:\n\n")

Phase II Monitoring with VAR Model:

# For Phase II, use VAR model to forecast and calculate residuals
# Then apply T² chart to residuals

# Refit VAR iteratively for Phase II (rolling window or expanding window)
# Here we use expanding window

all_T2 <- numeric(nrow(mts_data))
all_T2[1:m] <- T2_phase1

# For Phase II, use fixed Phase I model (simpler and more stable)
cat("Using Phase I VAR model for Phase II monitoring...\n\n")

Using Phase I VAR model for Phase II monitoring...

for (i in 151:nrow(mts_data)) {
  tryCatch({
    # Get test point as vector
    test_point <- as.numeric(mts_data[i, var_names])
    
    # Use Phase I model to forecast (expanding horizon)
    horizon <- i - 150
    forecast_i <- predict(var_model, n.ahead = horizon)
    predicted <- sapply(forecast_i$fcst, function(x) x[horizon, "fcst"])
    
    # Calculate residual (vector)
    residual_i <- test_point - predicted
    
    # Calculate T² (ensure residual_i is a vector)
    residual_centered <- as.numeric(residual_i - mu_resid)
    all_T2[i] <- t(residual_centered) %*% solve(Sigma_resid) %*% residual_centered
  }, error = function(e) {
    # If prediction fails, use mean as forecast
    test_point <- as.numeric(mts_data[i, var_names])
    predicted <- colMeans(X_phase1)
    residual_i <- test_point - predicted
    residual_centered <- as.numeric(residual_i - mu_resid)
    all_T2[i] <- t(residual_centered) %*% solve(Sigma_resid) %*% residual_centered
  })
}

# Phase II analysis
phase2_T2 <- all_T2[151:200]
ooc_phase2 <- sum(phase2_T2 > UCL_T2)

cat("Phase II out-of-control points:", ooc_phase2, "out of 50\n")

Phase II out-of-control points: 0 out of 50

cat("False alarm rate:", round(ooc_phase2 / 50 * 100, 2), "%\n\n")

False alarm rate: 0 %

# Full monitoring chart
monitoring_data <- tibble(
  Time = 1:length(all_T2),
  T2 = all_T2,
  UCL = UCL_T2,
  Phase = ifelse(Time <= 150, "Phase I", "Phase II"),
  Status = ifelse(T2 > UCL, "Out of Control", "In Control")
)

t2_monitoring_plot <- plot_ly(monitoring_data) %>%
  add_trace(
    x = ~Time,
    y = ~T2,
    type = "scatter",
    mode = "lines+markers",
    name = "T² Statistic",
    line = list(color = "blue"),
    marker = list(
      size = 4,
      color = ~ifelse(Status == "Out of Control", "red", "blue")
    )
  ) %>%
  add_trace(
    x = ~Time,
    y = ~UCL,
    type = "scatter",
    mode = "lines",
    name = "UCL",
    line = list(color = "red", dash = "dash", width = 3)
  ) %>%
  add_segments(
    x = 150.5, xend = 150.5,
    y = 0, yend = max(monitoring_data$T2) * 1.1,
    line = list(color = "gray", dash = "dot", width = 2),
    showlegend = FALSE
  ) %>%
  layout(
    title = list(
      text = "VAR-based T² Monitoring<br><sub>Full process monitoring (Phase I + Phase II)</sub>",
      font = list(size = 16, weight = "bold")
    ),
    xaxis = list(title = "Time"),
    yaxis = list(title = "T² Statistic"),
    annotations = list(
      list(
        x = 150.5,
        y = max(monitoring_data$T2),
        text = "Phase I | Phase II",
        showarrow = FALSE,
        font = list(size = 10, color = "gray")
      )
    )
  )

t2_monitoring_plot

4. Multivariate Time Series Cross-Validation

Rolling Origin Cross-Validation

cat("\n=== Multivariate Time Series Cross-Validation ===\n\n")
cat("Method: Rolling Origin (Expanding Window)\n")
cat("Forecasting all", p, "variables simultaneously\n\n")

# Rolling origin CV for multivariate forecasts
min_train <- 80
h <- 1  # One-step-ahead

# Storage for results
cv_errors <- matrix(0, nrow = 150 - min_train, ncol = p)
cv_T2 <- numeric(150 - min_train)
cv_success_count <- 0

cat("Running cross-validation...\n")

# Completely suppress ALL output
capture.output({
  suppressMessages({
    suppressWarnings({
      for (i in 1:(150 - min_train)) {
        success <- tryCatch({
          train_end <- min_train + i - 1
          test_point <- train_end + 1
          
          if (test_point <= 150 && train_end > optimal_lag + 10) {
            # Training data
            train_data <- X_phase1[1:train_end, ]
            test_data <- X_phase1[test_point, ]
            
            # Fit VAR
            var_cv <- VAR(train_data, p = optimal_lag)
            
            # Check if VAR fitting succeeded
            if (inherits(var_cv, "varest")) {
              # Forecast
              forecast_cv <- predict(var_cv, n.ahead = h)
              predicted <- sapply(forecast_cv$fcst, function(x) x[1, "fcst"])
              
              # Errors
              errors <- test_data - predicted
              cv_errors[i, ] <- errors
              
              # T² for forecast error  
              cv_T2[i] <- t(errors) %*% solve(Sigma_resid) %*% errors
              
              TRUE  # Success
            } else {
              FALSE  # VAR failed
            }
          } else {
            FALSE  # Not enough data
          }
        }, error = function(e) {
          FALSE  # Error occurred
        })
        
        if (success) cv_success_count <- cv_success_count + 1
      }
    })
  })
}, file = nullfile())

cat("✓ Cross-validation complete: ", cv_success_count, " successful iterations out of ", 
    150 - min_train, "\n\n", sep = "")

# Calculate metrics for each variable
# Remove rows where CV failed (all zeros)
valid_idx <- cv_T2 > 0
cv_errors_valid <- cv_errors[valid_idx, , drop = FALSE]

cat("CV Diagnostics:\n")

CV Diagnostics:

cat("  Total iterations:", 150 - min_train, "\n")

  Total iterations: 70 

cat("  Successful iterations:", cv_success_count, "\n")

  Successful iterations: 0 

cat("  Valid T² values:", sum(valid_idx), "\n")

  Valid T² values: 0 

cat("  Valid error rows:", nrow(cv_errors_valid), "\n")

  Valid error rows: 0 

if (sum(valid_idx) > 5) {  # Need at least 5 valid iterations
  cat("  Computing metrics from valid data...\n\n")
  
  cv_metrics <- tibble(
    Variable = var_names,
    MAE = colMeans(abs(cv_errors_valid)),
    RMSE = sqrt(colMeans(cv_errors_valid^2)),
    ME = colMeans(cv_errors_valid)
  )
  
  cat("Overall CV Performance:\n")
  cat("  Mean MAE across variables:", round(mean(cv_metrics$MAE, na.rm = TRUE), 4), "\n")
  cat("  Mean RMSE across variables:", round(mean(cv_metrics$RMSE, na.rm = TRUE), 4), "\n\n")
} else {
  cat("  [WARNING] Insufficient valid iterations (need >5, have", sum(valid_idx), ")\n\n")
  
  cv_metrics <- tibble(
    Variable = var_names,
    MAE = rep(NA, p),
    RMSE = rep(NA, p),
    ME = rep(NA, p)
  )
}

  [WARNING] Insufficient valid iterations (need >5, have 0 )

cv_metrics %>%
  gt() %>%
  tab_header(
    title = "Multivariate Cross-Validation Results",
    subtitle = "One-step-ahead forecast performance by variable"
  ) %>%
  fmt_number(
    columns = c(MAE, RMSE, ME),
    decimals = 4
  ) %>%
  sub_missing(
    columns = everything(),
    missing_text = "—"
  ) %>%
  tab_style(
    style = cell_fill(color = "lightblue"),
    locations = cells_column_labels()
  ) %>%
  tab_footnote(
    footnote = paste("Based on", cv_success_count, "successful CV iterations"),
    locations = cells_column_labels(columns = MAE)
  )

Multivariate Cross-Validation Results
One-step-ahead forecast performance by variable
Variable	MAE1	RMSE	ME
temperature	—	—	—
pressure	—	—	—
flow_rate	—	—	—
1 Based on 0 successful CV iterations

# Plot CV errors for first variable (only valid data)
if (sum(valid_idx) > 0) {
  cv_plot_data <- tibble(
    Fold = which(valid_idx),
    Error = cv_errors_valid[, 1],
    `Training Size` = min_train + which(valid_idx) - 1
  )

  cv_error_plot <- plot_ly(cv_plot_data) %>%
    add_trace(
      x = ~`Training Size`,
      y = ~Error,
      type = "scatter",
    mode = "lines+markers",
    name = var_names[1],
    line = list(color = "blue"),
    marker = list(size = 4)
  ) %>%
  add_segments(
    x = min(cv_plot_data$`Training Size`),
    xend = max(cv_plot_data$`Training Size`),
    y = 0, yend = 0,
    line = list(color = "black", width = 1),
    showlegend = FALSE
  ) %>%
  layout(
    title = list(
      text = paste0("CV Forecast Errors: ", var_names[1], 
                   "<br><sub>Rolling origin cross-validation (", cv_success_count, " valid iterations)</sub>"),
      font = list(size = 16, weight = "bold")
    ),
    xaxis = list(title = "Training Sample Size"),
    yaxis = list(title = "Forecast Error")
  )

  cv_error_plot

  # Plot T² for CV (multivariate measure) - only valid data
  cv_t2_valid <- cv_T2[valid_idx]
  cv_t2_data <- tibble(
    Fold = which(valid_idx),
    T2 = cv_t2_valid,
    `Training Size` = min_train + which(valid_idx) - 1
  )

  cv_t2_plot <- plot_ly(cv_t2_data) %>%
    add_trace(
      x = ~`Training Size`,
      y = ~T2,
      type = "scatter",
      mode = "lines+markers",
      name = "T² (CV)",
      line = list(color = "purple"),
      marker = list(size = 4)
    ) %>%
    layout(
      title = list(
        text = "Multivariate Forecast Error (T²)<br><sub>Combined measure across all variables</sub>",
        font = list(size = 16, weight = "bold")
      ),
      xaxis = list(title = "Training Sample Size"),
      yaxis = list(title = "T² Statistic")
    )

  cv_t2_plot
} else {
  cat("\n[INFO] Insufficient valid CV data for plotting\n")
}

[INFO] Insufficient valid CV data for plotting

Multivariate Block Bootstrap

cat("\nMultivariate Block Bootstrap:\n\n")

Multivariate Block Bootstrap:

cat("Method: Moving Block Bootstrap\n")

Method: Moving Block Bootstrap

cat("Preserves both autocorrelation AND cross-correlation\n\n")

Preserves both autocorrelation AND cross-correlation

# Block bootstrap for multivariate time series
block_length <- 10
n_bootstrap <- 500  # Reduced for speed
n_obs <- nrow(X_phase1)

# Storage
boot_means <- matrix(0, nrow = n_bootstrap, ncol = p)
boot_cors <- numeric(n_bootstrap)

set.seed(999)

for (b in 1:n_bootstrap) {
  # Generate bootstrap sample
  n_blocks <- ceiling(n_obs / block_length)
  boot_sample <- matrix(0, nrow = 0, ncol = p)
  
  for (i in 1:n_blocks) {
    start <- sample(1:(n_obs - block_length + 1), 1)
    block <- X_phase1[start:(start + block_length - 1), ]
    boot_sample <- rbind(boot_sample, block)
  }
  
  # Trim
  boot_sample <- boot_sample[1:n_obs, ]
  
  # Calculate statistics
  boot_means[b, ] <- colMeans(boot_sample)
  
  # Correlation between first two variables
  if (p >= 2) {
    boot_cors[b] <- cor(boot_sample[, 1], boot_sample[, 2])
  }
}

# Confidence intervals
alpha_boot <- 0.05
ci_means <- apply(boot_means, 2, quantile, probs = c(alpha_boot/2, 1-alpha_boot/2))

if (p >= 2) {
  ci_cor <- quantile(boot_cors, probs = c(alpha_boot/2, 1-alpha_boot/2))
}

cat("Bootstrap Results (", n_bootstrap, "iterations, block =", block_length, "):\n\n")

Bootstrap Results ( 500 iterations, block = 10 ):

# Results table
boot_results <- tibble(
  Variable = var_names,
  `Original Mean` = colMeans(X_phase1),
  `Bootstrap Mean` = colMeans(boot_means),
  `CI Lower` = ci_means[1, ],
  `CI Upper` = ci_means[2, ]
)

boot_results %>%
  gt() %>%
  tab_header(
    title = "Multivariate Block Bootstrap Results",
    subtitle = "95% confidence intervals for means"
  ) %>%
  fmt_number(
    columns = c(`Original Mean`, `Bootstrap Mean`, `CI Lower`, `CI Upper`),
    decimals = 3
  )

Multivariate Block Bootstrap Results
95% confidence intervals for means
Variable	Original Mean	Bootstrap Mean	CI Lower	CI Upper
temperature	98.072	97.985	96.538	99.560
pressure	198.243	198.180	196.659	199.660
flow_rate	48.368	48.313	46.523	50.207

if (p >= 2) {
  cat("\nCorrelation Bootstrap:\n")
  cat("Original:", round(cor(X_phase1[, 1], X_phase1[, 2]), 4), "\n")
  cat("Bootstrap mean:", round(mean(boot_cors), 4), "\n")
  cat("95% CI: [", round(ci_cor[1], 4), ",", round(ci_cor[2], 4), "]\n\n")
  
  # Histogram of bootstrap correlations
  boot_cor_data <- tibble(Correlation = boot_cors)
  
  cor_hist <- plot_ly(x = boot_cor_data$Correlation, type = "histogram",
                      marker = list(color = "lightblue",
                                   line = list(color = "darkblue", width = 1)),
                      nbinsx = 30) %>%
    add_segments(
      x = cor(X_phase1[, 1], X_phase1[, 2]),
      xend = cor(X_phase1[, 1], X_phase1[, 2]),
      y = 0, yend = 50,
      line = list(color = "red", width = 3),
      name = "Original"
    ) %>%
    add_segments(
      x = ci_cor[1], xend = ci_cor[1],
      y = 0, yend = 50,
      line = list(color = "green", dash = "dash", width = 2),
      name = "95% CI"
    ) %>%
    add_segments(
      x = ci_cor[2], xend = ci_cor[2],
      y = 0, yend = 50,
      line = list(color = "green", dash = "dash", width = 2),
      showlegend = FALSE
    ) %>%
    layout(
      title = list(
        text = paste0("Bootstrap Distribution of Correlation<br><sub>",
                     var_names[1], " vs ", var_names[2], "</sub>"),
        font = list(size = 16, weight = "bold")
      ),
      xaxis = list(title = "Correlation"),
      yaxis = list(title = "Frequency")
    )
  
  cor_hist
}

Correlation Bootstrap:
Original: 0.8437 
Bootstrap mean: 0.828 
95% CI: [ 0.6781 , 0.9138 ]

5. Dashboard and Summary

# Summary dashboard
dashboard_metrics <- tibble(
  Metric = c(
    "Total Observations",
    "Variables",
    "Phase I Samples",
    "Phase II Samples",
    "VAR Model Order",
    "Strongest Correlation",
    "Phase I T² OOC",
    "Phase II T² OOC",
    "CV Success Rate",
    "Mean CV MAE"
  ),
  Value = c(
    nrow(mts_data),
    p,
    150,
    50,
    optimal_lag,
    paste0(round(max(abs(cor_phase1[lower.tri(cor_phase1)])), 3)),
    sum(T2_phase1 > UCL_T2),
    ooc_phase2,
    paste0(cv_success_count, "/", 150 - min_train),
    ifelse(cv_success_count > 0 && !all(is.na(cv_metrics$MAE)),
           round(mean(cv_metrics$MAE, na.rm = TRUE), 4),
           "—")
  )
)

dashboard_metrics %>%
  gt() %>%
  tab_header(
    title = "Multivariate Time Series SPC Dashboard",
    subtitle = paste("Analysis Date:", Sys.Date())
  ) %>%
  tab_style(
    style = cell_fill(color = "lightblue"),
    locations = cells_column_labels()
  )

Multivariate Time Series SPC Dashboard
Analysis Date: 2025-12-22
Metric	Value
Total Observations	200
Variables	3
Phase I Samples	150
Phase II Samples	50
VAR Model Order	1
Strongest Correlation	0.851
Phase I T² OOC	0
Phase II T² OOC	0
CV Success Rate	0/70
Mean CV MAE	—

6. Recommendations

When to Use Multivariate Time Series SPC

Use VAR-based SPC when: - ✅ Multiple variables (2+) - ✅ Variables correlated with each other - ✅ Each variable autocorrelated - ✅ Lead-lag relationships between variables - ✅ Want to monitor system holistically

Key Advantages: 1. Models cross-correlations explicitly 2. Captures lead-lag dynamics 3. Single T² chart for all variables 4. Better detection of multivariate shifts 5. More efficient than separate ARIMA models

Implementation Guidelines

cat("\n### Analysis-Specific Recommendations:\n\n")

### Analysis-Specific Recommendations:

# Check if multivariate approach needed
max_cor <- max(abs(cor_phase1[lower.tri(cor_phase1)]))

if (max_cor > 0.5) {
  cat("[REQUIRED] Strong cross-correlations detected (max =", round(max_cor, 3), ")\n")
  cat("  Action: Use multivariate approach (VAR + T²)\n")
  cat("  Benefit: Captures variable interactions\n\n")
} else {
  cat("[OPTIONAL] Weak cross-correlations\n")
  cat("  Action: Could use separate ARIMA models\n")
  cat("  Alternative: VAR still provides system view\n\n")
}

[REQUIRED] Strong cross-correlations detected (max = 0.851 )
  Action: Use multivariate approach (VAR + T²)
  Benefit: Captures variable interactions

if (resid_p_value > 0.05) {
  cat("[OK] VAR residuals are white noise (p =", round(resid_p_value, 4), ")\n")
  cat("  Interpretation: Model adequately captures dynamics\n")
  cat("  Action: T² chart on residuals is valid\n\n")
} else {
  cat("[WARNING] VAR residuals still show autocorrelation\n")
  cat("  Action: Consider higher lag order or different model\n\n")
}

[WARNING] VAR residuals still show autocorrelation
  Action: Consider higher lag order or different model

if (ooc_phase2 > 5) {
  cat("[ALERT] Multiple Phase II signals (", ooc_phase2, "/50)\n")
  cat("  Action: Investigate process change\n")
  cat("  First signal at sample ~", min(which(all_T2[151:200] > UCL_T2)) + 150, "\n\n")
}

cat("Model Performance:\n")

Model Performance:

cat("  VAR order:", optimal_lag, "\n")

  VAR order: 1 

cat("  CV MAE:", round(mean(cv_metrics$MAE), 4), "\n")

  CV MAE: NA 

cat("  Bootstrap CI stable:", ifelse(all(ci_means[2,] - ci_means[1,] < 5), "Yes", "No"), "\n")

  Bootstrap CI stable: Yes 

---

Report Generated: 2025-12-22 17:03:56.241547
Analysis Method: VAR-based Multivariate Time Series SPC
Statistical Software: R with vars, MTS, and forecast packages
Methodology: Vector autoregression with Hotelling T² monitoring