5.3 计量分析3：一元线性回归模型

例子2.1

为了研究2015年中国各省农村居民人均消费支出（Cons）和可支配收入（DPI）之间的关系，我们以人均消费支出（Cons）为因变量，以人均可支配收入（DPI）为自变量,进行回归，估计如下一元线性回归模型：

数据来源：2016年《中国统计年鉴》 表6-31，表6-29

代码部分

# lec03
# 计量实验2-代码部分

# 导入库文件
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_white
from statsmodels.stats.diagnostic import het_breuschpagan
import ssl

# 导入本地excel数据（路径：当前目录/data/）
# df01 = pd.read_excel(r"./data/cons_rural_2020.xlsx")

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

# 导入在线数据
df01 = pd.read_excel(r"https://cdn.seit2019.xyz/data/econometrics/cons_rural_2020.xlsx")

# 展示数据
print("原始数据展示如下：")
print(df01, "\n")

# 简要描述性统计分析
print("原始数据的主要统计指标：")
print(df01.describe(), "\n")

# 绘图

# 设置绘图格式,支持汉字
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 绘制散点图

# 定义变量
x01 = df01["地区"]
y01 = df01["2020农村居民人均消费支出(元)"]
y02 = df01["2020农村居民人均可支配收入(元)"]

# 2020农村居民人均消费支出(元)散点图
plt.scatter(x01, y01, c="blue")
plt.xlabel("地区")
plt.xticks([])             # 关闭x轴刻度显示
plt.ylabel("农村居民人均消费支出(元)")
plt.title("2020农村居民人均消费支出散点图", fontsize="xx-large", fontweight="bold")
plt.show()

# 2020农村居民人均可支配收入(元)散点图
plt.scatter(x01, y02, c="green")
plt.xlabel("地区")
plt.xticks([])             # 关闭x轴刻度显示
plt.ylabel("农村居民人均可支配收入(元)")
plt.title("2020农村居民人均可支配收入散点图", fontsize="xx-large", fontweight="bold")
plt.show()

# 绘制XY散点图

# 定义变量
y1 = df01["2020农村居民人均消费支出(元)"]
x1 = df01["2020农村居民人均可支配收入(元)"]

# 消费支出与可支配收入散点图
plt.scatter(x1, y1, c="black")
plt.xlabel("农村居民人均可支配收入(元)")
plt.ylabel("农村居民人均消费支出(元)")
plt.title("消费支出与可支配收入散点图", fontsize="xx-large", fontweight="bold")
plt.show()

# 统计描述

# 重新定义数据:剔除第一列,修改列名
df02 = df01.iloc[:, 1:3]
df02 = df02.rename(columns={"2020农村居民人均消费支出(元)": "CONS", "2020农村居民人均可支配收入(元)": "DPI"})
print("回归分析的数据为：")
print(df02, "\n")

# 定义变量
y2 = df02["CONS"]
x2 = df02["DPI"]

# 展示数据的简要统计描述
print("数据的主要统计指标：")
print(df02.describe(), "\n")

# 展示统计指标
# 定义统计指标函数
def stat(x):
    return pd.Series([x.count(), x.min(), x.quantile(.25), x.median(), x.quantile(.75), x.max(), x.mean(), x.var(), x.std(), x.skew(), x.kurt(), x.max()-x.min()], index=['总数', '最小值', '下四分位数', '中位数', '上四分位数', '最大值', '平均数', '方差', '标准差', '偏度', '峰度', '极差'])

# 计算统计指标
s1 = pd.DataFrame(stat(y2))
s2 = pd.DataFrame(stat(x2))

# 合并显示统计指标
dfs = pd.DataFrame(columns=['CONS', 'DPI'])
dfs['CONS'] = s1
dfs['DPI'] = s2
print(dfs, "\n")
dfs.to_excel('./data/cons_rural_2020_stats.xlsx')

# 展示数据的协方差矩阵
print("数据的协方差矩阵：")
print(df02.cov(), "\n")

# 展示数据的相关系数矩阵
print("数据的相关系数矩阵：")
print(df02.corr(), "\n")

# 一元线性回归

# 添加常数项
X = sm.add_constant(x2)

# OLS回归
OLS_model = sm.OLS(y2, X).fit()

# 显示OLS回归结果
print(OLS_model.summary(), "\n")
dfc = pd.DataFrame(OLS_model.params).round(2)     # 获取回归参数,保留小数点后4位
dfc = dfc.astype("string")                        # 转换格式为字符串
beta0 = dfc.iloc[0, 0]                            # 定义beta0
beta1 = dfc.iloc[1, 0]                            # 定义beta1
OLS_equation = "Y=" + beta0 + "+" + beta1 + "*X"  # 定义回归方程，格式：字符串
print("估计的回归方程为：")
print(OLS_equation, "\n")

# 绘制趋势和散点图
y_fitted = OLS_model.predict()
plt.scatter(x2, y2, label="实际值")
plt.plot(x2, y_fitted, color="black", label="预测值")
plt.title("农村居民人均消费支出对人均可支配收入的回归", fontsize="xx-large", fontweight="bold")
plt.xlabel("农村居民人均可支配收入（元）")
plt.ylabel("农村居民人均消费支出（元）")
plt.text(27000, 17800, s=OLS_equation, fontsize="x-large")
plt.legend()
plt.show()

# 残差检验
model_resid = OLS_model.resid    # 生成残差
# 图示法
# 残差自身序列
plt.scatter(x01, model_resid)
plt.title("残差散点图", fontsize="xx-large", fontweight="bold")
plt.xlabel("地区")
plt.ylabel("残差")
plt.xticks([])             # 关闭x轴刻度显示
# plt.legend()
plt.show()

# 残差对X序列的散点图
plt.scatter(x2, model_resid)
plt.title("残差对X散点图", fontsize="xx-large", fontweight="bold")
plt.xlabel("农村居民人均可支配收入（元）")
plt.ylabel("残差")
# plt.xticks([])             # 关闭x轴刻度显示
# plt.legend()
plt.show()

# 残差对y_fitted序列的散点图
plt.scatter(y_fitted, model_resid)
plt.title("残差对y_fitted散点图", fontsize="xx-large", fontweight="bold")
plt.xlabel("农村居民人均消费支出拟合值")
plt.ylabel("残差")
# plt.xticks([])             # 关闭x轴刻度显示
# plt.legend()
plt.show()

# 异方差检验
# white test
white_output = het_white(model_resid, X)
white_result = pd.Series(white_output[0:4], index=['white_lm_statistic', 'white_lm_p_value', 'white_F_statistic', 'white_F_statistic_p_value'])
print("White异方差检验结果：")
print(white_result)
if white_result['white_lm_p_value'] < 0.05:
    print("White异方差检验表明：存在异方差", "\n")
else:
    print("White异方差检验表明：不存在异方差", "\n")

# Breuschpagan test
bp_output = het_breuschpagan(model_resid, X)
bp_result = pd.Series(bp_output[0:4], index=['bp_lm_statistic', 'bp_lm_p_value', 'bp_F_statistic', 'bp_F_statistic_p_value'])
print("BreuschPagan异方差检验结果：")
print(bp_result, "\n")
if bp_result['bp_lm_p_value'] < 0.05:
    print("BreuschPagan异方差检验表明：存在异方差", "\n")
else:
    print("BreuschPagan异方差检验表明：不存在异方差", "\n")