Placebo Distribution

# ============================================================
# Required package: R2DAG
# ============================================================
# The replication codes use functions from the author's R package
# R2DAG. The package is available from GitHub and can be installed
# with the command below.
#
# `dependencies = FALSE` avoids installing a large set of additional
# packages automatically, and `upgrade = "never"` prevents R from
# changing already installed package versions during replication.

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

# Load the package after installation. The main function used below is
# R2_network(), which returns the directed R2 network table.
library(R2DAG)

# Packages used for reading the prepared CSV file and handling
# time-indexed intraday observations.
library(readr)
library(xts)

# Replace [YOUR FOLDER] with the folder that contains the full
# one-minute log-return dataset created in the data-preparation notebook.
df <- read_csv("[YOUR FOLDER]/full_dataset.csv")

# The event windows in the paper are defined in Eastern Time, so the
# timestamp is explicitly converted to the America/New_York timezone.
df$datetime <- as.POSIXct(df$datetime, tz = "America/New_York")

# Store the return panel as an xts object. The first column is the
# timestamp, while all remaining columns are ticker return series.
all_data <- xts(df[, -1], order.by = df$datetime)

# ============================================================
# Ticker groups used in the placebo exercise
# ============================================================
# The AI-related basket is split into the groups used throughout the
# empirical analysis. The non-AI tickers are used only as placebo
# replacements for the exposed group.
default_tickers_m7 <- c("NVDA", "MSFT", "AMZN", "META", "GOOG", "TSLA", "AAPL")
default_tickers_computing <- c("AVGO", "AMD", "MU")
default_tickers_platform <- c("ORCL", "IBM", "PLTR", "CSCO")
default_tickers_adoption <- c("CRM", "NOW", "INTU", "WDAY", "ADBE")

# In this replication file, the disruption-exposed group is the same
# ticker set as the adoption group.
default_tickers_disruption <- default_tickers_adoption

# Non-AI comparison tickers used for the placebo replacement groups.
default_non_ai_tickers <- c("WMT", "COST", "CVX", "ABBV", "PG", "MRK", "PM", "MCD", "PEP", "T")

# Check that all columns needed for a given calculation are available
# in the input dataset. This produces a clear error message if the
# user runs the notebook with an incomplete CSV file.
check_columns_exist <- function(data, cols, obj_name = "data") {
  miss <- setdiff(cols, colnames(data))
  if (length(miss) > 0) {
    stop(
      paste0(
        "Missing column(s) from ", obj_name, ": ",
        paste(miss, collapse = ", ")
      )
    )
  }
}

# Select all intraday observations from a single calendar day in
# Eastern Time. This keeps the event-day selection consistent with
# the timestamps used in the paper.
get_day_data_et <- function(data, date, tz = "America/New_York") {
  idx_date <- format(index(data), tz = tz, format = "%Y-%m-%d")
  data[idx_date == date]
}

# Remove the market component from each stock return series by
# regressing the stock return on the S&P 500 return and keeping the
# residual. The directed network is then estimated from these
# demarketed returns.
remove_market_effect <- function(data_xts) {
  if (!"SP500" %in% colnames(data_xts)) {
    stop("SP500 column is missing from the data.")
  }

  df <- as.data.frame(data_xts)
  market <- df$SP500
  stock_names <- setdiff(colnames(df), "SP500")

  # Run one market-model regression for each ticker and collect the
  # residuals. These residuals are interpreted as firm-specific
  # intraday returns after removing the common market movement.
  resid_mat <- sapply(stock_names, function(col) {
    resid(lm(df[[col]] ~ market))
  })

  resid_mat <- as.matrix(resid_mat)
  colnames(resid_mat) <- stock_names

  xts(resid_mat, order.by = index(data_xts))
}

# Standardize group names so that the main placebo function accepts
# a small number of natural naming variants.
normalize_group_name <- function(group) {
  g <- tolower(group)

  if (g %in% c("m7")) return("m7")
  if (g %in% c("computing", "computation", "compution")) return("computing")
  if (g %in% c("platform")) return("platform")
  if (g %in% c("adoption")) return("adoption")
  if (g %in% c("disruption")) return("disruption")

  stop("Unknown group name. Allowed values: m7, computing, platform, adoption, disruption")
}

