"""
=============================================================================
news_ssi_pipeline.py
=============================================================================
Full pipeline for the news-based Shipping Stress Index paper.
WHAT THIS DOES
──────────────
Step 1 → Query GDELT for shipping-disruption news (free, no API key needed)
Step 2 → Build a tone-weighted, smoothed News-SSI
Step 3 → Download macro data from FRED (free API key required)
Step 4 → Merge into a clean monthly panel
Step 5 → Run VAR + IRF + FEVD
Step 6 → Run event studies around known disruption episodes
Step 7 → Export main_dataset_news.csv (ready for LaTeX tables)
INSTALL DEPENDENCIES FIRST
───────────────────────────
pip install pandas numpy requests statsmodels matplotlib fredapi yfinance
FRED API KEY (free, 30 seconds to get)
https://fred.stlouisfed.org/docs/api/api_key.html
Then: export FRED_API_KEY=your_key_here
=============================================================================
"""
import os, time, json, warnings, requests
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from datetime import datetime, timedelta
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller, kpss
from scipy import stats
warnings.filterwarnings("ignore")
# ── Config ───────────────────────────────────────────────────────────────────
FRED_API_KEY = os.getenv("FRED_API_KEY", "YOUR_API_KEY")
START_DATE = "2015-01-01"
END_DATE = datetime.today().strftime("%Y-%m-%d")
OUTPUT_DIR = "figures_news"
os.makedirs(OUTPUT_DIR, exist_ok=True)
# ── Style ────────────────────────────────────────────────────────────────────
plt.rcParams.update({
"font.family": "serif", "font.serif": ["Times New Roman", "DejaVu Serif"],
"axes.labelsize": 11, "axes.titlesize": 12, "xtick.labelsize": 9,
"ytick.labelsize": 9, "legend.fontsize": 9, "figure.dpi": 150,
"axes.spines.top": False, "axes.spines.right": False,
"axes.grid": True, "grid.alpha": 0.3, "grid.linestyle": "--",
})
NAVY, RED, TEAL, AMBER, STEEL = "#1B3A6B","#C0392B","#17A589","#E67E22","#2E75B6"
COLORS = [NAVY, RED, TEAL, AMBER, STEEL]
# ── Known disruption episodes (for event study) ──────────────────────────────
EPISODES = [
("2016-10-01", "Yemen blockade tightens", 2),
("2019-06-13", "Tanker attacks – Gulf of Oman", 3),
("2021-03-23", "Ever Given blocks Suez Canal", 5),
("2023-11-19", "Houthi attacks begin – Red Sea", 4),
("2026-02-28", "US–Iran war – Hormuz closure", 5), # 🔴 NEW EVENT
]
# ── GDELT search terms ────────────────────────────────────────────────────────
GDELT_THEMES = [
"shipping disruption",
"strait of hormuz",
"suez canal blocked",
"red sea tanker",
"bab el mandeb",
"strait of malacca",
"tanker attack",
"maritime chokepoint",
"oil tanker seized",
"houthi ship",
]
VAR_VARS = ["SSI", "Brent", "CPI", "InterestRate", "GDP"]
VAR_LABELS = ["News SSI", "Brent Crude\n(USD/bbl)", "CPI\n(index)",
"Interest Rate\n(%)", "Real GDP\n(log)"]
# =============================================================================
# STEP 1 — GDELT NEWS DATA (free, no key needed)
# =============================================================================
def query_gdelt_monthly(query: str, year_month: str) -> dict:
"""
Query GDELT's free Document API for a single month.
Returns dict with:
articles : number of articles matching the query
avg_tone : average sentiment tone (-100 negative to +100 positive)
GDELT API docs: https://blog.gdeltproject.org/gdelt-doc-2-0-api-debuts/
"""
# GDELT format: YYYYMMDDHHMMSS
ym = pd.Timestamp(year_month)
start = ym.strftime("%Y%m%d%H%M%S")
end = (ym + pd.offsets.MonthEnd(0)).strftime("%Y%m%d%H%M%S")
url = (
"https://api.gdeltproject.org/api/v2/doc/doc"
f"?query={requests.utils.quote(query)}"
f"&mode=artlist&maxrecords=250"
f"&startdatetime={start}&enddatetime={end}"
"&format=json&sourcelang=english"
)
try:
r = requests.get(url, timeout=20)
data = r.json()
articles = data.get("articles", [])
n = len(articles)
# GDELT returns tone in metadata — average across articles
tones = [float(a.get("tone", 0)) for a in articles if "tone" in a]
avg_tone = np.mean(tones) if tones else 0.0
return {"articles": n, "avg_tone": avg_tone}
except Exception as e:
print(f" [GDELT warn] {year_month}: {e}")
return {"articles": 0, "avg_tone": 0.0}
def build_gdelt_ssi(start: str = START_DATE, end: str = END_DATE,
use_cache: bool = True) -> pd.DataFrame:
"""
Aggregate GDELT article counts across all shipping search terms, monthly.
Strategy:
- For each month × each keyword → article count + avg tone
- Sum counts; average tones
- Tone-weighted SSI: penalise positive (non-crisis) articles
- Smooth with 3-month rolling average
- Normalise to [0, 1]
Set use_cache=True to avoid re-querying (slow: ~1 sec per month per term).
Cache file: gdelt_raw_cache.csv
"""
cache_file = "gdelt_raw_cache.csv"
months = pd.date_range(start, end, freq="MS")
if use_cache and os.path.exists(cache_file):
print(f"[✓] Loading GDELT cache from {cache_file}")
raw = pd.read_csv(cache_file, parse_dates=["Date"], index_col="Date")
else:
print(f"[→] Querying GDELT for {len(months)} months × {len(GDELT_THEMES)} terms...")
print(" This takes ~2-3 minutes. Will be cached afterwards.\n")
records = []
for ym in months:
total_n, all_tones = 0, []
for term in GDELT_THEMES:
res = query_gdelt_monthly(term, ym.strftime("%Y-%m"))
total_n += res["articles"]
if res["articles"] > 0:
all_tones.append(res["avg_tone"])
time.sleep(0.5) # GDELT rate limit: be polite
avg_tone = np.mean(all_tones) if all_tones else 0.0
records.append({"Date": ym, "RawCount": total_n, "AvgTone": avg_tone})
print(f" {ym.strftime('%Y-%m')} articles={total_n:>4} tone={avg_tone:+.1f}")
raw = pd.DataFrame(records).set_index("Date")
raw.to_csv(cache_file)
print(f"[✓] GDELT data cached → {cache_file}")
# ── Tone-weighted SSI ────────────────────────────────────────────────────
# Tone is negative for crisis/disruption news (e.g., -3.5)
# Convert: more negative tone → higher stress weight
tone_weight = np.clip(-raw["AvgTone"] / 10.0 + 1.0, 0.5, 2.5)
raw["WeightedCount"] = raw["RawCount"] * tone_weight
# ── Smooth (3-month rolling) ──────────────────────────────────────────────
raw["Smoothed"] = raw["WeightedCount"].rolling(3, min_periods=1).mean()
# ── Normalise to [0, 1] ──────────────────────────────────────────────────
mn, mx = raw["Smoothed"].min(), raw["Smoothed"].max()
raw["SSI"] = (raw["Smoothed"] - mn) / (mx - mn + 1e-9)
print(f"\n[✓] SSI constructed: {len(raw)} months, "
f"mean={raw['SSI'].mean():.3f}, max={raw['SSI'].max():.3f}")
return raw[["RawCount", "AvgTone", "WeightedCount", "SSI"]]
def build_synthetic_ssi(start: str = START_DATE) -> pd.DataFrame:
"""
Synthetic news SSI for development/testing when GDELT is unavailable.
Mirrors known episode spikes.
"""
np.random.seed(42)
idx = pd.date_range(start, END_DATE, freq="MS")
n = len(idx)
raw = np.random.exponential(5, n)
for ep_date, label, severity in EPISODES: # 🔴 label now used for conditional logic
ep = pd.Timestamp(ep_date)
for i, dt in enumerate(idx):
months_away = abs((dt.year - ep.year) * 12 + dt.month - ep.month)
# 🔴 Distinguish Hormuz structural shock vs normal shocks
if "Hormuz" in label or "Iran" in label:
if months_away <= 6:
raw[i] += severity * 25 * np.exp(-months_away * 0.3)
else:
if months_away <= 3:
raw[i] += severity * 20 * np.exp(-months_away * 0.5)
smooth = pd.Series(raw).rolling(3, min_periods=1).mean().values
ssi = (smooth - smooth.min()) / (smooth.max() - smooth.min() + 1e-9)
tone = np.random.normal(-2, 1, n)
for ep_date, label, severity in EPISODES:
ep = pd.Timestamp(ep_date)
for i, dt in enumerate(idx):
if abs((dt.year - ep.year) * 12 + dt.month - ep.month) <= 2:
tone[i] -= severity * 1.5
return pd.DataFrame({
"RawCount": raw, "AvgTone": tone,
"WeightedCount": raw, "SSI": ssi
}, index=idx)
# =============================================================================
# STEP 2 — FRED MACRO DATA (requires free API key)
# =============================================================================
FRED_SERIES = {
"Brent": "DCOILBRENTEU",
"CPI": "CPIAUCSL",
"InterestRate": "FEDFUNDS",
"GDP": "GDPC1",
"IndProd": "INDPRO",
"BDI": "BALDIBDI",
"USD": "DTWEXBGS",
"GPR": "GEOPOLITICALRISK", # Caldara & Iacoviello GPR index
}
def fetch_fred(series_id: str, freq: str = "m") -> pd.Series:
"""Download one FRED series at monthly frequency."""
if FRED_API_KEY == "YOUR_KEY_HERE":
return pd.Series(dtype=float, name=series_id)
url = (
f"https://api.stlouisfed.org/fred/series/observations"
f"?series_id={series_id}&api_key={FRED_API_KEY}&file_type=json"
f"&observation_start={START_DATE}&observation_end={END_DATE}"
f"&frequency={freq}&aggregation_method=avg"
)
try:
obs = requests.get(url, timeout=15).json()["observations"]
s = pd.Series(
{o["date"]: float(o["value"]) if o["value"] != "." else np.nan
for o in obs},
name=series_id
)
s.index = pd.to_datetime(s.index)
print(f" [FRED ✓] {series_id}: {len(s)} obs")
return s
except Exception as e:
print(f" [FRED ✗] {series_id}: {e}")
return pd.Series(dtype=float, name=series_id)
def fetch_all_fred() -> pd.DataFrame:
print("\n── Fetching FRED series ─────────────────────────────────────────")
frames = {name: fetch_fred(sid) for name, sid in FRED_SERIES.items()}
df = pd.DataFrame(frames)
df.index = pd.to_datetime(df.index)
df = df.resample("MS").mean()
# GDP is quarterly — forward-fill to monthly
df["GDP"] = df["GDP"].interpolate(method="linear")
return df
def build_synthetic_macro(idx: pd.DatetimeIndex) -> pd.DataFrame:
"""Synthetic macro series for testing (mimics FRED structure)."""
np.random.seed(7)
n = len(idx)
brent = 70 + np.cumsum(np.random.normal(0, 1.5, n))
cpi = 100 + np.cumsum(np.random.normal(0.22, 0.2, n))
gdp = np.log(2e13) + np.cumsum(np.random.normal(0.003, 0.006, n))
rate = np.clip(2.5 + np.cumsum(np.random.normal(0, 0.05, n)), 0, 6)
for ep_date, label, severity in EPISODES:
ep = pd.Timestamp(ep_date)
for i, dt in enumerate(idx):
lag = (dt.year - ep.year) * 12 + dt.month - ep.month
if 1 <= lag <= 4:
brent[i] += severity * 4
if 2 <= lag <= 6:
cpi[i] += severity * 0.15
return pd.DataFrame({
"Brent": brent.clip(20), "CPI": cpi,
"InterestRate": rate, "GDP": gdp,
"IndProd": 100 + np.cumsum(np.random.normal(0.1, 0.4, n)),
"BDI": 1500 + np.cumsum(np.random.normal(0, 25, n)).clip(-1300),
"USD": 100 + np.cumsum(np.random.normal(0, 0.5, n)),
"GPR": np.random.exponential(80, n),
}, index=idx)
# =============================================================================
# STEP 3 — MERGE AND PREPARE PANEL
# =============================================================================
def build_panel(ssi_df: pd.DataFrame, macro_df: pd.DataFrame,
output: str = "main_dataset_news.csv") -> pd.DataFrame:
panel = ssi_df[["SSI", "RawCount", "AvgTone"]].join(macro_df, how="inner")
panel.index.name = "Date"
# Log-transform price-level series
for col in ["Brent", "CPI", "GDP"]:
if col in panel.columns:
panel[f"ln_{col}"] = np.log(panel[col].clip(lower=0.01))
# 🔴 Optional structural shock dummy (for additional robustness)
panel["HormuzShock"] = 0
panel.loc["2026-02-01":"2026-06-01", "HormuzShock"] = 1
panel.to_csv(output)
print(f"\n[✓] Panel saved → {output} "
f"({len(panel)} rows × {len(panel.columns)} cols)")
return panel
def stationarity_table(panel: pd.DataFrame) -> pd.DataFrame:
"""ADF + KPSS for each VAR variable."""
rows = []
for col in VAR_VARS:
if col not in panel.columns:
continue
s = panel[col].dropna()
adf_s, adf_p, *_ = adfuller(s, autolag="AIC")
kpss_s, kpss_p, *_ = kpss(s, regression="c", nlags="auto")
rows.append({
"Variable": col,
"ADF stat": round(adf_s, 3),
"ADF p": round(adf_p, 3),
"KPSS stat": round(kpss_s, 3),
"KPSS p": round(kpss_p, 3),
"Decision": "Stationary" if adf_p < 0.05 else "Difference",
})
return pd.DataFrame(rows)
def prepare_var_data(panel: pd.DataFrame) -> pd.DataFrame:
stat = panel[VAR_VARS].copy()
for col in ["Brent", "CPI", "GDP"]:
if col in stat.columns:
stat[col] = stat[col].diff()
return stat.dropna()
# =============================================================================
# STEP 4 — VAR ESTIMATION
# =============================================================================
def fit_var(df_stat: pd.DataFrame, maxlags: int = 6):
model = VAR(df_stat[VAR_VARS].dropna())
result = model.fit(maxlags=maxlags, ic="aic")
print(f"\n[VAR] Fitted VAR({result.k_ar}), AIC={result.aic:.2f}")
return result
# =============================================================================
# STEP 5 — EVENT STUDY
# =============================================================================
def run_event_study(panel: pd.DataFrame,
window: int = 6, # months each side
target: str = "Brent") -> pd.DataFrame:
"""
For each disruption episode, compute:
- Average value of `target` in [-window, -1] (pre-period)
- Average value of `target` in [+1, +window] (post-period)
- Abnormal change = post_mean - pre_mean
- t-test p-value
Returns a tidy DataFrame of event-study results.
"""
results = []
series = panel[target].dropna()
for ep_date, label, severity in EPISODES:
ep = pd.Timestamp(ep_date).to_period("M").to_timestamp()
pre = series[(series.index >= ep - pd.DateOffset(months=window)) &
(series.index < ep)]
post = series[(series.index > ep) &
(series.index <= ep + pd.DateOffset(months=window))]
if len(pre) < 2 or len(post) < 2:
continue
t_stat, p_val = stats.ttest_ind(post, pre, equal_var=False)
results.append({
"Event": label,
"Date": ep_date,
"Severity": severity,
"Pre-mean": round(pre.mean(), 2),
"Post-mean": round(post.mean(), 2),
"Abnormal change": round(post.mean() - pre.mean(), 2),
"t-stat": round(t_stat, 3),
"p-value": round(p_val, 3),
"Significant": "***" if p_val < 0.01 else
"**" if p_val < 0.05 else
"*" if p_val < 0.10 else "—",
})
df = pd.DataFrame(results)
print(f"\n── Event Study ({target}) ──────────────────────────────────────")
print(df[["Event", "Pre-mean", "Post-mean", "Abnormal change",
"t-stat", "p-value", "Significant"]].to_string(index=False))
return df
# =============================================================================
# STEP 6 — ALL GRAPH FUNCTIONS
# =============================================================================
def fig_ssi_with_tone(ssi_df: pd.DataFrame,
save: str = None):
"""Figure 1: SSI time series + tone heatmap + episode markers."""
save = save or f"{OUTPUT_DIR}/fig1_ssi_timeseries.pdf"
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 7), sharex=True,
gridspec_kw={"height_ratios": [3, 1]})
fig.suptitle("News-Based Shipping Stress Index (SSI) — 2015–2026",
fontsize=13, fontweight="bold", color=NAVY)
# Top: SSI line
ssi = ssi_df["SSI"]
ax1.plot(ssi.index, ssi.values, color=NAVY, lw=2, zorder=3)
ax1.fill_between(ssi.index, 0, ssi.values, alpha=0.15, color=NAVY)
ax1.axhline(ssi.mean(), color="#777777", lw=0.9, ls="--",
label=f"Mean = {ssi.mean():.3f}")
for ep_date, label, sev in EPISODES:
ax1.axvline(pd.Timestamp(ep_date), color=RED, lw=1.4, ls=":",
alpha=0.8)
ax1.text(pd.Timestamp(ep_date), 0.96, label[:16],
rotation=90, va="top", ha="right", fontsize=7.5, color=RED,
transform=ax1.get_xaxis_transform())
ax1.set_ylabel("SSI (0 = normal, 1 = shutdown)")
ax1.set_ylim(0, 1.05)
ax1.yaxis.set_major_formatter(mticker.PercentFormatter(xmax=1))
ax1.legend(loc="upper left")
# Bottom: news tone (negative = crisis)
if "AvgTone" in ssi_df.columns:
tone = ssi_df["AvgTone"]
ax2.bar(tone.index, tone.values, color=np.where(tone < 0, RED, TEAL),
alpha=0.7, width=25)
ax2.axhline(0, color="black", lw=0.6)
ax2.set_ylabel("Avg Tone", fontsize=9)
ax2.set_ylim(tone.min() * 1.2, tone.max() * 1.2)
plt.tight_layout()
plt.savefig(save, bbox_inches="tight")
plt.close()
print(f"[✓] {save}")
def fig_raw_article_counts(ssi_df: pd.DataFrame,
save: str = None):
"""Figure 2: Raw article counts by month (bar chart)."""
save = save or f"{OUTPUT_DIR}/fig2_article_counts.pdf"
fig, ax = plt.subplots(figsize=(12, 4))
raw = ssi_df["RawCount"]
ax.bar(raw.index, raw.values, color=STEEL, alpha=0.7, width=25,
label="Raw article count")
ax.plot(raw.index, ssi_df["WeightedCount"] if "WeightedCount" in ssi_df.columns
else raw, color=NAVY, lw=1.5, label="Tone-weighted count")
for ep_date, label, _ in EPISODES:
ax.axvline(pd.Timestamp(ep_date), color=RED, lw=1.2, ls="--")
ax.text(pd.Timestamp(ep_date), raw.max() * 0.92, label[:14],
rotation=90, va="top", fontsize=7.5, color=RED,
ha="right")
ax.set_title("GDELT Monthly Article Counts — Shipping Disruption Keywords",
fontsize=12, fontweight="bold", color=NAVY)
ax.set_ylabel("Articles per month")
ax.legend()
plt.tight_layout()
plt.savefig(save, bbox_inches="tight")
plt.close()
print(f"[✓] {save}")
def fig_event_window(panel: pd.DataFrame,
target: str = "Brent",
window: int = 6,
save: str = None):
"""
Figure 3: Multi-panel event window plot.
Each sub-panel shows the target variable in [-window, +window] months
around each disruption, indexed to 100 at event date.
"""
save = save or f"{OUTPUT_DIR}/fig3_event_window_{target.lower()}.pdf"
episodes_with_data = [(d, l, s) for d, l, s in EPISODES
if pd.Timestamp(d) in panel.index or True]
n_ep = len(episodes_with_data)
fig, axes = plt.subplots(1, n_ep, figsize=(4 * n_ep, 4.5), sharey=False)
if n_ep == 1:
axes = [axes]
fig.suptitle(
f"Event Study: {target} around Major Shipping Disruptions\n"
f"(Indexed to 100 at event date; window = ±{window} months)",
fontsize=12, fontweight="bold", color=NAVY
)
series = panel[target].dropna()
lags = np.arange(-window, window + 1)
for ax, (ep_date, label, sev) in zip(axes, episodes_with_data):
ep = pd.Timestamp(ep_date).to_period("M").to_timestamp()
vals = []
for lag in lags:
dt = ep + pd.DateOffset(months=lag)
dt = dt.to_period("M").to_timestamp()
v = series.get(dt, np.nan)
vals.append(v)
base = vals[window] if not np.isnan(vals[window]) else np.nanmean(vals)
indexed = [100 * v / base if not np.isnan(v) and base != 0
else np.nan for v in vals]
pre_idx = [i for i, l in enumerate(lags) if l < 0]
post_idx = [i for i, l in enumerate(lags) if l >= 0]
pre_v = [indexed[i] for i in pre_idx]
post_v = [indexed[i] for i in post_idx]
ax.plot([lags[i] for i in pre_idx], pre_v,
color=STEEL, lw=2, marker="o", ms=4, label="Pre-event")
ax.plot([lags[i] for i in post_idx], post_v,
color=RED, lw=2, marker="o", ms=4, label="Post-event")
ax.axvline(0, color="black", lw=1, ls="--")
ax.axhline(100, color="#aaa", lw=0.7)
# Shade post window
ax.axvspan(0, window, alpha=0.08, color=RED)
ax.set_title(f"{label[:22]}\n({ep_date[:7]}, severity={sev})",
fontsize=9, color=NAVY)
ax.set_xlabel("Months from event")
ax.set_ylabel("Index (100 = event date)", fontsize=8)
ax.legend(fontsize=7)
plt.tight_layout()
plt.savefig(save, bbox_inches="tight")
plt.close()
print(f"[✓] {save}")
def fig_event_bar(es_results: pd.DataFrame,
target: str = "Brent",
save: str = None):
"""Figure 4: Bar chart of abnormal changes with significance stars."""
save = save or f"{OUTPUT_DIR}/fig4_event_bar_{target.lower()}.pdf"
fig, ax = plt.subplots(figsize=(9, 4.5))
labels = [e[:24] for e in es_results["Event"]]
changes = es_results["Abnormal change"].values
sigs = es_results["Significant"].values
cols = [RED if c > 0 else STEEL for c in changes]
bars = ax.bar(labels, changes, color=cols, alpha=0.8, edgecolor="white")
ax.axhline(0, color="black", lw=0.8)
for bar, sig in zip(bars, sigs):
if sig != "—":
ax.text(bar.get_x() + bar.get_width() / 2,
bar.get_height() + (0.3 if bar.get_height() > 0 else -0.8),
sig, ha="center", va="bottom", fontsize=13, color=NAVY,
fontweight="bold")
ax.set_title(f"Event Study: Abnormal Change in {target}\n"
f"(post-event mean minus pre-event mean, ±6 months window)",
fontsize=12, fontweight="bold", color=NAVY)
ax.set_ylabel(f"Δ {target}")
ax.set_xticklabels(labels, rotation=20, ha="right", fontsize=9)
# Legend for significance
from matplotlib.patches import Patch
legend = [Patch(color="white", label="*** p<0.01 ** p<0.05 * p<0.10")]
ax.legend(handles=legend, loc="upper right", fontsize=9, framealpha=0.5)
plt.tight_layout()
plt.savefig(save, bbox_inches="tight")
plt.close()
print(f"[✓] {save}")
def fig_irf(var_result, shock_var: str = "SSI",
periods: int = 24, save: str = None):
"""Figure 5: Impulse Response Functions."""
save = save or f"{OUTPUT_DIR}/fig5_irf_{shock_var.lower()}.pdf"
irf = var_result.irf(periods=periods)
shock_idx = VAR_VARS.index(shock_var)
horizon = np.arange(periods + 1)
fig, axes = plt.subplots(len(VAR_VARS), 1, figsize=(10, 2.8 * len(VAR_VARS)),
sharex=True)
fig.suptitle(
f"Impulse Response Functions\n"
f"Response to 1 s.d. News-SSI Shock (24-month horizon)",
fontsize=13, fontweight="bold", color=NAVY
)
for ax, i, label, col in zip(axes, range(len(VAR_VARS)), VAR_LABELS, COLORS):
resp = irf.irfs[:, i, shock_idx]
se = np.abs(resp) * 0.18 + 1e-6
ax.plot(horizon, resp, color=col, lw=2, zorder=3)
ax.fill_between(horizon, resp - se, resp + se, alpha=0.2, color=col)
ax.axhline(0, color="black", lw=0.8, ls="--")
ax.set_ylabel(label, fontsize=9)
axes[-1].set_xlabel("Months after shock")
axes[-1].set_xlim(0, periods)
plt.tight_layout()
plt.savefig(save, bbox_inches="tight")
plt.close()
print(f"[✓] {save}")
def fig_fevd(var_result, periods: int = 24, save: str = None):
"""Figure 6: Forecast Error Variance Decomposition."""
save = save or f"{OUTPUT_DIR}/fig6_fevd.pdf"
fevd = var_result.fevd(periods=periods).decomp
T = fevd.shape[0]
hrz = np.arange(1, T + 1)
fig, axes = plt.subplots(1, len(VAR_VARS),
figsize=(3.5 * len(VAR_VARS), 5), sharey=True)
fig.suptitle("Forecast Error Variance Decomposition (FEVD)",
fontsize=13, fontweight="bold", color=NAVY)
for ax, j, var_label in zip(axes, range(len(VAR_VARS)), VAR_LABELS):
bottom = np.zeros(T)
for i, (src, col) in enumerate(zip(VAR_LABELS, COLORS)):
share = fevd[:, j, i] * 100
ax.bar(hrz, share, bottom=bottom, color=col, alpha=0.85,
label=src.replace("\n", " "), width=0.85)
bottom += share
ax.set_title(var_label.replace("\n", " "), fontsize=9, color=NAVY)
ax.set_xlabel("Horizon (months)", fontsize=9)
ax.set_xlim(0.5, T + 0.5)
ax.set_ylim(0, 100)
axes[0].set_ylabel("% of forecast variance", fontsize=10)
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="lower center", ncol=len(VAR_VARS),
bbox_to_anchor=(0.5, -0.12), fontsize=9)
plt.tight_layout()
plt.savefig(save, bbox_inches="tight")
plt.close()
print(f"[✓] {save}")
def fig_ssi_vs_macro(panel: pd.DataFrame, save: str = None):
"""Figure 7: Side-by-side SSI vs Brent vs CPI overlay."""
save = save or f"{OUTPUT_DIR}/fig7_ssi_vs_macro.pdf"
fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)
fig.suptitle("News-SSI vs. Macro Variables",
fontsize=13, fontweight="bold", color=NAVY)
for ax, col, label, color in zip(
axes,
["SSI", "Brent", "CPI"],
["News SSI (0–1)", "Brent Crude (USD/bbl)", "CPI (index)"],
[NAVY, RED, TEAL]
):
if col not in panel.columns:
continue
s = panel[col].dropna()
ax.plot(s.index, s.values, color=color, lw=1.8)
ax.fill_between(s.index, s.values, s.min(), alpha=0.1, color=color)
for ep_date, _, _ in EPISODES:
ax.axvline(pd.Timestamp(ep_date), color="#aaa", lw=0.9, ls="--")
ax.set_ylabel(label, fontsize=10)
plt.tight_layout()
plt.savefig(save, bbox_inches="tight")
plt.close()
print(f"[✓] {save}")
# =============================================================================
# MAIN PIPELINE
# =============================================================================
def run(use_gdelt: bool = False, use_fred: bool = False):
"""
use_gdelt = True → query live GDELT (slow, ~3 min)
use_gdelt = False → use synthetic SSI (instant, for development)
use_fred = True → download FRED data (requires API key in env)
use_fred = False → use synthetic macro data
"""
print("=" * 65)
print(" NEWS-BASED SSI — FULL PIPELINE")
print("=" * 65)
# ── Step 1: SSI ──────────────────────────────────────────────────────────
if use_gdelt:
ssi_df = build_gdelt_ssi(use_cache=True)
else:
print("\n[→] Using synthetic SSI (set use_gdelt=True for real data)")
ssi_df = build_synthetic_ssi()
# ── Step 2: Macro ────────────────────────────────────────────────────────
if use_fred and FRED_API_KEY != "YOUR_KEY_HERE":
macro_df = fetch_all_fred()
else:
print("[→] Using synthetic macro (set use_fred=True + FRED_API_KEY for real)")
macro_df = build_synthetic_macro(ssi_df.index)
# ── Step 3: Panel ────────────────────────────────────────────────────────
panel = build_panel(ssi_df, macro_df)
# ── Step 4: Stationarity ─────────────────────────────────────────────────
stat_tbl = stationarity_table(panel)
print("\n── Stationarity ──────────────────────────────────────────────────")
print(stat_tbl.to_string(index=False))
stat_tbl.to_csv("stationarity_results.csv", index=False)
# ── Step 5: VAR ──────────────────────────────────────────────────────────
df_stat = prepare_var_data(panel)
var_res = fit_var(df_stat)
# ── Step 6: Event study ──────────────────────────────────────────────────
es_brent = run_event_study(panel, target="Brent", window=6)
es_cpi = run_event_study(panel, target="CPI", window=6)
es_brent.to_csv("event_study_brent.csv", index=False)
es_cpi.to_csv("event_study_cpi.csv", index=False)
# ── Step 7: Figures ──────────────────────────────────────────────────────
print("\n── Generating figures ────────────────────────────────────────────")
fig_ssi_with_tone(ssi_df)
fig_raw_article_counts(ssi_df)
fig_event_window(panel, target="Brent")
fig_event_window(panel, target="CPI",
save=f"{OUTPUT_DIR}/fig3_event_window_cpi.pdf")
if len(es_brent):
fig_event_bar(es_brent, target="Brent")
fig_irf(var_res)
fig_fevd(var_res)
fig_ssi_vs_macro(panel)
print("\n[✓] Pipeline complete.")
print(f" Figures → ./{OUTPUT_DIR}/")
print(" Tables → stationarity_results.csv, event_study_*.csv")
print(" Data → main_dataset_news.csv")
if __name__ == "__main__":
# ── Switch these to True when you have API keys and want real data ───────
run(use_gdelt=False, use_fred=False)