Network of Each Event

# ============================================================
# Required package: R2DAG
# ============================================================
# The replication code uses functions from the author's R package
# R2DAG. The package is available from GitHub.
#
# The arguments below are deliberately conservative:
#   - dependencies = FALSE avoids reinstalling packages that may already
#     be available in the user's local R environment;
#   - upgrade = "never" prevents automatic package upgrades that could
#     change the replication environment.
#
# If R2DAG is already installed, this installation cell can be skipped.

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

library(readr)
library(xts)

# Read the cleaned return dataset and convert it to an xts object.
# The first column is assumed to contain the timestamp. All remaining
# columns are kept as return series and ordered by time.
load_stocks_data <- function(path) {

  df <- read_csv(path, show_col_types = FALSE)

  # Use the first column as the date-time index.
  dt_col <- colnames(df)[1]
  df[[dt_col]] <- as.POSIXct(df[[dt_col]], tz = "UTC")

  # Convert the rectangular data frame to an xts time-series object.
  x <- xts(df[, -1, drop = FALSE], order.by = df[[dt_col]])

  # Ensure chronological ordering even if the CSV is not sorted.
  x <- x[order(index(x))]

  return(x)
}

# Local path to the cleaned one-minute return dataset.
# Replace this path if the replication folder is stored elsewhere.
path = "[YOUR FOLDER]/full_dataset.csv"

# Load the complete dataset as an xts object.
all_data <- load_stocks_data(path)