# Compute the aggregate effect from one focal node to a target group.
#
# Convention used here:
#   total_table[receiver, transmitter]
#
# Therefore, total_table[targets, node] measures transmission from the
# focal node to the target group, while total_table[node, targets]
# measures transmission from the target group back to the focal node.
# The net effect is outgoing transmission minus incoming transmission.
compute_node_to_targets <- function(total_table, node, targets, mode = c("net", "total")) {
  mode <- match.arg(mode)

  if (!node %in% rownames(total_table) || !node %in% colnames(total_table)) {
    stop("The selected node is not present in total_table.")
  }

  # A node should not be treated as its own target.
  targets <- setdiff(targets, node)

  if (length(targets) == 0) {
    stop("The target group is empty after removing the node itself.")
  }

  miss_r <- setdiff(targets, rownames(total_table))
  miss_c <- setdiff(targets, colnames(total_table))

  if (length(miss_r) > 0 || length(miss_c) > 0) {
    stop("Some target tickers are missing from total_table.")
  }

  # Total outgoing effect from the selected node to the target group.
  out_to_targets <- sum(total_table[targets, node, drop = FALSE])

  if (mode == "total") {
    return(out_to_targets)
  }

  # Net effect subtracts the transmission running in the opposite direction.
  in_from_targets <- sum(total_table[node, targets, drop = FALSE])
  out_to_targets - in_from_targets
}

# Estimate one directed R2 network for a selected intraday return panel.
# If keep_tickers is provided, only those tickers and SP500 are kept
# before the market component is removed.
run_one_network <- function(day_data, keep_tickers = NULL) {
  if (!"SP500" %in% colnames(day_data)) {
    stop("SP500 column is missing from day_data.")
  }

  if (is.null(keep_tickers)) {
    x <- day_data
  } else {
    use_cols <- colnames(day_data)[colnames(day_data) %in% c(keep_tickers, "SP500")]
    check_columns_exist(day_data, c(keep_tickers, "SP500"), "day_data")
    x <- day_data[, use_cols]
  }

  # Demarket the selected ticker set before network estimation.
  x_dem <- remove_market_effect(x)

  # Estimate the directed R2 network. The DAG is learned using LiNGAM,
  # as in the main empirical analysis.
  out <- R2_network(
    data = as.data.frame(x_dem),
    directed = TRUE,
    dag_method = "lingam",
    standardize_for_dag = FALSE
  )

  if (is.null(out$total_table)) {
    stop("The output of R2_network() does not contain total_table.")
  }

  out$total_table
}

