CI Training module: Port Infrastructure Resilience under Coastal Hazards

This notebook walks through an end-to-end cyberinfrastructure (CI) workflow that connects:

• Coastal hazard signals (NOAA CO-OPS water-level anomalies, plus NOAA NCEI storm exposure features), with
• Operational disruption metrics derived from AIS broadcast points near ports.

We build a daily panel dataset (port × date), identify hazard events (extreme water-level days and storm-day sequences), compute resilience metrics (drop and recovery), and run regression models to quantify how hazards relate to vessel activity.

Ports in this training run: Houston and Norfolk.

By the end of this notebook, you will be able to:

1. Build a daily panel of port × date (data combined for all hazard and port) dataset by integrating AIS activity with coastal hazard variables
2. Define hazard events and event windows and compute resilience metrics (viz. amount of drop in port operations and time taken to recover operations)
3. Visualize and interpret operation disruption patterns
4. Quantify hazard–activity relationships using regression based on individual hazard events and on daily panel
5. Demonstrate parallel processing for estimating uncertainty for correlated time-series data using resampling

---

SECTION 0: Library import and data path setup

We start by importing core scientific Python packages and defining a consistent folder layout:

• data/ contains precomputed daily datasets (AIS activity, water-level hazards, storm features) to be used for the training module.
• outputs/ will store intermediate panels, event-level metrics, and any plots/figures.

import os, re, time
from pathlib import Path
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from event_viz import make_event_slider_map # helper py file: event_viz.py
from IPython.display import IFrame

import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats import spearmanr

from joblib import Parallel, delayed
import warnings

DATA_DIR = Path("data")

OUT_DIR  = Path("outputs")
OUT_DIR.mkdir(parents=True, exist_ok=True)

---

SECTION 1: Water-level hazards (NOAA CO-OPS)

1.1) Load daily water-level hazards (NOAA CO-OPS)

waterlevel_daily_2019_2024.parquet contains per-port daily summaries including:

• wl_mean, wl_max, wl_min (daily water level stats)
• wl_rollmed_30d (30-day rolling median baseline)
• wl_anom (anomaly relative to baseline)
• port / port_name and date

haz = pd.read_parquet(DATA_DIR / "waterlevel_daily_2019_2024.parquet").copy()
haz["date"] = pd.to_datetime(haz["date"]).dt.normalize()

haz.head()

	date	wl_mean	wl_max	wl_min	wl_rollmed_30d	wl_anom	port_name
0	2019-01-01 00:00:00+00:00	0.096850	0.275	-0.191	NaN	NaN	Port of Houston (Galveston gauge)
1	2019-01-02 00:00:00+00:00	0.157825	0.434	-0.175	NaN	NaN	Port of Houston (Galveston gauge)
2	2019-01-03 00:00:00+00:00	0.227037	0.399	0.038	NaN	NaN	Port of Houston (Galveston gauge)
3	2019-01-04 00:00:00+00:00	-0.156642	0.284	-0.485	NaN	NaN	Port of Houston (Galveston gauge)
4	2019-01-05 00:00:00+00:00	-0.159746	0.213	-0.474	NaN	NaN	Port of Houston (Galveston gauge)

haz.shape

(4384, 7)

---

Next we:

1. normalize dates to day-level timestamps, and
2. map gauge/port naming to the AIS port identifiers used in the activity dataset.

# normalize dates
haz["date"] = pd.to_datetime(haz["date"], utc=True).dt.tz_convert(None).dt.normalize()

# map hazard gauge names -> AIS port short names
PORT_TO_HAZARD = {
    "Houston": "Port of Houston (Galveston gauge)",
    "Norfolk": "Port of Norfolk (Sewells Point gauge)"
}

inv = {v: k for k, v in PORT_TO_HAZARD.items()}
haz["port"] = haz["port_name"].map(inv)

print(haz["port"].value_counts(dropna=False).head(10))

port
Houston    2192
Norfolk    2192
Name: count, dtype: int64

haz2 = haz[haz["port"].isin(["Houston", "Norfolk"])].copy()
print(haz2.shape)

(4384, 8)

---

1.2) Define “extreme” water-level days (per port)

Ports can have different baseline variability in water level anomalies (e.g., due to local hydrodynamics and gauge behavior). To compare ports fairly, we define extremes relative to each port’s own distribution.

Here we compute a 95th percentile threshold of wl_anom for each port and create:

• wl_thr: per-port anomaly threshold
• extreme_wl: 1 if the day's anomaly exceeds the threshold, else 0

This produces ~5% “extreme” days within each port’s own time series, which becomes the “event trigger” for the water-level event study later.

thr = (haz2.dropna(subset=["wl_anom"])
          .groupby("port")["wl_anom"]
          .quantile(0.95)
          .rename("wl_thr")
          .reset_index())

haz2 = haz2.merge(thr, on="port", how="left")

haz2["extreme_wl"] = np.where(
    haz2["wl_anom"].notna() & (haz2["wl_anom"] >= haz2["wl_thr"]),
    1, 0
)

print(haz2.groupby("port")["extreme_wl"].sum())

port
Houston    110
Norfolk    110
Name: extreme_wl, dtype: int64

---

1.3) Load daily AIS activity near ports

AIS preprocessing note : Raw AIS vessel transit broadcast points are extremely large (~850GB) and were processed in a separate notebook. That upstream step identifies the two most extreme water-level events per port, subsets AIS to a ±7 day window around each event, applies a port buffer, and aggregates points into daily port-level activity metrics. For training and runtime feasibility, we begin here from the resulting parquet output.

daily_ais_port_activity_event_sample_2019_2024.parquet contains daily operational activity metrics derived from AIS broadcast points within a port buffer, such as:

• unique_vessels: number of distinct MMSIs observed near the port that day
• n_points: total AIS points observed (proxy for traffic intensity)
• port and date identifiers

Because this file is event-sampled (top 2 events per port-year, ±7 days), the same (port, date) can appear more than once when event windows overlap; hence presence of duplicates.

We de-duplicate by (port, date) and then merge hazards onto this AIS table to form a single daily panel.

daily_ais = pd.read_parquet(DATA_DIR / "daily_ais_port_activity_event_sample_2019_2024.parquet")
daily_ais = daily_ais.drop_duplicates(subset=["port","date"])
print("After dedupe:", daily_ais.shape)

After dedupe: (1294, 4)

daily_ais["date"] = pd.to_datetime(daily_ais["date"]).dt.normalize()

panel = daily_ais.merge(
    haz2[["date","port","wl_anom","extreme_wl"]],
    on=["date","port"],
    how="left"
)

print("Panel shape:", panel.shape)
print("Missing hazard rows:", panel["wl_anom"].isna().sum())
panel.head()

Panel shape: (1294, 6)
Missing hazard rows: 12

	date	port	unique_vessels	n_points	wl_anom	extreme_wl
0	2019-01-04	Houston	422	240476	NaN	0
1	2019-01-05	Houston	425	226553	NaN	0
2	2019-01-06	Houston	419	235795	NaN	0
3	2019-01-07	Houston	402	240075	NaN	0
4	2019-01-08	Houston	411	249350	NaN	0

---

Merge sanity checks

After merging hazards onto AIS, we quickly check for missing wl_anom values.
Absence of the values usually comes from:

• port naming mismatches (e.g., gauge name vs AIS port label),
• gaps in the hazard time series,
• differences in date normalization / time zones.

We inspect missing rows by port and date range before continuing.

missing = panel[panel["wl_anom"].isna()].sort_values(["port","date"])
print("Missing rows:", len(missing))
display(missing.head(20))

print("Date range of missing by port:")
display(missing.groupby("port")["date"].agg(["min","max","count"]))

Missing rows: 12

	date	port	unique_vessels	n_points	wl_anom	extreme_wl
0	2019-01-04	Houston	422	240476	NaN	0
1	2019-01-05	Houston	425	226553	NaN	0
2	2019-01-06	Houston	419	235795	NaN	0
3	2019-01-07	Houston	402	240075	NaN	0
4	2019-01-08	Houston	411	249350	NaN	0
5	2019-01-09	Houston	403	267216	NaN	0
647	2019-01-04	Norfolk	184	84394	NaN	0
648	2019-01-05	Norfolk	167	84544	NaN	0
649	2019-01-06	Norfolk	177	84492	NaN	0
650	2019-01-07	Norfolk	183	85176	NaN	0
651	2019-01-08	Norfolk	177	86807	NaN	0
652	2019-01-09	Norfolk	183	99482	NaN	0

Date range of missing by port:

	min	max	count
port
Houston	2019-01-04	2019-01-09	6
Norfolk	2019-01-04	2019-01-09	6

panel = panel.copy()
panel["log_unique_vessels"] = np.log1p(panel["unique_vessels"])
panel["log_points"] = np.log1p(panel["n_points"])

# for analyses requiring wl_anom, drop rows where wl_anom is NaN
panel_nonnull = panel.dropna(subset=["wl_anom"]).copy()
print("panel_nonnull:", panel_nonnull.shape)

panel_nonnull: (1282, 8)

---

1.4) Quick pooled extreme threshold (for a simple impact check of high water levels on port activity)

Note: “extreme_wl” in this section is based on a pooled 95% threshold for a quick check; the earlier per-port threshold is the main event-definition used in the hazard labeling.

thr = panel_nonnull["wl_anom"].quantile(0.95)
panel_nonnull["extreme_wl"] = (panel_nonnull["wl_anom"] >= thr).astype(int)

panel_nonnull["extreme_wl"].value_counts()

extreme_wl
0    1217
1      65
Name: count, dtype: int64

# 65 occurences of abnormal water level (according to 95% threshold)
# Same output as above, but with more stats
panel_nonnull.groupby("extreme_wl")["wl_anom"].describe()

	count	mean	std	min	25%	50%	75%	max
extreme_wl
0	1217.0	-0.012627	0.140086	-0.673113	-0.088065	0.005658	0.082604	0.245325
1	65.0	0.313541	0.075775	0.245417	0.263654	0.288963	0.341333	0.621317

# distribution of the anomalities
panel_nonnull.groupby("port")["extreme_wl"].sum()

port
Houston    22
Norfolk    43
Name: extreme_wl, dtype: int64

Out of the 65 anomalous water-level events, 22 occured in Houston and 43 in Norfolk

---

1.5) Event-window construction for an event-study view (±7 days)

Now that we have a daily port × date panel with wl_anom, AIS activity, and an extreme_wl flag, we reshape the data into an event-study format.

For each port:

1. Identify all dates where extreme_wl == 1 (our “event dates”).
2. Extract a ±7 day window around each event date.
3. Create a relative-time index tau (days from the event; 0 = event day, -7..+7 = pre/post days).

