Statistical Tests and Diagnostic Checks

# ============================================================
# Required package: R2DAG
# ============================================================
# The replication code uses functions from the author's R package
# R2DAG. The package is available from GitHub and can be installed
# with the command below.
#
# dependencies = FALSE and upgrade = "never" are used to avoid
# automatically modifying the user's existing R package environment.
# After installation, the package is loaded for the LiNGAM-based
# residual extraction used later in the notebook.

devtools::install_github(
  "espanm/R2DAG",
  dependencies = FALSE,
  upgrade = "never"
)

library(R2DAG)

library(readr)
library(xts)

# Read the one-minute log-return dataset created by the data-preparation notebook.
# Replace [YOUR FOLDER] with the folder where full_dataset.csv is stored.
df <- read_csv("[YOUR FOLDER]/full_dataset.csv")

# Event windows are defined in Eastern/New York trading time.
df$datetime <- as.POSIXct(df$datetime, tz = "America/New_York")

# These non-AI tickers were included in the broad data file but are not used in
# this notebook's event-day statistical tables.
non_ai_tickers <- c("WMT", "COST", "CVX", "ABBV", "PG", "MRK", "PM", "MCD", "PEP", "T")

# Drop only those non-AI tickers that are actually present in the data.
cols_to_drop <- intersect(non_ai_tickers, colnames(df))

df <- df[, !colnames(df) %in% cols_to_drop]

# Convert the data frame into an xts object for convenient intraday subsetting.
data_merged <- xts(df[, -1], order.by = df$datetime)

library(xts)
library(moments)
library(readr)

# The SP500 column is required because demarketing is done by regressing
# each stock return on the market return.
stopifnot("SP500" %in% colnames(data_merged))

# Event dates used in the empirical analysis.
events <- c(
  "ChatGPT launch" = "2022-11-30",
  "NVIDIA Hopper/H100" = "2022-03-22",
  "GitHub Copilot preview" = "2021-06-29",
  "Claude Opus 4.6" = "2026-02-05"
)

# Extract the regular U.S. trading session for a given event day.
get_event_window <- function(data, date_str,
                             start_hms = "09:30:00",
                             end_hms = "16:00:00") {
  from <- paste0(date_str, " ", start_hms)
  to <- paste0(date_str, " ", end_hms)
  data[paste0(from, "/", to)]
}

# Estimate the degrees-of-freedom parameter of a Student-t density by maximum
# likelihood. The transformation nu = 2 + exp(theta) keeps the fitted degrees
# of freedom above 2, so the variance of the fitted Student-t distribution exists.
fit_student_t_df <- function(z) {
  z <- na.omit(as.numeric(z))

  if (length(z) < 30) return(NA_real_)

  sdz <- sd(z)
  if (!is.finite(sdz) || sdz <= 0) return(NA_real_)

  negloglik_t <- function(par, x) {
    mu <- par[1]
    sigma <- exp(par[2])
    nu <- 2 + exp(par[3])

    ll <- dt((x - mu) / sigma, df = nu, log = TRUE) - log(sigma)

    if (any(!is.finite(ll))) return(1e12)
    -sum(ll)
  }

  # Start the optimization at the sample mean, sample standard deviation,
  # and an initial Student-t degrees-of-freedom value of 8.
  init <- c(
    mean(z),
    log(sdz + 1e-8),
    log(8 - 2)
  )

  fit <- tryCatch(
    optim(
      par = init,
      fn = negloglik_t,
      x = z,
      method = "BFGS",
      hessian = FALSE,
      control = list(maxit = 2000, reltol = 1e-10)
    ),
    error = function(e) NULL
  )

  if (is.null(fit)) return(NA_real_)
  if (!is.list(fit) || is.null(fit$convergence) || fit$convergence != 0) return(NA_real_)

  nu_hat <- 2 + exp(fit$par[3])

  if (!is.finite(nu_hat)) return(NA_real_)
  nu_hat
}

