О программировании, алгоритмах и не только

Random Forest с примерами на R

2012-05-16T19:09:00.001+04:00

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

В чем же идея такого замечательного алгоритма? Идея проста: допустим у нас есть какой-то очень слабый алгоритм, скажем, дерево принятия решений. Если мы сделаем очень много разных моделей с использованием этого слабого алгоритма и усредним результат их предсказаний, то итоговый результат будет существенно лучше. Это, так называемое, обучение ансамбля в действии. Алгоритм Random Forest потому и называется "Случайный Лес", для полученных данных он создает множество деревьев приятия решений и потом усредняет результат их предсказаний. Важным моментом тут является элемент случайности в создании каждого дерева. Ведь понятно, что если мы создадим много одинаковых деревьев, то результат их усреднения будет обладать точностью одного дерева.

Как он работает? Предположим, у нас есть некие данные на входе. Каждая колонка соответствует некоторому параметру, каждая строка соответствует некоторому элементу данных.

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

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

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

Этот удивительно простой алгоритм и самым сложным шагом в его реализации является построение дерева decision tree. И не смотря на свою простоту он дает очень хорошие результаты в реальных задачах. С практической точки зрения у него есть одно огромное преимущество: он почти не требует конфигурации. Если мы возьмем любой другой алгоритм машинного обучения, будь то регрессия или нейронная сеть, они все имеют кучу параметров и их надо уметь подбирать под конкретную задачу. Алгоритм random forest имеет по-сути лишь один параметр: размер случайного подмножества выбираемого на каждом шаге построения дерева. Этот параметр важен, но даже значения по-умолчанию позволяют получить весьма приемлемые результаты.

Перейдем теперь к примеру на R
В R random forest реализован в пакете под названием randomForest (удивительно). Для сегодняшнего примера я воспользуюсь теми же данными что и в статье про деревья принятия решений. Попробуем предсказать изменения в объемах торгов акциями Газпрома на следующий день.

#Где взять пакет rusquant и как его установить можно посмотреть тут http://www.algorithmist.ru/p/rusquant.html 
library(rusquant)
getSymbols('GAZP', from='2003-05-01', to='2012-05-09', src='Finam')

#для повторяемости результатов
set.seed(1)
label <- (Vo(GAZP) - lag(Vo(GAZP)))/lag(Vo(GAZP))
names(label) <- c('label')
data <- lag(GAZP)

data <-cbind(lag(label), lag(label, k=2), lag(label, k=3), lag(label, k=4), lag(label, k=5), RSI(Cl(data)), MACD(Cl(data)), volatility(data), factor(weekdays(index(label))))
names(data) <- c('Volume1', 'Volume2', 'Volume3', 'Volume4', 'Volume5', 'RSI', 'macd', 'signal', 'volatility', 'weekdays')
data <- cbind(label, data)
data <- data[complete.cases(data),]

# 80% данных используем для тренировки модели, 20% для тестирования
split <- runif(dim(data)[1]) > 0.2
train <- data[split,]
test <- data[!split,]

Код для тренировки random forest в R очень прост:

#random forest
library(randomForest)
rf <- randomForest(label ~ ., train)
predictions <- predict(rf, test)
print(sqrt(sum((as.vector(predictions - test$label))^2))/length(predictions))

В результате получается точность 0.02157 что полностью согласуется с результатами полученными с помощью ансамбля деревьев. Можно ли как-то подкрутить этот алгоритм? Если не считать вопросы скорости работы, есть только один параметр который имеет смысл менять mtry - он определяет, сколько столбцов тестируется при поиске следующего разбиения, для каждого дерева. Но его изменения не сказываются на результате принципиальным образом. Лучшее значение, которое мне удалось выжать, при mtry=3 было 0.02149. Подобрать оптимальное значение mtry поможет функция tuneRF Кроме того, random forest позволяет оценить качество используемых предикторов.

importance(rf)

В нашем примере получилось:

           IncNodePurity
Volume1        136.11990
Volume2         24.80163
Volume3         15.96727
Volume4         19.43955
Volume5         32.80988
RSI             18.90843
macd            18.68504
signal          18.35107
volatility      21.04426
weekdays        13.54460

Т.е. самым значимым предиктором оказалось изменение объема в предыдущий день, все остальные несли гораздо меньше информации. Это функция алгоритма random forest может помочь при выборе предикторов и провести простейший feature selection. Помимо перечисленных возможностей в пакете randomForest есть еще много интересных функций, часть из них позволяют лучше распределять работу при тренировке случайного леса, часть визуализировать данные, вставлять пропущенные значения итп. В следующей статье я расскажу о том, как мы боролись за скорость при использовании случайного леса в последнем соревновании.

Деревья принятия решений с примерами на R

2012-05-10T09:19:00.001+04:00

Деревья принятия решения это алгоритм из серии "для домохозяек", обычно в том или ином виде с ним встречались все. Мое начальство очень любит такого рода картинки, наверняка и вам доводилось видеть подобное.

Это и есть дерево принятия решений. В силу своей наглядности этот алгоритм часто используется для визуализации данных. И тем не менее, несмотря на свою простоту, алгоритм является чрезвычайно важным в data mining и лежит в основе многих продвинутых методов, таких как random forest, boosted trees итп.

Итак, как же работает дерево принятия решений? Каждая вершина представляет собой точку разбиения данных. Вершина определяет некоторый простой критерий по которому мы делим данные на части. Далее я буду рассматривать только бинарные деревья, где каждая вершина делит данные на две части, именно они используются на практике. Кроме того, я буду рассматривать только простые критерии разбиения данных на части, а именно разбиение по единственному параметру. Таким образом, в каждой вершине дерева мы разбиваем данные на две части согласно значению одного из параметров:

Дальше для каждой части мы повторяем операцию и в итоге получаем дерево. Для простоты реализации часто используют рекурсию. Здесь неочевидными остаются только три вещи:
1. Как выбирается параметр по которому происходит разбиение
2. Как выбирается пороговое значения параметра
3. Когда алгоритм должен остановиться

Ответы на эти вопросы зависят от реализации. Существует несколько популярных реализаций со странными названиями типа ID3, C4.5, CARD. Возьмем, к примеру, CARD, как он отвечает на эти вопросы? Параметр по которому происходит разбиение и пороговое значение этого параметра находятся методом полного перебора. Т.е. мы просто пробуем все возможные способы разбиения данных по каждому из параметров и оцениваем, насколько разбиение было удачным. А как оценить удачность разбиения? Очень просто, выполняя каждое разбиение мы получаем некоторое промежуточное дерево и его возможности классификации или регрессии нужно оценить. Вот тут то на помощь приходят понятия вроде информационной энтропии и коэффициентов Джини. В моей практике, степень удачности разбиения часто определялась условиями задачи, а именно, критерием, который мы пытаемся оптимизировать.

Когда алгоритму остановиться? Есть два варианта: когда дальнейшие разбиения не способны улучшить дерево, например если в текущей вершине все элементы имеют одинаковое значение. Либо когда мы достигли определенного уровня, например в вершине 5 или меньше элементов (этот критерий используется по-умолчанию в реализации random forest в R)

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

Мы плавно пришли к примеру. Для примера попробуем предсказывать изменения в объемах торгов акциями Газпрома на рынке ММВБ на следующий день (цены предсказывать дело не благодарное, поэтому объемы). В качестве предикторов будем использовать изменения в объемах торгов за предыдущие 5 дней, несколько стандартных технических индикаторов и день недели. Итак, скачаем и сформируем данные:

#Где взять пакет rusquant и как его установить можно посмотреть тут http://www.algorithmist.ru/p/rusquant.html 
library(rusquant)
getSymbols('GAZP', from='2003-05-01', to='2012-05-09', src='Finam')

#для повторяемости результатов
set.seed(1)
label <- (Vo(GAZP) - lag(Vo(GAZP)))/lag(Vo(GAZP))
names(label) <- c('label')
data <- lag(GAZP)

data <-cbind(lag(label), lag(label, k=2), lag(label, k=3), lag(label, k=4), lag(label, k=5), RSI(Cl(data)), MACD(Cl(data)), volatility(data), factor(weekdays(index(label))))
names(data) <- c('Volume1', 'Volume2', 'Volume3', 'Volume4', 'Volume5', 'RSI', 'macd', 'signal', 'volatility', 'weekdays')
data <- cbind(label, data)
data <- data[complete.cases(data),]

# 80% данных используем для тренировки модели, 20% для тестирования
split <- runif(dim(data)[1]) > 0.2
train <- data[split,]
test <- data[!split,]

Для начала проведем небольшую проверку на адекватность (sanity check) чтобы понимать, будет ли наш алгоритм вообще предсказывать что-то или нет. Предположим, что объемы торгов вообще не меняются и посмотрим, какую ошибку даст это предположение:

#бенчмарк
print(sqrt(sum((as.vector(test$label))^2))/length(test$label))

В результате получаем: 0.02397 Теперь посмотрим что способно предсказать одно единственное дерево:

library(rpart)

#тренируем дерево принятия решений
tree <- rpart(label ~ ., train)
predictions <- predict(tree, test)
print(sqrt(sum((as.vector(predictions - test$label))^2))/length(predictions))

Результат: 0.02442 - чуть-чуть хуже чем бенчмарк. Не очень то оптимистично. Как можно улучшить этот результат? Как я уже писал, тенденция последнего времени, если одна модель работает плохо, надо объединить много моделей. Поступим следующим образом, будем из наших данных выбирать подмножества: примерно 70% столбцов и примерно 70% рядов, и потом для этих подмножеств генерить деревья. Дальше, получив много таких деревьев для разных подмножеств мы просто посчитаем арифметическое среднее для результатов их предсказаний. Смысл генерации подмножеств в том, чтобы как-то рандомизировать деревья, ведь если просто сгенерить много одинаковых деревьев и усреднить их предсказания то результаты никак не улучшатся. Здесь есть очень важный момент в вызове функции rpart следует передавать параметр control=rpart.control(cp=.0), он необходим для того чтобы алгоритм строил полные деревья не удаляя вершины, которые не приводят к существенному улучшению результата. Без этой опции деревья будут получаться слишком похожими и усреднение ничего не даст.

ntrees <- 300
res <- numeric(ntrees)
prd <- numeric(nrow(test))
for(i in 1:ntrees){
    #выберем 70% рядов и колонок из множества данных
    x <- runif(nrow(train))>0.3;
    y <- runif(ncol(train))>0.3;
    #обязательно включим label
    y[1] <- TRUE
    traindata <- train[x,y]
    #генерим полное дерево не оптимизируя процесс
    tree <- rpart(label ~ ., traindata, control=rpart.control(cp=.0))
    #усредняем предсказание со всеми предыдущими деревьями
    prd <- prd + predict(tree, test)
    predictions <- prd / i
    #оцениваем ошибку
    res[i] <- sqrt(sum((as.vector(predictions - test$label))^2))/length(predictions)
    print(res[i])
}

plot(res, type='l')

В результате для 300 деревьев мы получаем ошибку: 0.02170 что намного лучше чем одно дерево и бенчмарк. Теперь результаты стали гораздо интересней. Посмотрим как изменяется ошибка в зависимости от количества деревьев:

На графике видно типичное поведение ошибки при тренировке ансамбля. Видно что где-то к 30 деревьям результат стабилизируется и дальнейшее добавление деревьев слабо влияет на ошибку. По-сути мы только что построили подобие алгоритма random forest для регрессии. Аналогичным образом можно построить алгоритм для классификации. Когда я запускал эти же алгоритмы на предварительно обработанных данных для kaggle.com, первый дал результат 0.85076, что принесло бы примерно 100-ое место, а второй всего лишь с десятью деревьями дал бы 0.73623 очков и 10-ое место. Что не сильно хуже полноценного random forest, но об этом в следующий раз.

Kaggle.com мы третьи!

2012-05-03T17:53:00.003+04:00

И вот, наконец, закончилось мое первое соревнование в kaggle.com с неожиданным для меня результатом: третье место. Это, пожалуй, были самые лучшие два месяца с точки зрения изучения data mining-а и я постараюсь в ближайшее время поделиться своими находками.

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

Сложности добавляла низкая ликвидность на рынке бондов. Некоторые контракты не торговались неделями и, понятно, что ни о каком High Frequency Trading тут и речи идти не может.

Моя первая идея была использовать random forest, благо что и пример организаторов использовал его и в целом было понятно что при такой неопределенности данных эта модель должна как-минимум работать неплохо. Поэтому, первым, что я сделал, это нормализовал цены относительно curve_based_price так, чтобы избавиться насовсем от абсолютных значений в долларах США. Это простое изменение в тот момент подняло меня на 5-ую строчку в рейтинге и мотивировало на дальнейшие исследования. Следующим существенным изменением было разбить тренировочные данные на несколько множеств (по типу сделки и по типу бонда). Основной эффект этого разбиения заключался в том, что эффективное время обучения random forest на большом множестве данных упало многократно. Это изменение подняло меня на первое место где я и находился около месяца. За это время успев получить приглашение от организаторов и посетить Stanford Conference of Quantitative Finance.

А дальше был творческий кризис, я пробовал добавлять различные предикторы, пробовал оптимизировать модель, удалять outliers и даже переписал реализацию random forest на Java надеясь что оптимизировав функцию оценки смогу получить лучший результат. Но тщетно, некоторые изменения давали мне небольшое улучшение, но большинство не влияли на результат либо влияли негативно. Постепенно участники стали формировать команды и я потерял свое первое место и рисковал оказаться в районе 5-го места в финале. Что было бы весьма обидно ведь денежный приз дается только за первые три места. Но в последние три дня нам удалось объединить усилия с участником VikP, посчитать арифметическое среднее наших моделей и выйти на третье место.

Немного разочаровало то, что модели победителей это просто смесь большого количества не самых лучших индивидуальных моделей. Т.е. рецепт победы в общем случае прост: натренируйте кучу разных моделей а потом с помощью линейной регресси найдите оптимальную комбинацию. Такой подход мне кажется немного ущербным. Ведь в основе часто лежат модели типа random forest и GBM которые уже сами по себе являются ансамблями большого числа еще более слабых моделей. Но, видимо, в настоящее время, этот подход является лучшим из известного, достаточно например посмотреть на описания решений победителей промежуточных туров Heritage Health Prize.

