*! ============================================================================
*! event_study_replication.do
*! Goodwin — "Read the Fine Print: What ExxonMobil's Proxy Actually Says
*!           About Texas Redomiciliation" (CLS Blue Sky Blog, April 2026)
*!
*! Replicates the five-specification market battery described in the Market
*! Reaction section and the placebo / ADH ratio tests underlying Figure 1.
*!
*! Specification 1 : Synthetic control (XOM vs 21 energy-sector donors,
*!                   non-negative simplex weights summing to 1)
*! Specification 2 : Market model (XOM vs SPY)
*! Specification 3 : Oil-augmented market model (XOM vs SPY + BNO)
*! Specification 4 : Matched-pair, both market-model adjusted (XOM - CVX)
*! Specification 5 : Raw differential (XOM - CVX)
*!
*! Plus: Patell z-test, Corrado rank test, bootstrap p, placebo permutation
*!       inference (gap-based and post/pre RMSPE ratio, ADH), multi-window
*!       CARs, and a TOST two-one-sided equivalence test at +/- 2 pp.
*!
*! Requires: Stata 15 or later, synth (ssc install synth), optionally
*!           mat2txt (ssc install mat2txt), estout (ssc install estout).
*!
*! Input   : daily_closes.csv — Date column plus Adj Close for
*!           XOM, CVX, COP, EOG, SLB, OXY, MPC, VLO, PSX, DVN, FANG, HAL,
*!           BKR, CTRA, OVV, APA, EQT, WMB, OKE, TRGP, LNG, ET, SPY, BNO.
*!
*! Author  : Shane Goodwin (Stata port auto-generated 2026-04-23)
*! ============================================================================

version 15
clear all
set more off
set seed 42

*-----------------------------------------------------------------------------*
* 0. USER SETTINGS — edit path and event date if needed.
*-----------------------------------------------------------------------------*
global DATA   "daily_closes.csv"        // relative or absolute path
global OUTDIR "."                       // where logs and results go
global EVENT   "10mar2026"              // event date T
global EST_N   240                      // estimation window length (trading days)
global TOST_BAND 0.02                   // equivalence band, +/- 2 percentage points

* 21 energy-sector donors (article spec)
global DONORS "CVX COP EOG SLB OXY MPC VLO PSX DVN FANG HAL BKR CTRA OVV APA EQT WMB OKE TRGP LNG ET"

cap log close
log using "$OUTDIR/event_study_replication.log", replace text

*-----------------------------------------------------------------------------*
* 1. LOAD PRICES, COMPUTE SIMPLE RETURNS
*-----------------------------------------------------------------------------*
import delimited "$DATA", varnames(1) case(upper) clear
gen date = date(DATE, "YMD")
format date %tdNN/DD/CCYY
drop DATE
order date