This produces event_panel, which contains multiple rows per (port, tau) because multiple events contribute windows.
We then average across all available windows to get an average disruption/recovery curve by port and relative day.

def build_event_windows_from_panel(df, window=7):
    rows = []
    for port, g in df.groupby("port"):
        g = g.sort_values("date")
        event_dates = g.loc[g["extreme_wl"] == 1, "date"].unique()
        for ed in event_dates:
            w = g[(g["date"] >= ed - pd.Timedelta(days=window)) &
                  (g["date"] <= ed + pd.Timedelta(days=window))].copy()
            w["event_date"] = pd.Timestamp(ed)
            w["tau"] = (w["date"] - pd.Timestamp(ed)).dt.days
            rows.append(w)
    if not rows:
        return pd.DataFrame()
    return pd.concat(rows, ignore_index=True)

event_panel = build_event_windows_from_panel(panel_nonnull, window=7)
print("event_panel:", event_panel.shape)

event_panel: (904, 10)

curve = (event_panel
         .groupby(["port","tau"])
         .agg(
             mean_log_unique=("log_unique_vessels","mean"),
             mean_unique=("unique_vessels","mean"),
             mean_wl_anom=("wl_anom","mean"),
             n_obs=("unique_vessels","size")
         )
         .reset_index()
         .sort_values(["port","tau"])
)

display(curve.head())

curve_path = OUT_DIR / "event_study_curve_extreme_wl_2019_2024.parquet"
curve.to_parquet(curve_path, index=False)
print("Saved:", curve_path)

	port	tau	mean_log_unique	mean_unique	mean_wl_anom	n_obs
0	Houston	-7	5.999700	404.190476	0.012912	21
1	Houston	-6	6.004539	406.136364	-0.046236	22
2	Houston	-5	6.021919	413.428571	-0.078628	21
3	Houston	-4	6.016222	410.150000	-0.045547	20
4	Houston	-3	6.019761	412.190476	0.004189	21

Saved: outputs/event_study_curve_extreme_wl_2019_2024.parquet

def per_event_metrics(event_df, recovery_frac=0.9):
    out = []
    for (port, event_date), g in event_df.groupby(["port","event_date"]):
        g = g.sort_values("tau")

        base = g[g["tau"].between(-7, -1)]["unique_vessels"].mean()
        evt  = g[g["tau"] == 0]["unique_vessels"].mean()

        if pd.isna(base) or base <= 0 or pd.isna(evt):
            continue

        drop = (base - evt) / base

        thresh = recovery_frac * base
        post = g[g["tau"] >= 0]
        rec = post[post["unique_vessels"] >= thresh]
        recovery_days = int(rec["tau"].iloc[0]) if len(rec) else np.nan

        wl0 = g.loc[g["tau"] == 0, "wl_anom"]
        wl0 = float(wl0.iloc[0]) if len(wl0) else np.nan

        out.append({
            "port": port,
            "event_date": pd.Timestamp(event_date),
            "baseline_mean": float(base),
            "event_tau0": float(evt),
            "drop_frac": float(drop),
            "recovery_days_to_90pct": recovery_days,
            "wl_anom_tau0": wl0
        })
    return pd.DataFrame(out)

event_metrics = per_event_metrics(event_panel, recovery_frac=0.9) # Recovery fraction after disruption = 90%, meaning 90% operation recovered
print("event_metrics:", event_metrics.shape)
display(event_metrics.head())

event_metrics_path = OUT_DIR / "resilience_metrics_per_event_extreme_wl_2019_2024.parquet"
event_metrics.to_parquet(event_metrics_path, index=False)
print("Saved:", event_metrics_path)

event_metrics: (65, 7)

	port	event_date	baseline_mean	event_tau0	drop_frac	recovery_days_to_90pct	wl_anom_tau0
0	Houston	2019-02-01	426.857143	393.0	0.079317	0	0.276608
1	Houston	2019-02-19	389.428571	388.0	0.003668	0	0.313856
2	Houston	2019-03-13	371.857143	373.0	-0.003073	0	0.256092
3	Houston	2020-01-23	373.714286	434.0	-0.161315	0	0.277871
4	Houston	2020-04-05	411.833333	425.0	-0.031971	0	0.250640

Saved: outputs/resilience_metrics_per_event_extreme_wl_2019_2024.parquet

bench = (event_metrics
         .groupby("port")
         .agg(
             n_events=("event_date","count"),
             median_drop=("drop_frac","median"),
             mean_drop=("drop_frac","mean"),
             median_recovery=("recovery_days_to_90pct","median"),
             mean_recovery=("recovery_days_to_90pct","mean"),
             median_wl_anom=("wl_anom_tau0","median")
         )
         .reset_index()
         .sort_values("median_drop", ascending=False)
)

display(bench)

bench_path = OUT_DIR / "benchmark_resilience_summary_extreme_wl_2019_2024.csv"
bench.to_csv(bench_path, index=False)
print("Saved:", bench_path)

	port	n_events	median_drop	mean_drop	median_recovery	mean_recovery	median_wl_anom
0	Houston	22	0.022645	0.017370	0.0	0.090909	0.283398
1	Norfolk	43	-0.012931	-0.011648	0.0	0.139535	0.295852

Saved: outputs/benchmark_resilience_summary_extreme_wl_2019_2024.csv

Inference from the tables above:

median_drop & mean_drop help us answer the question: "During extreme WL events, how big was the discrete drop?"

Houston median_drop = 0.0226 : On a typical extreme water day, Houston saw about 2.3% drop in vessel activity

Norflok median_drop = -0.0129 : Negative here likely means effectively very small average impact

median_recovery = 0 : Activity resumed normally almost immediately

median_wl_anom: Water level on extreme days were roughly +28–30 cm above normal tide

Key takeaway: In this sample, extreme water-level anomaly events correspond to modest disruption in AIS activity. Houston shows a typical median drop of ~2.3% (median_drop ≈ 0.0226), while Norfolk shows near-zero / inconsistent change (median_drop ≈ -0.0129), with median recovery time ≈ 0 days (activity returns to ≥90% baseline immediately).

event_metrics[["wl_anom_tau0","drop_frac"]].corr()

	wl_anom_tau0	drop_frac
wl_anom_tau0	1.000000	0.237936
drop_frac	0.237936	1.000000

This tells us that Higher water levels correspond to larger drops in vessel activity. However, 0.24 corr value is a weak - to - moderate positivie relationship, since it's not zero either.

Thus, through this section, we know have confirmed that Higher water level days coincide with larger vessel activity drop.

---

SECTION 2: Storm events (IBTrACS)

AIS + storm preprocessing note (computed upstream):
To keep this training notebook lightweight, we start from a precompiled AIS storm-sample parquet. In a separate preprocessing notebook, we:

1. Download NOAA IBTrACS (North Atlantic) track data and compute daily storm exposure per port centroid, including a storm_day_200km flag (storm track points within 200 km) plus min distance and max wind statistics.
2. Treat storm_day_200km == 1 as a storm event and build ±7-day event windows, deduplicating overlapping windows and restricting to the AIS date range (2019-01-04 to 2024-04-20).
3. For the resulting storm-window days, download/extract the relevant daily AIS broadcast-point files and aggregate near-port activity within a 20 km buffer around each port centroid into daily metrics (unique_vessels, n_points).

This produces daily_ais_port_activity_storm_sample_2019_2024.parquet, which we load below.

---

2.1) Load storm exposure features (IBTrACS)

Storm exposure is summarized daily per port using the IBTrACS tropical cyclone tracks. Typical columns include:

• storm_day_200km: indicator that a storm was within 200 km of the port on that date
• storm_min_dist_km: closest approach distance (km)
• storm_max_wind_kt: maximum wind speed (knots) within the exposure window

We load:

1. an AIS activity table aligned to the storm sampling window,
2. the water-level hazards, and
3. the storm daily exposure table.

These are then merged into a unified daily panel.

# Load AIS storm-sample daily activity
ais_storm_path = DATA_DIR / "daily_ais_port_activity_storm_sample_2019_2024.parquet"
daily_ais = pd.read_parquet(ais_storm_path).copy()
daily_ais["date"] = pd.to_datetime(daily_ais["date"], errors="coerce").dt.normalize()

daily_ais = daily_ais.sort_values(["port","date"])
daily_ais = daily_ais.drop_duplicates(subset=["date","port"], keep="last").reset_index(drop=True) # drop duplicates

print("AIS storm-sample rows:", daily_ais.shape)
print("AIS date range:", daily_ais["date"].min(), "->", daily_ais["date"].max())

AIS storm-sample rows: (416, 4)
AIS date range: 2019-08-30 00:00:00 -> 2023-10-01 00:00:00

# Load waterlevel hazards
haz_path = DATA_DIR / "waterlevel_daily_2019_2024.parquet"
haz = pd.read_parquet(haz_path).copy()
haz["date"] = pd.to_datetime(haz["date"], errors="coerce", utc=True).dt.tz_convert(None).dt.normalize()

# map hazard port_name -> AIS port
PORT_TO_HAZARD = {
    "Houston": "Port of Houston (Galveston gauge)",
    "Norfolk": "Port of Norfolk (Sewells Point gauge)",
}
inv = {v: k for k, v in PORT_TO_HAZARD.items()}
haz["port"] = haz["port_name"].map(inv)

haz = haz.dropna(subset=["port"]).copy()

# Load storm daily features you computed from IBTrACS
storm_path = DATA_DIR / "storm_daily_ibtracs_NA_v04r01_2019_2024.parquet"
storm = pd.read_parquet(storm_path).copy()
storm["date"] = pd.to_datetime(storm["date"], errors="coerce").dt.normalize()

def to_binary_flag(s: pd.Series) -> pd.Series:
    x = s.astype(str).str.strip().str.lower()
    return x.isin(["1","1.0","true","t","yes","y"]).astype(int)

# check storm_day_200km is truly 0/1 int
storm["storm_day_200km"] = to_binary_flag(storm["storm_day_200km"])

---

2.2) Merge: AIS + waterlevel + storm

After merging, each row represents a (port, date) with:

• operational activity (AIS),
• hydrologic stress (water-level anomaly),
• storm exposure (distance/wind/indicator).

This panel is the backbone for both:

• event study metrics (drop + recovery around hazards), and
• regression modeling (average effect sizes and uncertainty).

panel = (
    daily_ais
    .merge(haz[["date","port","wl_anom"]], on=["date","port"], how="left")
    .merge(storm[["date","port","storm_day_200km","storm_min_dist_km","storm_max_wind_kt_200km","storm_pts_200km"]],
           on=["date","port"], how="left")
)