# Compute all descriptive statistics for one stock on one event day.
compute_one_asset_one_event <- function(day_data, ticker) {
  stat_names <- c(
    "Mean",
    "Variance",
    "Skewness",
    "Excess Kurtosis",
    "Student-t density shape parameter"
  )

  # Return NA values if the requested ticker is not available.
  if (!(ticker %in% colnames(day_data))) {
    out <- rep(NA_real_, length(stat_names))
    names(out) <- stat_names
    return(out)
  }

  # Keep the stock and market return only, then drop missing observations.
  sub <- day_data[, c(ticker, "SP500")]
  sub <- na.omit(sub)

  if (NROW(sub) < 30) {
    out <- rep(NA_real_, length(stat_names))
    names(out) <- stat_names
    return(out)
  }

  y <- as.numeric(sub[, ticker])
  x <- as.numeric(sub[, "SP500"])

  # Demarket the stock return by taking residuals from y_t = alpha + beta * SP500_t + eps_t.
  fit_lm <- lm(y ~ x)
  eps <- resid(fit_lm)

  if (length(eps) < 30 || sd(eps) <= 0 || !is.finite(sd(eps))) {
    out <- rep(NA_real_, length(stat_names))
    names(out) <- stat_names
    return(out)
  }

  c(
    Mean = mean(eps),
    Variance = var(eps),
    Skewness = moments::skewness(eps),
    `Excess Kurtosis` = moments::kurtosis(eps) - 3,
    `Student-t density shape parameter` = fit_student_t_df(eps)
  )
}

# All stock columns except the market index are treated as assets.
asset_names <- setdiff(colnames(data_merged), "SP500")

event_asset_results <- list()
panel_results <- list()

for (ev_name in names(events)) {
  ev_date <- events[[ev_name]]
  day_data <- get_event_window(data_merged, ev_date)

  # Compute a stock-by-statistic table for the current event.
  # If the event window is unavailable, keep the table structure and fill it with NA.
  if (NROW(day_data) == 0) {
    stats_mat <- matrix(NA_real_, nrow = 5, ncol = length(asset_names))
    rownames(stats_mat) <- c(
      "Mean",
      "Variance",
      "Skewness",
      "Excess Kurtosis",
      "Student-t density shape parameter"
    )
    colnames(stats_mat) <- asset_names
  } else {
    stats_mat <- sapply(asset_names, function(tkr) {
      compute_one_asset_one_event(day_data, tkr)
    })
  }

  # Store the asset-level results in case the detailed table is needed later.
  event_asset_results[[ev_name]] <- as.data.frame(stats_mat)

  # Collapse the cross-section into median, minimum, and maximum values,
  # matching the summary format used in the paper.
  panel_mat <- t(apply(stats_mat, 1, function(z) {
    z <- as.numeric(z)
    c(
      median = if (all(is.na(z))) NA_real_ else median(z, na.rm = TRUE),
      min = if (all(is.na(z))) NA_real_ else min(z, na.rm = TRUE),
      max = if (all(is.na(z))) NA_real_ else max(z, na.rm = TRUE)
    )
  }))

  panel_results[[ev_name]] <- panel_mat
}

# Combine the four event-level summary tables into one publication-style table.
final_table <- cbind(
  panel_results[[1]],
  panel_results[[2]],
  panel_results[[3]],
  panel_results[[4]]
)

final_table <- as.data.frame(final_table)

colnames(final_table) <- c(
  "Event1_median", "Event1_min", "Event1_max",
  "Event2_median", "Event2_min", "Event2_max",
  "Event3_median", "Event3_min", "Event3_max",
  "Event4_median", "Event4_min", "Event4_max"
)

rownames(final_table) <- c(
  "Mean",
  "Variance",
  "Skewness",
  "Excess Kurtosis",
  "Student-t density shape parameter"
)

# Rename columns with readable event labels.
final_table_named <- final_table

colnames(final_table_named) <- c(
  "ChatGPT_median", "ChatGPT_min", "ChatGPT_max",
  "Hopper_median", "Hopper_min", "Hopper_max",
  "Copilot_median", "Copilot_min", "Copilot_max",
  "Claude_median", "Claude_min", "Claude_max"
)

# Format small mean and variance values in scientific notation, while keeping
# the other statistics in a compact four-decimal format.
format_final_table <- function(df) {
  out <- as.data.frame(matrix("", nrow = nrow(df), ncol = ncol(df)))
  rownames(out) <- rownames(df)
  colnames(out) <- colnames(df)

  sci_rows <- c("Mean", "Variance")

  for (i in seq_len(nrow(df))) {
    rn <- rownames(df)[i]

    if (rn %in% sci_rows) {
      out[i, ] <- sapply(df[i, ], function(x) {
        if (is.na(x)) {
          NA_character_
        } else {
          formatC(x, format = "e", digits = 4)
        }
      })
    } else {
      out[i, ] <- sapply(df[i, ], function(x) {
        if (is.na(x)) {
          NA_character_
        } else {
          sprintf("%.4f", round(x, 4))
        }
      })
    }
  }

  out
}