# ============================================================
# Plot one total connectedness matrix as a two-panel network
# ============================================================
# Input:
#   tab  - a square total connectedness matrix, where rows are
#          receivers and columns are transmitters/sources.
#   Name - optional output filename. If provided, the plot is saved
#          as a PNG file.
#
# Output:
#   The function invisibly returns the qgraph objects and the matrices
#   used for plotting. This is useful for checking the exact thresholded
#   matrices behind the figures.
plot_fixed_basket_network <- function(
  tab,
  Name = NULL,
  label = NULL,
  output_dir = ".",
  env = parent.frame(),

  pct_keep = 1,
  top_n_avg = FALSE,
  top = 250,
  limit = NULL,
  group_limit = 0,

  vsize_range = c(11, 22),
  group_vsize_range = c(18, 30),
  edge_size = 8,
  group_edge_size = 14,
  arrow_size = 8,
  pie_border = 0.2,

  png_width = 4200,
  png_height = 2000,
  middle_gap = 0.12,

  reorder_to_fixed_basket = TRUE,
  include_non_ai = TRUE,
  highlight_top_sources = TRUE,
  source_share_cutoff = 0.60,
  source_max_share = 0.10
) {

  # ============================================================
  # 0. Package checks
  # ============================================================

  needed_packages <- c("qgraph", "scales")
  missing_packages <- needed_packages[
    !vapply(needed_packages, requireNamespace, logical(1), quietly = TRUE)
  ]

  if (length(missing_packages) > 0) {
    stop(
      "Missing R package(s): ", paste(missing_packages, collapse = ", "),
      "\nInstall them with:\ninstall.packages(c(",
      paste(sprintf('"%s"', missing_packages), collapse = ", "),
      "))"
    )
  }

  # ============================================================
  # 1. Handle multiple matrices if provided
  # ============================================================

  # `tab` can also be a character vector containing object names.
  # In that case, the objects are retrieved from `env`.
  if (is.character(tab)) {
    obj_names <- tab

    missing_objects <- obj_names[!vapply(obj_names, exists, logical(1), envir = env)]
    if (length(missing_objects) > 0) {
      stop(
        "The following objects were not found in the specified environment: ",
        paste(missing_objects, collapse = ", ")
      )
    }

    tab <- mget(obj_names, envir = env)
    names(tab) <- obj_names
  }

  # If a list of matrices is provided, the function recursively
  # produces one plot per matrix.
  if (is.list(tab) && !is.data.frame(tab)) {
    if (is.null(names(tab)) || any(!nzchar(names(tab)))) {
      names(tab) <- paste0("matrix_", seq_along(tab))
    }

    safe_names <- gsub("[^A-Za-z0-9_\\-]+", "_", names(tab))

    if (is.null(Name)) {
      if (!dir.exists(output_dir)) {
        dir.create(output_dir, recursive = TRUE)
      }
      Names <- file.path(output_dir, paste0(safe_names, "_two_panel_network.png"))
    } else if (length(Name) == length(tab)) {
      Names <- Name
    } else {
      stop("If multiple matrices are provided, `Name` must be NULL or have the same length as the matrix list.")
    }

    out <- vector("list", length(tab))
    names(out) <- names(tab)

    for (i in seq_along(tab)) {
      out[[i]] <- plot_fixed_basket_network(
        tab = tab[[i]],
        Name = Names[i],
        label = names(tab)[i],
        output_dir = output_dir,
        env = env,
        pct_keep = pct_keep,
        top_n_avg = top_n_avg,
        top = top,
        limit = limit,
        group_limit = group_limit,
        vsize_range = vsize_range,
        group_vsize_range = group_vsize_range,
        edge_size = edge_size,
        group_edge_size = group_edge_size,
        arrow_size = arrow_size,
        pie_border = pie_border,
        png_width = png_width,
        png_height = png_height,
        middle_gap = middle_gap,
        reorder_to_fixed_basket = reorder_to_fixed_basket,
        include_non_ai = include_non_ai,
        highlight_top_sources = highlight_top_sources,
        source_share_cutoff = source_share_cutoff,
        source_max_share = source_max_share
      )
    }

    return(invisible(out))
  }

  # ============================================================
  # 2. Fixed stock basket and groups
  # ============================================================

  # The order below fixes the visual position of the same ticker
  # across all event-window figures.
  tickers_m7 <- c("NVDA", "MSFT", "AMZN", "META", "GOOG", "TSLA", "AAPL")
  tickers_computing <- c("AVGO", "AMD", "MU")
  tickers_platform <- c("ORCL", "IBM", "PLTR", "CSCO")
  tickers_disruption_adoption <- c("CRM", "NOW", "INTU", "WDAY", "ADBE")
  non_ai_tickers <- c("WMT", "COST", "CVX", "ABBV", "PG", "MRK", "PM", "MCD", "PEP", "T")

  fixed_order <- c(
    tickers_m7,
    tickers_computing,
    tickers_platform,
    tickers_disruption_adoption
  )

  if (include_non_ai) {
    fixed_order <- c(fixed_order, non_ai_tickers)
  }

  # ============================================================
  # 3. Internal helper functions
  # ============================================================

  # Convert a directed matrix into its pairwise net representation.
  # For each pair (i, j), only the stronger direction is kept and the
  # weaker direction is subtracted.
  net_pairwise <- function(X) {
    X <- as.matrix(X)
    Y <- X
    diag(Y) <- 0

    if (nrow(Y) >= 2) {
      for (i in seq_len(nrow(Y) - 1)) {
        for (j in (i + 1):ncol(Y)) {
          a_ij <- Y[i, j]
          a_ji <- Y[j, i]

          if (a_ij >= a_ji) {
            Y[i, j] <- a_ij - a_ji
            Y[j, i] <- 0
          } else {
            Y[j, i] <- a_ji - a_ij
            Y[i, j] <- 0
          }
        }
      }
    }

    Y
  }

  # Compute node sizes and pie-chart shares.
  # The size reflects total connectedness involvement, while the pie
  # share separates outgoing and incoming strength.
  size_pie_map_ir <- function(M, NAMES, rescale_ = c(7.2, 13.2)) {
    M <- as.matrix(M)
    colnames(M) <- NAMES
    rownames(M) <- NAMES
    diag(M) <- 0

    to_strength <- rowSums(M, na.rm = TRUE)
    from_strength <- colSums(M, na.rm = TRUE)

    total_strength <- to_strength + from_strength
    pie_map <- to_strength / total_strength
    pie_map[!is.finite(pie_map)] <- 0.5

    if (max(total_strength, na.rm = TRUE) == min(total_strength, na.rm = TRUE)) {
      vsize_map <- rep(mean(rescale_), length(total_strength))
    } else {
      vsize_map <- scales::rescale(total_strength, to = rescale_)
    }

    names(vsize_map) <- NAMES
    names(pie_map) <- NAMES

    list(
      vsize_map = vsize_map,
      pie_map = pie_map,
      to_strength = to_strength,
      from_strength = from_strength
    )
  }

  # Threshold edges before plotting. By default, the function keeps the
  # strongest share of edges determined by `pct_keep`; alternatively it
  # can use the average of the top `top` entries.
  thrsh_topmean <- function(data, top_n_avg = FALSE, top = 250, pct = 0.15) {
    data <- as.matrix(data)

    pct <- max(min(pct, 1), 0)

    if (top_n_avg) {
      vals <- sort(as.vector(data), decreasing = TRUE)
      vals <- vals[is.finite(vals) & vals > 0]

      if (length(vals) == 0) {
        limit <- 0
      } else {
        top <- min(top, length(vals))
        limit <- mean(vals[1:top])
      }
    } else {
      threshold <- 1 - pct
      pos <- ifelse(data == 0, NA, data)
      limit <- suppressWarnings(quantile(pos, threshold, na.rm = TRUE))
      if (!is.finite(limit)) limit <- 0
    }

    list(
      data_thrsh = ifelse(data < limit, 0, data),
      limit = limit
    )
  }

  # Build a fixed circular layout. Groups are placed on a larger circle,
  # and stocks within each group are arranged around the group center.
  make_layouts <- function(groups, n_nodes) {
    group_angles <- seq(pi / 2, pi / 2 + 2 * pi, length.out = length(groups) + 1)
    group_angles <- group_angles[-length(group_angles)]

    group_radius <- 5.7
    within_spacing <- 2.8

    node_layout <- matrix(NA_real_, nrow = n_nodes, ncol = 2)
    group_layout <- matrix(NA_real_, nrow = length(groups), ncol = 2)
    rownames(group_layout) <- names(groups)

    for (g in seq_along(groups)) {
      ids <- groups[[g]]
      k <- length(ids)

      cx <- group_radius * cos(group_angles[g])
      cy <- group_radius * sin(group_angles[g])

      group_layout[g, ] <- c(cx, cy)

      if (k == 1) {
        node_layout[ids, ] <- c(cx, cy)
      } else {
        radius_g <- max(within_spacing * k / (2 * pi), 2.2)

        node_angles <- seq(pi / 2, pi / 2 + 2 * pi, length.out = k + 1)
        node_angles <- node_angles[-length(node_angles)]

        node_layout[ids, 1] <- cx + radius_g * cos(node_angles)
        node_layout[ids, 2] <- cy + radius_g * sin(node_angles)
      }
    }

    list(
      node_layout = node_layout,
      group_layout = group_layout
    )
  }

  # Aggregate firm-level edges into group-to-group edges.
  # The convention is preserved: rows are receivers and columns are senders.
  make_aggregated_group_matrix <- function(M, groups) {
    group_names <- names(groups)

    G <- matrix(0, nrow = length(groups), ncol = length(groups))
    rownames(G) <- group_names
    colnames(G) <- group_names

    for (receiver_g in seq_along(groups)) {
      for (sender_g in seq_along(groups)) {
        if (receiver_g != sender_g) {
          receivers <- groups[[receiver_g]]
          senders <- groups[[sender_g]]

          G[receiver_g, sender_g] <- sum(
            M[receivers, senders, drop = FALSE],
            na.rm = TRUE
          )
        }
      }
    }

    diag(G) <- 0
    G
  }

  # ============================================================
  # 4. Matrix preparation
  # ============================================================

  # Convert to a numeric matrix and enforce a consistent ticker order.
  M <- as.matrix(tab)
  suppressWarnings(storage.mode(M) <- "numeric")

  if (nrow(M) != ncol(M)) {
    stop("The provided object is not a square matrix: ", nrow(M), " x ", ncol(M))
  }

  if (is.null(colnames(M)) && is.null(rownames(M))) {
    if (ncol(M) > length(fixed_order)) {
      stop("The matrix has no row or column names and is larger than the fixed stock basket.")
    }

    colnames(M) <- fixed_order[seq_len(ncol(M))]
    rownames(M) <- fixed_order[seq_len(nrow(M))]
  }

  if (is.null(colnames(M)) && !is.null(rownames(M))) {
    colnames(M) <- rownames(M)
  }

  if (!is.null(colnames(M)) && is.null(rownames(M))) {
    rownames(M) <- colnames(M)
  }

  if (!setequal(rownames(M), colnames(M))) {
    stop("The row names and column names of the matrix do not contain the same tickers.")
  }

  M <- M[colnames(M), colnames(M), drop = FALSE]

  if (reorder_to_fixed_basket) {
    ordered_names <- c(
      fixed_order[fixed_order %in% colnames(M)],
      setdiff(colnames(M), fixed_order)
    )

    M <- M[ordered_names, ordered_names, drop = FALSE]
  }

  # Missing or infinite values would break the plotting routine, so they
  # are replaced by zero after the matrix checks.
  if (any(!is.finite(M))) {
    warning("The matrix contained NA/NaN/Inf values. These values were replaced by zero.")
    M[!is.finite(M)] <- 0
  }

  names_ <- colnames(M)
  diag(M) <- 0

  # ============================================================
  # 5. Group definitions
  # ============================================================

  # These groups match the categories used in the empirical section.
  groups <- list(
    M7 = which(names_ %in% tickers_m7),
    Computing = which(names_ %in% tickers_computing),
    Platform = which(names_ %in% tickers_platform),
    `Disruption/Adoption` = which(names_ %in% tickers_disruption_adoption)
  )

  if (include_non_ai) {
    groups$`Non-AI` <- which(names_ %in% non_ai_tickers)
  }

  assigned <- sort(unique(unlist(groups)))
  other <- setdiff(seq_along(names_), assigned)

  if (length(other) > 0) {
    groups$Other <- other
  }

  groups <- groups[sapply(groups, length) > 0]

  group_colors <- c(
    M7 = "#197EC0",
    Computing = "#FED439",
    Platform =  "#8A9197",
    `Disruption/Adoption` = "#1A9993",
    `Non-AI` = "white",
    Other = "gray85"
  )

  group_colors <- group_colors[names(groups)]
  group_colors[is.na(group_colors)] <- "gray85"

  edge_group_colors <- group_colors
  edge_group_colors[edge_group_colors == "white"] <- "gray45"

  # ============================================================
  # 6. Pairwise-netted firm-level matrix
  # ============================================================

  # The firm-level plot is based on pairwise net flows. This avoids
  # displaying two arrows between the same pair of stocks.
  M_net <- net_pairwise(M)

  # qgraph uses the opposite orientation convention here:
  # M[receiver, sender], therefore A_plot[sender, receiver]
  A_plot <- t(M_net)
  diag(A_plot) <- 0

  if (is.null(limit)) {
    limit <- thrsh_topmean(
      A_plot,
      top_n_avg = top_n_avg,
      top = top,
      pct = pct_keep
    )$limit
  }

  A_thrsh <- ifelse(A_plot < limit, 0, A_plot)

  max_val <- max(A_thrsh, na.rm = TRUE)
  if (!is.finite(max_val) || max_val <= 0) max_val <- 1

  maps <- size_pie_map_ir(M, names_, vsize_range)

  node_color <- rep(NA_character_, length(names_))
  for (g in seq_along(groups)) {
    node_color[groups[[g]]] <- group_colors[g]
  }

  # Optionally highlight the strongest net source nodes in light pink.
  if (highlight_top_sources) {
    node_net <- maps$from_strength - maps$to_strength

    max_node_net <- max(node_net, na.rm = TRUE)
    max_source_nodes <- max(1, floor(length(node_net) * source_max_share + 0.5))

    hi_node <- rep(FALSE, length(node_net))

    if (is.finite(max_node_net) && max_node_net > 0) {
      candidate <- node_net > 0 & node_net >= source_share_cutoff * max_node_net
      candidate_order <- order(node_net, decreasing = TRUE)
      selected <- candidate_order[candidate[candidate_order]]
      selected <- head(selected, max_source_nodes)
      hi_node[selected] <- TRUE
    }

    node_color[hi_node] <- "lightpink1"
  }

  label_color_node <- rep("black", length(names_))

  layouts <- make_layouts(groups, length(names_))
  layout_custom <- layouts$node_layout
  group_layout <- layouts$group_layout

  edge_col_mat <- matrix(
    NA_character_,
    nrow = length(names_),
    ncol = length(names_)
  )
  rownames(edge_col_mat) <- names_
  colnames(edge_col_mat) <- names_

  for (g in seq_along(groups)) {
    edge_col_mat[groups[[g]], ] <- edge_group_colors[g]
  }

  diag(edge_col_mat) <- NA

  # ============================================================
  # 7. Pairwise-netted group-level matrix
  # ============================================================

  # The right panel aggregates the original firm-level matrix to group
  # level and then applies the same pairwise-netting logic.
  G_raw <- make_aggregated_group_matrix(M, groups)
  G_net <- net_pairwise(G_raw)

  A_group_plot <- t(G_net)
  diag(A_group_plot) <- 0

  if (is.null(group_limit)) {
    group_limit <- 0
  }

  A_group_thrsh <- ifelse(A_group_plot < group_limit, 0, A_group_plot)

  group_max_val <- max(A_group_thrsh, na.rm = TRUE)
  if (!is.finite(group_max_val) || group_max_val <= 0) group_max_val <- 1

  group_names <- names(groups)
  group_maps <- size_pie_map_ir(G_raw, group_names, group_vsize_range)

  group_groups <- as.list(seq_along(group_names))
  names(group_groups) <- group_names

  group_edge_col <- matrix(
    NA_character_,
    nrow = length(group_names),
    ncol = length(group_names)
  )
  rownames(group_edge_col) <- group_names
  colnames(group_edge_col) <- group_names

  for (i in seq_along(group_names)) {
    group_edge_col[i, ] <- edge_group_colors[i]
  }

  diag(group_edge_col) <- NA

  group_labels <- group_names
  group_labels[group_labels == "Computing"] <- "Comp"
  group_labels[group_labels == "Platform"] <- "Plat"
  group_labels[group_labels == "Disruption/Adoption"] <- "A/D"

  # ============================================================
  # 8. Plot construction
  # ============================================================

  if (!is.null(Name)) {
    png(Name, width = png_width, height = png_height, res = 180)
    on.exit(dev.off(), add = TRUE)
  }

  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par), add = TRUE)

  layout(
    matrix(c(1, 0, 2), nrow = 1),
    widths = c(1, middle_gap, 1)
  )

  par(mar = c(1, 1, 1, 1))

  # Left panel: firm-level network.
  qg_node <- qgraph::qgraph(
    A_thrsh,
    groups = groups,
    layout = layout_custom,
    layoutScale = c(1.08, 1.08),
    maximum = max_val,
    edge.labels = FALSE,
    label.font = rep(2, length(names_)),
    label.cex = rep(1.1, length(names_)),
    shape = "circle",
    labels = names_,
    curveAll = TRUE,
    color = node_color,
    edge.color = edge_col_mat,
    fade = FALSE,
    pieBorder = pie_border,
    label.fill.vertical = 0.9,
    label.fill.horizontal = 0.9,
    curve = 1.35,
    asize = arrow_size,
    arrows = TRUE,
    bidirectional = FALSE,
    vsize = maps$vsize_map,
    label.color = label_color_node,
    legend = FALSE,
    esize = edge_size,
    pie = maps$pie_map,
    pieColor = rep("green3", length(names_)),
    pieColor2 = rep("red2", length(names_))
  )

  # Right panel: group-level network.
  qg_group <- qgraph::qgraph(
    A_group_thrsh,
    groups = group_groups,
    layout = group_layout,
    layoutScale = c(0.85, 0.85),
    maximum = group_max_val,
    edge.labels = FALSE,
    label.font = rep(2, length(group_names)),
    label.cex = rep(1.1, length(group_names)),
    shape = "circle",
    labels = group_labels,
    curveAll = TRUE,
    color = group_colors,
    edge.color = group_edge_col,
    fade = FALSE,
    pieBorder = pie_border,
    label.fill.vertical = 0.9,
    label.fill.horizontal = 0.9,
    curve = 1.15,
    asize = arrow_size,
    arrows = TRUE,
    bidirectional = FALSE,
    vsize = group_maps$vsize_map,
    label.color = rep("black", length(group_names)),
    legend = FALSE,
    esize = group_edge_size,
    pie = group_maps$pie_map,
    pieColor = rep("green3", length(group_names)),
    pieColor2 = rep("red2", length(group_names))
  )

  if (!is.null(Name)) {
    message("Picture saved")
  }

  invisible(list(
    qgraph_object_firm = qg_node,
    qgraph_object_group = qg_group,
    original_matrix = M,
    pairwise_net_matrix = M_net,
    plot_matrix = A_plot,
    thresholded_plot_matrix = A_thrsh,
    group_raw_matrix = G_raw,
    group_net_matrix = G_net,
    group_plot_matrix = A_group_plot,
    group_thresholded_plot_matrix = A_group_thrsh,
    limit = limit,
    group_limit = group_limit,
    pie_map = maps$pie_map,
    group_pie_map = group_maps$pie_map,
    vsize_map = maps$vsize_map,
    group_vsize_map = group_maps$vsize_map,
    groups = groups,
    layout = layout_custom,
    group_layout = group_layout,
    label = label
  ))
}

