利用python实现Diebold-Mariano检验
文章目录
- DM检验的原理
- 代码实现
- 函数说明
- 实例
- 本文参考
DM检验的原理
给定两个预测的预测结果，我们希望比较他们的预测结果，以用于模型预测精度的比较。
Diebold-Mariano检验本质是一个t检验，用于检验产生预测的两个损失序列的平均值是否相等。即，它是一系列损失差的零均值的t检验。
- 原假设：DM统计量均值为0，即两个模型的预测效率一致。
- 备择假设：两个模型的预测效率不一致。
**注意：**在使用DM检验式时，其假设损失序列是平稳的。
另外,DM检验在小样本数据时往往会拒绝零假设。对于小样本数据，推荐Harvey, Leybourne and Newbold (HLN)检验【1】；
代码实现
该代码来源于：https://github.com/johntwk/Diebold-Mariano-Test
# Author : John Tsang
# Date : December 7th, 2017
# Purpose : Implement the Diebold-Mariano Test (DM test) to compare
# forecast accuracy
# Input : 1) actual_lst: the list of actual values
# 2) pred1_lst : the first list of predicted values
# 3) pred2_lst : the second list of predicted values
# 4) h : the number of stpes ahead
# 5) crit : a string specifying the criterion
# i) MSE : the mean squared error
# ii) MAD : the mean absolute deviation
# iii) MAPE : the mean absolute percentage error
# iv) poly : use power function to weigh the errors
# 6) poly : the power for crit power
# (it is only meaningful when crit is "poly")
# Condition: 1) length of actual_lst, pred1_lst and pred2_lst is equal
# 2) h must be an integer and it must be greater than 0 and less than
# the length of actual_lst.
# 3) crit must take the 4 values specified in Input
# 4) Each value of actual_lst, pred1_lst and pred2_lst must
# be numerical values. Missing values will not be accepted.
# 5) power must be a numerical value.
# Return : a named-tuple of 2 elements
# 1) p_value : the p-value of the DM test
# 2) DM : the test statistics of the DM test
##########################################################
# References:
#
# Harvey, D., Leybourne, S., & Newbold, P. (1997). Testing the equality of
# prediction mean squared errors. International Journal of forecasting,
# 13(2), 281-291.
#
# Diebold, F. X. and Mariano, R. S. (1995), Comparing predictive accuracy,
# Journal of business & economic statistics 13(3), 253-264.
#
##########################################################
def dm_test(actual_lst, pred1_lst, pred2_lst, h = 1, crit="MSE", power = 2):# Routine for checking errorsdef error_check():rt = 0msg = ""# Check if h is an integerif (not isinstance(h, int)):rt = -1msg = "The type of the number of steps ahead (h) is not an integer."return (rt,msg)# Check the range of hif (h < 1):rt = -1msg = "The number of steps ahead (h) is not large enough."return (rt,msg)len_act = len(actual_lst)len_p1 = len(pred1_lst)len_p2 = len(pred2_lst)# Check if lengths of actual values and predicted values are equalif (len_act != len_p1 or len_p1 != len_p2 or len_act != len_p2):rt = -1msg = "Lengths of actual_lst, pred1_lst and pred2_lst do not match."return (rt,msg)# Check range of hif (h >= len_act):rt = -1msg = "The number of steps ahead is too large."return (rt,msg)# Check if criterion supportedif (crit != "MSE" and crit != "MAPE" and crit != "MAD" and crit != "poly"):rt = -1msg = "The criterion is not supported."return (rt,msg) # Check if every value of the input lists are numerical valuesfrom re import compile as re_compilecomp = re_compile("^\d+?\.\d+?$") def compiled_regex(s):""" Returns True is string is a number. """if comp.match(s) is None:return s.isdigit()return Truefor actual, pred1, pred2 in zip(actual_lst, pred1_lst, pred2_lst):is_actual_ok = compiled_regex(str(abs(actual)))is_pred1_ok = compiled_regex(str(abs(pred1)))is_pred2_ok = compiled_regex(str(abs(pred2)))if (not (is_actual_ok and is_pred1_ok and is_pred2_ok)): msg = "An element in the actual_lst, pred1_lst or pred2_lst is not numeric."rt = -1return (rt,msg)return (rt,msg)# Error checkerror_code = error_check()# Raise error if cannot pass error checkif (error_code[0] == -1):raise SyntaxError(error_code[1])return# Import librariesfrom scipy.stats import timport collectionsimport pandas as pdimport numpy as np# Initialise listse1_lst = []e2_lst = []d_lst = []# convert every value of the lists into real valuesactual_lst = pd.Series(actual_lst).apply(lambda x: float(x)).tolist()pred1_lst = pd.Series(pred1_lst).apply(lambda x: float(x)).tolist()pred2_lst = pd.Series(pred2_lst).apply(lambda x: float(x)).tolist()# Length of lists (as real numbers)T = float(len(actual_lst))# construct d according to critif (crit == "MSE"):for actual,p1,p2 in zip(actual_lst,pred1_lst,pred2_lst):e1_lst.append((actual - p1)**2)e2_lst.append((actual - p2)**2)for e1, e2 in zip(e1_lst, e2_lst):d_lst.append(e1 - e2)elif (crit == "MAD"):for actual,p1,p2 in zip(actual_lst,pred1_lst,pred2_lst):e1_lst.append(abs(actual - p1))e2_lst.append(abs(actual - p2))for e1, e2 in zip(e1_lst, e2_lst):d_lst.append(e1 - e2)elif (crit == "MAPE"):for actual,p1,p2 in zip(actual_lst,pred1_lst,pred2_lst):e1_lst.append(abs((actual - p1)/actual))e2_lst.append(abs((actual - p2)/actual))for e1, e2 in zip(e1_lst, e2_lst):d_lst.append(e1 - e2)elif (crit == "poly"):for actual,p1,p2 in zip(actual_lst,pred1_lst,pred2_lst):e1_lst.append(((actual - p1))**(power))e2_lst.append(((actual - p2))**(power))for e1, e2 in zip(e1_lst, e2_lst):d_lst.append(e1 - e2) # Mean of d mean_d = pd.Series(d_lst).mean()# Find autocovariance and construct DM test statisticsdef autocovariance(Xi, N, k, Xs):autoCov = 0T = float(N)for i in np.arange(0, N-k):autoCov += ((Xi[i+k])-Xs)*(Xi[i]-Xs)return (1/(T))*autoCovgamma = []for lag in range(0,h):gamma.append(autocovariance(d_lst,len(d_lst),lag,mean_d)) # 0, 1, 2V_d = (gamma[0] + 2*sum(gamma[1:]))/TDM_stat=V_d**(-0.5)*mean_dharvey_adj=((T+1-2*h+h*(h-1)/T)/T)**(0.5)DM_stat = harvey_adj*DM_stat# Find p-valuep_value = 2*t.cdf(-abs(DM_stat), df = T - 1)# Construct named tuple for returndm_return = collections.namedtuple('dm_return', 'DM p_value')rt = dm_return(DM = DM_stat, p_value = p_value)return rt函数说明
Input Parameter Description：
- actual_lst The list of actual values
- pred1_lst The first list of predicted values
- pred2_lst The second list of predicted values
- h The number of steps ahead of the prediction. The default value is 1.
- crit A string specifying the criterion. The default value is MSE.
MSE : d = (e1)^2 - (e2)^2
MAD : d = abs(e1) - abs(e2)
MAPE: d = abs((e1 - actual)/(actual))
Poly: d = (e1)^power - (e2)^power - power The power for crit equal “poly”. It is only meaningful when crit is “poly”. The default value is 2. (i.e. E[d] = (e1)^2 - (e2)^2)
Output:
7. DM The DM test statistics
8. p-value The p-value of DM test statistics
实例
from dm_test import dm_test
import randomrandom.seed(123)
actual_lst = range(0,100)
pred1_lst = range(0,100)
pred2_lst = range(0,100)actual_lst = random.sample(actual_lst,100)
pred1_lst = random.sample(pred1_lst,100)
pred2_lst = random.sample(pred2_lst,100)rt = dm_test(actual_lst,pred1_lst,pred2_lst,h = 1, crit="MAD")
print rt
rt = dm_test(actual_lst,pred1_lst,pred2_lst,h = 1, crit="MSE")
print rt
rt = dm_test(actual_lst,pred1_lst,pred2_lst,h = 1, crit="poly", power=4)
print rt
#dm_return(DM=1.3275742446369585, p_value=0.18737195617455585)
#dm_return(DM=1.2112523589452902, p_value=0.22868210381769466)
#dm_return(DM=0.9124498079287283, p_value=0.36374861695187799)本文参考
【1】https://www.real-statistics.com/time-series-analysis/forecasting-accuracy/diebold-mariano-test/
【2】https://blog.csdn.net/weixin_43948357/article/details/112884095
【3】https://blog.csdn.net/qq_40283970/article/details/104644122
【4】 Harvey, D., Leybourne, S., & Newbold, P. (1997). Testing the equality of prediction mean squared errors. International Journal of forecasting,
13(2), 281-291.
【5】 Diebold, F. X. and Mariano, R. S. (1995), Comparing predictive accuracy, Journal of business & economic statistics 13(3), 253-264.
利用python实现Diebold-Mariano检验相关推荐
- 利用python进行单边T检验
可以利用 python 中的 scipy.stats.ttest_ind 做关于两组数据的双边 t 检验,结果比较简单.但是做 大于或者小于的单边检测的时候需要做一些处理,才能得到正确的结果. fro ...
- python怎么实现检验_[python skill]利用python实现假设性检验方法
[python skill]利用python实现假设性检验方法 刀尔東 2018-08-03 09:19:13 1244 收藏 2 版权 hello,大噶好,最近新学习了利用python实现假设性检验 ...
- 统计系列（四）利用Python进行假设检验
统计系列(四)利用Python进行假设检验 z检验 主要应用场景:在大样本量的总体比例检验 核心:两样本的总体比例差异 单样本比例检验 # 检验样本合格率与0.38是否有差异 import numpy ...
- python 白噪声检验-利用python实现平稳时间序列的建模方式
假如某个观察值序列通过序列预处理可以判定为平稳非白噪声序列,就可以利用ARMA模型对该序列进行建模.建模的基本步骤如下: (1)求出该观察值序列的样本自相关系数(ACF)和样本偏自相关系数(PACF) ...
- 利用python实现非参数方法（拟合优度检验）
简介 数据不符合正态分布时,无法直接用t分布和F分布及正态分布来检验数据的显著性.如果检验对于数据分布没有要求,这类检验统称为非参数检验.本文将介绍其中的一种拟合优度检验. 拟合优度检验-期望频数相等 ...
- json字段顺序读取 python_如何利用Python批量读取视频文件的时间长度？
本期的主题是利用Python来实现对视频文件时间长度的读取. 在学习编程语言时,相比较于通过书本来学习知识,我更喜欢通过观看学习视频的方式来进行学习,通过主讲老师的讲解,我能很直观且快速的了解一些知识 ...
- python逐步回归筛选变量_利用python实现逐步回归
逐步回归的基本思想是将变量逐个引入模型,每引入一个解释变量后都要进行F检验,并对已经选入的解释变量逐个进行t检验,当原来引入的解释变量由于后面解释变量的引入变得不再显著时,则将其删除.以确保每次引入新 ...
- 数据分析 python 用途-利用Python数据分析可以实现些什么功能呢？
随着大数据时代的来临和Python编程语言的火爆,Python数据分析早已成为现在职场人的必备核心技能.那么利用Python数据分析可以做什么呢?简单来说,可以做到的内容有很多,比如检查数据表.数据表 ...
- python 数据分析学什么-利用Python做数据分析 需要学习哪些知识
根据调查结果,十大最常用的数据工具中有八个来自或利用Python.Python广泛应用于所有数据科学领域,包括数据分析.机器学习.深度学习和数据可视化.不过你知道如何利用Python做数据分析吗?需要 ...
- python怎么判断一个文件是否存在-利用Python如何判断一个文件是否存在
通常在读写文件之前,需要判断文件或目录是否存在,不然某些处理方法可能会使程序出错.所以最好在做任何操作之前,先判断文件是否存在. 这里将介绍三种判断文件或文件夹是否存在的方法,分别使用os模块.Try ...
最新文章
- facerec = dlib.face_recognition_model_v1()面部识别器用法
- python判断能否组成三角形_牛鹭学院：学员笔记|Python： 输入三条边，判断是否可以成为三角形...
- VS Code 的常用快捷键和插件
- 【Thymeleaf】 循环固定次数/循环次数由变量控制
- 【Task】- JVM逃逸分析等待学习任务
- “80后”作家应扮演更重要的角色
- P4301-[CQOI2013]新Nim游戏【线性基】
- java文件学生_文件存储学生信息（JavaIO流）
- mysql locked myisam_MySql 事务 隔离级别 知识点
- C++求100以内中的所有素数
- 摩斯电码php源码,PHP实现基于文本的莫斯电码生成器
- python列表元素可以重复吗_Python列表中的元素重复
- pandas 排序 excel
- 1.《天空之城》- 尤克里里指弹入门版
- 小程序获取用户openid，php获取微信小程序openid的方法
- 元镜头——手机相机的下一场革命
- Cisco思科IPS签名策略配置引擎告警和日志动作
- 解决 plt.savefig 保存图像为全白或全黑图像方法
- Cygwin、Msys、MinGW、Msys2的区别与联系
- 总谐波失真80_如何将总谐波失真（THD）控制着10%以下？