final_table_named_fmt <- format_final_table(final_table_named)

print(final_table_named_fmt, quote = FALSE)

library(xts)
library(tseries)

# The market index is required to construct demarketed residuals.
stopifnot("SP500" %in% colnames(data_merged))

events <- c(
  "ChatGPT launch" = "2022-11-30",
  "NVIDIA Hopper/H100" = "2022-03-22",
  "GitHub Copilot preview" = "2021-06-29",
  "Claude Opus 4.6" = "2026-02-05"
)

# Extract the regular trading session for the selected event date.
get_event_window <- function(data, date_str,
                             start_hms = "09:30:00",
                             end_hms = "16:00:00") {
  from <- paste0(date_str, " ", start_hms)
  to <- paste0(date_str, " ", end_hms)
  data[paste0(from, "/", to)]
}

# Compute normality-test results for one stock on one event day.
compute_test_results <- function(day_data, ticker) {
  out_na <- c(
    `Jarque and Bera p-value` = NA_real_,
    `Jarque and Bera statistic` = NA_real_,
    `Shapiro - Wilk p-value` = NA_real_,
    `Shapiro - Wilk statistic` = NA_real_
  )

  if (!(ticker %in% colnames(day_data))) {
    return(out_na)
  }

  # Demarketing is done on the same observations used in the event-day test.
  sub <- day_data[, c(ticker, "SP500")]
  sub <- na.omit(sub)

  if (NROW(sub) < 8) {
    return(out_na)
  }

  y <- as.numeric(sub[, ticker])
  x <- as.numeric(sub[, "SP500"])

  fit_lm <- lm(y ~ x)
  eps <- resid(fit_lm)
  eps <- eps[is.finite(eps)]

  if (length(eps) < 8 || sd(eps) <= 0 || !is.finite(sd(eps))) {
    return(out_na)
  }

  # Jarque--Bera is based on skewness and kurtosis.
  jb_out <- tryCatch(
    tseries::jarque.bera.test(eps),
    error = function(e) NULL
  )

  # Shapiro--Wilk is an additional finite-sample normality diagnostic.
  sw_out <- tryCatch(
    shapiro.test(eps),
    error = function(e) NULL
  )

  c(
    `Jarque and Bera p-value` = if (is.null(jb_out)) NA_real_ else unname(jb_out$p.value),
    `Jarque and Bera statistic` = if (is.null(jb_out)) NA_real_ else unname(jb_out$statistic),
    `Shapiro - Wilk p-value` = if (is.null(sw_out)) NA_real_ else unname(sw_out$p.value),
    `Shapiro - Wilk statistic` = if (is.null(sw_out)) NA_real_ else unname(sw_out$statistic)
  )
}

asset_names <- setdiff(colnames(data_merged), "SP500")

event_test_results <- list()
panel_results <- list()

for (ev_name in names(events)) {
  ev_date <- events[[ev_name]]
  day_data <- get_event_window(data_merged, ev_date)

  # Compute stock-level test results. Missing event windows are represented by NA tables.
  if (NROW(day_data) == 0) {
    res_mat <- matrix(NA_real_, nrow = 4, ncol = length(asset_names))
    rownames(res_mat) <- c(
      "Jarque and Bera p-value",
      "Jarque and Bera statistic",
      "Shapiro - Wilk p-value",
      "Shapiro - Wilk statistic"
    )
    colnames(res_mat) <- asset_names
  } else {
    res_mat <- sapply(asset_names, function(tkr) {
      compute_test_results(day_data, tkr)
    })
  }

  event_test_results[[ev_name]] <- as.data.frame(res_mat)

  # Summarize the cross-section of tickers for the table in the paper.
  panel_mat <- t(apply(res_mat, 1, function(z) {
    z <- as.numeric(z)
    c(
      median = if (all(is.na(z))) NA_real_ else median(z, na.rm = TRUE),
      min = if (all(is.na(z))) NA_real_ else min(z, na.rm = TRUE),
      max = if (all(is.na(z))) NA_real_ else max(z, na.rm = TRUE)
    )
  }))

  panel_results[[ev_name]] <- panel_mat
}

