目录
5.1Variance Inflation Factor (VIF)差异性通货膨胀系数(VIF)
1.线性回归
1.1导入所需模块
import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
%matplotlib inline
import math
# imports
import pandas as pd
import numpy as np
# from sklearn.linear_model import LinearRegression # for ml
import statsmodels.api as sm
# statistical tests
from statsmodels.tsa.stattools import adfuller, kpss # stationarity
from statsmodels.stats.api import linear_harvey_collier # linearity
from statsmodels.stats.diagnostic import linear_rainbow # linearity
from statsmodels.stats.outliers_influence import variance_inflation_factor as vif # multicollinearity
from scipy.stats import shapiro, anderson # normality
from statsmodels.stats.stattools import durbin_watson # autocorrelation
from statsmodels.stats.diagnostic import acorr_ljungbox, acorr_breusch_godfrey # autocorrelation
from statsmodels.stats.diagnostic import het_breuschpagan # heteroscedasticity
from statsmodels.stats.api import het_goldfeldquandt # heteroscedasticity
# stat graphics
from statsmodels.graphics.regressionplots import plot_ccpr # partial residuals plot
from statsmodels.graphics.gofplots import qqplot # qq plot for normality of residuals
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # autocorrelation plots
1.2读取数据
# 读取数据,并打印前5行
data = pd.read_csv('data/regdat1.csv', sep=';', header=0, index_col=0)
data.head()
1.3模型参数估计
#建立模型
X = data.iloc[:, 1:]
y = data.iloc[:, 0]
reg_model = sm.OLS(y, X)
reg_model = reg_model.fit()
reg_summary = reg_model.summary()
reg_summary
2.以下一组统计检验用于模型有效性的定量验证。
回归模型和系数的显著性:t准则、F准则。
变量和残差的稳定性:增强Dickey-Fuller检验(ADF)、Phillips-Perron检验(PP)、Kwiatkowski-Phillips-Schmidt-Shin检验(KPSS)。
因变量和自变量之间的依赖线:部分残差图,Ramsey RESET检验。
无多线性:VIF和条件指数。
正常残差:Q-Q图、Shapiro-Wilk检验、Kolmogorov-Smirnov检验、Jarka-Ber检验和Anderson-Darling检验。
无自相关:ACF和PACF图、Darbin-Watson检验、Broisch-Godfrey检验和Young Box检验。
残余物的同质性:残余物与对应值的关系图,Broysch-Pagan试验和Goldfeldt-Quandt试验。
回归系数稳定性:前向和后向稳定性检验,滚动周检验。
下文将详细介绍所有这些测试及其结果。
# 阈值
p_value_threshold = 0.05
confidence_threshold = 0.01
removal_rate = 0.3 # 30%
number_of_out_of_sample_data = 9
2.1显著性/重要性
变量值--T-тест
reg_summary.tables[1]
回归系数在5%的水平上显著。
回归方程的显著性--F-тест
reg_summary.tables[0]
F统计值为71.83,P值=6.1e-16,因此,回归方程显著。
3.固定性
对于一个时间序列模型,所有的因变量、自变量和残差都必须使用ADF检验、PP检验和KPSS检验来检验其稳定性。每个变量必须至少通过2次固定测试。对于每个固定检验,必须对因变量、自变量和残差进行固定性评价。
零平均静止性(仅针对残留物);
单一平均静止性;以及
这是一种趋势性的静止。
确定的趋势静止性是指变量包含趋势( 𝛽𝑡 ),去趋势后过程变得静止。对于模型残差,通常只检查零均值静止性。
变量:
如果一个变量已经通过了单均值静止性检验,就可以直接用于模型中。
如果该变量只通过了趋势检验,那么就需要明确说明为什么可以在模型中使用这个变量。
如果方程中含有一个非稳态变量,但余数通过了稳态检验,则说明变量是共积分的,即有一些多变量的线性组合是稳态的。在这种情况下,只要有使用非稳态变量的理由,变量仍可用于建模。否则,就得重新构建模型。
如果方程中包含一个非稳态变量,且残差也是非稳态的,那么就需要对模型进行重构。
对于ADF和PP检验,零假设是数据序列有单位根,即数据序列是非平稳的;而KPSS检验的零假设是数据序列没有单位根,所以ADF和PP检验的p值宜小,KPSS检验的p值宜大。
def get_stationarity_test_results(series, regime, series_name, p_value_threshold=p_value_threshold):
"""
regime: "none", "mean" or "trend"
"""
if regime == "none":
ADFType = "nc"
KPSSType = "c"
elif regime == "mean":
ADFType = "c"
KPSSType = "c"
elif regime == "trend":
ADFType = "ct"
KPSSType = "ct"
else:
print("There is no such regime name")
return None
adf_result = adfuller(series, regression=ADFType)[1]
kpss_result = kpss(series, regression=KPSSType)[1]
outcome = 'Pass' if (adf_result < p_value_threshold and kpss_result > p_value_threshold) else 'Fail'
result = {'series_name': series_name, 'test_regime': regime, 'ADF p-value': adf_result,
'KPSS p-value': kpss_result, 'outcome': outcome}
return result
# residuals残差
residuals = reg_model.resid
stationarity_residuals = get_stationarity_test_results(residuals, 'none', 'residuals')
stationarity_table = pd.DataFrame([stationarity_residuals])
# data
regimes = ['mean', 'trend']
for column in data.columns:
for regime in regimes:
stationarity_result = get_stationarity_test_results(data[column], regime, column)
stationarity_table = pd.concat([stationarity_table, pd.DataFrame([stationarity_result])], axis=0)
# stationarity tests results table固定性测试结果表
columns_right_order = ['series_name', 'test_regime', 'ADF p-value', 'KPSS p-value', 'outcome']
stationarity_table = stationarity_table[columns_right_order]
stationarity_table.index = list(range(len(data.columns) * 2 +1))
stationarity_table
对模型的变量和残差进行静止性检验。
4.模型的线性度
OLS模型假设自变量和因变量之间具有线性关系。以下是两种验证线性的方法。
部分残差画图
plot_ccpr(reg_model, data.columns[1])
plt.show()
plot_ccpr(reg_model, data.columns[2])
plt.show()
你可以追踪图表上的线性依赖关系。
线性测试
# Rainbow test
rainbow_pvalue = linear_rainbow(reg_model)[1]
rainbow_outcome = 'Pass' if rainbow_pvalue > p_value_threshold else 'Fail'
print('p-value of Rainbow test {}, outcome: {}'.format(rainbow_pvalue, rainbow_outcome))
# Harvey-Collier test
hc_pvalue = linear_harvey_collier(reg_model)[1]
hc_outcome = 'Pass' if hc_pvalue > p_value_threshold else 'Fail'
print('p-value of Harvey-Collier test {}, outcome: {}'.format(hc_pvalue, hc_outcome))
根据得到的p值的值,所选模型的线性规范是正确的。
5.独立变量的多重关联性
多重线性是指在多元回归中,两个或两个以上的独立变量强烈相关的现象。当多线性存在时,受影响的输入变量通常有很大的标准误差。在这种情况下,检查系数的显著性可能会导致拒绝假的零假设(系数与0没有显著差异)。此外,由于标准误差较大,输入数据的微小变化可能会导致模型的巨大变化,甚至导致系数符号的变化,这可能会导致独立变量内多共线性模型的稳定性检验结果不理想。
5.1Variance Inflation Factor (VIF)差异性通货膨胀系数(VIF)
采用VIF检验法对模型中的所有输入变量进行多共线性检验。没有强调VIF这个具体数字。它应该用来指导后续对有无特定变量的模型的检验。VIF越高,说明多线性度越高。VIF超过5即为大。
X_cor = X.corr()
inv_cor_data = pd.DataFrame(np.linalg.inv(X.corr().values), index = X_cor.index, columns=X_cor.columns)
vif_result = pd.DataFrame(data=np.diag(inv_cor_data), index=X.columns, columns=['VIF'])
vif_result
5.2状态指数Condition Index
条件指数(Condition Index)等于最大特征环境值X'X(X--回归变量矩阵)与每个独立回归者的平方根之比的最大值。条件指数值应小于30。
eig_vals = np.linalg.eigvals(np.dot(np.transpose(X), X))
cond_ind_result = pd.DataFrame(data=np.sqrt(max(eig_vals) ** 2 / eig_vals), index=X.columns, columns=['Condition Index'])
cond_ind_result
在模型中没有观察到多线性的迹象。
5.3残留物规范性
在OLS回归中,假设残余物呈正态分布。如果打破了正态性的假设,比率仍然是无偏的,但是,比率重要性的检验可能是不可靠的。
Q-Q Plot for 残差
Quantile-Quantile Plot分位数-四分位数图是一种图形方法,通过构建它们的分位数与标准正态分布相互比较。如果残余物是正常的,那么Q-Q Plot应该大约位于45度的线上。
qqplot(reg_model.resid, line='s')
plt.show()
残留物正常性测试
_, shapiro_pvalue = shapiro(reg_model.resid)
shapiro_pvalue = round(shapiro_pvalue, 2)
shapiro_outcome = 'Pass' if shapiro_pvalue > p_value_threshold else 'Fail'
print('p-value of Shapiro Normality test {}, outcome: {}'.format(shapiro_pvalue, shapiro_outcome))
result = anderson(reg_model.resid)
print('Statistic: %.3f' % result.statistic)
for i in range(len(result.critical_values)):
sl, cv = result.significance_level[i], result.critical_values[i]
if result.statistic < result.critical_values[i]:
print('%.3f: %.3f, data looks normal (fail to reject H0)' % (sl, cv))
else:
print('%.3f: %.3f, data does not look normal (reject H0)' % (sl, cv))
根据检验结果,回归模型的遗迹呈正态分布。
6.残留物的自相关性
OLS假设残余物没有自相关。如果违反了这一假设,关于自变量的显著性的结论就可能不可靠。以下是几种估计自相关的方法。
6.1ACF and PACF plots
ACF和PACF图通常用于可视化给定滞后时的选择性自相关和部分自相关。它们也可以作为一个有用的工具,以确定一个潜在的模型(AR,MA或ARMA)和相应的滞后的建模残留物。图中还显示了一个近似检验的结果,即对于某个阵营来说,(部分)自相关为零。临界区域定义为±2/√n,其中n-数据点的数量。如果(部分)自相关在临界区域内,我们就得出结论,这个阵营中不存在(部分)自相关。如果一个分化行的PACF出现急剧收缩和/或滞后1自相关为正,即如果该行看起来有点 "无分化",我们可以考虑在模型中加入AR项。滞后,其中PACF被切断,是指定的AR成员数。如果一个分化行的ACF出现急剧收缩和/或滞后1自相关为负,即如果该行看起来有点 "过度分化",那么我们可以考虑在模型中加入MA项。ACF被切断的滞后时间是MA的指定成员数。
plot_acf(reg_model.resid, lags=30)
plt.show()
plot_pacf(reg_model.resid, lags=30)
plt.show()
6.2残余自相关检验
# Durbin-Watson test德宾-沃森测试
autocorrelation_results = pd.DataFrame(columns=['Autocorrelation test', 'lag', 'p-value / DW statistic value', 'outcome'])
autocorrelation_results.loc[0] = ['Durbin-Watson', 'NA', durbin_watson(reg_model.resid), 'Pass']
# Breusch-Godfrey test布鲁什-戈弗雷试验
for lag in range(1, 5):
tmp_res = acorr_breusch_godfrey(reg_model, nlags=lag)[3]
autocorrelation_results.loc[lag] = ['Breusch-Godfrey', lag, tmp_res, 'Pass' if tmp_res > p_value_threshold else 'Fail']
# Ljung-Box test
for lag in range(1, 5):
tmp_res = acorr_ljungbox(reg_model.resid, lags=lag)[1][0]
autocorrelation_results.loc[lag + 4] = ['Breusch-Godfrey', lag, tmp_res, 'Pass' if tmp_res > p_value_threshold else 'Fail']
autocorrelation_results
Darbin-Watson检验的检验统计量约为2*(1-r),其中r-采样自相关的残基。因此,对于r==0表示不存在序列相关性,检验统计量等于2。这个统计数字将永远在0到4之间。统计数据越接近0,越能证明序列正相关。越接近4,越能证明负的序列相关性。
因此,回归的残余物不进行自动校正。
7.残留物的异方差性
7.1残差与拟合值对比图
plt.scatter(reg_model.fittedvalues, reg_model.resid)
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.show()
该图对异方差性进行了可视化研究。如果模型残差是随机分散在零附近,而对于异方差没有系统的规律,则可以接受同方差。
# Breush-Pagan test
bp_pvalue = het_breuschpagan(reg_model.resid, X)[3]
bp_outcome = 'Pass' if bp_pvalue > p_value_threshold else 'Fail'
print('p-value of Breush-Pagan test {}, outcome: {}'.format(bp_pvalue, bp_outcome))
# Goldfeld-Quandt test
gq_pvalue = het_goldfeldquandt(reg_model.resid, X)[1]
gq_outcome = 'Pass' if gq_pvalue > p_value_threshold else 'Fail'
print('p-value of Goldfeld-Quandt test {}, outcome: {}'.format(gq_pvalue, gq_outcome))
根据图示和试验结果分析,可以认为残余物是同源的。
8.回归模型系数的稳定性
对于所有的回归方程,都需要检查每个系数的稳定性。核查的目的是为了确保所提出的建模方法是可靠的。测试开始时,先将原数据集的最后30%(前向稳定性)或前30%(后向稳定性)数据删除,作为测试集。在每次迭代时,在测试数据集中增加一个观测值,并重新评估与所研究模型输入变量完全相同的模型。可以通过比较模型系数的符号、数值和置信区间来估计模型的稳定性。预计回归系数变化不大,且会在相应的置信区间内。
8.1前进稳定性
removed_len = int(removal_rate * data.shape[0])
coeff_1_forward = pd.DataFrame(data=np.zeros((removed_len, 4)), columns=['lower', 'rolling', 'upper', 'actual'])
coeff_2_forward = pd.DataFrame(data=np.zeros((removed_len, 4)), columns=['lower', 'rolling', 'upper', 'actual'])
for idx in range(data.shape[0] - removed_len, data.shape[0]):
temp_reg_model = sm.OLS(y[:idx], X.iloc[:idx, :]).fit()
conf_int = temp_reg_model.conf_int(confidence_threshold)
coeff_1_forward.loc[idx - data.shape[0] + removed_len] = [conf_int.iloc[0, 0],
temp_reg_model.params[0], conf_int.iloc[0, 1], reg_model.params[0]]
coeff_2_forward.loc[idx - data.shape[0] + removed_len] = [conf_int.iloc[1, 0],
temp_reg_model.params[1], conf_int.iloc[1, 1], reg_model.params[1]]
coeff_1_forward.plot()
plt.ylabel(X.columns[0])
plt.title('Forward stability of ' + X.columns[0] + ' regression coefficient')
plt.show()
coeff_2_forward.plot()
plt.ylabel(X.columns[1])
plt.title('Forward stability of ' + X.columns[1] + ' regression coefficient')
plt.show()
8.2后向稳定性
coeff_1_backward = pd.DataFrame(data=np.zeros((removed_len, 4)), columns=['lower', 'rolling', 'upper', 'actual'])
coeff_2_backward = pd.DataFrame(data=np.zeros((removed_len, 4)), columns=['lower', 'rolling', 'upper', 'actual'])
for idx in range(removed_len, -1, -1):
temp_reg_model = sm.OLS(y[idx:], X.iloc[idx:, :]).fit()
conf_int = temp_reg_model.conf_int(confidence_threshold)
coeff_1_backward.loc[removed_len - idx] = [conf_int.iloc[0, 0],
temp_reg_model.params[0], conf_int.iloc[0, 1], reg_model.params[0]]
coeff_2_backward.loc[removed_len - idx] = [conf_int.iloc[1, 0],
temp_reg_model.params[1], conf_int.iloc[1, 1], reg_model.params[1]]
coeff_1_backward.plot()
plt.ylabel(X.columns[0])
plt.title('Backward stability of ' + X.columns[0] + ' regression coefficient')
plt.show()
coeff_2_backward.plot()
plt.ylabel(X.columns[1])
plt.title('Backward stability of ' + X.columns[1] + ' regression coefficient')
plt.show()
根据收到的图表,可以得出结论,最终的回归系数在置信区间内,证明了回归系数的稳定性。
9.性能测试
模型的性能测试是用来评估模型预测力的准确性,特别是当被分析的模型已经使用了所有可用数据时。在这种情况下,应该对样本内和样本外数据的模型性能进行评估。性能测试包括以下步骤:
从原始数据集中删除最后9个观测点(通常为最后9个季度的季度观测)。删除的数据--样本外数据。而其他数据--样本内数据。
用相同的变量覆盖所提出的模型方程,只针对样本内区间。
根据样本内和样本外数据计算模型关键性能指标(CERRPCT、MAE/MAR和RMSE/STD)。
对9个点的实际A值(𝑖)和模型预测的P值(𝑖)分别进行比较。试验期各期(季度)的实际结果与预测结果之和确定如下:
预测误差--每个时期(季度或月份)的误差定义如下:
𝐸𝑅𝑅(𝑖)=𝑃(𝑖)−𝐴(𝑖)
累计预测误差--逆向测试期间的累计误差,它是基本指标CERRPCT的分子。
𝐶𝐸𝑅𝑅(𝑁)=𝐶𝑃(𝑁)−𝐶𝐴(𝑁)
预测误差百分比--每个时期(季度或月份)的误差百分比确定如下:
𝐸𝑅𝑅𝑃𝐶𝑇(𝑖)=𝐸𝑅𝑅(𝑖)/𝐴(𝑖)
预测误差的累计百分比--测试期间的累计误差百分比是模型的主要性能指标,其定义如下。
𝐶𝐸𝑅𝑅𝑃𝐶𝑇(𝑖)=𝐶𝐸𝑅𝑅(𝑖)/𝐶𝐴(𝑖)
MAE/MAR和RMSE/SD计算方法如下:
其中n为样本内数据的样本量(对于样本外的n=9)。
in_sample_X = X.iloc[:X.shape[0] - number_of_out_of_sample_data, :]
in_sample_y = y[:len(y) - number_of_out_of_sample_data]
out_of_sample_X = X.iloc[X.shape[0] - number_of_out_of_sample_data:, :]
out_of_sample_y = y[len(y) - number_of_out_of_sample_data:]
# in-sample regression
reg_model_in_sample = sm.OLS(in_sample_y, in_sample_X).fit()
# get predictions
predictions = reg_model_in_sample.get_prediction(out_of_sample_X)
predictions = predictions.summary_frame(alpha=0.05)
performance_df = pd.DataFrame(index=data.index.values, columns=[data.columns[0], 'fitted',
'lower bound', 'upper_bound'])
performance_df.iloc[:, 0] = data.iloc[:, 0]
performance_df.iloc[:data.shape[0] - number_of_out_of_sample_data, 1] = reg_model_in_sample.fittedvalues
performance_df.iloc[data.shape[0] - number_of_out_of_sample_data:, 1] = predictions.iloc[:, 0]
performance_df.iloc[data.shape[0] - number_of_out_of_sample_data:, 2] = predictions.iloc[:, -2]
performance_df.iloc[data.shape[0] - number_of_out_of_sample_data:, 3] = predictions.iloc[:, -1]
performance_df.plot()
plt.show()
# in-sample performance metrics
n = data.shape[0] - number_of_out_of_sample_data
mae_in_sample = (1 / n) * sum([abs(in_sample_y[idx] - reg_model_in_sample.fittedvalues[idx]) for idx in range(n)])
mar_in_sample = (1 / (n - 1)) * sum([abs(in_sample_y[idx] - in_sample_y[idx - 1]) for idx in range(1, n)])
rmse_in_sample = math.sqrt((1 / n) * sum([abs(in_sample_y[idx] -
reg_model_in_sample.fittedvalues[idx])**2 for idx in range(n)]))
sd_in_sample = math.sqrt((1 / n) * sum([abs(in_sample_y[idx] -
np.mean(in_sample_y))**2 for idx in range(n)]))
ceerpct_in_sample = abs((reg_model_in_sample.fittedvalues.sum() - in_sample_y.sum()) / in_sample_y.sum())
print('MAE in sample: ', round(mae_in_sample, 3))
print('MAR in sample: ', round(mar_in_sample, 3))
print('MAE / MAR in sample: ', round(mae_in_sample / mar_in_sample, 3))
print('RMSE in sample: ', round(rmse_in_sample, 3))
print('SD in sample: ', round(sd_in_sample, 3))
print('RMSE / SD in sample: ', round(rmse_in_sample / sd_in_sample, 3))
print('CEERPCT in sample: ', round(ceerpct_in_sample, 3))
# out-of-sample performance metrics
n = number_of_out_of_sample_data
mae_out_of_sample = (1 / n) * sum([abs(out_of_sample_y[idx] - predictions.iloc[idx, 1]) for idx in range(n)])
mar_out_of_sample = (1 / (n - 1)) * sum([abs(out_of_sample_y[idx] - out_of_sample_y[idx - 1]) for idx in range(1, n)])
rmse_out_of_sample = math.sqrt((1 / n) * sum([abs(out_of_sample_y[idx] -
predictions.iloc[idx, 1])**2 for idx in range(n)]))
sd_out_of_sample = math.sqrt((1 / n) * sum([abs(out_of_sample_y[idx] -
np.mean(out_of_sample_y))**2 for idx in range(n)]))
ceerpct_out_sample = abs((predictions.iloc[:, 1].values.sum() - out_of_sample_y.sum()) / out_of_sample_y.sum())
print('MAE out of sample: ', round(mae_out_of_sample, 3))
print('MAR out of sample: ', round(mar_out_of_sample, 3))
print('MAE / MAR out of sample: ', round(mae_out_of_sample / mar_out_of_sample, 3))
print('RMSE out of sample: ', round(rmse_out_of_sample, 3))
print('SD out of sample: ', round(sd_out_of_sample, 3))
print('RMSE / SD inout ofsample: ', round(rmse_out_of_sample / sd_out_of_sample, 3))
print('CEERPCT out of sample: ', round(ceerpct_out_sample, 3))
Performance 模式是可以接受的。
10.输入数据检查
10.1视觉检测
for column in data.columns:
data[column].plot()
plt.legend(data.columns)
plt.show()
10.2基本统计
# descriptive statistics
data.describe()
10.3高影响点
data.plot.box()
plt.show()
fig, ax = plt.subplots(figsize=(12,8))
fig = sm.graphics.influence_plot(reg_model, ax=ax, criterion="cooks")
11.影响程度
from statsmodels.stats import outliers_influence
result = outliers_influence.OLSInfluence(reg_model).summary_frame()
result
以下为部分截取
11.1DFFITS
data.shape
result[['dffits_internal']].plot();
11.2DFBETA
result[['dfb_lags']].plot();
result[['dfb_O_M3_US']].plot();
11.3 Cook's distance
result[['cooks_d']].plot();
【参考文献】
1.Изучаем pandas. Урок 2. Структуры данных Series и DataFrame