Моя новая жизнь

2012-01-12T17:43:00.000+04:00

Вот собственно и всё. 17-ти часовой перелет позади, Россия осталась за океаном. А в окно уютной 2-ух спальной квартиры на нас смотрит Сан-Франциско, знаменитая Кремниевая долина, Калифорния, США. Да, это и есть та самая причина, по которой я практически не писал последнее время. Мы переехали.

Всё это началось еще в апреле 2011 года, когда я проходил телефонное интервью в компании Zynga. Тогда это все казалось какой-то игрой не имеющей отношения к реальности и я и представить себе не мог, во что это выльется. В июне 2011 года Zynga приехали в Москву и провели серию собеседований, рассматривалось около 60 кандидатов прошедших телефонное интервью и из них было отобрано около 15 человек (точное число не знаю, кто-то потом передумал, кто-то сразу отказался). Интервью оказалось неожиданно простым. Ни тебе задачек на программирование, ни заковыристых вопросов про форму люков, в основном проверялись способности болтать. А знания, на мой взгляд, оценивались лишь поверхностно.

А дальше началась канитель. Сначала мы ждали результатов, потом офера, потом одобрение LCA, потом одобрения петиции на визу, потом документы из США, потом очередь в посольстве, потом дополнительную проверку, потом визу. Временами мне казалось, что я готов все бросить и забить. Временами я сомневался, а нужна ли нам эта Америка ведь и в России не плохо. Весь процесс занял где-то около полугода, в итоге, в середине декабря мы получили визы и начали готовиться к отъезду.

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

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

Дивиденды rusquant 0.3.3

2011-11-29T06:57:00.000+04:00

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

library(rusquant)
divs <-  getAllDividends()

Все дивиденды, можно получить командой getAllDividends(allYears = TRUE)
Возвращаемый объект содержит следующие данные:

Symbol - id акции на бирже
Name - полное название эмитента
BoardDate - дата собрания совета директоров
RegistryDate - дата закрытия реестра
RecommendedDivs - рекомендуемые дивиденды
Divs - выплаченные дивиденды

Для примера использования, распечатаем дивидендную доходность российских компаний. В качестве базовой цены, берем цену закрытия акции в день закрытия реестра. Почему-то на сайте тройки этой информации нет, а она ведь гораздо интересней чем абсолютные величины дивидендов.
Внимание! Код выполняется довольно долго, т.к. для каждой акции требуется сделать запрос на сервера finam и получить её стоимость.

result <- NULL
    for(i in (1:length(divs[,1]))){
       d <- divs[i,]
       if (d$Divs>0){
          try({
              quotes <- getSymbols(d$Symbol, src="Finam", from="2010-01-01", auto.assign=FALSE)
             if (!is.nan(quotes)){ 
                 price <- Cl(quotes[d$RegistryDate])
                 if (length(price)>0){
                      dd <- d$Divs
                      result <- rbind(result, data.frame(d$Symbol, d$Name, d$RegistryDate, as.numeric(dd)/as.numeric(price), stringsAsFactors=FALSE))
                 }
             }
          }, silent=TRUE)
       }
    }
colnames(result) <- c('Symbol', 'Name', 'RegistryDate', 'Divs')
result[with(result, order(Divs, decreasing=TRUE)),c('Symbol','RegistryDate','Divs')]

Аналогично можно построить статистику для прошлых лет.

Регуляризация линейной регрессии с примерами на R

2011-11-28T06:00:00.000+04:00

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

Дабы сильно не усложнять задачу, сегодня я воспользуюсь синтетическими данными. Как обычно, пример будет на языке программирования R. Итак, данные:

"gen" <- function(x){
    return(3*x^2 + 2*x + 2 + rnorm(length(x))*0.5)
}

#Генерируем тренировочные данные
X <- runif(6)
Y <- gen(X)
#Данные для кросс-валидации
Xcv <- runif(50)
Ycv <- gen(Xcv)

Тут все просто, данные распределены согласно некой выдуманной формуле 3*x²+2*x+2 плюс некоторая погрешность, случайным образом распределенная. Тренировочных данных было умышленно взято очень мало, для наглядности (в реальных задачах их все равно, как правило, не хватает, так что это ничуть не мешает обобщению результатов).

Допустим, мы не знаем как на самом деле распределены наши данные, но мы хотим построить какую-нибудь модель. Простейшая линейная регрессия по переменной x нам не подойдет, т.к. линейная регрессия предполагает только линейные зависимости. Поэтому, мы попробуем построить полиномиальную модель. И для этого добавим в наши данные все возможные степени x от 1 до 5.

#Линейная регрессия по полиномиальному набору данных
train <- data.frame(Y, X, X^2, X^3, X^4, X^5)
colnames(train) <- c('Y', 'X', 'X2', 'X3', 'X4', 'X5')
simple <- lm(Y ~ X+X2+X3+X4+X5, train)
error <- sum((predict(simple, train)-Y)^2)/length(Y)
cat("Train error: ",error,"\n")

cv <- data.frame(Xcv, Xcv^2, Xcv^3, Xcv^4, Xcv^5)
colnames(cv) <- c('X', 'X2', 'X3', 'X4', 'X5')
error <- sum((predict(simple, cv)-Ycv)^2)/length(Ycv)
cat("Cross-validation error: ",error,"\n")

#Построим кривую, описывающую наше решение
x <- (1:100)/100
test = data.frame(x, x^2, x^3, x^4, x^5)
names(test) <- c('X', 'X2', 'X3', 'X4', 'X5')
y0 <- predict(simple, test)

plot(X, Y, ylim=range(y0), xlim=c(0,1))
lines(x,y0, col='red')

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

Train error:  4.951098e-13 
Cross-validation error:  6.751967

Это, так называемый, overfitting. И он абсолютно ожидаемый, ведь через 6 точек на плоскости всегда можно провести полином 5-ой степени.

Тут, читатель может заметить, а к чему, собственно, нам полином 5-ой степени, при таком количестве точек разумней было бы использовать меньше степеней. Согласен, случай выдуманный, но аналогичная проблема возникает в случае сильно-многомерных данных. Например, при обработке текстов. Где переменная X запросто может быть вектором с 10000 элементов.

Что делать?

Для начала, посмотрим на формулу, которой мы пытаемся описать наши данные:
y_t = a₀ + a₁*x + a₂*x + a₃*x + a₄*x + a₅*x
Мы пытаемся подобрать эти коэффициенты a_i таким образом, чтобы минимизировать ошибку т.е. минимизировать выражение:
L = Σ(y_t-y)²
А что если в это выражение добавить еще какой-нибудь член, чтобы как-то уменьшить величины коэффициентов a_i? Например, вместо L можно было бы минимизировать выражение:

L1 = Σ(y_t-y)² + λ*Σ|a_i|

или другой вариант

L2 = Σ(y_t-y)² + λ*Σ(a_i)²

Что мы здесь сделали? Мы добавили в наше выражение штраф за большие значения a_i. И величина этого штрафа пропорциональна величине параметра λ, с помощью которого мы теперь можем настраивать наш алгоритм.

Первый вариант носит название L1-регуляризация (в английской литературе LASSO regression), второй вариант L2-регуляризация или регуляризация Тихонова (в английской литературе ridge regression)

Оба варианта реализованы в R. К сожалению, реализация немного не последовательна в используемых входных параметрах и возвращаемых значениях. Тем не менее, мы можем посмотреть её использование чтобы сравнить результат работы.

L1-регуляризация или Lasso regression

L1 регуляризация реализована в пакете lars.

#Lasso regression
#Требуется установить пакет lars если его нет
#install.packages('lars')
library(lars)

train <- cbind(X, X^2, X^3, X^4, X^5)
colnames(train) <- c('X', 'X2', 'X3', 'X4', 'X5')
lasso <- lars(train, Y, type='lasso')

cv <- cbind(Xcv, Xcv^2, Xcv^3, Xcv^4, Xcv^5)
colnames(cv) <- c('X', 'X2', 'X3', 'X4', 'X5')

x <- (1:100)*(max(X)*1.1)/100
test <- cbind(x, x^2, x^3, x^4, x^5)
colnames(test) <- c('X', 'X2', 'X3', 'X4', 'X5')


stats <- NULL
lambda <- 0
plot(X, Y, ylim=c(-5, 10), pch=20, col='red')
points(Xcv, Ycv, pch=20, col='blue')
ls <- NULL
cs <- NULL
for(i in 1:15){
    Yp <- predict(lasso, train, s=lambda, type='fit', mode='lambda')
    trainError <- sum((Yp$fit-Y)^2)/length(Y)

    Yp <- predict(lasso, cv, s=lambda, type='fit', mode='lambda')
    cvError <- sum((Yp$fit-Ycv)^2)/length(Ycv)

    stats <- rbind(stats, c(lambda, trainError, cvError))
    if (i%%5==1){
        #Построим кривую, описывающую наше решение
        y0 <-  predict(lasso, test, lambda, type='fit', mode='lambda')
        lines(x,y0$fit, col=i)
        ls <- c(ls, paste("lambda=",lambda,sep=''))
        cs <- c(cs, i)
    }

    lambda <- ifelse(lambda==0,0.00000001, lambda*10)
}
legend("topleft", c("Train", "Cross-validation"), pch=20, col=c('red', 'blue'))
legend("bottomright",inset=0.05, legend=ls, pch=20, col=cs, text.col=cs, bty="n")



plot(stats[,2], ylim=range(stats[,2:3]), type='l', col='red', ylab="error", xlab="log(lambda)")
lines(stats[,3], col='blue')

#Оптимальное значение lamda=1
coef(lasso, s=1, mode='lambda')

На графике показано несколько решений, для разных значений λ. При λ=0 получается обычная линейная регрессия, которую я описал выше. При очень больших значениях λ получается горизонтальная линия, т.е. решение при котором коэффициенты a_i=0, для всех i>0. Оба крайних случая очевидно не являются оптимальными решениями. Как же найти оптимальное? Для этого построим график ошибки.

Здесь, на графике, красным цветом отмечена ошибка для тренировочного множества, а синим цветом ошибка для cross-validation множества. Видно, что ошибка тренировочного множества сначала очень мала (overfitting), но постепенно возрастает (underfitting). В то же время ошибка cross-validation имеет выраженный минимум, который для нас и будет оптимальным значением λ В нашем случае, это оптимальное значение оказалось равным 10.

У L1 регуляризации есть одна полезная особенность. Многие коэффициенты оказываются в точности равны 0. Например, в нашем случае:

X        X2        X3        X4        X5 
0.2989565 0.0000000 0.0000000 0.0000000 5.2679886

Т.е. L1 регуляризация может быть использована еще и как алгоритм feature selection.

Регуляризация Тихонова или Ridge regression

Этот алгоритм реализован в пакете MASS. К сожалению, реализация не совсем стандартная, например, пакет MASS не предоставляет нам функции predict которой так привычно пользоваться тем, кто знаком с R. Тем не менее, не представляет никакой сложности самостоятельная реализация этой функции через умножение матриц.

#Ridge regression
#Требуется установить пакет MASS если его нет
#install.packages('MASS')
library(MASS)

train <- data.frame(Y, X, X^2, X^3, X^4, X^5)
colnames(train) <- c('Y', 'X', 'X2', 'X3', 'X4', 'X5')

stats <- NULL
lambda <- 0
x <- (1:100)/100
test <- as.matrix(cbind(1, x, x^2, x^3, x^4, x^5))

plot(X, Y, ylim=c(-5,15), pch=20, col='red')
points(Xcv, Ycv, pch=20, col='blue')
ls <- NULL
cs <- NULL
for(i in (1:15)) {

    ridge <- lm.ridge(Y ~ X+X2+X3+X4+X5, train, lambda=lambda)
    Yp <- as.matrix(cbind(1, X, X^2, X^3, X^4, X^5)) %*% as.matrix(coef(ridge))
    trainError <- sum((Yp - Y)^2)/length(Y)

    cv <- as.matrix(cbind(1, Xcv, Xcv^2, Xcv^3, Xcv^4, Xcv^5))
    Yp <- cv %*% as.matrix(coef(ridge))
    cvError <- sum((Yp-Ycv)^2)/lentgh(Ycv)

    stats <- rbind(stats, c(lambda, trainError, cvError))
    if (i%%5==1){
        #Построим кривую, описывающую наше решение
        y0 <- test %*% as.matrix(coef(ridge))
        lines(x,y0, col=i)
        ls <- c(ls, paste("lambda=",lambda,sep=''))
        cs <- c(cs, i)
    }
    lambda <- ifelse(lambda==0,0.00000001, lambda*10)

}
legend("topleft", c("Train", "Cross-validation"), pch=20, col=c('red', 'blue'))
legend("bottomright",inset=0.05, legend=ls, pch=20, col=cs, text.col=cs, bty="n")


#Рисуем изменение ошибки CV
plot(stats[,2], ylim=range(stats[,2:3]), type='l', col='red', ylab="error", xlab="log(lambda)")
lines(stats[,3], col='blue')

#Оптимальное значение lamda=10
ridge <- lm.ridge(Y ~ X+X2+X3+X4+X5, train, lambda=10)
coef(ridge)

Аналогично случаю L1 регрессии, здесь мы тоже можем определить оптимальные значения параметра λ

Оптимальное значение λ в данном случае оказалось равным 1. И соответствующие коэффициенты, как видим, совсем не равны 0.

X        X2        X3        X4        X5 
3.0257674 0.7951164 0.7637522 0.8631815 1.0352908 1.2643954

Естественно методы регуляризации не ограничиваются L1 и L2-регуляризациями. Существует еще несколько популярных подходов, например алгоритм elastic net, включающий в себя как L1 так и L2.

[Text mining] Индекс настроений фондовых рынков

2011-09-30T07:14:00.000+04:00

В прошлом году я уже писал о нашумевшей работе американцев, которые умудрились предсказывать индекс S&P500 по сообщениям в twitter. Это было лишь начало, сейчас в интернете можно найти множество статей посвященных анализу текстовых данных для предсказания котировок на бирже:

