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

蒙特卡洛模拟和回溯检验

考虑一个随机过程，它会产生不同的结果，研究员想知道不同结果发生的可能性，以帮助解决实际问题，蒙特卡洛模拟就是常用的概率分析方法。

在开发策略的过程中，通过优化参数得到“最优”模型是否意味着回溯检验的结束？答案毫无疑问是否定的，这仅仅是第一阶段的结束，接下来需要对交易记录和净值曲线进行深度分析，主要目的是评估模型的预测能力，我们必须回答一个问题：回溯检验的优秀表现能否在交易中重现？

蒙特卡洛模拟可以评估策略运行面临的潜在风险，不管如何小心谨慎地优化参数，都不可避免地挖掘了样本数据，一旦部署策略，真实回撤可能显著大于回溯检验的结果，如何能够判断潜在回撤的概率分布？蒙特卡洛模拟解决了这个问题。首先我们要准备交易纪录，即每一笔交易的pnl，然后进行重复抽样，既可以有放回抽样，也可以是不放回抽样，然后计算每一条模拟的净值曲线的样本统计量，如最大回撤和年化收益率，根据中心极限定理，年化收益率的抽样分布渐进服从正态分布。最后进行概率推断，如计算抽样分布的置信区间，这些信息有助于判断策略运行的潜在风险。

案例分析

library(tidyverse)
library(quantmod)
library(PerformanceAnalytics)
library(Quandl)

从Quandl下载标普500指数期货的历史价格，计算日收益率。假设使用买入并持有策略。这里把日收益率作为日pnl。

stock <- Quandl("CHRIS/CME_ES1/6", type = "xts", collapse = "daily")
rets <- stock/lag.xts(stock) - 1
colnames(rets) <- "returns"
head(rets)

## returns
## 1997-09-09 NA
## 1997-09-10 -0.020342612
## 1997-09-11 -0.007650273
## 1997-09-12 0.017621145
## 1997-09-15 -0.002164502
## 1997-09-16 0.027114967

笔者创建了一个自定义函数：monte_carlo_simulation，接受3个参数：

• pnl: 要求是xts对象，代表pnl，很容易扩展到普通的数值向量，之所以使用xts，是因为笔者想使用PerformanceAnalytics包计算常用业绩指标
• n：正整数，模拟次数，默认值为100
• replace: 布尔值，TRUE代表有放回抽样，FALSE代表不放回抽样，不放回抽样会导致相同的收益率，因为只是对交易顺序重新排列，只有最大回撤是不同的，有放回抽样是一种更加严格的检验，因为只抽取了部分交易，收益率和最大回撤都不同。

monte_carlo_simulation计算了5个常用的业绩指标：累积收益率，年化收益率，年化标准差，年化夏普比率和最大回撤。返回结果是一个列表，方便后续计算。

monte_carlo_simulation <- function(pnl, n = 100, replace = FALSE) {
# pnl: xts, trading results
# n: positive integer, number of simulations
# replace: logical, sampling with replacement or without replacement
if(!is.xts(pnl)) {
stop("pnl must be xts object")
}
ex_na <- na.omit(pnl)
dt <- index(ex_na)
x <- coredata(ex_na)
size <- length(x)
if(factorial(size) < n) {
stop("n exceeds limit of possible sampling size")
}
m <- replicate(n = n, sample(x, replace = replace))
m_ts <- xts(m, order.by = dt)
equity_curve <- xts(apply(m, 2, function(x) cumprod(1 + x) - 1), order.by = dt)
cum_rets <- Return.cumulative(m_ts)[1, ] # result is matrix, need to convert it to numeric vector
ann_rets <- Return.annualized(m_ts, scale = 252)[1, ] # change frequency if needed
ann_stds <- StdDev.annualized(m_ts, scale = 252)[1, ]
ann_sharpe <- SharpeRatio.annualized(m_ts, scale = 252)[1, ]
mdd <- maxDrawdown(m_ts)[1, ]
out <- list(equity_curve = equity_curve, cum_rets = cum_rets,
ann_rets = ann_rets, ann_stds = ann_stds,
ann_sharpe = ann_sharpe, mdd = mdd)
return(out)
}
sim_results <- monte_carlo_simulation(pnl = rets, n = 500, replace = TRUE)

