评估数据源是否回溯_IAI Trade：用“决策树模型”评估信号好坏

IAI Trade致力于降低量化交易门槛，在IAI Trade用户可以使用“可视化策略生成器”：0代码生成EA策略，并一键接通模拟交易及真实交易。

研究目的

交易策略由三个部分组成，指标，信号和规则，指标可以是简单的技术指标，也可以是复杂的统计模型，但不管怎么样，最终必然会生成进场和离场的信号，这也是今天研究的重点。

交易信号的结果只有两种：盈利和亏损。在什么情况下某信号导致成功的交易？什么情况下会出现亏损？我们想知道哪些预期之外的因素影响了信号的结果，本文的研究目的是评估一些普遍使用的因子，例如趋势，震荡，波动性和时间对信号的影响。如果发现某些因子能够决定信号的好坏，就可以设计相应的筛子机制，嵌入到交易策略中剔除噪音，提升长期业绩。

笔者把它当成一个分类问题处理，响应变量是信号的结果，用1代表盈利，0代表亏损。预测变量则会选择一系列技术指标。可供选择的建模方法很多，我们想看到预测变量和响应变量的关系，决策树模型可能是最合适的，此外也可以尝试随机森林，它可以评估变量的相对重要性。

首先笔者构建了一个简单的交易策略：价格突破100天高点做多，回落至100天低点平仓；价格创100天低点做空，反弹创100天高点平仓。不使用任何的风控机制和资金管理。原始策略中也没有包含任何筛子机制，只有这样才能得到可信的结果。

准备数据

这个试验需要两方面的数据，第一是回溯检验的交易记录，第二是历史价格。上述策略在美元/日元的5分钟图进行检验，样本是2012年1月1号至2017年12月30号。

library(tidyverse)
library(lubridate)

回溯检验的交易记录。这份报告从MT4平台直接导出，关键变量是：

• time: 进场时间，字符串，精确到分钟
• order_type: 订单类型，字符串，代表买卖行为，只有四个可能的取值：buy, sell, close, close at stop
• pnl: 交易盈亏，数值型变量，单位是美元

trades <- read_csv("trades.csv", col_types = "icciddiidd")
head(trades)

## # A tibble: 6 x 10
## obs time order_type order_num lots entry stoploss
## <int> <chr> <chr> <int> <dbl> <dbl> <int>
## 1 1 2012.01.03 15:23 sell 1 0.01 76.747 0
## 2 2 2012.01.04 19:37 close 1 0.01 76.730 0
## 3 3 2012.01.04 19:37 buy 2 0.01 76.730 0
## 4 4 2012.01.06 15:23 close 2 0.01 77.125 0
## 5 5 2012.01.06 15:23 sell 3 0.01 77.125 0
## 6 6 2012.01.06 21:30 close 3 0.01 77.281 0
## # ... with 3 more variables: takeprofit <int>, pnl <dbl>, balance <dbl>

提取交易结果，按进场时间进行排列。

entry_time <- trades %>%
filter(order_type %in% c("buy", "sell")) %>%
select(entry_time = time)
outcome <- trades %>%
filter(order_type %in% c("close", "close at stop")) %>%
mutate(result = factor(ifelse(pnl > 0, 1, 0))) %>%
select(result)
results <- cbind(entry_time, outcome)
head(results)

## entry_time result
## 1 2012.01.03 15:23 1
## 2 2012.01.04 19:37 1
## 3 2012.01.06 15:23 0
## 4 2012.01.06 21:30 0
## 5 2012.01.07 02:54 1
## 6 2012.01.10 02:09 0

把进场时间转换成K线开盘时间，方便接下来跟包含预测变量的数据框进行合并。回溯检验时用tickData，所以进场时间可以精确到秒(并没有体现在trades中)，但价格数据的时间戳却是按K线的开盘时间排序的，所以需要对前者进行转换。例如第一笔交易在2012年1月3号15点23分进场，它对应的K线是15点20分开盘的K线。

bar_open_time <- function(x) {
minutes <- str_sub(x, 15, 16)
minutes_left <- str_sub(minutes, 1, 1)
minutes_right <- as.numeric(str_sub(minutes, 2, 2))
minutes_right_adjusted <- ifelse(minutes_right >= 0 & minutes_right <= 4, 0, 5)
minutes_adjusted <- str_c(minutes_left, as.character(minutes_right_adjusted))
out <- ymd_hm(str_c(str_sub(x, 1, 14), minutes_adjusted))
return(out)
}
results <- mutate(results, datetime = bar_open_time(entry_time))
head(results)

## entry_time result datetime
## 1 2012.01.03 15:23 1 2012-01-03 15:20:00
## 2 2012.01.04 19:37 1 2012-01-04 19:35:00
## 3 2012.01.06 15:23 0 2012-01-06 15:20:00
## 4 2012.01.06 21:30 0 2012-01-06 21:30:00
## 5 2012.01.07 02:54 1 2012-01-07 02:50:00
## 6 2012.01.10 02:09 0 2012-01-10 02:05:00