# Bind all event summaries into a single table.
final_table_named <- cbind(
  panel_results[["ChatGPT launch"]],
  panel_results[["NVIDIA Hopper/H100"]],
  panel_results[["GitHub Copilot preview"]],
  panel_results[["Claude Opus 4.6"]]
)

final_table_named <- as.data.frame(final_table_named)

colnames(final_table_named) <- c(
  "ChatGPT_median", "ChatGPT_min", "ChatGPT_max",
  "Hopper_median", "Hopper_min", "Hopper_max",
  "Copilot_median", "Copilot_min", "Copilot_max",
  "Claude_median", "Claude_min", "Claude_max"
)

rownames(final_table_named) <- c(
  "Jarque and Bera p-value",
  "Jarque and Bera statistic",
  "Shapiro - Wilk p-value",
  "Shapiro - Wilk statistic"
)

print(round(final_table_named, 6))

library(xts)
library(moments)
library(readr)

# The SP500 column is needed for the first-stage demarketing regression.
stopifnot("SP500" %in% colnames(data_merged))

events <- c(
  "ChatGPT launch" = "2022-11-30",
  "NVIDIA Hopper/H100" = "2022-03-22",
  "GitHub Copilot preview" = "2021-06-29",
  "Claude Opus 4.6" = "2026-02-05"
)

# Extract the regular U.S. trading session for each event day.
get_event_window <- function(data, date_str,
                             start_hms = "09:30:00",
                             end_hms   = "16:00:00") {
  from <- paste0(date_str, " ", start_hms)
  to   <- paste0(date_str, " ", end_hms)
  data[paste0(from, "/", to)]
}

# Fit a Student-t density to one residual series and return the estimated
# degrees-of-freedom parameter. Smaller values indicate heavier tails.
fit_student_t_df <- function(z) {
  z <- na.omit(as.numeric(z))

  if (length(z) < 30) return(NA_real_)

  sdz <- sd(z)
  if (!is.finite(sdz) || sdz <= 0) return(NA_real_)

  negloglik_t <- function(par, x) {
    mu    <- par[1]
    sigma <- exp(par[2])
    nu    <- 2 + exp(par[3])

    ll <- dt((x - mu) / sigma, df = nu, log = TRUE) - log(sigma)

    if (any(!is.finite(ll))) return(1e12)
    -sum(ll)
  }

  init <- c(
    mean(z),
    log(sdz + 1e-8),
    log(8 - 2)
  )

  fit <- tryCatch(
    optim(
      par = init,
      fn = negloglik_t,
      x = z,
      method = "BFGS",
      hessian = FALSE,
      control = list(maxit = 2000, reltol = 1e-10)
    ),
    error = function(e) NULL
  )

  if (is.null(fit)) return(NA_real_)
  if (!is.list(fit) || is.null(fit$convergence) || fit$convergence != 0) return(NA_real_)

  nu_hat <- 2 + exp(fit$par[3])

  if (!is.finite(nu_hat)) return(NA_real_)
  nu_hat
}

# Compute descriptive statistics for a single residual vector.
compute_stats_from_vector <- function(z) {
  stat_names <- c(
    "Mean",
    "Variance",
    "Skewness",
    "Excess Kurtosis",
    "Student-t density shape parameter"
  )

  z <- na.omit(as.numeric(z))

  if (length(z) < 30 || sd(z) <= 0 || !is.finite(sd(z))) {
    out <- rep(NA_real_, length(stat_names))
    names(out) <- stat_names
    return(out)
  }

  c(
    Mean = mean(z),
    Variance = var(z),
    Skewness = moments::skewness(z),
    `Excess Kurtosis` = moments::kurtosis(z) - 3,
    `Student-t density shape parameter` = fit_student_t_df(z)
  )
}

# Construct the demarketed return matrix for one event window.
# Each stock is regressed separately on SP500, and the residuals are collected.
get_demarketed_matrix <- function(day_data, asset_names) {
  needed <- c(asset_names, "SP500")
  sub <- day_data[, needed]
  sub <- na.omit(sub)

  if (NROW(sub) < 30) return(NULL)

  x_mkt <- as.numeric(sub[, "SP500"])

  eps_mat <- sapply(asset_names, function(tkr) {
    y <- as.numeric(sub[, tkr])
    fit <- tryCatch(lm(y ~ x_mkt), error = function(e) NULL)

    if (is.null(fit)) {
      rep(NA_real_, length(y))
    } else {
      resid(fit)
    }
  })

  eps_mat <- as.matrix(eps_mat)
  colnames(eps_mat) <- asset_names
  eps_mat <- eps_mat[complete.cases(eps_mat), , drop = FALSE]

  if (nrow(eps_mat) < 30) return(NULL)
  eps_mat
}