* Drop columns that are entirely missing (HES, KMI, MRO in some pulls)
foreach v of varlist _all {
    if "`v'" == "date" continue
    qui count if missing(`v')
    if r(N) == _N drop `v'
}

tsset date

* Simple returns for all tickers
ds, has(type numeric)
local tickers = r(varlist)
local tickers : list tickers - date
foreach t of local tickers {
    gen double r_`t' = (`t' - L.`t') / L.`t'
}
drop `tickers'

rename (r_*) (*)
drop if missing(XOM)                    // drop the first (no lag) row

* Event-date indexing
gen long t = _n
summarize t if date == date("$EVENT", "DMY"), meanonly
local EIDX = r(mean)
if missing(`EIDX') {
    di as err "Event date $EVENT not found in data. Check daily_closes.csv."
    exit 498
}
di as txt "Event date $EVENT is row `EIDX' of " _N "."

* Estimation window: T-EST_N ... T-1 (EST_N observations)
local EST_LO = `EIDX' - $EST_N
local EST_HI = `EIDX' - 1
gen byte est_win = inrange(t, `EST_LO', `EST_HI')

di as txt "Estimation window: rows `EST_LO'–`EST_HI' ($EST_N days)."

*=============================================================================*
* 2. SPECIFICATION 2 — MARKET MODEL (XOM on SPY)
*=============================================================================*
di _n as res "==== 2. Market model: XOM = alpha + beta*SPY ===="

reg XOM SPY if est_win, robust
scalar a_mm     = _b[_cons]
scalar b_mm     = _b[SPY]
predict double r_mm_hat, xb
predict double r_mm_res if est_win, resid
qui summ r_mm_res if est_win
scalar sigma_mm = r(sd) * sqrt((r(N)-1)/(r(N)-2))   // ddof=2
scalar r2_mm    = e(r2)

* Day-0 AR
qui summ XOM if t == `EIDX', meanonly
scalar r_xom0 = r(mean)
qui summ SPY if t == `EIDX', meanonly
scalar r_spy0 = r(mean)
scalar ar_mm0     = r_xom0 - (a_mm + b_mm * r_spy0)
scalar patell_z   = ar_mm0 / sigma_mm
scalar patell_p   = 2 * (1 - normal(abs(patell_z)))

di as txt "alpha="    %9.6f a_mm     "  beta="  %7.4f b_mm
di as txt "sigma="    %9.6f sigma_mm "  R^2="   %7.4f r2_mm
di as txt "AR_0="     %9.6f ar_mm0   "  Patell z=" %7.4f patell_z "  p=" %7.4f patell_p

*=============================================================================*
* 3. SPECIFICATION 3 — OIL-AUGMENTED (XOM on SPY + BNO)
*=============================================================================*
di _n as res "==== 3. Oil-augmented: XOM = alpha + b1*SPY + b2*BNO ===="

reg XOM SPY BNO if est_win, robust
scalar a_o   = _b[_cons]
scalar b_om  = _b[SPY]
scalar b_oo  = _b[BNO]
predict double r_oil_res if est_win, resid
qui summ r_oil_res if est_win
scalar sigma_o = r(sd) * sqrt((r(N)-1)/(r(N)-3))   // ddof=3
scalar r2_o    = e(r2)

qui summ BNO if t == `EIDX', meanonly
scalar r_bno0 = r(mean)
scalar ar_oil0   = r_xom0 - (a_o + b_om * r_spy0 + b_oo * r_bno0)
scalar oil_z     = ar_oil0 / sigma_o
scalar oil_p_pat = 2 * (1 - normal(abs(oil_z)))

di as txt "alpha="    %9.6f a_o   "  beta_mkt=" %7.4f b_om "  beta_oil=" %7.4f b_oo
di as txt "sigma="    %9.6f sigma_o "  R^2="   %7.4f r2_o
di as txt "AR_0="     %9.6f ar_oil0 "  Patell z=" %7.4f oil_z "  p=" %7.4f oil_p_pat

* ---- Corrado rank test (oil-augmented) ----
cap drop r_oil_all
gen double r_oil_all = XOM - (a_o + b_om*SPY + b_oo*BNO) if est_win | t == `EIDX'
egen double rk = rank(r_oil_all) if !missing(r_oil_all), unique
qui count if !missing(r_oil_all)
scalar nT = r(N)
scalar mean_rk = (nT + 1) / 2

qui summ r_oil_all if !missing(r_oil_all)
cap drop rkdev
gen double rkdev = rk - mean_rk if !missing(r_oil_all)
qui summ rkdev
scalar sd_rk = sqrt(r(sd)^2 * (r(N)-1)/r(N))    // population sd

qui summ rk if t == `EIDX'
scalar rk_event = r(mean)
scalar corrado_z = (rk_event - mean_rk) / sd_rk
scalar corrado_p = 2 * (1 - normal(abs(corrado_z)))
di as txt "Corrado: event rank=" %5.0f rk_event " of " %5.0f nT ///
          "  z=" %7.4f corrado_z "  p=" %7.4f corrado_p

* ---- Bootstrap p-value for oil-augmented AR_0 ----
preserve
    keep if est_win
    keep r_oil_res
    drop if missing(r_oil_res)
    tempfile bootres
    save `bootres'
    local n_boot = 10000
    scalar hits = 0
    forvalues b = 1/`n_boot' {
        use `bootres', clear
        bsample 1
        qui summ r_oil_res
        if abs(r(mean)) >= abs(scalar(ar_oil0)) scalar hits = scalar(hits) + 1
    }
    scalar boot_p = hits / `n_boot'
restore
di as txt "Bootstrap p (oil-augmented, B=" %5.0f `n_boot' "): " %7.4f boot_p

*=============================================================================*
* 4. SPECIFICATION 4 — MATCHED PAIR (XOM - CVX, both market-model adjusted)
*=============================================================================*
di _n as res "==== 4. Matched pair: MM-adjusted XOM - CVX ===="

reg CVX SPY if est_win
scalar a_cvx = _b[_cons]
scalar b_cvx = _b[SPY]
predict double r_cvx_res if est_win, resid
qui summ r_cvx_res if est_win
scalar sigma_cvx = r(sd) * sqrt((r(N)-1)/(r(N)-2))
scalar r2_cvx    = e(r2)

qui summ CVX if t == `EIDX', meanonly
scalar r_cvx0 = r(mean)
scalar ar_cvx0 = r_cvx0 - (a_cvx + b_cvx * r_spy0)
scalar delta0  = ar_mm0 - ar_cvx0

* SE on delta = sqrt(sigma_XOM^2 + sigma_CVX^2 - 2*cov(res))
qui corr r_mm_res r_cvx_res if est_win, cov
scalar cov_xc = r(cov_12)
scalar se_delta = sqrt(sigma_mm^2 + sigma_cvx^2 - 2*cov_xc)
scalar t_pair   = delta0 / se_delta
scalar p_pair   = 2 * (1 - normal(abs(t_pair)))

di as txt "AR_XOM=" %8.5f ar_mm0 "  AR_CVX=" %8.5f ar_cvx0 ///
          "  Delta=" %8.5f delta0
di as txt "SE_delta=" %8.5f se_delta "  t=" %7.4f t_pair "  p=" %7.4f p_pair

*=============================================================================*
* 5. SPECIFICATION 5 — RAW DIFFERENTIAL (XOM - CVX)
*=============================================================================*
di _n as res "==== 5. Raw differential XOM - CVX on Day 0 ===="
scalar raw_diff = r_xom0 - r_cvx0
di as txt "XOM=" %8.5f r_xom0 "  CVX=" %8.5f r_cvx0 "  raw diff=" %8.5f raw_diff

*=============================================================================*
* 6. SPECIFICATION 1 — SYNTHETIC CONTROL (simplex weights, energy donors)
*    Implemented in Mata: min ||y - X*w||^2  s.t. w>=0, sum(w)=1
*=============================================================================*
di _n as res "==== 1. Synthetic control: XOM vs 21 energy donors ===="

mata:
// ---- projection onto the simplex (Duchi, Shalev-Shwartz, Singer, Chandra 2008) ----
real colvector proj_simplex(real colvector v) {
    real colvector u, cssv, rho_idx, cond
    real scalar n, rho, tau
    n = rows(v)
    u = sort(v, -1)
    cssv = runningsum(u)
    cond = u :- (cssv :- 1) :/ (1..n)'
    rho  = colsum(cond :> 0)
    tau  = (cssv[rho] - 1) / rho
    return(editmissing((v :- tau) :* ((v :- tau) :> 0), 0))
}

// ---- projected-gradient solver for simplex-constrained LS ----
real colvector sc_weights(real colvector y, real matrix X) {
    real colvector w, w_new, g
    real scalar n_donor, step, tol, iter, change
    n_donor = cols(X)
    w = J(n_donor, 1, 1/n_donor)
    step = 2 / (2 * trace(X'X))
    tol = 1e-14
    for (iter=1; iter<=200000; iter++) {
        g = 2 * X' * (X*w - y)
        w_new = proj_simplex(w - step*g)
        change = max(abs(w_new - w))
        w = w_new
        if (change < tol) break
    }
    return(w)
}

// ---- pull data from Stata ----
st_view(Y=., ., "XOM", "est_win")
donor_names = tokens(st_global("DONORS"))
st_view(X=., ., donor_names, "est_win")

w = sc_weights(Y, X)
st_matrix("W_SC", w)
st_local("wsum", strofreal(sum(w)))

// pre-period fit
y_hat = X * w
resid = Y - y_hat
ss_res = sum(resid:^2)
ss_tot = sum((Y :- mean(Y)):^2)
pre_rmspe = sqrt(mean(resid:^2))
corr_pre  = correlation((Y, y_hat))[1,2]
r2_pre    = 1 - ss_res/ss_tot

st_numscalar("pre_rmspe", pre_rmspe)
st_numscalar("corr_pre",  corr_pre)
st_numscalar("r2_pre",    r2_pre)
end

matrix colnames W_SC = w
matrix rownames W_SC = $DONORS

di as txt "Simplex weight sum = `wsum'   (should be ~1)"
di as txt "Pre-period RMSPE=" %9.6f pre_rmspe "  corr=" %7.4f corr_pre "  R^2=" %7.4f r2_pre

di _n as txt "Top donor weights:"
mat li W_SC, format(%7.4f)

* Day-0 gap
mata:
st_view(Xrow=., ., donor_names, "t == `EIDX'")
w = st_matrix("W_SC")
synth0 = Xrow * w
st_numscalar("synth_r0", synth0)
end
scalar day0_gap = r_xom0 - synth_r0
di as txt "XOM Day-0 return=" %9.6f r_xom0 "   synthetic=" %9.6f synth_r0 ///
          "   gap=" %9.6f day0_gap

* Post-event [-1,+5] RMSPE (synthetic)
local POST_LO = `EIDX' - 1
local POST_HI = `EIDX' + 5
mata:
st_view(Ypost=., ., "XOM",       "inrange(t, `POST_LO', `POST_HI')")
st_view(Xpost=., ., donor_names, "inrange(t, `POST_LO', `POST_HI')")
w = st_matrix("W_SC")
post_resid = Ypost - Xpost * w
st_numscalar("post_rmspe", sqrt(mean(post_resid:^2)))
end
scalar rmspe_ratio = post_rmspe / pre_rmspe
di as txt "Post [-1,+5] RMSPE=" %9.6f post_rmspe ///
          "   post/pre ratio=" %7.4f rmspe_ratio

*=============================================================================*
* 7. PLACEBO / PERMUTATION INFERENCE (gap-based AND ADH ratio)
*    Re-fit synthetic control treating each donor as the "pseudo-treated"
*    and record day-0 gap and post/pre RMSPE ratio.
*=============================================================================*
di _n as res "==== 7. Placebo permutation inference (22-unit set) ===="

tempname PGAPS PRATIOS
matrix `PGAPS'   = J(22, 2, .)
matrix `PRATIOS' = J(22, 2, .)

* XOM is first row; donors fill rows 2..22
matrix `PGAPS'[1,1]   = day0_gap
matrix `PRATIOS'[1,1] = rmspe_ratio

local i = 2
foreach D of global DONORS {
    local OTHER : list global(DONORS) - D
    mata:
    Dname = st_local("D")
    other = tokens(st_local("OTHER"))
    st_view(Yp=., ., Dname, "est_win")
    st_view(Xp=., ., other, "est_win")
    wp = sc_weights(Yp, Xp)

    // day-0 gap
    st_view(Yp0=., ., Dname, "t == `EIDX'")
    st_view(Xp0=., ., other, "t == `EIDX'")
    gap = Yp0 - Xp0 * wp

    // pre/post RMSPE
    pre_r = sqrt(mean((Yp - Xp*wp):^2))
    st_view(Ypost2=., ., Dname, "inrange(t, `POST_LO', `POST_HI')")
    st_view(Xpost2=., ., other, "inrange(t, `POST_LO', `POST_HI')")
    post_r = sqrt(mean((Ypost2 - Xpost2*wp):^2))
    ratio = post_r / pre_r

    st_numscalar("g_`D'", gap)
    st_numscalar("r_`D'", ratio)
    end

    matrix `PGAPS'[`i',1]   = scalar(g_`D')
    matrix `PRATIOS'[`i',1] = scalar(r_`D')
    local ++i
}

