bartardl

Non-parametric ARDL modelling with Bayesian Additive Regression Trees (BART).

A faithful, self-contained Python implementation of

Mahdavi, P., Ehsani, M. A., Ahelegbey, D. F. and Mohammadpour, M. (2024). Measuring Causal Effect with ARDL-BART: A Macroeconomic Application. Applied Mathematics, 15, 292–312. doi:10.4236/am.2024.154018

bartardl fits an autoregressive distributed lag (ARDL) model whose conditional mean is estimated non-parametrically with BART instead of OLS, so it can capture the non-linear and interaction effects the paper argues dominate macroeconomic data. It benchmarks BART against LASSO, Elastic Net and a Bayesian-Network (spike-and-slab) selector, and ships publication-quality figures and tables.

The BART sampler is written from scratch in NumPy — no external BART engine — so the package installs with only numpy, scipy, pandas, scikit-learn and matplotlib.

Highlights

• 🌲 Faithful BART — sum-of-trees with the exact priors of Chipman, George and McCulloch (2010): tree prior α(1+d)^(−β), conjugate-normal leaves, and the (ν, q)-calibrated inverse-gamma error prior. Paper defaults (m = 200, k = 2, ν = 3, q = 0.9, α = 0.95, β = 2) are the library defaults.
• 🔁 ARDL layer — one equation per target, every variable entered at lags 1…p, with the paper's Table 3 stationarity transforms built in.
• 🧮 Competitors as drop-ins — Lasso, ElasticNet, BayesianNetwork share the estimator protocol, so any of them can be the ARDL conditional mean.
• 📊 Top-journal graphics — RMSE box plots, variable-inclusion bars with 95% whiskers, actual-vs-fitted overlays, LaTeX/HTML/console tables, and the MATLAB Parula palette (parula_colors, turbo_colors, …).
• 🧪 Reproduces the paper — the Friedman simulations and the eight-variable US-macro application, with an offline synthetic panel so every example runs without a network.

Installation

git clone https://github.com/merwanroudane/bartardl.git
cd bartardl
pip install -e .

# optional: live FRED download for the real macro data
pip install -e ".[fred]"

Requires Python ≥ 3.9.

Quick start

import bartardl as ba

data  = ba.load_macro()                     # stationary US macro panel (8 vars)
model = ba.ARDLBART(                         # ARDL(4) with a BART conditional mean
    n_lags=4,
    estimator=ba.BART,
    estimator_kwargs=dict(n_trees=100, n_burn=500, n_draws=500, random_state=0),
)
res = model.fit(data, target="GDP")

print("In-sample RMSE:", ba.rmse(res.y_train, res.fitted))
print(res.top_features(5))                  # 5 most influential lagged regressors

In-sample RMSE: 0.0019

  feature variable  lag  importance
  PC.lag1       PC    1    0.0892
  IP.lag1       IP    1    0.0760
 GDP.lag1      GDP    1    0.0601
   U.lag3        U    3    0.0402
 INF.lag3      INF    3    0.0401

(this is exactly python examples/quickstart.py)

The model

For a target y the ARDL-BART equation is

y_t = f( y_{t-1}, …, y_{t-p},  x_{1,t-1}, …, x_{k,t-p} ) + e_t ,   e_t ~ N(0, σ²)

and f is a sum of m regression trees (BART):

f(x) = Σ_{l=1}^{m} g(x ; T_l , M_l)

The priors (paper §2) regularise every tree to be a weak learner:

The ensemble is drawn by Bayesian back-fitting MCMC: each sweep proposes a single grow or prune move per tree (Metropolis–Hastings), then re-samples the leaf values and σ². Variable importance is the inclusion proportion — the share of all splitting rules that use each predictor.

Faithfulness notes. The paper's "ARDL" is the single-equation, distributed-lag regression sense (no bounds test / no error-correction term), and its "causal effect" framing is operationalised here as forecasting + variable importance, exactly as in the paper's own experiments. See docs/guide.md for a point-by-point reconciliation with the article.

Reproducing the paper

1 · Simulations (Tables 1–2, Figures 1–6)

python examples/01_simulation_linear.py and python examples/02_simulation_nonlinear.py.

Relative RMSE (RMSE ÷ noise σ; a perfect model scores 1.0), p = 10:

Non-linear DGP — BART wins (paper Table 2)

Linear DGP — the linear selectors win (paper Table 1)

2 · Empirical application (Table 4, Figure 7)

python examples/03_macro_forecasting.py — eight US macro variables, ARDL(4), hold-out forecast RMSE with per-target ranks:

BART ranks first for 6 of 8 targets, with the linear selectors taking the two (near-)linear series — the same qualitative pattern as the paper's Table 4.

