Kaggle: не можем ходить — будем бегать / forpes.ru

Главная
Kaggle: не можем ходить — будем бегать

Kaggle: не можем ходить — будем бегать +27

06.03.2019 05:22

niklyaev 13 6800 Источник

Насколько сложна тема машинного обучения? Если Вы неплохо математически подкованы, но объем знаний о машинном обучении стремится к нулю, как далеко Вы сможете зайти в серьезном конкурсе на платформе Kaggle?

О площадке и конкурсе

Kaggle — это сообщество интересующихся ML (от новичков до крутых профи) и площадка для проведения конкурсов (часто с внушительным призовым фондом).

Чтобы сразу погрузиться во все прелести ML, решил сразу выбрать серьезный конкурс. Такой как раз был в наличии: Two Sigma: Using News to Predict Stock Movements. Суть конкурса в двух словах — предсказывать цену акций различных компаний на основе состояния актива и новостей, связанных с данным активом. Призовой фонд конкурса составляет 100 000 $, которые распределят между участниками, занявшими первые 7 мест.

Конкурс является особенным по двум причинам:

это Kernels-only конкурс: обучать модели можно только в облаке Kaggle Kernels;
окончательное распределение мест будет известно только через полгода после завершения приема решений; в течение этого времени решения будут предсказывать цены на текущую дату.

О задаче

По условию мы должны предсказывать уверенность

$\hat{y}_{ti} \in [-1,1]$ в том, что доходность актива вырастет. Доходность актива считается относительно доходности рынка в целом. Целевая метрика является кастомной — это не более привычные RMSE или MAE, а коэффициент Шарпа, который в данном случае считается так:

$\text{score} = \frac{\bar{x}_t}{\sigma(x_t)},$

где

$x_t = \sum_i \hat{y}_{ti} r_{ti} u_{ti}$ ,

$r_{ti}$ — доходность актива i относительно рынка для дня t на 10-дневном горизонте,

$u_{ti}$ — булева переменная, показывающая, включается ли i-й актив в оценку для дня t,

$\bar{x}_t$ — среднее значение

$inline$ ,

$\sigma(x_t)$ — среднеквадратичное отклонение

$inline$ .

Коэффициент Шарпа — это доходность с поправкой на риск, значения коэффициента показывают эффективность трейдера:

меньше 1: плохая эффективность,
1 — 2: средняя, нормальная эффективность,
2 — 3: отличная эффективность,
больше 3: идеальная.

Данные по движению рынка

time(datetime64[ns, UTC]) — текущее время (в данных по движению рынка во всех строках 22:00 UTC)
assetCode(object) — идентификатор актива
assetName(category) — идентификатор группы активов для связи с данными новостей
universe(float64) — булево значение, показывающее, будет ли учитываться этот актив в расчете score
volume(float64) — объем торгов за день
close(float64) — цена закрытия для данного дня
open(float64) — цена открытия для данного дня
returnsClosePrevRaw1(float64) — доходность от закрытия до закрытия за предыдущий день
returnsOpenPrevRaw1(float64) — доходность от открытия до открытия за предыдущий день
returnsClosePrevMktres1(float64) — доходность от закрытия до закрытия за предыдущий день, скорректированная относительно движения рынка в целом
returnsOpenPrevMktres1(float64) — доходность от открытия до открытия за предыдущий день, скорректированная относительно движения рынка в целом
returnsClosePrevRaw10(float64) — доходность от закрытия до закрытия за предыдущие 10 дней
returnsOpenPrevRaw10(float64) — доходность от открытия до открытия за предыдущие 10 дней
returnsClosePrevMktres10(float64) — доходность от закрытия до закрытия за предыдущие 10 дней, скорректированная относительно движения рынка в целом
returnsOpenPrevMktres10(float64) — доходность от открытия до открытия за предыдущие 10 дней, скорректированная относительно движения рынка в целом
returnsOpenNextMktres10(float64) — доходность от открытия до открытия за следующие 10 дней, скорректированная относительно движения рынка в целом. Это значение мы и будем предсказывать.

Новостные данные

time(datetime64[ns, UTC]) — время в UTC доступности данных
sourceTimestamp(datetime64[ns, UTC]) — время в UTC публикации новости
firstCreated(datetime64[ns, UTC]) — время в UTC первой версии данных
sourceId(object) — идентификатор записи
headline(object) — заголовок
urgency(int8) — типы новостей (1: alert, 3: article)
takeSequence(int16) — не совсем понятный параметр, номер в некоторой последовательности
provider(category) — идентификатор провайдера новости
subjects(category) — список кодов тем новостей (может быть географический признак, событие, сектор промышленности и т.д.)
audiences(category) — список кодов аудиторий новости
bodySize(int32) — кол-во символов в теле новости
companyCount(int8) — кол-во компаний, явно упомянутых в новости
headlineTag(object) — некий тэг заголовка от компании Thomson Reuters
marketCommentary(bool) — признак, что новость касается общих условий рынка
sentenceCount(int16) — кол-во предложений в новости
wordCount(int32) — кол-во слов и знаков препинания в новости
assetCodes(category) — список активов, упомянутых в новости
assetName(category) — код группы активов
firstMentionSentence(int16) — предложение, в котором впервые упоминается актив:
relevance(float32) — число от 0 до 1, показывающее релевантность новости относительно актива
sentimentClass(int8) — класс тональности новости
sentimentNegative(float32) — вероятность, что тональность отрицательная
sentimentNeutral(float32) — вероятность, что тональность нейтральная
sentimentPositive(float32) — вероятность, что тональность положительная
sentimentWordCount(int32) — кол-во слов в тексте, которые относятся к активу
noveltyCount12H(int16) — «новизна» новости за 12 часов, рассчитывается относительно предыдущих новостей о данном активе
noveltyCount24H(int16) — то же, за 24 часа
noveltyCount3D(int16) — то же, за 3 дня
noveltyCount5D(int16) — то же, за 5 дней
noveltyCount7D(int16) — то же, за 7 дней
volumeCounts12H(int16) — объем новостей о данном активе за 12 часов
volumeCounts24H(int16) — то же, за 24 часа
volumeCounts3D(int16) — то же, за 3 дня
volumeCounts5D(int16) — то же, за 5 дней
volumeCounts7D(int16) — то же, за 7 дней