# Run the LiNGAM-based R2_network routine and extract the residuals.
# These residuals are the estimated structural shocks used in the appendix table.
get_lingam_residuals_from_R2_network <- function(X) {
  X <- as.matrix(X)

  if (nrow(X) < 30 || ncol(X) < 2) return(NULL)

  X_df <- as.data.frame(X)

  fit <- tryCatch(
    R2_network(
      data = X_df,
      directed = TRUE,
      amat = FALSE,
      dag_method = "lingam",
      standardize_for_dag = FALSE
    ),
    error = function(e) NULL
  )

  if (is.null(fit)) return(NULL)
  if (is.null(fit$lingam_residuals)) return(NULL)

  E <- as.matrix(fit$lingam_residuals)

  if (!all(dim(E) == dim(X))) return(NULL)

  colnames(E) <- colnames(X)
  E
}

asset_names <- setdiff(colnames(data_merged), "SP500")

event_asset_results <- list()
panel_results <- list()

for (ev_name in names(events)) {

  ev_date <- events[[ev_name]]
  day_data <- get_event_window(data_merged, ev_date)

  stat_names <- c(
    "Mean",
    "Variance",
    "Skewness",
    "Excess Kurtosis",
    "Student-t density shape parameter"
  )

  # Keep a fixed table shape even if a window, demarketing step, or LiNGAM fit fails.
  if (NROW(day_data) == 0) {
    stats_mat <- matrix(NA_real_, nrow = length(stat_names), ncol = length(asset_names))
    rownames(stats_mat) <- stat_names
    colnames(stats_mat) <- asset_names

  } else {
    demarketed_X <- get_demarketed_matrix(day_data, asset_names)

    if (is.null(demarketed_X)) {
      stats_mat <- matrix(NA_real_, nrow = length(stat_names), ncol = length(asset_names))
      rownames(stats_mat) <- stat_names
      colnames(stats_mat) <- asset_names

    } else {
      lingam_E <- get_lingam_residuals_from_R2_network(demarketed_X)

      if (is.null(lingam_E)) {
        stats_mat <- matrix(NA_real_, nrow = length(stat_names), ncol = length(asset_names))
        rownames(stats_mat) <- stat_names
        colnames(stats_mat) <- asset_names

      } else {
        stats_mat <- sapply(asset_names, function(tkr) {
          compute_stats_from_vector(lingam_E[, tkr])
        })
      }
    }
  }

  event_asset_results[[ev_name]] <- as.data.frame(stats_mat)

  # Summarize residual statistics across assets for each event.
  panel_mat <- t(apply(stats_mat, 1, function(z) {
    z <- as.numeric(z)
    c(
      median = if (all(is.na(z))) NA_real_ else median(z, na.rm = TRUE),
      min    = if (all(is.na(z))) NA_real_ else min(z, na.rm = TRUE),
      max    = if (all(is.na(z))) NA_real_ else max(z, na.rm = TRUE)
    )
  }))

  panel_results[[ev_name]] <- panel_mat
}

final_table <- cbind(
  panel_results[[1]],
  panel_results[[2]],
  panel_results[[3]],
  panel_results[[4]]
)

final_table <- as.data.frame(final_table)

colnames(final_table) <- c(
  "Event1_median", "Event1_min", "Event1_max",
  "Event2_median", "Event2_min", "Event2_max",
  "Event3_median", "Event3_min", "Event3_max",
  "Event4_median", "Event4_min", "Event4_max"
)

rownames(final_table) <- c(
  "Mean",
  "Variance",
  "Skewness",
  "Excess Kurtosis",
  "Student-t density shape parameter"
)

final_table_named <- final_table
colnames(final_table_named) <- c(
  "ChatGPT_median", "ChatGPT_min", "ChatGPT_max",
  "Hopper_median", "Hopper_min", "Hopper_max",
  "Copilot_median", "Copilot_min", "Copilot_max",
  "Claude_median", "Claude_min", "Claude_max"
)

