shapneurgc

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

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 ←

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
• Installation
• Sixty-second quickstart
• Detailed syntax — every call, with defaults
• Worked examples — eight complete, runnable recipes
• Documentation
• Applied study: US monetary policy transmission
• Benchmark results
• The selection rule matters more than the estimator
• Figures
• Example scripts
• Differences from the reference implementation
• Project layout
• Testing
• 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:

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

---

Installation

pip install shapneurgc

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

pip install "shapneurgc[all]"

From source, for development:

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:

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):

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

---

Sixty-second quickstart

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

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

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

# ---- 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.

---

Worked examples

Eight self-contained recipes. Each runs as written.

1. Minimal fit on your own data

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

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

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

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:

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

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

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

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

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	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	Every public signature, argument, default and return type.
THEORY.md	Each paper equation mapped to the line of code that implements it, with the numerical verification.
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.

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.

---

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

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:

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	fit → attribute → select → plot, in 30 lines	~1 min
02_simulation_var.py	VAR(3) benchmark, 4 penalties × 2 sample sizes; knockoff vs stability calibration	~5 min
03_monetary_policy.py	full applied study on real FRED data	~5 min
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.

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

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:

@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:

@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 · merwanroudane920@gmail.com

Released under the MIT 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.