Макроэкономическое прогнозирование

В этом руководстве описывается, как преобразовать данные, оценить простой базовый метод и стандартную многомерную регрессию временных рядов, а затем сравнить с регрессией гауссовского процесса (GP).

Макроэкономическое прогнозирование с гауссовскими процессами

Точное макроэкономическое прогнозирование необходимо для принятия деловых решений, финансовой торговли и в качестве основы для политических решений политиков и экономистов. Традиционно это делалось с помощью линейной регрессии данных временных рядов или тесно связанных методов. В этом руководстве мы попытаемся спрогнозировать ВВП для ряда разных стран, используя регрессию гауссовского процесса, популярный нелинейный метод, и посмотрим, получим ли мы лучшие результаты, чем традиционные методы.

Наиболее важной макроэкономической переменной, вероятно, является валовой внутренний продукт (ВВП), поскольку он измеряет стоимость всех товаров и услуг, проданных в регионе. Поэтому остановимся на прогнозировании ВВП.

Учебник организован следующим образом. Сначала мы загружаем и преобразуем некоторые данные, затем оцениваем простой базовый метод и стандартную многомерную регрессию временных рядов и, наконец, сравниваем с регрессией гауссовского процесса (GP).

Окружающая обстановка

В этом руководстве мы используем Python 3.5 (хотя код также работает с Python 2.7). Мы будем использовать пакеты pandas, statsmodels и scikit-learn, которые вам нужно будет установить на свой компьютер. Чтобы установить эти пакеты, вы можете ввести в свой терминал / консоль следующее…

    > pip3 install pandas
    > pip3 install statsmodels
    > pip3 install sklearn

Использование диспетчера пакетов Python pip3 (> - это приглашение терминала). Если вы используете дистрибутив Python anaconda, эти пакеты включены по умолчанию.

Конфигурация

Сначала мы указываем некоторую конфигурацию. Мы можем изменить конфигурацию, если хотим запускать код для разных стран, индикаторов или периодов времени.

Ввод:

# Forecasting years, number of years to use for out-of-sample 
# evaluation (We will create one forecast per country for each year
# which we compare to the actual GDP)
nyears = 10
# Number of lags to use for GP regression
lags = 5
# Indicator labels and names in the World Bank API
indicators  = {"gdp"        : "NY.GDP.MKTP.CD",
               "population" :    "SP.POP.TOTL",
               "inflation"  : "FP.CPI.TOTL.ZG"}
nindicators = len(indicators)
# The variable to forecast, should be one of the indicator labels
target_variable = "gdp"
# Countries to include in the data, specified as ISO country codes
countries  = ['au','ca','de','es','fr','gb','jp','us']
ncountries = len(countries)
# Start and end year for the data set
start_year = 1976
end_year   = 2016

Данные

Загрузка

Сначала мы загружаем некоторые макроэкономические показатели из HTTP API Всемирного банка с помощью удобного пакета Python requests, входящего в вашу установку Python.

Мы указываем URL-адрес шаблона, который мы заполняем разными странами, показателями и датами, которые мы указали в конфигурации выше. Когда мы запрашиваем URL-адрес, HTTP API ответит выбранными данными. Данные извлекаются в формате json, который мы перебираем и помещаем в pandas фрейм данных.

Ввод:

import requests
import pandas as pd
template_url = "http://api.worldbank.org/v2/countries/{0}/indi"
template_url +="cators/{1}?date={2}:{3}&format=json&per_page=999"
# Countries should be ISO identifiers separated by semi-colon
country_str = ';'.join(countries)
raw_data = pd.DataFrame()
for label, indicator in indicators.items():
# Fill in the template URL
    url = template_url.format(country_str, indicator, 
                                  start_year, end_year)
    
    # Request the data
    json_data = requests.get(url)
    
    # Convert the JSON string to a Python object
    json_data = json_data.json()
    
    # The result is a list where the first element is meta-data, 
    # and the second element is the actual data
    json_data = json_data[1]
    
    # Loop over all data points, pick out the values and append 
    # them to the data frame
    for data_point in json_data:
        
        country = data_point['country']['id']
        
        # Create a variable for each country and indicator pair
        item    = country + '_' + label
        
        year    = data_point['date']
        
        value   = data_point['value']
        
        # Append to data frame
        new_row  = pd.DataFrame([[item, year, value]],
                                columns=['item', 'year', 'value'])
        raw_data = raw_data.append(new_row)