library(xts)
library(lubridate)

# =========================================================
# 0. Event definitions
# =========================================================

# Each event has:
#   - a short name used in output object names;
#   - a calendar date used for the 1-day and 3-day windows;
#   - an exact Eastern Time timestamp used for the 24-hour window.
events_info <- list(
  "ChatGPT launch" = list(
    short_name = "ChatGPT",
    date = as.Date("2022-11-30"),
    datetime = as.POSIXct("2022-11-30 15:14:00", tz = "America/New_York")
  ),
  "NVIDIA Hopper/H100" = list(
    short_name = "Hopper_H100",
    date = as.Date("2022-03-22"),
    datetime = as.POSIXct("2022-03-22 11:27:02", tz = "America/New_York")
  ),
  "GitHub Copilot preview" = list(
    short_name = "Copilot",
    date = as.Date("2021-06-29"),
    datetime = as.POSIXct("2021-06-29 11:53:00", tz = "America/New_York")
  ),
  "Claude Opus 4.6" = list(
    short_name = "Claude_Opus_46",
    date = as.Date("2026-02-05"),
    datetime = as.POSIXct("2026-02-05 12:45:00", tz = "America/New_York")
  )
)


# =========================================================
# 1. Helper: remove the market effect via SP500 regression
# =========================================================