* Rank by absolute gap (two-sided)
preserve
    clear
    svmat double `PGAPS', names(col)
    rename c1 gap
    gen absgap = abs(gap)
    gen str8 unit = ""
    replace unit = "XOM" in 1
    local k = 2
    foreach D of global DONORS {
        replace unit = "`D'" in `k'
        local ++k
    }
    gsort -absgap
    gen rank = _n
    list rank unit gap, sepby(rank) noobs
    qui summ rank if unit == "XOM"
    scalar xom_rank_gap = r(mean)
    scalar p_gap = xom_rank_gap / _N
restore
di as txt "Placebo p (two-sided, gap-based): " %7.4f p_gap ///
          "   XOM |AR| rank=" %3.0f xom_rank_gap " of 22"

* Rank by post/pre RMSPE ratio (ADH)
preserve
    clear
    svmat double `PRATIOS', names(col)
    rename c1 ratio
    gen str8 unit = ""
    replace unit = "XOM" in 1
    local k = 2
    foreach D of global DONORS {
        replace unit = "`D'" in `k'
        local ++k
    }
    gsort -ratio
    gen rank = _n
    list rank unit ratio, sepby(rank) noobs
    qui summ rank if unit == "XOM"
    scalar xom_rank_ratio = r(mean)
    scalar p_ratio = xom_rank_ratio / _N
