Neural-GC: Deep Learning Framework for Granger Causality

Practical implementations and demonstrations of the approach are provided in the following Jupyter notebooks:

• cmlp_lagged_var_demo.ipynb: Demonstrates the cMLP model with lagged variables.
• clstm_lorenz_demo.ipynb: Demonstrates the cLSTM model using Lorenz system data.
• crnn_lorenz_demo.ipynb: Demonstrates the cRNN model using Lorenz system data.

When using Neural-GC, consider the following best practices:

• Regularization Strength: Selecting the optimal strength is difficult. To evaluate multiple strengths efficiently, consider running parallel training jobs or implementing a warm start strategy.
• Initialization and Pretraining: Training a model without regularization (pretraining) and then applying ISTA may yield different results than training directly with ISTA due to the non-convex nature of the objective function. Experiment with both random initialization and pretraining to find the best results for your data.
• Debiasing: To train a model that uses the learned sparsity pattern without the bias introduced by the penalty, use the cMLPSparse, cLSTMSparse, or cRNNSparse classes.

To install Neural-GC, clone the repository and ensure you have the following dependencies installed:

• Python 3
• PyTorch (>= 0.4.0)
• numpy
• scipy

# Clone the repository manually
# Ensure dependencies are installed via pip or conda
pip install torch numpy scipy

To perform Granger causality discovery, initialize a cMLP model and train it using the train_model_ista function. This method uses Iterative Soft-Thresholding Algorithm (ISTA) to induce sparsity in the learned interactions.

Key parameters for train_model_ista:

• lam: Sparsity regularization strength.
• lam_ridge: Ridge regularization strength.
• lr: Learning rate.
• penalty: The type of penalty to use (e.g., 'H').
• max_iter: Maximum number of training iterations.

Note: The model should be moved to the desired device (e.g., .cuda()) before training.

from models.cmlp import cMLP, train_model_ista

# Initialize model
cmlp = cMLP(input_dim, lag=5, hidden=[100]).cuda(device=device)

# Train with ISTA
train_loss_list = train_model_ista(
    cmlp, X, lam=0.002, lam_ridge=1e-2, lr=5e-2, penalty='H', max_iter=50000,
    check_every=100
)

After discovering the causal structure with a regularized cMLP, you can build a 'debiased' model that only includes the discovered interactions. This is done using cMLPSparse and training it with train_unregularized to remove the bias introduced by the sparsity penalty.

Workflow:

1. Extract the sparsity pattern from the trained cMLP using .GC().bool().
2. Initialize cMLPSparse using that sparsity pattern.
3. Train using train_unregularized with a standard learning rate.

from models.cmlp import cMLPSparse, train_unregularized

# 1. Get sparsity pattern from trained cMLP
sparsity = cmlp.GC().bool()

# 2. Create debiased model
cmlp_sparse = cMLPSparse(X.shape[-1], sparsity, lag=5, hidden=[100]).cuda(device=device)

# 3. Train without regularization
train_loss_list = train_unregularized(
    cmlp_sparse, X, lr=1e-3, max_iter=10000, check_every=100, verbose=1
)

If you have already learned a sparsity pattern and wish to train a debiased model using that pattern, use the specific 'Sparse' variants of the model classes instead of the standard versions:

• cMLPSparse
• cLSTMSparse
• cRNNSparse

A trained cMLP (or a debiased cMLPSparse) can be used as a forecasting model by passing the input time series X to the model instance. The model will output predicted values for the subsequent time steps.

# Assuming cmlp_sparse is a trained debiased model
X_pred = cmlp_sparse(X)

Use simulate_lorenz_96 from the synthetic module to generate time series data and the ground truth Granger Causality (GC) matrix.

Parameters:

• p: Number of variables in the system.
• F: Forcing parameter.
• T: Number of time points.

The function returns a tuple containing the simulated data X_np (as a NumPy array) and the true Granger Causality matrix GC.

from synthetic import simulate_lorenz_96

# p=10 variables, F=10 forcing, T=1000 time points
X_np, GC = simulate_lorenz_96(p=10, F=10, T=1000)

Once a cMLP model is trained, you can extract the estimated Granger causality (GC) matrix using the .GC() method. This matrix represents the discovered causal relationships between variables.

To inspect the lag-specific causality, use .GC(ignore_lag=False, threshold=False). This returns the causal influence mapped across different lags.

# Get the general Granger causality matrix
GC_est = cmlp.GC().cpu().data.numpy()

# Get lag-specific Granger causality
GC_est_lag = cmlp.GC(ignore_lag=False, threshold=False)

Use train_model_ista to train the cRNN model on your time series data X. This function implements Iterative Soft-Thresholding Algorithm (ISTA) for sparse regularization.

Arguments:

• crnn: The cRNN model instance.
• X: The input time series tensor (expected shape [1, T, variables]).
• context: Context window size.
• lam: Sparsity parameter (L1 penalty).
• lam_ridge: Ridge regularization parameter.
• lr: Learning rate.
• max_iter: Maximum number of training iterations.
• check_every: Frequency of loss logging/checking.

from models.crnn import train_model_ista

train_loss_list = train_model_ista(
    crnn, X, context=10, lam=10.0, lam_ridge=1e-2, lr=1e-3, max_iter=20000,
    check_every=50
)