Задача по сути является задачей бинарной классификации, то есть, мы предсказываем бинарный признак, будет доходность расти (1 класс) или уменьшаться (0 класс).

Об инструментах

Kaggle Kernels — облачная вычислительная платформа, поддерживающая совместную работу. Поддерживаются kernel’ы следующих типов:

Python Script
R Script
Jupyter Notebook
RMarkdown

Каждый kernel запускается в своем docker-контейнере. В контейнере установлено большое количество пакетов, список для питона можно посмотреть здесь. Технические спецификации такие:

CPU: 4 ядра,
RAM: 17 ГБ,
диск: 5 ГБ постоянных и 16 ГБ временных,
максимальное время работы скрипта: 9 часов (на момент старта конкурса было 6 часов).

Также в Kernels доступны и GPU, однако в данном конкурсе использование GPU было запрещено.

Keras — высокоуровневый нейросетевой фреймворк, работающий поверх TensorFlow, CNTK или Theano. Очень удобный и понятный API, и есть возможность добавлять свои топологии сетей, функции потерь и всякое другое, используя API бэкенда.

Scikit-learn — крупная библиотека алгоритмов машинного обучения. Полезный источник алгоритмов предобработки и анализа данных для использования вместе с более специализированными фреймворками.

Валидация модели

Прежде чем отправлять модель на оценку, нужно как-то проверить локально, насколько хорошо она работает, — то есть, придумать способ локальной валидации. Я попробовал следующие подходы:

перекрестная валидация vs простое пропорциональное разбиение на тренировочное/тестовое множества;
локальный расчет коэффициента Шарпа vs ROC AUC.

В итоге наиболее близкие к конкурсной оценке результаты, как ни странно, показало сочетание пропорционального разбиения (эмпирически подобрал разбиение 0.85/0.15) и AUC. Вероятно, перекрестная валидация не очень подходит, так как поведение рынка сильно отличается на ранних этапах тренировочных данных и на оценочном периоде. Почему AUC сработал лучше, чем коэффициент Шарпа — вообще не берусь сказать.

Первые попытки

Так как задача заключается в предсказании временного ряда, то первым было опробовано классическое решение — рекуррентная нейросеть (RNN), а точнее, её варианты LSTM и GRU.

Главный принцип работы рекуррентных сетей заключается в том, что для каждого выходного значения на вход подается не один семпл, а целая последовательность. Из этого следует, что:

нам нужна некоторая предобработка исходных данных — генерация этих самых последовательностей длиной t дней для каждого актива;
модель на основе рекуррентной сети не сможет предсказать выходное значение, если нет данных за предыдущие t дней.

Я генерировал последовательности для каждого дня, начиная с t, поэтому при довольно больших t (от 20) полный набор тренировочных семплов перестал помещаться в память. Проблема решилась использованием генераторов, так как Keras умеет использовать генераторы в качестве входных и выходных наборов данных как для тренировки, так и для предсказания.

Первоначальная подготовка данных была максимально наивной: рыночные данные берем целиком плюс добавляем пару фич (день недели, месяц, номер недели в году), а данные новостей пока не трогаем вообще.

Первая модель использовала t = 10 и выглядела так:

model = Sequential()
model.add(LSTM(256, activation=act.tanh, return_sequences=True, input_shape=(data.timesteps, data.features)))
model.add(LSTM(256, activation=act.relu))
model.add(Dense(data.assets, activation=act.relu))
model.add(Dense(data.assets))

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

Temporal Convolutional Networks

Более современным нейросетевым решением для предсказания временных рядов являются TCN. Суть этой топологии весьма проста: берем одномерную сверточную сеть и применяем её к нашей последовательности длиной t. Более продвинутые варианты используют несколько сверточных слоев с разным dilation. Реализация TCN была частично скопирована (иногда на уровне идей) отсюда (визуализация стека TCN взята из статьи о Wavenet).

Первым относительно удачным решением стала вот такая модель, включающая слой GRU поверх TCN:

model = Sequential()
model.add(Conv1D(512,3, activation=act.relu, padding='causal', input_shape=(data.timesteps, data.features)))
model.add(Conv1D(100,3, activation=act.relu, padding='causal', dilation_rate=2))
model.add(Conv1D(100,3, activation=act.relu, padding='causal', dilation_rate=4))
model.add(GRU(256))
model.add(Dense(data.assets, activation=act.relu))

Такая модель выдает score = 0.27668. Небольшим тюнингом (количество фильтров TCN, размер batch) и увеличением t до 100 получаем уже 0.41092:

batch_size = 512
model = Sequential()
model.add(Conv1D(8,3, activation=act.relu, padding='causal', input_shape=(data.timesteps, data.features)))
model.add(Conv1D(4,3, activation=act.relu, padding='causal', dilation_rate=2))
model.add(Conv1D(4,3, activation=act.relu, padding='causal', dilation_rate=4))
model.add(GRU(16))
model.add(Dense(1, activation=act.sigmoid))