# Main placebo routine.
#
# For a selected event day, focal node and exposed group, the function:
#   1. estimates the actual event-day network;
#   2. computes the actual node-to-group effect;
#   3. replaces the exposed group with all non-AI combinations of equal size;
#   4. re-estimates the network for each placebo combination;
#   5. summarizes the actual value relative to the placebo distribution.
placebo_group_test <- function(all_data,
                               day,
                               event_day,
                               node,
                               group,
                               mode = c("net", "total"),
                               tz = "America/New_York",
                               drop_original_group_in_placebo = TRUE,
                               tickers_m7 = NULL,
                               tickers_computing = NULL,
                               tickers_platform = NULL,
                               tickers_adoption = NULL,
                               tickers_disruption = NULL,
                               non_ai_tickers = NULL) {
  mode <- match.arg(mode)

  # Use default ticker groups unless custom groups are supplied.
  if (is.null(tickers_m7)) tickers_m7 <- default_tickers_m7
  if (is.null(tickers_computing)) tickers_computing <- default_tickers_computing
  if (is.null(tickers_platform)) tickers_platform <- default_tickers_platform
  if (is.null(tickers_adoption)) tickers_adoption <- default_tickers_adoption
  if (is.null(tickers_disruption)) tickers_disruption <- default_tickers_disruption
  if (is.null(non_ai_tickers)) non_ai_tickers <- default_non_ai_tickers

  group_name <- normalize_group_name(group)

  group_map <- list(
    m7 = tickers_m7,
    computing = tickers_computing,
    platform = tickers_platform,
    adoption = tickers_adoption,
    disruption = tickers_disruption
  )

  selected_group <- group_map[[group_name]]

  # The AI basket defines the original network universe for the event-day
  # calculation before placebo replacements are introduced.
  ai_basket <- unique(c(
    tickers_m7,
    tickers_computing,
    tickers_platform,
    tickers_adoption,
    tickers_disruption
  ))

  required_cols <- unique(c("SP500", ai_basket, non_ai_tickers))
  check_columns_exist(all_data, required_cols, "all_data")

  # Extract the selected event day.
  day_data <- get_day_data_et(all_data, day, tz = tz)

  if (nrow(day_data) == 0) {
    stop("No data found for the selected day.")
  }

  ai_tickers <- colnames(day_data)[colnames(day_data) %in% ai_basket]

  if (!node %in% ai_tickers) {
    stop("The selected node is not part of the AI basket.")
  }

  # If the focal node is itself part of the selected group, it is removed
  # from the target set. This avoids measuring a self-edge.
  node_in_group <- node %in% selected_group

  real_targets <- if (node_in_group) {
    setdiff(selected_group, node)
  } else {
    selected_group
  }

  n_replace <- length(real_targets)

  if (n_replace < 1) {
    stop("There is nothing to replace in the placebo construction.")
  }

  if (length(non_ai_tickers) < n_replace) {
    stop("There are fewer non-AI tickers than required replacements.")
  }

  # Generate all non-AI replacement groups with the same size as the
  # actual target group.
  combo_list <- combn(non_ai_tickers, n_replace, simplify = FALSE)

  # Actual event-day network and actual node-to-group statistic.
  real_total_table <- run_one_network(
    day_data = day_data,
    keep_tickers = ai_tickers
  )

  real_value <- compute_node_to_targets(
    total_table = real_total_table,
    node = node,
    targets = real_targets,
    mode = mode
  )

  # Define the baseline ticker set used in each placebo network.
  # With drop_original_group_in_placebo = TRUE, the original exposed
  # group is removed before the non-AI replacement tickers are added.
  if (drop_original_group_in_placebo) {
    if (node_in_group) {
      to_remove <- setdiff(selected_group, node)
    } else {
      to_remove <- selected_group
    }

    placebo_base_tickers <- ai_tickers[!(ai_tickers %in% to_remove)]
  } else {
    placebo_base_tickers <- ai_tickers
  }

  # Re-estimate the directed network for each placebo replacement group
  # and store the resulting node-to-group statistic.
  placebo_results <- lapply(seq_along(combo_list), function(i) {
    combo <- combo_list[[i]]
    placebo_tickers <- unique(c(placebo_base_tickers, combo))

    placebo_total_table <- run_one_network(
      day_data = day_data,
      keep_tickers = placebo_tickers
    )

    placebo_value <- compute_node_to_targets(
      total_table = placebo_total_table,
      node = node,
      targets = combo,
      mode = mode
    )

    row <- data.frame(
      combo_id = i,
      day = day,
      event_day = event_day,
      node = node,
      group = group_name,
      mode = mode,
      node_in_group = node_in_group,
      drop_original_group_in_placebo = drop_original_group_in_placebo,
      n_replaced = n_replace,
      placebo_value = placebo_value,
      stringsAsFactors = FALSE
    )

    # Store the identity of each replacement ticker explicitly. This makes
    # the placebo table auditable and helps identify which non-AI baskets
    # generate unusually high or low placebo values.
    for (j in seq_len(n_replace)) {
      row[[paste0("replacement_", j)]] <- combo[j]
    }

    row
  })

  placebo_df <- do.call(rbind, placebo_results)

  # Summary statistics comparing the actual value to the full placebo
  # distribution. The upper-tail p-value is the fraction of placebo values
  # at least as large as the actual event-day value.
  summary_df <- data.frame(
    day = day,
    event_day = event_day,
    node = node,
    group = group_name,
    mode = mode,
    node_in_group = node_in_group,
    drop_original_group_in_placebo = drop_original_group_in_placebo,
    real_targets = paste(real_targets, collapse = ", "),
    real_value = real_value,
    placebo_mean = mean(placebo_df$placebo_value, na.rm = TRUE),
    placebo_median = median(placebo_df$placebo_value, na.rm = TRUE),
    placebo_sd = sd(placebo_df$placebo_value, na.rm = TRUE),
    placebo_min = min(placebo_df$placebo_value, na.rm = TRUE),
    placebo_q05 = as.numeric(quantile(placebo_df$placebo_value, 0.05, na.rm = TRUE)),
    placebo_q25 = as.numeric(quantile(placebo_df$placebo_value, 0.25, na.rm = TRUE)),
    placebo_q50 = as.numeric(quantile(placebo_df$placebo_value, 0.50, na.rm = TRUE)),
    placebo_q75 = as.numeric(quantile(placebo_df$placebo_value, 0.75, na.rm = TRUE)),
    placebo_q95 = as.numeric(quantile(placebo_df$placebo_value, 0.95, na.rm = TRUE)),
    placebo_max = max(placebo_df$placebo_value, na.rm = TRUE),
    percentile = mean(placebo_df$placebo_value <= real_value, na.rm = TRUE),
    p_value_upper = mean(placebo_df$placebo_value >= real_value, na.rm = TRUE),
    n_placebo = nrow(placebo_df),
    stringsAsFactors = FALSE
  )

  # Return both the summary and the underlying objects so that the user can
  # inspect the actual network, the placebo combinations, or the event-day
  # data after running the function.
  list(
    summary = summary_df,
    placebo_df = placebo_df,
    real_value = real_value,
    real_targets = real_targets,
    real_total_table = real_total_table,
    day_data = day_data,
    placebo_base_tickers = placebo_base_tickers,
    all_data = all_data
  )
}