# fill storm flag missing as 0
panel["storm_day_200km"] = panel["storm_day_200km"].fillna(0).astype(int)

# log features
panel["log_unique_vessels"] = np.log1p(panel["unique_vessels"])
panel["log_points"] = np.log1p(panel["n_points"])

print("\nPanel shape:", panel.shape)
print("Missing wl_anom:", int(panel["wl_anom"].isna().sum()))
print("Storm days (by port):")
print(panel.groupby("port")["storm_day_200km"].sum())

Panel shape: (416, 11)
Missing wl_anom: 0
Storm days (by port):
port
Houston    10
Norfolk    13
Name: storm_day_200km, dtype: int64

# Sanity check: show a few storm-day rows that overlap AIS
storm_overlap = panel[panel["storm_day_200km"] == 1].sort_values(["port","date"]).head()
display(storm_overlap)

	date	port	unique_vessels	n_points	wl_anom	storm_day_200km	storm_min_dist_km	storm_max_wind_kt_200km	storm_pts_200km	log_unique_vessels	log_points
18	2019-09-17	Houston	418	263447	0.357877	1	28.934956	40.0	5	6.037871	12.481611
19	2019-09-18	Houston	411	251568	0.161696	1	33.003822	30.0	8	6.021023	12.435473
20	2019-09-19	Houston	405	229145	0.117319	1	158.594207	25.0	1	6.006353	12.342115
90	2020-08-27	Houston	360	242152	0.590460	1	189.346471	130.0	2	5.888878	12.397325
105	2020-09-21	Houston	292	163536	0.732352	1	193.146426	43.0	2	5.680173	12.004795

---

2.3) Visualising AIS activity vs water-level anomaly

def plot_daily_timeseries(panel, ports=None, date_min=None, date_max=None):
    df = panel.copy()
    if ports is not None:
        df = df[df["port"].isin(ports)]
    if date_min is not None:
        df = df[df["date"] >= pd.to_datetime(date_min)]
    if date_max is not None:
        df = df[df["date"] <= pd.to_datetime(date_max)]

    ports_list = sorted(df["port"].unique())
    n = len(ports_list)

    fig, axes = plt.subplots(n, 1, figsize=(14, 3.5*n), sharex=True)
    if n == 1:
        axes = [axes]

    for ax, port in zip(axes, ports_list):
        g = df[df["port"] == port].sort_values("date")

        ax.plot(g["date"], g["unique_vessels"], label="unique_vessels", color="orange")
        ax.set_title(f"{port}: Daily AIS activity vs Water-level anomaly (wl_anom)")
        ax.set_ylabel("unique_vessels")

        ax2 = ax.twinx()
        ax2.plot(g["date"], g["wl_anom"], linestyle="--", label="wl_anom")
        ax2.set_ylabel("wl_anom")

        # storm overlay: shade days where storm_day_200km == 1 (if present)
        if "storm_day_200km" in g.columns:
            storm_days = g[g["storm_day_200km"] == 1]["date"]
            for d in storm_days:
                ax.axvspan(d, d + pd.Timedelta(days=1), alpha=0.1, color="red")

        # combine legends
        lines, labels = ax.get_legend_handles_labels()
        lines2, labels2 = ax2.get_legend_handles_labels()
        ax.legend(lines + lines2, labels + labels2, loc="upper right")

    plt.tight_layout()
    plt.show()

plot_daily_timeseries(panel)

This plot overlays two daily signals for each port:

• Orange (left axis): unique_vessels = number of distinct AIS MMSIs observed within the port buffer (proxy for port activity).
• Blue dashed (right axis): wl_anom = water-level anomaly relative to the rolling baseline.
• Red shading: days where storm_day_200km == 1 (a storm track was within 200 km of the port centroid).

We use this figure as a quick visual check for

(i) data alignment across sources,
(ii) whether hazard spikes coincide with activity changes, and
(iii) whether storm days correspond to sharper disruptions. 

Note: AIS activity is event-window sampled, so the line can connect non-adjacent days; interpret trends within the dense clusters rather than across gaps.Hence, visual patterns here are suggestive; formal effect estimates are computed later using event windows and regression models.

From the graph, it's clear that as the water level increases above a certain level at the ports, the vessel activity decreseases, drastically falls down in case of the storm presence.

---

2.4) Water-level anomaly distribution and the “extreme” threshold (q=0.95)

We use this plot to:

• confirm the anomaly values are centered near 0 with a right tail (higher-than-normal water levels), and
• visualize how the “extreme” threshold differs by port based on each port’s anomaly distribution.

def plot_wl_anom_distribution(panel, q=0.95):
    df = panel.copy()
    ports = sorted(df["port"].unique())

    fig, axes = plt.subplots(len(ports), 1, figsize=(12, 3.5*len(ports)))
    if len(ports) == 1:
        axes = [axes]

    for ax, port in zip(axes, ports):
        g = df[df["port"] == port]
        x = g["wl_anom"].dropna().values

        thr = np.quantile(x, q) if len(x) else np.nan
        ax.hist(x, bins=40)
        ax.axvline(thr, linestyle="--", linewidth=2)
        ax.set_title(f"{port}: wl_anom distribution (q={q:.2f} threshold shown)")
        ax.set_xlabel("wl_anom")
        ax.set_ylabel("count")

    plt.tight_layout()
    plt.show()

plot_wl_anom_distribution(panel, q=0.95)

In this sample, wl_anom values for both ports are centered near ~0, which is consistent with the anomaly being computed relative to a rolling baseline (most days are “normal” after detrending). Both distributions are slightly right-skewed, meaning unusually high water-level anomalies occur and form a longer positive tail than unusually low-water days.

The dashed line marks the 95th-percentile (q=0.95) cutoff for each port (≈ the top 5% highest anomalies). In this plot the cutoff is around ~0.3 for both ports, so “extreme” days correspond to roughly wl_anom ≥ 0.3 in this dataset.

Houston shows a wider/heavier positive tail (more large outliers) than Norfolk, suggesting greater variability/magnitude of high-water anomalies in Houston for the days represented here.

Note: this histogram reflects the rows present in the merged panel (event-sampled AIS windows), not necessarily the full continuous 2019–2024 hazard record.

out_path = OUT_DIR / "panel_daily_ais_wl_storm_STORM_SAMPLE_2019_2024.parquet"
panel.to_parquet(out_path, index=False)
print("\nSaved:", out_path)

Saved: outputs/panel_daily_ais_wl_storm_STORM_SAMPLE_2019_2024.parquet

---

2.5) Storm event-study: from storm-day flags to disruption/recovery curves

In the daily panel, storm_day_200km is a binary daily indicator (1 = at least one IBTrACS storm track point within 200 km of the port centroid on that date). Because a single storm can affect a port across multiple consecutive days, we first convert this daily flag into storm events per port by collapsing consecutive storm days into one event and defining event_date as the first day of each run.

We then build an event-study panel around each storm event by generating a relative-time index tau for a fixed window (tau = -7…+14, where tau=0 is event_date). For each (port, event_date, tau), we merge in outcomes and covariates from the main panel, including AIS activity (unique_vessels, n_points), water-level anomaly (wl_anom), and storm intensity metrics (storm_min_dist_km, storm_max_wind_kt_200km, storm_pts_200km). Some tau-days may have missing AIS outcomes if they fall outside the precomputed AIS sample; these are expected and tracked.

Finally, we aggregate across all storm events to obtain an average disruption/recovery curve by port: grouping by (port, tau) and computing mean/median activity, log-scaled activity, n_obs (how many observations support each tau), and average hazard context. The resulting event-study panel and curve are saved as compact parquet outputs for downstream plotting and modeling.

panel["date"] = pd.to_datetime(panel["date"]).dt.normalize()

''' 
    Define storm events (per port)
        We collapse consecutive storm days into ONE event:
        Event_date = first day of each consecutive run
'''

def make_events_from_binary(df, date_col="date", flag_col="storm_day_200km", group_col="port"):
    df = df.sort_values([group_col, date_col]).copy()
    df["prev_flag"] = df.groupby(group_col)[flag_col].shift(1).fillna(0).astype(int)
    df["prev_date"] = df.groupby(group_col)[date_col].shift(1)

    # new event starts if flag==1 and either previous flag==0 OR date gap > 1 day
    is_new = (df[flag_col] == 1) & (
        (df["prev_flag"] == 0) |
        ((df[date_col] - df["prev_date"]).dt.days > 1)
    )

    events = df.loc[is_new, [group_col, date_col]].rename(columns={date_col: "event_date"})
    events = events.reset_index(drop=True)
    return events

events = make_events_from_binary(panel, flag_col="storm_day_200km")
print("Storm events:", events.shape)
display(events.head(20))

Storm events: (14, 2)

	port	event_date
0	Houston	2019-09-17
1	Houston	2020-08-27
2	Houston	2020-09-21
3	Houston	2020-10-09
4	Houston	2021-09-14
5	Norfolk	2019-09-06
6	Norfolk	2019-10-20
7	Norfolk	2020-07-09
8	Norfolk	2020-08-04
9	Norfolk	2020-08-14
10	Norfolk	2020-10-29
11	Norfolk	2021-06-21
12	Norfolk	2021-07-08
13	Norfolk	2023-09-23

# Build tau-window panel around each event (e.g., -7..+14)
PRE_DAYS  = 7
POST_DAYS = 14

taus = np.arange(-PRE_DAYS, POST_DAYS + 1)

event_panel_rows = []
for _, r in events.iterrows():
    port = r["port"]
    event_date = r["event_date"]
    for tau in taus:
        d = event_date + pd.Timedelta(days=int(tau))
        event_panel_rows.append({"port": port, "event_date": event_date, "tau": int(tau), "date": d})

event_panel = pd.DataFrame(event_panel_rows)

# Merge in outcomes + hazards from main panel
keep_cols = [
    "date","port","unique_vessels","n_points",
    "wl_anom",
    "storm_day_200km","storm_min_dist_km","storm_max_wind_kt_200km","storm_pts_200km"
]
event_panel = event_panel.merge(panel[keep_cols], on=["date","port"], how="left")

print("event_panel shape:", event_panel.shape)
print("Missing unique_vessels:", int(event_panel["unique_vessels"].isna().sum()))
display(event_panel.head())

event_panel_path = OUT_DIR / "event_study_panel_storm_2019_2024.parquet"
event_panel.to_parquet(event_panel_path, index=False)
print("Saved:", event_panel_path)

event_panel shape: (308, 11)
Missing unique_vessels: 57

	port	event_date	tau	date	unique_vessels	n_points	wl_anom	storm_day_200km	storm_min_dist_km	storm_max_wind_kt_200km	storm_pts_200km