# Pivot the data to get unique years along the columns,
# and variables along the rows
raw_data = raw_data.pivot('year', 'item', 'value')
# Let's look at the first few rows and columns
print('\n', raw_data.iloc[:10, :5], '\n')

Вывод:

Исследование

Сначала давайте исследуем данные с помощью графиков.

Ввод:

import matplotlib.pyplot as plt
for lab in indicators.keys():
    
    indicator = raw_data[[x for x in raw_data.columns 
                              if x.split("_")[-1] == lab]]
    indicator.plot(title=lab)
    plt.show()

Вывод:

Более значительный рост ВВП в США отчасти можно объяснить более быстрым ростом населения. Это должно быть зафиксировано нашими моделями.

Трансформация

Прежде чем использовать данные в наших моделях, мы преобразуем их в более подходящий формат.

Ввод:

import numpy as np
# Calculate rates of change instead of absolute levels
# (Runtime warning expected due to NaN)
data = np.log(raw_data).diff().iloc[1:,:]
# Set NaN to zero
data.fillna(0, inplace=True)
# Subtract the mean from each series
data = data - data.mean()
# Convert to date type
data.index = pd.to_datetime(data.index, format='%Y')
# Put the target variable into a separate data frame
target = data[[x for x in data.columns 
                   if x.split("_")[-1] == target_variable]]

Оценка

Все модели оцениваются одинаково. Для каждого года прогнозирования (данные за последние годы) мы обучаем модель, используя предыдущие годы, и создаем прогноз для каждой страны. Для каждого прогноза мы вычисляем квадрат ошибки по сравнению с истинным значением, а затем вычисляем среднее значение по странам и годам, прежде чем извлекать квадратный корень для получения среднеквадратичной ошибки (RMSE). Другими словами

где ŷ𝒾 - прогноз, 𝘺𝒾 - фактическое значение целевой переменной, а mm - количество создаваемых нами прогнозов (количество лет прогнозирования, умноженное на количество стран).

Базовые прогнозы

Сначала мы установим производительность простого базового метода, который будет использоваться для сравнения с последующими методами. Наш базовый метод просто предсказывает ту же скорость изменения, что и предыдущее наблюдение целевой переменной.

Только в том случае, если более сложный метод значительно превосходит наивный базовый прогноз, ему следует отдавать предпочтение в соответствии с принципом бритвы Оккама и принципом KISS в разработке программного обеспечения.

Ввод:

errors = target.iloc[-nyears:] - target.shift().iloc[-nyears:]
# Root mean squared error
rmse = errors.pow(2).sum().sum()/(nyears*ncountries)**.5
print('\n\t' + '-' * 18)
print("\t| Error: ", np.round(rmse, 4), '|')
print('\t' + '-' * 18 + '\n')

Вывод:

Многомерный временной ряд

Мы исследуем эффективность многомерной регрессии временных рядов, также называемой векторной авторегрессией (VAR) в эконометрике. Модель принимает вид

где 𝑑 - количество независимых переменных, 𝜀 - количество лагов, εtεt - член ошибки, а βi, jβi, j - параметры, которые необходимо оценить.

Чтобы узнать больше об анализе временных рядов, хорошей вводной книгой является «Введение в временные ряды и прогнозирование» Броквелла и Дэвиса, а еще две более сложные книги - «Временные ряды: теория и методы», также написанные Броквеллом и Дэвисом, и «Время. Анализ рядов »Гамильтона.

