Авторы: Наби Рагимов, Александра Грищенко, Грейсон Майкл Калберт
Защищено авторским правом.
В этом примере я расскажу, как модели векторной авторегрессии (VAR) и машинного обучения (ML) работают при прогнозировании инфляции в трех странах Балтии на основе теории кривой Филлипса. Прежде чем строить модели, давайте обсудим модель VAR.
Основным инструментом, используемым для моделирования и прогнозирования, является векторная авторегрессионная (VAR) модель. Эта модель многомерного временного ряда работного дома связывает текущие наблюдения переменной с прошлыми наблюдениями за самой собой и прошлыми наблюдениями за другими переменными в системе. Модели VAR характеризуются своим порядком, который относится к количеству более ранних периодов времени, которые будет использовать модель. Например, VAR 5-го порядка будет моделировать цену каждого года как линейную комбинацию цен за последние пять лет. Лаг — это значение переменной в предыдущий период времени. Таким образом, в общем, VAR p-го порядка относится к модели VAR, которая включает лаги для последних p периодов времени. VAR p-го порядка обозначается как «VAR(p)» и иногда называется «VAR с p-лагами».
В качестве инструмента для построения всех моделей и прогнозов использовался Python. Основные библиотеки, которые мы использовали: pandas, numpy, statsmodels, scipy, matplotlib, seaborn, sklearn и xgboost.
Нами были выполнены следующие шаги:
- Импортируйте наборы данных
- Визуализируйте временной ряд
- Проверка причинности с помощью теста причинности Грейнджера
- Коинтеграционный тест
- Разделите серию на обучающие и тестовые данные
- Проверить стационарность и сделать временной ряд стационарным
- Выбор порядка () модели VAR
- Обучите модель VAR выбранного заказа (p)
- Проверка последовательной корреляции остатков (погрешностей) с использованием статистики Дарбина-Уотсона.
- Как прогнозировать модель VAR с помощью статистических моделей
- Обучите модель VAR выбранного заказа (p)
- Инвертируйте преобразование, чтобы получить реальный прогноз
- График прогноза против фактов
- Оцените прогнозы
- Результаты и выводы
эмпирические результаты
Наши данные содержат реальный ВВП, Гармонизированные индексированные потребительские цены (HICP) и уровень безработицы в трех странах Балтии. Мы проанализируем три страны Балтии и выберем одну из них, для которой будем прогнозировать инфляцию, применяя различные модели. Мы построили график изменения этих переменных для стран Балтии за определенный период времени. Обугливание показано на рисунке 1. Очевидно, что каждый из рядов имеет довольно схожий тренд по годам. Поскольку экономики Эстонии, Латвии и Литвы связаны друг с другом и взаимосвязаны, разумно наблюдать схожие закономерности в макроэкономических переменных.
Для статистической проверки связи между этими переменными мы использовали корреляционную матрицу. Матрица корреляции представлена на рисунке 3. Между тремя странами матрица показывает высокие коэффициенты корреляции по всем трем переменным. Однако это не свидетельствует об их влиянии друг на друга (это трудно определить по имеющимся данным), но свидетельствует об их сходных тенденциях. Во всех трех странах существует отрицательная корреляция между инфляцией и уровнем безработицы. Более того, уровень безработицы и уровень ВВП имеют довольно сильную отрицательную корреляцию в каждой из прибалтийских стран. И наоборот, мы можем видеть положительную корреляцию между ВВП и индексом инфляции.
Мы также создали точечный график, чтобы визуализировать взаимосвязь между инфляцией, ВВП и уровнем безработицы в Эстонии. На следующих рисунках изображены графики рассеяния. Несмотря на некоторые выбросы, существует изогнутая отрицательная связь между уровнем безработицы и инфляцией, что дает некоторые визуальные доказательства существования кривой Филлипса. Однако четкой связи между инфляцией и уровнем ВВП, которую можно было бы увидеть на графике, нет. Запустив приведенный ниже код, вы можете получить результат.
def scatter_plot(data, col_1, col_2, color):
"""
Creates scatter plots for inflation, unemployment rate, GDP
"""
fig, ax = plt.subplots(figsize=(7,4))
ax.scatter(x=data[col_1], y=data[col_2], alpha=0.4, color=color)
ax.set_xlabel("Inflation", fontname="Verdana")
if col_2 == "EE_Un":
ax.set_title("Inflation and unemployment rate EE", fontsize=16, fontname="Verdana", loc="left")
ax.set_ylabel("Unemployment rate", fontname="Verdana")
else:
ax.set_title("Inflation and GDP EE", fontsize=16, fontname="Verdana", loc="left")
ax.set_ylabel("GDP", fontname="Verdana")
for f,c in dict(zip(["EE_Un", "EE_GDP"], ["navy", "firebrick"])).items():
scatter_plot(data=baltic_data_ee, col_1="EE_HICP", col_2=f, color=c)
Проверка причинности с помощью теста причинности Грейнджера
В основе векторной авторегрессии лежит то, что каждый из временных рядов в системе влияет друг на друга. То есть мы можем предсказать ряд с прошлыми значениями самого себя вместе с другими рядами в системе. Используя тест причинности Грейнджера, можно проверить эту связь перед построением модели. Тест причинности Грейнджера проверяет нулевую гипотезу о том, что коэффициенты прошлых значений в уравнении регрессии равны нулю. Проще говоря, прошлые значения временного ряда (X) не вызывают другой ряд (Y). Итак, если значение p, полученное в результате теста, меньше уровня значимости 5%, то можно смело отвергать нулевую гипотезу. В таблице 1 мы видим, что инфляция и безработица в Эстонии по Грейнджеру обусловливают уровень ВВП в Эстонии с p-значениями 0 и 0,0367 соответственно. Хотя инфляция в Эстонии по Грейнджеру вызывает безработицу, обратное неверно, исходя из матрицы причинно-следственных связей Грейнджера.
# Testing Causation using Granger’s Causality Test
from statsmodels.tsa.stattools import grangercausalitytests
maxlag=12
test = 'ssr_chi2test'
def grangers_causation_matrix(data, variables, test='ssr_chi2test', verbose=False):
"""Check Granger Causality of all possible combinations of the Time series.
The rows are the response variable, columns are predictors. The values in the table
are the P-Values. P-Values lesser than the significance level (0.05), implies
the Null Hypothesis that the coefficients of the corresponding past values is
zero, that is, the X does not cause Y can be rejected.
data : pandas dataframe containing the time series variables
variables : list containing names of the time series variables.
"""
df = pd.DataFrame(np.zeros((len(variables), len(variables))), columns=variables, index=variables)
for c in df.columns:
for r in df.index:
test_result = grangercausalitytests(data[[r, c]], maxlag=maxlag, verbose=False)
p_values = [round(test_result[i+1][0][test][1],4) for i in range(maxlag)]
if verbose: print(f'Y = {r}, X = {c}, P Values = {p_values}')
min_p_value = np.min(p_values)
df.loc[r, c] = min_p_value
df.columns = [var + '_x' for var in variables]
df.index = [var + '_y' for var in variables]
return df
grangers_causation_matrix(baltic_data_ee, variables = baltic_data_ee.columns)
Коинтеграционный тест
В модели векторной авторегрессии (VAR) коинтеграция относится к существованию линейной комбинации переменных в модели. Другими словами, это означает, что переменные в модели находятся в долгосрочном равновесии и не являются независимыми друг от друга. Чтобы проверить коинтеграцию в модели VAR, мы можем использовать тест на коинтеграцию, такой как тест Энгла-Грейнджера или тест Йохансена. Эти тесты включают регрессию переменных в модели друг относительно друга и проверку того, являются ли остатки этой регрессии стационарными. Если остатки стационарны, это предполагает, что переменные коинтегрированы. В тесте на коинтеграцию для векторной авторегрессионной (VAR) модели нулевая гипотеза состоит в том, что переменные в модели не коинтегрированы. Другими словами, нулевая гипотеза состоит в том, что переменные в модели не имеют долгосрочных равновесных отношений и не зависят друг от друга. Основываясь на данных Эстонии, в следующей таблице 2 мы видим, что инфляция и безработица не коинтегрированы с ВВП и не имеют долгосрочной связи.
Тестирование на стационарность
Поскольку модель VAR требует, чтобы временные ряды, которые вы хотите спрогнозировать, были стационарными, принято проверять все временные ряды в системе на стационарность. Стационарный временной ряд — это ряд, характеристики которого, такие как среднее значение и дисперсия, не меняются с течением времени. Существует набор тестов, называемых тестами единичного корня. Давайте воспользуемся тестом ADF для нашей цели. Нулевая гипотеза об остатке единичного корня подразумевает, что единичный корень существует и временной ряд не является стационарным. Когда p-значение больше 5%, мы можем отклонить нулевую гипотезу и подтвердить, что временной ряд является стационарным. До первого дифференцирования ВВП и уровень инфляции в Эстонии не являются стационарными, за исключением безработицы. После первой разности ВВП в Эстонии все еще был нестационарным. Все ряды в модели VAR должны иметь одинаковое количество наблюдений. Вот почему мы взяли вторую разность, чтобы сделать временной ряд ВВП также стационарным.
# Check for Stationarity and Make the Time Series Stationary
def adfuller_test(series, signif=0.05, name='', verbose=False):
"""Perform ADFuller to test for Stationarity of given series and print report"""
r = adfuller(series, autolag='AIC')
output = {'test_statistic':round(r[0], 4), 'pvalue':round(r[1], 4), 'n_lags':round(r[2], 4), 'n_obs':r[3]}
p_value = output['pvalue']
def adjust(val, length= 6): return str(val).ljust(length)
# Print Summary
print(f' Augmented Dickey-Fuller Test on "{name}"', "\n ", '-'*47)
print(f' Null Hypothesis: Data has unit root. Non-Stationary.')
print(f' Significance Level = {signif}')
print(f' Test Statistic = {output["test_statistic"]}')
print(f' No. Lags Chosen = {output["n_lags"]}')
for key,val in r[4].items():
print(f' Critical value {adjust(key)} = {round(val, 3)}')
if p_value <= signif:
print(f" => P-Value = {p_value}. Rejecting Null Hypothesis.")
print(f" => Series is Stationary.")
else:
print(f" => P-Value = {p_value}. Weak evidence to reject the Null Hypothesis.")
print(f" => Series is Non-Stationary.")
for name, column in data_train.iteritems():
adfuller_test(column, name=column.name)
print('\n')
На рисунке 9 мы видим сравнение исторических тенденций до и после преобразования. После преобразования ряд становится стационарным, а это означает, что среднее значение и дисперсия не меняются со временем.
Выбор порядка (p) модели VAR
Чтобы выбрать правильный порядок модели VAR, мы итеративно подбираем возрастающие порядки модели VAR и выбираем порядок, который дает модель с наименьшим AIC. Хотя обычной практикой является рассмотрение AIC, вы также можете проверить другие наиболее подходящие сравнительные оценки BIC, FPE и HQIC. На рисунке 10 мы описали изменение значения AIC через количество лагов. Из графика видно, что AIC продолжает снижаться до 12-го лага и начинает расти на 13-м лаге. Итак, для автоматического выбора лагового порядка VAR со встроенной функцией мы будем использовать 12-лаговый диапазон для нашей модели VAR.
# Select the Order (P) of VAR model
forecasting_model = VAR(data_diff)
results_aic = []
for p in range(1,15):
results = forecasting_model.fit(p)
results_aic.append(results.aic)
np.argmin(results_aic)+1
sns.set()
plt.plot(list(np.arange(1,15,1)), results_aic)
plt.xlabel("Order")
plt.ylabel("AIC")
plt.show()
Тест Дарбина-Ватсона для последовательной корреляции остатков
В статистике серийная корреляция (также известная как автокорреляция) относится к корреляции между значениями временного ряда в разные моменты времени. Когда остатки (наблюдаемые значения минус подогнанные значения из статистической модели) временного ряда демонстрируют последовательную корреляцию, это означает, что остатки в один момент времени коррелируют с остатками в предыдущие моменты времени. Это может указывать на то, что модель указана неправильно или что есть пропущенные переменные, влияющие на модель. Мы будем использовать тест Дарбина-Уотсона, нулевая гипотеза теста состоит в том, что в остатках нет последовательной корреляции. Статистика теста Дарбина-Ватсона представляет собой значение от 0 до 4, при этом значения ближе к 2 указывают на отсутствие автокорреляции. Тест сравнивает наблюдаемое значение статистики Дарбина-Ватсона с критическим значением, и если наблюдаемое значение выходит за пределы критического значения, нулевая гипотеза об отсутствии автокорреляции отклоняется. Для наших переменных последовательная корреляция выглядит хорошо, и у нас нет проблемы сезонности: [EE_GDP: 1,97, EE_HICP: 2,27, EE_Un: 1,82].
# Durbin-Watson test for serial correlation in residuals
def adjust(val, length= 6):
return str(val).ljust(length)
from statsmodels.stats.stattools import durbin_watson
out = durbin_watson(model_fitted.resid)
for col, val in zip(baltic_data_ee.columns, out):
"""
Durbin-Watson:The null hypothesis of the test is that there is no serial correlation in the residuals.
The Durbin-Watson test statistic is a value between 0 and 4, with values closer to 2 indicating the absence
of autocorrelation. Regarding our variables, test statistic values are closer to 2, which implies that there is
no serial correlation among residuals.
"""
print(adjust(col), ':', round(val, 2))
После приведения ряда к стационарному уровню мы повторно использовали критерий причинности и коинтеграции Грейнджера. Тест причинности Грейнджера дал несколько лучшие результаты. Уровень безработицы не влияет на ВВП по Грейнджеру по сравнению с предыдущим результатом. Хотя безработица по-прежнему не вызывает инфляцию по Грейнджеру, p-значение снизилось до 0,09, что является значимым на уровне 10%. Коинтеграционный тест показывает, что между всеми временными рядами существует долгосрочная связь.
Прогнозирование
Для прогнозирования все имеющиеся наблюдения были разделены на обучающие и тестовые данные. Период тестовых данных: с 01.04.2020 по 01.07.2022. Как видно далее, из-за большого количества непредсказуемых изменений за указанный период, таких как пандемия COVID-19, энергетический кризис и война в Украине, модель не способна сделать достаточно точный прогноз. Для сравнения мы использовали три метода прогнозирования и три модели:
- VAR: традиционный подход
- VAR: метод лассо
- VAR: подход XGB
Как упоминалось выше, для выбора p-порядка использовалась автоматическая функция, которая сравнивает значения AIC.
Оценка прогнозов
В этом разделе мы оценим эффективность традиционных подходов VAR и машинного обучения (ML) для прогнозирования инфляции в Эстонии.
Мы построили фактические и подобранные значения инфляции после прогнозирования и оценили эффективность прогноза. Для оценки прогнозов мы рассчитали комплексный набор показателей, а именно MAPE, ME, MAE, MPE и RMSE. Во-первых, на рисунке 11 мы сначала предсказали инфляцию в Эстонии с помощью нашей модели VAR. Вы можете видеть, что визуализации показывают, что VAR на самом деле не работает хорошо в течение прогнозируемого периода. Во-вторых, на Рисунке 12 мы применили регрессию Лассо, чтобы определить, лучше ли характеристики, выбранные с помощью VAR, предсказывают инфляцию в Эстонии. Можно видеть, что Лассо немного лучше, чем традиционная модель VAR, предсказывает инфляцию на основе исторических данных.
start_date = "2020-07-01"
end_date = "2022-07-01"
def var_create(columns, data):
"""
Creates vector autoregressive model given data and list of selected variables
Returns the MSE between forecasted and actual values
Also returns the concatenated dataset for visualization purposes
"""
data = data[columns]
data = data.dropna(axis=0)
data.index.to_period("M")
# Split dataset and run VAR on the trained part
data_train = data.loc["2000-01":"2020-04", :]
var_train = VAR(data_train)
results = var_train.fit(12)
lag_order = results.k_ar
forecasted = pd.DataFrame(results.forecast(data_train.values[-lag_order:], 24)) # Forecast 24 months
# Rename forecasted columns
forecasted_names = list(forecasted.columns.values)
data_train_names = list(data_train.columns.values)
var_dict = dict(zip(forecasted_names, data_train_names))
for f,t in var_dict.items():
forecasted = forecasted.rename(columns={f:t + "_fcast"})
# forecasted.index = pd.DatetimeIndex(start=start_date, periods=forecasted.shape[0], freq="MS")
forecasted.index = pd.date_range(start=start_date, periods=forecasted.shape[0], freq="MS")
forecasted.index.names = ["date"]
# Parse together forecasted data with original dataset
final_data = pd.merge(forecasted, data, left_index=True, right_index=True)
final_data = final_data.sort_index(axis=0, ascending=True)
final_data = pd.concat([data_train, final_data], sort=True, axis=0)
final_data = final_data.sort_index(axis=0, ascending=True)
var_mse = metrics.mean_squared_error(final_data.loc[start_date:end_date,"EE_HICP_diff_fcast"],
final_data.loc[start_date:end_date,"EE_HICP_diff"])
return var_mse, final_data
def plot_cpi(final_data, var_mse, approach):
"""
Plots the actual values against forecast
"""
fig, ax = plt.subplots(figsize=(14,6))
colors = sns.color_palette("deep", 8)
final_data["EE_HICP_diff_fcast"].plot(ax=ax, legend=True, linewidth=2.5, linestyle="dashed")
final_data["EE_HICP_diff"].plot(ax=ax, legend=True, alpha=0.6, linestyle="solid")
if approach=="traditional":
ax.set_title("VAR in-sample forecast, traditional approach", fontsize=16,
fontweight="bold", fontname="Verdana", loc="left")
elif approach=="lasso":
ax.set_title("VAR in-sample forecast, Lasso approach", fontsize=16,
fontweight="bold", fontname="Verdana", loc="left")
elif approach=="XGBoost":
ax.set_title("VAR in-sample forecast, XGBoost approach", fontsize=16,
fontweight="bold", fontname="Verdana", loc="left")
else:
ax.set_title("VAR in-sample forecast", fontsize=16,
fontweight="bold", fontname="Verdana", loc="left")
ax.set_ylabel("First differences", fontname="Verdana")
ax.legend([f"VAR Forecast, MSE={var_mse}", "EE_HICP Real Values"])
plot_cpi(final_data=df1, var_mse=mse1, approach="traditional")
plot_cpi(final_data=df2, var_mse=mse2, approach="lasso")
Следует понимать, что в данных временных рядов Лассо выполняет k-кратность на скользящей основе, поскольку данные складываются по наблюдениям во времени. На рисунке 13 мы видим сравнение фактических и подогнанных значений на основе Лассо. Фактические и подогнанные значения не перекрываются.
Наконец, на рисунке 14 мы использовали XGBoost для выбора функций для VAR. Оказывается, XGBoost и другие древовидные модели лучше обрабатывают тренды по сравнению с VAR и Lasso в наборах данных временных рядов. Мы обучили и протестировали XGBoost на наборе данных об инфляции в Эстонии. XGBoost немного превосходит Lasso и традиционный подход VAR. На рисунке 15 мы построили график, на котором оценка производительности трех моделей была описана на основе среднеквадратичной ошибки (MSE).
В целом, мы можем сказать, что все модели работают более или менее одинаково, но XGBoost делает прогноз немного лучше.
Заключение
Можно сделать вывод, что ни одна из этих моделей не смогла предсказать инфляцию с высокой точностью. Это показывает неэффективность предлагаемых моделей для стран Балтии, особенно для Эстонии. Трудно сказать, является ли плохой результат специфичным для страны и региона или нет, поскольку он требует дальнейшего анализа.
