5.6 计量分析6：异方差的诊断和处理

于由管理员发布

例子

以课本表11-5（第396页）中的数据为样本，建立如下回归模型：

img

代码部分

# lec06
# 计量实验4-代码部分 异方差的诊断与处理

# 导入程序包
import pandas as pd
import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.stats.stattools import jarque_bera
import ssl

# 禁用SSL证书校验
ssl._create_default_https_context = ssl._create_unverified_context

# 导入(在线)数据
df = pd.read_excel(r"https://cdn.seit2019.xyz/data/econometrics/table11.5.xlsx", index_col="Industry")

# 展示数据
print("数据展示如下：")
print(df, "\n")

# 简要统计描述
print("简要统计描述：")
print(df.describe(), "\n")

# OLS回归
# 定义变量
y = df["RD"]
x = df["Sales"]
# 添加常数项
X = sm.add_constant(x)

# OLS回归
OLS_model = sm.OLS(y, X).fit()
print("OLS回归结果：")
print(OLS_model.summary(), "\n")

# 生成残差序列
df["e"] = OLS_model.resid
df["e^2"] = df["e"]**2
df["Log(e^2)"] = df["e^2"].apply(np.log)
df["Abs(e)"] = df["e"].abs()

# 残差正态性
df_jb = pd.DataFrame()
df_jb["Stats"] = ["Jarque_Bera", "P_value", "Skew", "Kurt"]
df_jb["Value"] = jarque_bera(OLS_model.resid)
df_jb["Value"] = df_jb["Value"].round(4)
print("残差正态性：")
print(df_jb.to_string(index=False), "\n")
# 绘图：设置绘图格式,支持汉字
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘图：残差直方图
plt.hist(OLS_model.resid, bins=7)
plt.title("残差直方图", fontsize="xx-large", fontweight="bold")
plt.show()

# 绘图：残差散点图
# 绘图：设置绘图格式,支持汉字
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘图：定义变量
x = df.index
y = df["e"]
# 绘图：残差散点图
plt.scatter(x, y, c="blue")
plt.title("残差散点图", fontsize="xx-large", fontweight="bold")
plt.xlabel("Industry", fontsize="xx-large", fontweight="bold")
plt.ylabel("e", fontsize="xx-large", fontweight="bold")
plt.xticks([])
plt.show()

# 绘图：残差与解释变量的散点图
# 绘图：定义变量
x = df["Sales"]
y = df["e"]
# 绘图：残差散点图
plt.scatter(x, y, c="blue")
plt.title("残差与解释变量的散点图", fontsize="xx-large", fontweight="bold")
plt.xlabel("Sales", fontsize="xx-large", fontweight="bold")
plt.ylabel("e", fontsize="xx-large", fontweight="bold")
plt.show()

# 设置检验结果汇总表
df_test = pd.DataFrame(columns=["Test_Name", "F_stat", "P_value", "Sig"])
df_test["Test_Name"] = ["Park test", "Glejser Test", "White Test", "B-P-G Test"]

# 异方差诊断

# 帕克检验(Park test)
# 定义变量
y1 = df["Log(e^2)"]
df["Log(Sales)"] = df["Sales"].apply(np.log)
x1 = df["Log(Sales)"]
# 辅助回归
X1 = sm.add_constant(x1)
OLS_aux1 = sm.OLS(y1, X1).fit()
print("辅助回归1结果：")
print(OLS_aux1.summary(), "\n")
df_test.iloc[0, 1] = OLS_aux1.fvalue.round(4)
df_test.iloc[0, 2] = OLS_aux1.f_pvalue.round(4)
# F统计量显著性
if OLS_aux1.f_pvalue < 0.05:
    df_test.iloc[0, 3] = "显著"
    print("结论：帕克检验(Park test)的F统计量显著，OLS回归模型中存在异方差。")
else:
    df_test.iloc[0, 3] = "不显著"
    print("结论：帕克检验(Park test)的F统计量不显著，OLS回归模型中不存在异方差。")

# 格莱泽检验(Glejser Test)
# 定义变量
y2 = df["Abs(e)"]
x2 = df["Sales"]
# 添加常数项
X2 = sm.add_constant(x2)
# 辅助回归
OLS_aux2 = sm.OLS(y2, X2).fit()
print("辅助回归2结果：")
print(OLS_aux2.summary(), "\n")
df_test.iloc[1, 1] = OLS_aux2.fvalue.round(4)
df_test.iloc[1, 2] = OLS_aux2.f_pvalue.round(4)
# F统计量显著性
if OLS_aux2.f_pvalue < 0.05:
    df_test.iloc[1, 3] = "显著"
    print("结论：格莱泽检验(Glejser Test)的F统计量显著，OLS回归模型中存在异方差。")