# Copilot preview event.
#
# The actual statistic is the net effect from MSFT to the
# adoption/disruption-exposed group on the event day. The original
# target group is dropped and replaced by equally sized non-AI
# placebo groups.
res1 <- placebo_group_test(
  all_data = all_data,
  day = "2021-06-29",
  event_day = "2021-06-29",
  node = "MSFT",
  group = "disruption",
  mode = "net",
  drop_original_group_in_placebo = TRUE
)

# Hopper/H100 announcement event.
#
# The actual statistic is the net effect from NVDA to the
# computing-exposed group. The placebo distribution replaces the
# computing group with non-AI ticker combinations of the same size.
res2 <- placebo_group_test(
  all_data = all_data,
  day = "2022-03-22",
  event_day = "2022-03-22",
  node = "NVDA",
  group = "computing",
  mode = "net",
  drop_original_group_in_placebo = TRUE
)

# Claude Opus 4.6 event.
#
# The actual statistic is the net effect from CRM to the
# adoption/disruption-exposed group. The placebo distribution is
# built from non-AI replacement groups.
res3 <- placebo_group_test(
  all_data = all_data,
  day = "2026-02-05",
  event_day = "2026-02-05",
  node = "CRM",
  group = "disruption",
  mode = "net",
  drop_original_group_in_placebo = TRUE
)

# ChatGPT release event.
#
# This specification measures the net effect from MSFT to the
# adoption/disruption-exposed group on the ChatGPT release day.
res4 <- placebo_group_test(
  all_data = all_data,
  day = "2022-11-30",
  event_day = "2022-11-30",
  node = "MSFT",
  group = "disruption",
  mode = "net",
  drop_original_group_in_placebo = TRUE
)

# ChatGPT release event.
#
# This specification measures the net effect from NVDA to the
# computing-exposed group on the ChatGPT release day.
res5 <- placebo_group_test(
  all_data = all_data,
  day = "2022-11-30",
  event_day = "2022-11-30",
  node = "NVDA",
  group = "computing",
  mode = "net",
  drop_original_group_in_placebo = TRUE
)

library(ggplot2)