0	Houston	2019-09-17	-7	2019-09-10	426.0	226997.0	0.124567	0.0	4725.128012	NaN	0.0
1	Houston	2019-09-17	-6	2019-09-11	422.0	228951.0	0.140450	0.0	5568.762390	NaN	0.0
2	Houston	2019-09-17	-5	2019-09-12	404.0	232822.0	0.198871	0.0	2187.560086	NaN	0.0
3	Houston	2019-09-17	-4	2019-09-13	444.0	227440.0	0.182317	0.0	2060.820653	NaN	0.0
4	Houston	2019-09-17	-3	2019-09-14	446.0	215787.0	0.183025	0.0	1785.763568	NaN	0.0

Saved: outputs/event_study_panel_storm_2019_2024.parquet

# Obtain event-study curve (mean/median by tau, per port)
curve = (
    event_panel
    .groupby(["port","tau"], as_index=False)
    .agg(
        mean_unique=("unique_vessels","mean"),
        median_unique=("unique_vessels","median"),
        mean_log_unique=("unique_vessels", lambda x: np.log1p(x).mean() if x.notna().any() else np.nan),
        n_obs=("unique_vessels","count"),
        mean_wl_anom=("wl_anom","mean"),
        storm_flag_rate=("storm_day_200km","mean"),
        mean_min_dist=("storm_min_dist_km","mean"),
        mean_max_wind=("storm_max_wind_kt_200km","mean"),
    )
)

curve_path = OUT_DIR / "event_study_curve_storm_2019_2024.parquet"
curve.to_parquet(curve_path, index=False)
print("Saved:", curve_path)
display(curve.head())

Saved: outputs/event_study_curve_storm_2019_2024.parquet

	port	tau	mean_unique	median_unique	mean_log_unique	n_obs	mean_wl_anom	storm_flag_rate	mean_min_dist	mean_max_wind
0	Houston	-7	418.6	426.0	6.037779	5	0.034572	0.0	2636.572784	NaN
1	Houston	-6	428.2	422.0	6.058990	5	0.072347	0.0	2014.176808	NaN
2	Houston	-5	427.2	420.0	6.056828	5	0.116661	0.0	1244.089963	NaN
3	Houston	-4	421.0	417.0	6.044026	5	0.110689	0.0	1329.546649	NaN
4	Houston	-3	421.2	416.0	6.044666	5	0.122203	0.0	1583.706888	NaN

---

2.6) Visualizing average disruption and recovery in AIS activity

This figure plots the event-study curve of port activity around storm events.
For each port, we stack ±7 to +14 day windows around each event_date (the first day of a consecutive storm run) and compute mean_unique = the average number of unique vessels observed near the port for each relative day tau.

def plot_event_study_curve(curve, y_col="mean_unique"):
    df = curve.copy()
    ports = sorted(df["port"].unique())

    fig, axes = plt.subplots(1, len(ports), figsize=(6*len(ports), 4), sharey=True)
    if len(ports) == 1:
        axes = [axes]

    for ax, port in zip(axes, ports):
        g = df[df["port"] == port].sort_values("tau")
        ax.plot(g["tau"], g[y_col], marker="o")
        ax.axvline(0, linestyle="--")
        ax.set_title(f"{port}: Event-study curve ({y_col})")
        ax.set_xlabel("tau (days from event)")
        ax.grid(True, alpha=0.3)

    axes[0].set_ylabel(y_col)
    plt.tight_layout()
    plt.show()

plot_event_study_curve(curve, y_col="mean_unique")

• x-axis (tau): days relative to the storm event start (tau=0 is the event start day)
• y-axis (mean_unique): average daily unique vessels (AIS activity proxy)
• vertical dashed line: tau=0 (storm event start)

Inference from this plot:

• Houston: there is a clear drop in activity at tau=0 (storm onset), followed by a recovery over the next few days toward the pre-event baseline. This pattern is consistent with short-term operational disruption around storm start and gradual normalization afterward.
• Norfolk: there is also a dip at tau=0, but the magnitude is smaller and the curve is flatter overall, suggesting weaker or less consistent disruption in this sample (and/or fewer storm events contributing).
• Post-event fluctuations (e.g., later bumps/dips) likely reflect limited sample size and overlapping operational factors; we treat these curves as descriptive and quantify effects later with summary metrics/regression.

---

2.7) Per-storm-event resilience metrics (drop magnitude + recovery time)

This table summarizes event-level disruption and recovery for each storm event (per port), using the event-study panel.

For each (port, event_date):

• baseline_mean: mean unique_vessels over the pre-event baseline window tau = -7…-1
• event_tau0: unique_vessels on the event start day tau = 0
• drop_frac: fractional drop at event onset = (baseline_mean - event_tau0) / baseline_mean
• recovery_days_to_90pct: first tau ≥ 0 day when activity reaches ≥ 90% of baseline_mean
• wl_anom_tau0: water-level anomaly on the event day (tau=0) for context
• storm_min_dist_km_tau0, storm_max_wind_kt_tau0: storm proximity and intensity on tau=0 (IBTrACS-derived context)

These per-event metrics turn the full time-series window into a compact “impact + recovery” summary that we can compare across ports and relate to storm intensity.

# Per-event resilience metrics (drop + recovery time)
def per_event_metrics(event_df, recovery_frac=0.9):
    out = []
    for (port, event_date), g in event_df.groupby(["port","event_date"]):
        g = g.sort_values("tau")

        base = g.loc[g["tau"].between(-7, -1), "unique_vessels"].mean()
        evt  = g.loc[g["tau"] == 0, "unique_vessels"].mean()

        if pd.isna(base) or base <= 0 or pd.isna(evt):
            continue

        drop = (base - evt) / base

        thresh = recovery_frac * base
        post = g.loc[g["tau"] >= 0].copy()
        rec = post.loc[post["unique_vessels"] >= thresh]
        recovery_days = int(rec["tau"].iloc[0]) if len(rec) else np.nan

        wl0 = g.loc[g["tau"] == 0, "wl_anom"]
        wl0 = float(wl0.iloc[0]) if len(wl0) else np.nan

        mind0 = g.loc[g["tau"] == 0, "storm_min_dist_km"]
        mind0 = float(mind0.iloc[0]) if len(mind0) else np.nan

        wind0 = g.loc[g["tau"] == 0, "storm_max_wind_kt_200km"]
        wind0 = float(wind0.iloc[0]) if len(wind0) else np.nan

        out.append({
            "port": port,
            "event_date": pd.Timestamp(event_date),
            "baseline_mean": float(base),
            "event_tau0": float(evt),
            "drop_frac": float(drop),
            "recovery_days_to_90pct": recovery_days,
            "wl_anom_tau0": wl0,
            "storm_min_dist_km_tau0": mind0,
            "storm_max_wind_kt_tau0": wind0,
        })

    return pd.DataFrame(out)

storm_event_metrics = per_event_metrics(event_panel, recovery_frac=0.9)
print("storm_event_metrics:", storm_event_metrics.shape)
display(storm_event_metrics.sort_values(["port","event_date"]).head(20))

metrics_path = OUT_DIR / "resilience_metrics_per_event_storm_2019_2024.parquet"
storm_event_metrics.to_parquet(metrics_path, index=False)
print("Saved:", metrics_path)

storm_event_metrics: (14, 9)

	port	event_date	baseline_mean	event_tau0	drop_frac	recovery_days_to_90pct	wl_anom_tau0	storm_min_dist_km_tau0	storm_max_wind_kt_tau0
0	Houston	2019-09-17	427.857143	418.0	0.023038	0	0.357877	28.934956	40.0
1	Houston	2020-08-27	417.714286	360.0	0.138167	1	0.590460	189.346471	130.0
2	Houston	2020-09-21	413.285714	292.0	0.293467	3	0.732352	193.146426	43.0
3	Houston	2020-10-09	438.857143	365.0	0.168294	1	0.304992	186.222669	90.0
4	Houston	2021-09-14	380.714286	298.0	0.217261	1	0.488708	26.224295	65.0
5	Norfolk	2019-09-06	235.428571	167.0	0.290655	1	0.583879	178.540097	85.0
6	Norfolk	2019-10-20	281.714286	239.0	0.151623	1	0.003883	40.080863	40.0
7	Norfolk	2020-07-09	254.857143	231.0	0.093610	0	0.143919	191.482176	40.0
8	Norfolk	2020-08-04	223.714286	191.0	0.146232	1	-0.109342	104.735133	60.0
9	Norfolk	2020-08-14	221.000000	229.0	-0.036199	0	0.066823	188.328815	35.0
10	Norfolk	2020-10-29	305.428571	299.0	0.021048	0	-0.125456	199.196618	45.0
11	Norfolk	2021-06-21	287.428571	272.0	0.053678	0	-0.012352	96.365054	40.0
12	Norfolk	2021-07-08	267.142857	249.0	0.067914	0	-0.148431	190.049225	45.0
13	Norfolk	2023-09-23	292.428571	236.0	0.192965	1	0.653879	120.015615	45.0

Saved: outputs/resilience_metrics_per_event_storm_2019_2024.parquet

Inference from the above event-metrics table

In this sample there are 14 storm events total (Houston: 5; Norfolk: 9).
Most events show positive drops at tau=0, meaning vessel activity decreases on storm onset relative to the pre-event baseline. The drop magnitude varies substantially across events (from near-zero to ~0.29–0.30), indicating that storms do not all disrupt port activity equally.

Recovery is often fast since many events return to ≥90% of baseline within 0–1 days, while a few take longer (e.g., up to ~3 days). We also see at least one case with negative drop_frac (activity higher at tau=0 than baseline), which can happen due to operational scheduling noise, limited AIS sampling around events, or storms that were present within 200 km but did not materially disrupt near-port traffic.

Focussing on these two rows as example:

Houston 2020-09-21:

- baseline: 413 vessels
- storm day: 292 vessels
~29% drop in activity & Recovery took 3 days => Strong disruption

Norfolk 2019-09-06:

- baseline: 235 vessels
- storm day: 167 vessels
Also ~29% drop, but recovered in 1 day

def plot_recovery_days(storm_event_metrics):
    df = storm_event_metrics.copy()
    ports = sorted(df["port"].unique())

    fig, axes = plt.subplots(1, len(ports), figsize=(6*len(ports), 4), sharey=True)
    if len(ports) == 1:
        axes = [axes]

    for ax, port in zip(axes, ports):
        g = df[df["port"] == port]["recovery_days_to_90pct"].dropna()
        ax.hist(g, bins=np.arange(-0.5, 5.5, 1))
        ax.set_title(f"{port}: recovery_days_to_90pct")
        ax.set_xlabel("days")
        ax.grid(True, alpha=0.3)

    axes[0].set_ylabel("count")
    plt.tight_layout()
    plt.show()