Чтобы соответствовать модели, мы используем пакет statsmodels, VAR реализация которого использует метод наименьших квадратов для оценки параметров. Мы обучаем модель с использованием одного лага, то есть регрессии целевой переменной по последним предыдущим значениям независимых переменных, но автоматические методы выбора количества лагов также поддерживаются в пакете statsmodels, например, информационный критерий Акаике (AIC). или байесовский информационный критерий (BIC).

Ввод:

from statsmodels.tsa.api import VAR
# Sum of squared errors
sse = 0
for t in range(nyears):
    
    # Create a VAR model
    model = VAR(target.iloc[t:-nyears+t], freq='AS')
    
    # Estimate the model parameters
    results = model.fit(maxlags=1)
    
    actual_values = target.values[-nyears+t+1]
    
    forecasts = results.forecast(target.values[:-nyears+t], 1)
    forecasts = forecasts[0,:ncountries]
sse += ((actual_values - forecasts)**2).sum()
# Root mean squared error
rmse = (sse / (nyears * ncountries))**.5
print('\n\t' + '-' * 18)
print("\t| Error: ", np.round(rmse, 4), '|')
print('\t' + '-' * 18 + '\n')

Вывод:

Мы видим, что ошибка примерно такая же, как и в базовой модели, которая выступает против использования более сложной модели VAR для макроэкономического прогнозирования с текущими данными.

Регрессия гауссовского процесса

Теперь мы попытаемся спрогнозировать ВВП, используя регрессию по Гауссу. Наша модель похожа на предыдущую и дается формулой

где 𝘺𝑡 - это наблюдения целевой переменной, каждый - вектор наблюдений независимой (предикторной) переменной, 𝑓 (𝓍) - функция регрессии, а шумовой член ε𝑡 имеет нормальное (гауссово) распределение с дисперсией α. Для каждого прогноза вне выборки мы поместим nn предыдущих наблюдений целевой переменной 𝘺 в вектор-столбец, который мы обозначим 𝘺. Мы также будем обозначать 𝑋 матрицу nn предыдущих наблюдений переменной-предиктора, где каждое наблюдение занимает строку в матрице. Здесь каждое наблюдение будет включать запаздывающие значения наших переменных по сравнению с предыдущими (страны и индикаторы).

Вместо функции регрессии 𝑓 (𝓍), являющейся линейной функцией с неизвестными параметрами, функция теперь является случайной, так что ее значение в каждой точке 𝓍𝑡 является нормально распределенной случайной величиной, и такой, что ковариация между двумя случайными величинами 𝑓 (𝓍𝑡 ) и 𝑓 (𝓍𝑠) задается так называемой ковариационной функцией (или ядром) 𝑘 (𝓍𝑡, 𝓍𝑠), которую необходимо указать. Функция 𝑓 (𝓍) называется гауссовским процессом. Поскольку член ошибки имеет нормальное распределение, целевые значения также будут нормально распределены с ковариацией, равной (𝓍𝑡, 𝓍𝑠) + 𝑎, если 𝑡 = 𝑠, и 𝑘 (𝓍𝑡, 𝓍𝑠) в противном случае.

Предположим, мы хотим оценить модель, используя nn точек данных для 𝘺𝑡 и 𝓍𝑡, и что у нас есть новая контрольная точка 𝓍𝑛 + ₁, на основе которой мы хотим сделать прогноз ŷ 𝑛 + ₁. Условное распределение f (𝓍𝑛 + ₁) при заданных 𝘺, 𝑋 и 𝓍𝑛 + ₁ тогда также нормально распределено, и мы можем принять наш прогноз ŷ𝑛 + ₁ как условное ожидаемое значение 𝑓 (𝓍𝑛 + ₁). По известной формуле многомерного гауссовского распределения оно дается выражением

