【编者按】【作者简介】陈东阳，LTR首席基金经理，拥有多年量化交易经验，尤其擅长R语言，并精通外汇、商品期货和美股ETF等市场运行机制。坚信数学和
【作者简介】陈东阳，LTR首席基金经理，拥有多年量化交易经验，尤其擅长R语言，并精通外汇、商品期货和美股ETF等市场运行机制。坚信数学和编程是打败市场的唯一方式，擅长时间序列分析和统计建模，使用自动化模型进行分析和交易，定期测试模型的预测能力和平稳性，及时进行订单管理及风险管理。
本文是机器学习系列的开篇之作，简单来说，笔者想利用高阶统计模型预测价格涨跌，通过交易赚取财富。毫无疑问这是非常困难的事情，但也正因为有巨大的挑战性，才使得过程更加丰富精彩。我们会主要关注外汇和贵金属市场，因为经验之顾，笔者最熟悉这两类市场环境。
研究目的
预测货币对的趋势。
文献和思路
笔者阅读了几篇相关文献，最有启发的是《Stock Trend Prediction with Technical Indicators with SVM》。
作者的研究主题是预测股票价格的趋势变化，预测变量使用了技术指标，个人认为这是非常好的思路，比很多学术文献专注于宏观经济变量实际得多。预测价格趋势通常视为分类问题，所以要用分类模型解决，作者最终选择了支持向量机(SVM)。我们会尝试不同类型的模型，包括逻辑回归，决策树，支持向量机等，然后选出最优的结果。
金融市场高度有效，准确预测下一个交易日的涨跌幅非常困难，很多模型的预测误差(MAPE)都接近10%，这样的结果是无法接受的。预测趋势的变化相对简单，这里指的是未来k个交易日的收盘价变化(k是正整数)。部分文献指出，当响应变量是未来3-10个交易日的涨跌(上涨记为1，下跌记为0)时，预测精度可超过70%。我们还是从最简单的情形开始，即预测未来一个交易日的涨跌概率。
模型设定
响应变量
响应变量是未来一个交易日的波动(更准确来说，是上涨优势比的自然对数)，上涨记为1，下跌记为0。计算时使用开盘价而不是收盘价，这更符合实际情况。首先尝试经典的逻辑回归。
预测变量
预测变量全部由技术指标构成，总共包含17个变量，分属3个类别：趋势，震荡和波动性。它们从不同方面描述了价格波动。指标名称后面的数字代表参数(回溯期)，没有写数值的指标使用默认参数，如MACD，计算时使用12(快速均线)，26(慢速均线)，9(均线差异的平滑窗口)的组合。
整理数据
我们从ducascopy下载了2010年1月至2017年6月的及时价格，然后将频率转化成日，转化时把时区从世界协调时转换成北京时间。ducascopy是最高质量的免费数据源。
library(tidyverse)
library(lubridate)
library(quantmod)
library(TTR)
选择EURUSD作为研究对象。
get_data <- function(file) {
# file: string, file path
# return a dataframe contains OHLC
if(!file.exists(file)) {
stop("file doesn't exist.")
}
df <- read_csv(file, col_types = "c_dddd_",
col_names = c("date", "open", "high", "low", "close"),
skip = 1) %>%
mutate(date = ymd(date))
return(df)
}
eurusd <- get_data("../daily forex/EURUSD.csv")
head(eurusd)
## # A tibble: 6 x 5
## date open high low close
## <date> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-01 1.43283 1.43354 1.43181 1.43337
## 2 2010-01-02 1.43337 1.43356 1.43225 1.43335
## 3 2010-01-04 1.43024 1.44556 1.42559 1.44369
## 4 2010-01-05 1.44366 1.44834 1.43864 1.44145
## 5 2010-01-06 1.44130 1.44237 1.42807 1.43865
## 6 2010-01-07 1.43869 1.44432 1.42976 1.43263
计算技术指标和响应变量。首先计算开盘价的一阶差分，然后转化成分类变量，若open(t+1) > open(t)，记为1，代表上涨，若open(t+1) < open(t)，记为0，代表下跌。注意，必须把结果变量滞后一期，因为我们用当前的数据预测未来涨跌。剔除周末的数据。如果用中国大陆的服务器，周六会有几个小时的交易，这部分低波动行情应该剔除。
结果存储在df数据框。
cal_variables <- function(ohlc) {
# calculate predictive variables and response variable
# ohlc: dataframe that have standard OHLC data
# return: a dataframe
# xts and TTR provide better mechanism to manipulate time series data
ts <- timetk::tk_xts(ohlc, select = open:close, date_var = date)
if(!is.OHLC(ts)) {
stop("df must contain open/high/low/close.")
}
# technical indicators
ema3 <- EMA(Cl(ts), n = 3)
ema10 <- EMA(Cl(ts), n = 10)
ema100 <- EMA(Cl(ts), n = 100)
macd <- MACD(Cl(ts), nFast = 12, nSlow = 26, nSig = 9, maType = "EMA")
macd_hist <- macd$macd
macd_diff <- macd_hist - macd$signal
adx20 <- ADX(HLC(ts), n = 20, maType = "EMA")
pdimdi_diff <- adx20$DIp - adx20$DIn
adx20 <- adx20$ADX
sto <- stoch(Cl(ts), nFastK = 14, nFastD = 3, nSlowD = 3, maType = "EMA")
sto_fastd <- sto$fastD
sto_diff <- sto$fastD - sto$slowD
cci3 <- CCI(Cl(ts), n = 3, maType = "EMA")
cci10 <- CCI(Cl(ts), n = 10, maType = "EMA")
mom3 <- momentum(Cl(ts), n = 3)
mom10 <- momentum(Cl(ts), n = 10)
atr3 <- ATR(HLC(ts), n = 3, maType = "EMA")$atr
atr10 <- ATR(HLC(ts), n = 10, maType = "EMA")$atr
bbands3 <- BBands(HLC(ts), n = 3, maType = "EMA", sd = 2)
bbands3_diff <- bbands3$up - bbands3$dn
bbands10 <- BBands(HLC(ts), n = 10, maType = "EMA", sd = 2)
bbands10_diff <- bbands10$up - bbands10$dn
op <- Op(ts)
# merge data, convert xts to dataframe, df is better at modeling
joined <- cbind(ema3, ema10, ema100, macd_hist, macd_diff, adx20, pdimdi_diff, sto_fastd,
sto_diff, cci3, cci10, mom3, mom10, atr3, atr10, bbands3_diff,
bbands10_diff, op)
colnames(joined) <- c("ema3", "ema10", "ema100", "macd_hist", "macd_diff", "adx20",
"pdimdi_diff", "sto_fastd", "sto_diff", "cci3", "cci10", "mom3",
"mom10", "atr3", "atr10", "bbands3_diff", "bbands10_diff", "op")
out <- joined %>%
timetk::tk_tbl(preserve_index = T, rename_index = "date") %>%
mutate(dop = c(NA, diff(op)),
dop_lag = lag.xts(dop, k = -1), # shift changes in open by 1 period
trend = factor(ifelse(dop_lag > 0, 1, 0))) %>%
filter(wday(date, label = F, week_start = 1) != 6)
return(out)
}
df <- cal_variables(eurusd)
查看变量和数据类型。
str(df)
## Classes 'tbl_df', 'tbl' and 'data.frame': 1956 obs. of 22 variables:
## $ date : Date, format: "2010-01-01" "2010-01-04" ...
## $ ema3 : num NA 1.44 1.44 1.44 1.44 ...
## $ ema10 : num NA NA NA NA NA ...
## $ ema100 : num NA NA NA NA NA NA NA NA NA NA ...
## $ macd_hist : num NA NA NA NA NA NA NA NA NA NA ...
## $ macd_diff : num NA NA NA NA NA NA NA NA NA NA ...
## $ adx20 : num NA NA NA NA NA NA NA NA NA NA ...
## $ pdimdi_diff : num NA NA NA NA NA NA NA NA NA NA ...
## $ sto_fastd : num NA NA NA NA NA NA NA NA NA NA ...
## $ sto_diff : num NA NA NA NA NA NA NA NA NA NA ...
## $ cci3 : num NA 100 36.69 -6.27 -53.44 ...
## $ cci10 : num NA NA NA NA NA ...
## $ mom3 : num NA NA 0.00808 0.0053 -0.01106 ...
## $ mom10 : num NA NA NA NA NA ...
## $ atr3 : num NA NA 0.0103 0.0123 0.0134 ...
## $ atr10 : num NA NA NA NA NA ...
## $ bbands3_diff : num NA 0.01 0.0159 0.0108 0.013 ...
## $ bbands10_diff: num NA NA NA NA NA ...
## $ op : num 1.43 1.43 1.44 1.44 1.44 ...
## $ dop : num NA -0.00313 0.01342 -0.00236 -0.00261 ...
## $ dop_lag : num 0.00054 0.01342 -0.00236 -0.00261 -0.00612 ...
## $ trend : Factor w/ 2 levels "0","1": 2 2 1 1 1 1 2 1 1 1 ...
划分训练集和检验集，两者比例约为4：1。
训练集：2015年1月1号至2016年12月30号
检验集：2017年1月1号至2017年6月30号
training <- df %>%
filter(date >= ymd("2015-01-01") & date <= ymd("2016-12-30")) %>%
select(-date, -op, -dop, -dop_lag)
head(training)
## # A tibble: 6 x 18
## ema3 ema10 ema100 macd_hist macd_diff adx20 pdimdi_diff
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.211746 1.217811 1.257811 -0.6621183 -0.1014720 9.584787 -5.624275
## 2 1.207443 1.215144 1.256729 -0.7192700 -0.1268990 11.010923 -10.865014
## 3 1.198748 1.209016 1.254383 -0.8538427 -0.1798587 14.861747 -18.608190
## 4 1.195324 1.205904 1.253145 -0.9176341 -0.1949200 17.097913 -17.619765
## 5 1.187782 1.201237 1.251702 -1.0346916 -0.2495821 19.741232 -20.886877
## 6 1.183496 1.197232 1.250266 -1.1221991 -0.2696717 22.453051 -22.434944
## # ... with 11 more variables: sto_fastd <dbl>, sto_diff <dbl>, cci3 <dbl>,
## # cci10 <dbl>, mom3 <dbl>, mom10 <dbl>, atr3 <dbl>, atr10 <dbl>,
## # bbands3_diff <dbl>, bbands10_diff <dbl>, trend <fctr>
建立模型
fit <- glm(trend ~ ., data = training, family = binomial("logit"))
summary(fit)
##
## Call:
## glm(formula = trend ~ ., family = binomial("logit"), data = training)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.14777 -0.38034 -0.05511 0.34232 2.66700
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.985e+00 7.561e+00 1.321 0.186639
## ema3 -2.937e+02 1.544e+02 -1.903 0.057051 .
## ema10 2.865e+02 1.594e+02 1.798 0.072229 .
## ema100 1.904e+00 1.549e+01 0.123 0.902169
## macd_hist 1.172e+00 9.147e-01 1.281 0.200219
## macd_diff 2.850e-01 3.473e+00 0.082 0.934606
## adx20 -9.994e-03 2.159e-02 -0.463 0.643376
## pdimdi_diff -6.306e-02 3.105e-02 -2.031 0.042304 *
## sto_fastd -7.094e+00 1.567e+00 -4.528 5.94e-06 ***
## sto_diff -6.835e+00 4.297e+00 -1.591 0.111650
## cci3 5.439e-02 6.311e-03 8.620 < 2e-16 ***
## cci10 2.990e-02 7.580e-03 3.944 8.00e-05 ***
## mom3 -7.221e+01 3.239e+01 -2.230 0.025764 *
## mom10 7.573e+01 2.245e+01 3.373 0.000744 ***
## atr3 -1.790e+02 1.110e+02 -1.613 0.106780
## atr10 1.880e+02 1.426e+02 1.318 0.187519
## bbands3_diff -1.905e+01 2.155e+01 -0.884 0.376497
## bbands10_diff 6.236e+00 1.286e+01 0.485 0.627826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 723.52 on 521 degrees of freedom
## Residual deviance: 302.34 on 504 degrees of freedom
## AIC: 338.34
##
## Number of Fisher Scoring iterations: 6
简约模型
从单个变量的统计显著性来看，pdimdi_diff，sto_fastd, sto_diff, cci3, cci10, mom3, mom10是最强的预测因子。
尝试拟合简约模型。必须指出，这是最简单但可能是无效的变量选择方法。为了避免过度挖掘，笔者不会使用逐步回归。由于变量数目不是很多，也没有必要使用PCA来降维，接下来也会考虑正则化。
fit_reduced <- glm(trend ~ pdimdi_diff + sto_fastd + sto_diff + cci3 +
cci10 + mom3 + mom10, data = training,
family = binomial("logit"))
summary(fit_reduced)
##
## Call:
## glm(formula = trend ~ pdimdi_diff + sto_fastd + sto_diff + cci3 +
## cci10 + mom3 + mom10, family = binomial("logit"), data = training)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.09945 -0.42486 -0.06932 0.37159 2.98779
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.498e+00 6.890e-01 6.528 6.66e-11 ***
## pdimdi_diff -5.234e-02 2.055e-02 -2.546 0.01088 *
## sto_fastd -8.637e+00 1.353e+00 -6.386 1.71e-10 ***
## sto_diff -8.984e+00 3.154e+00 -2.849 0.00439 **
## cci3 5.102e-02 5.813e-03 8.777 < 2e-16 ***
## cci10 3.397e-02 6.855e-03 4.955 7.25e-07 ***
## mom3 -1.173e+02 2.413e+01 -4.859 1.18e-06 ***
## mom10 3.699e+01 1.490e+01 2.483 0.01304 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 723.52 on 521 degrees of freedom
## Residual deviance: 324.52 on 514 degrees of freedom
## AIC: 340.52
##
## Number of Fisher Scoring iterations: 6
选择5%的显著性水平，所有参数估计值仍然是统计显著的，除mom3和mom10外，其它参数估计值的符号和取值变化不大，但AIC却上升。所以仍然使用第一个模型，简约模型没有优势。
拟合优度
空偏差(Null deviance)是不使用任何预测变量而使用常数概率来解释响应变量的偏差，这里的常数概率是训练集中事件发生的频率，空偏差类似回归模型的总离差平方和(TSS)。残余偏差(Residual deviance)是没有被预测变量解释的变异，类似于回归模型的残差平方和(RSS)。
仿照线性回归模型的
的计算方法，计算逻辑回归的伪
。
1 - fit$deviance/fit$null.deviance
## [1] 0.582134
预测变量“解释了”58.2%的响应变量的变异。
空偏差和残余偏差的差值近似服从自由度为k(预测变量个数)的卡方分布，接下来计算差值的p值，看是否显著异于零。
null_deviance <- fit$null.deviance
resi_deviance <- fit$deviance
diff_deviance <- null_deviance - resi_deviance
diff_df <- fit$df.null - fit$df.residual
pchisq(diff_deviance, df = diff_df, lower.tail = F)
## [1] 6.833464e-79
p值接近零，说明空偏差和残余偏差的差值是显著异于零的，基于技术指标的逻辑回归模型的确比空模型更好。
检验模型的预测能力
先生成检验集，检验模型在样本外的预测结果，再计算混淆矩阵。阈值选择0.5。
testing <- df %>%
filter(date >= ymd("2017-01-01") & date <= ymd("2017-06-30")) %>%
select(-date, -op, -dop, -dop_lag)
prob <- predict(fit, newdata = testing, type = "response")
pred <- factor(ifelse(prob > 0.5, 1, 0))
cm <- table(testing$trend, pred, dnn = c("actual", "forecast"))
cm
## forecast
## actual 0 1
## 0 49 16
## 1 4 61
精度指标。
eval_accuracy <- function(cm) {
if(nrow(cm) != 2 | ncol(cm) != 2) {
stop("invalid input, confusion matrix must be 2*2 table")
}
tn <- cm[1, 1] # actual = 0, pred = 0
fp <- cm[1, 2] # actual = 0, pred = 1
fn <- cm[2, 1] # actual = 1, pred = 0
tp <- cm[2, 2] # actual = 1, pred = 1
acc <- (tp + tn) / sum(cm)
sensitivity <- tp / (tp + fn)
specificity <- tn / (tn + fp)
ppp <- tp / (tp + fp)
npp <- tn / (tn + fn)
out <- data.frame(indicator = c("accuracy", "sensitivity", "specificity",
"positive predictive power", "negative predictive power"),
value = c(round(acc, 2), round(sensitivity, 2), round(specificity, 2),
round(ppp, 2), round(npp, 2)))
return(out)
}
eval_accuracy(cm)
## indicator value
## 1 accuracy 0.85
## 2 sensitivity 0.94
## 3 specificity 0.75
## 4 positive predictive power 0.79
## 5 negative predictive power 0.92
最主要的精度衡量指标是accuracy，达到85%，结果令人惊讶(too good to be true ?)，很多同类模型的预测精度只能达到60-70%。接下来衡量预测模型的平稳程度。
评估预测模型的平稳性
在回溯检验交易策略时，平稳性是非常重要的概念。如果一个模型能够在不同的货币对和不同的时间周期取得相似的结果，模型被认为是平稳的。这里进行简单的多货币检验，除了几个流动性最强的直盘外，还选择了EURGBP和GBPJPY两个交叉盘，以及黄金/美元(CFD)，它们具有完全不同的属性。
currencys <- c("AUDUSD", "EURGBP", "EURUSD", "GBPJPY", "GBPUSD", "USDCAD", "USDJPY", "XAUUSD")
files <- paste("../daily forex/", currencys, ".csv", sep = "")
results <- vector("list", length(currencys))
for(i in seq_along(results)) {
tmp_data <- get_data(files[i])
tmp_df <- cal_variables(tmp_data)
tmp_training <- tmp_df %>%
filter(date >= ymd("2015-01-01") & date <= ymd("2016-12-30")) %>%
select(-date, -op, -dop, -dop_lag)
tmp_fit <- glm(trend ~ ., data = tmp_training, family = binomial("logit"))
tmp_testing <- tmp_df %>%
filter(date >= ymd("2017-01-01") & date <= ymd("2017-06-30")) %>%
select(-date, -op, -dop, -dop_lag)
tmp_prob <- predict(tmp_fit, newdata = tmp_testing, type = "response")
tmp_pred <- factor(ifelse(tmp_prob > 0.5, 1, 0))
tmp_cm <- table(tmp_testing$trend, tmp_pred, dnn = c("actual", "prediction"))
results[[i]] <- eval_accuracy(tmp_cm)
}
names(results) <- currencys
results_df <- bind_cols(results)
results_df %>%
select(-num_range(prefix = "indicator", range = 1:7)) %>%
rename(AUDUSD = value, EURGBP = value1, EURUSD = value2, GBPJPY = value3,
GBPUSD = value4, USDCAD = value5, USDJPY = value6, XAUUSD = value7)
## indicator AUDUSD EURGBP EURUSD GBPJPY GBPUSD USDCAD
## 1 accuracy 0.85 0.85 0.85 0.88 0.85 0.89
## 2 sensitivity 0.94 0.80 0.94 0.90 0.92 0.84
## 3 specificity 0.77 0.89 0.75 0.85 0.78 0.94
## 4 positive predictive power 0.81 0.88 0.79 0.88 0.79 0.93
## 5 negative predictive power 0.92 0.82 0.92 0.88 0.91 0.86
## USDJPY XAUUSD
## 1 0.87 0.88
## 2 0.87 0.91
## 3 0.87 0.83
## 4 0.88 0.86
## 5 0.86 0.89
8个货币对的准确率都超过85%，初步探索表明，逻辑回归不仅得到非常高的预测精度，还具备平稳性，接下来需要优化模型，并构建交易策略和回溯检验。
欢迎加入【和高手一起玩赚EA】社群（群主QQ：3157002952，加好友时请注明暗号：66），在这里聚集了”外汇EA交易者、EA技术大神、EA爱好者、计算机宽客、金融极客、数据大师、草根经济学家”分享各种想法、算法、策略，共享思维碰撞带来的灵感。不仅有料，更有圈子！
本社群有哪些资源？
1）外汇资管业务，传统外汇资管业务，传统拨头皮EA+对冲EA+马丁EA
2）外汇培训课程+外汇线下活动+外汇线下展会。
3）EA编程，EA测试，DLL加密等MT4+MT5业务。
4）海外私募挂靠
5）策略资金对接