"Trading Strategies to Exploit Blog and News Sentiment" (Zhang, Skiena 2010)
"The Predictive Power of Financial Blogs" (Frisbee 2010)
"An analysis of verbs in financial news articles and their impact on stock price" (Schumaker 2010)

Но это все там, за океаном. А что же у нас, в России?

Искать правду в твиттере для российского рынка, мне показалось, бесполезным. Во-первых твиттером у нас пользуются по-прежнему единицы, а во-вторых у американцев, так или иначе, в фондовый рынок вовлечено почти все население, у нас же это развлечение для избранных. Зато в России есть раскрученная социальная сеть трейдеров comon.ru. На этом ресурсе публикуется несколько десятков статей ежедневно и, что самое главное, они по теме. Что может быть лучше для задачи анализа текста?

Итак, сказано - сделано. Я решил проверить связано ли как-то настроение блогеров comon.ru с движениями на рынке и, самое интересное, могут ли эти данные быть использованы для предсказания движений рынка. Первым делом, было выкачано около 10 тыс. статей от 100 авторов с comon.ru, вместе с датой публикации и именем автора. Потом текст статей был обработан: удалены все элементы HTML-разметки и стоп-слова а у всех оставшихся слов была выделена основа.

Дальше надо было научиться как-то характеризовать текст. Для начала я выбрал всего две категории: негативную и позитивную. Под позитивной понимаются тексты в которых идет речь о росте, восходящем тренде и повышении. Негативные тексты, напротив, изобилуют словами типа дефолт, кризис, обвал и понижение. Интересно, что негативных слов получилось вдвое больше чем позитивных.

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

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

Когда я прогнал полученный алгоритм по набору текстов, стало понятно, что результаты крайне нестабильны. И понятно почему. Одни авторы используют много негативных слов, другие мало, одни оптимисты, другие пессимисты, одни пишут часто и регулярно, другие редко, у одних тексты состоят из нескольких тысяч слов, у других они больше похоже на посты в твиттере. Чтобы сгладить влияние отдельных авторов, было проведено нормирование индекса настроений, по каждому автору, предположив, что среднее настроение автора за весь период жизни на комоне нейтрально. Попутно, кстати, получилась интересная метрика, показывающая насколько одни авторы более пессимистичны или оптимистичны чем другие.

Далее, после усреднения данных по дням, получился некоторый индекс, назовем его индексом настроений comon.ru. Это число, которое характеризует среднее настроение блогеров comon.ru в этот день. Интересным оказался график изменения индекса настроения:

Из графика видно, что вначале года, несмотря на ряд провалов, среднее значение индекса настроений большую часть времени больше нуля. Во второй половине года настроения изменились и стали все больше негативными. Что в это время происходило на бирже, комментировать не надо. Кроме того, индекс настроений успешно предсказал два самых крупных обвала года.

Следующим шагом, я проверил а сколько можно было бы заработать если покупать вечером того дня когда индекс настроений больше нуля и закрывать позицию на следующий день. И аналогично для отрицательных значений индекса настроений. Оказалось, что среднее дневное изменение индекса ММВБ в день, следующий за днем с позитивным настроением составляет: 0.13%. А в дни следующие за днями с негативным настроением -0.34%. Результат мне показался вдохновляющим. И было решено проверить математически гипотезу о том, что индекс настроения комона предсказывает значения ММВБ. Для этого, я воспользовался методом описанным в оригинальной работе, а именно, причинным тестом Грейнджера:

Model 1: MICEX ~ Lags(MICEX, 1:1) + Lags(SIGNAL, 1:1)
Model 2: MICEX ~ Lags(MICEX, 1:1)
  Res.Df Df      F  Pr(>F)  
1    178                    
2    179 -1 2.9539 0.08741 .

P-value теста составило 0.087, что меньше 0.1. Значит наша гипотеза о том что индекс настроений комона предсказывает изменения индекса ММВБ на следующий день верна с вероятностью 90%. Что является неплохим результатом, вполне сравнимым с результатом в оригинальной работе с твиттером.

Большим недостатком полученного индекса является его высокая волатильность. По-идее, добавление других источников данных и улучшение алгоритма анализа настроений должны её сгладить. С другой стороны, увеличение количества источников данных может снизить точность результатов.

Что использовалось для данной работы?
1. Граббер текстов с comon.ru и предварительная обработка были написан на Java
2. Вся математика, анализ и построение графиков выполнялись на R с использованием библиотек lmtest и rusquant

Какие варианты развития?
1. В отличии от твиттера, количество авторов комона которые пишут регулярно не велико. Можно ввести что-то вроде предсказательного рейтинга автора, тем самым улучшив общий прогноз.
2. Учитывать оценки других пользователей и комментарии к каждой статье.
3. Добавить анализ морфологии русского языка, учитывать время в котором написаны отдельные предложения статьи.

Этот безумный R

2011-09-20T06:55:00.000+04:00

Сначала мне казалось, что это простой и примитивный язык, для небольших скриптов, с уклоном в статистику. Потом мне стало казаться, что это самый обычный язык программирования, с уклоном в статистику. Сейчас мне кажется, что это гремучая смесь классических императивных языков вроде C++ или Java и матерой функциональщины вроде Haskell или Lisp. И всё это R.

На самом деле, действительно, R многое заимствовал из Lisp. Функции в R являются полноценными сущностями. Их можно не только вызывать и определять но и передавать в качестве параметров или составлять из них сложные выражения. В R доступны замыкания и отложенные вычисления. При этом, в отличии, например, от Haskell, R не страдает от функциональной чистоты и допускает функции с побочными эффектами. Я уже не говорю о том, что в R еще и объектно-ориентированный язык. Правда, признаюсь честно, реализация объектно-ориентированного программирования в R больше похожа на помойку. В R есть несколько ООП-систем и каждая живет своей жизнью. Все это многообразие, делает работу сложней и интересней. Нередко, разобравшись в концепции, начинаешь понимать, как тот или иной код можно написать короче и эффективней. Некоторыми такими концепциями я бы и хотел сегодня поделиться.

Векторизация

Векторизация это первое с чем сталкиваются все R разработчики не всегда, правда, осознанно. Векторы, в R повсюду. Даже когда мы просто объявляем в консоли число 1, мы создаем вектор с единственным элементом.

> 1
[1] 1

Как следствие этого, большинство функций в R работают с векторами. И обычно об этом даже не нужно думать. Например:

> c(1,2,3)+c(4,5,6)
[1] 5 7 9
> c(1,2,3)*c(4,5,6)
[1]  4 10 18

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

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

> a <- c(1,-2,-3,4)
> ifelse(a>0, 1, 0)
[1] 1 0 0 1

Recycling rule

Эта концепция важна в сочетании с векторизацией. Что будет, если в примере выше два вектора будут разной длины? Оказывается, что если размер более длинного вектора кратен размеру короткого, то короткий вектор будет автоматически повторен нужное количество раз. Если же не кратен, то произойдет то же самое, но будет показано предупреждение (которое, правда, можно подавить)

> c(1,2)+c(1,2,3,4)
[1] 2 4 4 6
> c(1,2,3)+c(1,2,3,4)
[1] 2 4 6 5
Предупреждение
In c(1, 2, 3) + c(1, 2, 3, 4) :
  длина большего объекта не является произведением длины меньшего объекта

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

> c(1,2,3)+1
[1] 2 3 4

Мы ведь помним, что единица это тоже вектор, с длиной 1.

Всё является функцией