restore
di as txt "Placebo p (ADH post/pre RMSPE ratio): " %7.4f p_ratio ///
          "   XOM ratio rank=" %3.0f xom_rank_ratio " of 22"

*=============================================================================*
* 8. MULTI-WINDOW CARs (synthetic, market model, oil)
*=============================================================================*
di _n as res "==== 8. Multi-window CARs ===="

foreach win in "0 0" "-1 1" "-1 3" "-1 5" "-2 5" {
    tokenize "`win'"
    local LO = `1'
    local HI = `2'
    local ALO = `EIDX' + `LO'
    local AHI = `EIDX' + `HI'

    mata:
    st_view(Yw=., ., "XOM",       "inrange(t, `ALO', `AHI')")
    st_view(Xw=., ., donor_names, "inrange(t, `ALO', `AHI')")
    w = st_matrix("W_SC")
    st_numscalar("syn_car", sum(Yw - Xw*w))
    end

    qui gen double _tmp_mm  = XOM - (a_mm + b_mm*SPY)                   if inrange(t,`ALO',`AHI')
    qui gen double _tmp_oil = XOM - (a_o  + b_om*SPY + b_oo*BNO)       if inrange(t,`ALO',`AHI')
    qui summ _tmp_mm
    scalar mm_car = r(sum)
    qui summ _tmp_oil
    scalar oil_car = r(sum)
    drop _tmp_mm _tmp_oil

    di as txt "[`LO',`HI']  Synth=" %+7.3f syn_car*100 "%  " ///
              "MM=" %+7.3f mm_car*100 "%  " ///
              "Oil=" %+7.3f oil_car*100 "%"
}

*=============================================================================*
* 9. TOST EQUIVALENCE TEST at +/- $TOST_BAND (two one-sided tests)
*    Applied to market-model AR_0.  Using standardized AR ~ N(0,1) under null.
*=============================================================================*
di _n as res "==== 9. TOST equivalence test (market-model AR_0) ===="
scalar lo  = -$TOST_BAND
scalar hi  =  $TOST_BAND
scalar t_lo = (ar_mm0 - lo) / sigma_mm
scalar t_hi = (hi - ar_mm0) / sigma_mm
scalar p_lo = 1 - normal(t_lo)
scalar p_hi = 1 - normal(t_hi)
scalar p_tost = max(p_lo, p_hi)
di as txt "Band=[" %6.3f lo ", " %6.3f hi "]   AR_0=" %8.5f ar_mm0 ///
          "   sigma=" %8.5f sigma_mm
di as txt "t_lo=" %7.4f t_lo " (p=" %6.4f p_lo ")   " ///
          "t_hi=" %7.4f t_hi " (p=" %6.4f p_hi ")"
di as txt "TOST p (max of two one-sided) = " %7.4f p_tost ///
          "   reject equivalence = " cond(p_tost<0.05, "NO (cannot reject)", "YES — governance effect is within +/- 2pp")

*=============================================================================*
* 10. WRITE SUMMARY RESULTS TABLE
*=============================================================================*
di _n as res "==== 10. Summary ===="
tempname R
matrix `R' = J(5, 4, .)
matrix rownames `R' = "1_SyntheticControl" "2_MarketModel" "3_OilAugmented" "4_MatchedPair" "5_RawDifferential"
matrix colnames `R' = "AR_day0" "test_stat" "p_value" "note"

matrix `R'[1,1] = day0_gap
matrix `R'[1,2] = xom_rank_gap
matrix `R'[1,3] = p_gap

matrix `R'[2,1] = ar_mm0
matrix `R'[2,2] = patell_z
matrix `R'[2,3] = patell_p

matrix `R'[3,1] = ar_oil0
matrix `R'[3,2] = oil_z
matrix `R'[3,3] = oil_p_pat

matrix `R'[4,1] = delta0
matrix `R'[4,2] = t_pair
matrix `R'[4,3] = p_pair

matrix `R'[5,1] = raw_diff
matrix `R'[5,2] = .
matrix `R'[5,3] = .

matrix list `R', format(%9.5f) title("Event-study battery — XOM, 2026-03-10")

di _n as txt "Placebo-ADH: XOM post/pre RMSPE ratio rank = " %3.0f xom_rank_ratio " of 22, p=" %6.4f p_ratio
di as txt    "TOST at +/- 2pp: p = " %6.4f p_tost

log close
*! END of event_study_replication.do