plot_recovery_days(storm_event_metrics)

Norfolk experienced a higher number of storm events, but recovered faster compared to Houston. This could also be due to Houston having larger port activity than Norfolk (as we'll see from the data in next part).

---

2.8) Using Correlation heatmap to understand how storm intensity relates to disruption

This heatmap shows Pearson correlations across per-event storm metrics:

• drop_frac: fractional activity drop at storm onset (tau=0) relative to the pre-event baseline
• storm_max_wind_kt_tau0: storm wind intensity (knots) at tau=0 (within 200 km)
• storm_min_dist_km_tau0: closest storm distance (km) at tau=0
• wl_anom_tau0: water-level anomaly at tau=0

def corr_heatmap(df, cols, title):
    corr = df[cols].corr()
    fig, ax = plt.subplots(figsize=(6,5))
    im = ax.imshow(corr.values)
    ax.set_xticks(range(len(cols)))
    ax.set_yticks(range(len(cols)))
    ax.set_xticklabels(cols, rotation=45, ha="right")
    ax.set_yticklabels(cols)

    for i in range(len(cols)):
        for j in range(len(cols)):
            ax.text(j, i, f"{corr.values[i,j]:.2f}", ha="center", va="center")

    ax.set_title(title)
    plt.colorbar(im, ax=ax)
    plt.tight_layout()
    plt.show()

corr_heatmap(
    storm_event_metrics,
    cols=["drop_frac","storm_max_wind_kt_tau0","storm_min_dist_km_tau0","wl_anom_tau0"],
    title="Storm-event metrics correlation"
)

• drop_frac is positively correlated with wl_anom_tau0 (~0.69) and moderately with storm_max_wind_kt_tau0 (~0.38), suggesting larger drops tend to coincide with higher water anomalies and stronger winds.
• drop_frac shows near-zero correlation with storm_min_dist_km_tau0 (~ -0.01) here, which may reflect small sample size, the 200 km thresholding, and/or limited variation in distances among these events.

---

def corr_heatmap(df, cols, title):
    corr = df[cols].corr()
    fig, ax = plt.subplots(figsize=(6,5))
    im = ax.imshow(corr.values)
    ax.set_xticks(range(len(cols)))
    ax.set_yticks(range(len(cols)))
    ax.set_xticklabels(cols, rotation=45, ha="right")
    ax.set_yticklabels(cols)

    for i in range(len(cols)):
        for j in range(len(cols)):
            ax.text(j, i, f"{corr.values[i,j]:.2f}", ha="center", va="center")

    ax.set_title(title)
    plt.colorbar(im, ax=ax)
    plt.tight_layout()
    plt.show()

corr_heatmap(
    storm_event_metrics,
    cols=["drop_frac","storm_max_wind_kt_tau0","storm_min_dist_km_tau0","wl_anom_tau0"],
    title="Storm-event metrics correlation"
)

• drop_frac is positively correlated with wl_anom_tau0 (~0.69) and moderately with storm_max_wind_kt_tau0 (~0.38), suggesting larger drops tend to coincide with higher water anomalies and stronger winds.
• drop_frac shows near-zero correlation with storm_min_dist_km_tau0 (~ -0.01) here, which may reflect small sample size, the 200 km thresholding, and/or limited variation in distances among these events.

2.9) Summary statistics by port for storm events

We summarize the per-event resilience metrics by port to compare:

• how often storms occur in the sample (events per port)
• typical disruption magnitude (median drop_frac)
• typical recovery speed (median recovery_days_to_90pct)

# summaries
print("\nEvents per port:")
print(storm_event_metrics.groupby("port").size())

print("\nMedian drop_frac per port:")
print(storm_event_metrics.groupby("port")["drop_frac"].median())

print("\nMedian recovery days per port:")
print(storm_event_metrics.groupby("port")["recovery_days_to_90pct"].median())

Events per port:
port
Houston    5
Norfolk    9
dtype: int64

Median drop_frac per port:
port
Houston    0.168294
Norfolk    0.093610
Name: drop_frac, dtype: float64

Median recovery days per port:
port
Houston    1.0
Norfolk    0.0
Name: recovery_days_to_90pct, dtype: float64

Interpretation:

• Norfolk has more storm events (9) than Houston (5).
• Median disruption is larger in Houston (~0.168) than Norfolk (~0.094), i.e., Houston shows a bigger typical activity drop at storm onset.
• Median recovery time to 90% baseline is ~1 day for Houston and ~0 days for Norfolk (many Norfolk events return to ≥90% immediately at tau=0).

These summaries suggest that, in this dataset, storms are more frequent for Norfolk but (on average) produce stronger immediate shocks and slightly slower recovery for Houston.

While water-level anomaly events caused modest disruption (~2%), storm events cause significant shocks (10-30% drop).

---

SECTION 3: Hazard Event Comparison and Interactive Validation

3.1) Compare extreme water-level vs storm events, and “compound” overlap

At this point we have two per-event resilience tables:

• Extreme water-level events (from resilience_metrics_per_event_extreme_wl_2019_2024.parquet)
• Storm events (from storm_event_metrics)

Here we align them by (port, event_date) using an outer merge and label each row as:

• wl_only: has a water-level event metric but no storm metric on that event_date
• storm_only: has a storm metric but no water-level metric on that event_date
• compound: both metrics exist on the same event_date (overlap)

We then compare typical disruption magnitudes across event types by taking medians of:

• drop_frac_wl for wl_only
• drop_frac_storm for storm_only

# Load extreme WL event metrics
extreme_path = OUT_DIR / "resilience_metrics_per_event_extreme_wl_2019_2024.parquet"
extreme = pd.read_parquet(extreme_path).copy()

extreme["event_date"] = pd.to_datetime(extreme["event_date"]).dt.normalize()
storm_event_metrics["event_date"] = pd.to_datetime(storm_event_metrics["event_date"]).dt.normalize()

# Merge on port + event_date
combined = extreme.merge(
    storm_event_metrics,
    on=["port","event_date"],
    how="outer",
    suffixes=("_wl","_storm")
)

# Label event type
def classify(row):
    has_wl = not pd.isna(row.get("drop_frac_wl"))
    has_storm = not pd.isna(row.get("drop_frac_storm"))

    if has_wl and has_storm:
        return "compound"
    elif has_wl:
        return "wl_only"
    elif has_storm:
        return "storm_only"
    else:
        return "unknown"

combined["event_type"] = combined.apply(classify, axis=1)

print(combined["event_type"].value_counts())
display(combined.sort_values(["port","event_date"]).head())

event_type
wl_only       65
storm_only    14
Name: count, dtype: int64

	port	event_date	baseline_mean_wl	event_tau0_wl	drop_frac_wl	recovery_days_to_90pct_wl	wl_anom_tau0_wl	baseline_mean_storm	event_tau0_storm	drop_frac_storm	recovery_days_to_90pct_storm	wl_anom_tau0_storm	storm_min_dist_km_tau0	storm_max_wind_kt_tau0	event_type
0	Houston	2019-02-01	426.857143	393.0	0.079317	0.0	0.276608	NaN	NaN	NaN	NaN	NaN	NaN	NaN	wl_only
1	Houston	2019-02-19	389.428571	388.0	0.003668	0.0	0.313856	NaN	NaN	NaN	NaN	NaN	NaN	NaN	wl_only
2	Houston	2019-03-13	371.857143	373.0	-0.003073	0.0	0.256092	NaN	NaN	NaN	NaN	NaN	NaN	NaN	wl_only
3	Houston	2019-09-17	NaN	NaN	NaN	NaN	NaN	427.857143	418.0	0.023038	0.0	0.357877	28.934956	40.0	storm_only
4	Houston	2020-01-23	373.714286	434.0	-0.161315	0.0	0.277871	NaN	NaN	NaN	NaN	NaN	NaN	NaN	wl_only

combined.groupby("event_type")[
    ["drop_frac_wl","drop_frac_storm"]
].median()

	drop_frac_wl	drop_frac_storm
event_type
storm_only	NaN	0.1422
wl_only	-0.004274	NaN

The median wl_only drop is close to zero (small/inconsistent impact), while storm_only events show a clearly larger median drop (≈0.14 here). This supports the narrative that storm exposure produces the dominant operational shock, whereas water-level anomaly extremes alone are comparatively modest in this sample.

Storm exposure is the primary driver of operational disruption, not just anomalous water levels alone.

---

Next, we compute correlations within storm events (drop_frac vs storm wind, distance, and wl_anom_tau0) as an exploratory check to see which hazard signals tend to co-occur with larger drops.

storm_event_metrics[[
    "drop_frac",
    "storm_max_wind_kt_tau0",
    "storm_min_dist_km_tau0",
    "wl_anom_tau0"
]].corr()

	drop_frac	storm_max_wind_kt_tau0	storm_min_dist_km_tau0	wl_anom_tau0
drop_frac	1.000000	0.375428	-0.012868	0.688380
storm_max_wind_kt_tau0	0.375428	1.000000	0.244178	0.423879
storm_min_dist_km_tau0	-0.012868	0.244178	1.000000	0.010604
wl_anom_tau0	0.688380	0.423879	0.010604	1.000000

Inference from above table:

storm_max_wind_kt_tau0: Moderate positive correlation => Stronger winds -> larger vessel activity drop. Makes intuitive sense.

storm_min_dist_km_tau0: Almost zero => Distance alone is not explaining disruption.

wl_anom_tau0: Strong positive correlation. Larger surge anomaly -> larger activity drop.

This suggests that Storm wind alone matters, but surge/water-level anomaly is the strongest predictor

---

summary = combined.groupby("event_type").agg(
    median_drop_wl=("drop_frac_wl","median"),
    median_drop_storm=("drop_frac_storm","median"),
    median_rec_wl=("recovery_days_to_90pct_wl","median"),
    median_rec_storm=("recovery_days_to_90pct_storm","median"),
)
print(summary)

            median_drop_wl  median_drop_storm  median_rec_wl  median_rec_storm
event_type
storm_only             NaN             0.1422            NaN               1.0
wl_only          -0.004274                NaN            0.0               NaN

---

Finally, we run a Spearman rank correlation test for drop_frac vs wl_anom_tau0 to assess whether the monotonic relationship remains evident even with a small number of storm events.

from scipy.stats import spearmanr

spearmanr(storm_event_metrics["drop_frac"],
          storm_event_metrics["wl_anom_tau0"])