The bundled panel is a reproducible synthetic data set (built so its Table-3 transforms carry genuine, economically-motivated non-linearity) so the example runs offline. For the real series use ba.load_macro(source="fred") (needs internet + pandas-datareader).

Syntax guide

A complete reference for every public entry point: parameters, defaults, return types and a runnable snippet for each. Import once as import bartardl as ba.

Contents: The estimator protocol · Data preparation · ARDLBART · ARDLResult · BART · Competitors · Simulation · Metrics & tables · Datasets · Visualisation · Colours · End-to-end recipes

0 · The estimator protocol

Every conditional-mean model in the package — BART, Lasso, ElasticNet, BayesianNetwork — implements the same three-part interface, which is why any of them can be dropped into ARDLBART(estimator=…):

est.fit(X, y)        # X: (n, p) array, y: (n,) array          -> self
est.predict(Xnew)    # Xnew: (m, p) array                       -> (m,) array
est.importance_      # (p,) non-negative weights; BART exposes  inclusion_proportion_

To plug in your own model, satisfy those three and pass the class.

1 · Data preparation

transform_series(x, code) -> np.ndarray

Apply one FRED-MD stationarity transformation (paper Table 3). Leading values lost to differencing come back as NaN.

transform_frame(df, codes) -> pd.DataFrame

Transform each column by its code (codes is a {column: code} dict) and drop the rows lost to differencing.

make_lag_matrix(data, target, n_lags, include_target_lags=True) -> (X, y)

Build one ARDL equation's design matrix.

Returns X (columns named "<var>.lag<k>") and y (aligned target). With 8 variables and n_lags=4 → 32 regressors.

data = ba.load_macro()
X, y = ba.make_lag_matrix(data, target="GDP", n_lags=4)
X.shape            # (251, 32)
list(X.columns)[:3]   # ['GDP.lag1', 'GDP.lag2', 'GDP.lag3']

2 · ARDLBART

ARDLBART(n_lags=4, estimator=BART, include_target_lags=True, estimator_kwargs=None)

Single-equation ARDL with a pluggable conditional mean.

Methods

model = ba.ARDLBART(
    n_lags=4,
    estimator=ba.BART,
    estimator_kwargs=dict(n_trees=200, n_draws=1000, random_state=0),
)
res  = model.fit(data, target="GDP")            # -> ARDLResult
pred = model.predict(data, target="GDP")        # -> pandas Series (aligned to dates)

3 · ARDLResult

The object returned by ARDLBART.fit.

res.top_features(5)
res.importance_frame().query("variable == 'M2'")   # all M2 lags

4 · BART

BART(config=None, **overrides)          # or BART(n_trees=200, n_draws=1000, …)

A self-contained Bayesian Additive Regression Trees regressor. Pass options either as a BARTConfig or directly as keyword arguments.

Methods & fitted attributes

m = ba.BART(n_trees=200, k=2.0, nu=3.0, q=0.90,
            alpha=0.95, beta=2.0, random_state=0).fit(X, y)
m.predict(Xnew)
lo, hi = np.percentile(m.yhat_train_draws_, [2.5, 97.5], axis=0)   # 95% band

Match the paper exactly: ba.BART(n_trees=200, k=2, nu=3, q=0.9, alpha=0.95, beta=2, n_draws=20000, n_burn=250).

5 · Competitors

All three share the estimator protocol plus a .coef_.

Lasso(cv=5, random_state=None, max_iter=20000)

Cross-validated L1 regression (Tibshirani, 1996). .importance_ = |coef| normalised.

ElasticNet(cv=5, l1_ratio=None, random_state=None, max_iter=20000)