# Use the same formatting convention as in the demarketed-return table.
format_final_table <- function(df) {
  out <- as.data.frame(matrix("", nrow = nrow(df), ncol = ncol(df)))
  rownames(out) <- rownames(df)
  colnames(out) <- colnames(df)

  sci_rows <- c("Mean", "Variance")

  for (i in seq_len(nrow(df))) {
    rn <- rownames(df)[i]

    if (rn %in% sci_rows) {
      out[i, ] <- sapply(df[i, ], function(x) {
        if (is.na(x)) {
          NA_character_
        } else {
          formatC(x, format = "e", digits = 4)
        }
      })
    } else {
      out[i, ] <- sapply(df[i, ], function(x) {
        if (is.na(x)) {
          NA_character_
        } else {
          sprintf("%.4f", round(x, 4))
        }
      })
    }
  }

  out
}

final_table_named_fmt <- format_final_table(final_table_named)

print(final_table_named_fmt, quote = FALSE)

library(xts)
library(readr)
library(tseries)

# SP500 is required for the first-stage demarketing step.
stopifnot("SP500" %in% colnames(data_merged))

events <- c(
  "ChatGPT launch" = "2022-11-30",
  "NVIDIA Hopper/H100" = "2022-03-22",
  "GitHub Copilot preview" = "2021-06-29",
  "Claude Opus 4.6" = "2026-02-05"
)

# Extract the regular trading session for one event day.
get_event_window <- function(data, date_str,
                             start_hms = "09:30:00",
                             end_hms   = "16:00:00") {
  from <- paste0(date_str, " ", start_hms)
  to   <- paste0(date_str, " ", end_hms)
  data[paste0(from, "/", to)]
}

# First-stage demarketing: regress each stock on SP500 and collect residuals.
get_demarketed_matrix <- function(day_data, asset_names) {
  needed <- c(asset_names, "SP500")
  sub <- day_data[, needed]
  sub <- na.omit(sub)

  if (NROW(sub) < 8) return(NULL)

  x_mkt <- as.numeric(sub[, "SP500"])

  eps_mat <- sapply(asset_names, function(tkr) {
    y <- as.numeric(sub[, tkr])

    fit <- tryCatch(
      lm(y ~ x_mkt),
      error = function(e) NULL
    )

    if (is.null(fit)) {
      rep(NA_real_, length(y))
    } else {
      resid(fit)
    }
  })

  eps_mat <- as.matrix(eps_mat)
  colnames(eps_mat) <- asset_names
  eps_mat <- eps_mat[complete.cases(eps_mat), , drop = FALSE]

  if (nrow(eps_mat) < 8) return(NULL)
  eps_mat
}

# Run the LiNGAM-based network routine and extract its residual matrix.
get_lingam_residuals_from_R2_network <- function(X) {
  X <- as.matrix(X)

  if (nrow(X) < 8 || ncol(X) < 2) return(NULL)

  X_df <- as.data.frame(X)

  fit <- tryCatch(
    R2_network(
      data = X_df,
      directed = TRUE,
      amat = FALSE,
      dag_method = "lingam",
      standardize_for_dag = FALSE
    ),
    error = function(e) NULL
  )

  if (is.null(fit)) return(NULL)
  if (is.null(fit$lingam_residuals)) return(NULL)

  E <- as.matrix(fit$lingam_residuals)

  if (!all(dim(E) == dim(X))) return(NULL)

  colnames(E) <- colnames(X)
  E
}

# Apply the two normality tests to one residual vector.
compute_test_results_from_vector <- function(z) {
  z <- as.numeric(z)
  z <- z[is.finite(z)]

  out <- c(
    `Jarque-Bera statistic` = NA_real_,
    `Jarque-Bera p-value`   = NA_real_,
    `Shapiro-Wilk W`        = NA_real_,
    `Shapiro-Wilk p-value`  = NA_real_
  )

  if (length(z) < 8 || sd(z) <= 0 || !is.finite(sd(z))) {
    return(out)
  }

  jb_res <- tryCatch(
    tseries::jarque.bera.test(z),
    error = function(e) NULL
  )

  if (!is.null(jb_res)) {
    out["Jarque-Bera statistic"] <- as.numeric(jb_res$statistic)
    out["Jarque-Bera p-value"]   <- jb_res$p.value
  }

  sw_res <- tryCatch(
    shapiro.test(z),
    error = function(e) NULL
  )

  if (!is.null(sw_res)) {
    out["Shapiro-Wilk W"]       <- as.numeric(sw_res$statistic)
    out["Shapiro-Wilk p-value"] <- sw_res$p.value
  }

  out
}

asset_names <- setdiff(colnames(data_merged), "SP500")