价格数据。数据源是ducascopy，这跟回溯检验的数据集是同一个。

usdjpy <- read_csv("usdjpy_m5.csv", col_types = "ccdddd_",
col_names = c("date", "time", "open", "high", "low", "close"),
progress = FALSE)
head(usdjpy)

## # A tibble: 6 x 6
## date time open high low close
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2003.05.05 05:00 118.940 118.972 118.940 118.946
## 2 2003.05.05 05:05 118.952 118.956 118.926 118.926
## 3 2003.05.05 05:10 118.929 119.016 118.929 118.951
## 4 2003.05.05 05:15 118.947 118.961 118.941 118.943
## 5 2003.05.05 05:20 118.944 118.960 118.933 118.944
## 6 2003.05.05 05:25 118.951 118.975 118.949 118.949

准备预测变量

我们用一系列技术指标代表价格波动的3个特征：趋势，震荡，波动性。

趋势指标：

• 短期趋势：5/20EMA
• 中期趋势：20/100EMA
• 长期趋势：100/500EMA

震荡指标：

• RSI(10), RSI(50), RSI(100)
• ROC(10), RSI(50), ROC(100)

波动性：

• ATR(10)，ATR(50)，ATR(100)

用双均线差异表现趋势是常用方法，为了具有普遍适用性，将双均线差异转化成分类变量，如果短期均线大于长期均线，用+1表示，代表上涨趋势，反之用-1表示，代表下跌趋势。

library(quantmod)
dt <- ymd_hm(paste(usdjpy$date, usdjpy$time, sep = " "))
ts <- xts(usdjpy[, c("open", "high", "low", "close")], order.by = dt)
ema5 <- EMA(Cl(ts), n = 5)
ema20 <- EMA(Cl(ts), n = 20)
ema100 <- EMA(Cl(ts), n = 100)
ema500 <- EMA(Cl(ts), n = 500)
rsi10 <- RSI(Cl(ts), n = 10, maType = "EMA")
rsi50 <- RSI(Cl(ts), n = 50, maType = "EMA")
rsi100 <- RSI(Cl(ts), n = 100, maType = "EMA")
roc10 <- ROC(Cl(ts), n = 10)
roc50 <- ROC(Cl(ts), n = 50)
roc100 <- ROC(Cl(ts), n = 100)
atr10 <- ATR(HLC(ts), n = 10, maType = "EMA")$atr
atr50 <- ATR(HLC(ts), n = 50, maType = "EMA")$atr
atr100 <- ATR(HLC(ts), n = 100, maType = "EMA")$atr
indicators <- cbind(ema5, ema20, ema100, ema500, rsi10, rsi50, rsi100, roc10, roc50,
roc100, atr10, atr50, atr100)
colnames(indicators) <- c("ema5", "ema20", "ema100", "ema500", "rsi10",
"rsi50", "rsi100", "roc10", "roc50", "roc100",
"atr10", "atr50", "atr100")
indicators_df <- indicators %>%
timetk::tk_tbl(preserve_index = T, rename_index = "datetime") %>%
mutate(short_trend = factor(ifelse(ema5 - ema20 > 0, 1, -1)),
medium_trend = factor(ifelse(ema20 - ema100 > 0, 1, -1)),
long_trend = factor(ifelse(ema100 - ema500 > 0, 1, -1))) %>%
select(-starts_with("ema"))
tail(indicators_df)

## # A tibble: 6 x 13
## datetime rsi10 rsi50 rsi100 roc10
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2018-01-31 12:30:00 26.44211 46.92540 49.48928 -0.0001929447
## 2 2018-01-31 12:35:00 50.46841 50.15935 51.00424 0.0001929110
## 3 2018-01-31 12:40:00 36.85777 47.37992 49.60036 -0.0002113126
## 4 2018-01-31 12:45:00 32.01000 46.15682 48.96857 -0.0006063502
## 5 2018-01-31 12:50:00 44.15566 48.04543 49.85027 -0.0002847864
## 6 2018-01-31 12:55:00 49.86659 49.01655 50.30727 -0.0002112776
## # ... with 8 more variables: roc50 <dbl>, roc100 <dbl>, atr10 <dbl>,
## # atr50 <dbl>, atr100 <dbl>, short_trend <fctr>, medium_trend <fctr>,
## # long_trend <fctr>

合并数据集

将results(包含响应变量)和indicators_df(包含预测变量)两个数据框合并，合并需要用到datetime变量，注意避免前视偏误，在15点23分进场交易时只能得到15点15分这根K线对应的指标值，所以要先把results的datetime减去5分钟。此外还会新增两个预测变量：工作日和时间，代表每周和日内的季节性效应。

合并数据后，时间序列就转变成横截面数据，时间戳不再是必须的，可以剔除。此外合并数据集中包含21个缺失值，由于比例太小而不会影响结果，直接删除。