Cross-validated L1/L2 (Zou & Hastie, 2005). l1_ratio (the paper's α) is a grid searched by CV; default leans toward the LASSO end.

BayesianNetwork(n_iter=2000, burn=1000, pi=0.5, tau=0.1, c=100.0, random_state=None)

Spike-and-slab (SSVS) Gibbs sampler. A non-zero coefficient is read as a present network edge (Ahelegbey–Billio–Casarin, 2016).

Extra attribute: .inclusion_prob_ (posterior mean of the edge indicators).

for est, kw in [(ba.Lasso, {}), (ba.ElasticNet, {}),
                (ba.BayesianNetwork, dict(n_iter=2000, burn=1000))]:
    r = ba.ARDLBART(n_lags=4, estimator=est, estimator_kwargs=kw).fit(data, "GDP")
    print(est.__name__, ba.rmse(r.y_train, r.fitted))

6 · Simulation

friedman(n=100, p=10, kind="nonlinear", sigma=1.0, random_state=None) -> SimData

Draw a Friedman (1991) data set; only x1…x5 are relevant.

SimData fields: .X (n,p), .y (n,), .f (noise-free signal), .relevant (array([0..4])), .kind.

d = ba.friedman(n=200, p=100, kind="nonlinear", sigma=1.0, random_state=0)
m = ba.BART(random_state=0).fit(d.X, d.y)

7 · Metrics & tables

• rmse_draws is {method: array_of_rmse_per_replication}.
• rmse_by_target is {target: {method: rmse}}.

tab = ba.simulation_table({"BART": bart_draws, "LASSO": lasso_draws})
rank = ba.ranking_table({"GDP": {"BART": 0.0032, "LASSO": 0.0057}, ...})
ba.best_method_counts(rank)["BART"]

8 · Datasets

source ∈ {"auto", "bundled", "synthetic", "fred"}. transform=True applies the Table-3 codes and drops differencing NaNs.

data = ba.load_macro()                    # offline, transformed, ready to model
levels = ba.load_macro(transform=False)   # raw levels
real = ba.load_macro(source="fred")       # live series (pip install bartardl[fred])
ba.transform_table()

9 · Visualisation

from bartardl import viz — every function returns a Matplotlib Axes.

from bartardl import viz
ax = viz.rmse_boxplot(rmse_draws, title="Non-linear DGP (p = 10)")
viz.inclusion_plot(res.importance, res.feature_names, top=12,
                   title="BART inclusion proportions for GDP")
print(viz.journal_table(rank, fmt="latex", caption="Table 4"))

10 · Colours

from bartardl import colors

from bartardl.colors import parula_colors, get_cmap
parula_colors(8)                 # ['#352a87', '#0268e1', …]
plt.imshow(Z, cmap=get_cmap("Parula"))

11 · End-to-end recipes

Compare all four estimators on one target and rank them

import bartardl as ba
data = ba.load_macro(); H = 40; train = data.iloc[:-H]
methods = {
    "BART":         (ba.BART,            dict(n_trees=200, n_draws=1000, random_state=0)),
    "LASSO":        (ba.Lasso,           {}),
    "Elastic-Net":  (ba.ElasticNet,      {}),
    "Bayesian-Net": (ba.BayesianNetwork, dict(n_iter=2000, burn=1000, random_state=0)),
}
rmse_by_target = {}
for tgt in data.columns:
    rmse_by_target[tgt] = {}
    for name, (est, kw) in methods.items():
        m = ba.ARDLBART(n_lags=4, estimator=est, estimator_kwargs=kw).fit(train, tgt)
        pred = m.predict(data, tgt).iloc[-H:]
        rmse_by_target[tgt][name] = ba.rmse(data[tgt].loc[pred.index], pred)
print(ba.ranking_table(rmse_by_target))

Run the non-linear Monte-Carlo and summarise (Table 2)

draws = {"BART": [], "LASSO": []}
for rep in range(100):
    tr = ba.friedman(n=100, p=10, kind="nonlinear", random_state=rep)
    te = ba.friedman(n=100, p=10, kind="nonlinear", random_state=1000 + rep)
    for name, est in [("BART", ba.BART(random_state=0)), ("LASSO", ba.Lasso())]:
        est.fit(tr.X, tr.y)
        draws[name].append(ba.relative_rmse(te.y, est.predict(te.X), 1.0))
print(ba.simulation_table(draws))

Documentation

• docs/guide.md — a detailed, from-the-ground-up guide to how the code is written: the BART sampler internals (grow/prune acceptance ratios, back-fitting, σ calibration), the ARDL design, how the competitors and the SSVS sampler work, how to extend the package, and a full reconciliation with the paper.
• examples/ — runnable scripts for the quick start, both simulations, and the macro application.
• docs/make_figures.py — regenerates every figure and table in this README.

Testing

pip install -e ".[dev]"
pytest -q

Citing

If you use this software, please cite the paper and the package:

@article{mahdavi2024ardlbart,
  author  = {Mahdavi, Pegah and Ehsani, Mohammad Ali and
             Ahelegbey, Daniel Felix and Mohammadpour, Mehrnaz},
  title   = {Measuring Causal Effect with {ARDL-BART}: A Macroeconomic Application},
  journal = {Applied Mathematics},
  year    = {2024},
  volume  = {15},
  pages   = {292--312},
  doi     = {10.4236/am.2024.154018}
}

@software{roudane_bartardl,
  author = {Merwan Roudane},
  title  = {bartardl: Non-parametric ARDL modelling with BART},
  url    = {https://github.com/merwanroudane/bartardl}
}

License

MIT © Merwan Roudane. See LICENSE.