event_test_results <- list()
panel_test_results <- list()

for (ev_name in names(events)) {

  ev_date <- events[[ev_name]]
  day_data <- get_event_window(data_merged, ev_date)

  test_names <- c(
    "Jarque-Bera statistic",
    "Jarque-Bera p-value",
    "Shapiro-Wilk W",
    "Shapiro-Wilk p-value"
  )

  # Preserve the output structure even when the window, demarketing, or LiNGAM step fails.
  if (NROW(day_data) == 0) {
    test_mat <- matrix(NA_real_, nrow = length(test_names), ncol = length(asset_names))
    rownames(test_mat) <- test_names
    colnames(test_mat) <- asset_names

  } else {
    demarketed_X <- get_demarketed_matrix(day_data, asset_names)

    if (is.null(demarketed_X)) {
      test_mat <- matrix(NA_real_, nrow = length(test_names), ncol = length(asset_names))
      rownames(test_mat) <- test_names
      colnames(test_mat) <- asset_names

    } else {
      lingam_E <- get_lingam_residuals_from_R2_network(demarketed_X)

      if (is.null(lingam_E)) {
        test_mat <- matrix(NA_real_, nrow = length(test_names), ncol = length(asset_names))
        rownames(test_mat) <- test_names
        colnames(test_mat) <- asset_names

      } else {
        test_mat <- sapply(asset_names, function(tkr) {
          compute_test_results_from_vector(lingam_E[, tkr])
        })
      }
    }
  }

  event_test_results[[ev_name]] <- as.data.frame(test_mat)

  # Summarize the cross-section for each event and each test statistic.
  panel_mat <- t(apply(test_mat, 1, function(z) {
    z <- as.numeric(z)
    c(
      median = if (all(is.na(z))) NA_real_ else median(z, na.rm = TRUE),
      min    = if (all(is.na(z))) NA_real_ else min(z, na.rm = TRUE),
      max    = if (all(is.na(z))) NA_real_ else max(z, na.rm = TRUE)
    )
  }))

  panel_test_results[[ev_name]] <- panel_mat
}

final_test_table <- cbind(
  panel_test_results[[1]],
  panel_test_results[[2]],
  panel_test_results[[3]],
  panel_test_results[[4]]
)

final_test_table <- as.data.frame(final_test_table)

colnames(final_test_table) <- c(
  "ChatGPT_median", "ChatGPT_min", "ChatGPT_max",
  "Hopper_median", "Hopper_min", "Hopper_max",
  "Copilot_median", "Copilot_min", "Copilot_max",
  "Claude_median", "Claude_min", "Claude_max"
)

rownames(final_test_table) <- c(
  "Jarque-Bera statistic",
  "Jarque-Bera p-value",
  "Shapiro-Wilk W",
  "Shapiro-Wilk p-value"
)

print(round(final_test_table, 6))

library(xts)

events <- c(
  "ChatGPT launch" = "2022-11-30",
  "NVIDIA Hopper/H100" = "2022-03-22",
  "GitHub Copilot preview" = "2021-06-29",
  "Claude Opus 4.6" = "2026-02-05"
)

# The market index is needed for the demarketing step before the VAR comparison.
stopifnot("SP500" %in% colnames(data_merged))

stock_names <- setdiff(colnames(data_merged), "SP500")

# For the BIC comparison, use the full event calendar day rather than only
# the regular trading-session subset.
get_one_day <- function(x, date_string) {
  x[paste0(date_string, " 00:00:00/", date_string, " 23:59:59")]
}

# Construct the demarketed return matrix for one event day.
get_demarketed_day <- function(day_data, stock_names) {
  needed <- c(stock_names, "SP500")
  sub <- day_data[, needed]
  sub <- na.omit(sub)

  if (NROW(sub) < 5) return(NULL)

  market_vec <- as.numeric(sub[, "SP500"])

  eps_mat <- sapply(stock_names, function(nm) {
    y <- as.numeric(sub[, nm])

    fit <- tryCatch(
      lm(y ~ market_vec),
      error = function(e) NULL
    )

    if (is.null(fit)) {
      rep(NA_real_, length(y))
    } else {
      resid(fit)
    }
  })

  eps_mat <- as.matrix(eps_mat)
  colnames(eps_mat) <- stock_names
  eps_mat <- eps_mat[complete.cases(eps_mat), , drop = FALSE]

  if (nrow(eps_mat) < 5 || ncol(eps_mat) < 2) return(NULL)

  eps_mat
}