SignificanceResult(statistic=np.float64(0.6483516483516484), pvalue=np.float64(0.012144222573105738))

rho = 0.648 p = 0.012

Strong monotonic relationship

Statistically significant (p < 0.05)

Even with only 14 storm events (which is a small sample size)

Therefore:

Surge (water-level anomaly) explains disruption better than distance

Wind explains some variation

Distance alone explains almost nothing

WL-only anomalies do not disrupt ports much

---

3.2) Interactive event map to inspect individual events spatially

This cell generates an HTML event slider to explore events one-by-one. The map allows you to:

• choose a port and event type (WL extreme vs storm),
• step through events using a slider (each step is an event_date),
• view an AIS activity window around the event and basic impact metrics (e.g., drop_frac).

Under the hood, the slider uses:

• the merged daily panel (AIS + hazards),
• per-event metrics for storms (storm_event_metrics) and water-level extremes (wl_event_metrics),
• the daily IBTrACS-derived storm exposure table (storm_daily_*),
• fixed port coordinates (PORT_COORDS).

This visualization is meant as a qualitative companion to the tables/curves: it helps validate that event dates and windows look reasonable and lets us sanity-check a few high-impact or surprising events by looking at their spatial/temporal context.

from event_viz import make_event_slider_map
from IPython.display import IFrame

wl_event_metrics = pd.read_parquet(OUT_DIR/"resilience_metrics_per_event_extreme_wl_2019_2024.parquet")
wl_event_metrics["event_date"] = pd.to_datetime(wl_event_metrics["event_date"]).dt.normalize()

storm_daily = pd.read_parquet(DATA_DIR/"storm_daily_ibtracs_NA_v04r01_2019_2024.parquet")

PORT_COORDS = {
    "Houston":  (29.73, -95.27),
    "Norfolk":  (36.85, -76.29),
}

html_path = make_event_slider_map(
    panel=panel,
    storm_event_metrics=storm_event_metrics,
    wl_event_metrics=wl_event_metrics,
    storm_daily=storm_daily,
    out_dir=OUT_DIR,
    port_coords=PORT_COORDS,
    window_days=7,
    dots_cap=160,
    dot_scale_divisor=5.0,
)

IFrame(str(html_path), width=1100, height=650)

The separately generated interactive event-slider map is not included in this hosted notebook.

---

---

SECTION 4: Regression Modeling

We fit regressions for two different “response modes”:

1. Shock response (event-level):
   Use storm event metrics where the outcome is the immediate fractional drop in activity on the event date.

2. Gradual stress response (daily panel):
   Use the daily panel to estimate how activity changes with day-to-day water-level anomalies, controlling for port effects.

4.1) Shock response (Sudden collapse) regression

Here we quantify the sudden-collapse (“shock”) response to storms by modeling the event-level disruption metric:

• Outcome (y): drop_frac = fractional decline in AIS activity at storm onset (tau=0) relative to the pre-event baseline (tau=-7…-1).
• Predictors (X) (measured at tau=0):
  • storm_max_wind_kt_tau0 (storm intensity proxy)
  • storm_min_dist_km_tau0 (storm proximity)
  • wl_anom_tau0 (water-level anomaly / surge context)

We fit an OLS model with an intercept to ask: which hazard signals are most associated with larger immediate drops?

import statsmodels.api as sm

X = storm_event_metrics[
    ["storm_max_wind_kt_tau0",
     "storm_min_dist_km_tau0",
     "wl_anom_tau0"]
]

X = sm.add_constant(X)
y = storm_event_metrics["drop_frac"]