где K - матрица размером 𝑛 ×, где запись в строке 𝑡 и столбце задается как 𝑘 (𝓍𝑠, 𝓍𝑡), 𝐼 - это единичная матрица, а 𝐤 - вектор-столбец длины 𝑛 с элементом 𝑡, заданным как 𝑘 (𝓍𝑡 , 𝓍𝑛 + ₁).

Чтобы узнать больше о методах гауссовских процессов, которые также можно использовать для классификации, см., Например, книгу Расмуссена и Вильямса «Гауссовские процессы в машинном обучении».

Для реализации гауссовской регрессии процесса мы используем готовый функционал в пакете scikit-learn (см. Эту ссылку для документации). Для каждого года прогнозирования мы создаем вектор 𝘺 прошлых наблюдений целевой переменной и матрицу данных 𝑋 наблюдений независимой переменной. Мы вводим их в класс GaussianProcessRegressor для оценки модели, а затем делаем прогноз на следующий год. Как и раньше, мы вычисляем ошибку ряда прогнозов, которую суммируем в среднеквадратичную ошибку.

Мы будем использовать общую гауссовскую ковариационную функцию, которая задается формулой

где ∥v∥ - длина вектора vv, а σσ - параметр.

Ввод:

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
# We set the parameter of the covariance function to
# approximately equal to the median distance between data points,
# a common heuristic for this covariance function. The 'alpha' 
# argument is the noise variance, which we set equal to the 
# covariance parameter.
gpr = GaussianProcessRegressor(kernel=RBF(0.1), alpha=0.1)
# Number of data points for estimation/fitting for each forecast
ndata = target.shape[0] - nyears - lags
# Sum of squared errors
sse = 0
for t in range(nyears):
    
    # Observations for the target variables
    y = np.zeros((ndata, ncountries))
# Observations for the independent variables
    X = np.zeros((ndata, lags*ncountries*nindicators))
    
    for i in range(ndata):
        
        y[i] = target.iloc[t+i+1]
        X[i] = data.iloc[t+i+2:t+i+2+lags].values.flatten()
        
    gpr.fit(X, y)
    
    x_test   = np.expand_dims(data.iloc[t+1:t+1+lags].values.flatten(), 0)
    forecast = gpr.predict(x_test)
    
    sse += ((target.iloc[t].values - forecast)**2).sum()
    
rmse = (sse / (nyears * ncountries))**.5
print('\n\t' + '-' * 18)
print("\t| Error: ", np.round(rmse, 4), '|')
print('\t' + '-' * 18 + '\n')

Вывод:

С помощью гауссовской регрессии процесса мы получаем лучшие на данный момент результаты с небольшим улучшением базовой и линейной регрессионных моделей.

об авторе

Фредрик - аспирант 4-го курса по статистике в UCL, работает в области машинного обучения. Перед тем как получить докторскую степень, он несколько лет работал в хедж-фонде, а до этого изучал физику и экономику. Его исследовательские интересы в основном связаны с ядерными методами, но он также интересуется нейронными сетями, числовой линейной алгеброй и статистическим оцениванием. Одно из его хобби - сбор подписей за книгами, среди которых в настоящее время есть два нобелевских лауреата и медалист Филдса.

Связаться с Фредриком

Спасибо за внимание! Если вам понравился пост, мы будем благодарны за вашу поддержку, нажав кнопку хлопка (👏🏼) ниже или поделившись этой статьей, чтобы ее могли найти другие.

Если вы хотите первым узнавать о нашем новом содержании, включая новые руководства, тематические исследования и другие полезные ресурсы, нажмите здесь!

Дополнительные уроки для отработки навыков:

- Развертывание модели машинного обучения в сети
- Шесть проектов Data Science для расширения ваших навыков и знаний
- Модульное тестирование с PySpark
- Настройка гиперпараметров в XGBoost
- Начало работы с XGBoost