# The network is estimated on demarketed stock returns. For each stock,
# this function estimates a simple regression on the SP500 return and
# returns the residual series.
remove_market_effect <- function(data) {

  if (!xts::is.xts(data)) {
    stop("data must be an xts object.")
  }

  if (!"SP500" %in% colnames(data)) {
    stop("The input data must contain an SP500 column.")
  }

  # Convert to a data frame because lm() works naturally with vectors.
  df <- as.data.frame(data)
  market <- df$SP500

  stocks <- df[, colnames(df) != "SP500", drop = FALSE]

  resid_mat <- sapply(names(stocks), function(col) {

    y <- stocks[[col]]

    # Keep only finite observations for the stock and the market return.
    ok <- is.finite(y) & is.finite(market)

    if (sum(ok) < 3) {
      stop(paste("Too few valid observations for ticker:", col))
    }

    out <- rep(NA_real_, length(y))
    out[ok] <- resid(lm(y[ok] ~ market[ok]))

    return(out)
  })

  # Put the residual matrix back into an xts object with the original index.
  resid_xts <- xts(resid_mat, order.by = index(data))
  resid_xts <- resid_xts[complete.cases(resid_xts), ]

  return(resid_xts)
}


# =========================================================
# 2. Helper: slice a full trading day
# =========================================================

