Машинное обучение в трейдинге: теория, модели, практика и алготорговля - страница 213

 

Еще один пример на симуляцию.

Строится 20 000 линейных моделей (везде 1000 наблюдений, количество предикторов от 1 до 20 (по 1000 модели на каждое количество), плюс одна независимая переменная). Данные i.i.d., N(0,1).

 Цель симуляции - убедиться, что при построении МНК регрессии на независимых данных (не содержащих зависимостей), удовлетворяющих требованиям лин.модели, F-статистика не превышает критическое значение. Значит, может использоваться как индикатор обученности модели.

 

############### simulate lm f-stats with random vars


rm(list=ls());gc()


library(data.table)

library(ggplot2)


start <- Sys.time()


set.seed(1)


x <- as.data.table(matrix(rnorm(21000000, 0, 1), ncol = 21))

x[, sampling:= sample(1000, nrow(x), replace = T)]


lm_models <- x[, 

 {

  lapply(c(1:20), function(x) summary(lm(data = .SD[, c(1:x, 21), with = F], formula = V21 ~ . -1))$'fstatistic'[[1]])

 }

 , by = sampling

]


lm_models_melted <- melt(lm_models, measure.vars = paste0('V', c(1:20)))


crtitical_f_stats <- qf(p = 0.99, df1 = c(1:20), df2 = 1000, lower.tail = TRUE, log.p = FALSE)


boxplot(data = lm_models_melted, value ~ variable); lines(crtitical_f_stats, type = 's', col = 'red')


Sys.time() - start


gc()

Время работы кода: 1.35 минуты. 

 

Полезный код. Визуализируются последовательности сделок в трех ипостасях.

 

 

##########################


rm(list=ls());gc()


library(data.table)

library(ggplot2)

library(gridExtra)

library(tseries)


start <- Sys.time()


set.seed(1)


x <- as.data.table(matrix(rnorm(1000000, 0.1, 1), ncol = 1)) #random normal value with positive expectation

x[, variable:= rep(1:1000, times = 1000)]

x[, trade:= 1:.N, by = variable]


x.cast = dcast.data.table(x, variable ~ trade, value.var = 'V1', fun.aggregate = sum)

x_cum <- x.cast[, as.list(cumsum(unlist(.SD))), by = variable]


monte_trades <- melt(x_cum, measure.vars = names(x_cum)[-1], variable.name = "trade", value.name = 'V1')

setorder(monte_trades, variable, trade)

monte_trades_last <- as.data.table(monte_trades[trade == '1000', V1])

quantile_trade <- monte_trades[, quantile(V1, probs = 0.05), by = trade]

RF_last <- monte_trades[, V1[.N] / maxdrawdown(V1)[[1]], by = variable]



p1 <- ggplot(data = monte_trades, aes(x = trade, y = V1, group = variable)) +

geom_line(size = 2, color = 'blue', alpha = 0.01) + 

geom_line(data = quantile_trade, aes(x = trade, y = V1, group = 1), size = 2, alpha = 0.5, colour = 'blue') +

ggtitle('Simulated Trade Sequences of Length 1000')


p2 <- ggplot(data = monte_trades_last, aes(V1)) +

geom_density(alpha = 0.1, size = 1, color = 'blue', fill = 'blue') + 

scale_x_continuous(limits = c(min(monte_trades$V1), max(monte_trades$V1))) +

coord_flip() +

ggtitle('Cumulative Profit Density')


p3 <- ggplot(data = RF_last, aes(V1)) +

geom_density(alpha = 0.1, size = 1, color = 'blue', fill = 'blue') + 

geom_vline(xintercept = mean(RF_last$V1), colour = "blue", linetype = 2, size = 1) +

geom_vline(xintercept = median(RF_last$V1), colour = "red", linetype = 2, size = 1) +

ggtitle('Recovery Factor Density + Mean (blue) and Median (red)')


grid.arrange(p1, p2, p3, ncol = 3)


Sys.time() - start


gc()

 Работает примерно 45 сек. Отрисовывает примерно 1.5 минуты.

 
Alexey Burnakov:


Красиво, спасибо.
 
Alexey Burnakov:

Цель симуляции - убедиться, что при построении МНК регрессии на независимых данных (не содержащих зависимостей), удовлетворяющих требованиям лин.модели, F-статистика не превышает критическое значение. Значит, может использоваться как индикатор обученности модели.