Далее добавляем нормализацию и dropout:

Код

batch_size = 512
dropout_rate = 0.05

def channel_normalization(x):
	max_values = K.max(K.abs(x), 2, keepdims=True) + 1e-5
	out = x / max_values
	return out

model = Sequential()
if(data.timesteps > 1):
	model.add(Conv1D(16,2, activation=act.relu, padding='valid', input_shape=(data.timesteps, data.features)))
	model.add(Lambda(channel_normalization))
	model.add(SpatialDropout1D(dropout_rate))
	model.add(Conv1D(16,1, padding='valid'))
	for i in range(1, 6):
		model.add(Conv1D(16,2, activation=act.relu, padding='valid', dilation_rate=2**i))
		model.add(Lambda(channel_normalization))
		model.add(SpatialDropout1D(dropout_rate))
		model.add(Conv1D(16,1, padding='valid'))
	model.add(Flatten())
else:
	model.add(Flatten(input_shape=(data.timesteps, data.features)))
	model.add(Dense(256, activation=act.relu))
model.add(Dense(1, activation=act.sigmoid))

Применяя эту модель в том числе и на ранних шагах (с t = 1), получаем score = 0.53578.

Gradient Boosting Machines

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

Gradient Boosting — это метод ML, на выходе которого мы получаем ансамбль простых моделей (чаще всего деревьев решений). За счет большого количества таких простых моделей и происходит оптимизация функции потерь. Более подробно о Gradient Boosting можно почитать, например, здесь.

В качестве реализации GBM использовали lightgbm — достаточно известный фреймворк от Microsoft.

Взятые отсюда модель и предобработка данных сразу дают score около 0.64:

Код

def prepare_data(marketdf, newsdf):
	# a bit of feature engineering
	marketdf['time'] = marketdf.time.dt.strftime("%Y%m%d").astype(int)
	marketdf['bartrend'] = marketdf['close'] / marketdf['open']
	marketdf['average'] = (marketdf['close'] + marketdf['open'])/2
	marketdf['pricevolume'] = marketdf['volume'] * marketdf['close']
	newsdf['time'] = newsdf.time.dt.strftime("%Y%m%d").astype(int)
	newsdf['assetCode'] = newsdf['assetCodes'].map(lambda x: list(eval(x))[0])
	newsdf['position'] = newsdf['firstMentionSentence'] / newsdf['sentenceCount']
	newsdf['coverage'] = newsdf['sentimentWordCount'] / newsdf['wordCount']
	# filter pre-2012 data, no particular reason
	marketdf = marketdf.loc[marketdf['time'] > 20120000]
	# get rid of extra junk from news data
	droplist = ['sourceTimestamp','firstCreated','sourceId','headline','takeSequence','provider','firstMentionSentence', 'sentenceCount','bodySize','headlineTag','marketCommentary','subjects','audiences','sentimentClass', 'assetName', 'assetCodes','urgency','wordCount','sentimentWordCount']
	newsdf.drop(droplist, axis=1, inplace=True)
	marketdf.drop(['assetName', 'volume'], axis=1, inplace=True)
	# combine multiple news reports for same assets on same day
	newsgp = newsdf.groupby(['time','assetCode'], sort=False).aggregate(np.mean).reset_index()
	# join news reports to market data, note many assets will have many days without news data
	return pd.merge(marketdf, newsgp, how='left', on=['time', 'assetCode'], copy=False)

import lightgbm as lgb
print ('Training lightgbm')
# money
params = {
	"objective" : "binary",
	"metric" : "binary_logloss",
	"num_leaves" : 60,
	"max_depth": -1,
	"learning_rate" : 0.01,
	"bagging_fraction" : 0.9, # subsample
	"feature_fraction" : 0.9, # colsample_bytree
	"bagging_freq" : 5, # subsample_freq
	"bagging_seed" : 2018,
	"verbosity" : -1
}
lgtrain, lgval = lgb.Dataset(Xt, Yt[:,0]), lgb.Dataset(Xv, Yv[:,0])
lgbmodel = lgb.train(params, lgtrain, 2000, valid_sets=[lgtrain, lgval], early_stopping_rounds=100, verbose_eval=200)

Предобработка здесь уже включает данные новостей, объединяя их с рыночными данными (однако делая это достаточно наивно, учитывается только один код актива из всех, которые упоминаются в новости). Этот вариант предобработки я взял за основу во всех последующих решениях.

Немного добавив фич (firstMentionSentence, marketCommentary, sentimentClass), а также заменив метрику на ROC AUC, получим score 0.65389.

Ансамбль

Следующим удачным решением было использование ансамбля, состоящего из нейросетевой модели и GBM (хотя «ансамбль» — это громкое название для двух моделей). Результирующее предсказание получаем усреднением предсказаний двух моделей, применяя, таким образом, механизм soft voting. Такое решение позволило получить score 0.66879.

Exploratory Data Analysis и Feature Engineering

Еще одна вещь, с которой нужно было начать, это EDA. Прочитав, что важно понять корреляцию между фичами, строим такую картинку (картинки в этом разделе кликабельны):

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

Видно, что после 10-дневного периода зависимость сильно падает. Вероятно, это и обуславливает неплохую работу GBM, учитывающих только фичи с 10-дневной задержкой (которые уже есть в изначальном наборе данных).

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

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

market_x = market_data.loc[:,features]
scaler = StandardScaler()
scaler.fit(market_x)
market_x = scaler.transform(market_x)