# The one-day window is the full regular trading session on the event date.
# The default market hours are 09:30--16:00 Eastern Time.
slice_full_trading_day <- function(data,
                                   event_date,
                                   tz = "America/New_York",
                                   market_open = "09:30:00",
                                   market_close = "16:00:00") {

  if (!xts::is.xts(data)) {
    stop("data must be an xts object.")
  }

  # Build the start and end timestamps in Eastern Time.
  start_time <- as.POSIXct(
    paste(event_date, market_open),
    tz = tz
  )

  end_time <- as.POSIXct(
    paste(event_date, market_close),
    tz = tz
  )

  # Keep all observations inside the regular trading session.
  out <- data[index(data) >= start_time & index(data) <= end_time]

  if (nrow(out) == 0) {
    stop(paste("No data found for full trading day:", event_date))
  }

  cat("Window type:    full trading day\n")
  cat("Date ET:        ", as.character(event_date), "\n")
  cat("Window start:   ", format(start_time, tz = tz, usetz = TRUE), "\n")
  cat("Window end:     ", format(end_time, tz = tz, usetz = TRUE), "\n")
  cat("Actual obs:     ", nrow(out), "\n")
  cat("Data start ET:  ", format(first(index(out)), tz = tz, usetz = TRUE), "\n")
  cat("Data end ET:    ", format(last(index(out)), tz = tz, usetz = TRUE), "\n\n")

  return(out)
}


