# 依赖: numpy, statsmodels
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

rng = np.random.default_rng(0)
X = rng.normal(size=(60, 3))
y = 1.5 + X @ np.array([2.0, -1.0, 0.5]) + rng.normal(scale=0.3, size=60)

Xd = sm.add_constant(X)
model = sm.OLS(y, Xd).fit()
print(model.summary().tables[1])
print('VIF =', [round(variance_inflation_factor(Xd, i), 2) for i in range(Xd.shape[1])])
print('残差 JB p =', round(sm.stats.stattools.jarque_bera(model.resid)[1], 4))