Действительно ли ваша новая модель лучше?
Как специалисты по данным и экономисты, мы часто сталкиваемся с ситуациями, когда нам нужно сравнить точность двух моделей прогнозирования. Тест Diebold Mariano, также известный как тест DM, представляет собой статистический инструмент, который позволяет нам сделать именно это. В этой статье мы рассмотрим основы теста DM и посмотрим, как его реализовать на Python.
Но прежде чем мы углубимся в детали, давайте сначала разберемся с проблемой, которую пытается решить тест DM. Предположим, у нас есть две модели прогнозирования, модель A и модель B, которые используются для прогнозирования будущей стоимости временного ряда. Мы хотим определить, какая из этих моделей является более точной. Один из способов сделать это — сравнить ошибки прогноза двух моделей. Если ошибки модели А меньше ошибок модели Б, можно сделать вывод, что модель А более точна.
Но как измерить ошибки прогноза? Один из вариантов — использовать среднеквадратичную ошибку (MSE), которая определяется как среднее квадратов разностей между истинными значениями и прогнозируемыми значениями. Чем меньше MSE, тем лучше прогноз. Однако у этого подхода есть существенный недостаток: MSE учитывает только величину ошибок, но не их направление.
Здесь на помощь приходит тест DM. Это вариант классического парного t-критерия, который сравнивает средние значения двух парных выборок. В нашем случае выборки представляют собой ошибки прогноза двух моделей. Тест DM вычисляет тестовую статистику, называемую статистикой DM, которая измеряет стандартизированный дифференциал потерь двух моделей. Если статистика DM значительно отличается от нуля, мы можем сделать вывод, что одна модель точнее другой. Формально,
Теперь, когда у нас есть общее представление о тесте DM, давайте посмотрим, как его реализовать на Python.
import numpy as np import pandas as pd from statsmodels.tsa.ar_model import AutoReg as AR import pandas_datareader as web import datetime as dt from scipy.stats import t import collections START = dt.datetime(1960, 1, 1) END = dt.datetime(2022, 8, 1) # US inflation time series cpi = web.DataReader(["CPIAUCSL"], "fred", START, END) # Ensure stationarity cpi = np.log(cpi).diff(1) cpi.dropna(inplace=True) forecast_A = [] forecast_B = [] h = 1 # Horizon of prediction num_oos = 5 # Testset length # Recursivly estimate the model and predict the next observation for i in range(-num_oos, 0, h): # Fit the model model_A = AR(cpi[:i], lags=1).fit() model_B = AR(cpi[:i], lags=3).fit() # Predict the next value pred_A = model_A.forecast(1) pred_B = model_B.forecast(1) forecast_A.append(pred_A) forecast_B.append(pred_B) # Testset oos_set = cpi.iloc[-num_oos:].squeeze().to_list() def dm_test(real_values, pred1, pred2, h=1, harvey_adj=True): e1_lst = [] e2_lst = [] d_lst = [] real_values = pd.Series(real_values).apply(lambda x: float(x)).tolist() pred1 = pd.Series(pred1).apply(lambda x: float(x)).tolist() pred2 = pd.Series(pred2).apply(lambda x: float(x)).tolist() # Length of forecasts T = float(len(real_values)) # Construct loss differential according to error criterion (MSE) for real, p1, p2 in zip(real_values, pred1, pred2): e1_lst.append((real - p1)**2) e2_lst.append((real - p2)**2) for e1, e2 in zip(e1_lst, e2_lst): d_lst.append(e1 - e2) # Mean of loss differential mean_d = pd.Series(d_lst).mean() # Calculate autocovariance def autocovariance(Xi, N, k, Xs): autoCov = 0 T = float(N) for i in np.arange(0, N-k): autoCov += ((Xi[i+k])-Xs)*(Xi[i]-Xs) return (1/(T))*autoCov # Calculate the denominator of DM stat gamma = [] for lag in range(0, h): gamma.append(autocovariance(d_lst, len(d_lst), lag, mean_d)) # 0, 1, 2 V_d = (gamma[0] + 2*sum(gamma[1:]))/T # Calculate DM stat DM_stat = V_d**(-0.5)*mean_d # Calculate and apply Harvey adjustement # It applies a correction for small sample if harvey_adj is True: harvey_adj = ((T+1-2*h+h*(h-1)/T)/T)**(0.5) DM_stat = harvey_adj*DM_stat # Calculate p-value p_value = 2*t.cdf(-abs(DM_stat), df=T - 1) dm_return = collections.namedtuple('dm_return', 'DM p_value') result = dm_return(DM=DM_stat, p_value=p_value) return result results = dm_test(oos_set, forecast_A, forecast_B, h=1, harvey_adj=True) print(results)
В нашем примере мы получаем следующий результат:
dm_return (DM = 2,090864056688617, p_value = 0,10472839911640372)
Поскольку значение p превышает 0,05 (учитывая уровень значимости 95%), мы принимаем нулевую гипотезу о равной прогностической точности. Другими словами, мы не можем сделать вывод, что один прогноз значительно точнее другого.
В заключение, тест Дибольда-Мариано является полезным инструментом для сравнения точности двух моделей прогнозирования. Это позволяет нам определить, какая модель лучше, если таковая имеется, и может помочь нам принимать более обоснованные решения при выборе метода прогнозирования.
Спасибо за прочтение!
Полный скрипт доступен в этом репозитории. Полная реализация теста Дибольда-Мариано и исходный код доступны здесь.
Рекомендации
Диболд, Ф.Х. и Мариано, Р. (1995), «Сравнение точности прогнозов», Journal of Business and Economic Statistics, 13, 253–265.
Фрэнсис X. Диболд (2015) Сравнение точности прогнозов, двадцать лет спустя: личный взгляд на использование и злоупотребление тестами Дибольда-Мариано, Журнал деловой и экономической статистики, 33: 1, 1, DOI: 10.1080/07350015.2014.983236