Да, R действительно функциональный язык программирования. Функцией является всё, даже скобочки {, [, операторы присвоения <- и условные операторы if, while и т.п. Даже ключевое слово, объявляющее функции function - тоже является функцией. И главное, их можно переопределить!

> `if`
.Primitive("if")
> `if` <- function(...) {"Hello"}
> if (TRUE) 1 else 0
[1] "Hello"
> `if` <- .Primitive("if")
> if (TRUE) 1 else 0
[1] 1

Для чего это может понадобиться не берусь даже представить.

Отложенные вычисления

R, как и Haskell, использует концепцию отложенных (lazy) вычислений. Это значит, что аргументы функций не вычисляются до тех пор пока их значения не потребуются. На самом деле, незаметно для программиста, R передает в функции специальные объекты "обещания" (promise), эти объекты содержат в себе всё, что необходимо, чтобы по требованию вычислить значение. Если значение было вычислено, то оно запоминается в объекте promise и используется в дальнейшем по мере необходимости.

> fun <- function(x,a,b) {if (x) a else b}
> fun(TRUE, 1, factorial(10000))
[1] 1
> fun(FALSE, 1, factorial(10000))
[1] Inf
Предупреждение
In factorial(10000) : значение в 'gammafn' -- за пределами

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

Окружение

Окружение (environment) это часть механизма замыканий (closures), реализованного в R. Любая переменная доступная в окружении на момент вызова функции сохраняется и передается вместе с окружением в вызываемую функцию. Но самое неожиданное, это то, что с помощью функции assign можно определить некую переменную не только в текущем окружении, но и в родительском, например так:

> foo <- function(){assign("a", 123, envir=parent.frame())}
> bar <- function(){foo(); message(a)}
> bar()
123

Зачем это может быть нужно, мне сложно представить. В любом случае, полагаю, что использовать эту возможность следует с большой осторожностью. В целом, язык программирования R оказался намного больше, разнообразней и сложнее чем казалось в начале. Здесь я перечислил далеко не все его возможности. Более полное описание языка можно найти в его спецификации.

Тиковые данные rusquant v0.3

2011-09-19T06:26:00.000+04:00

Очередное обновление rusquant. Теперь стали доступны тиковые данные.

Как обычно, обновиться, тем у кого уже есть, или установить тем, у кого еще нет, можно командой в консоли R:

install.packages("rusquant", repos="http://R-Forge.R-project.org")

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

Для получения тиковых данных можно воспользоваться знакомой функцией getSymbols()

> library(rusquant)
>  getSymbols("SBER", from=Sys.Date()-1, src="Finam", period="tick")
[1] "SBER"
>  SBER[1:5]
                    SBER.Close SBER.Volume
2011-09-14 09:59:59      79.24          30
2011-09-14 09:59:59      79.24          30
2011-09-14 09:59:59      79.24          10
2011-09-14 09:59:59      79.24          20
2011-09-14 09:59:59      79.24         500

Сразу отмечу, что вызов выполняется достаточно долго, т.к. объем получаемых данных очень большой. Для удобства, сохранена система именования OHLC данных. Т.е. цены можно получить всё той же функций Cl, а объем Vo.

Если же есть необходимость преобразовать тиковые данные в OHLC, то нет ничего проще, в модуле xts для этого есть ряд удобных функций типа to.minutes, to.hourly, to.minutes10 и более общая функция to.period

> to.period(SBER, "mins", k=5)
                    SBER.Open SBER.High SBER.Low SBER.Close SBER.Volume
2011-09-14 09:59:59     79.24     79.24    79.24      79.24      287520
2011-09-14 10:04:58     79.47     79.68    79.11      79.52     4042460
2011-09-14 10:09:59     79.52     79.66    79.34      79.35     3095280

Доступные значения периодов: “secs”, “seconds”, “mins”, “minutes”, “hours”, “days”, “weeks”, “months”, “quarters”, и “years”. Хотя, я сомневаюсь, что последними кто-то будет пользоваться. k=5 - показывает какое количество периодов объединять в группу. В данном примере мы получили пятиминутные данные.

Геомагнитные бури на фондовом рынке

2011-09-13T06:50:00.000+04:00

Сегодняшняя тема немного несерьезная. А может быть наоборот очень серьезная. Так сложилось, что во времена моей прыщавой юности, когда рухнул совок и было не во что больше верить, люди начали слушать всяких Чумаков, Кашпировских и прочих телепузиков. Кажется именно тогда особой популярностью пользовались геомагнитные бури и бабушки у подъезда с радостью списывали очередной приступ остеохондроза на повышенную геомагнитную активность. Наверное, именно такое отношение и привело в последствии к предвзятому взгляду на прогнозы геомагнитной активности (я то её влияние никогда не чувствовал). Тем не менее, явление то существует и объективно регистрируется.

Есть масса исследований, утверждающих, что геомагнитные бури влияют на настроение людей. Люди склонны пессимистичней оценивать происходящее вокруг них во время и после бури, нежели в спокойные дни. Теперь подумаем, если люди глядят на вещи с пессимизом, причем это люди по всей земле. Значит, это должно непременно сказываться на фондовом рынке. Логично предположить, что магнитные бури должны приводить к сокращению вложений в рискованные активы, такие как акции. А раз так, то этим можно было бы воспользоваться для построения собственной торговой стратегии. И вот, в статье "Playing the Field: Geomagnetic Storms and the Stock Market" такое исследование уже проводилось и результаты мне показались достаточно интересными чтобы проверить их на практике.

Итак, сила магнитной бури измеряется двумя величинами Kp и Ap. Какой за ними тайный смысл я не вникал, создалось ощущение, что в любом случае они очень сильно взаимосвязаны. И не важно, какую из них взять за основу. Данные по геомагнитным бурям аж с 1932 года доступны всем желающим на публичном ftp сервере: ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/INDICES/KP_AP/ Но мы, не для того столько времени посвятили языку программирования R, чтобы качать вручную и разбирать их экселем. Напишем маленькую программку, которая эти данные соберет и представит нам в удобоваримом формате:

"getGeomagnetic" <-
function(year){
    require(xts)
    url <- paste("ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/INDICES/KP_AP/", year, sep='')
    tmp <- tempfile()
    download.file(url, destfile=tmp, quiet=TRUE)
    fr <- readLines(con = tmp, warn=FALSE)
    unlink(tmp)
    date <- as.Date(strptime(substr(fr,1,6), "%y%m%d"))
    apIndex <- as.numeric(substr(fr, 56, 58))
    kpIndex <- as.numeric(substr(fr, 29, 31))
    result <- xts(as.matrix(cbind(apIndex, kpIndex)), date)
    result <- result[complete.cases(result)]
    return(result)
}

Совместим, для начала, графики геомагнитной обстановки и индекса RTS (нам понадобится модуль rusquant о котором я уже писал в одной из прошлых статей):

gm <- c(getGeomagnetic(2010),getGeomagnetic(2011))

library(rusquant)
getSymbols("RTSI", from="2010-01-01", src="Finam")
data <- cbind(RTSI, gm)

chartSeries(data, TA=NULL, name="RTSI")
addTA(data$kpIndex, legend="Геомагнитная обстановка")

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

Значит ли это, что обвалы можно спрогнозировать по геомагнитным бурям или нет? Для начала повторим эксперимент, проведенный в статье, в применении к Российскому рынку. А именно, посмотрим на среднее изменение величины индекса РТС в день геомагнитной бури и в течении 6 дней после нее. Будем считать магнитной бурей такие дни, когда значение величины Ap больше 29 (именно такое значение использовалось в оригинальной работе).

library(rusquant)
getSymbols("RTSI", from="2002-01-01", src="Finam")
gm <- getGeomagnetic(2002)
for(i in (2003:2011)){
	gm <- c(gm, getGeomagnetic(i))
}

GMS <- ifelse(gm$apIndex>29, 1, 0)
GMS <- apply(Lag(GMS, k=(0:6))>0, 1, any)
data <- cbind(GMS, ROC(Cl(RTSI)))
names(data) <- c("GMS", "RTSI")
mean(data[data$GMS==1]$RTSI, na.rm=TRUE)
mean(data[data$GMS==0]$RTSI, na.rm=TRUE)

Получается, что в случае спокойной геомагнитной обстановки среднее увеличение индекса РТС составляет 0.092% в день, в случае же магнитной бури 0.01% в день, т.е. в 9 раз меньше!

Попробуем теперь подтвердить или опровергнуть гипотезу о влиянии геомагнитных бурь математически.

Когда-то давно я рассказывал о причинном тесте Грейнджера и показал тогда, как его можно реализовать на Java. В этот раз мы проведем аналогичный тест, но на R и увидим, насколько это проще.

Причинный тест Грейнджера это, по сути, авторегрессия. Т.е. линейная регрессия временного ряда на его собственных значениях в прошлом. С точки зрения математической формулы, это выглядит так:
Y_i = μ_i + ∑α_kY_i-k + ∑β_kX_i-k + ε_i
где, Y_i - значение исследуемой переменной Y в момент времени i, в нашем случае, это индекс РТС
X_i - значение переменной X, влияние которой мы хотим проверить, в момент времени i. В нашем случае это сила геомагнитной бури.
k - временная задержка. Она определяет, влияние какой давности мы хотим проверить.

Мы будем проверять влияние "на следующий день" после магнитной бури. Реализация теста Грейнджера доступна в пакете lmtest.

#install.packages("lmtest")
library(lmtest)

GMS <- ifelse(gm$apIndex>29, 1, 0)
data <- cbind(GMS, ROC(Cl(RTSI)))
names(data) <- c("GMS", "RTSI")
grangertest(RTSI~GMS, data=data, order=1)

Получаем

Granger causality test

Model 1: RTSI ~ Lags(RTSI, 1:1) + Lags(GMS, 1:1)
Model 2: RTSI ~ Lags(RTSI, 1:1)
  Res.Df Df      F Pr(>F)
1   2374                 
2   2375 -1 1.3522  0.245

Значение p-value=0.245 не позволяет нам утверждать о наличии зависимости между индексом RTS и геомагнитными бурями. Однако, если немного увеличить порог критерия магнитной бури, то получим:

GMS <- ifelse(gm$apIndex>33, 1, 0)
data <- cbind(GMS, ROC(Cl(RTSI)))
names(data) <- c("GMS", "RTSI")
grangertest(RTSI~GMS, data=data, order=1)

результат:

Granger causality test

Model 1: RTSI ~ Lags(RTSI, 1:1) + Lags(GMS, 1:1)
Model 2: RTSI ~ Lags(RTSI, 1:1)
  Res.Df Df      F  Pr(>F)  
1   2374                    
2   2375 -1 2.9949 0.08366

В данном случае, значение p-value=0.08366 меньше 0.1 и значит, что с вероятностью 90% мы не можем отвергать гипотезу о влиянии геомагнитных бурь на фондовый рынок. Выходит, что все-таки влияет. Но влияют только очень сильные магнитные бури, которых с 2002 года было всего 137. Можно ли это как-то использовать? Мне тоже интересно. Буду ждать ближайшей бури, чтобы проверить.

Обновление rusquant 0.2

2011-09-06T06:57:00.000+04:00

Небольшая библиотека для работы с Российскими фондовыми рынками в R пополнилась новыми возможностями. Теперь, помимо доступа к данным с Финама, прямо из R вы можете получить доступ к котировкам фьючерсов и опционов Forts, доступным на сайте rts.ru

Итак, для тех что еще не установил себе rusquant - ставим, для тех, кто уже установил - обновляемся командой:

install.packages("rusquant", repos="http://R-Forge.R-project.org")

Что новенького?

1. Стала доступной доска опционов:

library(rusquant)
rts <- getOptionChain("RTS-9.11", "2011-09-15", src="Forts", session = "MAIN")

В качестве параметров передается название базового инструмента, дата экспирации и источник данных src="Forts", кроме того, можно указать сессию session = "MAIN" - основную или session = "EVENING" - вечернюю. Для примера, построим улыбку волатильности. Замечу, что ожидаемая волатильность в данном примере считается биржей и точной методики расчета мне найти не удалось.

plot(cbind(rts$calls[,'Strike'], rts$calls[,'IV']), type='l', main="Implied Volatility", ylab="Volatility", xlab="Strike")

Графики открытого интереса:

plot(cbind(rts$calls[,'Strike'], rts$calls[,'OI']), type='l', col='blue', main="Open Interest", ylab="R", xlab="Strike")
lines(cbind(rts$puts[,'Strike'], rts$puts[,'OI']), type='l', col='red')

2. Стали доступны дневные исторические данные по опционам:

getSymbols("RTS-9.11M150911CA 185000", from="2011-03-21", src="Forts")
chartSeries(RTS911M150911CA185000)

Планы на будущее: добавить дивиденды, статистику и тиковые данные.

И снова о коинтеграции временных рядов

2011-09-01T10:57:00.000+04:00

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

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

Что такое стационарный временной ряд? Если совсем просто, то это временной ряд свойства которого не меняются во времени. Т.е. если мы возьмем некий его отрезок, посчитаем стандартные статистические характеристики, такие как математическое ожидание или дисперсию, мы должны получить одинаковые величины, в пределах погрешности, естественно.

Рассмотрим пример на языке программирования R:

data <- rnorm(50)
plot(data, type='l')

Стационарный временной ряд имеет порядок интегрирования 0. Обратное не верно, но нам не очень то и нужно. Возьмем теперь стационарный ряд X₁, X₂... X_n, и преобразуем его в следующий ряд: Y₁ = X₁, Y₂ = Y₁ + X₂, ... Y_n = Y_n-1 + X_n получившийся ряд имеет порядок интегрирования 1. Если вернуться к определению, то временной ряд имеет порядок интегрирования 1 если ряд его приращений имеет порядок интегрирования 0. Аналогично можно определить более высокие порядки интегрирования. Рассмотрим пример ряда, интегрированого с порядком 1:

data <- cumsum(rnorm(100))
plot(data, type='l')

Классическое случайное блуждание как раз является примером такого ряда. В некотором приближении, цены на акции тоже можно считать примером временного ряда с порядком интегрирования 1. Теперь мы можем переформулировать определение коинтеграции в применении для цен. Два временных ряда цен коинтегрированы если мы можем найти такую их линейную комбинацию, которая является стационарным временным рядом. Классическим примером коинтеграции является соотношение цен на золото GLD и индекс акций золотодобывающих компаний GDX. Оба доступны на yahoo.finance.com для скачивания, но мы воспользуемся R чтобы упростить себе жизнь

library(tseries)
library(quantmod)

#сначала скачаем данные
from <- "2006-05-23"
to <- "2008-07-23"
getSymbols("GLD", from=from, to=to)
getSymbols("GDX", from=from, to=to)

#Очень важный момент в исследованиях на коинтеграцию, выплаты дивидендов, ребалансировки индекса все должно быть учетно. В случае GLD все просто, золото оно и есть золото. В случае GDX следует взять adjusted данные.
data <- cbind(Cl(GLD), GDX$GDX.Adjusted)
data <- data[complete.cases(data)]

names(data) <- c("GLD", "GDX")

#проводим линейную регрессию для определения правильного соотношения
model <- lm(GLD~GDX+0, data)

#вычисляем разницу цен (спред)
spread <- data$GLD - coef(model)[1]*data$GDX

#проводим тест Дики-Фуллера на стационарность
adf.test(as.vector(spread), k=0)

Тест Дики-Фуллера это далеко не единственный тест на стационарность, но в данном примере мы используем именного его. Результат:

Augmented Dickey-Fuller Test

data:  as.vector(spread) 
Dickey-Fuller = -3.8556, Lag order = 0, p-value = 0.01622
alternative hypothesis: stationary

Видим, что. p-value мал, и отсюда можем утверждать, что действительно GLD и GDX коинтегрированы. Что подтверждается визуально на графике:

plot(as.vector(Cl(GLD)), col='blue', type='l', main='Коинтеграция GLD и GDX', ylab='price')
lines(as.vector(GDX$GDX.Adjusted*coef(model)[1]), col='red', type='l')

Красиво? Безусловно, если не знать о том, что в 2008, когда цены на нефть взлетели, на мир обрушился финансовый кризис, этот спред развалился и более не является стационарным. Изменим в нашей программе значение переменной to на "2011-09-01" и получим:

Augmented Dickey-Fuller Test

data:  as.vector(spread) 
Dickey-Fuller = -2.9022, Lag order = 0, p-value = 0.1964
alternative hypothesis: stationary

Большое значение p-value уже не позволяет нам говорить о стационарности спреда, а значит и о коинтеграции двух цен, и график только подтверждает это:

Тем не менее, никто не отменял возможности восстановления стационарности спреда по мере стабилизации на рынке в будущем. И никто не отменял более сложных комбинаций инструментов, например золото, золотодобывающие компании и нефть.

Сравнение временных рядов

2011-08-29T06:13:00.001+04:00

Огромное количество данных в data mining вообще и в финансах в частности приходит к нам в виде временных рядов. Это не удивительно, ведь очень часто нас интересуют какие-то события или показатели изменяющиеся во времени. При этом, огромный пласт классической математики веками создавался для работы с множествами чисел. В результате, одним из самых популярных подходов к работе с такого рода данными был отказ от оси времени как полноценной оси координат и переход к простым множествам данных. Действительно, все алгоритмы, что я до сих пор рассматривал в серии data mining с примерами на R работают с наборами данных и игнорируют время.

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

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

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

Рассмотрим пример:

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

mean <- 0
std <- 0.01 

res1 <- cumsum(rnorm(1000, mean=mean, sd=std))     
res2 <- cumsum(rnorm(1000, mean=mean, sd=std))     

plot(res1, type="l", col="red", ylim=range(res1, res2), main="Два случайных блуждания", ylab='Y')      
lines(res2, col="blue")

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

library(MASS)   

mean <- 0
std <- 0.01 
cor <- 0.8  

covmat <- (std^2) * matrix(c(1, cor, cor, 1), nrow=2)
res <- mvrnorm(1000, c(mean, mean), covmat)  

res1 <- cumsum(res[,1])
res2 <- cumsum(res[,2])   

plot(res1, type="l", col="red", ylim=range(res1, res2), main="Два случайных блуждания с корреляцией", ylab='Y')      
lines(res2, col="blue")

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

Какая здесь аналогия с финансами? Непосредственная, предположим что мы нашли две акции, изменения цен которых сильно коррелируют. Значит ли это, что если вдруг, в какой-то момент времени одна стала стоить существенно дешевле чем другая, то мы можем купить дешевую, продать дорогую и ждать когда цена сойдется? Нет, не значит. Цена может не сойтись никогда. Вернемся к нашим алкоголикам, предположим теперь, что они знают друг друга и на каждый шаг кричат: эй, дружище, ты где? Услышав такой крик, они оценивают направление и расстояние друга до друга и стараются шагнуть в нужном направлении.

library(xts)
mean <- 0
std <- 0.01
coint <- 0.02

x <- 0 
y <- 0
for (i in 1:1000) {
  px <- last(x)
  py <- last(y)  
  x <- c(x,  px + (coint * (py - px)) + rnorm(1, mean=mean, sd=std))  
  y <- c(y,  py + (coint * (px - py)) + rnorm(1, mean=mean, sd=std))  
}  

ylim <- range(x, y)  
plot(x, ylim=ylim, type="l", col="red", main="Два коинтегрированных случайных блуждания")  
lines(y, col="blue")

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

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

Куплю будущее, недорого

2011-08-25T06:41:00.000+04:00

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

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

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

А важны нам сейчас фьючерсы. О них и поговорим.

Фьючерс, от английского future-будущее, уже своим названием как-бы намекает, что это некий контракт касательно будущего. А именно, фьючерс это право и обязательство купить или продать некий базовый актив по заранее определенной цене в фиксированный момент времени в будущем. Фьючерс может быть выпущен на акции, индексы, нефть, валютные пары и много чего еще. Сегодня мы будем рассматривать только фьючерсы на акции.

Например, я хочу купить акции сбербанка на новый год, но я не знаю, сколько они будут стоить к новому году поэтому я хочу договориться о цене сейчас. Иду на биржу и говорю, за сколько вы мне продадите акции сбербанка на новый год? Именно такую возможность мне предоставляет фьючерс РТС на акции сбербанка. Когда-то, для меня стало неожиданным открытием, что согласно современной теории, цена на сбербанковский фьючерс не зависит от перспектив сбербанка. Цена на фьючерс зависит лишь от текущей стоимости акций и процентной ставки. Это объясняется довольно простой теорией арбитража.

Суть такая: предположим, что мы можем брать неограниченное количество денег в долг, а также отдавать деньги на временное хранение по одной и той же процентной ставке 10% годовых. Предположим, что акция сбербанка стоит 100 рублей. Сколько должен стоит фьючерс на акцию сбербанка с датой окончания (экспирации) ровно через год?

Предположим, что есть человек, желающий купить фьючерс за 115 рублей, тогда мы можем взять в долг 100 рублей, купить на них акцию сбербанка и, одновременно, продать желающему фьючерс за 115. Через год мы отдадим акцию сбербанка покупателю фьючерса, получим за нее 115 рублей, отдадим 100+10% = 110 рублей в банк и останемся с чистой и гарантированной прибылью 5 рублей. Заметьте, нам абсолютно не важно сколько будет стоить акция сбербанка через год. И для получения прибыли нам не понадобилось ни копейки собственных денег.

Предположим теперь, что есть человек, желающий продать фьючерс за 105 рублей и у нас есть акция сбербанка. Тогда мы можем продать эту акцию, положить деньги в банк и одновременно купить фьючерс. Через год мы заберем из банка 100+10% = 110 рублей и выполним обязательства по фьючерсу купив обратно акцию сбербанка за 105 рублей. Снова в итоге 5 рублей гарантированной прибыли. И заметьте, нам опять абсолютно не важно сколько будет стоит акция сбербанка через год.

На рынке много игроков как первого типа так и второго. Они создают противоположные силы. И эти силы толкают цену фьючерса к некоторому значению, в нашем случае 110 рублей, что равно цене базового актива + процентная ставка. Конечно, формула на самом деле несколько сложнее, ведь проценты могут капитализироваться, акция может приносить дивиденды а брокер может брать комиссию за проведение операции. Все это можно было бы учесть, но для дальнейшего рассмотрения, оно не очень интересно.

Ключевой момент, который тут надо запомнить, это то, что цена на фьючерс определяется лишь стоимостью акции сейчас и процентной ставкой. Что будет если, например, все знают, что компания через год разорится? Будет ли это влиять на стоимость её фьючерсов? Ответ: будет, но через цену на акции. Если все знают, что компания через год разорится, её акции упадут до нуля и вместе с ними упадут фьючерсы.

Перейдем от теории к практике. Я тут несколько раз упомянул некую гипотетическую процентную ставку. Исходя из вывода формулы можно предположить, что это процентная ставка, по которой крупные игроки, банки, могут занимать и давать в долг деньги. Попробуем, для начала, подтвердить теорию а потом и оценить чему равна эта процентная ставка на практике. Для запуска примера понадобится расширение к модулю quantmod для работы на Российском рынке.

library(rusquant)
from <- "2010-09-01"
#июньский фьючерс на акции сбербанка
future <- getSymbols("SRM1", src="Finam", period="hour", from=from, auto.assign=FALSE)
stock <- getSymbols("SBER", src="Finam", period="hour", from=from, auto.assign=FALSE)
#стоимость фьючерса равна стоимости 100 акций
data <- cbind(100*Cl(stock), Cl(future))
data <- data[complete.cases(data)]
names(data) <- c("Stock", "Future")
chartSeries(future)
addTA(data$Future-data$Stock, legend="Spread", col='blue', yrange=c(-50, 300))

На картинке, в верхней части цена на фьючерс, в средней объемы торгов, а в нижней части разница цен фьючерса и акций (спред).

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

Попробуем найти угол наклона прямой в правой части графика с помощью уже известной нам линейной регрессии

#отбираем линейную часть данных
data <- data['2011-05-03/']
#вычисляем спред
spread <- data$Future-data$Stock
#дней до экспирации
exptime <- as.double(difftime(last(index(spread)), index(spread), units="days"))

model <- cbind(exptime, spread)
names(model) <- c("ExpTime", "Spread")
lm(Spread~ExpTime-1, model)

Результат:
ExpTime
1.468

Т.е. за каждый день держания сбербанковского фьючерса вы неявно платите 1 рубль, 47 копеек. Что без рекапитализации, в годовом эквиваленте дает приблизительно 5.5%. Это и есть та самая процентная ставка.

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

Feature selection с примерами на R

2011-08-18T06:40:00.000+04:00

В серии постов о data mining я уже описал несколько различных методов и алгоритмов. Они все разные, но у них есть одна неизбежная общая черта: на вход подается некоторый набор "входных параметров" (в англоязычной литературе используется слово feature) и уже дальше начинаются все чудеса и преобразования. В случае с финансовыми данными, этими параметрами могут быть цены, объемы торгов, день недели, технические индикаторы и много чего еще. Но откуда берутся эти параметры и почему они именно такие? До сих пор, этот вопрос я не рассматривал.

Часто в различных обсуждениях проскакивает точка зрения, что в модель нужно напихать как можно больше различных параметров, а дальше мощные алгоритмы data mining сами разберутся, что им нужно, а что нет. На самом деле, это конечно чушь. Неограниченное увеличение количества входных параметров не только снижает общее понимание модели, но и приводит к оверифитингу. Кроме того, скорость работы некоторых моделей экспоненциально убывает с ростом количества параметров. Поэтому, слишком много параметров, это плохо. Слишком мало параметров, очевидно, тоже не хорошо. Ведь чем меньше параметров, тем меньше информации о системе мы имеем. Получается, что необходимо каким-то образом определять, когда параметров слишком много, а когда мало. Но даже этого не достаточно. Необходимо еще каким-то образом выбирать из имеющегося множества параметров те, которые наиболее нам важны. Те, которые несут в себе больше всего уникальной информации. Как это сделать? Об этом и поговорим.

Проблема отбора параметров решается в большой области знаний, под названием Dimensionality Reduction (уменьшение размерности). В свою очередь, эта область знаний делится на две части Feature Selection (выбор параметров) и Feature Extraction (преобразование параметров). Сегодня мы поговорим о Feature Selection.

Что же из себя представляет Feature Selection? Предположим, что у нас есть некие данные, скажем цены на все акции на бирже и мы хотим по ним спрогнозировать движение рынка в целом (ничего себе задачка...). Что с этим делать? Понятно, что у нас может быть сотни тысяч акций и никакой метод опорных векторов или нейронная сеть не переварит это, а если даже и переварит, то неизбежно произойдет оверфитинг. Стало быть, нам надо часть данных отбросить. Какие именно? Можно, конечно, воспользоваться знаниями предметной области и вручную отобрать значимые параметры, которые, как нам кажется, должны играть наибольшую роль. На самом деле, этот подход является одним из лучших, если в наличии есть специалист глубоко понимающий предметную область. В случае финансов, это проблематично. В других областях, вроде анализа ДНК или обработки текста ситуация аналогична. Стало быть, процесс отбора параметров следует автоматизировать.

Решим задачу на коленке

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

Очевидно, что по каждому из параметров значения не разделимы. Но при этом, если рассматривать связку (x,y), то данные очень легко разделить на две группы. Т.е. даже в системах где параметры по одному не несут никакой дополнительной информации, в совокупности они могут привести к хорошим результатам.

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

Таких методов можно придумать очень много. Все они классифицируются по двум характеристикам. Рассматриваем ли мы параметры по одному или группами и каким образом мы оцениваем удачность подборки параметров.

Классификация методов отбора

По количеству рассматриваемых параметров, методы делятся на:

1) Методы ранжирования параметров.

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

2) Методы работающие с множествами параметров.

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

