bayespecon Replication of panelg_manual and sdemo_programs¶

import numpy as np
import pandas as pd
from libpysal.graph import Graph
from scipy.sparse import csr_matrix

from bayespecon import SDEM, SDM, SDEMPanelFE, dgp

def make_knn_w(xcoord: np.ndarray, ycoord: np.ndarray, k: int) -> np.ndarray:
    """Build row-standardized k-NN weights (no self-neighbors)."""
    coords = np.column_stack([xcoord, ycoord])
    n = coords.shape[0]
    d = np.sqrt(((coords[:, None, :] - coords[None, :, :]) ** 2).sum(axis=2))
    np.fill_diagonal(d, np.inf)
    nn = np.argpartition(d, kth=k - 1, axis=1)[:, :k]

    W = np.zeros((n, n), dtype=float)
    for i in range(n):
        W[i, nn[i]] = 1.0

    rs = W.sum(axis=1, keepdims=True)
    rs[rs == 0.0] = 1.0
    return W / rs


def to_graph(W: np.ndarray) -> Graph:
    return Graph.from_sparse(csr_matrix(W)).transform("r")


def summarize_recovery(df: pd.DataFrame, abs_tol: float = 0.25) -> pd.DataFrame:
    out = df.copy()
    out["abs_error"] = (out["estimate"] - out["truth"]).abs()
    out["relative_error_pct"] = (
        100.0 * out["abs_error"] / np.maximum(np.abs(out["truth"]), 1e-12)
    )
    out["within_tol"] = out["abs_error"] <= abs_tol
    return out

rng = np.random.default_rng(10203040)
n, t, k = 100, 20, 2
lam_true = 0.7
beta_true = np.ones(k)
theta_true = -np.ones(k)
sige = 0.1

x = rng.standard_normal((n * t, k))
Wn = make_knn_w(rng.random(n), rng.random(n), k=5)
Wbig = np.kron(np.eye(t), Wn)

SFE = np.kron(np.ones(t), np.arange(1, n + 1) / n)
TFE = np.kron(np.arange(1, t + 1) / t, np.ones(n))
u = np.linalg.solve(
    np.eye(n * t) - lam_true * Wbig, rng.standard_normal(n * t) * np.sqrt(sige)
)
y = x @ beta_true + (Wbig @ x) @ theta_true + SFE + TFE + u

m_panel = SDEMPanelFE(
    y=y,
    X=pd.DataFrame(x, columns=["x1", "x2"]),
    W=to_graph(Wn),
    N=n,
    T=t,
    model=3,
    priors={"lam_lower": -0.95, "lam_upper": 0.95},
)
idata_panel = m_panel.fit(
    draws=160, tune=160, chains=1, target_accept=0.92, progressbar=False, random_seed=11
)
display(m_panel.summary(round_to=4))

beta_hat = idata_panel.posterior["beta"].mean(("chain", "draw")).to_numpy()
lam_hat = float(idata_panel.posterior["lam"].mean(("chain", "draw")).to_numpy())
faith_panel = pd.DataFrame(
    {
        "term": ["lam", "x1", "x2", "W*x1", "W*x2"],
        "truth": [lam_true, 1.0, 1.0, -1.0, -1.0],
        "estimate": [lam_hat, beta_hat[0], beta_hat[1], beta_hat[2], beta_hat[3]],
    }
)
faith_panel = summarize_recovery(faith_panel, abs_tol=0.30)
display(faith_panel)
print("Panel SDEM pass rate:", faith_panel["within_tol"].mean())

rng = np.random.default_rng(30203040)
n, k = 1000, 2
rho_true = 0.5
intercept_true = 2.0
beta_true = np.ones(k)
theta_true = -0.5 * np.ones(k)

Wn = make_knn_w(rng.random(n), rng.random(n), k=6)
sdm_data = dgp.simulate_sdm(
    W=Wn,
    rho=rho_true,
    beta1=np.array([intercept_true, *beta_true]),
    beta2=theta_true,
    sigma=np.sqrt(0.5),
    rng=rng,
)
y = sdm_data["y"]
X_arr = sdm_data["X"]

X = pd.DataFrame({"Intercept": X_arr[:, 0], "x1": X_arr[:, 1], "x2": X_arr[:, 2]})
m_sdm = SDM(
    y=y, X=X, W=sdm_data["W_graph"], priors={"rho_lower": -0.95, "rho_upper": 0.95}
)
idata_sdm = m_sdm.fit(
    draws=160, tune=160, chains=1, target_accept=0.92, progressbar=False, random_seed=22
)
display(m_sdm.summary(round_to=4))

post = idata_sdm.posterior
rho_hat = float(post["rho"].mean(("chain", "draw")).to_numpy())
beta_hat = post["beta"].mean(("chain", "draw")).to_numpy()