# =========================================================
# 3. Helper: slice event day plus the next two available trading days
# =========================================================

# The three-day window uses trading days rather than calendar days.
# This avoids including weekends or holidays when markets are closed.
slice_three_trading_days <- function(data,
                                     event_date,
                                     tz = "America/New_York",
                                     market_open = "09:30:00",
                                     market_close = "16:00:00") {

  if (!xts::is.xts(data)) {
    stop("data must be an xts object.")
  }

  event_date <- as.Date(event_date)

  # Identify available trading dates in the data starting from the event date.
  data_dates <- as.Date(index(data), tz = tz)
  available_dates <- sort(unique(data_dates[data_dates >= event_date]))

  if (length(available_dates) < 3) {
    stop(
      paste0(
        "Fewer than three available trading dates found from event date: ",
        event_date
      )
    )
  }

  # Use the event trading day and the next two available trading days.
  selected_dates <- available_dates[1:3]

  keep <- rep(FALSE, nrow(data))

  for (d in as.character(selected_dates)) {

    start_time <- as.POSIXct(
      paste(d, market_open),
      format = "%Y-%m-%d %H:%M:%S",
      tz = tz
    )

    end_time <- as.POSIXct(
      paste(d, market_close),
      format = "%Y-%m-%d %H:%M:%S",
      tz = tz
    )

    keep <- keep | (index(data) >= start_time & index(data) <= end_time)
  }

  out <- data[keep]

  if (nrow(out) == 0) {
    stop(paste("No data found for three-day trading window:", event_date))
  }

  cat("Window type:    event day + next two available trading days\n")
  cat("Dates ET:       ", paste(as.character(selected_dates), collapse = ", "), "\n")
  cat("Actual obs:     ", nrow(out), "\n")
  cat("Data start ET:  ", format(first(index(out)), tz = tz, usetz = TRUE), "\n")
  cat("Data end ET:    ", format(last(index(out)), tz = tz, usetz = TRUE), "\n\n")

  return(out)
}