По способу отбора алгоритмы делятся на:

1) Фильтры.

Название говорит за себя. Очень распространены в совокупности с ранжированием параметров. Идея состоит в том, чтобы отобрать интересующие нас параметры по какому-то критерию. Часто используются стандартные статистические тесты. Этот подход полностью независимы от модели, которая используется для классификации.

2) Обертки.

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

3) Встроенные методы.

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

Пример

Перейдем к примеру, думаю, сразу все станет понятней. Допустим мы хотим узнать какие внешние факторы влияют на поведение индекса РТС. Не секрет, что Россия сырьевая страна и, наверное, можно было бы предположить, что влияет на нас более всего цена на нефть и курс доллара. Но как это проверить? Для проверки, предлагается взять значения индекса РТС и некоторое количество потенциальных параметров: курс доллара, цена на нефть brent, цена на нефть light, цена на золото, американский индекс s&p500, европейский индекс ftse100, цена на серебро и просто некоторая случайная величина, для проверки на вшивость.

Итак, запускаем R, собираем данные (нам потребуется последняя версия модуля rusquant (как его можно установить написано тут). Т.к. все эти параметры торгуются на разных рынках в разных часовых поясах, придется немного поработать над тем, чтобы выровнять их значения по времени. Собственно, я взял все значения в 18 вечера по Москве.

library(rusquant)
"select.hours" <-
function(data, hour){
    return(data[format(index(data), format="%H")==hour])
}


from <- "2007-10-01"
hour <- "18"
period <- "hour"

RTSI <- select.hours(getSymbols("RTSI", src="Finam", from=from, auto.assign=FALSE, period=period), hour)
BRENT <- select.hours(getSymbols("ICE.BRN", src="Finam", from=from, auto.assign=FALSE, period=period), hour)
LIGHT <- select.hours(getSymbols("NYMEX.CL", src="Finam", from=from, auto.assign=FALSE, period=period), hour)
USDRUB <- select.hours(getSymbols("USDRUB", src="Finam", from=from, auto.assign=FALSE, period=period), hour)
GOLD <- select.hours(getSymbols("comex.GC", src="Finam", from=from, auto.assign=FALSE, period=period), hour)
SP500 <- select.hours(getSymbols("SP500", src="Finam", from=from, auto.assign=FALSE, period=period), hour)
FTSE <- select.hours(getSymbols("FTSE", src="Finam", from=from, auto.assign=FALSE, period=period), hour)
SILVER <- select.hours(getSymbols("comex.SI", src="Finam", from=from, auto.assign=FALSE, period=period), hour)

data <- cbind(Cl(RTSI), Cl(USDRUB), Cl(BRENT), Cl(LIGHT), Cl(GOLD), Cl(SP500), Cl(FTSE), Cl(SILVER))

data <- ROC(data[complete.cases(data)])
data <- data[complete.cases(data)][-1]

#добавим случайное значение
data <- cbind(data, rnorm(length(data$RTSI.Close), mean(data$RTSI.Close), sd(data$RTSI.Close)))

names(data) <- c("RTS", "SI", "BRENT", "LIGHT", "GOLD", "SP500", "FTSE", "SILV", "RND")

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

panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y, use="complete.obs"))
    txt <- format(c(r, 0.123456789), digits=digits)[1]
    txt <- paste(prefix, txt, sep="")
    if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex * r)
}

pairs(~RTS+SI+BRENT+LIGHT+GOLD+SP500+FTSE+RND, data=data, lower.panel=panel.smooth, upper.panel=panel.cor)

Видно, что все наши параметры, кроме случайной величины, сильно коррелируют между собой и с индексом РТС. Какие-то параметры коррелируют сильнее, какие-то слабее. В целом, картина кажется заслуживающей дальнейшего исследования. Множество алгоритмов feature selection реализовано в пакете FSelector. Попробуем, для начала, самый простой фильтр: посчитаем корреляцию между каждым из параметров и индексом РТС.

library(FSelector)
linear.correlation(RTS~SI+BRENT+LIGHT+GOLD+SP500+FTSE+SILV+RND, data=data)

Результат:

SI         0.54385755
BRENT      0.53846199
LIGHT      0.43151621
GOLD       0.26078314
SP500      0.61501880
FTSE       0.64300959
SILV       0.36046822
RND        0.06929931

В этом тесте наибольшую значимость показали европейский ftse100, американский s&p500, доллар и нефть. Вполне ожидаемо, что наименьшую значимость имеет случайное число. В принципе, на этом можно было бы и остановиться, выбрав 3-4 самых значимых параметров и передав их в модель. Но кто сказал, что это будет оптимальное решение? Другой метод сортировки параметров основан на энтропии. Все кто изучали физику помнят, что энтропия системы это мера беспорядка. В данном случае, мы говорим о информационной энтропии, которая, как ни странно, тоже является, своего рода, мерой беспорядка в системе.

information.gain(RTS~SI+BRENT+LIGHT+GOLD+SP500+FTSE+SILV+RND, data=data)

Результат:

SI         0.32154920
BRENT      0.21683641
LIGHT      0.09690079
GOLD       0.03985278
SP500      0.35216883
FTSE       0.38601392
SILV       0.07476571
RND        0.00000000

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

evaluator <- function(subset) {
    #k-fold cross validation
    k lt;- 5
    splits lt;- runif(nrow(data))
    results = sapply(1:k, function(i) {
      test.idx lt;- (splits >= (i - 1) / k) & (splits lt; i / k)
      train.idx lt;- !test.idx
      test lt;- data[test.idx, , drop=FALSE]
      train lt;- data[train.idx, , drop=FALSE]
      model lt;- glm(as.simple.formula(subset, "RTS"), data=train, family = gaussian)
      delta lt;- test$RTS - predict(model, test)
      error.rate = sum(delta*delta) / nrow(test)
      return(1 - error.rate)
    })
    return(mean(results))
}

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

exhaustive.search(names(data)[-1], evaluator)

Результат: SI,BRENT,FTSE Здесь у нас уже нет никаких цифр, нам выдается лучший из найденных наборов параметров. И т.к. мы решали методом полного перебора, то мы можем быть уверенными, что это вообще лучшее возможное решение для данной модели. Как видим, оно вполне соответствует нашим изначальным прогнозам. Индекс РТС более всего определяется значениями курса доллара, цены на нефть марки brent и европейским фондовым индексом. Честно говоря, я предполагал что здесь еще будет индекс S&P500, но его не оказалось. Возможно, отчасти потому, что наши биржи слишком сильно разнесены по времени. Возможно если бы мы взяли время отсечки не 18 часов, а 23 часа, то S&P500 играл бы существенно большую роль. Но вернемся к нашим алгоритмам. В данном случае, на входе было всего 8 параметров и нам пришлось ждать ответа несколько минут. Что если параметров было бы 100 или 1000 или еще больше? В таком случае мы можем воспользоваться алгоритмом градиентного спуска или одним из жадных алгоритмов:

> hill.climbing.search(names(data)[-1], evaluator)
[1] "SI"    "BRENT" "SP500" "FTSE"  "SILV" 
> forward.search(names(data)[-1], evaluator)
[1] "SI"    "BRENT" "GOLD"  "FTSE"  "SILV" 
> backward.search(names(data)[-1], evaluator)
[1] "SI"    "BRENT" "SP500" "FTSE"  "SILV"  "RND"