pca = PCA(.95)
pca.fit(market_x)
market_pca = pca.transform(market_x)

Посмотрим, какие фичи генерирует PCA:

Видим, что метод работает не очень хорошо на наших данных, так как итоговая корреляция новых фич с целевым значением невелика.

Тонкая настройка и нужна ли она

Многие модели ML имеют достаточно большое количество гиперпараметров, то есть, «настроек» самого алгоритма. Их можно подбирать вручную, но есть также и механизмы автоматического подбора. Для последнего существует библиотека hyperopt, которая реализует два алгоритма подбора — случайный поиск и Tree-structured Parzen Estimator (TPE). Пробовал оптимизировать:

параметры lightgbm (тип алгоритма, кол-во листьев, learning rate и другие),
параметры нейросетевых моделей (кол-во фильтров TCN, кол-во блоков памяти GRU, dropout rate, learning rate, тип решателя).

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

Финальное решение

По правилам соревнования участники могут выбрать два решения для финального этапа. Мои финальные решения практически одинаковы и содержат ансамбль из двух моделей — GBM и многослойного GRU. Отличия только в том, что одно решение совсем не использует новостные данные, а другое использует, но только для нейросетевой модели.

Решение с новостными данными:

Импорты

import numpy as np
import pandas as p
import itertools
import functools
from kaggle.competitions import twosigmanews
from sklearn.preprocessing import StandardScaler, LabelEncoder
import tensorflow as tf
from keras.models import Sequential, Model
from keras.layers import Dense, GRU, LSTM, Conv1D, Reshape, Flatten, SpatialDropout1D, Lambda, Input, Average
from keras.optimizers import Adam, SGD, RMSprop
from keras import losses as ls
from keras import activations as act
import keras.backend as K
import lightgbm as lgb

Предобработка данных

# fix random
from numpy.random import seed
seed(42)
from tensorflow import set_random_seed
set_random_seed(42)

env = twosigmanews.make_env()
(market_train_df, news_train_df) = env.get_training_data()

def cleanData(market_data, news_data):
    market_data = market_data[(market_data['returnsOpenNextMktres10'] <= 1) & (market_data['returnsOpenNextMktres10'] >= -1)]
    return market_data, news_data

def prepareData(marketdf, newsdf, scaler=None):
    print('Preparing data...')
    
    print('...preparing features...')
    marketdf = marketdf.copy()
    newsdf = newsdf.copy()
    # a bit of feature engineering
    marketdf['time'] = marketdf.time.dt.strftime("%Y%m%d").astype(int)
    marketdf['bartrend'] = marketdf['close'] / marketdf['open']
    marketdf['average'] = (marketdf['close'] + marketdf['open'])/2
    marketdf['pricevolume'] = marketdf['volume'] * marketdf['close']
    
    newsdf['time'] = newsdf.time.dt.strftime("%Y%m%d").astype(int)
    newsdf['position'] = newsdf['firstMentionSentence'] / newsdf['sentenceCount']
    newsdf['coverage'] = newsdf['sentimentWordCount'] / newsdf['wordCount']

    # filter pre-2012 data, no particular reason
    marketdf = marketdf.loc[marketdf['time'] > 20120000]
    
    # get rid of extra junk from news data
    droplist = ['sourceTimestamp','firstCreated','sourceId','headline','takeSequence','provider',
                'sentenceCount','bodySize','headlineTag', 'subjects','audiences',
                'assetName', 'wordCount','sentimentWordCount', 'companyCount',
                 'coverage']
    newsdf.drop(droplist, axis=1, inplace=True)
    marketdf.drop(['assetName', 'volume'], axis=1, inplace=True)
    
    # unstack news
    newsdf['assetCodes'] = newsdf['assetCodes'].apply(lambda x: x[1:-1].replace("'", ""))
    codes = []
    indices = []
    for i, values in newsdf['assetCodes'].iteritems():
        explode = values.split(", ")
        codes.extend(explode)
        repeat_index = [int(i)]*len(explode)
        indices.extend(repeat_index)
    index_df = p.DataFrame({'news_index': indices, 'assetCode': codes})
    newsdf['news_index'] = newsdf.index.copy()

    # Merge news on unstacked assets
    news_unstack = index_df.merge(newsdf, how='left', on='news_index')
    news_unstack.drop(['news_index', 'assetCodes'], axis=1, inplace=True)
    
    # combine multiple news reports for same assets on same day
    newsgp = news_unstack.groupby(['time','assetCode'], sort=False).aggregate(np.mean).reset_index()
    
    # join news reports to market data, note many assets will have many days without news data
    res = p.merge(marketdf, newsgp, how='left', on=['time', 'assetCode'], copy=False) #, right_on=['time', 'assetCodes'])
    res.marketCommentary = res.marketCommentary.astype(float)
    
    targetcol = 'returnsOpenNextMktres10'
    target_presented = targetcol in res.columns
    features = [col for col in res.columns if col not in ['time', 'assetCode', 'universe', targetcol]]
    
    print('...scaling...')
    if(scaler == None):
        scaler = StandardScaler()
        scaler = scaler.fit(res[features])
    res[features] = scaler.transform(res[features])

    print('...done.')
    return type('', (object,), {
        'scaler': scaler,
        'data': res,
        'x': res[features],
        'y': (res[targetcol] > 0).astype(int).values if target_presented else None,
        'features': features,
        'samples': len(res),
        'assets': res['assetCode'].unique(),
        'target_presented': target_presented
    })

