Tuesday, November 29, 2011

Дивиденды rusquant 0.3.3

Выкатил очередное обновление 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')]


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

Monday, November 28, 2011

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

Последнее время немного выпал из реальности, погрузившись в изучение 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*x2+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 элементов.

Что делать?

Для начала, посмотрим на формулу, которой мы пытаемся описать наши данные:
yt = a0 + a1*x + a2*x + a3*x + a4*x + a5*x
Мы пытаемся подобрать эти коэффициенты ai таким образом, чтобы минимизировать ошибку т.е. минимизировать выражение:
L = Σ(yt-y)2
А что если в это выражение добавить еще какой-нибудь член, чтобы как-то уменьшить величины коэффициентов ai? Например, вместо L можно было бы минимизировать выражение:

L1 = Σ(yt-y)2 + λ*Σ|ai|

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

L2 = Σ(yt-y)2 + λ*Σ(ai)2

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

Первый вариант носит название 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 получается обычная линейная регрессия, которую я описал выше. При очень больших значениях λ получается горизонтальная линия, т.е. решение при котором коэффициенты ai=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.