Результаты другие, тем не менее, они не далеки от оптимальных.

Конечно же это далеко не все возможные подходы, и далеко не все, что можно сделать. Здесь еще можно рассмотреть много интересных вопросов, и к теме feature selection я еще обязательно вернусь.

Исследуем Российский рынок ценных бумаг с помощью R

2011-08-08T22:43:00.001+04:00

В прошлых статьях, посвященных среде разработки R, я многократно упоминал о расширении quantmod, пожалуй, одном из самых полезных расширений для желающих строить математические модели финансовых рынков в R. Единственный недостаток этого расширения - отсутствие поддержки российских рынков, и этот недостаток я решил исправить. Для этого, было написано маленькое расширеньице, оно позволяет автоматически выкачивать данные для российских площадок с сайта Finam. По мере возможности, постараюсь добавить еще данные с сайта РТС.

Установить это расширение просто. Если у вас R последней версии (2.13 на момент написания текста) то достаточно запустить его и выполнить команду:

install.packages("rusquant", repos="http://R-Forge.R-project.org")

Если R более ранней версии, но вы используете Windows, то можно поступить так:

install.packages("rusquant", contriburl="http://r-forge.r-project.org/bin/windows/contrib/latest/")

И третий вариант установки. Можно зайти на http://r-forge.r-project.org/projects/rusquant/ и вручную скачать либо исходники либо готовую версию.

Вот и всё, после сообщения об успешной установке этого и всех зависимых пакетов, вы сможете приступить к работе.

Например, можно построить дневной график фьючерса на индекс РТС:

library(rusquant)
getSymbols("SPFB.RTS", from="2011-01-01", src="Finam")
chartSeries(SPFB.RTS)

Кроме того, в отличии от оригинального quantmod, здесь есть маленький бонус, а именно, нам доступны не только дневные но и внутри-дневные данные. Например, можно посмотреть часовик фьючерса на золото:

library(rusquant)
getSymbols("SPFB.GOLD", from=Sys.Date()-5, src="Finam", period="hour")
chartSeries(SPFB.GOLD)

Или, даже, 5 минутку (будьте, однако осторожны выбирая диапазон, очень легко накачать реально много данных)

library(rusquant)
getSymbols("SPFB.GOLD", from=Sys.Date()-5, src="Finam", period="5min")
chartSeries(SPFB.GOLD)

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

Всего доступны для загрузки 9 временных интервалов: 1min, 5min, 10min, 15min, 30min, hour, day, week, month. Доступные инструменты и их коды можно смотреть на сайте источника в поле имя контракта.

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

Topcoder. String compression.

2011-07-29T10:42:00.001+04:00

В прошлый раз, я писал о своем желании провести детальный разбор решений лучших участников топкодера с целью постичь их тайное дао и примкнуть к их рядам. Итак отчитываюсь.
Первым делом я решил разобрать решение победителя последнего онлайн тура Topcoder Open 2011, участника из Китая, с ником ACRush. Его решение было написано на C++ и состояло из более чем 1800 строк кода. Для того чтобы это переварить, я решил переписать его решение на Java.

Для тех кто не в курсе, расскажу условия задачи. На вход вам давалась некая строка, длинной до 100000 символов. Вам нужно было её сжать, представив в виде нескольких строк вида:

1: abc2a3sd4
2: qwelwesdfdf
3: erwe2rew4dfds
4: qweqwe2dferw

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

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

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

Более точно условие задачи можно прочитать тут.

В марафона часто возникают задачи, в которых собственный результат можно оценить. В таких задачах, методика перебора нескольких способов решения и выбора лучшего, очень распространена среди лидеров. И решение ACRush не стало исключением. Я бы выделил в его решении 5 основных подходов.

Итак, решение.

Важный момент о решении, который следовало бы сразу заметить:

Нам важно правильно расставить цифры. Если цифры расставлены правильно, то нет ничего сложного в том, чтобы определить оптимальную комбинацию букв. Достаточно просто развернуть сжатое представление, сравнить с исходной строкой и посчитать какая буква и сколько раз попадает в ту или иную позицию в сжатом представлении. После чего, для каждой позиции, выбрать наиболее часто встречаемую букву. Это, вроде бы, простое умозаключение пришло в мою голову лишь на третий или четвертый день. В решении ACRush оно использовалось повсеместно.

Первый подход. Жадный алгоритм.

Этот алгоритм дает минимальное решение для любых начальных данных. По сути, мы перебираем строки вида:

xxxxxxxxx1
xxxxxxx2
xxxxx333333
xxxxxxxxxx

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

Второй подход. Поиск в глубину.

Этот подход дает идеальное решение для случаев когда процент ошибок в исходной строке мал. Он основан на том, что мы ищем наиболее вероятных кандидатов исходя из частоты вхождений символов. Идея простая, если посмотреть на мой самый первый пример, то понятно, что в разжатой строке чаще всего будет встречаться подстрока номер 2: qwelwesdfdf
И это свойство можно применить для поиска строк, использовавшихся для генерации. В своем решении, я сначала искал подстроки по частотам одного символа. Но точность этого подхода стремительно падала с увеличением ошибки. Потом, я стал искать на основе частот парных комбинаций символов. ACRush пошел еще на шаг дальше и искал наиболее часто встречаемые тройки последовательных символов. Таким образом, на каждом шаге поиска в глубину он выбирал подстроки из исходной, в которых тройки последовательных символов (триплеты) встречаются наиболее часто. Таких подстрок он выбирал сотню штук и для каждой из них оценивал, какое количество очков она даст, если ей воспользоваться (исходя из предполагаемого процента ошибки). Выбирал ту, которая дает максимальное количество очков и шел на следующий цикл поиска в глубину.

Для коротких исходных строк ( < 500 символов) частоты триплетов не считались, просто перебирались все возможные подстроки и для них считалось максимально возможное количество очков.

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

Третий подход. Динамическое программирование.

Здесь, на вход нам подавались частичные решения найденные на предыдущих этапах. Частичные они были в том смысле, что мы могли найти, скажем, 2 или 3 предполагаемых строки из четырех. Поиск в глубину, во втором подходе, идет снизу вверх, т.е. из нашего примера мы сначала ищем строку № 2, выполняем замену, потом ищем строку № 4 и т.д. Слабым местом такого подхода является как раз замена. Ведь мы не знаем процент ошибки в данном конкретном месте строки. Можно, конечно, брать ошибку с запасом, но увеличение ошибки приводит к ложным срабатываниям. Таким образом лучше всего было бы, перебрать все возможные варианты расстановки и выбрать из них оптимальный.
Это и реализовал ACRush использовав динамическое программирование. Все это проводилось для сотни наилучших частичных решений найденных к этому моменту.

Четвертый и пятый подходы.

По моим ощущениям, 99% успеха являются результатом второго и третьего подхода. Однако, видимо для того чтобы совсем добить конкурентов, в решении было реализовано еще два подхода. Случайное блуждание и метод грубой силы.
Для случайного блуждания выбиралось наилучшее имеющееся решение и проводились всевозможные замены цифр в каждой из строк (я уже говорил, что расстановка цифр однозначно определяет максимальное количество очков).
Метод грубой силы выполнялся все оставшееся время решения. Заключался он в том, чтобы перебрать все возможные порядки строк и все возможные расстановки цифр, до тех пор пока время не кончится.

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

Немного вечерней философии...

2011-07-13T22:52:00.000+04:00

Сегодня для меня завершился очередной Topcoder Open, не буду скрывать, в глубине души я все-таки надеялся попасть в финал. Но, надежды надеждами, а результаты так себе. И вроде бы с одной стороны можно успокаивать себя прогрессом, два года назад я вообще не прошел в третий тур, год назад я был на 66-ом месте, а в этом году предварительно на 52-ом, но уж больно медленный этот прогресс. И особенно обидно осознавать, читая описания решений лидеров, что идея то решения была правильной. Как всегда, не хватило каких-то деталей реализации, каких-то мелочей, из которых, конечно же, в нашем деле состоит всё.

И вот сейчас я бы хотел полностью пересмотреть свой подход к решению подобных задач. Нутром чую, что что-то делаю не так. Конечно было бы здорово поучиться у лидеров лично, но где они? Кто из них пишет блог или читает лекции? Пока мне видится единственный вариант: детальный разбор лучших решений, вылизывание и оттачивание собственных, до тех пор пока они не окажутся столь же успешным как и у лидеров.

Я постараюсь, да да, знаю, я уже много чего наобещал за последнее время, но все же, я постараюсь разбирать решения и публиковать тут детали реализации. Зачем? Прежде всего чтобы убедиться, что я сам всё понял, ведь лучший способ разобраться в проблеме это рассказать её другим людям, правда?

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

Метод опорных векторов с примерами на R

2011-07-12T21:26:00.001+04:00

Друзья, после продолжительного затишься, я продолжаю тему "data mining с примерами на R". И эта статья посвящается методу опорных векторов (в английской литературе support vector machines или, сокращенно, svm). В интернете масса информации по этому методу и, как всегда, она изобилует сложными математическими формулами, за которыми теряется общая идея. Поэтому я решил написать эту статью без единой формулы, в конце концов, если кому-то нужны формулы, он всегда может обратиться в википедию, скриншот из которой вы можете увидеть слева.

Итак, в чем же состоит идея метода опорных векторов? Давайте, сначала рассмотрим очень простой случай. Предположим, что у нас есть множество точек на плоскости, часть которых относится к классу A, а другая часть к классу B и есть точки класс которых нужно определить. Задачу такого рода, по-научному можно определить как задачу классификации, причем, для её решения, нужен алгоритм обучения "с учителем". Метод опорных векторов как раз подходит именно для таких задач.

Рассмотрим, для начала, упрощенный случай. Пусть точки, принадлежащие разным классам, можно разделить с помощью прямой линии:

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

Какую из них выбрать? Интуитивно понятно, что нам бы хотелось прямую где-нибудь по центру. Такая вот прямая, очевидно, не является лучшим выбором:

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

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

В простейшем случае, именно такую прямую и находит метод опорных векторов. Как именно он её находит? Чтобы ответить на этот вопрос пришлось бы привести кучу математических выкладок, которые в итоге сводятся к некоторой системе уравнений, которая, в свою очередь решается методами квадратичного программирования.

А если всё не так гладко?

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

Очевидно, что красная точка справа является выбросом и было бы лучше её игнорировать. Что делать в таком случае? В этом случае на помощь к нам приходит, так называемый, алгоритм с мягким зазором (soft marging). Математически это выражается введением дополнительных параметров в систему уравнений, а по сути, мы просто назначаем некий штраф, за каждую точку, оказавшуюся на чужой стороне. Используя такой подход мы можем работать даже с линейно неразделимыми данными.

Естественно, прямыми на плоскости все дело не ограничивается. Вместо плоскости могут быть сколь-угодно многомерные пространства, а вместо прямых гиперплоскости в этих пространствах. Сам алгоритм, при этом не меняется.

А если все-таки нет прямой?

Если бы это было всё, метод опорных векторов не был бы таким мощным. Действительно, в чем радость разделять данные линейно? С этим отлично могла бы справиться и простая линейная регрессия. На практике, мы не знаем структуру данных и очень редко их можно разделить прямой линией или плоскостью. Что если у нас на входе, например, такая картинка?

Оказывается, что система уравнений составляющая метод опорных векторов содержит только скалярные произведения координат различных точек. И с точки зрения математики, можно заменить каждое скалярное произведение на некоторую функцию (ядро) отвечающую определенным требованиям. Замена скалярного произведения ядром позволяет перевести проблему в пространство с другими измерениями. И в этом пространстве данные могут оказаться разделяемыми. Для примера, предлагаю посмотреть ролик, в котором все наглядно показано.

Таким образом, метод опорных векторов позволяет работать не только с линейно-разделяемыми данными, но и с более общими случаями. Кроме того, на практике, в отличии от, например, нейронных сетей, метод опорных векторов меньше страдает от перетренировки(overfitting): ситуаций, в которых алгоритм выявляет какие-то частные закономерности в тренировочных данных, в результате чего хорошо работает с тренировочным набором, но показывает очень плохой результат на новых и неизвестных данных.

Пример:

Уже по сложившейся традиции, я приведу пример из области финансовых рынков. Идея напрашивается сама собой, раз уж метод опорных векторов так хорошо классифицирует данные и способен работать как черный ящик с очень большим объемом входных параметров, то давайте скормим ему на вход как можно больше технических индикаторов и посмотрим, как хорошо он нам предскажет изменение цены.
Для анализа, мы возьмем индекс S&P500, его значения могут быть загружены автоматически в программу на R с помощью модуля quantmod. А реализацию метода опорных векторов мы найдем в пакете e1071.

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

Сначала нам предстоит провести подготовительную работу с данными, а именно, преобразовать их в удобоваримую форму:

library(quantmod)
library(e1071)

#Будем исследовать индекс S&P500
getSymbols("^GSPC", from='2005-01-01')
data <- GSPC

#Накидаем в нашу модель побольше разных индикаторов
atr <- ATR(data)[,'atr']
ema7 <- EMA(Cl(data), n=7)[,1]
ema14 <- EMA(Cl(data), n=14)[,1]
ema50 <- EMA(Cl(data), n=50)[,1]
macd <- MACD(Cl(data))[,'macd']
macd.signal <- MACD(Cl(data))[,'signal']
rsi <- RSI(Cl(data))
cmo <- CMO(Cl(data))
dc.high <- DonchianChannel(Cl(data))[,'high']
dc.low <- DonchianChannel(Cl(data))[,'low']
adx <- ADX(data)[,'ADX']

#определим модель, здесь важно сдвинуть либо предсказываемые значения на день вперед либо предикторы на день назад
model <- specifyModel(Next(Delt(Cl(data))) ~ atr + ema7 + ema14 + ema50 + macd + macd.signal + rsi + cmo + dc.high + dc.low + adx)

#разделим данные на две части - тренировочные и тестовые
train.window <- c('2005-01-01','2009-01-01')
validate.window <- c('2009-01-01', '2011-06-29')

train <- na.omit(as.data.frame(modelData(model, data.window=train.window)))
validate <- na.omit(as.data.frame(modelData(model, data.window=validate.window)))

