1 数据说明
“car_sales.xlsx”数据中包含了汽车的多项参数和销量信息，我们将根据该数据建立汽车参数对企业销量影响的多元线性回归模型。
2 初步回归
# 引入模块
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
# 读取外部数据
data = pd.read_excel('D:\\Desktop\\car_sales.xlsx', sheet_name='Sheet1', na_values='n/a')
data_balance = data.dropna(axis=0, how='any', thresh=None, subset=None, inplace=False)
y=data_balance.loc[:, 'lnsales']
X=data_balance.loc[:,'type':'mpg']
X=sm.add_constant(X)
ols_model = sm.OLS(y, X)
ols_results = ols_model.fit()
print(ols_results.summary()) # 输出模型的结果
初步回归结果显示，大部分变量不太显著，且Notes：的第2条显示模型中的多重共线问题比较严重。应当进行变量筛选！
3 自变量的筛选
3.1 方差膨胀因子法（VIF）
def calculate_vif(df):
vif = pd.DataFrame()
vif['index'] = df.columns
vif['VIF'] = [variance_inflation_factor(df.values,i) for i in range(df.shape[1])]
return vif
X_vif=X.iloc[:,1:]
vif = calculate_vif(X_vif)
while (vif['VIF'] > 10).any():
remove = vif.sort_values(by='VIF',ascending=False)['index'][:1].values[0]
X_vif.drop(remove,axis=1,inplace=True)
vif = calculate_vif(X_vif)
Selected_var = vif.iloc[:-1,:]['index'].values
print(Selected_var)
X_selected = X[Selected_var]
X_selected=sm.add_constant(X_selected)
ols_model_vif = sm.OLS(y, X_selected)
ols_results_vif = ols_model_vif.fit()
print(ols_results_vif.summary())
3.2 逐步回归方法
statsmodels没有可以直接调用的逐步回归方法，这里简单根据系数的显著性水平是否大于0.1来进行逐步回归。一般来说，在多元线性回归中，即使常数项也不显著，也应该加入模型中。代码如下：
X_step = X
pv_var = pd.DataFrame()
pv_var['index'] = X_step.columns
print(X_step.columns)
pv_var['pv'] = pd.Index(ols_results.pvalues)
while (pv_var['pv'] > 0.1).any():
remove = pv_var.sort_values(by='pv',ascending=False)['index'][:1].values[0]
X_step.drop(remove,axis=1,inplace=True)
ols_model_step = sm.OLS(y, X_step)
ols_results_step = ols_model_step.fit()
pv_var = pd.DataFrame()
pv_var['index'] = X_step.columns
pv_var['pv'] = pd.Index(ols_results_step.pvalues)
if 'const' not in X_step.columns:
X_step = sm.add_constant(X_step)
ols_model_step = sm.OLS(y, X_step)
ols_results_step = ols_model_step.fit()
print(ols_results_step.summary())
3.3 根据理论分析来筛选变量
在计量经济分析中，通过将自变量分为主要解释变量和控制变量，主要解释变量必须加入模型中，而控制变量只要不会对主要解释变量的系数产生影响，即使不显著也可以放入模型中，所以可以手动进行筛选。这里不再演示！
由于方差膨胀因子方法过于严格，所以我们选择逐步回归确定的自变量，接着进行分析！
4 绘制y的真实值与预测值图
prd_ols = ols_results_step.get_prediction()
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(range(0, len(y)), y, "o", label="data")
ax.plot(range(0, len(y)), ols_results_step.fittedvalues, "r--.", label="OLS")
ax.legend(loc="best")
plt.show()
5 残差正态性检验
5.1 Jarque-Bera test:
name = ["Jarque-Bera", "Chi^2 two-tail prob.", "Skew", "Kurtosis"]
test = sms.jarque_bera(ols_results_step.resid)
print(lzip(name, test))
运行结果如下，其中P值小于0.05，说明残差不服从正态分布
[('Jarque-Bera', 60.148856708878256), ('Chi^2 two-tail prob.', 8.686437901525318e-14), ('Skew', -0.9972366732963988), ('Kurtosis', 5.89141887088964)]
5.2 Omni test:
name = ["Chi^2", "Two-tail probability"]
test = sms.omni_normtest(ols_results_step.resid)
print(lzip(name, test))
运行结果如下，其中P值小于0.05，说明残差不服从正态分布
[('Chi^2', 28.90701709031857), ('Two-tail probability', 5.283491276420808e-07)]
6 加权最小二乘法
由于残差不服从正态分布，可能存在异方差能问题，所以这里采用加权最小二乘法进一步进行估计。
w = np.std(ols_results_step.resid)
mod_wls = sm.WLS(y, X_step, weights=1.0 / w)
res_wls = mod_wls.fit()
print(res_wls.summary())
加权最小二乘法的结果即为本案例最终选择的模型形式。price和wheelbas两个变量对汽车销量有显著的负面影响。
当然，由于type和fuel_cap变量不显著，所以也可以继续剔除这两个变量后再次进行加权最小二乘估计。同时也可以继续对残差的正态性进行检验，看结果是否有所好转！