09_模型设定与数据问题
第9章 模型设定与数据问题
如果模型设定不当，会带来设定误差（specification error）
- [[#9.1 遗漏变量|9.1 遗漏变量]]
- [[#9.2 无关变量|9.2 无关变量]]
- [[#9.3 建模策略：“由小到大”还是“由大到小”？|9.3 建模策略：“由小到大”还是“由大到小”？]]
- [[#9.4 解释变量个数的选择|9.4 解释变量个数的选择]]
- [[#9.5 对函数形式的检验|9.5 对函数形式的检验]]
另外，数据本事也可能存在问题: - [[#9.6 多重共线性|9.6 多重共线性]]
- [[#9.7 极端数据|9.7 极端数据]]
- [[#9.8 虚拟变量|9.8 虚拟变量]]
- [[#9.9 经济结构变动的检验|9.9 经济结构变动的检验]]
- [[#9.10 缺失数据与线性插值|9.10 缺失数据与线性插值]]
- [[#9.11 变量单位的选择|9.11 变量单位的选择]]
9.1 遗漏变量
由于某些数据难以获取，遗漏变量现象几乎难以避免。遗漏变量是否一定导致不一致的估计？
- 遗漏变量\(x_2\)与解释变量\(x_1\)不相关：OLS可一致估计，但扰动项方差增大
- 遗漏变量\(x_2\)与解释变量\(x_1\)相关：OLS估计不一致
解决遗漏变量的方法：
- 加入尽可能多的控制变量
- 随机实验与自然实验
- 工具变量法
- 面板数据
9.2 无关变量
与遗漏变量相反，加入了与被解释变量无关的变量。
- OLS仍然一致
- 估计量的方差会增大
9.3 建模策略
“由小到大”还是“由大到小”？
- 折中，凭感觉来。
9.4 解释变量个数的选择
好的经济理论应该能用简洁的模型很好的描述复杂的经济现实。但解释力（增大拟合优度）和简洁性（parismony）是两个矛盾的目标，需要如下方法进行权衡：
校正可决系数
选择K，使 \(\overline R^2\) 最小
赤池信息准则
（AIC，Akaike Information Criterion）
- 选择解释变量的个数，使目标函数最小$$\min_{K}AIC\equiv\ln(\frac{SSR}{n}+\frac{2}{n}K)$$
贝叶斯信息准则
（BIC，Bayesian Information Criterion）/施瓦茨信息准则（SIC）
- 选择解释变量的个数K，使目标函数最小$$\min_{K}BIC\equiv\ln(\frac{SSR}{n}+\frac{\ln n}{n}K)$$
由大到小的序贯t规则
常用于时间序列模型
- 指定一个最大滞后期\(P_{max}\)，
- 令\(\hat P \equiv P_{max}\)，进行估计，对最后一阶系数做t检验。如接受系数为0，则
- 同时还可以观察AIC和BIC的变化
- 令\(\hat P \equiv P_{max}-1\)，进行估计，对最后一阶系数做t检验。如显著，则停止，
- 否则依次类推。
9.5 对函数形式的检验
很多经济关系是非线性的。如果存在非线性项，但遗漏了，是模型设定误差（specification error）的一种形式。可使用如下检验方法：
Ramsey's RESET检验
Regression Equation Specification Error Test
基本思想：如怀疑非线性项被遗漏，则加入非线性项，并检验其系数是否显著。
- 既可以接受被解释变量X高次项
- 也可以接受解释变量y的高次项
Python实现
![[statsmodel_docs#函数形式检验（RESET检验）]]
提供一个一次性检测被解释变量、解释变量和所有解释变量的2次项的函数：
import pandas as pd
import statsmodels.api as sm
def reset(dataset, X_cols, y_col):
"""
对拟合值、全体解释变量、不同的解释变量进行RESET检测.
并显示可识别检测对象的判断结果
"""
X = dataset[X_cols]
y = dataset[y_col]
X = sm.add_constant(X)
Results_y = sm.OLS(y, X).fit()
reset_y = sm.stats.diagnostic.linear_reset(Results,
power=[2,3,4],
use_f=True)
reset_X = sm.stats.diagnostic.linear_reset(Results,
test_type='exog',
power=[2,3,4],
use_f=True)
reset_prin = sm.stats.diagnostic.linear_reset(Results,
test_type='princomp',
power=[2,3,4],
use_f=True)
print("被解释变量的RESET检测结果:",reset_y)
print("全体解释变量的RESET检测结果:",reset_X)
print('解释变量主成分的RESET检测结果:',reset_prin)
x2_list = []
for i in X_cols:
i2 = i+'2'
dataset[i2] = dataset[i]**2
X = dataset[X_cols+[i2]]
# print(X_cols+[i2])
Results_new = sm.OLS(y, sm.add_constant(X)).fit()
print(Results_new.params[i2])
reset_new = sm.stats.diagnostic.linear_reset(Results_new,
test_type='princomp',
power=[2,3,4],
use_f=True)
print(f'加入{i}的2次项的RESET检测结果:',reset_new.pvalue)
if reset_new.pvalue < 0.05 and Results_new.params[i2] > 0.05:
x2_list.append(i2)
print(f'加入2次项后,可考虑引入{x2_list}:')
代码: [[Chapter_09.ipynb]]
9.6 多重共线性
定义
定义 严格多重共线性
某解释变量可由其他解释变量线性表出。数据矩阵X不满秩，\(X'X\)不可逆。
定义 多重共线性
将第k个解释变量对其余解释变量进行回归，可决系数 \(R^2\) 较高
- OLS仍是BLUE
- 表现：
- 回归方程 \(R^2\) 较大，F检验也很显著，但是对单个系数的t检验却不显著
- 增减解释变量使得系数估计值发生较大变化
检验方法
定义 方差膨胀因子（VIF，Variance Inflation Factor）
\[VIF_k \equiv \frac{1}{1-R^2}
\]
可以证明$$Var(\hat\beta_k|X)=\frac{1}{1-R2}·\frac{\sigma2}{S_k}=VIF_k\frac{\sigma^2}{S_k}$$
- \(\sigma^2\equiv Var(\epsilon)\) 扰动项方差
- \(S_k \equiv \sum_{i=1}^n(x_{ik}-\overline x_k)^2\) 为\(x_k\)的离差平方和
title:经验
$VIF_k$ 越大，说明$x_k$ 的多重共线性越严重。在判断是否存在多重共线性的一个经验规则是：
- $\{VIF_1,\cdots,VIF_k \}$ 的最大值不应超过10
处理方法
- 如果不关系具体的回归系数，只关心整个方程预测能力
- 不必理会多重共线性，因为方程是显著的
- 如果关系具体的回归系数，但多重共线性不影响所关心的那个变量
- 不必理会
- 如果多重共线性影响到所关心变量的显著性，则应设法处理
- 增大样本容量
- 剔出导致严重共线性的变量
- 将变量标准化：减去均值，除以标准差
- 修改模型设定
解释变量之间存在相关性是普遍存在的，在一定程度上也是允许的。最常见的处理方法是“无为而治”。
python实现
![[statsmodel_docs#statsmodels.stats.outliers_influence.variance_inflation_factor]]
因为python只提供了一次计算一个变量的vif值，不便于观察，特写如下函数。
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
def vif(exog, criterion=5):
'''vif df格式输出
计算解释变量的方差膨胀因子和容忍度,其中解释变量不含有常数项
返回：
解释变量vif值、容忍度的dataframe
Arguments:
exog -- 解释变量:dataframe
criterion -- 方差膨胀因子阈值,默认5
'''
multicol = pd.DataFrame()
multicol['变量'] = exog.columns
multicol['方差膨胀因子'] = [variance_inflation_factor(exog.values, i) for i in range(exog.shape[1])]
multicol['容忍度'] = 1/multicol['方差膨胀因子']
_ = []
for i in multicol['方差膨胀因子']:
if i > criterion:
_.append('**是**')
else:
_.append('否')
multicol['是否多重共线'] = _
return multicol
案例 [[Chapter_09.ipynb]]
9.7 极端数据
定义
定义 极端值、离群值、高影响力数据
如果样本中有少数观测值离大多数观测值很远，可能对OLS的回归系数产生很大影响。
处理方法
- 检查数据是否正确录入
- 看极端值的背景，必要时直接删除
- 同时报告含或不含极端值的回归结果
9.8 虚拟变量
定义 虚拟变量
取值为0或1的变量
- 用于定性数据
- 或分类数据
例子：
- 性别：$$D =
\begin{cases}
1, \ 男\
0, \ 女
\end{cases}$$ - 五大洲：则需要4个虚拟变量
\[\begin{array}{rcl}
D_1 =
\begin{cases}
1, \ 亚洲\\
0, \ 其他
\end{cases}
,\quad D_2 =
\begin{cases}
1, \ 美洲\\
0, \ 其他
\end{cases}
,\quad D_3 =
\begin{cases}
1, \ 欧洲\\
0, \ 其他
\end{cases}
,\quad D_4 =
\begin{cases}
1, \ 非洲\\
0, \ 其他
\end{cases}
\end{array}\]
- 当$D_1=D_2=D_3=D_4=0$， 则为大洋洲
定义 虚拟变量陷阱
在有常数项的模型中，如果定性指标有M类，最多只能放入（M-1）个虚拟变量。
引入虚拟变量相当于在不同时期给予不同的截距项
引入虚拟变量和交互项相当于在不同时期给予不同的截距项和斜率
9.9 经济结构变动的检验
对于时间序列，模型系数的稳定性很重要。如果存在结构变动，但没考虑，也是一种模型设定误差。
邹检验 和 虚拟变量法 结果一致，虚拟变量法：
- 方便
- 可在存在异方差和自相关的情况下使用，只要在估计时，使用HAC标准误即可。
- 还可提供截距项和斜率变动的信息
1. 邹检验
（1）画图
import pandas as pd
import statsmodels.api as sm
import seaborn as sns
import matplotlib.pyplot as plt
# Load the data
consump = pd.read_stata('../2_Data/Data-2e/consumption.dta')
## 画图
sns.set_theme(style="darkgrid")
sns.lineplot(x='year', y='y', data=consump,marker='^',markers=True)
sns.lineplot(x='year', y='c', data=consump,color='r',marker='o',markers=True)
plt.axvline(x=1992, color='g', linestyle='-')
breakpoint = consump[consump['year']==1992].index[0]
for i in consump.columns:
print(consump[i].astype('float64')**2)
![[9-9-1邹检验-画图.png]]
（2）检验
import pandas as pd
import statsmodels.api as sm
import seaborn as sns
import scipy.stats as stats
def chow_test(dataset, X_cols, y_col, breakpoint):
'''chow_test 对某个回归进行邹志庄检验.
原假设：分段与总体的回归系数一致，不存在结构变化
Arguments:
dataset -- dataframe,变量数据
X_cols -- list,自变量列名
y_col -- str,因变量列名
breakpoint -- float,用于指定分界点.
Returns:
F -- float,邹志庄检验的F值.
p_value -- float,邹志庄检验的p值.
result -- str,邹志庄检验的结果.
'''
dataset1 = dataset[:breakpoint]
dataset2 = dataset[breakpoint:]
ssr_list = []
model_res = []
for d in [dataset, dataset1, dataset2]:
X = d[X_cols].astype('float64')
y = d[y_col].astype('float64')
X = sm.add_constant(X)
model = sm.OLS(y, X)
results = model.fit()
ssr_list.append(results.ssr)
model_res.append(results)
n = dataset.shape[0]
k = len(X_cols) + 1 # 加入常数项的自由度
F = ((n-2*k)/k)*((ssr_list[0]-ssr_list[1]-ssr_list[2])/(ssr_list[1]+ssr_list[2]))
p_value = 1 - stats.f.cdf(F, k, n-2*k)
if p_value < 0.05:
result = '拒绝原假设，存在结构变化'
else:
result = '接受原假设，不存在结构变化'
return F,p_value,result
2.引入虚拟变量
（1）F检验进行联合显著性检验
在使用f_test()时，需要将原假设设置正确，否则自由度会出现偏差，导致计算结果出错。
import numpy as np
# (4) rain 和 coast 联合显著性
hypotheses = '(d = 0), (yd = 0)'
f_test = results.f_test(hypotheses)
f_test.summary()
if f_test.pvalue < 0.05:
print("因p值为{:.4f}小于显著性水平,拒绝原假设".format(f_test.pvalue))
else:
print("因p值为{:.4f}大于显著性水平,不能拒绝原假设".format(f_test.pvalue))
print(f_test.summary())
此结果需要在同方差和自相关的情况下成立，因此还要使用：
- 异方差：White检验
- 自相关：BG检验
（2）异方差检验：white检验
from statsmodels.stats.diagnostic import het_white
result_white = het_white(resid=results.resid, exog=X)
result_white
拒绝原假设：存在异方差
（3）自相关检验：BG检验
from statsmodels.stats.diagnostic import acorr_breusch_godfrey
lag = 1 # 假设我们检验1阶自相关
bg_result = acorr_breusch_godfrey(results, nlags=lag)
print("BG检验的LM统计量:{:.4f}".format(bg_result[0]))
print("BG检验的LM统计量P值:{:.4f}".format(bg_result[1]))
if bg_result[1]<0.05:
print("拒绝原假设,模型存在自相关.")
else:
print("接受原假设，模型不存在自相关.")
BG检验的LM统计量:21.4828
BG检验的LM统计量P值:0.0000
拒绝原假设,模型存在自相关.
（4）NW法
OLS+异方差自相关稳健的标准误（HAC）
nw_cov_type = 'HAC'
p = int((len(consp)** 0.25 // 1) + 1)
nw_kwargs = {'maxlags': p,'use_correction':True}
nw_res = results.get_robustcov_results(cov_type=nw_cov_type, use_t=True, **nw_kwargs)
print(nw_res.summary())
（5）对新回归进行联合显著性检验
hypotheses = '(d = 0), (yd = 0)'
f_test = nw_res.f_test(hypotheses)
if f_test.pvalue < 0.05:
print("因p值为{:.4f}小于显著性水平,拒绝原假设".format(f_test.pvalue))
else:
print("因p值为{:.4f}大于显著性水平,不能拒绝原假设".format(f_test.pvalue))
print(f_test.summary())
因p值为0.0000小于显著性水平,拒绝原假设
F test: F=73.05469132974434, p=1.17872061607885e-12, df_denom=32, df_num=2
9.10 缺失数据与线性插值
中间缺的数据，用两边的直接算平均值替代。
9.11 变量单位的选择
多个变量的单位，应该尽可能的避免变量间数量级差别过大。
直接取对数。