Я про fstatistic не до конца понял. Данные тут случайные, но модель чему-то обучилась, значит можно сделать вывод о подгонке и переобученности модели. А значит оценка модели должна быть плохой. Т.е. я ожидал негативный fstatistiс, или какие-то другие признаки того что всё плохо на графике.
Как правильно интерпретировать результат из того примера?
Я так понял что первый предиктор можно считать качественне чем первый+второй. А 1+2 качественнее чем 1+2+3. Как-то так? Целесообразно ли подбирать генетикой набор предикторов, которые дадут наибольший fstatistic?
 
Dr.Trader:
Я про fstatistic не до конца понял. Данные тут случайные, но модель чему-то обучилась, значит можно сделать вывод о подгонке и переобученности модели. А значит оценка модели должна быть плохой. Т.е. я ожидал негативный fstatistiс, или какие-то другие признаки того что всё плохо на графике.
Как правильно интерпретировать результат из того примера?
Я так понял что первый предиктор можно считать качественне чем первый+второй. А 1+2 качественнее чем 1+2+3. Как-то так? Целесообразно ли подбирать генетикой набор предикторов, которые дадут наибольший fstatistic?

Посмотрите таблицу F распределения. http://www.socr.ucla.edu/applets.dir/f_table.html

 F-статистика - это величина, зависящая от степеней свободы. Она всегда положительна, так как имеем одностороннее распределение.

Но модель ничему не учится, так как обученная модель должна иметь высокую F-статистику (больше или равно critical при заданной alpha - как это звучит при проверке нулевой гипотезы).

Во всех случаях не превзошли критическое значение при alpha = 0.01, а можно задать 0.0001, например.

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

F-Distribution Tables
  • Ivo Dinov: www.SOCR.ucla.edu
  • www.socr.ucla.edu
Statistics Online Computational Resource
 
Alexey Burnakov:

Полезный код. Визуализируются последовательности сделок в трех ипостасях.

 

Относительно приведенного кода. Пожалуйста, пишите хотя бы краткие комментарии в коде. Особенно когда применяете сложные выражения. Не все знают и пользуют пакет "data.table", не лишне пояснить, что делает dcast.data.table, melt  что есть .N, .SD. Вы же публикуете код не для того, что бы показать насколько Вы глубоко в теме. По моему мнению публикуемый код должен помочь остальным пользователям (даже с начальным уровнем подготовки) разобраться в скрипте. 

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

Несколько предложений по коду:

- промежуточные переменные x, x.cast, x.cum не нужны при расчетах и только занимают память. Все  вычисления  не требующие сохранения промежуточных результатов желательно выполнять через pipe

Например

#---variant-------------
rm(list=ls());gc()
library(data.table)
library(ggplot2)
library(gridExtra)
library(tseries)
#----
require(magrittr)
require(dplyr)
start <- Sys.time()
monte_trades <- as.data.table(matrix(rnorm(1000000, 0.1, 1), ncol = 1)) %>%
        .[, variable := rep(1:1000, times = 1000)]%>%
        .[, trade := 1:.N, by = variable] %>%
        dcast.data.table(., variable ~ trade, value.var = 'V1', fun.aggregate = sum)%>%
        .[, as.list(cumsum(unlist(.SD))), by = variable]%>%
        melt(., measure.vars = names(.)[-1], variable.name = "trade", value.name = 'V1')%>%
        setorder(., variable, trade)
monte_trades_last <- as.data.table(monte_trades[trade == '1000', V1])
quantile_trade <- monte_trades[, quantile(V1, probs = 0.05), by = trade]
RF_last <- monte_trades[, V1[.N] / maxdrawdown(V1)[[1]], by = variable]
Sys.time() - start
#Time difference of 2.247022 secs

А строит графики конечно очень долго.

Не критика.

Удачи
 

 

 
Dr.Trader:
Я про fstatistic не до конца понял. Данные тут случайные, но модель чему-то обучилась, значит можно сделать вывод о подгонке и переобученности модели. А значит оценка модели должна быть плохой. Т.е. я ожидал негативный fstatistiс, или какие-то другие признаки того что всё плохо на графике.
Как правильно интерпретировать результат из того примера?
Я так понял что первый предиктор можно считать качественне чем первый+второй. А 1+2 качественнее чем 1+2+3. Как-то так? Целесообразно ли подбирать генетикой набор предикторов, которые дадут наибольший fstatistic?

А вот пример, когда мы подразумеваем, что полностью обученная модель будет включать 20 переменных, причем вес их возрастает (1 переменная - вес 1, 20ая переменная - вес 20). И посмотрим, как будет меняться распределение F статистики после добавления последовательно в модель предикторов:

############### simulate lm f-stats with non-random vars


rm(list=ls());gc()


library(data.table)