# Compute BIC values for VAR(0) and VAR(1) using the same multivariate sample.
compute_var0_var1_bic <- function(X, type = "const") {
  X <- as.matrix(X)

  if (nrow(X) < 5 || ncol(X) < 2) return(NULL)

  K <- ncol(X)
  N <- nrow(X)

  # VAR(1) uses X_{t-1} to explain X_t. The same effective sample is used
  # for both VAR(0) and VAR(1), which makes the BIC comparison consistent.
  Y <- X[-1, , drop = FALSE]
  Xlag <- X[-N, , drop = FALSE]
  T_eff <- nrow(Y)

  if (T_eff < 3) return(NULL)

  # Fit a multivariate system equation-by-equation using least squares.
  fit_system <- function(Ymat, Xreg = NULL, type = "const") {
    if (is.null(Xreg)) {
      if (type == "const") {
        Z <- matrix(1, nrow = nrow(Ymat), ncol = 1)
      } else if (type == "none") {
        Z <- NULL
      } else {
        stop("Only type = 'const' or type = 'none' is implemented.")
      }
    } else {
      if (type == "const") {
        Z <- cbind(1, Xreg)
      } else if (type == "none") {
        Z <- Xreg
      } else {
        stop("Only type = 'const' or type = 'none' is implemented.")
      }
    }

    if (is.null(Z)) {
      U <- Ymat
      n_star <- 0
    } else {
      Bhat <- tryCatch(qr.solve(Z, Ymat), error = function(e) NULL)
      if (is.null(Bhat)) return(NULL)

      U <- Ymat - Z %*% Bhat
      n_star <- ncol(Z)
    }

    # Residual covariance and log-determinant term of the Gaussian system likelihood.
    Sigma_u <- crossprod(U) / nrow(Ymat)

    det_out <- tryCatch(
      determinant(Sigma_u, logarithm = TRUE),
      error = function(e) NULL
    )
    if (is.null(det_out) || !is.finite(as.numeric(det_out$modulus))) return(NULL)

    logdet <- as.numeric(det_out$modulus)

    # BIC per observation for a K-variable system with n_star regressors per equation.
    BIC_val <- logdet + (log(nrow(Ymat)) / nrow(Ymat)) * K * n_star

    list(
      T_eff = nrow(Ymat),
      K = K,
      n_star = n_star,
      logdet = logdet,
      BIC = BIC_val
    )
  }

  var0 <- fit_system(Ymat = Y, Xreg = NULL, type = type)
  var1 <- fit_system(Ymat = Y, Xreg = Xlag, type = type)

  if (is.null(var0) || is.null(var1)) return(NULL)

  list(
    var0 = var0,
    var1 = var1
  )
}

bic_results <- lapply(names(events), function(ev_name) {
  d <- events[[ev_name]]

  raw_day <- get_one_day(data_merged, d)
  if (NROW(raw_day) == 0) {
    return(data.frame(
      event = ev_name,
      date = d,
      T_eff = NA_integer_,
      K = NA_integer_,
      BIC_VAR0 = NA_real_,
      BIC_VAR1 = NA_real_,
      preferred_BIC = NA_character_
    ))
  }

  demarketed_day <- get_demarketed_day(raw_day, stock_names)
  if (is.null(demarketed_day)) {
    return(data.frame(
      event = ev_name,
      date = d,
      T_eff = NA_integer_,
      K = NA_integer_,
      BIC_VAR0 = NA_real_,
      BIC_VAR1 = NA_real_,
      preferred_BIC = NA_character_
    ))
  }

  bic <- compute_var0_var1_bic(demarketed_day, type = "const")
  if (is.null(bic)) {
    return(data.frame(
      event = ev_name,
      date = d,
      T_eff = NA_integer_,
      K = NA_integer_,
      BIC_VAR0 = NA_real_,
      BIC_VAR1 = NA_real_,
      preferred_BIC = NA_character_
    ))
  }

  data.frame(
    event = ev_name,
    date = d,
    T_eff = bic$var0$T_eff,
    K = bic$var0$K,
    BIC_VAR0 = bic$var0$BIC,
    BIC_VAR1 = bic$var1$BIC,
    preferred_BIC = ifelse(bic$var1$BIC < bic$var0$BIC, "VAR(1)", "VAR(0)")
  )
})

bic_results <- do.call(rbind, bic_results)
print(bic_results, digits = 6)