Metadata-Version: 2.4
Name: shapneurgc
Version: 0.1.0
Summary: Shapley Regularized Neural Granger Causality: global feature attribution for temporal causal discovery in multivariate time series
Author-email: Merwan Roudane <merwanroudane920@gmail.com>
Maintainer-email: Merwan Roudane <merwanroudane920@gmail.com>
License-Expression: MIT
Project-URL: Homepage, https://github.com/merwanroudane/shapneurgc
Project-URL: Repository, https://github.com/merwanroudane/shapneurgc
Project-URL: Researcher's guide, https://github.com/merwanroudane/shapneurgc/blob/main/docs/GUIDE.md
Project-URL: Syntax reference, https://github.com/merwanroudane/shapneurgc/blob/main/docs/SYNTAX.md
Project-URL: Theory to code map, https://github.com/merwanroudane/shapneurgc/blob/main/docs/THEORY.md
Project-URL: Issues, https://github.com/merwanroudane/shapneurgc/issues
Keywords: granger causality,causal discovery,time series,shapley value,feature attribution,neural networks,econometrics,macroeconomics
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Financial and Insurance Industry
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Artificial Intelligence
Classifier: Topic :: Scientific/Engineering :: Information Analysis
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: torch>=2.0
Requires-Dist: numpy>=1.22
Requires-Dist: pandas>=1.5
Requires-Dist: scikit-learn>=1.1
Requires-Dist: scipy>=1.9
Requires-Dist: matplotlib>=3.6
Provides-Extra: plots
Requires-Dist: networkx>=3.0; extra == "plots"
Provides-Extra: econ
Requires-Dist: statsmodels>=0.13; extra == "econ"
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: pytest-cov; extra == "dev"
Requires-Dist: ruff; extra == "dev"
Provides-Extra: all
Requires-Dist: networkx>=3.0; extra == "all"
Requires-Dist: statsmodels>=0.13; extra == "all"
Requires-Dist: pytest>=7.0; extra == "all"
Requires-Dist: tabulate; extra == "all"
Dynamic: license-file

# shapneurgc

**Shapley Regularized Neural Granger Causality — a complete, documented Python library for temporal causal discovery.**