#мы хотим покупать или продавать, когда индекс изменяется как минимум на 0.5% в ту или иную сторону, остальные случаи нас не интересуют
signals<-function(x) {

	if(x >= -0.005 && x<=0.005) { result<-"none" } else
	if(x > 0.005) { result<-"buy" } else
	if(x < -0.005) { result<-"sell" }

	result
}

#классифицируем имеющиеся данные, не забыв при этом выбросить из них сами цены, чтобы алгоритм не подсматривал ;)
class<-sapply(train$Next.Delt.Cl.data,signals)
train<-cbind(train,class)
train$Next.Delt.Cl.data <- NULL

class<-sapply(validate$Next.Delt.Cl.data,signals)
validate<-cbind(validate,class)
validate$Next.Delt.Cl.data <- NULL

Теперь когда у нас готова модель данных мы можем запустить наш алгоритм для анализа. Тут надо заметить, что метод опорных векторов, как впрочем и любой алгоритм data mining имеет множество параметров. На данном шаге мы выберем некоторые значения параметров произвольно:

#подготовка данных закончилась, можно запускать SVM на тренировочных данных
sv<-svm(class~.,train,gamma=0.01,cost=5,kernel="radial")

#а теперь посмотрим, насколько точно наш алгоритм предсказывает тестовые данные
table(actual=validate$class,predicted=predict(sv,newdata=validate,type="class"))

Полученные результаты оказываются лишь немного лучше, чем бросание монетки:

predicted
actual  buy none sell
  buy   60  136    3
  none  39  227    3
  sell  52  102    4

Что делать? Попробуем оптимизировать параметры, в конце концов, в документации к svm английским по белому сказано "Parameters of SVM-models usually must be tuned to yield sensible results! " Попробуем подобрать более оптимальные значения. К счастью для нас, нет никакой необходимости делать это вручную, модуль e1071 содержит метод tune который сделает всю работу. Для этого он применит кросс-валидацию. В процессе кросс-валидации тренировочные данные разбиваются на две части, модель тренируется на одной части, а тестируется на другой и так многократно. В итоге, выбираются лучшие параметры и модель с лучшими параметрами тренируется снова уже на полном наборе тренировочных данных. Все просто. Приступим:

#попробуем подобрать оптимальные параметры для алгоритма методом кросс-валиадации, эта строчка может выполняться несколько часов перебирая все возможные 

комбинации параметров
tuned <- tune(svm, class~., data = train,ranges = list(gamma = c(0.0001,0.001,0.05,0.1,0.2,0.3), cost = c(1,5,10,20,50,100,120,130)),tunecontrol = tune.control(sampling = "cross"),cross=3)

tuned

Дальше мы просто копируем найденные оптимальные значения и запускаем наш алгоритм снова

#подготовка данных закончилась, можно запускать SVM на тренировочных данных
sv<-svm(class~.,train,gamma=0.05,cost=10,kernel="radial")

#а теперь посмотрим, насколько точно наш алгоритм предсказывает тестовые данные
table(actual=validate$class,predicted=predict(sv,newdata=validate,type="class"))

Результаты, конечно не идеальны, но очевидно, что стали намного лучше.

predicted
actual buy none sell
  buy   73  123    3
  none  46  215    8
  sell  48   99   11

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

Фильтр Калмана, руководство полного идиота

2011-06-16T06:44:00.003+04:00

Фильтру Калмана посвящено множество статей и книг, совсем недавно была опубликована статья на хабре. Однако, чтение подобных статей у меня всегда опасение: а не идиот ли я? К сожалению, разобраться в чем же собственно суть дела, за всеми этими многоэтажными индексами и закорючками довольно сложно. Да и нужно ли? Здесь я не буду в очередной раз расписывать алгоритм и его математические основы, об этом можно почитать в книжках, википедии или где угодно еще. Сегодня я постараюсь изложить практическую сторону вопроса.

Итак, что же такое фильтр Калмана. Предположим у вас есть некая величина X изменяющаяся во времени. У вас есть некий прибор, измеряющий величину X с некоторой погрешностью. Вы хотите оценить, чему на самом деле равна величина X. Например, в случае акций X_k может быть истинной ценой акции в момент времени k (она определяется доходами компании, курсом валют, ожиданиями инвесторов итп). У нас есть только один способ измерить эту величину - пойти на биржу и посмотреть котировки. Но котировки на бирже подвержены влиянию массы факторов, цены колеблются в широком диапазоне и измеренная величина Y_k может отличаться от X_x. Алгоритм Калмана позволяет нам оценить истинное значение X_k на основе истории изменения измеряемой величины Y_k. При этом, естественно, делается ряд предположений о связи между Y и X и характере изменения X.

Предположение номер 1

X_k = A_k * X_k-1 + B_k * U_k + W_k-1

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

Что такое управляющий сигнал? Вернемся к нашему примеру с акциями. Предположим, что ФРС США активно печатает деньги, например 600 млрд долларов, эти деньги льются широкой рекой на рынок и неизбежно двигают цены акций вверх. Это и есть управляющий сигнал.

Но вернемся к нашей формуле, нам потребуется сделать еще одно важное уточнение. Шум w должен иметь нормальное распределение с дисперсией dW. На практике, часто так оно и есть, здесь на руку нам играет центральная предельная теорема, ведь шум это, как правило, совокупность большого множества неизвестных факторов.

Предположение номер 2

Y_k = H_k * X_k + V_k

Т.е. наша измеряемая величина Y это линейная комбинация реальной величины X и погрешности измерения V. При этом, делается дополнительное предположение что погрешность измерения тоже имеет нормальное распределение со своей дисперсией dV. Очень часто на практике H = 1.

Что мы только что сделали? Мы построили модель, описывающую нашу систему, при этом, заметьте, пока никакого фильтра Калмана не упоминалось.

Что мы хотим сделать?

Мы хотим получить некое уравнение

X_k = F(X_{k-1, k-2...}, Y_{k, k-1...}, dW, dV)

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

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

X_k = K_k * Y_k + (1 - K_k) * X_k-1

Где K_k - это Калмановский коэффициент усиления. Именно эту величину, в процессе работы и определяет алгоритм фильтра Калмана. При этом фильтр Калмана находит оптимальную величину коэффициента усиления при заданных дисперсиях dW и dV. Чем лучше вы оцените эти дисперсии, тем лучше будет работать фильтр Калмана.

А теперь, представьте, что X и Y это не одномерные величины, а вектор. A, B и H - матрицы. Погрешность и шум V и W тоже векторы. Оказывается, что при этом, все уравнения будут работать точно так же, и алгоритм фильтра Калмана легко обобщается на многомерные величины.

А теперь, по традиции, пример на R

Фильтр Калмана можно найти в нескольких модулях R, для примера я буду использовать реализацию из модуля dlm. В данном примере мы применим фильтр Калмана к значениям индекса S&P500.

Предположим, что управляющего сигнала нет, это конечно, на самом деле, не так, но оценка управляющего сигнала это отдельная задача. Кроме того, отдельными задачами являются оценки параметров дисперсий dV и dW. Для нашего примера я просто подобрал некоторые значения эмпирическим путем. Если же подойти к этой задаче основательно, то их вполне можно оценить исходя из волатильности на рынке.

Итак, пример:

library(dlm)
library(quantmod)

getSymbols("^GSPC", from="2010-09-01")
y <- as.ts(Cl(GSPC))
x <- dlmFilter(y, dlmModPoly(1, dV=0.1, dW=0.01))

plot(y, type='l')
lines(x$f[-1], col = "red")
#Добавим EMA для сравнения
lines(EMA(y), col = "green")

Для сравнения, на график, зеленым цветом, было добавлено 10-дневное экспоненциальное скользящее среднее Получилась картинка:

Большое спасибо читателю по имени Kirill за найденную в коде примера ошибку.

Пока Дойче. Привет Одноклассники!

2011-05-20T06:54:00.001+04:00

Как некоторым уже стало известно, затишье в моем блоге последние пару недель было вызвано сменой работы. Да, отныне я более не сотрудник Deutsche Bank, а работаю в компании Одноклассники.ру Это, однако, отнюдь не значит что финансовые рынки мне более не интересны. Скорее наоборот, теперь я освобожден от ограничений, накладываемых контрактом инвестиционного банка и могу свободно спекулировать на бирже. Так что, рассчитываю в ближайшее время опробовать на практике множество идей, о которых пока мог только теоретизировать. Надеюсь так же, что и для графоманской деятельности вскоре найдется время и мой блог пополнится множеством новых замечательных статей.

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

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

Кластерный анализ с примерами на R

2011-05-03T06:37:00.003+04:00

Доброго дня, уважаемый читатель! Перед вами очередной опус из серии о data mining. В прошлый раз я рассказал о методе ближайших соседей. Сегодня, как логическое продолжение, поговорим о кластерном анализе или кластеризации. С устоявшейся терминологией тут проблемы, т.к. большинство публикаций на английском и приходится придумывать русский эквивалент английским терминам. Потому и я иногда буду тоже скатываться на англицизмы.

Почему это логическое продолжение? Потому что идеи лежащие в основе этих подходов очень похожи. Напомню, что суть метода ближайших соседей состоит в том, что для каждого объекта мы ищем ближайших к нему соседей и на основании имеющихся данных о соседях делаем вывод об исходном объекте, на языке data mining, это называется обучением с учителем. Для того, чтобы этот подход работал, нужно иметь набор тренировочных данных.

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

Предварительная работа

Рассмотрим, к примеру, набор данных:

Имя	Возраст	Город	Доход
Петя	25	Москва	120000
Вася	34	Киров	45000
Маша	27	Самара	15000
Света	36	Нью-Йорк	150000
Катя	24	Москва	30000

Каким образом можно разбить эти данные на кластеры? Можно по зарплате: Петя и Света в один кластер, Вася и Катя во второй, а Маша в третий. Можно по возрасту: Петя, Маша и Катя в один, Вася и Света во второй. Можно еще разбить по месту жительства или полу. И всякий раз результат будет получаться разный. Какой результат лучше? Ответить можно только зная задачу, для чего именно проводится анализ. Данные сами по себе никак не определяют возможные группировки. Чтобы получить однозначный результат нам необходимо ввести понятие расстояния. Причем, мы можем использовать сразу несколько характеристик объектов для определения расстояния между ними, например зарплату, возраст и пол.

Существуют различные меры расстояний, но для нас интуитивно понятным является классическое Евклидово расстояние и, справедливости ради, замечу, что его в большинстве задач более чем достаточно. Важно лишь правильно нормализовать данные. Как в нашем примере: зарплаты измеряются в тысячах, а возраст в годах. Непосредственно сравнивать эти две величины нельзя. А вот если предварительно привести их к диапазону от 0 до 1 то они станут вполне сравнимыми.

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

Алгоритмы Кластерного Анализа

Известно множество алгоритмов кластерного анализа. Далеко не полная и не единственная возможная классификация представлена на рисунке ниже.

Итак, первым делом, все алгоритмы кластерного анализа делятся на две группы: иерархические и неиерархические. Первые строят не просто разбиение на классы, а иерархию разбиений. Результатом их работы, как правило, является дендрограмма, на основе которой, пользователь может сам выбрать желаемое разбиение. Неиерархические алгоритмы, напротив, в результате работы выдают некоторое конкретное разбиение и имеют ряд параметров, позволяющих настраивать алгоритм для имеющихся данных. Естественно, вторые работают быстрее чем первые.

Все методы кластеризации работают с данными в виде векторов в многомерном пространстве. Каждый вектор определяется значениями нескольких направлений, а направления и есть известные нам характеристики (пол, возраст, образование). Характеристики могут быть как количественные, так и качественные и искусство специалиста по data mining состоит в том, чтобы правильно отобрать и нормализовать эти характеристики, а потом выбрать подходящую меру расстояний. И только после этого в игру вступают алгоритмы кластеризации.

К числу наиболее популярных, неиерархических алгоритмов относится алгоритм k-средних. Он особенно популярен в силу простоты реализации и скорости работы. Главным его недостатком является сходимость к локальному минимуму и зависимость результата от начального распределения. Кроме того, требуется заранее знать предполагаемое число кластеров k.

Итак, сам алгоритм:
1) Выбрать k случайных центров в пространстве, содержащем исходные данные
2) Приписать каждый объект из множества исходных данных кластеру исходя из того, какой центр к нему ближе
3) Пересчитать центры кластеров используя полученное распределение объектов
4) Если алгоритм не сошелся, то перейти к п. 2. Типичные критерии схождения алгоритма это либо среднеквадратичная ошибка, либо отсутствие перемещений объектов из кластера в кластер.

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

Иерархическая кластеризация

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

Таким образом, алгоритмы с полной связью склонны находить более компактные кластеры. Алгоритмы с одиночной связью, напротив, склонны выявлять сильно вытянутые и сложны формы. Но часто страдают из-за шумов.
На рисунке ниже приведен классический пример работы алгоритма одиночной и полной связи. Нам даны два кластера, соединенных "дорожкой шумов". Слева результат работы алгоритма с одиночной связью, а справа алгоритма с полной связью:

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

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

Алгоритмы основанные на теории графов

На практике, в чистом виде, применяются редко, т.к. их вычислительная сложность не позволяет обрабатывать большие графы. Чаще используются в сочетании с другими алгоритмами, такими как k-средних. Классический пример, это алгоритм минимального покрывающего дерева. Идея состоит в том, чтобы представить весь набор данных в виде графа, вершины которого это наши элементы данных, а вес каждого ребра равен расстоянию между соответствующими элементами (помните, здесь тоже нет никаких чудес, сначала придется определить понятие расстояния). Тогда, мы можем построить минимальное покрывающее дерево на данном графе, а затем последовательно удалять ребра с наибольшим весом. Кластером считается множество элементов, соединенных "остатком" дерева. С каждым убранным ребром, количество кластеров увеличивается.

Статистические алгоритмы

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

Пример:

Для примера будет использоваться программная среда R. Что это такое, где взять и как настроить можно почитать тут.