else:
    df_test.iloc[1, 3] = "不显著"
    print("结论：格莱泽检验(Glejser Test)的F统计量不显著，OLS回归模型中不存在异方差。")

# 怀特检验(White Test)
# 定义变量
y3 = df["e^2"]
df["Sales^2"] = df["Sales"]**2
x3 = df["Sales^2"]
# 添加常数项
X3 = sm.add_constant(x3)
# 辅助回归
OLS_aux3 = sm.OLS(y3, X3).fit()
print("辅助回归3结果：")
print(OLS_aux3.summary(), "\n")
df_test.iloc[2, 1] = OLS_aux3.fvalue.round(4)
df_test.iloc[2, 2] = OLS_aux3.f_pvalue.round(4)
# F统计量显著性
if OLS_aux3.f_pvalue < 0.05:
    df_test.iloc[2, 3] = "显著"
    print("结论：怀特检验(White Test)的F统计量显著，OLS回归模型中存在异方差。")
else:
    df_test.iloc[2, 3] = "不显著"
    print("结论：怀特检验(White Test)的F统计量不显著，OLS回归模型中不存在异方差。")

# B-P-G检验（Breusch-Pagan-Godfrey Test）
# 定义变量
y4 = df["e^2"]
x4 = df["Sales"]
# 添加常数项
X4 = sm.add_constant(x4)
# 辅助回归
OLS_aux4 = sm.OLS(y4, X4).fit()
print("辅助回归4结果：")
print(OLS_aux4.summary(), "\n")
df_test.iloc[3, 1] = OLS_aux4.fvalue.round(4)
df_test.iloc[3, 2] = OLS_aux4.f_pvalue.round(4)
# F统计量显著性
if OLS_aux4.f_pvalue < 0.05:
    df_test.iloc[3, 3] = "显著"
    print("结论：B-P-G检验（Breusch-Pagan-Godfrey Test）的F统计量显著，OLS回归模型中存在异方差。")
else:
    df_test.iloc[3, 3] = "不显著"
    print("结论：B-P-G检验（Breusch-Pagan-Godfrey Test）的F统计量不显著，OLS回归模型中不存在异方差。")
print("\n")

# 异方差检验结果汇总
print("异方差检验结果汇总：")
# print(df_test, "\n")
print(df_test.to_string(index=False), "\n")
# 检验结果综合判断
if df_test["P_value"].min() < 0.05:
    print("最终结论：OLS模型存在异方差。")
    num = df_test["P_value"].values.argmin()
    if df_test.iloc[num, 0] == "B-P-G Test":
        print("异方差形式：误差正比于X.", "\n")
    else:
        print("异方差形式：误差正比于X平方.", "\n")
else:
    print("最终结论：OLS模型不存在异方差。", "\n")

# 异方差处理-WLS估计
# 异方差检验：B-P-G test；异方差形式：误差正比于x；误差处理：WLS,w=sqrt(sales)
# 定义变量
y = df["RD"]
x = df["Sales"]
# 定义权重
# w = df["Sales"]       # White test:误差方差正比于X平方
w = df["Sales"]**0.5    # B-P-G test:误差方差正比于X
# 添加常数项
X = sm.add_constant(x)
# WLS回归
WLS_model = sm.WLS(y, X, weights=1.0/(w**2)).fit()
print("WLS回归结果：")
print(WLS_model.summary())

# WLS估计结果的异方差检验
# 对WLS估计结果的异方差检验需要留意辅助回归的解释变量形式。Eviews软件中默认采用原始OLS回归的解释变量，这是错误的。
# 生成残差
df["e_wls"] = WLS_model.resid

# B-P-G检验（Breusch-Pagan-Godfrey Test）
# 定义变量
df["e^2_wls"] = df["e_wls"]**2
y = df["e^2_wls"]
df["Sales_c"] = df["Sales"]**(-0.5)
df["Sales_x"] = df["Sales"]**0.5
x = df[["Sales_c", "Sales_x"]]
# 添加常数项
X = sm.add_constant(x)
# 辅助回归
OLS_aux = sm.OLS(y, X).fit()
print("辅助回归结果：")
print(OLS_aux.summary(), "\n")
# F统计量显著性
if OLS_aux.f_pvalue < 0.05:
    print("结论：对WLS估计结果的B-P-G检验表明，F统计量显著，WLS回归模型中存在异方差。")