[![Python](https://img.shields.io/badge/python-3.9%2B-blue)](https://www.python.org/)
[![PyTorch](https://img.shields.io/badge/PyTorch-2.0%2B-ee4c2c)](https://pytorch.org/)
[![License](https://img.shields.io/badge/license-MIT-green)](LICENSE)
[![Tests](https://img.shields.io/badge/tests-58%20passing-brightgreen)](tests/)

An implementation of the Information-Theoretic Shapley value (Info-Shap) for
Granger-causal discovery in multivariate time series, following Yang, Zhu, Tian,
Gao and Li (2026), *Shapley Regularized Neural Granger Causality*, ICML,
PMLR 306 — extended with the tooling an applied study actually needs:
publication-quality tables and figures, a validated edge-selection rule, a real
macroeconomic application, and a long-form guide.

### → **[Read the full researcher's guide](https://github.com/merwanroudane/shapneurgc/blob/main/docs/GUIDE.md)** ←

*A step-by-step walkthrough of how to write a complete SRNGC application, from
raw data to submission-ready tables — including the validation checks that
should pass before you interpret anything.*

---

## Contents

- [Why this method](#why-this-method)
- [Installation](#installation)
- [Sixty-second quickstart](#sixty-second-quickstart)
- [Detailed syntax](#detailed-syntax) — every call, with defaults
- [Worked examples](#worked-examples) — eight complete, runnable recipes
- [Documentation](#documentation)
- [Applied study: US monetary policy transmission](#applied-study-us-monetary-policy-transmission)
- [Benchmark results](#benchmark-results)
- [The selection rule matters more than the estimator](#the-selection-rule-matters-more-than-the-estimator)
- [Figures](#figures)
- [Example scripts](#example-scripts)
- [Differences from the reference implementation](#differences-from-the-reference-implementation)
- [Project layout](#project-layout)
- [Testing](#testing)
- [Citation](#citation)

---

## Why this method

Neural Granger causality rests on three interlocking pieces — a **predictive
model**, a **sparsity penalty**, and an **importance measure** that reads the
causal graph off the fitted model. The literature has invested heavily in the
first and treated the third as an afterthought, reading causality off
coefficient matrices, first-layer weights, or the input–output Jacobian.

Those proxies give only *necessary* conditions. For a self-explaining model the
true derivative is

```
∂X^i_t / ∂X^j_{t-k}  =  [g(X_{t-k})]_ij  +  Σ_m ∂[g(X_{t-k})]_im/∂X^j_{t-k} · X^m_{t-k}
                        └ the "interpretable" part ┘   └──── what it misses ────┘
```

so sparsity in the interpretable component does not imply the output is
invariant to the input. A Jacobian escapes the architectural constraint but is
still *local* and blind to feature interactions.

**Info-Shap** replaces the proxy with an axiomatic global measure:

```
φ_ij  =  E[J_ij²]  +  ½ Σ_{k≠j} E[|J_ij · J_ik|]
         └ individual ┘   └──── interaction ────┘
```

Two properties make it work. It is **zero exactly under Granger
non-causality** (all terms are non-negative and carry a factor `J_ij`), and it
has a **closed form** — because the underlying coalition game is 2-additive, the
exponential Shapley sum collapses to these two terms. That makes it cheap enough
to use as a training penalty *and* as the inference rule, aligning the two.

In this library the split is exposed, so you can see how much of an attribution
a local measure would have missed:

```python
model.individual_    # E[J²]
model.interaction_   # ½ Σ E|J J|
```

---

## Installation

```bash
pip install shapneurgc
```

With the optional extras (network figures, ADF tests, Markdown tables):

```bash
pip install "shapneurgc[all]"
```

From source, for development:

```bash
git clone https://github.com/merwanroudane/shapneurgc.git
cd shapneurgc
pip install -e ".[all]"
```

| extra | pulls in | needed for |
|---|---|---|
| *(core)* | `torch`, `numpy`, `pandas`, `scikit-learn`, `scipy`, `matplotlib` | everything except the items below |
| `plots` | `networkx` | `plot_network` |
| `econ` | `statsmodels` | `adf_report` |
| `dev` | `pytest`, `ruff` | the test suite |
| `all` | all of the above + `tabulate` | `tables.to_markdown` |

Python 3.9+. CPU is fine — every number in this README was produced on CPU.

The US macro panel **ships inside the package**, so this works immediately after
install, with no repository checkout:

```python
from shapneurgc import load_us_macro
panel = load_us_macro()          # MacroPanel(T=795, d=10, 1960-02-01 to 2026-06-01)
```

Refresh it from FRED at any time (no API key needed):

```bash
python -m shapneurgc.fetch_fred --start 1960-01-01
```

---

## Sixty-second quickstart

```python
import numpy as np
from shapneurgc import SRNGC, make_var
from shapneurgc.metrics import binary_scores

# a sparse stationary VAR(3) whose true graph we know
sim = make_var(dim=8, T=800, sparsity=0.3, max_lag=3, seed=7)

model = SRNGC(lag=3, penalty="shap", hidden_dim=64, layers=2,
              lam_ind=1e-3, max_epochs=400, patience=30).fit(sim.X)

print(model.score(sim.graph))        # {'auroc': 1.0, 'auprc': 1.0, ...}
graph = model.discover(q=0.10)       # binary [d, d] graph
print(model.summary(top=10))
```

```
==========================================================================
  Shapley Regularized Neural Granger Causality
==========================================================================
  Backbone        : ResidualMLP (27,528 parameters)
  Penalty         : shap   lam_ind=0.001  lam_int=0
  Series          : d = 8,  lag K = 3,  effective N = 797
  Validation MSE  : 0.395820
  Info-Shap split : individual 3.866, interaction 15.25 (79.8% interaction)
  Selection       : knockoff+ at q = 0.1, threshold 0.05739
  Edges selected  : 22 of 64 (density 0.344)
--------------------------------------------------------------------------
  Strongest 10 lagged effects
--------------------------------------------------------------------------
         x8 -> x8       lag 1   1.02110  *   [interaction 49%]
         x2 -> x2       lag 1   0.70177  *   [interaction 63%]
         x7 -> x7       lag 1   0.62489  *   [interaction 67%]
         x6 -> x2       lag 1   0.56572  *   [interaction 68%]
         x8 -> x7       lag 2   0.43321  *   [interaction 74%]
   (* = retained by the knockoff filter)
==========================================================================
```

`python examples/01_quickstart.py` runs exactly this and writes the figures.

---

## Detailed syntax

Complete signatures with every default. Conventions used throughout:
`graph[i, j]` means **source `j` → target `i`** (row = effect, column = cause);
the design vector is `[X_{t-1}, …, X_{t-K}]`, so series `j` at lag `k` sits at
flat index `(k-1)·d + j`; a trailing underscore marks an attribute learned
during `fit`.

### The estimator

```python
from shapneurgc import SRNGC

model = SRNGC(
    lag           = 3,             # int   K, lag order
    penalty       = "shap",        # 'shap'|'fshap'|'jacob_l1'|'jacob_f'|'layer_weight'|'none'
    model_type    = "ResidualMLP", # 'ResidualMLP'|'MLP'|'cMLP'|'LSTM'|'cLSTM'
    hidden_dim    = 64,            # int
    layers        = 2,             # int   0 allowed for ResidualMLP (linear map)
    dropout       = 0.0,           # float
    lam_ind       = 1e-3,          # float weight on the individual term
    lam_int       = 0.0,           # float weight on the interaction term (fshap only)
    num_proj      = 1,             # int   R random projections; -1 = exact
    lr            = 1e-3,          # float
    weight_decay  = 1e-5,          # float
    batch_size    = -1,            # int   -1 = full batch
    max_epochs    = 2000,          # int
    patience      = 50,            # int   early-stopping patience
    standardize   = True,          # bool
    train_frac    = 0.7,           # float chronological split
    sync_dropout  = True,          # bool  hold the dropout mask across projections
    summarise     = "max",         # 'max'|'sum'|'mean' collapse over lags
    seed          = 2025,          # int
    device        = "cpu",         # 'cpu'|'cuda'|'cuda:0'
    verbose       = 1,             # 0 silent, 1 periodic, 2 every epoch
)
```

> `penalty="shap"` applies **one** λ to `individual + interaction` (the paper's
> convention), so `lam_int` is ignored there. `penalty="fshap"` uses both.

**Methods**

```python
model.fit(X, names=None)                    # X: [T, d] or [B, T, d] replicates
model.discover(q=0.1, offset=1, per_lag=False, refit=True,
               method="block", block=None)  # -> [d, d] binary graph
model.score(truth, ignore_diagonal=False)   # -> dict of AUROC/AUPRC/...
model.to_frame(top=None, include_lags=True) # -> tidy edge DataFrame
model.summary(top=15)                       # -> printable report string
model.graph_frame(binary=False)             # -> labelled [d, d] DataFrame
model.predict(X=None)                       # -> one-step-ahead forecasts
```

**Attributes after `fit`**

| attribute | shape | meaning |
|---|---|---|
| `importance_` | `[d, K·d]` | raw Info-Shap scores |
| `importance_by_lag_` | `[K, d, d]` | reshaped by lag |
| `summary_graph_` | `[d, d]` | lag-collapsed continuous scores |
| `individual_` | `[d, K·d]` | individual component `E[J²]` |
| `interaction_` | `[d, K·d]` | interaction component |
| `val_loss_` | float | best unpenalised validation MSE |
| `history_` | dict | `train`, `val`, `penalty` per epoch |
| `knockoff_` | object | populated by `.discover()` |
| `model_`, `dataset_`, `names_` | — | the fitted network, data, labels |

### Choosing the penalty

| `penalty` | measure | estimation | cost | pick it when |
|---|---|---|---|---|
| `"shap"` | global | exact | `O(d·C_g)` | default; small-to-moderate `d`, small samples |
| `"fshap"` | global | stochastic | `O(R·C_g)` | `d` large, or exact Jacobians too slow |
| `"jacob_l1"` | local | exact | `O(d·C_g)` | JRNGC baseline |
| `"jacob_f"` | local | stochastic | `O(R·C_g)` | JRNGC baseline, large `d` |
| `"none"` | — | — | — | reference point; expect AUPRC to collapse |

### Lower-level building blocks

```python
# ---- data ----------------------------------------------------------------
from shapneurgc import TimeSeriesDataset, lag_index, unflatten_importance, \
                       summarise_over_lags

ds = TimeSeriesDataset(X, lag=3, standardize=True, device="cpu")
train, val = ds.split(train_frac=0.7, shuffle=False)   # chronological
ds.inputs, ds.outputs        # [N, K*d], [N, d]
lag_index(lag_k=2, series_j=1, dim=10)                 # -> flat column index
unflatten_importance(phi, dim=10, lag=3)               # [d, K*d] -> [K, d, d]
summarise_over_lags(by_lag, how="max")                 # [K, d, d] -> [d, d]

# ---- measures ------------------------------------------------------------
from shapneurgc import InfoShapley, JacobianMeasure, full_jacobian

J   = full_jacobian(y, x, create_graph=False)          # [batch, d_out, d_in]
phi = InfoShapley().compute(model, x)                  # [d_out, d_in]
phi, ind, inter = InfoShapley(return_parts=True).compute(model, x)
loc = JacobianMeasure().compute(model, x)              # E|J_ij|, JRNGC baseline

# ---- penalties -----------------------------------------------------------
from shapneurgc import ShapleyPenalty, build_penalty

pen = ShapleyPenalty(kind="shap", num_proj=1, sync_dropout=True,
                     reuse_jvp=True, dist="rademacher", device="cpu")
individual, interaction = pen(net, x)                  # two scalar tensors
build_penalty("Fast_Shap")            # paper names work as aliases

# ---- models --------------------------------------------------------------
from shapneurgc import build_model, MODEL_NAMES, count_parameters

net = build_model("ResidualMLP", dim=10, lag=3, hidden_dim=40,
                  layers=2, dropout=0.0, device="cpu")   # [B, K*d] -> [B, d]

# ---- training ------------------------------------------------------------
from shapneurgc import TrainConfig, train_model, set_seed

res = train_model(net, train, val, pen,
                  TrainConfig(lam_ind=1e-3, tie_lambdas=True, lr=1e-3,
                              max_epochs=2000, patience=50))
res.model, res.best_val_loss, res.history, res.epochs_run

# ---- selection -----------------------------------------------------------
from shapneurgc import grid_search, ParamGrid, default_grid, select_lag

select_lag(X, candidate_lags=(1, 2, 3, 4, 6))          # -> DataFrame, best first
grid_search(ds, penalty_name="shap", grid=ParamGrid({...}),
            select_by="val_loss",                       # never on a graph metric
            num_workers=1, worker_index=1)

# ---- edge selection ------------------------------------------------------
from shapneurgc import select_graph, knockoff_threshold, make_knockoff_series
from shapneurgc import stability_selection, expected_false_positives

stability_selection(X, fit_fn, n_subsamples=30, subsample_frac=0.75,
                    top_m=None, density=0.25, pi=0.7,
                    exclude_diagonal=False, seed=0)

# ---- metrics -------------------------------------------------------------
from shapneurgc import graph_scores, binary_scores, confusion_map

graph_scores(truth, scores, ignore_diagonal=False)     # .auroc .auprc
binary_scores(truth, graph)                            # .precision .recall .f1 .fdr .shd

# ---- simulation ----------------------------------------------------------
from shapneurgc import make_var, make_lorenz96, VAR3Generator, Lorenz96Generator

make_var(dim=10, T=1000, sparsity=0.3, max_lag=3, seed=0)
make_lorenz96(dim=20, T=1000, F=10.0, seed=0)

# ---- economic data -------------------------------------------------------
from shapneurgc import load_us_macro, fetch_fred, adf_report

load_us_macro(path=None, start=None, end=None, columns=None, drop_covid=False)
fetch_fred(["INDPRO", "UNRATE"], start="1990-01-01")
adf_report(panel.data)
```

Every signature, argument and return type is in
[SYNTAX.md](docs/SYNTAX.md).

---

## Worked examples

Eight self-contained recipes. Each runs as written.

### 1. Minimal fit on your own data

```python
import numpy as np, pandas as pd
from shapneurgc import SRNGC

df = pd.read_csv("mydata.csv", parse_dates=["date"]).set_index("date").dropna()
X, names = df.to_numpy(float), list(df.columns)

model = SRNGC(lag=3, penalty="shap", lam_ind=1e-3).fit(X, names=names)
print(model.summary(top=20))
print(model.graph_frame().round(4))       # labelled [target, source] matrix
```

### 2. Choose the lag order before anything else

```python
from shapneurgc.selection import select_lag

tbl = select_lag(X, candidate_lags=(1, 2, 3, 4, 6, 12),
                 hidden_dim=64, layers=2, max_epochs=400)
print(tbl)              # sorted best-first by validation MSE
K = int(tbl.iloc[0]["lag"])
```

Selected with **no penalty**: the lag order is a property of the data, not of
the regularizer.

### 3. Tune hyperparameters honestly

```python
from shapneurgc.selection import grid_search, ParamGrid

grid = ParamGrid({
    "lr":           [5e-4, 1e-3, 5e-3],
    "hidden_dim":   [2 * d, 4 * d],
    "layers":       [1, 2, 3],
    "dropout":      [0.0, 0.1],
    "lam_ind":      [1e-4, 1e-3, 1e-2, 1e-1],
    "lam_int":      [0.0],
    "weight_decay": [1e-5],
})

res = grid_search(ds, penalty_name="shap", grid=grid, select_by="val_loss")
print(res.best)         # winning configuration
print(res.top(10))      # the ten best rows
```

Selection is on **validation MSE**. Passing `select_by="auroc"` still works for
reproducing published numbers but requires a ground-truth graph and warns that
it is oracle tuning.

### 4. Turn scores into a graph — the rule that calibrates

```python
from shapneurgc.stability import stability_selection

def fit_fn(Xs):
    m = SRNGC(lag=K, penalty="shap", hidden_dim=4 * d, layers=2,
              lam_ind=1e-3, max_epochs=300, patience=25,
              seed=7, verbose=0).fit(Xs)
    return m.summary_graph_, m.importance_by_lag_

stab = stability_selection(X, fit_fn, n_subsamples=30, subsample_frac=0.75,
                           density=0.25, pi=0.7)
graph = stab.graph
print(stab)                              # includes the Meinshausen-Bühlmann bound
print(stab.top_edges(names, top=20))     # ranked by selection frequency
```

For the knockoff filter of the original paper — and a check on whether it is
usable on your data:

```python
graph_ko = model.discover(q=0.10, offset=1)     # offset=1 -> knockoff+
kr = model.knockoff_
print(f"{(kr.W < 0).sum()} of {kr.W.size} statistics negative")
# roughly half should be. Far fewer means the threshold is not trustworthy.
```

### 5. Compare the global measure against the local baseline

```python
from scipy.stats import spearmanr

glob = SRNGC(lag=K, penalty="shap",     lam_ind=1e-3, verbose=0).fit(X)
loc  = SRNGC(lag=K, penalty="jacob_l1", lam_ind=1e-3, verbose=0).fit(X)

rho = spearmanr(glob.summary_graph_.ravel(), loc.summary_graph_.ravel()).statistic
share = glob.interaction_.sum() / (glob.individual_.sum() + glob.interaction_.sum())
print(f"Spearman(global, local) = {rho:.3f} | interaction share = {share:.1%}")
```

A high interaction share with a low correlation is the case where the global
measure earns its cost. Report both numbers either way.

### 6. Validate the pipeline on a graph you already know

```python
from shapneurgc import make_var
from shapneurgc.metrics import graph_scores, binary_scores

sim = make_var(dim=d, T=len(X), sparsity=0.3, max_lag=K, seed=0)
chk = SRNGC(lag=K, penalty="shap", lam_ind=1e-3, verbose=0).fit(sim.X)
print(graph_scores(sim.graph, chk.summary_graph_))
assert graph_scores(sim.graph, chk.summary_graph_).auroc > 0.8
```

If it cannot recover a known graph at your `d` and `T`, it will not recover an
unknown one.

### 7. Publication figures and tables

```python
from shapneurgc.plots import (set_style, save_figure, plot_importance_heatmap,
                              plot_lag_panels, plot_network, plot_decomposition,
                              plot_error_map, plot_roc_pr)
from shapneurgc import tables

set_style()                                   # once per script

save_figure(plot_importance_heatmap(model.summary_graph_, names),
            "fig_heatmap")                    # writes fig_heatmap.pdf and .png
save_figure(plot_lag_panels(model.importance_by_lag_, names,
                            summary=model.summary_graph_), "fig_lags")
save_figure(plot_network(graph, names, weights=model.summary_graph_,
                         groups=groups), "fig_network")
save_figure(plot_decomposition(model.individual_, model.interaction_,
                               names, dim=d, lag=K), "fig_decomp")

open("tab_edges.tex", "w").write(
    tables.edge_table(model.to_frame(), top=20,
                      caption="Strongest Granger-causal effects.",
                      label="tab:edges"))
```

All LaTeX output is `booktabs`; add `\usepackage{booktabs}` and
`\usepackage{multirow}`.

### 8. Bring your own network

```python
import torch.nn as nn
from shapneurgc import TimeSeriesDataset, ShapleyPenalty, TrainConfig, train_model
from shapneurgc.measures import InfoShapley

class MyNet(nn.Module):
    def __init__(self, in_dim, out_dim, hidden=128):
        super().__init__()
        self.f = nn.Sequential(nn.Linear(in_dim, hidden), nn.GELU(),
                               nn.Linear(hidden, hidden), nn.GELU(),
                               nn.Linear(hidden, out_dim))
    def forward(self, x):
        return self.f(x)

ds = TimeSeriesDataset(X, lag=3)
tr, va = ds.split(0.7)
net = MyNet(ds.input_dim, ds.output_dim)
res = train_model(net, tr, va, ShapleyPenalty("shap"),
                  TrainConfig(lam_ind=1e-3, tie_lambdas=True))
phi = InfoShapley().compute(res.model, ds.full_tensors()[0])
```

Anything mapping `[B, K·d] → [B, d]` and twice differentiable works — that is
the model-agnostic claim in practice. Avoid non-differentiable operations
(hard thresholds, `argmax`): they zero the Jacobian and with it the attribution.

---

## Documentation

| document | what it covers |
|---|---|
| **[GUIDE.md](docs/GUIDE.md)** | **Start here.** Writing a full application end to end: data preparation, lag choice, backbone, penalty, honest tuning, selection, validation, reporting, interpretation, a complete copy-paste study, troubleshooting, and a pre-submission checklist. |
| [SYNTAX.md](docs/SYNTAX.md) | Every public signature, argument, default and return type. |
| [THEORY.md](docs/THEORY.md) | Each paper equation mapped to the line of code that implements it, with the numerical verification. |
| [REPRODUCIBILITY.md](docs/REPRODUCIBILITY.md) | Where this library departs from the authors' reference release, and why — with measurements. |

---

## Applied study: US monetary policy transmission

`examples/03_monetary_policy.py` is a complete study on **real data**: ten
monthly FRED series, February 1960 to June 2026 (`T = 795`, `d = 10`), reduced
to stationary form with McCracken–Ng transformation codes.

| label | series | transform | group |
|---|---|---|---|
| `IP` | Industrial production | `100·Δlog` | Real activity |
| `EMP` | Nonfarm payroll employment | `100·Δlog` | Real activity |
| `UNRATE` | Unemployment rate | `Δ` | Real activity |
| `CPI` | Consumer price index | `100·Δlog` | Prices |
| `PPI` | Producer price index | `100·Δlog` | Prices |
| `M2` | M2 money stock | `100·Δlog` | Money |
| `FFR` | Effective federal funds rate | `Δ` | Policy & rates |
| `GS10` | 10-year Treasury yield | `Δ` | Policy & rates |
| `HOUST` | Housing starts | `100·Δlog` | Housing |
| `CREDSPR` | Baa–Aaa credit spread | level | Credit |

The script runs the whole pipeline: descriptive statistics, ADF tests, lag
selection by validation loss, the penalised fit, a JRNGC baseline for
comparison, both selection rules, a post-1985 subsample check, and every table
and figure.

```bash
python examples/03_monetary_policy.py
```

### What it finds

Lag order `K = 2`, chosen by validation MSE. **81.3%** of total Info-Shap mass
sits in the *interaction* term — this is a system where a purely local measure
would rank edges differently, and the strongest cross-variable attributions are
80–85% interaction-driven.

Own-lag dynamics dominate the raw scores, as they should in macro data, so the
script also ranks cross-variable channels separately (`exclude_diagonal=True`)
and reports each channel's **selection frequency** across 20 subsamples rather
than a bare yes/no.

Cross-variable edges surviving at `π = 0.7`:

| edge | frequency |
|---|---|
| `GS10 → FFR` | 1.00 |
| `GS10 → M2` | 1.00 |
| `HOUST → GS10` | 0.95 |
| `FFR → M2` | 0.85 |
| `UNRATE → M2` | 0.80 |

Textbook transmission channels, with their frequencies:

| channel | Info-Shap | frequency | selected |
|---|---|---|---|
| unemployment → policy rate | 0.164 | 0.35 | no |
| policy rate → consumer prices | 0.098 | 0.30 | no |
| inflation → policy rate | 0.084 | 0.25 | no |
| credit spread → employment | 0.050 | 0.65 | no |
| policy rate → real activity | 0.030 | 0.00 | no |

**None of the classic channels survives a strict stability threshold.** The
robust structure at monthly frequency is concentrated in the rates-and-money
block; the macro-transmission edges rank respectably on raw Info-Shap but are
not stable across subsamples. That is the honest result, and it is the kind of
result this library is built to let you report rather than paper over. Two
supporting diagnostics: the knockoff filter selects 65 of 100 edges here with
only 41 of 200 statistics negative — the degenerate regime described below —
and the importance matrix correlates 0.66 (Spearman) between the full sample
and the post-1985 subsample, so some of the instability is genuine parameter
drift rather than estimation noise.

The global and local measures agree closely on this dataset
(Spearman 0.937), which is itself worth reporting: here the interaction term
reorders edges without overturning the picture.

> **Granger causality is incremental predictability, not structural causation.**
> An edge from the policy rate to output says the rate helps forecast output
> beyond output's own past. It does not identify a monetary policy shock.
> Anticipation, omitted variables and temporal aggregation can all generate
> edges with no structural reading. `GUIDE.md` §12 spells out what you may and
> may not claim.

---

## Benchmark results

All numbers below were produced by the scripts in `examples/`, on CPU, and are
reproducible with `python examples/02_simulation_var.py`.

### VAR(3), `d = 10`, 3 replicates, fresh graph per replicate

| Metric | Penalty | T = 500 | T = 1000 |
|---|---|---|---|
| AUROC | **Shap** | **0.980** (0.014) | **0.995** (0.006) |
| | Jacob-ℓ₁ | 0.974 (0.018) | 0.988 (0.008) |
| | Jacob-F | 0.966 (0.021) | 0.977 (0.012) |
| | F-Shap | 0.958 (0.025) | 0.970 (0.012) |
| AUPRC | **Shap** | **0.970** (0.019) | **0.992** (0.009) |
| | Jacob-ℓ₁ | 0.958 (0.025) | 0.978 (0.015) |
| | Jacob-F | 0.941 (0.032) | 0.962 (0.022) |
| | F-Shap | 0.935 (0.033) | 0.955 (0.021) |

A VAR is *linear*, so its Jacobian characterises the dependence structure
completely and the four penalties should be close — they are. This benchmark is
a wiring test, not a demonstration of superiority. The exact ℓ₁ forms edge out
the random-projection forms, consistent with the paper.

### Lorenz-96, `d = 20`, `T = 1000`, 3 replicates

This is the regime the method is built for: bilinear chaotic dynamics, exactly
four parents per coordinate, and a local Jacobian that is a poor summary of
dependence averaged over the attractor. Settings follow Appendix D.1 of the
paper (`dt = 0.01`, observation noise sd 1.0).

| Metric | Penalty | F = 10 (weakly chaotic) | F = 40 (strongly chaotic) |
|---|---|---|---|
| AUROC | **Shap** | **0.722** (0.015) | **0.988** (0.011) |
| | Jacob-ℓ₁ | 0.719 (0.031) | 0.981 (0.017) |
| | Jacob-F | 0.671 (0.026) | 0.748 (0.017) |
| | F-Shap | 0.670 (0.024) | 0.744 (0.018) |
| AUPRC | **Shap** | **0.553** (0.022) | **0.977** (0.019) |
| | Jacob-ℓ₁ | 0.547 (0.020) | 0.964 (0.027) |
| | Jacob-F | 0.512 (0.009) | 0.553 (0.004) |
| | F-Shap | 0.511 (0.010) | 0.553 (0.006) |

Two things to read here, and the script prints both directly.

**The global measure leads the local one, and the gap widens with the
nonlinearity.** Holding the estimation scheme fixed, Shap − Jacob-ℓ₁ is
`+0.003` AUROC / `+0.006` AUPRC at `F = 10` and `+0.007` / `+0.013` at
`F = 40`. Small, but consistent in sign across both metrics and both regimes,
and moving in the direction the theory predicts: stronger forcing means more
of the dependence lives in the interaction term.

**The exact ℓ₁ forms dominate the random-projection forms by far more than the
global-vs-local gap** — 0.988 against 0.748 AUROC at `F = 40`. On this problem
the estimator matters more than the measure, which is worth knowing before
reaching for `fshap` on a hard nonlinear system.

---

## The selection rule matters more than the estimator

Turning continuous scores into a 0/1 graph is the step most often reported
without validation — and on dependent data the knockoff filter of the original
paper does not calibrate.

**Diagnosis** (VAR(3), `d=10`, `T=1000`, 300 candidate triples): the knockoff
copies receive only **31%** of the attribution the originals do, so
`W = φ̂ − φ̂^ko` is positive almost everywhere. Only **25 of 300** statistics
come out negative, against roughly 125 under exchangeability. The filter
estimates its false-discovery count from that tail, so the threshold collapses:

| nominal `q` | selected | precision | recall | **realised FDP** |
|---|---|---|---|---|
| 0.05 | 93 | 0.323 | 1.000 | 0.677 |
| 0.10 | 100 | 0.300 | 1.000 | 0.700 |
| 0.20 | 100 | 0.300 | 1.000 | 0.700 |
| 0.30 | 100 | 0.300 | 1.000 | 0.700 |

Precision 0.300 is exactly the base rate — every one of the 100 possible edges
was selected. The cause is structural: a knockoff copy is decoupled from the
targets, but a *null-yet-correlated* original is not, so a flexible model gives
it more attribution than the knockoff either way.

**`shapneurgc` therefore ships `stability_selection`** — refit on many
contiguous subsamples, keep the edges that survive in at least a fraction `π` of
them. Same data, same model, same importance measure:

| `π` | selected | precision | recall | **realised FDP** |
|---|---|---|---|---|
| 0.6 | 29 | 0.966 | 0.933 | **0.034** |
| 0.7 | 29 | 0.966 | 0.933 | **0.034** |
| 0.8 | 27 | 1.000 | 0.900 | **0.000** |
| 0.9 | 27 | 1.000 | 0.900 | **0.000** |

(true graph: 30 edges). FDP `0.034` against `0.700`.

Both rules are implemented, `examples/02` runs them side by side, and the
knockoff threshold emits a warning whenever its negative tail is too thin to
trust. Details in [REPRODUCIBILITY.md §5](docs/REPRODUCIBILITY.md).

---

## Figures

Ten publication-ready figure types, all vector (PDF) plus raster preview,
serif-set, Parula-coloured, and colour-vision-deficiency safe.

| figure | function |
|---|---|
| Importance heatmap | `plot_importance_heatmap` |
| Per-lag panels + summary | `plot_lag_panels` |
| TP/TN/FP/FN confusion map | `plot_error_map` |
| Original vs knockoff attributions | `plot_knockoff_scores` |
| Knockoff statistic and threshold | `plot_w_statistics` |
| ROC and precision–recall curves | `plot_roc_pr` |
| Directed network, grouped nodes | `plot_network` |
| Loss and penalty trajectories | `plot_training_history` |
| Individual vs interaction split | `plot_decomposition` |
| Grouped bars with error bars | `plot_method_comparison` |

```python
from shapneurgc.plots import set_style, plot_lag_panels, save_figure
set_style()                                   # once per script
fig = plot_lag_panels(model.importance_by_lag_, names,
                      summary=model.summary_graph_)
save_figure(fig, "results/figures/lags")      # writes .pdf and .png
```

Tables are emitted as `booktabs` LaTeX and as Markdown from the same object:

```python
from shapneurgc import tables
open("tab.tex", "w").write(
    tables.comparison_table(df, metrics=("auroc", "auprc"),
                            caption="Performance.", label="tab:perf"))
```

---

## Example scripts

| script | what it does | runtime (CPU) |
|---|---|---|
| [`01_quickstart.py`](examples/01_quickstart.py) | fit → attribute → select → plot, in 30 lines | ~1 min |
| [`02_simulation_var.py`](examples/02_simulation_var.py) | VAR(3) benchmark, 4 penalties × 2 sample sizes; knockoff vs stability calibration | ~5 min |
| [`03_monetary_policy.py`](examples/03_monetary_policy.py) | full applied study on real FRED data | ~5 min |
| [`04_lorenz96_nonlinear.py`](examples/04_lorenz96_nonlinear.py) | chaotic nonlinear system; where global beats local | ~13 min |

Every script takes `--quick` for a fast smoke run.

---

## Differences from the reference implementation

The estimators are numerically identical to the authors' release — verified
term by term in `tests/test_equations.py`. The differences are in the
scaffolding, and each is measured in
[REPRODUCIBILITY.md](docs/REPRODUCIBILITY.md).

| topic | reference | here |
|---|---|---|
| Model selection | 3 of 4 aggregation notebooks select on the test metric | validation MSE by default; oracle mode warns |
| Dropout + random projections | each JVP draws a fresh mask, breaking Eq. (51) | RNG synchronised across projections |
| Architecture dispatch | `TSMixer`/`SimpleMamba` silently build a 0-layer linear MLP | every name builds a distinct model; unknown names raise |
| Individual-term estimator | separate VJP (2 JVP + 1 VJP) | reuses the JVPs per Algorithm 1; both available |
| Exact F-Shap | reachable | reachable via `num_proj=-1`, verified against Eq. (30) |
| Knockoff threshold | §3.4 form | `offset` parameter; default knockoff+ per Algorithm 2 |
| Edge selection | knockoff only | knockoff **and** validated stability selection |
| Train/validation split | by index | chronological, documented |

---

## Project layout

```
shapneurgc/
├── shapneurgc/
│   ├── measures.py     Info-Shap and the Jacobian baseline
│   ├── penalties.py    Shap / F-Shap / Jacob-ℓ₁ / Jacob-F
│   ├── models.py       Residual MLP, MLP, LSTM (joint and component-wise)
│   ├── data.py         lagged design construction, scaling, splitting
│   ├── train.py        penalised training with early stopping
│   ├── selection.py    grid search and lag choice on validation loss
│   ├── knockoff.py     residual-bootstrap and block knockoff filter
│   ├── stability.py    stability selection
│   ├── metrics.py      AUROC/AUPRC and thresholded graph metrics
│   ├── simulate.py     VAR(3) and Lorenz-96 generators
│   ├── economics.py    FRED loaders, transformation codes, ADF
│   ├── tables.py       booktabs LaTeX and Markdown builders
│   ├── plots.py        ten publication figure types
│   ├── api.py          the SRNGC estimator
│   ├── fetch_fred.py   refresh the bundled panel from FRED
│   └── datasets/       bundled FRED panel, shipped with the wheel
├── docs/               GUIDE, SYNTAX, THEORY, REPRODUCIBILITY
├── examples/           four runnable studies
├── tests/              58 tests
└── results/            figures and tables written by the examples
```

---

## Testing

```bash
pytest -q                       # 58 tests
pytest tests/test_equations.py  # numerical equivalence to the paper
```

`test_equations.py` reconstructs every formula from its definition — explicit
loops where that is clearest — and compares against the optimised code:
Info-Shap against Eq. (28) to `7.5e-09`, the exact penalty against Eq. (29) and
against `Σφ`, F-Shap against Eq. (30), and both stochastic estimators for
unbiasedness. `test_pipeline.py` covers the lag layout, replicate isolation,
model dispatch, metrics, simulators, both selection rules and the API.

---

## Citation

If you use this library, please cite the method paper:

```bibtex
@inproceedings{yang2026shapley,
  title     = {Shapley Regularized Neural Granger Causality},
  author    = {Yang, Maolin and Zhu, Zhoufan and Tian, Yuanhe and
               Gao, Kun and Li, Muyi},
  booktitle = {Proceedings of the 43rd International Conference on Machine Learning},
  series    = {Proceedings of Machine Learning Research},
  volume    = {306},
  year      = {2026}
}
```

and, if the software itself was useful:

```bibtex
@software{roudane2026shapneurgc,
  author  = {Roudane, Merwan},
  title   = {shapneurgc: Shapley Regularized Neural Granger Causality in Python},
  year    = {2026},
  version = {0.1.0},
  url     = {https://github.com/merwanroudane/shapneurgc}
}
```

---

## Author

**Dr Merwan Roudane** — [github.com/merwanroudane](https://github.com/merwanroudane)
· merwanroudane920@gmail.com

Released under the [MIT License](LICENSE). The method is due to Yang et al.
(2026); this is an independent implementation with additional tooling, and any
errors in it are mine.
