Python之 【模型建立和测试-模型测试模板】

目录

1.线性回归

1.1导入所需模块

1.2读取数据

1.3模型参数估计

2.以下一组统计检验用于模型有效性的定量验证。

2.1显著性/重要性

3.固定性

4.模型的线性度

5.独立变量的多重关联性

5.1Variance Inflation Factor (VIF)差异性通货膨胀系数(VIF)

5.2状态指数Condition Index

5.3残留物规范性

6.残留物的自相关性

6.1ACF and PACF plots

6.2残余自相关检验

7.残留物的异方差性

7.1残差与拟合值对比图

8.回归模型系数的稳定性

8.1前进稳定性

8.2后向稳定性

9.性能测试

10.输入数据检查

10.1视觉检测

10.2基本统计

10.3高影响点

11.影响程度

11.1DFFITS

11.2DFBETA

11.3 Cook's distance


 

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值