Как всегда, пример у меня из области финансовых рынков. Допустим, у нас есть акции большого количества разных компаний. Одни из них растут, другие падают и мы накопили большую историю изменения их цен. Можем ли мы на основе этих данных объединить акции в группы? Логично было бы предположить, что акции из одного сектора рынка растут или падают вместе. И было бы логично если полученные группы это как-то отразили.

Для начала, я выкачал котировки примерно 20 российских компаний из разных областей, за два года с периодом 1 день. Для удобства, все они собраны в одном файле и доступны здесь. Попробуем проанализировать их средствами кластерного анализа. Т.к. цены разных акций отличаются очень сильно, и абсолютная величина цены нам совершенно не интересна, а интересно относительное изменение, приведем все цены к общему знаменателю: а именно, вместо самой цены будем считать логарифмическое изменение цены:
X_i = log(C_i) - log(C_i-1) ,где C_i - цена закрытия в i-ый день

x <- read.table('C:\\prj\\CNSVN\\stats\\stocks.csv', sep=',', header=TRUE) #читаем файл с ценами 
x <- log(x) #логарифмируем цены
x <- x[-1] #отбросим колонку с датой
x <- apply(x, 2, diff) #считаем разницу между последовательными элементами
x <- t(x) #транспонируем таблицу
kmeans(x, 5, 1000000) #будем разбивать на 5 кластеров, максимум 1000000 иттераций

Получаем результат:

1	2	3	4	5
RASP	AFLT, PMTL, PLZL, MMBM, MTSI, VZRZ	GAZP, LKOH, VTBR, URKA, TRNFP, SIBN, SBER, ROSN, GMKN, SNGS	OGK1, OGK2, OGK5	KMAZ, AVAZ

Сразу заметна тенденция к группировке в кластеры компаний из одного сектора: автопроизводители попали в пятый кластер, нефтедобывающие компании и крупные банки в третий, все генерирующие компании в четвертый. Если запустить этот алгоритм снова, то результат может получится другой, но тенденция, тем не менее, сохранится. Почему было решено разбивать именно на 5 кластеров, а не 6 или 4? К сожалению, решение вопроса о количестве кластеров не алгоритмизируется, в конечном счете это всегда должен решать пользователь. Но с использованием иерархической кластеризации это решение можно отложить на потом:

x <- read.table('C:\\prj\\CNSVN\\stats\\stocks.csv', sep=',', header=TRUE) #читаем файл с ценами 
x <- log(x) #логарифмируем цены
x <- x[-1] #отбросим колонку с датой
x <- apply(x, 2, diff) #считаем разницу между последовательными элементами
x <- t(x) #транспонируем таблицу

hc <- hclust(dist(x))

plot(hc)

Вот за что я люблю R! Нам понадобилось заменить всего одну строчку в прежнем примере и получилось следующая картинка:

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

Метод поиска ближайшего соседа с примерами на R

2011-04-14T06:55:00.002+04:00

Продолжаю серию про методы и алгоритмы data mining. Предыдущая статья была посвящена линейной регрессии. Эта статья посвящается методу ближайшего соседа. Идея метода ближайшего соседа очень проста на интуитивном уровне. Давайте рассмотрим простой пример.

Простой пример метода ближайшего соседа

Рассмотрим зарплаты людей, в зависимости от их места жительства. Предположим, что вы живете в центре Москвы и что рядом с вами живут люди с зарплатой 150 тыс рублей. Тогда, можно предположить, что и ваша зарплата близка к этой величине. Если же вы проживаете в поселке Новотагилка Челябинской области, и зарплаты ваших соседей около 3 тыс рублей, то скорее всего и ваша зарплата близка к 3 тыс рублей. Естественно, что здесь нет 100% точной зависимости, кто-то может жить в центре Москвы на пенсию в 6 тыс, а кто-то в Новотагилке зарабатывать 50. Но если единственная информация, которой мы владеем о человеке является его место жительства, то мы можем воспользоваться приведенным подходом для оценки его доходов.

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

А в попугаях я гораздо длиннее

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

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

А если нет элемента?

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

Обычно эта проблема решается следующим образом, выбирается некоторое последовательное число котировок. Например 10, первые 9 из них объявляются предикторами, а последнее искомым значением. Получается что-то вроде движущегося окна. Например, если у нас есть 100 последовательных значений, мы можем сначала взять значения 1-9 для тренировки, 10 — как искомое. Потом значения 2-10 для тренировки, 11 — искомое. Естественно, что при таком разбиении часть данных будет перекрываться, но искомое значение, каждый раз новое.

Метод k ближайших соседей

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

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

Очевидный способ улучшить положение - учитывать мнение нескольких ближайших соседей вместо одного. Это улучшение приводит нас к другому методу: метод k ближайших соседей. Идея состоит в том, чтобы взять не просто ближайшего соседа, а k ближайших соседей и предсказывать значения искомой величины на основе взвешенной оценки.

Метод k ближайших соседей, так же, позволяет оценить точность прогноза. Если все k ближайших соседей имеют одно и то же значение, то шансы что проверяемый объект будет иметь такое же значение очень высоки. Если же среди этих объектов 40% имеют одно, а 60% другое значение, то вероятность правильного прогноза существенно снижается. Таким образом мы можем использовать распределение голосов между соседями для оценки точности прогноза.

Пример

Пример сегодня будет на основе... биржевых котировок. Что может быть жизненней? Не считать же нам пестики и тычинки ирисов – кстати, если серьезно, по ссылке хороший пример применения метода ближайшего соседа.

Итак, задача: предсказать движения индекса S&P500 по его изменениям за прошлые 4 дня. Изменения считаем от закрытия к закрытию. К счастью для нас, ничего кроме R и нескольких дополнительных модулей к нему не нужно:

#раскоментировать если библиотека knnflex или quanmod не установлена
#install.packages("knnflex")
#install.packages("quantmod")

library(knnflex)
library(quantmod)

#получаем данные по индексу S&P500 с 2007 года
getSymbols("^GSPC", from="2007-10-01")

#считаем доходность по дням от закрытия до закрытия
roc <- ROC(Cl(GSPC)) 

# записываем в каждый ряд x изменения за прошедшие 4 дня
x <- cbind(lag(roc), lag(roc, k=2), lag(roc, k=3), lag(roc, k=4))

# в у записываем то, что нам надо предсказать
y <- roc

#тренировочные данные - это все данные за исключением последних 30 дней
train <- 1:(length(roc)-30)

#проверочные данные это все остальные данные, т.е. данные за последние 30 дней.
test <- (1:(length(roc)))[-train]

#Это самая важная строчка! 
#Считаем расстояние между каждой парой точек, в нашем случае это просто Евклидово расстояние 
#в 4-х мерном пространстве 
kdist <- knn.dist(x)

#считаем методом k-ближайших соседей результат
#первый параметр – индексы тренировочных данных
#второй параметр – индексы тестовых данных
#третий параметр – вектор с результатами для тренировочных данных
#четвертый параметр – матрица расстояний посчитаная строчкой выше
#k=3 – используем 3 ближайших соседа
preds <- knn.predict(train,test,y,kdist, k=3)

#связываем результат с исходными данными
results <- cbind(preds, y[test])

#считаем доходность за месяц
last(cumprod(1+ results[,2] * ifelse(results[,1]< 0, -1, 1)))

Результат работы:

2011-04-12      1.024020

Т.е. доходность за 30 дней составила 2.4%.

Вот такой простой и интуитивно понятный метод позволяет иногда давать достаточно точные прогнозы.

Линейная регрессия с примерами на R

2011-04-12T06:51:00.005+04:00

Меня давно преследует идея пройтись по ключевым алгоритмам data mining, систематизировать имеющиеся знания и составить некий обзор с примерами. Собственно, этой статьей и ознаменую начало данной серии. Для примеров будет использоваться программная среда R, о которой я уже писал в одном из прошлых постов.

Линейная регрессия

Строго говоря, линейная регрессия не является алгоритмом data mining. Это один из методов пришедших из статистики. В статистике, под регрессией, обычно подразумевают прогнозирование, в той или иной форме. Существует множество различных типов регрессий, но в основе любого из них лежит одна и та же идея: построить модель, связывающую предсказываемое значение с исходными данными (предикторами), минимизируя ошибку.

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

На языке математики, мы можем описать нашу линейную модель в виде уравнения: Y = a+b*X, где X – это предиктор или исходные данные, а Y – это предсказываемая величина. А саму задачу переформулировать в виде: найти коэффициенты a и b минимизирующие величину ошибки.
Когда модель построена, а коэффициенты найдены мы можем использовать полученное уравнение для предсказания неизвестных нам значений.

Пример:

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

Цена Акции = K * Цена барреля нефти + Некая константа

Если мы можем определить константу и коэффициент, то мы можем по цене на нефть предсказывать цену акции. Если возникает отклонение, то можно купить спред и на этом заработать, в теории, естественно.

А теперь практика, возьмем R, и проведем линейную регрессию цен акций компании Роснефть по цене на нефть. Исходные данные можно взять здесь:

data <- read.csv('C:\\ rosn.txt', sep='\t') #читаем данные из файла
fit <- lm(data$ROSN ~ data$BRN ) #строим линейную регрессию
summary(fit) #печатаем результат

Результат:

На рисунке, интересующие нас коэффициенты, обозначены (1). Подставив их в уравнение регрессии, получим выражение: ROSN = 61.405 + 1.944*Brent Помимо значений, собственно, коэффициентов, R показывает нам величины ошибок, или стандартного отклонения, для каждого коэффициента. Но это не все... Например, нам может быть интересно, объясняют ли вообще хоть что-нибудь наши коэффициенты. Чтобы проверить это, мы, выдвигаем нулевую гипотезу, что, к примеру, коэффициент 61.405 является лишь результатом погрешности и его значением можно пренебречь. Для проверки такой гипотезы, используется t-критерий Стьюдента. Здесь R за нас делает всю работу, вычисляя как саму величину t так и степень значимости нашей гипотезы Pr(>|t|). Так, в нашем случае величина (2) 0.0036 означает что мы на 100*(1-0.0036) = 99.64% уверены в том, что свободный член в нашем выражении отличен от нуля. Далее мы можем проверить, насколько точно наша модель описывает данные. Для этого используются коэффициенты R² (3). Чем ближе величина этих значений к 1, тем лучше. 1 это идеальный результат, означающий, что модель на 100% описывает данные. И, наконец, последнее, что мы можем проверить, это то, насколько предсказываемая величина зависит от предикторов. Для этого выдвигается нулевая гипотеза, что предсказываемая величина вообще не зависит от предикторов. Для этой гипотезы определяется p-значение (4). В нашем случае, оно получилось равным 2.65*10^-8. Т.е. мы можем быть уверенны на 99.99999735%, что предсказываемая величина действительно зависит от предикторов. Обычно, имеет смысл смотреть на этот параметр в первую очередь, ведь он определяет, насколько вообще наша модель адекватна. Да, кстати, график чуть выше, как раз показывает наши данные и результат линейной регрессии.

Что если мои данные не описываются прямой линией?

Линейная регрессия с одним предиктором это простейший вариант регрессии, на практике возможен целый ряд усложнений позволяющих решать гораздо более сложные задачи: 1) добавление нескольких предикторов – дополнительные предикторы могут добавить в модель больше информации и, тем самым, улучшить результаты. Уравнение линейной регрессии, в таком случае может иметь следующий вид: Y = a + b1*X1 + … bn*Xn 2)предварительное преобразование предикторов - Возведение предикторов в степень либо извлечение корней является одним из способов увеличить возможности линейной регрессии. На практике, нередко приходится проводить большое количество различных преобразований с целью выяснить, какое из них даст наилучшую оценку. 3)в уравнение модели могут быть добавлены члены содержащие произведение предикторов 4)часто возникает задача предсказания вероятности возникновения события. Т.е. предсказываемая величина принимает лишь два значения 1 и 0 или да и нет. Для решения такого рода задач применяется логистическая регрессия. К счастью, все эти и множество других дополнений реализованы в R и готовы к использованию. Информацию о том, как ими пользоваться можно почерпнуть набрав

> ?lm

в командном интерфейсе R.

Зачем Goldman Sachs крысы или генетические алгоритмы и нейронные сети на службе в инвестиционных банках

2011-04-01T08:03:00.000+04:00

Алгоритмический трейдинг в наше время составляет более 80% сделок практически на всех биржах. Это значит, что среднестатистический трейдер это уже не безумный человек, бегающий с телефоном, и кричащий "покупаем - продаем", это тихо жужащий в углу сервер с толстым каналом в интернет и надежной системой охлаждения. И это уже давно ни для кого не секрет. Но так было до недавнего времени. Последние годы ситуация стремительно меняется...

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

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

Существуют отточенные методики тренировки: крыса сажается перед торговым терминалом, ей показываются текущие курсы, а перед ней ставится две кнопки: купить и продать. Когда крыса совершат удачую сделку, ей автоматически выдается кусочек еды. Когда крыса ошибается её бьет током. Поверьте, крыса, мотивированная ударами тока, это гораздо дешевле и эффективней чем трейдеры, мотивированные годовым бонусом.

Таким образом, вместо ненадежных искуственных нейронных сетей, банки стали использовать естественные нейронные сети в мозгу крыс. И это уже дает свои результаты. Прибыль Goldman Sachs в прошлом году превзошла 8 млрд долларов.

Но технологии, однако, не стоят на месте. Вы слышали про генетические алгоритмы? Генетические алгоритмы были придуманы природой и заимствованы человеком. Но нет ничего эффективней самой природы. Следующим шагом в алгоритмическом трейдинге стала генетическая селекция крыс. Лучших крыс-трейдеров принуждают спариваться и приносить потомство. Среди потомства снова отбираются лучшие трейдеры и спариваются и так до бесконечности. Таким образом уже через 5-6 поколений получаются крысы заточенные под тот или иной финансовый инструмент. Одни способны безупречно торговать индексом S&P500, другие фьючерсами на золото, третьи валютами на форексе.

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

Коллеги, не останемся же и мы равнодушными, почему Россия опять позади планеты всей, это при том, что по показателю ККНД (количество крыс на душу населения) мы лидируем? Пока еще не позно, начинаем собирать крыс и сажать их за наши компьютеры. Будущие Уоррен Баффеты сейчас мечутся по подворотням и подвалам, не упустите свой шанс!