faith_sdm = pd.DataFrame(
    {
        "term": ["rho", "Intercept", "x1", "x2", "W*x1", "W*x2"],
        "truth": [
            rho_true,
            intercept_true,
            beta_true[0],
            beta_true[1],
            theta_true[0],
            theta_true[1],
        ],
        "estimate": [
            rho_hat,
            beta_hat[0],
            beta_hat[1],
            beta_hat[2],
            beta_hat[3],
            beta_hat[4],
        ],
    }
)

faith_sdm = summarize_recovery(faith_sdm, abs_tol=0.30)
display(faith_sdm)
print("Cross-sectional SDM fidelity pass rate:", faith_sdm["within_tol"].mean())

rng = np.random.default_rng(221010)
n, k = 3000, 4
lam_true = 0.8
intercept_true = 1.0
beta_true = np.ones(k)
theta_true = 0.5 * np.ones(k)

Wn = make_knn_w(rng.standard_normal(n), rng.standard_normal(n), k=5)
sdem_data = dgp.simulate_sdem(
    W=Wn,
    lam=lam_true,
    beta1=np.array([intercept_true, *beta_true]),
    beta2=theta_true,
    sigma=1.0,
    rng=rng,
)
y = sdem_data["y"]
X_arr = sdem_data["X"]

X = pd.DataFrame({"Intercept": X_arr[:, 0]})
for j in range(k):
    X[f"x{j + 1}"] = X_arr[:, j + 1]

m_sdem = SDEM(
    y=y, X=X, W=sdem_data["W_graph"], priors={"lam_lower": -0.95, "lam_upper": 0.95}
)
idata_sdem = m_sdem.fit(
    draws=160, tune=160, chains=1, target_accept=0.92, progressbar=False, random_seed=33
)
display(m_sdem.summary(round_to=4))

post = idata_sdem.posterior
lam_hat = float(post["lam"].mean(("chain", "draw")).to_numpy())
beta_hat = post["beta"].mean(("chain", "draw")).to_numpy()

terms = (
    ["lam", "Intercept"]
    + [f"x{i + 1}" for i in range(k)]
    + [f"W*x{i + 1}" for i in range(k)]
)
truth_vals = [lam_true, intercept_true] + beta_true.tolist() + theta_true.tolist()
est_vals = (
    [lam_hat, beta_hat[0]]
    + beta_hat[1 : 1 + k].tolist()
    + beta_hat[1 + k : 1 + 2 * k].tolist()
)

faith_sdem = pd.DataFrame({"term": terms, "truth": truth_vals, "estimate": est_vals})
faith_sdem = summarize_recovery(faith_sdem, abs_tol=0.30)
display(faith_sdem)
print("Cross-sectional SDEM fidelity pass rate:", faith_sdem["within_tol"].mean())

from bayespecon import SpatialProbit

rng = np.random.default_rng(707)

# semip-like setup: m regions with mobs observations each
m = 48
mobs = np.full(m, 60, dtype=int)
n = int(mobs.sum())
k = 3

# Region-level kNN weights and spatial random effects
Wm = make_knn_w(rng.random(m), rng.random(m), k=6)
rho_true = 0.45
sigma_a_true = 1.2

u = rng.standard_normal(m) * sigma_a_true
a = np.linalg.solve(np.eye(m) - rho_true * Wm, u)

region_ids = np.repeat(np.arange(m), mobs)
X = rng.standard_normal((n, k))
X = np.column_stack([np.ones(n), X])
feature_names = ["Intercept", "x1", "x2", "x3"]

beta_true = np.array([0.4, 0.8, -0.6, 0.5])
eta = X @ beta_true + a[region_ids]
y = (eta + rng.standard_normal(n) > 0.0).astype(float)

sp = SpatialProbit(
    y=y,
    X=pd.DataFrame(X, columns=feature_names),
    W=to_graph(Wm),
    region_ids=region_ids,
    priors={"rho_lower": -0.95, "rho_upper": 0.95, "beta_sigma": 5.0},
)
idata_sp = sp.fit(
    draws=220,
    tune=220,
    chains=1,
    target_accept=0.92,
    progressbar=False,
    random_seed=707,
)

display(sp.summary(var_names=["beta", "rho", "sigma_a"], round_to=4))

beta_hat = idata_sp.posterior["beta"].mean(("chain", "draw")).to_numpy()
rho_hat = float(idata_sp.posterior["rho"].mean(("chain", "draw")).to_numpy())

a_mean_hat = sp.random_effects_mean().to_numpy()
a_rmse = float(np.sqrt(np.mean((a_mean_hat - a) ** 2)))

faith_sp = pd.DataFrame(
    {
        "term": ["rho"] + feature_names,
        "truth": [rho_true] + beta_true.tolist(),
        "estimate": [rho_hat] + beta_hat.tolist(),
    }
)
faith_sp = summarize_recovery(faith_sp, abs_tol=0.30)

print(f"Regional-effects RMSE: {a_rmse:.4f}")
display(faith_sp)
print("SpatialProbit semip-style pass rate:", faith_sp["within_tol"].mean())