# =========================================================
# 4. Helper: slice exact calendar hours from the event timestamp
# =========================================================

# The 24-hour window starts at the exact event timestamp and ends exactly
# 24 calendar hours later. Since the data contain trading minutes only,
# the actual number of observations depends on market hours.
slice_exact_calendar_hours <- function(data,
                                       event_datetime_et,
                                       hours_after = 24,
                                       tz = "America/New_York") {

  if (!xts::is.xts(data)) {
    stop("data must be an xts object.")
  }

  # Convert the event timestamp to Eastern Time and define the endpoint.
  event_time <- as.POSIXct(event_datetime_et, tz = tz)
  end_time <- event_time + hours_after * 3600

  out <- data[index(data) >= event_time & index(data) <= end_time]

  if (nrow(out) == 0) {
    stop(paste("No data in exact calendar-hour event window after:", event_datetime_et))
  }

  cat("Window type:    exact", hours_after, "calendar hours from event time\n")
  cat("Event time ET:  ", format(event_time, tz = tz, usetz = TRUE), "\n")
  cat("Window end ET:  ", format(end_time, tz = tz, usetz = TRUE), "\n")
  cat("Actual obs:     ", nrow(out), "\n")
  cat("Data start ET:  ", format(first(index(out)), tz = tz, usetz = TRUE), "\n")
  cat("Data end ET:    ", format(last(index(out)), tz = tz, usetz = TRUE), "\n\n")

  return(out)
}


# =========================================================
# 5. Main function: compute one event-window network
# =========================================================

# This function contains the full workflow for one event-window pair:
#   1. select the required ticker columns and SP500;
#   2. slice the requested event window;
#   3. remove the market effect;
#   4. estimate the directed R2 network and return the total table.
get_total_table_event_window <- function(
    data,
    event_date,
    event_datetime,
    window_type = c("1d", "3d", "24h"),

    tickers_m7 = c("NVDA", "MSFT", "AMZN", "META", "GOOG", "TSLA", "AAPL"),
    tickers_computing = c("AVGO", "AMD", "MU"),
    tickers_platform = c("ORCL", "IBM", "PLTR", "CSCO"),
    tickers_disruption_adoption = c("CRM", "NOW", "INTU", "WDAY", "ADBE"),
    non_ai_tickers = c(),

    dag_method = "lingam",
    tz = "America/New_York"
) {

  window_type <- match.arg(window_type)

  if (!xts::is.xts(data)) {
    stop("data must be an xts object.")
  }

  # The network itself is estimated only on stock returns.
  selected_tickers <- c(
    tickers_m7,
    tickers_computing,
    tickers_platform,
    tickers_disruption_adoption,
    non_ai_tickers
  )

  # SP500 is added only because it is needed for demarketing.
  needed_cols <- c(selected_tickers, "SP500")
  missing_cols <- setdiff(needed_cols, colnames(data))

  if (length(missing_cols) > 0) {
    stop(
      paste0(
        "These required columns are missing from data: ",
        paste(missing_cols, collapse = ", ")
      )
    )
  }

  data_sub <- data[, needed_cols]

  # Select the appropriate event window.
  if (window_type == "1d") {

    window_data <- slice_full_trading_day(
      data = data_sub,
      event_date = event_date,
      tz = tz
    )

  } else if (window_type == "3d") {

    window_data <- slice_three_trading_days(
      data = data_sub,
      event_date = event_date,
      tz = tz
    )

  } else if (window_type == "24h") {

    window_data <- slice_exact_calendar_hours(
      data = data_sub,
      event_datetime_et = event_datetime,
      hours_after = 24,
      tz = tz
    )
  }

  # Remove the market component before network estimation.
  demarketed_data <- remove_market_effect(window_data)

  demarketed_data <- demarketed_data[, selected_tickers]

  # Estimate the directed R2 connectedness matrix.
  total_table <- R2_network(
    demarketed_data,
    dag_method = dag_method
  )$total_table

  return(total_table)
}