# Plot the placebo distribution and mark the actual event-day value.
#
# Inputs:
#   placebo_df: data frame returned by placebo_group_test()$placebo_df
#   real_value: actual event-day node-to-group effect
#   bins: number of histogram bins
#   save_path: optional file path for saving the figure
plot_placebo_distribution <- function(placebo_df,
                                      real_value = attr(placebo_df, "real_value"),
                                      event_name = NULL,
                                      event_day = NULL,
                                      bins = 10,
                                      xlab = "Placebo NET value",
                                      ylab = "Relative frequency",
                                      title_text = NULL,
                                      subtitle_text = NULL,
                                      title_prefix = "Placebo Distribution of the",
                                      net_label = "Net Effect",
                                      total_label = "Total Effect",
                                      from_prefix = "From",
                                      to_prefix = "",
                                      event_prefix = "",
                                      date_prefix = "",
                                      distribution_note = "",
                                      subtitle_sep = "",
                                      original_effect_prefix = "Actual edge =",
                                      digits = 4,
                                      unknown_node = "Unknown node",
                                      unknown_group = "Unknown group",
                                      unknown_event = "Unknown event",
                                      save_path = NULL,
                                      width = 8,
                                      height = 5,
                                      dpi = 300) {

  # Basic input checks. These make the plotting function fail with a
  # clear message if it is called with an incomplete placebo object.
  if (!"placebo_value" %in% names(placebo_df)) {
    stop("The placebo_df object does not contain a 'placebo_value' column.")
  }

  if (is.null(real_value) || length(real_value) != 1 || !is.finite(real_value)) {
    stop("The real_value is missing. Provide it explicitly or store it as an attribute.")
  }

  # Keep only finite placebo values before constructing the histogram.
  x <- placebo_df$placebo_value
  x <- x[is.finite(x)]

  if (length(x) == 0) {
    stop("There are no valid placebo_value observations to plot.")
  }

  # Use metadata from the placebo table when it is available. These labels
  # can also be overwritten manually through the function arguments.
  node <- if ("node" %in% names(placebo_df)) unique(placebo_df$node)[1] else unknown_node
  group <- if ("group" %in% names(placebo_df)) unique(placebo_df$group)[1] else unknown_group
  mode <- if ("mode" %in% names(placebo_df)) unique(placebo_df$mode)[1] else "effect"

  if (is.null(event_day) && "event_day" %in% names(placebo_df)) {
    event_day <- unique(placebo_df$event_day)[1]
  }

  if (is.null(event_name)) {
    event_name <- unknown_event
  }

  mode_label <- if (tolower(mode) == "net") net_label else total_label

  if (is.null(title_text)) {
    title_text <- paste(title_prefix, mode_label)
  }

  if (is.null(subtitle_text)) {
    subtitle_parts <- c(
      paste(from_prefix, node, to_prefix, group),
      paste(event_prefix, event_name),
      if (!is.null(event_day)) paste(date_prefix, event_day) else NULL,
      distribution_note
    )
    subtitle_text <- paste(subtitle_parts, collapse = subtitle_sep)
  }

  # Compute histogram bins manually so that the y-axis can be expressed
  # as relative frequency instead of raw counts.
  h <- hist(x, breaks = bins, plot = FALSE)

  hist_df <- data.frame(
    xmin = head(h$breaks, -1),
    xmax = tail(h$breaks, -1),
    xmid = h$mids,
    count = h$counts
  )

  hist_df$rel_freq <- hist_df$count / sum(hist_df$count)

  ymax <- max(hist_df$rel_freq, na.rm = TRUE)

  # Place the label on the side of the vertical line where it is less
  # likely to overlap with the plot boundary.
  x_center <- mean(range(c(x, real_value), na.rm = TRUE))
  hjust_val <- if (real_value <= x_center) -0.05 else 1.05

  line_label <- paste0(original_effect_prefix, " ", round(real_value, digits))

  # Histogram of placebo values with the actual event-day effect shown
  # as a vertical reference line.
  p <- ggplot(hist_df, aes(x = xmid, y = rel_freq)) +
    geom_col(
      width = diff(h$breaks)[1],
      fill = "grey70",
      color = "black",
      linewidth = 0.45
    ) +
    geom_segment(
      aes(
        x = real_value,
        xend = real_value,
        y = 0,
        yend = ymax
      ),
      color = "red",
      linewidth = 1.35,
      inherit.aes = FALSE
    ) +
    annotate(
      "text",
      x = real_value,
      y = ymax * 0.95,
      label = line_label,
      color = "red",
      hjust = hjust_val,
      vjust = 1,
      size = 5,
      fontface = "bold"
    ) +
    labs(
      title = "",
      subtitle = "",
      x = xlab,
      y = ylab
    ) +
    theme_minimal(base_size = 15) +
    theme(
      axis.title = element_text(face = "bold", size = 16),
      axis.text = element_text(size = 14, face = "bold", color = "black"),
      axis.line = element_line(color = "black", linewidth = 0.8),
      axis.ticks = element_line(color = "black", linewidth = 0.8),
      axis.ticks.length = grid::unit(0.18, "cm"),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(linewidth = 0.35),
      plot.title = element_text(face = "bold", size = 16),
      plot.subtitle = element_text(size = 12)
    )

  # Save the figure only when an output path is supplied.
  if (!is.null(save_path)) {
    ggsave(
      filename = save_path,
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }

  return(p)
}

# Suppress non-critical warnings during batch figure generation.
options(warn = -1)

# Create and save Figure E.3 for the Copilot preview placebo exercise.
plot_placebo_distribution(
  placebo_df = res1$placebo_df,
  real_value = res1$real_value,
  event_name = "",
  event_day = "",
  xlab = "NET value",
  ylab = "Relative frequency",
  title_prefix = "",
  net_label = "",
  total_label = "",
  #from_prefix = "From MSFT to disruption",
  bins=10,
  #to_prefix = "Receiver group:",
  #event_prefix = "Event",
  #date_prefix = "Event day",
  #distribution_note = "Non-AI placebo distribution",
  original_effect_prefix = "Actual edge = ",
   save_path = "[YOUR FOLDER]/placebo_net_distribution1.png"
)

# Create and save Figure E.4 for the Hopper/H100 placebo exercise.
plot_placebo_distribution(
  placebo_df = res2$placebo_df,
  real_value = res2$real_value,
  event_name = "",
  event_day = "",
  xlab = "NET value",
  ylab = "Relative frequency",
  title_prefix = "",
  net_label = "",
  total_label = "",
  #from_prefix = "From MSFT to disruption",
  bins=10,
  #to_prefix = "Receiver group:",
  #event_prefix = "Event",
  #date_prefix = "Event day",
  #distribution_note = "Non-AI placebo distribution",
  original_effect_prefix = "Actual edge = ",
   save_path = "[YOUR FOLDER]/placebo_net_distribution2.png"
)

# Create and save Figure E.5 for the Opus 4.6 placebo exercise.
plot_placebo_distribution(
  placebo_df = res3$placebo_df,
  real_value = res3$real_value,
  event_name = "",
  event_day = "",
  xlab = "NET value",
  ylab = "Relative frequency",
  title_prefix = "",
  net_label = "",
  total_label = "",
  #from_prefix = "From MSFT to disruption",
  bins=10,
  #to_prefix = "Receiver group:",
  #event_prefix = "Event",
  #date_prefix = "Event day",
  #distribution_note = "Non-AI placebo distribution",
  original_effect_prefix = "Actual edge = ",
   save_path = "[YOUR FOLDER]/placebo_net_distribution3.png"
)

# Create and save Figure E.1 for the ChatGPT/MSFT placebo exercise.
plot_placebo_distribution(
  placebo_df = res4$placebo_df,
  real_value = res4$real_value,
  event_name = "",
  event_day = "",
  xlab = "NET value",
  ylab = "Relative frequency",
  title_prefix = "",
  net_label = "",
  total_label = "",
  #from_prefix = "From MSFT to disruption",
  bins=10,
  #to_prefix = "Receiver group:",
  #event_prefix = "Event",
  #date_prefix = "Event day",
  #distribution_note = "Non-AI placebo distribution",
  original_effect_prefix = "Actual edge = ",
   save_path = "[YOUR FOLDER]/placebo_net_distribution4.png"
)

# Create and save Figure E.2 for the ChatGPT/NVDA placebo exercise.
plot_placebo_distribution(
  placebo_df = res5$placebo_df,
  real_value = res5$real_value,
  event_name = "",
  event_day = "",
  xlab = "NET value",
  ylab = "Relative frequency",
  title_prefix = "",
  net_label = "",
  total_label = "",
  #from_prefix = "From MSFT to disruption",
  bins=10,
  #to_prefix = "Receiver group:",
  #event_prefix = "Event",
  #date_prefix = "Event day",
  #distribution_note = "Non-AI placebo distribution",
  original_effect_prefix = "Actual edge = ",
   save_path = "[YOUR FOLDER]/placebo_net_distribution5.png"
)