model = sm.OLS(y, X).fit()
print(model.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:              drop_frac   R-squared:                       0.485
Model:                            OLS   Adj. R-squared:                  0.330
Method:                 Least Squares   F-statistic:                     3.133
Date:                Tue, 24 Feb 2026   Prob (F-statistic):             0.0743
Time:                        00:25:41   Log-Likelihood:                 17.620
No. Observations:                  14   AIC:                            -27.24
Df Residuals:                      10   BIC:                            -24.68
Df Model:                           3
Covariance Type:            nonrobust
==========================================================================================
                             coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------------
const                      0.0648      0.064      1.020      0.332      -0.077       0.206
storm_max_wind_kt_tau0     0.0004      0.001      0.446      0.665      -0.002       0.003
storm_min_dist_km_tau0 -7.096e-05      0.000     -0.204      0.843      -0.001       0.001
wl_anom_tau0               0.2012      0.079      2.538      0.029       0.025       0.378
==============================================================================
Omnibus:                        1.149   Durbin-Watson:                   1.545
Prob(Omnibus):                  0.563   Jarque-Bera (JB):                0.974
Skew:                          -0.499   Prob(JB):                        0.615
Kurtosis:                       2.181   Cond. No.                         592.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Output Inference:

    Drop = y0 + y1.Wind + y2.Distance + y3.Water_anomaly + e

Where:

Drop = fraction decline in vessel activity during the storm
Wind = storm max wind (kt)
Distance = storm distance (km)
WaterAnomaly = tide gauge anomaly on event day

A positive coefficient means larger values of that predictor are associated with a larger drop_frac (bigger disruption).

Wind and distance are statistically not significant, whereas increase in water anomaly (stats. sig. = 0.029) by 1m leads to increase in drop by 20% (coef=0.2012).

---

Port-control sensitivity check

Because ports can differ in baseline activity and operational patterns, we rerun the regression with a simple port dummy (port_dummy = 1 for Houston, 0 for Norfolk). This checks whether the estimated sensitivity to hazards changes materially once we control for which port the event occurred at.

storm_event_metrics["port_dummy"] = (storm_event_metrics["port"] == "Houston").astype(int)

X = storm_event_metrics[
    ["storm_max_wind_kt_tau0",
     "wl_anom_tau0",
     "port_dummy"]
]
X = sm.add_constant(X)

y = storm_event_metrics["drop_frac"]

model = sm.OLS(y, X).fit()
print(model.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:              drop_frac   R-squared:                       0.513
Model:                            OLS   Adj. R-squared:                  0.367
Method:                 Least Squares   F-statistic:                     3.515
Date:                Tue, 24 Feb 2026   Prob (F-statistic):             0.0569
Time:                        00:25:52   Log-Likelihood:                 18.022
No. Observations:                  14   AIC:                            -28.04
Df Residuals:                      10   BIC:                            -25.49
Df Model:                           3
Covariance Type:            nonrobust
==========================================================================================
                             coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------------
const                      0.0526      0.051      1.024      0.330      -0.062       0.167
storm_max_wind_kt_tau0     0.0006      0.001      0.632      0.541      -0.001       0.003
wl_anom_tau0               0.2377      0.088      2.696      0.022       0.041       0.434
port_dummy                -0.0456      0.057     -0.797      0.444      -0.173       0.082
==============================================================================
Omnibus:                        1.445   Durbin-Watson:                   1.791
Prob(Omnibus):                  0.486   Jarque-Bera (JB):                0.977
Skew:                          -0.356   Prob(JB):                        0.613
Kurtosis:                       1.919   Cond. No.                         281.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Inference:

Addition of port dummy doesn't have strong impact on storm sensitivity.

---

Caveats:

• This analysis uses a small number of storm events (n ≈ 14), so uncertainty is high and p-values should be interpreted cautiously.
• The goal is exploratory: identify which variables are most informative for shock magnitude and motivate more robust models (e.g., more ports/years, robust SEs, or hierarchical models) in future work.

---

4.2) Daily panel Regression (Gradual stress response)

In contrast to the event-level “shock” model (where storm onset implied immediate drop), this section estimates a gradual-stress relationship using the full daily panel.

Model setup:

• Outcome: log_vessels = log(1 + unique_vessels) (log transform stabilizes variance and makes coefficients interpretable as approximate % changes).
• Key predictor: wl_anom (daily water-level anomaly).
• Port fixed effects: C(port) controls for baseline differences in activity between ports (e.g., Houston has higher typical vessel counts than Norfolk).

We fit three related specifications:

1. Base model with robust (HC3) standard errors:
   log_vessels ~ wl_anom + C(port)
2. Interaction model (port-specific sensitivity):
   log_vessels ~ wl_anom * C(port)
   tests whether the water-anomaly effect differs between Houston and Norfolk.
3. HAC-robust model (time-series autocorrelation):
   log_vessels ~ wl_anom + C(port) with HAC errors (7 lags)
   checks whether inference changes once we account for serial correlation in daily data.

panel

	date	port	unique_vessels	n_points	wl_anom	storm_day_200km	storm_min_dist_km	storm_max_wind_kt_200km	storm_pts_200km	log_unique_vessels	log_points
0	2019-08-30	Houston	439	228740	-0.000163	0	2473.236388	NaN	0	6.086775	12.340346
1	2019-08-31	Houston	442	236449	0.059942	0	2082.890406	NaN	0	6.093570	12.373492
2	2019-09-01	Houston	425	230437	0.097958	0	1786.073218	NaN	0	6.054439	12.347737
3	2019-09-02	Houston	448	228706	0.214625	0	1668.947950	NaN	0	6.107023	12.340197
4	2019-09-03	Houston	433	251193	0.281608	0	693.009619	NaN	0	6.073045	12.433981
...	...	...	...	...	...	...	...	...	...	...	...
411	2023-09-27	Norfolk	276	131754	0.376508	0	3031.664722	NaN	0	5.624018	11.788699
412	2023-09-28	Norfolk	286	131375	0.335575	0	2909.072264	NaN	0	5.659482	11.785819
413	2023-09-29	Norfolk	304	137144	0.226552	0	2911.190774	NaN	0	5.720312	11.828794
414	2023-09-30	Norfolk	278	118295	0.142190	0	2848.628526	NaN	0	5.631212	11.680945
415	2023-10-01	Norfolk	277	126435	0.149069	0	2284.706411	NaN	0	5.627621	11.747492

416 rows × 11 columns

panel = panel.copy()
panel["log_vessels"] = np.log(panel["unique_vessels"] + 1)

import statsmodels.formula.api as smf
model2 = smf.ols(
    "log_vessels ~ wl_anom + C(port)",
    data=panel
).fit(cov_type="HC3")

print(model2.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:            log_vessels   R-squared:                       0.858
Model:                            OLS   Adj. R-squared:                  0.857
Method:                 Least Squares   F-statistic:                     1266.
Date:                Tue, 24 Feb 2026   Prob (F-statistic):          6.46e-177
Time:                        00:26:05   Log-Likelihood:                 357.28
No. Observations:                 416   AIC:                            -708.6
Df Residuals:                     413   BIC:                            -696.5
Df Model:                           2
Covariance Type:                  HC3
======================================================================================
                         coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------
Intercept              6.0430      0.005   1138.560      0.000       6.033       6.053
C(port)[T.Norfolk]    -0.5029      0.010    -49.769      0.000      -0.523      -0.483
wl_anom               -0.1073      0.042     -2.542      0.011      -0.190      -0.025
==============================================================================
Omnibus:                        5.397   Durbin-Watson:                   0.647
Prob(Omnibus):                  0.067   Jarque-Bera (JB):                6.659
Skew:                          -0.112   Prob(JB):                       0.0358
Kurtosis:                       3.578   Cond. No.                         6.89
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity robust (HC3)

Model Interpretation:

log_vessels ~ wl_anom + C(port)

Vessel activity depends on water-level anomaly and port fixed effect.

For wl_anom, p-value = 0.011. Because dependent variable is log(vessels), a 1 meter increase in water level anomaly leads to ≈ 10.7% decrease in vessel activity.

So, water anomalies significantly reduce port activity (Statistically significant (p < 0.05)). Also, Norfolk baseline vessel counts are lower than Houston (port dummy coef is -0.5029).

---

model3 = smf.ols(
    "log_vessels ~ wl_anom * C(port)",
    data=panel
).fit(cov_type="HC3")

print(model3.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:            log_vessels   R-squared:                       0.858
Model:                            OLS   Adj. R-squared:                  0.857
Method:                 Least Squares   F-statistic:                     876.0
Date:                Tue, 24 Feb 2026   Prob (F-statistic):          2.38e-178
Time:                        00:26:13   Log-Likelihood:                 358.56
No. Observations:                 416   AIC:                            -709.1
Df Residuals:                     412   BIC:                            -693.0
Df Model:                           3
Covariance Type:                  HC3
==============================================================================================
                                 coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------------------
Intercept                      6.0444      0.005   1163.793      0.000       6.034       6.055
C(port)[T.Norfolk]            -0.5057      0.010    -50.742      0.000      -0.525      -0.486
wl_anom                       -0.1494      0.053     -2.820      0.005      -0.253      -0.046
wl_anom:C(port)[T.Norfolk]     0.0977      0.088      1.110      0.267      -0.075       0.270
==============================================================================
Omnibus:                        4.932   Durbin-Watson:                   0.652
Prob(Omnibus):                  0.085   Jarque-Bera (JB):                6.127
Skew:                          -0.086   Prob(JB):                       0.0467
Kurtosis:                       3.569   Cond. No.                         15.4
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity robust (HC3)

Model summary Interpretation:

log_vessels ~ wl_anom * C(port)

Now here we allow water impact to differ by port.

Base WL effect (Houston) = -0.1494

p = 0.005

So in Houston: 1m rise -> ~15% decrease

Interaction term (Norfolk adjustment) = +0.0977

p = 0.267 (NOT significant)

So Norfolk effect = -0.1494 + 0.0977 = -0.0517 (~5%)

But since interaction p-value > 0.05:

We cannot statistically say Norfolk behaves differently.

---

model_hac = smf.ols(
    "log_vessels ~ wl_anom + C(port)",
    data=panel
).fit(cov_type="HAC", cov_kwds={"maxlags":7})

print(model_hac.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:            log_vessels   R-squared:                       0.858
Model:                            OLS   Adj. R-squared:                  0.857
Method:                 Least Squares   F-statistic:                     277.9
Date:                Tue, 24 Feb 2026   Prob (F-statistic):           3.50e-77
Time:                        00:26:16   Log-Likelihood:                 357.28
No. Observations:                 416   AIC:                            -708.6
Df Residuals:                     413   BIC:                            -696.5
Df Model:                           2
Covariance Type:                  HAC
======================================================================================
                         coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------
Intercept              6.0430      0.010    595.349      0.000       6.023       6.063
C(port)[T.Norfolk]    -0.5029      0.022    -23.306      0.000      -0.545      -0.461
wl_anom               -0.1073      0.055     -1.938      0.053      -0.216       0.001
==============================================================================
Omnibus:                        5.397   Durbin-Watson:                   0.647
Prob(Omnibus):                  0.067   Jarque-Bera (JB):                6.659
Skew:                          -0.112   Prob(JB):                       0.0358
Kurtosis:                       3.578   Cond. No.                         6.89
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 7 lags and without small sample correction

HAC with 7 lags adjusts for time-series autocorrelation. So, effect size same (~10%) but significance slightly weaker (just above 5%).

A 1-meter positive water-level anomaly is associated with approximately 10–15% reduction in daily vessel activity across major U.S. ports, with no statistically significant difference in sensitivity between Houston and Norfolk.

---

Summary:

• In the base HC3 model, the wl_anom coefficient is negative and statistically significant (p ≈ 0.011), indicating that higher water-level anomalies are associated with lower daily AIS activity after controlling for port.
  • Because the outcome is log(1 + vessels), the coefficient can be interpreted approximately as a percent change: a +1 unit increase in wl_anom corresponds to about a 10–15% decrease in vessel activity (depending on the model).
• The port fixed effect confirms Norfolk has a lower baseline activity level than Houston (large negative coefficient on C(port)[T.Norfolk]).
• In the interaction model, the wl_anom × port interaction term is not statistically significant (p > 0.05), so we do not have evidence (in this sample) that Norfolk’s sensitivity to wl_anom differs from Houston’s.
• With HAC standard errors (7 lags), the estimated wl_anom effect size is similar but the p-value becomes weaker (≈0.053), reflecting the fact that daily time series can be autocorrelated and effective sample size is smaller than the raw row count.

Conclusion:

Daily water-level anomalies appear to act as a gradual stressor that is associated with reduced near-port vessel activity, while differences in baseline activity are largely captured by port fixed effects. Port-to-port differences in sensitivity are not statistically distinguishable here, especially after accounting for autocorrelation.

---

SECTION 5: Scaling up with Cyberinfrastructure (parallel resampling)

A single OLS/HAC regression fit is fast for the datasets used in this training notebook.
However, the number of storm events is small and single-fit coefficients can be sensitive to individual observations, so we use resampling (jackknife and bootstrap) to assess stability and uncertainty.

Resampling requires refitting the model many times; we then use joblib.Parallel to speed up these repeated refits (not the base regression itself). This demonstrates a portable CI pattern:

• split work into independent refits,
• run them across available CPU cores,
• collect coefficient distributions & compare runtimes (1 worker vs many workers).

5.1) Storm regression: parallel jackknife confidence intervals (LOO refits)

We revisit the storm “shock-response” regression where the outcome is drop_frac (immediate activity drop at storm onset).
Because we only have ~14 storm events, coefficients can be overly influenced by any single event. To check stability, we run a leave-one-out jackknife:

• for each i in 1…n, refit the OLS model after dropping event i,
• collect the parameter vector each time,
• summarize variability using empirical quantiles (2.5%, 50%, 97.5%).

Since each refit is independent, we can leverage parallelization by distributing the refits across CPU cores. We time the same computation with 1 worker vs n_jobs workers and report the speedup.

import statsmodels.api as sm
from joblib import Parallel, delayed
import warnings

def effective_n_jobs(max_jobs=8):
    return max(1, min(os.cpu_count() or 1, max_jobs))

def time_run(func, *args, **kwargs):
    t0 = time.perf_counter()
    out = func(*args, **kwargs)
    t1 = time.perf_counter()
    return out, (t1 - t0)

def fit_ols(df, y_col, x_cols):
    X = sm.add_constant(df[x_cols], has_constant="add")
    y = df[y_col].astype(float)
    return sm.OLS(y, X).fit()

def jackknife_params(df, y_col, x_cols, n_jobs=None):
    """
    Leave-one-out refits. Returns DataFrame of parameter estimates (n rows).
    """
    if n_jobs is None:
        n_jobs = effective_n_jobs()

    df = df.dropna(subset=[y_col] + x_cols).reset_index(drop=True)
    n = len(df)
    if n < (len(x_cols) + 2):
        raise ValueError(f"Too few rows for jackknife (n={n}).")

    def one(i):
        dfi = df.drop(index=i).reset_index(drop=True)
        m = fit_ols(dfi, y_col, x_cols)
        return m.params

    params = Parallel(n_jobs=n_jobs, backend="loky")(delayed(one)(i) for i in range(n))
    return pd.DataFrame(params)

# Regression setup
df_reg = storm_event_metrics.copy()
y_col = "drop_frac"
x_cols = ["storm_max_wind_kt_tau0", "storm_min_dist_km_tau0", "wl_anom_tau0"]

df_clean = df_reg.dropna(subset=[y_col] + x_cols).reset_index(drop=True)
print("n rows used in regression:", len(df_clean))

model_jk = fit_ols(df_clean, y_col, x_cols)
print(model_jk.summary())


# CI demo using Jackknife (parallel)

nj = effective_n_jobs(max_jobs=8)

jk1, t1 = time_run(jackknife_params, df_clean, y_col, x_cols, n_jobs=1)
jkp, tp = time_run(jackknife_params, df_clean, y_col, x_cols, n_jobs=nj)

print(f"\nJackknife refits = {len(df_clean)}")
print(f"1 worker:  {t1:.2f} sec")
print(f"{nj} workers: {tp:.2f} sec")
print(f"Speedup:   {t1/tp:.2f}×" if tp > 0 else "Speedup: N/A")

ci = jkp.quantile([0.025, 0.5, 0.975]).T
ci.columns = ["p2.5", "median", "p97.5"]
ci["ols_point_est"] = model.params.reindex(ci.index)

display(ci)

n rows used in regression: 14
                            OLS Regression Results
==============================================================================
Dep. Variable:              drop_frac   R-squared:                       0.485
Model:                            OLS   Adj. R-squared:                  0.330
Method:                 Least Squares   F-statistic:                     3.133
Date:                Mon, 23 Feb 2026   Prob (F-statistic):             0.0743
Time:                        14:01:08   Log-Likelihood:                 17.620
No. Observations:                  14   AIC:                            -27.24
Df Residuals:                      10   BIC:                            -24.68
Df Model:                           3
Covariance Type:            nonrobust
==========================================================================================
                             coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------------
const                      0.0648      0.064      1.020      0.332      -0.077       0.206
storm_max_wind_kt_tau0     0.0004      0.001      0.446      0.665      -0.002       0.003
storm_min_dist_km_tau0 -7.096e-05      0.000     -0.204      0.843      -0.001       0.001
wl_anom_tau0               0.2012      0.079      2.538      0.029       0.025       0.378
==============================================================================
Omnibus:                        1.149   Durbin-Watson:                   1.545
Prob(Omnibus):                  0.563   Jarque-Bera (JB):                0.974
Skew:                          -0.499   Prob(JB):                        0.615
Kurtosis:                       2.181   Cond. No.                         592.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Jackknife refits = 14
1 worker:  0.03 sec
8 workers: 0.01 sec
Speedup:   2.02×

	p2.5	median	p97.5	ols_point_est
const	-0.005948	0.065251	0.114150	0.052621
storm_max_wind_kt_tau0	0.000055	0.000386	0.001872	0.000590
storm_min_dist_km_tau0	-0.000354	-0.000068	0.000096	NaN
wl_anom_tau0	0.138627	0.199233	0.237536	0.237656

Inference:
Parallel execution yields a runtime improvement even for a small dataset of 14 events. In large-scale resilience studies across hundreds of ports and decades of storms, this scaling becomes essential.

The jackknife coefficient distribution provides a quick stability check. In this run, the wl_anom_tau0 coefficient remains consistently positive across refits, while wind and distance are closer to zero / less stable: supporting the earlier conclusion that surge context is the most informative predictor of storm-time drops in this sample.

---

5.2) Daily panel regression: parallel moving-block bootstrap for HAC models

Daily port activity is a time series, so observations are autocorrelated. To quantify uncertainty more robustly for the daily HAC regression, we use a moving-block bootstrap:

• within each port, sample contiguous blocks of days (e.g., 14-day blocks) with replacement,
• concatenate blocks until reaching the original series length per port,
• fit the same HAC regression on each bootstrap sample,
• repeat B times (e.g., B=400) to obtain coefficient distributions and confidence intervals.

This is computationally expensive because it involves hundreds of refits, so we parallelize the bootstrap iterations using joblib.

import statsmodels.formula.api as smf
from joblib import Parallel, delayed

def effective_n_jobs(max_jobs=8):
    return max(1, min(os.cpu_count() or 1, max_jobs))

def time_run(func, *args, **kwargs):
    t0 = time.perf_counter()
    out = func(*args, **kwargs)
    t1 = time.perf_counter()
    return out, (t1 - t0)

PANEL_GLOBAL = None
FORMULA_GLOBAL = None
HAC_LAGS_GLOBAL = None
BLOCK_DAYS_GLOBAL = None
SEED_GLOBAL = None

def _init_bootstrap_globals(panel, formula, hac_lags, block_days, seed):
    global PANEL_GLOBAL, FORMULA_GLOBAL, HAC_LAGS_GLOBAL, BLOCK_DAYS_GLOBAL, SEED_GLOBAL
    PANEL_GLOBAL = panel
    FORMULA_GLOBAL = formula
    HAC_LAGS_GLOBAL = hac_lags
    BLOCK_DAYS_GLOBAL = block_days
    SEED_GLOBAL = seed

def _one_bootstrap_sample(rng):
    """
    Moving-block bootstrap over time within each port:
    - sample contiguous blocks of days with replacement
    - concatenate until we reach original length per port
    """
    df = PANEL_GLOBAL
    block_days = BLOCK_DAYS_GLOBAL

    out_parts = []
    for port, g in df.groupby("port"):
        g = g.sort_values("date")
        dates = g["date"].values
        n = len(g)
        if n == 0:
            continue

        # start positions for blocks (indices)
        max_start = max(1, n - block_days)
        starts = rng.integers(0, max_start, size=int(np.ceil(n / block_days)))

        idxs = []
        for s in starts:
            idxs.extend(range(s, min(s + block_days, n)))
            if len(idxs) >= n:
                break
        idxs = idxs[:n]

        out_parts.append(g.iloc[idxs])

    boot = pd.concat(out_parts, ignore_index=True)
    return boot

def _fit_hac_on_bootstrap(boot_df):
    m = smf.ols(FORMULA_GLOBAL, data=boot_df).fit(
        cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS_GLOBAL}
    )
    return m.params, m.bse

def _run_bootstrap_chunk(chunk_id, B_chunk):
    rng = np.random.default_rng(SEED_GLOBAL + chunk_id)
    params_list = []
    bse_list = []
    for _ in range(B_chunk):
        boot = _one_bootstrap_sample(rng)
        params, bse = _fit_hac_on_bootstrap(boot)
        params_list.append(params)
        bse_list.append(bse)
    return pd.DataFrame(params_list), pd.DataFrame(bse_list)

def bootstrap_hac_parallel(panel, formula, hac_lags=7, block_days=14, B=200, n_jobs=None, seed=42):
    """
    Parallel moving-block bootstrap for HAC regression.
    Returns:
      params_df: B x coefficients
      bse_df:    B x coefficients (SEs)
    """
    if n_jobs is None:
        n_jobs = effective_n_jobs()

    # Date required for block bootstrap
    if "date" not in panel.columns:
        raise ValueError("panel must have a 'date' column for block bootstrap.")

    # initialize globals once in parent (workers will receive on spawn)
    _init_bootstrap_globals(panel, formula, hac_lags, block_days, seed)

    # chunking: each worker does many iterations
    total_chunks = n_jobs  # keep simple
    Bs = [B // total_chunks] * total_chunks
    for i in range(B - sum(Bs)):
        Bs[i] += 1

    results = Parallel(n_jobs=n_jobs, backend="loky")(
        delayed(_run_bootstrap_chunk)(cid, Bc)
        for cid, Bc in enumerate(Bs) if Bc > 0
    )

    params_parts = [r[0] for r in results]
    bse_parts = [r[1] for r in results]

    return pd.concat(params_parts, ignore_index=True), pd.concat(bse_parts, ignore_index=True)

# CI demo

formula = "log_vessels ~ wl_anom + C(port)"
B = 400
block_days = 14 # 2-week blocks (for daily autocorrelation)
hac_lags = 7
nj = effective_n_jobs(max_jobs=8)

# Base HAC model (single full-sample fit)
base_model = smf.ols(
    formula,
    data=panel
).fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags})