# =========================================================
# 6. Run all events and create the twelve replication objects
# =========================================================

# The loop computes 4 events x 3 windows = 12 matrices.
# If assign_to_globalenv = TRUE, the matrices are also saved as
# objects named total_[event]_[window], which are used by the plotting cells.
run_event_total_tables_12 <- function(
    data,
    events_info,
    dag_method = "lingam",
    assign_to_globalenv = TRUE
) {

  out <- list()

  for (event_name in names(events_info)) {

    cat("\n========================================\n")
    cat("Event:", event_name, "\n")
    cat("========================================\n")

    event_date <- events_info[[event_name]]$date
    event_datetime <- events_info[[event_name]]$datetime
    short_name <- events_info[[event_name]]$short_name

    out[[short_name]] <- list()

    # Compute the three windows for the current event.
    for (window_type in c("1d", "3d", "24h")) {

      cat("\nRunning window:", window_type, "\n")

      total_table <- get_total_table_event_window(
        data = data,
        event_date = event_date,
        event_datetime = event_datetime,
        window_type = window_type,
        dag_method = dag_method
      )

      out[[short_name]][[window_type]] <- total_table

      # Object names are standardized so that the plotting cells below
      # can call them directly.
      object_name <- paste0("total_", short_name, "_", window_type)

      if (assign_to_globalenv) {
        assign(
          x = object_name,
          value = total_table,
          envir = .GlobalEnv
        )
      }

      cat("Created object:", object_name, "\n")
    }
  }

  if (assign_to_globalenv) {
    assign(
      x = "event_total_tables_12",
      value = out,
      envir = .GlobalEnv
    )

    cat("\nCreated list object: event_total_tables_12\n")
  }

  return(invisible(out))
}


# =========================================================
# 7. Execute the full event-window loop
# =========================================================

# Running this line creates:
#   - the nested list event_total_tables_12;
#   - twelve individual matrices in the global environment.
event_total_tables_12 <- run_event_total_tables_12(
  data = all_data,
  events_info = events_info,
  dag_method = "lingam",
  assign_to_globalenv = TRUE
)

plot_and_show_network <- function(tab, Name, ...) {
  plot_fixed_basket_network(
    tab = tab,
    Name = Name,
    ...
  )

  if (requireNamespace("IRdisplay", quietly = TRUE)) {
    IRdisplay::display_png(file = Name)
  } else {
    warning("Install IRdisplay to show PNG files inside the notebook: install.packages('IRdisplay')")
  }
}

plot_and_show_network(
  tab = total_ChatGPT_1d,
  Name = "ChatGPT_1d_network.png"
)

# Save the ChatGPT 3-day network figure.
plot_and_show_network(
  tab = total_ChatGPT_3d,
  Name = "ChatGPT_3d_network.png"
)

# Save the ChatGPT exact 24-hour network figure.
plot_and_show_network(
  tab = total_ChatGPT_24h,
  Name = "ChatGPT_24h_network.png"
)

# Save the Hopper/H100 1-day network figure.
plot_and_show_network(
  tab = total_Hopper_H100_1d,
  Name = "Hopper_H100_1d_network.png"
)

# Save the Hopper/H100 3-day network figure.
plot_and_show_network(
  tab = total_Hopper_H100_3d,
  Name = "Hopper_H100_3d_network.png"
)

# Save the Hopper/H100 exact 24-hour network figure.
plot_and_show_network(
  tab = total_Hopper_H100_24h,
  Name = "Hopper_H100_24h_network.png"
)

# Save the Copilot 1-day network figure.
plot_and_show_network(
  tab = total_Copilot_1d,
  Name = "Copilot_1d_network.png"
)

# Save the Copilot 3-day network figure.
plot_and_show_network(
  tab = total_Copilot_3d,
  Name = "Copilot_3d_network.png"
)

# Save the Copilot exact 24-hour network figure.
plot_and_show_network(
  tab = total_Copilot_24h,
  Name = "Copilot_24h_network.png"
)

# Save the Claude Opus 4.6 1-day network figure.
plot_and_show_network(
  tab = total_Claude_Opus_46_1d,
  Name = "Claude_Opus_46_1d_network.png"
)

# Save the Claude Opus 4.6 3-day network figure.
plot_and_show_network(
  tab = total_Claude_Opus_46_3d,
  Name = "Claude_Opus_46_3d_network.png"
)

# Save the Claude Opus 4.6 exact 24-hour network figure.
plot_and_show_network(
  tab = total_Claude_Opus_46_24h,
  Name = "Claude_Opus_46_24h_network.png"
)