library(ggplot2)


start <- Sys.time()


set.seed(1)


x <- as.data.table(matrix(rnorm(20000000, 0, 1), ncol = 20))

x[, (paste0('coef', c(1:20))):= lapply(1:20, function(x) rnorm(.N, x, 1))]

x[, output:= Reduce(`+`, Map(function(x, y) (x * y), .SD[, (1:20), with = FALSE], .SD[, (21:40), with = FALSE])), .SDcols = c(1:40)]

x[, sampling:= sample(1000, nrow(x), replace = T)]


lm_models <- x[, 

 {

  lapply(c(1:20), function(x) summary(lm(data = .SD[, c(1:x, 41), with = F], formula = output ~ . -1))$'fstatistic'[[1]])

 }

 , by = sampling

]


lm_models_melted <- melt(lm_models, measure.vars = paste0('V', c(1:20)))


crtitical_f_stats <- qf(p = 0.99, df1 = c(1:20), df2 = 1000, lower.tail = TRUE, log.p = FALSE)


boxplot(data = lm_models_melted, value ~ variable, log = 'y'); lines(crtitical_f_stats, type = 's', col = 'red')


summary(lm(data = x[sample(1000000, 1000, replace = T), c(1:20, 41), with = F], formula = output ~ . -1))


Sys.time() - start


gc()

 График с логарифмической y-осью:

 

Видно, да...

> summary(lm(data = x[sample(1000000, 1000, replace = T), c(1:20, 41), with = F], formula = output ~ . -1))


Call:

lm(formula = output ~ . - 1, data = x[sample(1e+06, 1000, replace = T), 

    c(1:20, 41), with = F])


Residuals:

     Min       1Q   Median       3Q      Max 

-19.6146  -2.8252   0.0192   3.0659  15.8853 


Coefficients:

    Estimate Std. Error t value Pr(>|t|)    

V1    0.9528     0.1427   6.676  4.1e-11 ***

V2    1.7771     0.1382  12.859  < 2e-16 ***

V3    2.7344     0.1442  18.968  < 2e-16 ***

V4    4.0195     0.1419  28.325  < 2e-16 ***

V5    5.2817     0.1479  35.718  < 2e-16 ***

V6    6.2776     0.1509  41.594  < 2e-16 ***

V7    6.9771     0.1446  48.242  < 2e-16 ***

V8    7.9722     0.1469  54.260  < 2e-16 ***

V9    9.0349     0.1462  61.806  < 2e-16 ***

V10  10.1372     0.1496  67.766  < 2e-16 ***

V11  10.8783     0.1487  73.134  < 2e-16 ***

V12  11.9129     0.1446  82.386  < 2e-16 ***

V13  12.8079     0.1462  87.588  < 2e-16 ***

V14  14.2017     0.1487  95.490  < 2e-16 ***

V15  14.9080     0.1458 102.252  < 2e-16 ***

V16  15.9893     0.1428 111.958  < 2e-16 ***

V17  17.4997     0.1403 124.716  < 2e-16 ***

V18  17.8798     0.1448 123.470  < 2e-16 ***

V19  18.9317     0.1470 128.823  < 2e-16 ***

V20  20.1143     0.1466 137.191  < 2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 4.581 on 980 degrees of freedom

Multiple R-squared:  0.9932, Adjusted R-squared:  0.993 

F-statistic:  7123 on 20 and 980 DF,  p-value: < 2.2e-16

 
Vladimir Perervenko:

Спасибо! Это я еще не умел. Но действительно, надо по-максимуму держать расчеты в памяти. Будет быстрее. Хорошее кун-фу...

С учетом ваших доработок:

#---variant-------------

rm(list=ls());gc()

library(data.table)

library(ggplot2)

library(gridExtra)

library(tseries)

#----

require(magrittr)

require(dplyr)

start <- Sys.time()

monte_trades <- as.data.table(matrix(rnorm(1000000, 0.1, 1), ncol = 1)) %>%

.[, variable := rep(1:1000, times = 1000)]%>%

.[, trade := 1:.N, by = variable] %>%

dcast.data.table(., variable ~ trade, value.var = 'V1', fun.aggregate = sum)%>%

.[, as.list(cumsum(unlist(.SD))), by = variable]%>%

melt(., measure.vars = names(.)[-1], variable.name = "trade", value.name = 'V1')%>% 

setorder(., variable, trade)

monte_trades_last <- as.data.table(monte_trades[trade == '1000', V1])

quantile_trade <- monte_trades[, quantile(V1, probs = 0.05), by = trade]