joined <- results %>%
mutate(weekday = wday(datetime, label = FALSE, week_start = 1)) %>%
mutate(hour = hour(datetime)) %>%
mutate(datetime = datetime - period(5, units = "minute")) %>%
left_join(indicators_df, by = "datetime") %>%
select(-entry_time, -datetime) %>%
na.omit()
head(joined)

## result weekday hour rsi10 rsi50 rsi100 roc10
## 1 1 2 15 56.74986 42.42667 41.60538 2.605218e-04
## 2 1 3 19 72.02399 57.19172 53.48414 5.085973e-04
## 3 0 5 15 18.42304 39.76033 46.32272 -9.587356e-04
## 4 0 5 21 64.38439 53.61241 51.13853 8.039003e-04
## 5 1 6 2 28.62550 42.57676 45.02211 -6.489335e-05
## 6 0 2 2 68.60085 54.53831 52.14809 3.381762e-04
## roc50 roc100 atr10 atr50 atr100 short_trend
## 1 -9.373170e-04 -0.0016007187 0.01478793 0.01228483 0.01202077 -1
## 2 5.738207e-04 0.0004433780 0.02508215 0.02052653 0.01896610 1
## 3 -3.240000e-04 0.0004537764 0.02195081 0.01814913 0.01785781 -1
## 4 2.592218e-05 -0.0009974934 0.02431289 0.02307899 0.02243724 1
## 5 -1.374838e-03 0.0002855808 0.01146411 0.01989305 0.02377123 -1
## 6 6.764668e-04 0.0001040420 0.01869754 0.02111537 0.02104921 1
## medium_trend long_trend
## 1 -1 -1
## 2 -1 -1
## 3 1 1
## 4 -1 1
## 5 -1 1
## 6 1 -1

决策树模型

决策树模型将预测变量和响应变量的关系映射为一系列if-then-else的规则。我们会使用rpart包的rpart函数拟合模型。

随机划分训练集和检验集，两者比例为4：1。

set.seed(1234)
sample_vector <- sample(nrow(joined), size = as.integer(nrow(joined) * 0.8),
replace = FALSE)
training <- joined[sample_vector, ]
testing <- joined[-sample_vector, ]

拟合模型。

library(rpart)
fit <- rpart(result ~ ., data = training, method = "class",
control = rpart.control(minsplit = 2, minbucket = 1, cp = 0.001))

初始模型的结果非常复杂，这里不打印结果，它往往会过度挖掘数据，需要进行剪枝，先查看复杂度参数表。CP即复杂度参数，nsplit是分裂节点总和，rel error是训练集中各种树对应的误差，xerror是基于训练集的10折交叉验证误差，xstd是10折交叉验证误差的标准差。令人失望的是，这里根本不存在建模的必要，当nsplit = 0(不进行分裂)时10折交叉验证误差是最小的。

fit$cptable

## CP nsplit rel error xerror xstd
## 1 0.006961507 0 1.00000000 1.000000 0.02815957
## 2 0.006142506 3 0.97911548 1.045455 0.02843075
## 3 0.004914005 6 0.96068796 1.045455 0.02843075
## 4 0.004299754 7 0.95577396 1.047912 0.02844444
## 5 0.003685504 38 0.77641278 1.067568 0.02855036
## 6 0.003071253 57 0.69778870 1.121622 0.02880953
## 7 0.002457002 82 0.60687961 1.157248 0.02895506
## 8 0.002149877 127 0.48157248 1.191646 0.02907685
## 9 0.002047502 131 0.47297297 1.208845 0.02913091
## 10 0.001842752 153 0.41769042 1.224816 0.02917708
## 11 0.001719902 203 0.31818182 1.227273 0.02918384
## 12 0.001638002 211 0.30343980 1.235872 0.02920677
## 13 0.001228501 229 0.27395577 1.242015 0.02922246
## 14 0.001000000 376 0.08968059 1.261671 0.02926885

复杂度参数表也证实了这一点，我们倾向于使用低于虚线的点对应的复杂度参数。

plotcp(fit)

将复杂度参数修改为0.007，拟合剪枝模型。

library(rpart.plot)
fit_pruned <- rpart(result ~ ., data = training, method = "class",
control = rpart.control(cp = 0.007))
summary(fit_pruned)

## Call:
## rpart(formula = result ~ ., data = training, method = "class",
## control = rpart.control(cp = 0.007))
## n= 2296
##
## CP nsplit rel error xerror xstd
## 1 0.006961507 0 1 0 0
##
## Node number 1: 2296 observations
## predicted class=0 expected loss=0.3545296 P(node) =1
## class counts: 1482 814
## probabilities: 0.645 0.355

结论

这是一个失败的实验。上述模型无法有效区分成功和失败的信号，可能是预测变量的选择太糟糕。在有监督机器学习领域，好的预测变量比复杂的算法更加重要。这个实验暂时告一段落，接下来我们继续用分类模型研究价格涨跌。