第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 變數單位的選擇
多個變數的單位，應該儘可能的避免變數間數量級差別過大。
直接取對數。