def generateTimeSeries(data, n_timesteps=1):
    
    data.data[data.features] = data.data[data.features].fillna(data.data[data.features].mean())
    #data.data[data.features] = data.data[data.features].fillna(0)
    assets = data.data.groupby('assetCode', sort=False)
    
    def grouper(n, iterable):
        it = iter(iterable)
        while True:
           chunk = list(itertools.islice(it, n))
           if not chunk:
               return
           yield chunk
    
    def sample_generator():
        while True:
            for assetCode, days in assets:
                x = days[data.features].values
                y = (days['returnsOpenNextMktres10'] > 0).astype(int).values if data.target_presented else None
                for i in range(0, len(days) - n_timesteps + 1):
                    yield (x[i: i + n_timesteps], y[i + n_timesteps - 1] if data.target_presented else 0)
    
    def batch_generator(batch_size):
        for batch in grouper(batch_size, sample_generator()):
            yield tuple([np.array(t) for t in zip(*batch)])
    
    n_samples = functools.reduce(lambda x,y : x + y, map(lambda t : 0 if len(t[1]) + 1 <= n_timesteps else len(t[1]) - n_timesteps + 1, assets))

    return type('', (object,), {
        'gen': batch_generator,
        'timesteps': n_timesteps,
        'features': len(data.features),
        'samples': n_samples,
        'assets': list(map(lambda x: x[0], filter(lambda t : len(t[1]) + 1 > n_timesteps, assets)))
    })

Нейросетевая модель