RF_last <- monte_trades[, V1[.N] / maxdrawdown(V1)[[1]], by = variable]


p1 <- ggplot(data = monte_trades, aes(x = trade, y = V1, group = variable)) +

geom_line(size = 2, color = 'blue', alpha = 0.01) + 

geom_line(data = quantile_trade, aes(x = trade, y = V1, group = 1), size = 2, alpha = 0.5, colour = 'blue') +

ggtitle('Simulated Trade Sequences of Length 1000')


p2 <- ggplot(data = monte_trades_last, aes(V1)) +

geom_density(alpha = 0.1, size = 1, color = 'blue', fill = 'blue') + 

scale_x_continuous(limits = c(min(monte_trades$V1), max(monte_trades$V1))) +

coord_flip() +

ggtitle('Cumulative Profit Density')


p3 <- ggplot(data = RF_last, aes(V1)) +

geom_density(alpha = 0.1, size = 1, color = 'blue', fill = 'blue') + 

geom_vline(xintercept = mean(RF_last$V1), colour = "blue", linetype = 2, size = 1) +

geom_vline(xintercept = median(RF_last$V1), colour = "red", linetype = 2, size = 1) +

ggtitle('Recovery Factor Density + Mean (blue) and Median (red)')


grid.arrange(p1, p2, p3, ncol = 3)


Sys.time() - start

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

 
Alexey Burnakov:

Спасибо! Это я еще не умел. Но действительно, надо по-максимуму держать расчеты в памяти. Будет быстрее. Хорошее кун-фу...

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

У меня расчет занимает 

#-время выполнения в секундах

# min            lq             mean         median       uq           max       neval

# 2.027561   2.253354   2.254134   2.275785     2.300051 2.610649  100 

 Но это не так важно. Речь шла о читаемости кода.

Удачи 

ПС. И распараллеливайте вычисления lm(). Это как раз тот случай когда нужно 

 
Vladimir Perervenko:

У меня расчет занимает 

#-время выполнения в секундах

# min            lq             mean         median       uq           max       neval

# 2.027561   2.253354   2.254134   2.275785     2.300051 2.610649  100 

 Но это не так важно. Речь шла о читаемости кода.

Удачи 

ПС. И распараллеливайте вычисления lm(). Это как раз тот случай когда нужно 

Нее. Вы приводите время для части кода до графиков. Я же указал, вместе с графиками.

До графиков у меня: 1.5 секунды. Вашим методом - 1.15 секунды.

 

rm(list=ls());gc()


library(data.table)

library(ggplot2)

library(gridExtra)

library(tseries)


start <- Sys.time()


set.seed(1)


x <- as.data.table(matrix(rnorm(1000000, 0.1, 1), ncol = 1)) #random normal value with positive expectation

x[, variable:= rep(1:1000, times = 1000)]

x[, trade:= 1:.N, by = variable]


x.cast = dcast.data.table(x, variable ~ trade, value.var = 'V1', fun.aggregate = sum)

x_cum <- x.cast[, as.list(cumsum(unlist(.SD))), by = variable]


monte_trades <- melt(x_cum, measure.vars = names(x_cum)[-1], variable.name = "trade", value.name = 'V1')

setorder(monte_trades, variable, trade)

monte_trades_last <- as.data.table(monte_trades[trade == '1000', V1])

quantile_trade <- monte_trades[, quantile(V1, probs = 0.05), by = trade]

RF_last <- monte_trades[, V1[.N] / maxdrawdown(V1)[[1]], by = variable]


Sys.time() - start

 

rm(list=ls());gc()

library(data.table)

library(ggplot2)

library(gridExtra)

library(tseries)

#----

require(magrittr)

require(dplyr)


start <- Sys.time()


monte_trades <- as.data.table(matrix(rnorm(1000000, 0.1, 1), ncol = 1)) %>%

.[, variable := rep(1:1000, times = 1000)]%>%

.[, trade := 1:.N, by = variable] %>%

dcast.data.table(., variable ~ trade, value.var = 'V1', fun.aggregate = sum)%>%

.[, as.list(cumsum(unlist(.SD))), by = variable]%>%

melt(., measure.vars = names(.)[-1], variable.name = "trade", value.name = 'V1')%>% 

setorder(., variable, trade)

monte_trades_last <- as.data.table(monte_trades[trade == '1000', V1])

quantile_trade <- monte_trades[, quantile(V1, probs = 0.05), by = trade]

RF_last <- monte_trades[, V1[.N] / maxdrawdown(V1)[[1]], by = variable]


Sys.time() - start

 

Получается, все же побыстрее у вас... 

Причина обращения: