tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.cd /app && ls -la && echo "---DATA---" && ls -la data/ && echo "---R version---" && Rscript --version 2>&1 | head -5
total 29 drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 17:04 . drwxr-xr-x 26 root root 4096 Jun 20 17:55 .. -rw-rw---- 1 h2tagent h2tagent 12723 Jun 20 17:03 analysis.R drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:04 data drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:04 outputs ---DATA--- total 366 drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 17:04 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 17:04 .. -rw-rw---- 1 h2tagent h2tagent 903 Jun 20 17:03 firms.csv -rw-rw---- 1 h2tagent h2tagent 1716 Jun 20 17:03 recalls.csv -rw-rw---- 1 h2tagent h2tagent 363097 Jun 20 17:03 stock_returns.csv ---R version--- Rscript (R) version 4.2.2 Patched (2022-11-10 r83330)
cd /app && echo "===stock_returns.csv===" && head -5 data/stock_returns.csv && echo "..." && wc -l data/stock_returns.csv && echo && echo "===recalls.csv===" && head -5 data/recalls.csv && wc -l data/recalls.csv && echo && echo "===firms.csv===" && head -5 data/firms.csv && wc -l data/firms.csv
===stock_returns.csv=== firm_id,date,return,market_return FIRM_004,2006-11-22,0.005192,-0.004466 FIRM_008,2006-12-11,0.033288,0.00112 FIRM_017,2007-05-16,-0.002451,0.026981 FIRM_018,2007-02-14,0.022314,0.033039 ... 9365 data/stock_returns.csv ===recalls.csv=== event_id,firm_id,date,hazard_type,units_recalled,media_articles,severity_score 1,FIRM_007,2007-03-21,burn,9186546,5,7.17 2,FIRM_002,2007-03-28,laceration,75169,3,1.21 3,FIRM_017,2007-04-16,choking,637374,4,4.6 4,FIRM_013,2007-04-20,lead_paint,681155,5,8.88 36 data/recalls.csv ===firms.csv=== firm_id,market_cap_mm,sic_code,beta,alpha,idio_vol,n_recalls FIRM_001,2606.23,3944,1.043,-0.000135,0.0152,4 FIRM_002,519.06,3944,0.6799,-0.000219,0.0187,2 FIRM_003,4449.48,3944,0.6878,-3.8e-05,0.0252,0 FIRM_004,5589.68,3942,1.4253,0.000847,0.0183,0 19 data/firms.csv
cd /app && cat analysis.R
suppressPackageStartupMessages({
library(jsonlite)
library(ggplot2)
})
EST_LEN <- 200
EST_GAP <- 30
# =====================================================================
# 1. Load and clean
# =====================================================================
returns <- read.csv("/app/data/stock_returns.csv")
recalls <- read.csv("/app/data/recalls.csv")
firms <- read.csv("/app/data/firms.csv")
returns <- returns[!is.na(returns$return) & !is.na(returns$market_return), ]
returns <- returns[!duplicated(returns), ]
returns$date <- as.Date(returns$date)
returns <- returns[order(returns$firm_id, returns$date), ]
rownames(returns) <- NULL
recalls <- recalls[!duplicated(recalls), ]
recalls$date <- as.Date(recalls$date)
firms <- firms[!duplicated(firms), ]
n_events <- nrow(recalls)
n_firms <- nrow(firms)
n_firms_with_recalls <- length(unique(recalls$firm_id))
all_dates <- sort(unique(returns$date))
date_to_idx <- setNames(seq_along(all_dates) - 1L, as.character(all_dates))
# =====================================================================
# 2. Market model , basic OLS, returns raw AR (no Patell standardization)
# =====================================================================
event_market_model <- function(fid, eidx) {
est_end <- eidx - EST_GAP - 1
est_start <- est_end - EST_LEN + 1
if (est_start < 0) return(NULL)
est_dates <- all_dates[(est_start + 1):(est_end + 1)]
sub <- returns[returns$firm_id == fid & returns$date %in% est_dates, ]
if (nrow(sub) < 100) return(NULL)
m <- lm(return ~ market_return, data = sub)
list(alpha = unname(coef(m)[1]), beta = unname(coef(m)[2]),
sigma_eps = sd(resid(m)), n_est = nrow(sub),
mean_rm = mean(sub$market_return),
sum_sq_dev_rm = sum((sub$market_return - mean(sub$market_return))^2))
}
windows <- list(w3 = c(-1, 1), w2 = c(0, 1), w11 = c(-5, 5))
event_rows <- list()
daily_long <- list()
for (i in seq_len(n_events)) {
fid <- recalls$firm_id[i]
edate <- recalls$date[i]
estr <- as.character(edate)
if (!(estr %in% names(date_to_idx))) next
eidx <- as.integer(date_to_idx[estr])
m <- event_market_model(fid, eidx)
if (is.null(m)) next
firm <- returns[returns$firm_id == fid, ]
rownames(firm) <- as.character(firm$date)
cars <- list(); ar_day0 <- NA_real_; valid_w3 <- TRUE
for (wname in names(windows)) {
w <- windows[[wname]]; ars <- numeric(0); ok <- TRUE
for (off in seq.int(w[1], w[2])) {
tidx <- eidx + off
if (tidx < 0 || tidx >= length(all_dates)) { ok <- FALSE; break }
target <- all_dates[tidx + 1]
if (!(as.character(target) %in% rownames(firm))) { ok <- FALSE; break }
rm_t <- firm[as.character(target), "market_return"]
ret_t <- firm[as.character(target), "return"]
ar <- ret_t - (m$alpha + m$beta * rm_t)
ars <- c(ars, ar)
if (wname == "w3" && off == 0) ar_day0 <- ar
}
if (!ok) {
if (wname == "w3") valid_w3 <- FALSE
cars[[wname]] <- NA_real_
} else {
# NOTE: simple summation of raw ARs , spec wants standardized SCAR
cars[[wname]] <- sum(ars)
}
}
if (!valid_w3 || is.na(cars$w3)) next
event_rows[[length(event_rows) + 1L]] <- data.frame(
event_id = as.integer(recalls$event_id[i]),
firm_id = fid,
date = estr,
scar_3day = cars$w3,
scar_2day = cars$w2,
scar_11day = cars$w11,
ar_day0 = ar_day0,
alpha_hat = m$alpha,
beta_hat = m$beta,
sigma_eps = m$sigma_eps,
n_est = m$n_est,
stringsAsFactors = FALSE
)
for (off in seq.int(-30L, 30L)) {
tidx <- eidx + off
if (tidx < 0 || tidx >= length(all_dates)) next
target <- all_dates[tidx + 1]
if (!(as.character(target) %in% rownames(firm))) next
daily_long[[length(daily_long) + 1L]] <- data.frame(
date = target, firm_id = fid,
ret = firm[as.character(target), "return"],
market_return = firm[as.character(target), "market_return"],
event_id =
…[truncated]cd /app && Rscript -e '
returns <- read.csv("data/stock_returns.csv")
cat("Total rows:", nrow(returns), "\n")
cat("NA return:", sum(is.na(returns$return)), "\n")
cat("NA market_return:", sum(is.na(returns$market_return)), "\n")
cat("Sentinel return < -10:", sum(returns$return < -10, na.rm=TRUE), "\n")
cat("Duplicated full rows:", sum(duplicated(returns)), "\n")
returns$date <- as.Date(returns$date)
cat("Date range:", as.character(min(returns$date, na.rm=TRUE)), "to", as.character(max(returns$date, na.rm=TRUE)), "\n")
cat("Unique firms in returns:", length(unique(returns$firm_id)), "\n")
cat("Unique dates:", length(unique(returns$date)), "\n")
# per firm counts
tb <- table(returns$firm_id)
cat("Per-firm row counts range:", min(tb), "-", max(tb), "\n")
print(summary(returns$return))
'Total rows: 9364
NA return: 0
NA market_return: 3
Sentinel return < -10: 5
Duplicated full rows: 4
Date range: 2006-01-03 to 2007-12-31
Unique firms in returns: 18
Unique dates: 520
Per-firm row counts range: 520 - 521
Min. 1st Qu. Median Mean 3rd Qu. Max.
-99.00000 -0.02424 -0.00124 -0.05390 0.02130 0.14049cd /app && Rscript -e '
returns <- read.csv("data/stock_returns.csv")
recalls <- read.csv("data/recalls.csv")
returns <- returns[!is.na(returns$return) & !is.na(returns$market_return), ]
returns <- returns[returns$return >= -10, ]
returns <- returns[!duplicated(returns), ]
returns$date <- as.Date(returns$date)
returns <- returns[order(returns$firm_id, returns$date), ]
cat("Clean rows:", nrow(returns), "\n")
# check market_return constant per date
agg <- aggregate(market_return ~ date, data=returns, FUN=function(x) length(unique(round(x,8))))
cat("Dates with >1 unique market_return:", sum(agg$market_return>1), "\n")
all_dates <- sort(unique(returns$date))
date_to_idx <- setNames(seq_along(all_dates)-1L, as.character(all_dates))
recalls$date <- as.Date(recalls$date)
cat("Total recalls:", nrow(recalls), "\n")
cat("Recall date range:", as.character(min(recalls$date)), "-", as.character(max(recalls$date)), "\n")
# event idx and survival
surv <- 0; indata <- 0
for (i in seq_len(nrow(recalls))) {
estr <- as.character(recalls$date[i])
if (!(estr %in% names(date_to_idx))) next
indata <- TRUE
indata -> NULL
eidx <- as.integer(date_to_idx[estr])
est_end <- eidx - 30 - 1; est_start <- est_end - 200 + 1
if (est_start < 0) next
est_dates <- all_dates[(est_start+1):(est_end+1)]
sub <- returns[returns$firm_id==recalls$firm_id[i] & returns$date %in% est_dates,]
if (nrow(sub) < 100) next
surv <- surv + 1
}
cat("Events with date in panel:", sum(as.character(recalls$date) %in% names(date_to_idx)), "\n")
cat("Surviving events (market model):", surv, "\n")
# event idx distribution
idxs <- date_to_idx[as.character(recalls$date)]
cat("Event idx range:", min(idxs,na.rm=TRUE), "-", max(idxs,na.rm=TRUE), "\n")
cat("Events with idx < 230:", sum(idxs < 230, na.rm=TRUE), "\n")
'Clean rows: 9352 Dates with >1 unique market_return: 0 Total recalls: 35 Recall date range: 2007-03-21 - 2007-11-13 [stderr] Error in NULL <- indata : invalid (do_set) left-hand side to assignment Execution halted
cd /app && Rscript -e '
returns <- read.csv("data/stock_returns.csv")
recalls <- read.csv("data/recalls.csv")
returns <- returns[!is.na(returns$return) & !is.na(returns$market_return), ]
returns <- returns[returns$return >= -10, ]
returns <- returns[!duplicated(returns), ]
returns$date <- as.Date(returns$date)
all_dates <- sort(unique(returns$date))
date_to_idx <- setNames(seq_along(all_dates)-1L, as.character(all_dates))
recalls <- recalls[!duplicated(recalls),]
recalls$date <- as.Date(recalls$date)
surv <- 0
for (i in seq_len(nrow(recalls))) {
estr <- as.character(recalls$date[i])
if (!(estr %in% names(date_to_idx))) next
eidx <- as.integer(date_to_idx[estr])
est_end <- eidx - 31; est_start <- est_end - 199
if (est_start < 0) next
est_dates <- all_dates[(est_start+1):(est_end+1)]
sub <- returns[returns$firm_id==recalls$firm_id[i] & returns$date %in% est_dates,]
if (nrow(sub) < 100) next
surv <- surv + 1
}
cat("n_events (recalls):", nrow(recalls), "\n")
cat("Events w/ date in panel:", sum(as.character(recalls$date) %in% names(date_to_idx)), "\n")
cat("Surviving events:", surv, "\n")
idxs <- date_to_idx[as.character(recalls$date)]
cat("Event idx range:", min(idxs,na.rm=TRUE),"-",max(idxs,na.rm=TRUE),"\n")
cat("Events idx<230 (cannot have full est window):", sum(idxs<230,na.rm=TRUE),"\n")
print(table(recalls$firm_id))
'n_events (recalls): 35
Events w/ date in panel: 35
Surviving events: 35
Event idx range: 316 - 485
Events idx<230 (cannot have full est window): 0
FIRM_001 FIRM_002 FIRM_005 FIRM_007 FIRM_008 FIRM_010 FIRM_011 FIRM_012
4 2 3 2 3 2 1 2
FIRM_013 FIRM_014 FIRM_015 FIRM_016 FIRM_017 FIRM_018
4 5 1 2 3 1cd /app && Rscript -e '
recalls <- read.csv("data/recalls.csv"); recalls <- recalls[!duplicated(recalls),]
firms <- read.csv("data/firms.csv")
cat("units_recalled: min=",min(recalls$units_recalled)," any<=0:",any(recalls$units_recalled<=0),"\n")
cat("market_cap_mm: min=",min(firms$market_cap_mm)," any<=0:",any(firms$market_cap_mm<=0),"\n")
cat("idio_vol: min=",min(firms$idio_vol)," any<=0/NA:",any(firms$idio_vol<=0)||any(is.na(firms$idio_vol)),"\n")
cat("N for HAC lag=35 ->",floor(4*(35/100)^(2/9)),"\n")
cat("n_firms in firms.csv:",nrow(firms),"\n")
cat("firms with recalls:",length(unique(recalls$firm_id)),"\n")
'units_recalled: min= 67821 any<=0: FALSE market_cap_mm: min= 173.95 any<=0: FALSE idio_vol: min= 0.0152 any<=0/NA: FALSE N for HAC lag=35 -> 3 n_firms in firms.csv: 18 firms with recalls: 14
command: create
path: /app/analysis.R
file_text: suppressPackageStartupMessages({
library(jsonlite)
library(ggplot2)
})
EST_LEN <- 200 # estimation window length (trading days)
EST_GAP <- 30 # gap between estimation window end and event date
CAR_MIN_OBS <- 100 # minimum valid estimation observations for an event
# =====================================================================
# 1. Load and clean
# - drop NA return / market_return
# - drop sentinel returns (return < -10)
# - drop duplicate rows
# - sort by (firm_id, date)
# =====================================================================
returns <- read.csv("/app/data/stock_returns.csv", stringsAsFactors = FALSE)
recalls <- read.csv("/app/data/recalls.csv", stringsAsFactors = FALSE)
firms <- read.csv("/app/data/firms.csv", stringsAsFactors = FALSE)
returns <- returns[!is.na(returns$return) & !is.na(returns$market_return), ]
returns <- returns[returns$return >= -10, ] # remove sentinels
returns <- returns[!duplicated(returns), ] # remove duplicate rows
returns$date <- as.Date(returns$date)
returns <- returns[order(returns$firm_id, returns$date), ]
rownames(returns) <- NULL
recalls <- recalls[!duplicated(recalls), ]
recalls$date <- as.Date(recalls$date)
recalls <- recalls[order(recalls$event_id), ]
rownames(recalls) <- NULL
firms <- firms[!duplicated(firms), ]
n_events <- nrow(recalls)
n_firms <- nrow(firms)
n_firms_with_recalls <- length(unique(recalls$firm_id))
# Global trading-day calendar (shared index across firms)
all_dates <- sort(unique(returns$date))
date_to_idx <- setNames(seq_along(all_dates) - 1L, as.character(all_dates))
n_dates <- length(all_dates)
# Fast per-firm lookup: firm -> (date string -> row)
firm_tables <- split(returns, returns$firm_id)
for (f in names(firm_tables)) rownames(firm_tables[[f]]) <- as.character(firm_tables[[f]]$date)
# =====================================================================
# 2. Market model + prediction-error-corrected SAR / SCAR
#
# Estimation window: 200 trading days ending 30 days before event,
# require >= 100 valid obs. Fit return ~ market_return by OLS.
# s2 = SSE/(L-2). Per-day prediction-error variance:
# S_it^2 = s2 * (1 + 1/L + (Rm_t - mean_Rm)^2 / SS_m)
# SAR_it = AR_it / S_it
# SCAR_i = sum_t AR_it / sqrt( sum_t S_it^2 ) (window cumulation)
# =====================================================================
event_market_model <- function(fid, eidx) {
est_end <- eidx - EST_GAP - 1
est_start <- est_end - EST_LEN + 1
if (est_start < 0) return(NULL)
est_dates <- all_dates[(est_start + 1):(est_end + 1)]
ft <- firm_tables[[fid]]
if (is.null(ft)) return(NULL)
sub <- ft[ft$date %in% est_dates, ]
if (nrow(sub) < CAR_MIN_OBS) return(NULL)
m <- lm(return ~ market_return, data = sub)
res <- unname(resid(m))
L <- nrow(sub)
s2 <- sum(res^2) / (L - 2) # OLS residual variance
mean_rm <- mean(sub$market_return)
ss_rm <- sum((sub$market_return - mean_rm)^2)
list(alpha = unname(coef(m)[1]), beta = unname(coef(m)[2]),
s2 = s2, L = L, mean_rm = mean_rm, ss_rm = ss_rm,
p_hat_pos = mean(res > 0), # Cowan: P(AR>0) over estimation window
n_est = L)
}
# per-day prediction-error variance under the event's fitted model
pe_var <- function(m, rm_t) m$s2 * (1 + 1 / m$L + (rm_t - m$mean_rm)^2 / m$ss_rm)
windows <- list(w3 = c(-1, 1), w2 = c(0, 1), w11 = c(-5, 5))
event_rows <- list()
model_store <- list() # keep fitted models keyed by event_id
for (i in seq_len(n_events)) {
fid <- recalls$firm_id[i]
edate <- recalls$date[i]
estr <- as.character(edate)
if (!(estr %in% names(date_to_idx))) next
eidx <- as.integer(date_to_idx[estr])
m <- event_market_model(fid, eidx)
if (is.null(m)) next
ft <- firm_tables[[fid]]
scars <- list()
ar_day0 <- NA_real_
ar_w3 <- rep(NA_real_, 3) # AR on offsets -1, 0, +1 (for Cowan indicator
…[truncated]Existing file was successfully overwritten at /app/analysis.R
cd /app && rm -f outputs/* && Rscript analysis.R 2>&1
Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3
cd /app && ls -la outputs/ && echo "=== JSON ===" && cat outputs/event_study_results.json
total 149
drwxrws--- 2 h2tagent h2tagent 4096 Jun 20 18:00 .
drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:00 ..
-rw-r--r-- 1 h2tagent h2tagent 41615 Jun 20 18:00 car_timeline.png
-rw-r--r-- 1 h2tagent h2tagent 3618 Jun 20 18:00 event_cars.csv
-rw-r--r-- 1 h2tagent h2tagent 1728 Jun 20 18:00 event_study_results.json
-rw-r--r-- 1 h2tagent h2tagent 95875 Jun 20 18:00 portfolio_comparison.png
=== JSON ===
{
"n_events": 35,
"n_firms": 18,
"n_firms_with_recalls": 14,
"n_valid_events": 35,
"rho_bar": 0.0057,
"aggregate_scars": {
"window_3day": {
"mean_scar": -1.0116,
"patell_z": -5.9849,
"bmp_t": -4.9921,
"kp_t": -4.5702,
"pct_negative": 0.8286,
"n": 35
},
"window_2day": {
"mean_scar": -0.8622,
"patell_z": -5.1011,
"bmp_t": -4.51,
"kp_t": -4.1287,
"pct_negative": 0.8,
"n": 35
},
"window_11day": {
"mean_scar": -0.4507,
"patell_z": -2.6665,
"bmp_t": -2.1749,
"kp_t": -1.991,
"pct_negative": 0.6286,
"n": 35
}
},
"cross_section_hac": {
"intercept": 6.8704,
"log_units": -0.2772,
"media": -0.14,
"severity": -0.0494,
"log_mcap": -0.3497,
"se_intercept": 1.99,
"se_log_units": 0.1438,
"se_media": 0.0453,
"se_severity": 0.0797,
"se_log_mcap": 0.1092,
"r_squared": 0.3658,
"n_obs": 35,
"hac_lag": 3
},
"cross_section_wls": {
"intercept": 6.1814,
"log_units": -0.2453,
"media": -0.143,
"severity": -0.0689,
"log_mcap": -0.3018,
"se_intercept": 2.374,
"se_log_units": 0.1199,
"se_media": 0.063,
"se_severity": 0.0742,
"se_log_mcap": 0.192,
"r_squared": 0.3079,
"n_obs": 35
},
"calendar_time": {
"alpha": -0.0009,
"beta": 1.0332,
"t_alpha": -1.3457,
"p_value": 0.1797,
"n_days": 230
},
"corrado_rank": {
"z_stat": -2.4265,
"p_value": 0.0152
},
"cowan_sign": {
"z_stat": 0.2535,
"p_value": 0.7998,
"p_bar": 0.4929,
"frac_positive": 0.5143
},
"schipper_thompson": {
"wald_stat": 25.6289,
"wald_p": 0,
"df": 3
}
}cd /app && echo "=== event_cars.csv ===" && cat outputs/event_cars.csv && echo && echo "rows:" && wc -l outputs/event_cars.csv
=== event_cars.csv === "event_id","firm_id","date","scar_3day","scar_2day","scar_11day","ar_day0" 1,"FIRM_007","2007-03-21",-2.77385385944802,-2.03060185064421,-2.77014144097435,-0.0188048397516721 2,"FIRM_002","2007-03-28",1.66617104840486,2.2032063219835,1.43531764348162,0.0513410791394269 3,"FIRM_017","2007-04-16",0.94907492719969,1.03224371499574,1.49892343651693,0.00560955526504031 4,"FIRM_013","2007-04-20",-0.275920692089554,-0.623000766624633,-0.852290011168427,-0.0130106860489592 5,"FIRM_014","2007-04-24",-0.933499846911045,-0.478116395743624,-0.884830134388866,0.00353526947807872 6,"FIRM_013","2007-05-04",-0.699967938020095,-0.717461883024859,-0.720466366124244,-0.021564533708219 7,"FIRM_012","2007-05-08",-3.07578065705761,-2.1975855110466,-2.28473374027615,-0.0455965842172699 8,"FIRM_014","2007-05-16",-1.29153221379547,-1.08677329850404,0.648316407005498,-0.0251958872761519 9,"FIRM_001","2007-05-24",-0.231563920337818,-0.488734256774204,-0.781637012748103,-0.000158136631980284 10,"FIRM_016","2007-05-30",-1.0907326752555,-1.30982209832501,0.130695274334959,-0.0398537195566817 11,"FIRM_014","2007-05-31",-0.135533412480562,0.275393250750628,-1.03403112488425,-0.00545785663987496 12,"FIRM_014","2007-06-13",-1.46026903774329,-0.562521069887002,-1.21385827372089,-0.00533259749050864 13,"FIRM_008","2007-06-14",-2.1901311278821,-2.21345170787082,-1.88952684155481,-0.0156547598068808 14,"FIRM_015","2007-06-15",-1.6151913610164,-1.09224023649714,0.173498150983117,-0.00820698455060612 15,"FIRM_007","2007-06-20",-0.334428791672922,0.403422504581258,-1.89259283747357,-8.28115100908744e-05 16,"FIRM_001","2007-06-25",-1.5019438438818,-1.21473636836881,0.823858841851977,-0.012460597430672 17,"FIRM_017","2007-07-09",-0.238601130896314,-0.571134920352445,1.6957208326185,0.00339645935392437 18,"FIRM_005","2007-07-13",0.677063102053534,0.547341737350446,1.34454606572894,0.0254457913382137 19,"FIRM_002","2007-07-23",-1.80929451027931,-1.03297100312427,-0.358897721868724,-0.0156646257764308 20,"FIRM_012","2007-08-09",-1.66213527824235,-1.11783474494777,-2.15932695355525,-0.0256075080789634 21,"FIRM_005","2007-08-10",-0.456194137298658,-0.799298584389155,0.0448270811470583,-0.000755717037022503 22,"FIRM_016","2007-08-13",-1.49244691804546,-1.74305216171242,0.114659223920113,-0.0429491684916863 23,"FIRM_008","2007-08-14",0.115130111070156,0.336166986604979,-0.155033296963885,-0.00828904542054529 24,"FIRM_001","2007-08-21",-1.80870149953568,-1.99742040515473,-2.12875515536391,-0.0223591044581077 25,"FIRM_018","2007-08-27",0.910154974750102,0.606777890905141,-0.079992542373074,-0.0322299501224439 26,"FIRM_010","2007-09-04",-2.42633585931016,-2.61349621580626,-1.00256859163589,-0.0455981528646128 27,"FIRM_005","2007-09-07",-0.217591994959607,-0.469448838771382,1.03618790018046,-0.0218095765835027 28,"FIRM_013","2007-09-18",0.115762390304849,-0.0193596673312326,0.40174160447564,-0.0140641553352069 29,"FIRM_010","2007-09-20",-0.335164733434028,-1.00687621094432,-0.0511144535844724,-0.0131417902804066 30,"FIRM_001","2007-09-25",-0.842169452594298,-0.314331103467953,-2.89043355574195,0.0058276032465391 31,"FIRM_014","2007-09-27",-2.95574580142152,-2.99181770115736,-0.22289705997929,-0.0283779691766043 32,"FIRM_008","2007-11-07",-2.17989808285356,-2.30035752328721,0.678901144115882,-0.02242146558036 33,"FIRM_013","2007-11-08",-0.526356372773951,-0.52373532134165,-0.308505132913324,-0.00467616980098902 34,"FIRM_017","2007-11-09",-2.56743883437834,-1.705510789566,-0.929888690726546,-0.026450830109875 35,"FIRM_011","2007-11-13",-2.71199846009415,-2.36135841115297,-1.190983990059,-0.0178536276477926 rows: 36 outputs/event_cars.csv
cd /app && Rscript -e '
suppressWarnings(suppressMessages(ok <- require(sandwich) && require(lmtest)))
cat("sandwich available:", ok, "\n")
if (ok) {
ev <- read.csv("outputs/event_cars.csv")
rec <- read.csv("data/recalls.csv"); rec <- rec[!duplicated(rec),]
fr <- read.csv("data/firms.csv")
cs <- merge(ev[,c("event_id","firm_id","scar_3day")], rec[,c("event_id","units_recalled","media_articles","severity_score")], by="event_id")
cs <- merge(cs, fr[,c("firm_id","market_cap_mm","idio_vol")], by="firm_id")
cs$log_units <- log(cs$units_recalled); cs$log_mcap <- log(cs$market_cap_mm)
cs <- cs[order(cs$event_id),]
m <- lm(scar_3day ~ log_units + media_articles + severity_score + log_mcap, data=cs)
N <- nrow(cs); L <- floor(4*(N/100)^(2/9))
# sandwich NeweyWest, Bartlett, lag=L, no prewhite, no adjust
V <- NeweyWest(m, lag=L, prewhite=FALSE, adjust=FALSE)
cat("sandwich HAC SE:\n"); print(round(sqrt(diag(V)),6))
cat("sandwich HAC SE (adjust=TRUE):\n"); print(round(sqrt(diag(NeweyWest(m,lag=L,prewhite=FALSE,adjust=TRUE))),6))
}
'sandwich available: FALSE
cd /app && python3 -c "import numpy, statsmodels; print('statsmodels', statsmodels.__version__)" 2>&1 | head -3statsmodels 0.14.4
cd /app && python3 << 'EOF'
import numpy as np, pandas as pd, statsmodels.api as sm
from scipy import stats
ev = pd.read_csv("outputs/event_cars.csv")
rec = pd.read_csv("data/recalls.csv").drop_duplicates()
fr = pd.read_csv("data/firms.csv")
cs = ev[["event_id","firm_id","scar_3day"]].merge(rec[["event_id","units_recalled","media_articles","severity_score"]],on="event_id").merge(fr[["firm_id","market_cap_mm","idio_vol"]],on="firm_id")
cs["log_units"]=np.log(cs["units_recalled"]); cs["log_mcap"]=np.log(cs["market_cap_mm"])
cs=cs.sort_values("event_id").reset_index(drop=True)
X=sm.add_constant(cs[["log_units","media_articles","severity_score","log_mcap"]].values)
y=cs["scar_3day"].values
N=len(y); L=int(np.floor(4*(N/100)**(2/9)))
m=sm.OLS(y,X).fit()
hac=m.get_robustcov_results(cov_type='HAC',maxlags=L,use_correction=False)
print("N",N,"L",L)
print("coef ",np.round(m.params,6))
print("HAC se ",np.round(hac.bse,6))
print("R2 ",round(m.rsquared,6))
# WLS
w=1/cs["idio_vol"].values**2
mw=sm.WLS(y,X,weights=w).fit()
print("WLS coef ",np.round(mw.params,6))
print("WLS se ",np.round(mw.bse,6))
print("WLS R2 ",round(mw.rsquared,6))
EOFN 35 L 3 coef [ 6.870415 -0.277206 -0.140035 -0.049356 -0.349665] HAC se [1.989976 0.143768 0.045345 0.079685 0.109201] R2 0.365807 WLS coef [ 6.181414 -0.245301 -0.142987 -0.068882 -0.301826] WLS se [2.373976 0.119944 0.063018 0.074244 0.192047] WLS R2 0.307908
cd /app && python3 << 'EOF'
import numpy as np, pandas as pd
from scipy import stats
r = pd.read_csv("data/stock_returns.csv")
r = r.dropna(subset=["return","market_return"])
r = r[r["return"]>=-10]
r = r.drop_duplicates()
r["date"]=pd.to_datetime(r["date"])
r = r.sort_values(["firm_id","date"]).reset_index(drop=True)
rec = pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id")
rec["date"]=pd.to_datetime(rec["date"])
all_dates = np.sort(r["date"].unique())
d2i = {pd.Timestamp(d):i for i,d in enumerate(all_dates)}
nD=len(all_dates)
ftab={f:g.set_index(g["date"].astype(str)) for f,g in r.groupby("firm_id")}
EST_LEN,EST_GAP=200,30
def fit(fid,eidx):
ee=eidx-EST_GAP-1; es=ee-EST_LEN+1
if es<0: return None
ed=set(pd.Timestamp(all_dates[k]) for k in range(es,ee+1))
g=ftab[fid]; sub=g[g["date"].isin(ed)]
if len(sub)<100: return None
x=sub["market_return"].values; y=sub["return"].values
b1,b0=np.polyfit(x,y,1) # slope,intercept
res=y-(b0+b1*x); L=len(sub); s2=np.sum(res**2)/(L-2)
return dict(a=b0,b=b1,s2=s2,L=L,mrm=x.mean(),ssrm=np.sum((x-x.mean())**2),phat=np.mean(res>0),n=L)
wins={"w3":(-1,1),"w2":(0,1),"w11":(-5,5)}
rows=[]
for _,row in rec.iterrows():
fid=row["firm_id"]; ts=pd.Timestamp(row["date"])
if ts not in d2i: continue
eidx=d2i[ts]; m=fit(fid,eidx)
if m is None: continue
g=ftab[fid]; out={"event_id":int(row["event_id"]),"firm_id":fid}
okw3=True
for wn,(lo,hi) in wins.items():
asum=0;vsum=0;ok=True
for off in range(lo,hi+1):
ti=eidx+off
if ti<0 or ti>=nD: ok=False;break
tgt=str(pd.Timestamp(all_dates[ti]))
if tgt not in g.index: ok=False;break
rm=g.loc[tgt,"market_return"]; rt=g.loc[tgt,"return"]
ar=rt-(m["a"]+m["b"]*rm); asum+=ar
vsum+=m["s2"]*(1+1/m["L"]+(rm-m["mrm"])**2/m["ssrm"])
out["scar_"+wn]= (asum/np.sqrt(vsum)) if ok else np.nan
if wn=="w3" and not ok: okw3=False
if not okw3 or np.isnan(out["scar_w3"]): continue
rows.append(out)
ev=pd.DataFrame(rows)
N=len(ev)
def agg(v):
v=v.dropna().values;n=len(v)
return dict(mean=v.mean(),patell=v.sum()/np.sqrt(n),bmp=v.mean()/(v.std(ddof=1)/np.sqrt(n)),pctneg=np.mean(v<0),n=n)
for wn in ["w3","w2","w11"]:
a=agg(ev["scar_"+wn]); print(wn,{k:round(val,4) if isinstance(val,float) else val for k,val in a.items()})
# Schipper-Thompson
S=ev[["scar_w3","scar_w2","scar_w11"]].dropna().values
mu=S.mean(0); cov=np.cov(S,rowvar=False,ddof=1)
W=S.shape[0]*mu@np.linalg.inv(cov)@mu
print("Schipper Wald",round(W,4),"p",stats.chi2.sf(W,3))
# Calendar-time
efi=[(row["firm_id"],d2i[pd.Timestamp(row["date"])]) for _,row in rec.iterrows() if pd.Timestamp(row["date"]) in d2i]
mbyd={str(pd.Timestamp(d)):None for d in all_dates}
for d,grp in r.groupby(r["date"].astype(str)):
mbyd[d]=grp["market_return"].iloc[0]
prows=[]
for g in range(nD):
dstr=str(pd.Timestamp(all_dates[g]))
fw=set(f for f,ei in efi if abs(ei-g)<=30)
if not fw: continue
rr=[ftab[f].loc[dstr,"return"] for f in fw if dstr in ftab[f].index]
if not rr: continue
prows.append((all_dates[g],np.mean(rr),mbyd[dstr]))
pdf=pd.DataFrame(prows,columns=["date","pr","mr"]).sort_values("date")
import statsmodels.api as sm
mm=sm.OLS(pdf["pr"].values,sm.add_constant(pdf["mr"].values)).fit()
print("Calendar alpha",round(mm.params[0],6),"beta",round(mm.params[1],4),"t_a",round(mm.tvalues[0],4),"p",round(mm.pvalues[0],4),"ndays",len(pdf))
EOF[stderr]
Traceback (most recent call last):
File "<stdin>", line 59, in <module>
File "/usr/local/lib/python3.11/site-packages/pandas/core/frame.py", line 4102, in __getitem__
indexer = self.columns.get_loc(key)
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.11/site-packages/pandas/core/indexes/range.py", line 417, in get_loc
raise KeyError(key)
KeyError: 'scar_w3'cd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd, sys
from scipy import stats
r = pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r = r[r["return"]>=-10].drop_duplicates()
r["date"]=pd.to_datetime(r["date"]); r=r.sort_values(["firm_id","date"]).reset_index(drop=True)
rec = pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id")
rec["date"]=pd.to_datetime(rec["date"])
all_dates=np.sort(r["date"].unique()); d2i={pd.Timestamp(d):i for i,d in enumerate(all_dates)}; nD=len(all_dates)
ftab={f:g.set_index(g["date"].astype(str)) for f,g in r.groupby("firm_id")}
def fit(fid,eidx):
ee=eidx-31; es=ee-199
if es<0: return None
ed=set(pd.Timestamp(all_dates[k]) for k in range(es,ee+1))
g=ftab[fid]; sub=g[g["date"].isin(ed)]
if len(sub)<100: return None
x=sub["market_return"].values; y=sub["return"].values
b1,b0=np.polyfit(x,y,1); res=y-(b0+b1*x); L=len(sub)
return dict(a=b0,b=b1,s2=np.sum(res**2)/(L-2),L=L,mrm=x.mean(),ssrm=np.sum((x-x.mean())**2),phat=np.mean(res>0),n=L)
wins={"scar_3day":(-1,1),"scar_2day":(0,1),"scar_11day":(-5,5)}
rows=[]
for _,row in rec.iterrows():
fid=row["firm_id"]; ts=pd.Timestamp(row["date"])
if ts not in d2i: continue
eidx=d2i[ts]; m=fit(fid,eidx)
if m is None: continue
g=ftab[fid]; out={"event_id":int(row["event_id"]),"firm_id":fid}; okw3=True
for wn,(lo,hi) in wins.items():
asum=0;vsum=0;ok=True
for off in range(lo,hi+1):
ti=eidx+off
if ti<0 or ti>=nD: ok=False;break
tgt=str(pd.Timestamp(all_dates[ti]))
if tgt not in g.index: ok=False;break
rm=g.loc[tgt,"market_return"]; rt=g.loc[tgt,"return"]
asum+=rt-(m["a"]+m["b"]*rm); vsum+=m["s2"]*(1+1/m["L"]+(rm-m["mrm"])**2/m["ssrm"])
out[wn]=(asum/np.sqrt(vsum)) if ok else np.nan
if wn=="scar_3day" and not ok: okw3=False
if not okw3 or np.isnan(out["scar_3day"]): continue
rows.append(out)
ev=pd.DataFrame(rows); print("N events",len(ev),"cols",list(ev.columns))
def agg(v):
v=v.dropna().values;n=len(v)
return dict(mean=round(v.mean(),4),patell=round(v.sum()/np.sqrt(n),4),bmp=round(v.mean()/(v.std(ddof=1)/np.sqrt(n)),4),pctneg=round(np.mean(v<0),4),n=n)
for wn in wins: print(wn,agg(ev[wn]))
S=ev[list(wins)].dropna().values; mu=S.mean(0); cov=np.cov(S,rowvar=False,ddof=1)
W=S.shape[0]*mu@np.linalg.inv(cov)@mu; print("Schipper Wald",round(W,4),"p",stats.chi2.sf(W,3))
EOFN events 0 cols []
Traceback (most recent call last):
File "<stdin>", line 44, in <module>
File "/usr/local/lib/python3.11/site-packages/pandas/core/frame.py", line 4102, in __getitem__
indexer = self.columns.get_loc(key)
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.11/site-packages/pandas/core/indexes/range.py", line 417, in get_loc
raise KeyError(key)
KeyError: 'scar_3day'cd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd
r = pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r = r[r["return"]>=-10].drop_duplicates()
r["date"]=pd.to_datetime(r["date"]); r=r.sort_values(["firm_id","date"]).reset_index(drop=True)
rec = pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id")
rec["date"]=pd.to_datetime(rec["date"])
all_dates=np.sort(r["date"].unique()); d2i={pd.Timestamp(d):i for i,d in enumerate(all_dates)}; nD=len(all_dates)
ftab={f:g.set_index(g["date"].astype(str)) for f,g in r.groupby("firm_id")}
row=rec.iloc[0]; fid=row["firm_id"]; ts=pd.Timestamp(row["date"])
print("fid",fid,"ts",ts,"in d2i",ts in d2i)
eidx=d2i[ts]; print("eidx",eidx)
ee=eidx-31; es=ee-199; print("es",es,"ee",ee)
ed=set(pd.Timestamp(all_dates[k]) for k in range(es,ee+1)); print("len ed",len(ed))
g=ftab[fid]; print("g date dtype",g["date"].dtype,"sample",list(g["date"].head(2)))
sub=g[g["date"].isin(ed)]; print("len sub",len(sub))
# try alternative: compare types
print("sample ed elem type", type(next(iter(ed))))
EOFfid FIRM_007 ts 2007-03-21 00:00:00 in d2i True
eidx 316
es 86 ee 285
len ed 200
g date dtype datetime64[ns] sample [Timestamp('2006-01-03 00:00:00'), Timestamp('2006-01-04 00:00:00')]
len sub 200
sample ed elem type <class 'pandas._libs.tslibs.timestamps.Timestamp'>cd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd
from scipy import stats
import statsmodels.api as sm
r = pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r = r[r["return"]>=-10].drop_duplicates()
r["date"]=pd.to_datetime(r["date"]); r["ds"]=r["date"].dt.strftime("%Y-%m-%d")
r=r.sort_values(["firm_id","date"]).reset_index(drop=True)
rec = pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id"); rec["date"]=pd.to_datetime(rec["date"])
all_dates=np.sort(r["date"].unique()); ds_list=[pd.Timestamp(d).strftime("%Y-%m-%d") for d in all_dates]
d2i={ds_list[i]:i for i in range(len(ds_list))}; nD=len(all_dates)
ftab={f:g.set_index("ds") for f,g in r.groupby("firm_id")}
def fit(fid,eidx):
ee=eidx-31; es=ee-199
if es<0: return None
ed=set(ds_list[k] for k in range(es,ee+1))
g=ftab[fid]; sub=g[g.index.isin(ed)]
if len(sub)<100: return None
x=sub["market_return"].values;y=sub["return"].values
b1,b0=np.polyfit(x,y,1); res=y-(b0+b1*x); L=len(sub)
return dict(a=b0,b=b1,s2=np.sum(res**2)/(L-2),L=L,mrm=x.mean(),ssrm=np.sum((x-x.mean())**2),phat=np.mean(res>0),n=L)
wins={"scar_3day":(-1,1),"scar_2day":(0,1),"scar_11day":(-5,5)}
rows=[]; models={}
for _,row in rec.iterrows():
fid=row["firm_id"]; ts=row["date"].strftime("%Y-%m-%d")
if ts not in d2i: continue
eidx=d2i[ts]; m=fit(fid,eidx)
if m is None: continue
g=ftab[fid]; out={"event_id":int(row["event_id"]),"firm_id":fid}; okw3=True; ar3=[]
for wn,(lo,hi) in wins.items():
asum=0;vsum=0;ok=True
for off in range(lo,hi+1):
ti=eidx+off
if ti<0 or ti>=nD: ok=False;break
tgt=ds_list[ti]
if tgt not in g.index: ok=False;break
rm=g.loc[tgt,"market_return"];rt=g.loc[tgt,"return"]; ar=rt-(m["a"]+m["b"]*rm); asum+=ar
vsum+=m["s2"]*(1+1/m["L"]+(rm-m["mrm"])**2/m["ssrm"])
if wn=="scar_3day": ar3.append(ar)
out[wn]=(asum/np.sqrt(vsum)) if ok else np.nan
if wn=="scar_3day" and not ok: okw3=False
if not okw3 or np.isnan(out["scar_3day"]): continue
out["phat"]=m["phat"]; out["ar3"]=ar3; out["n"]=m["n"]; rows.append(out); models[out["event_id"]]=(m,fid,eidx)
ev=pd.DataFrame(rows); N=len(ev); print("N events",N)
def agg(v):
v=np.array(v);v=v[~np.isnan(v)];n=len(v)
return dict(mean=round(v.mean(),4),patell=round(v.sum()/np.sqrt(n),4),bmp=round(v.mean()/(v.std(ddof=1)/np.sqrt(n)),4),pctneg=round(float(np.mean(v<0)),4),n=n)
for wn in wins: print(wn,agg(ev[wn]))
# rho_bar
resid={}
for f,g in ftab.items():
if len(g)<30: continue
x=g["market_return"].values;y=g["return"].values;b1,b0=np.polyfit(x,y,1)
resid[f]=pd.Series(y-(b0+b1*x),index=g.index)
fs=list(resid); pc=[]
for i in range(len(fs)):
for j in range(i+1,len(fs)):
a=resid[fs[i]];b=resid[fs[j]];com=a.index.intersection(b.index)
if len(com)<20: continue
pc.append(np.corrcoef(a[com],b[com])[0,1])
rho=np.mean(pc); print("rho_bar",round(rho,6),"npairs",len(pc))
for wn in wins:
v=ev[wn].dropna().values;n=len(v);kp=v.mean()/np.sqrt(np.var(v,ddof=1)*(1+(n-1)*rho)/n)
print(wn,"kp_t",round(kp,4))
# Schipper
S=ev[list(wins)].dropna().values;mu=S.mean(0);cov=np.cov(S,rowvar=False,ddof=1)
W=S.shape[0]*mu@np.linalg.inv(cov)@mu;print("Schipper W",round(W,4),"p",stats.chi2.sf(W,3))
# Cowan
pbar=ev["phat"].mean(); I=np.array([1 if any(a>0 for a in row) else 0 for row in ev["ar3"]]); fp=I.mean()
cz=(fp-pbar)/np.sqrt(pbar*(1-pbar)/N); print("Cowan z",round(cz,4),"p",round(2*stats.norm.sf(abs(cz)),4),"pbar",round(pbar,4),"fp",round(fp,4))
# Corrado
K=[]
for _,row in rec.iterrows():
fid=row["firm_id"];ts=row["date"].strftime("%Y-%m-%d")
if ts not in d2i: continue
eidx=d2i[ts];g=ftab[fid];offs=[];rts=[]
for off in range(-30,31):
ti=eidx+off
if ti<0 or ti>=nD: continue
tgt=ds_list[ti]
if tgt not in g.index: continue
offs.append(off);rt
…[truncated]N events 35
scar_3day {'mean': -1.0116, 'patell': -5.9849, 'bmp': -4.9921, 'pctneg': 0.8286, 'n': 35}
scar_2day {'mean': -0.8622, 'patell': -5.1011, 'bmp': -4.51, 'pctneg': 0.8, 'n': 35}
scar_11day {'mean': -0.4507, 'patell': -2.6665, 'bmp': -2.1749, 'pctneg': 0.6286, 'n': 35}
rho_bar 0.005682 npairs 153
scar_3day kp_t -4.5702
scar_2day kp_t -4.1287
scar_11day kp_t -1.991
Schipper W 25.6289 p 1.1405555227583067e-05
Cowan z 0.2535 p 0.7998 pbar 0.4929 fp 0.5143
Corrado z -2.4265 p 0.0152 N 35
Calendar alpha -0.000858 beta 1.0332 t -1.3457 p 0.1797 nd 230cd /app && Rscript -e '
W <- 25.6289
p <- pchisq(W, 3, lower.tail=FALSE)
cat("wald_p =", p, "\n")
cat("round6 =", round(p,6), "\n")
library(jsonlite)
cat("toJSON round6:", toJSON(round(p,6), auto_unbox=TRUE), "\n")
cat("toJSON raw:", toJSON(p, auto_unbox=TRUE), "\n")
cat("toJSON digits:", toJSON(round(p,6), auto_unbox=TRUE, digits=10), "\n")
'wald_p = 1.140538e-05 round6 = 1.1e-05 toJSON round6: 0 toJSON raw: 0 toJSON digits: 1.1e-05
command: str_replace
path: /app/analysis.R
old_str: dir.create("/app/outputs", showWarnings = FALSE, recursive = TRUE)
rnd <- function(x) if (is.numeric(x) && !is.integer(x)) round(x, 6) else x
results <- list(
n_events = as.integer(n_events),
n_firms = as.integer(n_firms),
n_firms_with_recalls = as.integer(n_firms_with_recalls),
n_valid_events = as.integer(n_valid_events),
rho_bar = round(rho_bar, 6),
aggregate_scars = list(
window_3day = lapply(agg_3, rnd),
window_2day = lapply(agg_2, rnd),
window_11day = lapply(agg_11, rnd)
),
cross_section_hac = lapply(cs_hac, rnd),
cross_section_wls = lapply(cs_wls, rnd),
calendar_time = lapply(calendar, rnd),
corrado_rank = lapply(corrado, rnd),
cowan_sign = lapply(cowan, rnd),
schipper_thompson = lapply(schipper, rnd)
)
write(toJSON(results, auto_unbox = TRUE, pretty = TRUE),
"/app/outputs/event_study_results.json")
new_str: dir.create("/app/outputs", showWarnings = FALSE, recursive = TRUE)
results <- list(
n_events = as.integer(n_events),
n_firms = as.integer(n_firms),
n_firms_with_recalls = as.integer(n_firms_with_recalls),
n_valid_events = as.integer(n_valid_events),
rho_bar = rho_bar,
aggregate_scars = list(
window_3day = agg_3,
window_2day = agg_2,
window_11day = agg_11
),
cross_section_hac = cs_hac,
cross_section_wls = cs_wls,
calendar_time = calendar,
corrado_rank = corrado,
cowan_sign = cowan,
schipper_thompson = schipper
)
# digits = NA preserves full double precision (default digits = 4 would
# truncate small values such as p-values and near-zero coefficients to 0).
write(toJSON(results, auto_unbox = TRUE, pretty = TRUE, digits = NA),
"/app/outputs/event_study_results.json")The file /app/analysis.R has been edited successfully.
cd /app && rm -f outputs/* && Rscript analysis.R 2>&1 && echo "=== JSON ===" && cat outputs/event_study_results.json
Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3
=== JSON ===
{
"n_events": 35,
"n_firms": 18,
"n_firms_with_recalls": 14,
"n_valid_events": 35,
"rho_bar": 0.0056821459330568,
"aggregate_scars": {
"window_3day": {
"mean_scar": -1.0116304539979,
"patell_z": -5.98488647686484,
"bmp_t": -4.99213776186089,
"kp_t": -4.57015805533409,
"pct_negative": 0.828571428571429,
"n": 35
},
"window_2day": {
"mean_scar": -0.862242761104183,
"patell_z": -5.10109696709245,
"bmp_t": -4.50996935631404,
"kp_t": -4.12874679471694,
"pct_negative": 0.8,
"n": 35
},
"window_11day": {
"mean_scar": -0.450723180620522,
"patell_z": -2.66651429664343,
"bmp_t": -2.17487494689865,
"kp_t": -1.99103525023878,
"pct_negative": 0.628571428571429,
"n": 35
}
},
"cross_section_hac": {
"intercept": 6.87041540570889,
"log_units": -0.277206000865995,
"media": -0.140035223934644,
"severity": -0.0493558669041088,
"log_mcap": -0.349664887829838,
"se_intercept": 1.989976224197,
"se_log_units": 0.143768062566243,
"se_media": 0.045344608521185,
"se_severity": 0.0796849034968555,
"se_log_mcap": 0.109201320305086,
"r_squared": 0.365806981865842,
"n_obs": 35,
"hac_lag": 3
},
"cross_section_wls": {
"intercept": 6.18141370419048,
"log_units": -0.245301252974036,
"media": -0.142987453304122,
"severity": -0.0688815591352763,
"log_mcap": -0.301826107135229,
"se_intercept": 2.37397582785115,
"se_log_units": 0.119944251522637,
"se_media": 0.0630179237202238,
"se_severity": 0.0742443322545675,
"se_log_mcap": 0.192046874929762,
"r_squared": 0.307908153604217,
"n_obs": 35
},
"calendar_time": {
"alpha": -0.000857770921066497,
"beta": 1.03320264843221,
"t_alpha": -1.34565994573033,
"p_value": 0.179748847934933,
"n_days": 230
},
"corrado_rank": {
"z_stat": -2.42651806653002,
"p_value": 0.015244491285374
},
"cowan_sign": {
"z_stat": 0.253546652207721,
"p_value": 0.799845815310467,
"p_bar": 0.492859296482412,
"frac_positive": 0.514285714285714
},
"schipper_thompson": {
"wald_stat": 25.6288672616647,
"wald_p": 1.14055552275834e-05,
"df": 3
}
}cd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd
r = pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r = r[r["return"]>=-10].drop_duplicates(); r["date"]=pd.to_datetime(r["date"]); r["ds"]=r["date"].dt.strftime("%Y-%m-%d")
rec = pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id"); rec["date"]=pd.to_datetime(rec["date"])
all_dates=np.sort(r["date"].unique()); ds_list=[pd.Timestamp(d).strftime("%Y-%m-%d") for d in all_dates]
d2i={ds_list[i]:i for i in range(len(ds_list))}; nD=len(all_dates)
ftab={f:g.set_index("ds") for f,g in r.groupby("firm_id")}
def fit(fid,eidx):
ee=eidx-31; es=ee-199
ed=set(ds_list[k] for k in range(es,ee+1)); g=ftab[fid]; sub=g[g.index.isin(ed)]
x=sub["market_return"].values;y=sub["return"].values; b1,b0=np.polyfit(x,y,1)
res=y-(b0+b1*x); L=len(sub)
return dict(a=b0,b=b1,s2=np.sum(res**2)/(L-2),L=L,mrm=x.mean(),ssrm=np.sum((x-x.mean())**2))
wins={"3day":(-1,1),"2day":(0,1),"11day":(-5,5)}
resA={w:[] for w in wins}; resB={w:[] for w in wins}; resC={w:[] for w in wins}
for _,row in rec.iterrows():
fid=row["firm_id"]; ts=row["date"].strftime("%Y-%m-%d"); eidx=d2i[ts]; m=fit(fid,eidx); g=ftab[fid]
for wn,(lo,hi) in wins.items():
ars=[];vs=[];devs=[];ok=True
for off in range(lo,hi+1):
tgt=ds_list[eidx+off]
if tgt not in g.index: ok=False;break
rm=g.loc[tgt,"market_return"];rt=g.loc[tgt,"return"]; ars.append(rt-(m["a"]+m["b"]*rm))
vs.append(m["s2"]*(1+1/m["L"]+(rm-m["mrm"])**2/m["ssrm"])); devs.append(rm-m["mrm"])
if not ok: continue
ars=np.array(ars);vs=np.array(vs);T=len(ars)
resA[wn].append(np.sum(ars/np.sqrt(vs))/np.sqrt(T))
resB[wn].append(np.sum(ars)/np.sqrt(np.sum(vs)))
varC=m["s2"]*(T + T**2/m["L"] + (np.sum(devs))**2/m["ssrm"])
resC[wn].append(np.sum(ars)/np.sqrt(varC))
for wn in wins:
for name,res in [("A sum(SAR)/sqrtT",resA),("B sumAR/sqrt(sumV)",resB),("C exact-PE",resC)]:
v=np.array(res);print(f"{wn:6s} {name:20s} mean={v.mean():.4f} patell={v.sum()/np.sqrt(len(v)):.4f}")
print()
EOFTraceback (most recent call last):
File "<stdin>", line 33, in <module>
File "/usr/local/lib/python3.11/site-packages/numpy/core/_methods.py", line 131, in _mean
ret = ret / rcount
~~~~^~~~~~~~
TypeError: unsupported operand type(s) for /: 'dict' and 'int'cd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd
r = pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r = r[r["return"]>=-10].drop_duplicates(); r["date"]=pd.to_datetime(r["date"]); r["ds"]=r["date"].dt.strftime("%Y-%m-%d")
rec = pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id"); rec["date"]=pd.to_datetime(rec["date"])
all_dates=np.sort(r["date"].unique()); ds_list=[pd.Timestamp(d).strftime("%Y-%m-%d") for d in all_dates]
d2i={ds_list[i]:i for i in range(len(ds_list))}; ftab={f:g.set_index("ds") for f,g in r.groupby("firm_id")}
def fit(fid,eidx):
ee=eidx-31; es=ee-199; ed=set(ds_list[k] for k in range(es,ee+1)); g=ftab[fid]; sub=g[g.index.isin(ed)]
x=sub["market_return"].values;y=sub["return"].values; b1,b0=np.polyfit(x,y,1); res=y-(b0+b1*x); L=len(sub)
return dict(a=b0,b=b1,s2=np.sum(res**2)/(L-2),L=L,mrm=x.mean(),ssrm=np.sum((x-x.mean())**2))
wins={"3day":(-1,1),"2day":(0,1),"11day":(-5,5)}
out={"A":{w:[] for w in wins},"B":{w:[] for w in wins},"C":{w:[] for w in wins}}
for _,row in rec.iterrows():
fid=row["firm_id"]; eidx=d2i[row["date"].strftime("%Y-%m-%d")]; m=fit(fid,eidx); g=ftab[fid]
for wn,(lo,hi) in wins.items():
ars=[];vs=[];devs=[];ok=True
for off in range(lo,hi+1):
tgt=ds_list[eidx+off]
if tgt not in g.index: ok=False;break
rm=g.loc[tgt,"market_return"];rt=g.loc[tgt,"return"]; ars.append(rt-(m["a"]+m["b"]*rm))
vs.append(m["s2"]*(1+1/m["L"]+(rm-m["mrm"])**2/m["ssrm"])); devs.append(rm-m["mrm"])
if not ok: continue
ars=np.array(ars);vs=np.array(vs);T=len(ars)
out["A"][wn].append(np.sum(ars/np.sqrt(vs))/np.sqrt(T))
out["B"][wn].append(np.sum(ars)/np.sqrt(np.sum(vs)))
out["C"][wn].append(np.sum(ars)/np.sqrt(m["s2"]*(T+T**2/m["L"]+(np.sum(devs))**2/m["ssrm"])))
for wn in wins:
for name in ["A","B","C"]:
v=np.array(out[name][wn]); lab={"A":"sum(SAR)/sqrtT","B":"sumAR/sqrt(sumV)","C":"exact-PE"}[name]
print(f"{wn:6s} {lab:18s} mean={v.mean():+.4f} patell={v.sum()/np.sqrt(len(v)):+.4f} bmp={v.mean()/(v.std(ddof=1)/np.sqrt(len(v))):+.4f}")
print()
EOF3day sum(SAR)/sqrtT mean=-1.0118 patell=-5.9862 bmp=-4.9935 3day sumAR/sqrt(sumV) mean=-1.0116 patell=-5.9849 bmp=-4.9921 3day exact-PE mean=-1.0075 patell=-5.9606 bmp=-4.9996 2day sum(SAR)/sqrtT mean=-0.8627 patell=-5.1038 bmp=-4.5132 2day sumAR/sqrt(sumV) mean=-0.8622 patell=-5.1011 bmp=-4.5100 2day exact-PE mean=-0.8609 patell=-5.0933 bmp=-4.5148 11day sum(SAR)/sqrtT mean=-0.4521 patell=-2.6747 bmp=-2.1807 11day sumAR/sqrt(sumV) mean=-0.4507 patell=-2.6665 bmp=-2.1749 11day exact-PE mean=-0.4396 patell=-2.6007 bmp=-2.1751
cd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd
r = pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r = r[r["return"]>=-10].drop_duplicates(); r["date"]=pd.to_datetime(r["date"]); r["ds"]=r["date"].dt.strftime("%Y-%m-%d")
rec = pd.read_csv("data/recalls.csv").drop_duplicates();
ftab={f:g.set_index("ds") for f,g in r.groupby("firm_id")}
event_firms=set(rec["firm_id"])
def rho(firm_subset):
resid={}
for f,g in ftab.items():
if firm_subset is not None and f not in firm_subset: continue
if len(g)<30: continue
x=g["market_return"].values;y=g["return"].values;b1,b0=np.polyfit(x,y,1); resid[f]=pd.Series(y-(b0+b1*x),index=g.index)
fs=list(resid); pc=[]
for i in range(len(fs)):
for j in range(i+1,len(fs)):
a=resid[fs[i]];b=resid[fs[j]];com=a.index.intersection(b.index)
if len(com)<20: continue
pc.append(np.corrcoef(a[com],b[com])[0,1])
return np.mean(pc),len(pc)
print("all 18 firms: rho_bar=%.6f npairs=%d"%rho(None))
print("14 event firms: rho_bar=%.6f npairs=%d"%rho(event_firms))
# impact on kp for 3day (bmp=-4.9921, N=35)
for rb in [rho(None)[0], rho(event_firms)[0]]:
N=35; bmp=-4.9921; kp=bmp/np.sqrt(1+(N-1)*rb); print("rho=%.6f -> kp_3day=%.4f"%(rb,kp))
EOFall 18 firms: rho_bar=0.005682 npairs=153 14 event firms: rho_bar=0.006605 npairs=91 rho=0.005682 -> kp_3day=-4.5701 rho=0.006605 -> kp_3day=-4.5112
command: view path: /app/outputs/car_timeline.png
<output_image> …[truncated]
cd /app && Rscript -e 'cat("ggplot2 version:", as.character(packageVersion("ggplot2")), "\n")'ggplot2 version: 3.4.1
command: str_replace
path: /app/analysis.R
old_str: x = "Trading days relative to recall announcement (day 0)",
y = "Cumulative abnormal return") +
theme_minimal()
ggsave("/app/outputs/car_timeline.png", p1, width = 8, height = 5, dpi = 120)
new_str: x = "Trading days relative to recall announcement (day 0)",
y = "Cumulative abnormal return") +
theme_minimal(base_size = 13) +
theme(plot.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA))
ggsave("/app/outputs/car_timeline.png", p1, width = 8, height = 5, dpi = 120, bg = "white")The file /app/analysis.R has been edited successfully.
command: str_replace
path: /app/analysis.R
old_str: x = "Date", y = "Cumulative simple return", color = "Series") +
theme_minimal() +
theme(legend.position = "bottom")
ggsave("/app/outputs/portfolio_comparison.png", p2, width = 10, height = 5, dpi = 120)
new_str: x = "Date", y = "Cumulative simple return", color = "Series") +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom",
plot.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA))
ggsave("/app/outputs/portfolio_comparison.png", p2, width = 10, height = 5, dpi = 120, bg = "white")The file /app/analysis.R has been edited successfully.
cd /app && rm -f outputs/*.png && Rscript analysis.R 2>&1 | tail -1
Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3
command: view path: /app/outputs/car_timeline.png
<output_image> …[truncated]
command: view path: /app/outputs/portfolio_comparison.png
<output_image> …[truncated]
command: str_replace
path: /app/analysis.R
old_str: labs(title = "Average Cumulative Abnormal Return around Recall Announcements",
subtitle = "Mean CAR across events with 95% confidence bands",
new_str: labs(title = "Average Cumulative Abnormal Return Around Recalls",
subtitle = "Mean CAR across events with 95% confidence bands",The file /app/analysis.R has been edited successfully.
cd /app && rm -f outputs/car_timeline.png && Rscript analysis.R 2>&1 | tail -1
Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3
command: view path: /app/outputs/car_timeline.png
<output_image> …[truncated]
cd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd, json
r = pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r = r[r["return"]>=-10].drop_duplicates(); r["date"]=pd.to_datetime(r["date"]); r["ds"]=r["date"].dt.strftime("%Y-%m-%d")
rec = pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id"); rec["date"]=pd.to_datetime(rec["date"])
all_dates=np.sort(r["date"].unique()); ds_list=[pd.Timestamp(d).strftime("%Y-%m-%d") for d in all_dates]
d2i={ds_list[i]:i for i in range(len(ds_list))}; ftab={f:g.set_index("ds") for f,g in r.groupby("firm_id")}
def fit(fid,eidx):
ee=eidx-31; es=ee-199; ed=set(ds_list[k] for k in range(es,ee+1)); g=ftab[fid]; sub=g[g.index.isin(ed)]
x=sub["market_return"].values;y=sub["return"].values; b1,b0=np.polyfit(x,y,1); return b0,b1
wins={"3day":(-1,1),"2day":(0,1),"11day":(-5,5)}
raw={w:[] for w in wins}
for _,row in rec.iterrows():
fid=row["firm_id"]; eidx=d2i[row["date"].strftime("%Y-%m-%d")]; a,b=fit(fid,eidx); g=ftab[fid]
for wn,(lo,hi) in wins.items():
s=0
for off in range(lo,hi+1):
tgt=ds_list[eidx+off]; rm=g.loc[tgt,"market_return"];rt=g.loc[tgt,"return"]; s+=rt-(a+b*rm)
raw[wn].append(s)
for wn in wins:
v=np.array(raw[wn]); print(f"{wn}: mean raw CAR = {v.mean()*100:+.2f}% (median {np.median(v)*100:+.2f}%)")
# cross-section t-stats (HAC)
j=json.load(open("outputs/event_study_results.json"))
h=j["cross_section_hac"]
for nm in ["log_units","media","severity","log_mcap","intercept"]:
print(f"HAC {nm}: coef={h[nm]:+.4f} se={h['se_'+nm]:.4f} t={h[nm]/h['se_'+nm]:+.2f}")
EOF3day: mean raw CAR = -2.92% (median -2.77%) 2day: mean raw CAR = -2.05% (median -2.22%) 11day: mean raw CAR = -2.50% (median -2.13%) HAC log_units: coef=-0.2772 se=0.1438 t=-1.93 HAC media: coef=-0.1400 se=0.0453 t=-3.09 HAC severity: coef=-0.0494 se=0.0797 t=-0.62 HAC log_mcap: coef=-0.3497 se=0.1092 t=-3.20 HAC intercept: coef=+6.8704 se=1.9900 t=+3.45
command: create path: /app/outputs/research_memo.md file_text: # The Stock-Market Impact of Product-Recall Announcements in the Toy Industry **Prepared for:** Non-technical readers (regulators, investors, corporate boards) **Subject:** How share prices respond when toy manufacturers announce product recalls **Sample:** 35 recall events involving 14 of 18 listed toy manufacturers during 2007 --- ## 1. Executive summary When a toy manufacturer announces a product recall, its share price falls sharply and almost immediately. Across the 35 recalls we studied, the average company lost roughly **2.9% of its market value in the three trading days centered on the announcement** (the day before, the day of, and the day after). This drop is large, consistent across events, and highly unlikely to be a coincidence. The market, in other words, treats a recall as genuinely bad news about a company's future profitability , not merely a public-relations nuisance. ## 2. Magnitude and statistical significance of the abnormal returns We measure the "abnormal" return , the portion of the stock move that cannot be explained by the overall market that day , using a market model estimated over the 200 trading days ending one month before each recall. We then standardize each company's cumulative abnormal return (CAR) by its own prediction error, so that large, volatile firms and small, stable firms are placed on a comparable footing. The results are strong and stable across event windows: | Window | Avg. raw CAR | Mean standardized CAR | Patell Z | BMP t | Kolari-Pynnonen t | % negative | |--------|-------------:|----------------------:|---------:|------:|------------------:|-----------:| | [-1,+1] (3-day) | -2.9% | -1.01 | -5.98 | -4.99 | -4.57 | 83% | | [0,+1] (2-day) | -2.0% | -0.86 | -5.10 | -4.51 | -4.13 | 80% | | [-5,+5] (11-day)| -2.5% | -0.45 | -2.67 | -2.17 | -1.99 | 63% | All three test statistics are well beyond the conventional significance threshold (about +/-2) for the 2- and 3-day windows. The Patell Z assumes each standardized return is an independent standard-normal draw; the Boehmer-Musumeci-Poulsen (BMP) t-statistic relaxes that by using the actual cross-sectional spread of the returns (guarding against an "event-induced" jump in volatility); and the Kolari-Pynnonen (KP) t-statistic further corrects for the possibility that firms' returns move together. Because the average cross-firm residual correlation is tiny (rho-bar = 0.006), the KP correction barely changes the conclusion. A joint Schipper-Thompson Wald test that all three windows are simultaneously zero is decisively rejected (chi-square = 25.6, p < 0.001). The effect is therefore not an artifact of one particular window or one particular test. The tightest, most informative window is the 3-day window. The 11-day window is also significant but weaker per day, which tells us the reaction is concentrated right at the announcement rather than spread out over two weeks. ## 3. What drives the differences across recalls? Not all recalls are equal. We regressed each event's 3-day standardized return on four characteristics, using heteroskedasticity- and autocorrelation-consistent (Newey-West) standard errors, and separately re-estimated the model weighting each firm by the inverse of its idiosyncratic variance (WLS). The two approaches agree. The model explains a meaningful share of the variation (R-squared ~ 0.37 for OLS, 0.31 for WLS). - **Media coverage is the most robust driver** (coefficient -0.14, t ~ -3.1). Each additional news article is associated with a larger price decline. Press attention both reflects and amplifies how widely investors, retailers, and consumers learn of the defect, raising the expected hit to sales and reputation. - **Firm size matters** (log market-cap coefficient -0.35, t ~ -3.2). Larger toy makers show *more* negative standardized reactions. Part of this is mechanical , bigger firms have lower idiosyncratic volatility, so a …[truncated]
New file created successfully at /app/outputs/research_memo.md
cd /app && echo "Word count: $(wc -w < outputs/research_memo.md)" && echo "---checklist of required topics---" && for kw in "significance" "severity" "media" "firm size" "calendar-time" "Corrado" "Cowan" "distribution" "regulator" "investor" "predictab"; do printf "%-14s: %s\n" "$kw" "$(grep -ic "$kw" outputs/research_memo.md)"; done
Word count: 1265 ---checklist of required topics--- significance : 2 severity : 4 media : 6 firm size : 1 calendar-time : 5 Corrado : 1 Cowan : 1 distribution : 1 regulator : 3 investor : 4 predictab : 1
cd /app && rm -rf outputs && Rscript analysis.R && echo "=== exit code: $? ===" && echo "--- files ---" && ls -la outputs/
Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3 === exit code: 0 === --- files --- total 154 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:07 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:07 .. -rw-r--r-- 1 h2tagent h2tagent 48651 Jun 20 18:07 car_timeline.png -rw-r--r-- 1 h2tagent h2tagent 3618 Jun 20 18:07 event_cars.csv -rw-r--r-- 1 h2tagent h2tagent 2273 Jun 20 18:07 event_study_results.json -rw-r--r-- 1 h2tagent h2tagent 92721 Jun 20 18:07 portfolio_comparison.png
command: str_replace
path: /app/analysis.R
old_str: scars <- list()
ar_day0 <- NA_real_
ar_w3 <- rep(NA_real_, 3) # AR on offsets -1, 0, +1 (for Cowan indicator)
valid_w3 <- TRUE
for (wname in names(windows)) {
w <- windows[[wname]]
ar_sum <- 0; var_sum <- 0; ok <- TRUE
new_str: scars <- list()
raw_car <- list()
ar_day0 <- NA_real_
ar_w3 <- rep(NA_real_, 3) # AR on offsets -1, 0, +1 (for Cowan indicator)
valid_w3 <- TRUE
for (wname in names(windows)) {
w <- windows[[wname]]
ar_sum <- 0; var_sum <- 0; ok <- TRUEThe file /app/analysis.R has been edited successfully.
command: str_replace
path: /app/analysis.R
old_str: if (!ok) {
if (wname == "w3") valid_w3 <- FALSE
scars[[wname]] <- NA_real_
} else {
scars[[wname]] <- ar_sum / sqrt(var_sum) # standardized CAR
}
}
new_str: if (!ok) {
if (wname == "w3") valid_w3 <- FALSE
scars[[wname]] <- NA_real_
raw_car[[wname]] <- NA_real_
} else {
scars[[wname]] <- ar_sum / sqrt(var_sum) # standardized CAR
raw_car[[wname]] <- ar_sum # raw (unstandardized) CAR
}
}The file /app/analysis.R has been edited successfully.
command: str_replace
path: /app/analysis.R
old_str: scar_11day = scars$w11,
ar_day0 = ar_day0,
ar_m1 = ar_w3[1], ar_0 = ar_w3[2], ar_p1 = ar_w3[3],
new_str: scar_11day = scars$w11,
car_3day = raw_car$w3,
car_2day = raw_car$w2,
car_11day = raw_car$w11,
ar_day0 = ar_day0,
ar_m1 = ar_w3[1], ar_0 = ar_w3[2], ar_p1 = ar_w3[3],The file /app/analysis.R has been edited successfully.
command: str_replace
path: /app/analysis.R
old_str: cat("Done. n_valid_events =", n_valid_events, "| rho_bar =", round(rho_bar, 4),
"| hac_lag =", hac_lag, "\n")
new_str: # =====================================================================
# Research memo (interpolated with the computed results)
# =====================================================================
mcar3 <- mean(event_df$car_3day, na.rm = TRUE) * 100
mcar2 <- mean(event_df$car_2day, na.rm = TRUE) * 100
mcar11 <- mean(event_df$car_11day, na.rm = TRUE) * 100
t_units <- cs_hac$log_units / cs_hac$se_log_units
t_media <- cs_hac$media / cs_hac$se_media
t_sev <- cs_hac$severity / cs_hac$se_severity
t_mcap <- cs_hac$log_mcap / cs_hac$se_log_mcap
f <- function(x, d = 2) formatC(x, format = "f", digits = d)
memo <- paste0(
"# The Stock-Market Impact of Product-Recall Announcements in the Toy Industry
**Prepared for:** Non-technical readers (regulators, investors, corporate boards)
**Subject:** How share prices respond when toy manufacturers announce product recalls
**Sample:** ", n_valid_events, " recall events involving ", n_firms_with_recalls,
" of ", n_firms, " listed toy manufacturers during 2007
---
## 1. Executive summary
When a toy manufacturer announces a product recall, its share price falls sharply and
almost immediately. Across the ", n_valid_events, " recalls we studied, the average
company lost roughly **", f(abs(mcar3)), "% of its market value in the three trading days
centered on the announcement** (the day before, the day of, and the day after). This
drop is large, consistent across events, and highly unlikely to be a coincidence. The
market treats a recall as genuinely bad news about future profitability, not merely a
public-relations nuisance.
## 2. Magnitude and statistical significance of the abnormal returns
We measure the *abnormal* return -- the part of a stock's move not explained by the
market that day -- using a market model estimated over the 200 trading days ending one
month before each recall. We then standardize each firm's cumulative abnormal return
(CAR) by its own prediction error, placing large and small firms on a comparable footing.
| Window | Avg. raw CAR | Mean std. CAR | Patell Z | BMP t | KP t | % negative |
|--------|-------------:|--------------:|---------:|------:|-----:|-----------:|
| [-1,+1] (3-day) | ", f(mcar3), "% | ", f(agg_3$mean_scar), " | ", f(agg_3$patell_z), " | ", f(agg_3$bmp_t), " | ", f(agg_3$kp_t), " | ", f(agg_3$pct_negative*100,0), "% |
| [0,+1] (2-day) | ", f(mcar2), "% | ", f(agg_2$mean_scar), " | ", f(agg_2$patell_z), " | ", f(agg_2$bmp_t), " | ", f(agg_2$kp_t), " | ", f(agg_2$pct_negative*100,0), "% |
| [-5,+5] (11-day) | ", f(mcar11),"% | ", f(agg_11$mean_scar), " | ", f(agg_11$patell_z), " | ", f(agg_11$bmp_t), " | ", f(agg_11$kp_t), " | ", f(agg_11$pct_negative*100,0), "% |
All three statistics sit far beyond the usual significance threshold (about +/-2) for the
2- and 3-day windows. The Patell Z treats each standardized return as an independent
standard-normal draw; the Boehmer-Musumeci-Poulsen (BMP) t uses the actual cross-sectional
spread to absorb any event-induced jump in volatility; and the Kolari-Pynnonen (KP) t
further corrects for firms moving together. Because the average cross-firm residual
correlation is tiny (rho-bar = ", f(rho_bar, 4), "), the KP correction barely changes the
verdict. A joint Schipper-Thompson Wald test that all three windows are simultaneously
zero is decisively rejected (chi-square = ", f(schipper$wald_stat), ", p < 0.001). The
effect is not an artifact of any single window or test. The reaction is concentrated right
at the announcement: the 11-day window is significant but weaker per day than the 3-day
window, implying the news is absorbed quickly rather than leaking out over two weeks.
## 3. What drives the differences across recalls?
Not all recalls are equal. We regressed each event's 3-day standardized return on four
char
…[truncated]The file /app/analysis.R has been edited successfully.
cd /app && rm -rf outputs && Rscript analysis.R && echo "=== exit: $? ===" && ls -la outputs/ && echo "memo words: $(wc -w < outputs/research_memo.md)"
Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3 === exit: 0 === total 161 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:08 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:08 .. -rw-r--r-- 1 h2tagent h2tagent 48651 Jun 20 18:08 car_timeline.png -rw-r--r-- 1 h2tagent h2tagent 3618 Jun 20 18:08 event_cars.csv -rw-r--r-- 1 h2tagent h2tagent 2273 Jun 20 18:08 event_study_results.json -rw-r--r-- 1 h2tagent h2tagent 92721 Jun 20 18:08 portfolio_comparison.png -rw-r--r-- 1 h2tagent h2tagent 7418 Jun 20 18:08 research_memo.md memo words: 1149
cd /app && sed -n '1,45p' outputs/research_memo.md
# The Stock-Market Impact of Product-Recall Announcements in the Toy Industry **Prepared for:** Non-technical readers (regulators, investors, corporate boards) **Subject:** How share prices respond when toy manufacturers announce product recalls **Sample:** 35 recall events involving 14 of 18 listed toy manufacturers during 2007 --- ## 1. Executive summary When a toy manufacturer announces a product recall, its share price falls sharply and almost immediately. Across the 35 recalls we studied, the average company lost roughly **2.92% of its market value in the three trading days centered on the announcement** (the day before, the day of, and the day after). This drop is large, consistent across events, and highly unlikely to be a coincidence. The market treats a recall as genuinely bad news about future profitability, not merely a public-relations nuisance. ## 2. Magnitude and statistical significance of the abnormal returns We measure the *abnormal* return -- the part of a stock's move not explained by the market that day -- using a market model estimated over the 200 trading days ending one month before each recall. We then standardize each firm's cumulative abnormal return (CAR) by its own prediction error, placing large and small firms on a comparable footing. | Window | Avg. raw CAR | Mean std. CAR | Patell Z | BMP t | KP t | % negative | |--------|-------------:|--------------:|---------:|------:|-----:|-----------:| | [-1,+1] (3-day) | -2.92% | -1.01 | -5.98 | -4.99 | -4.57 | 83% | | [0,+1] (2-day) | -2.05% | -0.86 | -5.10 | -4.51 | -4.13 | 80% | | [-5,+5] (11-day) | -2.50% | -0.45 | -2.67 | -2.17 | -1.99 | 63% | All three statistics sit far beyond the usual significance threshold (about +/-2) for the 2- and 3-day windows. The Patell Z treats each standardized return as an independent standard-normal draw; the Boehmer-Musumeci-Poulsen (BMP) t uses the actual cross-sectional spread to absorb any event-induced jump in volatility; and the Kolari-Pynnonen (KP) t further corrects for firms moving together. Because the average cross-firm residual correlation is tiny (rho-bar = 0.0057), the KP correction barely changes the verdict. A joint Schipper-Thompson Wald test that all three windows are simultaneously zero is decisively rejected (chi-square = 25.63, p < 0.001). The effect is not an artifact of any single window or test. The reaction is concentrated right at the announcement: the 11-day window is significant but weaker per day than the 3-day window, implying the news is absorbed quickly rather than leaking out over two weeks. ## 3. What drives the differences across recalls?
cd /app && sed -n '46,95p' outputs/research_memo.md && echo "=== JSON validity & key values ===" && python3 -c "
import json
d=json.load(open('outputs/event_study_results.json'))
print('JSON valid. top keys:',list(d.keys()))
assert d['n_events']==35 and d['n_valid_events']==35 and d['n_firms']==18 and d['n_firms_with_recalls']==14
assert d['schipper_thompson']['df']==3
assert d['cross_section_hac']['hac_lag']==3 and d['cross_section_hac']['n_obs']==35
print('calendar p_value:',d['calendar_time']['p_value'],'| wald_p:',d['schipper_thompson']['wald_p'])
print('all structural assertions passed')
"Not all recalls are equal. We regressed each event's 3-day standardized return on four characteristics using Newey-West (HAC) standard errors, and re-estimated the model weighting firms by the inverse of their idiosyncratic variance (WLS). The two agree, and the model explains a meaningful share of the variation (R-squared ~ 0.37 for OLS, 0.31 for WLS). - **Media coverage is the most robust driver** (coefficient -0.14, t ~ -3.09). Each extra news article is linked to a larger decline. Press attention both reflects and amplifies how widely investors, retailers, and consumers learn of the defect, raising the expected hit to sales and reputation. - **Firm size matters** (log market-cap coefficient -0.35, t ~ -3.20). Larger makers show *more* negative standardized reactions. Part is mechanical -- bigger firms have lower idiosyncratic volatility, so a given percentage drop is more standard deviations -- but it also fits the idea that bigger, more visible brands have more reputational capital at stake. - **The scale of the recall** (log units recalled, coefficient -0.28, t ~ -1.93) has the expected negative sign and is marginally significant: bigger recalls hurt more, but noisily. - **The engineered severity score is *not* statistically significant** (t ~ -0.62). Once media coverage and recall scale are accounted for, the hazard-severity label adds little -- the market prices *observable attention and breadth* more than a catalogued rating, plausibly because media volume already impounds public perceptions of danger. ## 4. Calendar-time portfolio versus the short-window event study To test whether the damage persists or is priced immediately, we built a calendar-time portfolio (Jaffe-Mandelker / Fama 1998): each day we held an equally-weighted basket of all firms within 30 trading days of a recall and regressed its return on the market. The portfolio alpha is **-0.086% per day with a t-statistic of just -1.35 (p = 0.18)** over 230 days, with a market beta near 1.03 -- economically negative but *not* statistically significant. This contrast is the key methodological lesson. The event study finds a large, significant effect because it looks exactly where the news lands. The calendar-time portfolio dilutes that one-time drop across a 61-day holding window full of ordinary noise, so the daily signal is small relative to volatility and the test loses power. The gap does **not** mean the event-study result is spurious; it confirms the response is fast and concentrated -- just what we expect in a reasonably efficient market, leaving little predictable post-event drift to exploit. ## 5. Do the non-parametric tests agree? Two distribution-free checks guard against outliers and non-normal returns: - **The Corrado (1989) rank test corroborates the parametric findings**: announcement-window ranks of raw daily returns are abnormally low (z = -2.43, p = 0.015). Because ranks ignore a few extreme values, this agreement shows the negative reaction is a feature of the whole sample, not a handful of crashes. - **The Cowan (1992) generalized sign test does not reject** (z = 0.25, p = 0.80): about 51% of events had at least one positive abnormal-return day in the 3-day window, essentially the 49% baseline. This is not a contradiction. The sign test asks whether positive *days* become rarer, whereas the magnitude tests ask whether *losses* are larger. Their divergence implies the recall effect works through the **size of the down moves (a fat, negative left tail), not a wholesale change in the direction of daily returns** -- investors should expect occasional large drops rather than a steady grind of small losses. ## 6. Practical conclusions === JSON validity & key values === JSON valid. top keys: ['n_events', 'n_firms', 'n_firms_with_recalls', 'n_valid_events', 'rho_bar', 'aggregate_scars', 'cross_section_hac', 'cross_section_wls', 'calendar_time', 'corrado_rank', 'cowan_sign', 'schipper_thompson'] calendar p_value: 0.179748847934933 | wald_p: …[truncated]
cd /app && cp outputs/event_study_results.json /tmp/run1.json && cp outputs/event_cars.csv /tmp/run1.csv && cp outputs/research_memo.md /tmp/run1.md && rm -rf outputs && Rscript analysis.R >/dev/null 2>&1 && echo "JSON diff:" && diff <(cat /tmp/run1.json) outputs/event_study_results.json && echo " (identical)" && echo "CSV diff:" && diff /tmp/run1.csv outputs/event_cars.csv && echo " (identical)" && echo "MEMO diff:" && diff /tmp/run1.md outputs/research_memo.md && echo " (identical)" && echo && echo "CSV header:" && head -1 outputs/event_cars.csv && echo "CSV rows (excl header): $(($(wc -l < outputs/event_cars.csv)-1))"
JSON diff: (identical) CSV diff: (identical) MEMO diff: (identical) CSV header: "event_id","firm_id","date","scar_3day","scar_2day","scar_11day","ar_day0" CSV rows (excl header): 35
cd /app && sed -n '1,130p' analysis.R
suppressPackageStartupMessages({
library(jsonlite)
library(ggplot2)
})
EST_LEN <- 200 # estimation window length (trading days)
EST_GAP <- 30 # gap between estimation window end and event date
CAR_MIN_OBS <- 100 # minimum valid estimation observations for an event
# =====================================================================
# 1. Load and clean
# - drop NA return / market_return
# - drop sentinel returns (return < -10)
# - drop duplicate rows
# - sort by (firm_id, date)
# =====================================================================
returns <- read.csv("/app/data/stock_returns.csv", stringsAsFactors = FALSE)
recalls <- read.csv("/app/data/recalls.csv", stringsAsFactors = FALSE)
firms <- read.csv("/app/data/firms.csv", stringsAsFactors = FALSE)
returns <- returns[!is.na(returns$return) & !is.na(returns$market_return), ]
returns <- returns[returns$return >= -10, ] # remove sentinels
returns <- returns[!duplicated(returns), ] # remove duplicate rows
returns$date <- as.Date(returns$date)
returns <- returns[order(returns$firm_id, returns$date), ]
rownames(returns) <- NULL
recalls <- recalls[!duplicated(recalls), ]
recalls$date <- as.Date(recalls$date)
recalls <- recalls[order(recalls$event_id), ]
rownames(recalls) <- NULL
firms <- firms[!duplicated(firms), ]
n_events <- nrow(recalls)
n_firms <- nrow(firms)
n_firms_with_recalls <- length(unique(recalls$firm_id))
# Global trading-day calendar (shared index across firms)
all_dates <- sort(unique(returns$date))
date_to_idx <- setNames(seq_along(all_dates) - 1L, as.character(all_dates))
n_dates <- length(all_dates)
# Fast per-firm lookup: firm -> (date string -> row)
firm_tables <- split(returns, returns$firm_id)
for (f in names(firm_tables)) rownames(firm_tables[[f]]) <- as.character(firm_tables[[f]]$date)
# =====================================================================
# 2. Market model + prediction-error-corrected SAR / SCAR
#
# Estimation window: 200 trading days ending 30 days before event,
# require >= 100 valid obs. Fit return ~ market_return by OLS.
# s2 = SSE/(L-2). Per-day prediction-error variance:
# S_it^2 = s2 * (1 + 1/L + (Rm_t - mean_Rm)^2 / SS_m)
# SAR_it = AR_it / S_it
# SCAR_i = sum_t AR_it / sqrt( sum_t S_it^2 ) (window cumulation)
# =====================================================================
event_market_model <- function(fid, eidx) {
est_end <- eidx - EST_GAP - 1
est_start <- est_end - EST_LEN + 1
if (est_start < 0) return(NULL)
est_dates <- all_dates[(est_start + 1):(est_end + 1)]
ft <- firm_tables[[fid]]
if (is.null(ft)) return(NULL)
sub <- ft[ft$date %in% est_dates, ]
if (nrow(sub) < CAR_MIN_OBS) return(NULL)
m <- lm(return ~ market_return, data = sub)
res <- unname(resid(m))
L <- nrow(sub)
s2 <- sum(res^2) / (L - 2) # OLS residual variance
mean_rm <- mean(sub$market_return)
ss_rm <- sum((sub$market_return - mean_rm)^2)
list(alpha = unname(coef(m)[1]), beta = unname(coef(m)[2]),
s2 = s2, L = L, mean_rm = mean_rm, ss_rm = ss_rm,
p_hat_pos = mean(res > 0), # Cowan: P(AR>0) over estimation window
n_est = L)
}
# per-day prediction-error variance under the event's fitted model
pe_var <- function(m, rm_t) m$s2 * (1 + 1 / m$L + (rm_t - m$mean_rm)^2 / m$ss_rm)
windows <- list(w3 = c(-1, 1), w2 = c(0, 1), w11 = c(-5, 5))
event_rows <- list()
model_store <- list() # keep fitted models keyed by event_id
for (i in seq_len(n_events)) {
fid <- recalls$firm_id[i]
edate <- recalls$date[i]
estr <- as.character(edate)
if (!(estr %in% names(date_to_idx))) next
eidx <- as.integer(date_to_idx[estr])
m <- event_market_model(fid, eidx)
if (is.null(m)) next
ft <- firm_tables[[fid]]
scars <- list()
raw_car <- list()
ar_day0 <- NA_real_
ar_w3 <- rep(NA_real_, 3) # AR on offsets -1, 0, +1 (for Cowan indicator)
valid_w3 <- TRUE
for (
…[truncated]cd /app && sed -n '285,360p' analysis.R
data = cs, weights = w)
wls_sum <- summary(wls_fit)
wc <- wls_sum$coefficients
cs_wls <- list(
intercept = unname(wc["(Intercept)", "Estimate"]),
log_units = unname(wc["log_units", "Estimate"]),
media = unname(wc["media_articles", "Estimate"]),
severity = unname(wc["severity_score", "Estimate"]),
log_mcap = unname(wc["log_mcap", "Estimate"]),
se_intercept = unname(wc["(Intercept)", "Std. Error"]),
se_log_units = unname(wc["log_units", "Std. Error"]),
se_media = unname(wc["media_articles", "Std. Error"]),
se_severity = unname(wc["severity_score", "Std. Error"]),
se_log_mcap = unname(wc["log_mcap", "Std. Error"]),
r_squared = wls_sum$r.squared,
n_obs = as.integer(nrow(cs))
)
# =====================================================================
# 6. Calendar-time portfolio (Jaffe-Mandelker / Fama 1998)
# Each trading day: equal-weight firms with >=1 recall within
# +/- 30 trading-day index distance. Regress port ret on market ret.
# =====================================================================
event_firm_idx <- data.frame(firm_id = recalls$firm_id,
eidx = as.integer(date_to_idx[as.character(recalls$date)]),
stringsAsFactors = FALSE)
event_firm_idx <- event_firm_idx[!is.na(event_firm_idx$eidx), ]
market_by_date <- tapply(returns$market_return, as.character(returns$date),
function(z) z[1])
port_rows <- list()
for (g in 0:(n_dates - 1)) {
dstr <- as.character(all_dates[g + 1])
# firms that have a recall within +/-30 index distance of this day
in_win <- unique(event_firm_idx$firm_id[abs(event_firm_idx$eidx - g) <= 30])
if (length(in_win) == 0) next
rets <- c()
for (fid in in_win) {
ft <- firm_tables[[fid]]
if (!is.null(ft) && dstr %in% rownames(ft)) rets <- c(rets, ft[dstr, "return"])
}
if (length(rets) == 0) next
port_rows[[length(port_rows) + 1L]] <- data.frame(
date = all_dates[g + 1],
port_ret = mean(rets),
market_return = as.numeric(market_by_date[dstr]),
stringsAsFactors = FALSE
)
}
port_df <- do.call(rbind, port_rows)
port_df <- port_df[order(port_df$date), ]
cal_fit <- lm(port_ret ~ market_return, data = port_df)
cal_sum <- summary(cal_fit)
calendar <- list(
alpha = unname(coef(cal_fit)["(Intercept)"]),
beta = unname(coef(cal_fit)["market_return"]),
t_alpha = unname(cal_sum$coefficients["(Intercept)", "t value"]),
p_value = unname(cal_sum$coefficients["(Intercept)", "Pr(>|t|)"]),
n_days = as.integer(nrow(port_df))
)
# =====================================================================
# 7. Corrado (1989) non-parametric rank test
# Rank raw returns over [-30,+30]; K_i = mean over [-1,+1] of
# standardized rank deviation; z = mean(K)/(sd(K)/sqrt(N)).
# =====================================================================
K_vals <- c()
for (i in seq_len(n_events)) {
fid <- recalls$firm_id[i]
estr <- as.character(recalls$date[i])
if (!(estr %in% names(date_to_idx))) next
eidx <- as.integer(date_to_idx[estr])
ft <- firm_tables[[fid]]
if (is.null(ft)) nextcd /app && python3 << 'EOF' 2>&1
import numpy as np, pandas as pd, json
from scipy import stats
import statsmodels.api as sm
J=json.load(open("outputs/event_study_results.json"))
r=pd.read_csv("data/stock_returns.csv").dropna(subset=["return","market_return"])
r=r[r["return"]>=-10].drop_duplicates(); r["date"]=pd.to_datetime(r["date"]); r["ds"]=r["date"].dt.strftime("%Y-%m-%d")
rec=pd.read_csv("data/recalls.csv").drop_duplicates().sort_values("event_id"); rec["date"]=pd.to_datetime(rec["date"])
fr=pd.read_csv("data/firms.csv").drop_duplicates()
ad=np.sort(r["date"].unique()); ds=[pd.Timestamp(d).strftime("%Y-%m-%d") for d in ad]; d2i={ds[i]:i for i in range(len(ds))}; nD=len(ad)
ft={f:g.set_index("ds") for f,g in r.groupby("firm_id")}
def fit(fid,ei):
ee=ei-31; es=ee-199; ed=set(ds[k] for k in range(es,ee+1)); g=ft[fid]; s=g[g.index.isin(ed)]
x=s["market_return"].values;y=s["return"].values;b1,b0=np.polyfit(x,y,1);res=y-(b0+b1*x);L=len(s)
return dict(a=b0,b=b1,s2=np.sum(res**2)/(L-2),L=L,mrm=x.mean(),ssrm=np.sum((x-x.mean())**2),phat=np.mean(res>0),n=L)
wins={"scar_3day":(-1,1),"scar_2day":(0,1),"scar_11day":(-5,5)}
rows=[]
for _,row in rec.iterrows():
fid=row["firm_id"];ts=row["date"].strftime("%Y-%m-%d")
if ts not in d2i:continue
ei=d2i[ts];m=fit(fid,ei);g=ft[fid];o={"event_id":int(row["event_id"]),"firm_id":fid};ok3=True;ar3=[]
for wn,(lo,hi) in wins.items():
a=0;v=0;ok=True
for off in range(lo,hi+1):
tg=ds[ei+off]
if tg not in g.index:ok=False;break
rm=g.loc[tg,"market_return"];rt=g.loc[tg,"return"];ar=rt-(m["a"]+m["b"]*rm);a+=ar
v+=m["s2"]*(1+1/m["L"]+(rm-m["mrm"])**2/m["ssrm"])
if wn=="scar_3day":ar3.append(ar)
o[wn]=a/np.sqrt(v) if ok else np.nan
if wn=="scar_3day" and not ok:ok3=False
if not ok3 or np.isnan(o["scar_3day"]):continue
o["phat"]=m["phat"];o["ar3"]=ar3;rows.append(o)
ev=pd.DataFrame(rows);N=len(ev)
# rho
resid={}
for f,g in ft.items():
if len(g)<30:continue
x=g["market_return"].values;y=g["return"].values;b1,b0=np.polyfit(x,y,1);resid[f]=pd.Series(y-(b0+b1*x),index=g.index)
fs=list(resid);pc=[]
for i in range(len(fs)):
for j in range(i+1,len(fs)):
a=resid[fs[i]];b=resid[fs[j]];c=a.index.intersection(b.index)
if len(c)<20:continue
pc.append(np.corrcoef(a[c],b[c])[0,1])
rho=np.mean(pc)
def chk(name,got,exp,tol=1e-6):
ok=abs(got-exp)<=tol*max(1,abs(exp)); print(("OK " if ok else "XX ")+f"{name}: json={got} recomputed={exp}")
return ok
allok=[]
allok.append(chk("rho_bar",J["rho_bar"],rho))
for wn,key in [("scar_3day","window_3day"),("scar_2day","window_2day"),("scar_11day","window_11day")]:
v=ev[wn].dropna().values;n=len(v)
allok.append(chk(key+".mean",J["aggregate_scars"][key]["mean_scar"],v.mean()))
allok.append(chk(key+".patell",J["aggregate_scars"][key]["patell_z"],v.sum()/np.sqrt(n)))
allok.append(chk(key+".bmp",J["aggregate_scars"][key]["bmp_t"],v.mean()/(v.std(ddof=1)/np.sqrt(n))))
allok.append(chk(key+".kp",J["aggregate_scars"][key]["kp_t"],v.mean()/np.sqrt(np.var(v,ddof=1)*(1+(n-1)*rho)/n)))
allok.append(chk(key+".pctneg",J["aggregate_scars"][key]["pct_negative"],float(np.mean(v<0))))
# cross section
cs=ev[["event_id","firm_id","scar_3day"]].merge(rec[["event_id","units_recalled","media_articles","severity_score"]],on="event_id").merge(fr[["firm_id","market_cap_mm","idio_vol"]],on="firm_id").sort_values("event_id")
X=sm.add_constant(np.c_[np.log(cs["units_recalled"]),cs["media_articles"],cs["severity_score"],np.log(cs["market_cap_mm"])]);y=cs["scar_3day"].values
L=int(np.floor(4*(N/100)**(2/9)));m=sm.OLS(y,X).fit();h=m.get_robustcov_results(cov_type="HAC",maxlags=L,use_correction=False)
H=J["cross_section_hac"]
for i,nm in enumerate(["intercept","log_units","media","severity","log_mcap"]):
allok.append(chk("hac."+nm,H[nm],m.params[i]));allok.append(chk("hac.se_"+nm,H["se_"+nm],h.bse[i]))
allok.append(chk("
…[truncated]OK rho_bar: json=0.0056821459330568 recomputed=0.005682145933056799 OK window_3day.mean: json=-1.0116304539979 recomputed=-1.0116304539978969 OK window_3day.patell: json=-5.98488647686484 recomputed=-5.984886476864844 OK window_3day.bmp: json=-4.99213776186089 recomputed=-4.992137761860886 OK window_3day.kp: json=-4.57015805533409 recomputed=-4.570158055334095 OK window_3day.pctneg: json=0.828571428571429 recomputed=0.8285714285714286 OK window_2day.mean: json=-0.862242761104183 recomputed=-0.8622427611041829 OK window_2day.patell: json=-5.10109696709245 recomputed=-5.101096967092448 OK window_2day.bmp: json=-4.50996935631404 recomputed=-4.509969356314038 OK window_2day.kp: json=-4.12874679471694 recomputed=-4.128746794716938 OK window_2day.pctneg: json=0.8 recomputed=0.8 OK window_11day.mean: json=-0.450723180620522 recomputed=-0.4507231806205222 OK window_11day.patell: json=-2.66651429664343 recomputed=-2.666514296643428 OK window_11day.bmp: json=-2.17487494689865 recomputed=-2.1748749468986492 OK window_11day.kp: json=-1.99103525023878 recomputed=-1.991035250238784 OK window_11day.pctneg: json=0.628571428571429 recomputed=0.6285714285714286 OK hac.intercept: json=6.87041540570889 recomputed=6.870415405709055 OK hac.se_intercept: json=1.989976224197 recomputed=1.9899762241969783 OK hac.log_units: json=-0.277206000865995 recomputed=-0.27720600086599434 OK hac.se_log_units: json=0.143768062566243 recomputed=0.14376806256623817 OK hac.media: json=-0.140035223934644 recomputed=-0.14003522393464252 OK hac.se_media: json=0.045344608521185 recomputed=0.045344608521184576 OK hac.severity: json=-0.0493558669041088 recomputed=-0.04935586690410798 OK hac.se_severity: json=0.0796849034968555 recomputed=0.07968490349685442 OK hac.log_mcap: json=-0.349664887829838 recomputed=-0.34966488782986666 OK hac.se_log_mcap: json=0.109201320305086 recomputed=0.10920132030508992 OK hac.r2: json=0.365806981865842 recomputed=0.36580698186584226 OK hac.lag: json=3 recomputed=3 OK wls.intercept: json=6.18141370419048 recomputed=6.181413704190462 OK wls.se_intercept: json=2.37397582785115 recomputed=2.373975827851146 OK wls.log_units: json=-0.245301252974036 recomputed=-0.24530125297403496 OK wls.se_log_units: json=0.119944251522637 recomputed=0.11994425152263676 OK wls.media: json=-0.142987453304122 recomputed=-0.14298745330412244 OK wls.se_media: json=0.0630179237202238 recomputed=0.06301792372022384 OK wls.severity: json=-0.0688815591352763 recomputed=-0.06888155913527738 OK wls.se_severity: json=0.0742443322545675 recomputed=0.07424433225456747 OK wls.log_mcap: json=-0.301826107135229 recomputed=-0.3018261071352276 OK wls.se_log_mcap: json=0.192046874929762 recomputed=0.19204687492976175 OK wls.r2: json=0.307908153604217 recomputed=0.3079081536042162 OK schipper.wald: json=25.6288672616647 recomputed=25.628867261664794 OK schipper.p: json=1.14055552275834e-05 recomputed=1.1405555227583052e-05 OK cowan.z: json=0.253546652207721 recomputed=0.25354665220772143 OK cowan.p: json=0.799845815310467 recomputed=0.7998458153104671 OK cowan.pbar: json=0.492859296482412 recomputed=0.49285929648241206 OK cowan.fp: json=0.514285714285714 recomputed=0.5142857142857142 OK corrado.z: json=-2.42651806653002 recomputed=-2.426518066530017 OK corrado.p: json=0.015244491285374 recomputed=0.015244491285373996 OK cal.alpha: json=-0.000857770921066497 recomputed=-0.0008577709210664968 OK cal.beta: json=1.03320264843221 recomputed=1.0332026484322134 OK cal.t: json=-1.34565994573033 recomputed=-1.345659945730334 OK cal.p: json=0.179748847934933 recomputed=0.17974884793493284 OK cal.nd: json=230 recomputed=230 ==> ALL MATCH
cd /app && rm -f /tmp/run1.* && echo "=== Forbidden package check ===" && grep -iE "eventstudies|estudy2|EventStudy|RcppEventStudy" analysis.R || echo "none found (good)" && echo && echo "=== library() calls ===" && grep -n "library(" analysis.R && echo && echo "=== Final clean run ===" && rm -rf outputs && Rscript analysis.R && echo "exit=$?" && echo && echo "=== Deliverables ===" && ls -la outputs/ && echo && echo "memo words: $(wc -w < outputs/research_memo.md) | csv rows: $(($(wc -l < outputs/event_cars.csv)-1))"=== Forbidden package check === none found (good) === library() calls === 2: library(jsonlite) 3: library(ggplot2) === Final clean run === Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3 exit=0 === Deliverables === total 161 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:10 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:10 .. -rw-r--r-- 1 h2tagent h2tagent 48651 Jun 20 18:10 car_timeline.png -rw-r--r-- 1 h2tagent h2tagent 3618 Jun 20 18:10 event_cars.csv -rw-r--r-- 1 h2tagent h2tagent 2273 Jun 20 18:10 event_study_results.json -rw-r--r-- 1 h2tagent h2tagent 92721 Jun 20 18:10 portfolio_comparison.png -rw-r--r-- 1 h2tagent h2tagent 7418 Jun 20 18:10 research_memo.md memo words: 1149 | csv rows: 35
No file changes, the agent left the workspace untouched.
# SyncVals normalized verifier summary # task: product-recall-stock-price-event # attempt: 2 status: completed raw_score: 0.0 # reward rule: binary: score 0.0 != 1.0 # => reward = 0.0 (FAIL) # # Note: these tasks score via an external hidden grader; no pytest/verifier # block is embedded in the agent transcript. The block below is the last # sizable tool-output captured in the run (supporting evidence). ---------------------------------------------------------------------- === Forbidden package check === none found (good) === library() calls === 2: library(jsonlite) 3: library(ggplot2) === Final clean run === Done. n_valid_events = 35 | rho_bar = 0.0057 | hac_lag = 3 exit=0 === Deliverables === total 161 drwxr-sr-x 2 h2tagent h2tagent 4096 Jun 20 18:10 . drwxrws--- 4 h2tagent h2tagent 4096 Jun 20 18:10 .. -rw-r--r-- 1 h2tagent h2tagent 48651 Jun 20 18:10 car_timeline.png -rw-r--r-- 1 h2tagent h2tagent 3618 Jun 20 18:10 event_cars.csv -rw-r--r-- 1 h2tagent h2tagent 2273 Jun 20 18:10 event_study_results.json -rw-r--r-- 1 h2tagent h2tagent 92721 Jun 20 18:10 portfolio_comparison.png -rw-r--r-- 1 h2tagent h2tagent 7418 Jun 20 18:10 research_memo.md memo words: 1149 | csv rows: 35
Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_85046dd0d2d4466a. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_85046dd0d2d4466a · verifier authoritative; classifier explanatory.