print("\n=== Base HAC Model Summary (Full Sample) ===\n")
print(base_model.summary())

# Sequential (1 worker)
(params1, bse1), t1 = time_run(
    bootstrap_hac_parallel, panel, formula,
    hac_lags=hac_lags, block_days=block_days, B=B, n_jobs=1, seed=42
)

# Parallel
(paramsp, bsep), tp = time_run(
    bootstrap_hac_parallel, panel, formula,
    hac_lags=hac_lags, block_days=block_days, B=B, n_jobs=nj, seed=42
)

print(f"Bootstrap-HAC (B={B}, block_days={block_days}, maxlags={hac_lags})")
print(f"1 worker:  {t1:.2f} sec")
print(f"{nj} workers: {tp:.2f} sec")
print(f"Speedup:   {t1/tp:.2f}×" if tp > 0 else "Speedup: N/A")

# Show uncertainty of wl_anom coefficient
ci = paramsp["wl_anom"].quantile([0.025, 0.5, 0.975]).to_frame().T
ci.index = ["wl_anom coef CI"]
ci.columns = ["p2.5", "median", "p97.5"]
display(ci)

=== Base HAC Model Summary (Full Sample) ===

                            OLS Regression Results
==============================================================================
Dep. Variable:            log_vessels   R-squared:                       0.858
Model:                            OLS   Adj. R-squared:                  0.857
Method:                 Least Squares   F-statistic:                     277.9
Date:                Mon, 23 Feb 2026   Prob (F-statistic):           3.50e-77
Time:                        14:22:22   Log-Likelihood:                 357.28
No. Observations:                 416   AIC:                            -708.6
Df Residuals:                     413   BIC:                            -696.5
Df Model:                           2
Covariance Type:                  HAC
======================================================================================
                         coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------
Intercept              6.0430      0.010    595.349      0.000       6.023       6.063
C(port)[T.Norfolk]    -0.5029      0.022    -23.306      0.000      -0.545      -0.461
wl_anom               -0.1073      0.055     -1.938      0.053      -0.216       0.001
==============================================================================
Omnibus:                        5.397   Durbin-Watson:                   0.647
Prob(Omnibus):                  0.067   Jarque-Bera (JB):                6.659
Skew:                          -0.112   Prob(JB):                       0.0358
Kurtosis:                       3.578   Cond. No.                         6.89
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 7 lags and without small sample correction
Bootstrap-HAC (B=400, block_days=14, maxlags=7)
1 worker:  2.28 sec
8 workers: 0.31 sec
Speedup:   7.31×

	p2.5	median	p97.5
wl_anom coef CI	-0.216477	-0.12749	-0.024843

Inference:
The parallel bootstrap shows substantial speedup (7.3x) with 8 workers; again illustrating the importance of CI once we scale to multiple ports/years. Moreover, we use HAC errors and a moving-block bootstrap that resamples contiguous day blocks, so uncertainty isn’t underestimated by treating days as independent.

The bootstrapped CI for the wl_anom coefficient remains negative [-0.22, -0.02], reinforcing the daily-panel result that higher water-level anomalies are associated with lower near-port activity in this sample.

---

Summary & Key Findings:

• Built a daily port × date dataset linking AIS activity (unique vessels) with hazards (water-level anomalies + storm exposure) and derived event-level resilience metrics (drop_frac, recovery_days_to_90pct).
• Using a pooled 95th percentile threshold, found 65 extreme water-level days total (Houston: 22, Norfolk: 43).
• Extreme water-level days showed small/near-zero disruption overall (median drop_frac: Houston ~0.023, Norfolk ~-0.013) with fast recovery (median = 0 days for both).
• Identified 14 storm events (Houston: 5, Norfolk: 9) with larger drops (median drop_frac: Houston ~0.168, Norfolk ~0.094; median recovery Houston 1 day, Norfolk 0 days).
• Storm-day disruption aligns most with water-level anomaly (surge proxy): corr(drop_frac, wl_anom_tau0) ~0.69, vs wind ~0.38, distance ~0.0. In event regression, wl_anom_tau0 is significant (coef ~0.201, p ~0.029).
• In the daily panel, higher wl_anom is associated with lower activity (wl_anom coef ~ -0.107; HC3 p ~0.011; bootstrap CI roughly [-0.216, -0.025]).
• Same-day matching produced no compound events: 65 wl_only vs 14 storm_only.

Future Work:

• Scale to more ports + more years and run as a CI pipeline (job arrays / Dask) for stable cross-port comparisons.
• Improve event definitions: storms as multi-day runs, and compound hazards using overlap windows (e.g., storm ±1–3 days with extreme WL).
• Add confounder controls in the daily model (seasonality, day-of-week, long-term trends, lagged activity).
• Expand hazard features (precip, discharge, waves) and validate disruptions with port closure / advisory data where available.