else:
    print("结论：对WLS估计结果的B-P-G检验表明，F统计量不显著，WLS回归模型中不存在异方差。")
print("\n")

# white检验（White Test）
# 定义变量
df["e^2_wls"] = df["e_wls"]**2
y = df["e^2_wls"]
df["Sales_c"] = df["Sales"]**(-1)
df["Sales_x"] = df["Sales"]
x = df[["Sales_c", "Sales_x"]]
# 添加常数项
X = sm.add_constant(x)
# 辅助回归
OLS_aux = sm.OLS(y, X).fit()
print("辅助回归结果：")
print(OLS_aux.summary(), "\n")
# F统计量显著性
if OLS_aux.f_pvalue < 0.05:
    print("结论：对WLS估计结果的White检验表明，F统计量显著，WLS回归模型中存在异方差。")
else:
    print("结论：对WLS估计结果的White检验表明，F统计量不显著，WLS回归模型中不存在异方差。")
print("\n")

# 对数模型 vs 线性模型
# 定义变量
df["Log(RD)"] = df["RD"].apply(np.log)
df["Log(Sales)"] = df["Sales"].apply(np.log)
y0 = df["RD"]
y1 = df["Log(RD)"]
x0 = df["Sales"]
x1 = df["Log(Sales)"]
# 添加常数项
X0 = sm.add_constant(x0)
X1 = sm.add_constant(x1)
# OLS回归
OLS_linear = sm.OLS(y0, X0).fit()
OLS_log = sm.OLS(y1, X1).fit()
# 显示回归结果
print("线性OLS回归结果：")
print(OLS_linear.summary(), "\n")
print("对数线性OLS回归结果：")
print(OLS_log.summary(), "\n")

# 异方差检验
# 生成残差
df["e_linear"] = OLS_linear.resid
df["e_log"] = OLS_log.resid
df["e^2_linear"] = df["e_linear"]**2
df["e^2_log"] = df["e_log"]**2
# 定义变量
y0 = df["e^2_linear"]
y1 = df["e^2_log"]
x0 = df["Sales"]
x1 = df["Log(Sales)"]
# 添加常数项
X0 = sm.add_constant(x0)
X1 = sm.add_constant(x1)
# OLS辅助回归
OLS_linear_aux = sm.OLS(y0, X0).fit()
OLS_log_aux = sm.OLS(y1, X1).fit()
# 显示回归结果
print("线性OLS-辅助回归结果：")
print(OLS_linear_aux.summary(), "\n")
print("对数线性OLS-辅助回归结果：")
print(OLS_log_aux.summary(), "\n")
# F统计量显著性

# Linear
if OLS_linear_aux.f_pvalue < 0.05:
    print("结论：B-P-G检验表明，F统计量显著，线性OLS回归模型中存在异方差。")
else:
    print("结论：B-P-G检验表明，F统计量不显著，线性OLS回归模型中不存在异方差。")
print("\n")

# Log
if OLS_log_aux.f_pvalue < 0.05:
    print("结论：B-P-G检验表明，F统计量显著，对数OLS回归模型中存在异方差。")
else:
    print("结论：B-P-G检验表明，F统计量不显著，对数OLS回归模型中不存在异方差。")
print("\n")

# Linear vs Log
if OLS_linear_aux.f_pvalue < OLS_log_aux.f_pvalue:
    print("结论：相比线性OLS回归，对数OLS回归能够减缓异方差的影响。")
else:
    print("结论：相比线性OLS回归，对数OLS回归不能减缓异方差的影响。")
print("\n")
分类： 教程笔记

管理员

管理员

相关文章

教程 笔记

5.10 计量分析10：文本分析

例子 对古典文学名著《红楼梦》进行文本分析，提取其中频数大于500的词 阅读更多…

教程 笔记

5.9 计量分析9：定性响应模型

例子：Probit模型与Logit模型的估计 代码部分 # lec09 阅读更多…

教程 笔记

5.8 计量分析8：虚拟变量模型

例子 以课本表9.1的数据为样本，建立虚拟变量模型： 代码部分 # l 阅读更多…

SEIT资源站