def buildRNN(timesteps, features):
    i = Input(shape=(timesteps, features))
    x1 = Lambda(lambda x: x[:,:,:13])(i)
    x1 = Conv1D(16,1, padding='valid')(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = GRU(10)(x1)
    x1 = Dense(1, activation=act.sigmoid)(x1)
    x2 = Lambda(lambda x: x[:,:,13:])(i)
    x2 = Conv1D(16,1, padding='valid')(x2)
    x2 = GRU(10, return_sequences=True)(x2)
    x2 = GRU(10, return_sequences=True)(x2)
    x2 = GRU(10, return_sequences=True)(x2)
    x2 = GRU(10)(x2)
    x2 = Dense(1, activation=act.sigmoid)(x2)
    x = Average()([x1, x2])
    model = Model(inputs=i, outputs=x)
    return model

def train_model_time_series(model, data, val_data=None):
    print('Building model...')
    batch_size = 4096
    
    optimizer = RMSprop()
    
    # define roc_callback, inspired by https://github.com/keras-team/keras/issues/6050#issuecomment-329996505
    def auc_roc(y_true, y_pred):
        value, update_op = tf.metrics.auc(y_true, y_pred)
        metric_vars = [i for i in tf.local_variables() if 'auc_roc' in i.name.split('/')[1]]
        for v in metric_vars:
            tf.add_to_collection(tf.GraphKeys.GLOBAL_VARIABLES, v)
        with tf.control_dependencies([update_op]):
            value = tf.identity(value)
            return value
    
    model.compile(loss=ls.binary_crossentropy, optimizer=optimizer, metrics=['binary_accuracy', auc_roc])
    
    print(model.summary())
    
    print('Training model...')
    
    if(val_data == None):
        model.fit_generator(data.gen(batch_size),
            epochs=8,
            steps_per_epoch=int(data.samples / batch_size),
            verbose=1)
    else:
        model.fit_generator(data.gen(batch_size),
            epochs=8,
            steps_per_epoch=int(data.samples / batch_size),
            validation_data=val_data.gen(batch_size),
            validation_steps=int(val_data.samples / batch_size),
            verbose=1)

    return type('', (object,), {
        'predict': lambda x: model.predict_generator(x, steps=1)
    })

GBM модель

def train_model(data, val_data=None):
    print('Building model...')
    
    params = {
        "objective" : "binary",
        "metric" : "auc",
        "num_leaves" : 60,
        "max_depth": -1,
        "learning_rate" : 0.01,
        "bagging_fraction" : 0.9,  # subsample
        "feature_fraction" : 0.9,  # colsample_bytree
        "bagging_freq" : 5,        # subsample_freq
        "bagging_seed" : 2018,
        "verbosity" : -1 }
    
    ds, val_ds = lgb.Dataset(data.x.iloc[:,:13], data.y), lgb.Dataset(val_data.x.iloc[:,:13], val_data.y)
    print('...training...')
    model = lgb.train(params, ds, 2000, valid_sets=[ds, val_ds], early_stopping_rounds=100, verbose_eval=100)
    print('...done.')
    
    return type('', (object,), {
        'model': model,
        'predict': lambda x: model.predict(x.iloc[:,:13], num_iteration=model.best_iteration)
    })

Обучение

n_timesteps = 30

market_data, news_data = cleanData(market_train_df, news_train_df)
dates = market_data['time'].unique()
train = range(len(dates))[:int(0.85*len(dates))]
val = range(len(dates))[int(0.85*len(dates)):]

train_data_prepared = prepareData(market_data.loc[market_data['time'].isin(dates[train])], news_data.loc[news_data['time'] <= max(dates[train])])
val_data_prepared = prepareData(market_data.loc[market_data['time'].isin(dates[val])], news_data.loc[news_data['time'] > max(dates[train])], scaler=train_data_prepared.scaler)

model_gbm = train_model(train_data_prepared, val_data_prepared)

train_data_ts = generateTimeSeries(train_data_prepared, n_timesteps=n_timesteps)
val_data_ts = generateTimeSeries(val_data_prepared, n_timesteps=n_timesteps)
rnn = buildRNN(train_data_ts.timesteps, train_data_ts.features)
model_rnn = train_model_time_series(rnn, train_data_ts, val_data_ts)

Предсказание

def make_predictions(data, template, model):
    if(hasattr(data, 'gen')):
        prediction = (model.predict(data.gen(data.samples)) * 2 - 1)[:,-1]
    else:
        prediction = model.predict(data.x) * 2 - 1
    predsdf = p.DataFrame({'ast':data.assets,'conf':prediction})
    template['confidenceValue'][template['assetCode'].isin(predsdf.ast)] = predsdf['conf'].values
    return template

day = 1
days_data = p.DataFrame({})
days_data_len = []
days_data_n = p.DataFrame({})
days_data_n_len = []
for (market_obs_df, news_obs_df, predictions_template_df) in env.get_prediction_days():
    print(f'Predicting day {day}')
    days_data = p.concat([days_data, market_obs_df], ignore_index=True, copy=False, sort=False)
    days_data_len.append(len(market_obs_df))
    days_data_n = p.concat([days_data_n, news_obs_df], ignore_index=True, copy=False, sort=False)
    days_data_n_len.append(len(news_obs_df))
    data = prepareData(market_obs_df, news_obs_df, scaler=train_data_prepared.scaler)
    predictions_df = make_predictions(data, predictions_template_df.copy(), model_gbm)
    if(day >= n_timesteps):
        data = prepareData(days_data, days_data_n, scaler=train_data_prepared.scaler)
        data = generateTimeSeries(data, n_timesteps=n_timesteps)
        predictions_df_s = make_predictions(data, predictions_template_df.copy(), model_rnn)
        predictions_df['confidenceValue'] = (predictions_df['confidenceValue'] + predictions_df_s['confidenceValue']) / 2
        days_data = days_data[days_data_len[0]:]
        days_data_n = days_data_n[days_data_n_len[0]:]
        days_data_len = days_data_len[1:]
        days_data_n_len = days_data_n_len[1:]
    env.predict(predictions_df)
    day += 1

env.write_submission_file()

Решение без новостных данных:

Код (только отличающийся метод)

def buildRNN(timesteps, features):
    i = Input(shape=(timesteps, features))
    x1 = Lambda(lambda x: x[:,:,:13])(i)
    x1 = Conv1D(16,1, padding='valid')(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = GRU(10)(x1)
    x1 = Dense(1, activation=act.sigmoid)(x1)
    model = Model(inputs=i, outputs=x1)
    return model

Оба решения выдавали похожий результат (около 0.69) на первом этапе конкурса, что соответствовало 566 месту из 2927. После первого месяца новых данных позиции в списке участников сильно перемешались, и решение с новостными данными оказалось на 65 месте из оставшихся 697 команд с результатом 3.19251, а что будет в течение следующих пяти месяцев, не знает никто.

Что еще попробовал

Кастомные метрики

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

Метрика для lightgbm:

def sharpe_metric(y_pred, train_data):
	y_true = train_data.get_label() * 2 - 1
	std =  np.std(y_true * y_pred)
	mean = np.mean(y_true * y_pred)
	sharpe = np.divide(mean, std, out=np.zeros_like(mean), where=std!=0)
	return "sharpe", sharpe, True

Проверка показала, что такая метрика работает в данной задаче хуже, чем AUC.

Механизм внимания

Механизм внимания позволяет нейросети сфокусироваться на «наиболее важных» признаках в исходных данных. Технически, внимание представляется вектором весов (чаще всего получаемых при помощи полносвязного слоя с активацией softmax), которые перемножаются на выход другого слоя. Я использовал реализацию, в которой внимание применяется к временной оси:

def buildRNN(timesteps, features):
    
    def attention_3d_block(inputs):
        a = Permute((2, 1))(inputs)
        a = Dense(timesteps, activation=act.softmax)(a)
        a = Permute((2, 1))(a)
        mul = Multiply()([inputs, a])
        return mul
    
    i = Input(shape=(timesteps, features))
    x1 = Lambda(lambda x: x[:,:,:13])(i)
    x1 = Conv1D(16,1, padding='valid')(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = attention_3d_block(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = attention_3d_block(x1)
    x1 = GRU(10, return_sequences=True)(x1)
    x1 = attention_3d_block(x1)
    x1 = GRU(10)(x1)
    x1 = Dense(1, activation=act.sigmoid)(x1)
    model = Model(inputs=i, outputs=x1)
    return model

Эта модель выглядит достаточно красиво, однако прироста score такой подход не дал, получилось около 0.67.

Что не успел сделать

Несколько направлений, которые выглядят перспективными:

более предметно разобраться с механизмом внимания,
попробовать применить автоэнкодеры,
попробовать онлайновое обучение,
внимательнее разобраться с объединением новостных и рыночных данных, а также с предобработкой новостей.

Выводы

Наше приключение подошло к концу, можно попробовать подвести итоги. Конкурс оказался сложным, но мы смогли не ударить в грязь лицом. Это намекает, что порог входа в ML не так уж высок, но, как в любом деле, настоящая магия (а её в машинном обучении предостаточно) доступна уже профессионалам.

Итоги в цифрах:

Максимальный score на первом этапе: ~0.69 против ~1.5 у первого места. Что-то вроде среднего по больнице, значение 0.7 преодолели немногие, максимальный score публичного решения так же ~0.69, чуть больше моего.
Место на первом этапе: 566 из 2927.
Score на втором этапе: 3.19251 после первого месяца.
Место на втором этапе: 65 из 697 после первого месяца.

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

Ссылки

Финальное решение с использованием новостей
Two Sigma: Using News to Predict Stock Movements — страница конкурса
Keras — нейросетевой фреймворк
LightGBM — GBM фреймворк
Scikit-learn — библиотека алгоритмов машинного обучения
Hyperopt — библиотека для оптимизации гиперпараметров
Статья о WaveNet

Комментарии (13)

roryorangepants
06.03.2019 09:45
#19841658
Вероятно, перекрестная валидация не очень подходит, так как поведение рынка сильно отличается на ранних этапах тренировочных данных и на оценочном периоде.

Почитайте что ли про time series cross-validation.

а также заменив метрику на ROC AUC

Если вы говорите не про таргет для оптимизации, а просто про метрику, по которой вы отбирали модель, то, что ROC AUC у вас лучше коррелировал с ЛБ, скорее всего означает, что у вас была неудачная валидация.

Таким образом, для этой задачи тонкая настройка не очень подходит, так как приводит к переобучению модели.

Или, опять таки, вы плохо сделали валидацию.
Ну, хотя настраивать сеть hyperopt-ом — это и правда в целом не лучшее решение, наверное.

Конкурс оказался сложным, но мы смогли не ударить в грязь лицом.

максимальный score публичного решения так же ~0.69, чуть больше моего.

Пардон, но если у вас скор ниже лучшего паблик кернела, я бы не называл это «удалось не ударить в грязь лицом».
1. niklyaev Автор
  06.03.2019 10:00
  #19841740
  Почитайте что ли про time series cross-validation.
  
  Спасибо, почитаю!
  
  Или, опять таки, вы плохо сделали валидацию.
  
  Это вот очень вероятно. Пока более-менее разобрался с тем, как надо делать, конкурс уже подходил к концу :)
  
  я бы не называл это «удалось не ударить в грязь лицом».
  
  По поводу этого пункта — да, спорное утверждение, согласен. В качестве небольшого оправдания: очень-очень не хотелось целиком тащить это конкретное решение. В нем также использовался только lightgbm, но были получше подобраны фичи и гиперпараметры, однако именно с такими фичами и гиперпараметрами мои модельки не заработали нормально.
1. markmipt
  06.03.2019 12:21
  #19842608
  +2
  Про time series cross-validation плюсуюсь.
  И как раз о ней же: в лучшем публичном решении была использована обычная валидация на основе train_test_split, которая не подходит для данной задачи. Там явный оверфиттинг, но насколько я понял людей это не останавливало и многие брали за основу это решение, не думая поменять валидацию на time series. Поэтому про «удалось не ударить в грязь лицом» не согласен.
  
  А вообще, добавлю от себя для сравнения чисел из статьи: свое решение на lightgbm я смог написать на 0.78. Что еще интереснее, публичный скор 0.7 мне удалось получить вообще не используя машинное обучение, а просто написав фиксированную формулу для расчета торговых весов. И эта стратегия хорошо вела себя за весь backtest период в 10 лет (в среднем получалось 0.65 +- 0.1 std)
  1. roryorangepants
    06.03.2019 14:22
    #19843418
    Там явный оверфиттинг
    
    Ну, там тоже плохая валидация, само собой, но ничем не хуже, чем у автора статьи. А по результатам текущего ЛБ куда лучше даже.
    
    niklyaev Автор
    06.03.2019 14:32
    #19843498
    По результатам текущего ЛБ у того решения сейчас 1.43030, насколько я понимаю. Хотя, это ничего особо не значит, судя по всему.
    
    roryorangepants
    06.03.2019 14:49
    #19843566
    Конкретно у того, может, 1.43. Я не знаю, какой конкретно кернел был тогда топовым, к сожалению.
    Но вот, например, паблик кернел на топ-21.

superstarstas306
06.03.2019 10:26
#19841878
Раз уж зашел разговор про Кагл
Уже неделю пытаюсь ворваться в соревнование по землетрясениям, используя генетическое программирование, каменный цветок не выходит. Если у кого-то из господ комментаторов есть пример использования библиотеки DEAP для проблемы регрессии, буду чрезмерно благодарен.

slava_k
06.03.2019 18:56
#19844810
Интересно, почему все так стремятся прогнозировать именно цену, а не, к примеру, диапазон волатильности такой цены в не очень отдаленном будущем (в рамках месяца). Вероятность построить модель, дающую точный прогноз именно по ожиданию цены на основе только ценового ряда прошлого — довольно низкая, и всегда будет убывать с ростом временнОго интервала. В то же время, оценивая свойства волатильности в определенные интервалы времени, в том числе выявляя длительность таких интервалов с волатильностью выше определенного уровня — уже это даст более интересный и практичный результат. Сезонность волатильности будет еще ближе к полезным данным, которые можно использовать для улучшения торговых систем или вообще брать за основу принятия решений.

К примеру, нет особого смысла входить в рынок, когда волатильность по активу достаточно низкая. Также если есть возможность оценить объемы торгов в определенный момент (доступ к листингу сделок/«стакан»), нет особого смысла влезать в тонкий рынок и брать повышенный риск выноса цены за ожидаемые рамки оценки рисков. В итоге, анализируя не просто цену, а то, как она изменяется во времени, в определенные периоды (часы/дни), сопоставляя с параметрами торговой системы и планируемых объемов сделок — можно ограничиться наилучшими периодами времени для торгов под нужную торговую стратегию.

Если инструмент как-то связан с сезонностью потребления (к примеру, фьючерс на зерно, кофе, пшеницу, бензин), то можно найти интересные периодичности в волатильности такого инструмента. Также практически по любому инструменту, у которого есть рынок фьючерсных и опционных контрактов — стоит обратить внимание на поведение цен таких деривативов и почему цены довольно интересно себя ведут близ уровней, кратных 50% (100%, 150% и далее) от среднедневной волатильности, и какая за таким поведением может быть реальная логика.

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

С точки зрения организаторов подобных конкурсов, я вижу стремление задешево ($100k) найти интересное и нестандартное решение/подход к прогнозированию, и уже сам этот подход потом эксклюзивно использовать в качестве торговой модели/робота или перепродать какому-либо инвестфонду/банку (уже на порядки дороже $100k).
1. eugene_bb
  06.03.2019 19:16
  #19844890
  +1
  Может я неправ, но мне кажется что в Kaggle (по крайней мере в последнее время) никто не пытается решить реальные задачи.
  
  Мне например это очень не понятно и не ясно в чём смысл соревнований. Реклама фирмы? Поиск кандидатов на позиции? Продажа дополнительных сервисов?
  
  В Two Sigma конкурсе было какое-то дурацкое ограничение на закачку данных на локальные компьютеры (я думаю всё равно все их закачали хотя и не было прямого способа)
  
  Данные не очень хорошего качества и не полные и т.п. Мне кажется столько реально интересных задач можно было бы решить, а тысячи людей занимаются выжиманием последних капелек за счёт стэкинга.
  
  Как площадка начального обучения Kaggle очень хорош, но потом лучше своё время тратить на что-то реально полезное.
  1. slava_k
    06.03.2019 19:35
    #19844962
    Реклама да, одна из очевидных причин. Ну и headhunting под задачи компании-организатора. На схожих с kaggle площадках рунета аналогичная ситуация.
    
    С учетом предыдущего ознакомления с данными по тем же конкурсам с нефтяными вышками — базовый набор данных безнадежно бессмысленен для практического применения, т.к. по такому сильно ограниченному срезу данных невозможно даже оценить рамки уровней добычи/отдачи с определенных вышек. На прогнозирование отдачи влияет вагон различных факторов, начиная от геологии и динамики ее изменения, температуры, где вышка находится, что за команда работает на кустах (опыт), проводятся ли близ вышки специфические работы (гидроразрыв, закачка пустот водой для подъема пластов). Консервация соседних вышек может приводить как к росту отдачи, так и к снижению из-за изменения давления подземных слоев. И все эти факторы существенно влияют на итоговый прогноз, без них достоверность и качество прогнозирования будет крайне низкой и бессмысленной.
    
    Ну и еще момент. Конкурс может так ничем и не закончиться, уровни сложности могут повысить, изменить условия и прочее.
    
    Не хочу ни в коем случае бросить тень на сервисы таких конкурсов, но, при конкурсах с ограничением kernel-only у организатора довольно широкие возможности для анализа кода и результатов участников. Теоретически, может быть получен интересный и нестандартный подход с высокой стоимостью для организатора конкурса и далее повышен уровень призеров так, чтобы никто его не прошел.
    
    С вашим выводом в конце про пользу полностью согласен и поддерживаю.
  1. markmipt
    06.03.2019 19:46
    #19845016
    +1
    Я согласен с вами частично. Из того, что я наблюдал по остальным соревнованиям — выигравают те, кто делают стэкинг и тому подобные вещи.
    
    Но конкретно в этом соревновании был большой простор для создания различных стратегий и оно непохоже на остальные, поэтому я и принял участие впервые на Kaggle.
    
    Но насчет реальных задач — я тоже сомневаюсь что Two Sigma будет использовать какое-либо из решений, скорее всего пиар/ищут сотрудников себе.
  1. roryorangepants
    06.03.2019 21:36
    #19845432
    Есть конкурсы и с реальными задачами. Например, в Human Protein Atlas был отдельный призовой трек для «быстрого» решения в продакшен.

slava_k
06.03.2019 19:21
#19844908
Предлагаю посмотреть на использование новостей под иным углом: не столь важно содержание новости, сколько «стоимость» источника публикации такой новости и ранее построенная несложная модель линейного предсказания влияния новости на рост волатильности связанного с новостью инструмента.

К примеру, возьмем за источник новости — твиттер Трампа и содержимое новости как "… нефть ...". Из текста новости можно легко выделить инструмент через прямое упоминание слова или его наиболее стОящих в плане цитирования синонимов. Далее, зная время публикации такого твита, смотрим на поведение целевого инструмента за 5 минут до публикации и 5 минут после публикации. Оцениваем рост волатильности инструмента как реакцию на новость и, если на интервале 2 лет (существование источника новостей с повышенной «ценностью») имеется закономерность в росте волатильности на подобные новости — уже на этом можно строить модель принятия решения по технике торговли волатильностью (схема покупки/продажи базового актива и обратного «хеджа» через дериватив на такой базовый актив).

Да, подобные стратегии не особо очевидные, но такой подход существует и работает с момента появления крупных новостных агентств и их популярных социальных преемников типа твиттера.

Для улучшения качества прогнозирования цены или диапазона ее волатильности можно использовать описанный выше механизм. Важно правильно оценить стоимость источников новостей, он должен быть максимально влиятельным на актив, иметь наивысшую «стоимость». К примеру, твиттер Маска и TSLA. И далее уже изучать статистические закономерности именно изменения волатильности актива и новостей. Не особо важен даже сентимент новости, плохо это или хорошо, важен сам факт воздействия и его размер. Оценку сентимента можно проводить по словарю наиболее однозначных слов и их синонимов, но это не всегда качественно работает. Это вообще по сути будет подстройка алгоритма под язык и манеру речи источника новостей, потому оценка сентимента новостей вряд ли когда-то станет статистически надежной.