## Warning in factorial(size): 'gammafn'里的值在范围外

这里进行了500次不放回抽样，出于示范目的，500次已经足够，在计算能力充足的情况下可以进行500-10000次模拟。部分研究员指出，1000次模拟的结果已经非常逼近正态分布。让我们绘制500条模拟的净值曲线的波动轨迹。

plot(sim_results$equity_curve, main = "Simulations of Equity Curve")

接下来看业绩指标的抽样分布，首先看年化收益率。

ar <- sim_results$ann_rets
ggplot(NULL, aes(x = ar, y = ..density..)) +
geom_histogram(bins = 30) +
geom_line(aes(x = ar, y = dnorm(ar, mean = mean(ar), sd = sd(ar))), col = "red",
size = 1.0) +
labs(title = "Sampling Distribution of Annualized Returns",
x = "annualized returns", y = "density") +
theme_bw()

年化标准差。

as <- sim_results$ann_stds
ggplot(NULL, aes(x = as, y = ..density..)) +
geom_histogram(bins = 30) +
geom_line(aes(x = as, y = dnorm(as, mean = mean(as), sd = sd(as))), col = "red",
size = 1.0) +
labs(title = "Sampling Distribution of Annualized Standard Deviations",
x = "annualized stdev", y = "density") +
theme_bw()

年化夏普比率。

sr <- sim_results$ann_sharpe
ggplot(NULL, aes(x = sr, y = ..density..)) +
geom_histogram(bins = 30) +
geom_line(aes(x = sr, y = dnorm(sr, mean = mean(sr), sd = sd(sr))), col = "red",
size = 1.0) +
labs(title = "Sampling Distribution of Annualized Sharpe Ratio",
x = "annualized sharpe ratio", y = "density") +
theme_bw()

最大回撤。

mdd <- sim_results$mdd
ggplot(NULL, aes(x = mdd, y = ..density..)) +
geom_histogram(bins = 30) +
geom_line(aes(x = mdd, y = dnorm(mdd, mean = mean(mdd), sd = sd(mdd))), col = "red",
size = 1.0) +
labs(title = "Sampling Distribution of Maximum Drawdown",
x = "max drawdown", y = "density") +
theme_bw()

除了最大回撤外，年化收益率，年化标准差和年化夏普比率都非常接近正态分布。最大回撤呈现温和右偏分布，有极大值的存在。

接下来计算样本统计量的置信区间。

sim_results[[1]] <- NULL
ci <- sapply(sim_results, function(x) quantile(x, probs = c(0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 1.0)))
stats <- data.frame(cum_rets = Return.cumulative(rets)[1, ],
ann_rets = Return.annualized(rets, 252)[1, ],
ann_stds = StdDev.annualized(rets, 252)[1, ],
ann_sharpe = SharpeRatio.annualized(rets, scale = 252)[1, ],
mdd = maxDrawdown(rets))
rownames(stats) <- "backtest"
rbind(stats, ci)

## cum_rets ann_rets ann_stds ann_sharpe mdd
## backtest 1.941649 0.05332191 0.1959935 0.2720595 0.5711340
## 50% 1.982086 0.05401410 0.1958034 0.2768144 0.5077284
## 60% 2.674193 0.06465941 0.1972487 0.3253984 0.5367205
## 70% 3.768695 0.07810906 0.1986608 0.3964761 0.5748254
## 80% 5.176337 0.09161906 0.2000416 0.4645394 0.6277824
## 90% 7.060499 0.10569998 0.2025783 0.5500235 0.6940598
## 95% 9.670987 0.12073953 0.2042243 0.6174111 0.7312092
## 99% 21.347121 0.16134220 0.2071967 0.8207452 0.8178700
## 100% 45.884982 0.20352417 0.2098272 1.0446744 0.8918892

结果显示，有95%的可能性年化收益率会小于等于13.1%，有95%的可能性最大回撤会小于等于71.6%，显著高于回溯结果的57%。蒙特卡洛模拟更多用于评估潜在风险，所以我们应该更多地关注最大回撤和类似的风险度量指标，一般情况下，最大回撤的95%置信上限不应该超过原始策略的2倍，否则认为策略无法通过检验，当然这仅仅是经验之谈，并不是客观准则。