Comparing Flexible and Elastic Asset Allocation

So recently, I tried to combine Flexible and Elastic Asset Allocation. The operative word being–tried. Essentially, I saw Flexible Asset Allocation as an incomplete algorithm — namely that although it was an excellent method for selecting securities, that there had to have been a better way to weigh stocks than a naive equal-weight scheme.

It turns out, the methods outlined in Elastic Asset Allocation weren’t doing the trick (that is, a four month cumulative return raised to the return weight multiplied by the correlation to a daily-rebalanced equal-weight index of the selected securities with cumulative return greater than zero). Either I managed a marginally higher return at the cost of much higher volatility and protracted drawdown, or I maintained my Sharpe ratio at the cost of much lower returns. Thus, I scrapped all of it, which was a shame as I was hoping to be able to combine the two methodologies into one system that would extend the research I previously blogged on. Instead, after scrapping it, I decided to have a look as to why I was running into the issues I was.

In any case, here’s the quick demo I did.

require(quantmod)
require(PerformanceAnalytics)
require(IKTrading)

symbols <- c("VTSMX", "FDIVX", "VEIEX", "VBMFX", "VFISX", "VGSIX", "QRAAX")

getSymbols(symbols, from="1990-01-01")
prices <- list()
for(i in 1:length(symbols)) {
  prices[[i]] <- Ad(get(symbols[i]))  
}
prices <- do.call(cbind, prices)
colnames(prices) <- gsub("\\.[A-z]*", "", colnames(prices))
ep <- endpoints(prices, "months")
adPrices <- prices
prices <- prices[ep,]
prices <- prices["1997-03::"]
adPrices <- adPrices["1997-04::"]

eaaOffensive <- EAA(monthlyPrices = prices, returnWeights = TRUE, cashAsset = "VBMFX", bestN = 3)
eaaOffNoCrash <- EAA(monthlyPrices = prices, returnWeights = TRUE, cashAsset ="VBMFX", 
                     bestN = 3, enableCrashProtection = FALSE)
faa <- FAA(prices = adPrices, riskFreeName = "VBMFX", bestN = 3, returnWeights = TRUE, stepCorRank = TRUE)
faaNoStepwise <- FAA(prices = adPrices, riskFreeName = "VBMFX", bestN = 3, returnWeights = TRUE, stepCorRank = FALSE)

eaaOffDaily <- Return.portfolio(R = Return.calculate(adPrices), weights = eaaOffensive[[1]])
eaaOffNoCrash <- Return.portfolio(R = Return.calculate(adPrices), weights = eaaOffNoCrash[[1]])
charts.PerformanceSummary(cbind(faa[[2]], eaaDaily))

comparison <- cbind(eaaOffDaily, eaaOffNoCrash, faa[[2]], faaNoStepwise[[2]])
colnames(comparison) <- c("Offensive EAA", "Offensive EAA (no crash protection)", "FAA (stepwise)", "FAA (no stepwise)")
charts.PerformanceSummary(comparison)

rbind(table.AnnualizedReturns(comparison), maxDrawdown(comparison))

Essentially, I compared FAA with the stepwise correlation rank algorithm, without it, and the offensive EAA with and without crash protection. The results were disappointing.

Here are the equity curves:

In short, the best default FAA variant handily outperforms any of the EAA variants.

And here are the statistics:

                          Offensive EAA Offensive EAA (no crash protection) FAA (stepwise) FAA (no stepwise)
Annualized Return             0.1247000                           0.1305000      0.1380000          0.131400
Annualized Std Dev            0.1225000                           0.1446000      0.0967000          0.089500
Annualized Sharpe (Rf=0%)     1.0187000                           0.9021000      1.4271000          1.467900
Worst Drawdown                0.1581859                           0.2696754      0.1376124          0.130865

Note of warning: if you run EAA, it seems it’s unwise to do it without crash protection (aka decreasing your stake in everything but the cash/risk free asset by a proportion of the number of negative return securities). I didn’t include the defensive variant of EAA since that gives markedly lower returns.

Not that this should discredit EAA, but on a whole, I feel that there should probably be a way to beat the (usually) equal-weight weighting scheme (sometimes the cash asset gets a larger value due to a negative momentum asset making it into the top assets by virtue of the rank of its volatility and correlation, and ends up getting zeroed out) that FAA employs, and that treating FAA as an asset selection mechanism as opposed to a weighting mechanism may yield some value. However, I have not yet found it myself.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

A New Harry Long Strategy and A Couple of New PerfA Functions

So, Harry Long (of Houston) came out with a new strategy on SeekingAlpha involving some usual mix of SPXL (3x leveraged SPY), TMF (3x leveraged TLT), and some volatility indices (in this case, ZIV and TVIX). Now, since we’ve tread this path before, expectations are rightfully set. It’s a strategy that’s going to look good in the sample he used, it’s going to give some performance back in the crisis, and it’ll ultimately prove to be a simple-to-implement, simple-to-backtest strategy with its own set of ups and downs that does outperform the usual buy-and-hold indices.

Once again, a huge thanks to Mr. Helmuth Vollmeier for the long history volatility data.

So here’s the code (I’ll skip a lot of the comparing equity curves of my synthetic instruments to the Yahoo-finance variants, as you’ve all seen that before) to get to the initial equity curve comparison.

require(downloader)
require(PerformanceAnalytics)

download("https://dl.dropboxusercontent.com/s/950x55x7jtm9x2q/VXXlong.TXT",
         destfile="longVXX.txt")
VXX &amp;amp;lt;- xts(read.zoo("longVXX.txt", sep=",", header=TRUE))
TVIXrets &amp;amp;lt;- Return.calculate(Cl(VXX))*2
getSymbols("TVIX", from="1990-01-01")
TVIX &amp;amp;lt;- TVIX[-which(index(TVIX)=="2014-12-30"),] #trashy Yahoo data, removing obvious bad print
compare &amp;amp;lt;- merge(TVIXrets, Return.calculate(Ad(TVIX)), join='inner')
charts.PerformanceSummary(compare)
charts.PerformanceSummary(compare["2012::"])
charts.PerformanceSummary(compare["2013::"])
charts.PerformanceSummary(compare["2014::"])
charts.PerformanceSummary(compare["2015::"]) #okay we're good

download("https://www.dropbox.com/s/jk3ortdyru4sg4n/ZIVlong.TXT",
         destfile="longZIV.txt")
ZIV &amp;amp;lt;- xts(read.zoo("longZIV.txt", sep=",", header=TRUE))
ZIVrets &amp;amp;lt;- Return.calculate(Cl(ZIV))

getSymbols("SPY", from="1990-01-01")
SPXLrets &amp;amp;lt;- Return.calculate(Ad(SPY))*3

getSymbols("TMF", from="1990-01-01")
TMFrets &amp;amp;lt;- Return.calculate(Ad(TMF))
getSymbols("TLT", from="1990-01-01")
TLTrets &amp;amp;lt;- Return.calculate(Ad(TLT))
tmf3TLT &amp;amp;lt;- merge(TMFrets, 3*TLTrets, join='inner')
charts.PerformanceSummary(tmf3TLT)
discrepancy &amp;amp;lt;- as.numeric(Return.annualized(tmf3TLT[,2]-tmf3TLT[,1]))
tmf3TLT[,2] &amp;amp;lt;- tmf3TLT[,2] - ((1+discrepancy)^(1/252)-1)
charts.PerformanceSummary(tmf3TLT)
modifiedTLT &amp;amp;lt;- 3*TLTrets - ((1+discrepancy)^(1/252)-1)
TMFrets &amp;amp;lt;- modifiedTLT

components &amp;amp;lt;- cbind(SPXLrets, ZIVrets, TMFrets, TVIXrets)
components &amp;amp;lt;- components["2004-03-29::"]
stratRets &amp;amp;lt;- Return.portfolio(R = components, weights = c(.4, .2, .35, .05), rebalance_on="years")
charts.PerformanceSummary(stratRets)
SPYrets &amp;amp;lt;- Return.calculate(Ad(SPY))
compare &amp;amp;lt;- merge(stratRets, SPYrets, join='inner')
charts.PerformanceSummary(compare["2011::"])

With the following equity curve display:

So far, so good. Let’s look at the full backtest performance.

charts.PerformanceSummary(compare)

With the resultant equity curve:

Which, given what we’ve seen before, isn’t outside the realm of expectation.

For those interested in the log equity curves, here you go:

compare[is.na(compare)] &amp;amp;lt;- 0
plot(log(cumprod(1+compare)), legend.loc="topleft")

And for fun, let’s look at the outperformance equity curve.

diff &amp;amp;lt;- compare[,1] - compare[,2]
charts.PerformanceSummary(diff, main="relative performance")

And the result:

Now this is somewhere in the ballpark of what you’d love to see from your strategy against a benchmark — aside from a couple of spikes which do a number on the corresponding drawdowns chart, it looks like a steady outperformance.

However, the new features I’d like to introduce in this blog post are a quicker way of generating the usual statistics table I display, and a more in-depth drawdown analysis.

Here’s how:

rbind(table.AnnualizedReturns(compare), maxDrawdown(compare))

Which gives us the following result:

&amp;amp;gt; rbind(table.AnnualizedReturns(compare), maxDrawdown(compare))
                          portfolio.returns SPY.Adjusted
Annualized Return                 0.2181000    0.0799000
Annualized Std Dev                0.2159000    0.1982000
Annualized Sharpe (Rf=0%)         1.0100000    0.4030000
Worst Drawdown                    0.4326138    0.5518552

Since this saves me typing, I’ll be using this format from now on. And as a bonus, it displays annualized standard deviation. While I don’t particularly care for that statistic as I believe that max drawdown captures the notion of “here’s the pain on the other end of your returns” better than “here’s how much your strategy wiggles from day to day”, the fact that it’s thrown in and is a statistic that a lot of other people (particularly portfolio managers, pension fund managers, etc.) are interested in, so much the better.

Now, moving onto a more in-depth analysis of drawdown, PerformanceAnalytics has the following functionality:

dd &amp;amp;lt;- table.Drawdowns(compare[,1], top=100)
dd &amp;amp;lt;- dd[dd$Depth &amp;amp;lt; -.05,]
dd
sum(dd$"To Trough")/nrow(compare)

This brings up the following table (it seems that with multiple return streams, it’ll just default to the first one), and a derived statistic.

&amp;amp;gt; dd
         From     Trough         To   Depth Length To Trough Recovery
1  2008-12-19 2009-03-09 2010-03-16 -0.4326    310        53      257
2  2007-10-30 2008-10-15 2008-12-04 -0.3275    278       243       35
3  2013-05-22 2013-06-24 2013-09-18 -0.1617     83        23       60
4  2004-04-02 2004-05-10 2004-09-17 -0.1450    116        26       90
5  2010-04-26 2010-07-02 2010-09-21 -0.1230    104        49       55
6  2006-03-20 2006-06-19 2006-09-20 -0.1229    129        64       65
7  2007-05-08 2007-08-15 2007-10-01 -0.1229    102        70       32
8  2005-07-29 2005-10-27 2005-12-14 -0.1112     97        64       33
9  2011-06-01 2011-08-11 2011-09-06 -0.1056     68        51       17
10 2005-02-09 2005-03-22 2005-05-31 -0.1051     77        29       48
11 2010-11-05 2010-11-17 2011-02-07 -0.1022     64         9       55
12 2011-09-23 2011-10-27 2011-12-19 -0.0836     61        25       36
13 2013-09-19 2013-10-09 2013-10-17 -0.0815     21        15        6
14 2012-05-02 2012-05-18 2012-06-29 -0.0803     42        13       29
15 2012-10-18 2012-11-15 2012-11-29 -0.0721     28        19        9
16 2014-09-02 2014-10-13 2014-11-05 -0.0679     47        30       17
17 2008-12-05 2008-12-08 2008-12-16 -0.0589      8         2        6
18 2011-02-18 2011-03-16 2011-04-01 -0.0580     30        18       12
19 2014-07-23 2014-08-06 2014-08-15 -0.0536     18        11        7
20 2012-04-03 2012-04-10 2012-04-26 -0.0517     17         5       12

What I did was I simply wanted to query the table for all drawdowns that were more than 5%, or the 100 biggest drawdowns (though considering that we have about 20 drawdowns over about a decade, it seems the rule is 2 drawdowns over 5% per year, give or take, and this is a pretty volatile strategy). Lastly, I wanted to know the proportion of the time that someone watching the strategy will be feeling the pain of watching the strategy go to those depths, so I took the sum of the “To Trough” column and divided it by the amount of days of the backtest. This is the result:

&amp;amp;gt; sum(dd$"To Trough")/nrow(compare)
[1] 0.3005505

I’m fairly certain some individuals more seasoned than I am would do something different given this information and functionality, and if so, feel free to leave a comment, but this is just a licked-finger-in-the-air calculation I did. So, 30% of the time, whoever is investing real money into this will want to go and grab a few more drinks than usual.

Let’s do the same analysis for the relative performance.

tmp &amp;amp;lt;- rbind(table.AnnualizedReturns(diff), maxDrawdown(diff))
rownames(tmp)[4] &amp;amp;lt;- "Worst Drawdown"
tmp

When running only one set of returns, apparently the last row will just simply be called “4”, so I had to manually rename that row. Here’s the result:

&amp;amp;gt; tmp
                          portfolio.returns
Annualized Return                 0.1033000
Annualized Std Dev                0.2278000
Annualized Sharpe (Rf=0%)         0.4535000
Worst Drawdown                    0.3874082

Far from spectacular, but there it is for what it’s worth.

Now the drawdowns.

dd &amp;amp;lt;- table.Drawdowns(diff, top=100)
dd &amp;amp;lt;- dd[dd$Depth &amp;amp;lt; -.05,]
dd
sum(dd$"To Trough")/nrow(diff)

With the following result:

&amp;amp;gt; dd
         From     Trough         To   Depth Length To Trough Recovery
1  2008-11-21 2009-06-22 2013-04-04 -0.3874   1097       145      952
2  2008-10-28 2008-11-04 2008-11-19 -0.2328     17         6       11
3  2007-11-27 2008-06-12 2008-10-07 -0.1569    218       137       81
4  2013-05-03 2013-06-25 2013-12-26 -0.1386    165        37      128
5  2008-10-13 2008-10-13 2008-10-24 -0.1327     10         1        9
6  2005-06-08 2006-06-28 2006-11-08 -0.1139    360       267       93
7  2004-03-29 2004-05-10 2004-09-20 -0.1102    121        30       91
8  2007-03-08 2007-06-12 2007-11-26 -0.0876    183        67      116
9  2005-02-10 2005-03-28 2005-05-31 -0.0827     76        31       45
10 2014-09-02 2014-09-17 2014-10-14 -0.0613     31        12       19

In short, the spikes in outperformance gave us some pretty…interesting…drawdown statistics, which just essentially meant that the strategy wasn’t roaring at the same exact time that the SPY had its bounce from the bottom. And for interest, my finger in the air pain statistic:

&amp;amp;gt; sum(dd$"To Trough")/nrow(diff)
[1] 0.2689908

So approximately 27% of the time, the strategy underperforms its benchmark–meaning that 73% of the time, you’re fairly happy–assuming you’ve chosen the correct benchmark.

In short, overall, more freebies from Harry Long with a title to attract readers. The strategy is what it is–something that boasts strong absolute returns and definitely outperforms SPY, though like anything else, has its moments that it’ll burn you (no pain, no gain, as they say). However, the quicker statistics table functionality combined with the more in-depth drawdown analysis is something that I am definitely happy to have stumbled upon.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

An Introduction to Change Points (packages: ecp and BreakoutDetection)

A forewarning, this post is me going out on a limb, to say the least. In fact, it’s a post/project requested from me by Brian Peterson, and it follows a new paper that he’s written on how to thoroughly replicate research papers. While I’ve replicated results from papers before (with FAA and EAA, for instance), this is a first for me in terms of what I’ll be doing here.

In essence, it is a thorough investigation into the paper “Leveraging Cloud Data to Mitigate User Experience from ‘Breaking Bad’”, and follows the process from the aforementioned paper. So, here we go.

*********************

Twitter Breakout Detection Package
Leveraging Cloud Data to Mitigate User Experience From ‘Breaking Bad’

Summary of Paper

Introduction: in a paper detailing the foundation of the breakout detection package (arXiv ID 1411.7955v1), James, Kejariwal, and Matteson demonstrate an algorithm that detects breakouts in twitter’s production-level cloud data. The paper begins by laying the mathematical foundation and motivation for energy statistics, the permutation test, and the E-divisive with medians algorithm, which create a fast way of detecting a shift in median between two nonparametric distributions that is robust to the presence of anomalies. Next, the paper demonstrates a trial run through some of twitter’s production cloud data, and compares the non-parametric E-divisive with medians to an algorithm called PELT. For the third topic, the paper discusses potential applications, one of which is quantitative trading/computational finance. Lastly, the paper states its conclusion, which is the addition of the E-divisive with medians algorithm to the existing literature of change point detection methodologies.

The quantitative and computational methodologies for the paper use a modified variant of energy statistics more resilient against anomalies through the use of robust statistics (viz. median). The idea of energy statistics is to compare the distances of means of two random variables contained within a larger time series. The hypothesis test to determine if this difference is statistically significant is called the permutation test, which permutes data from the two time series a finite number of times to make the process of comparing permuted time series computationally tractable. However, the presence of anomalies, such as in twitter’s production cloud data, would limit the effectiveness of using this process when using simple means. To that end, the paper proposes using the median, and due to the additional computational time resulting from the weaker distribution assumptions to extend the generality of the procedure, the paper devises the E-divisive with medians algorithms, one of which works off of distances between observations, and one works with the medians of the observations themselves (as far as I understand). To summarize, the E-divisive with medians algorithms exist as a way of creating a computationally tractable procedure for determining whether or not a new chunk of time series data is considerably different from the previous through the use of advanced distance statistics robust to anomalies such as those present in twitter’s cloud data.

To compare the performance of the E-divisive with medians algorithms, the paper compares the algorithms to an existing algorithm called PELT (which stands for Pruned Extract Linear Time) in various quantitative metrics, such as “Time To Detect”, meaning the exact moment of the breakout to when the algorithms report it (if at all), along with precision, recall, and the F-measure, defined as the product of precision and recall over their respective sum. Comparing PELT to the E-divisive with medians algorithm showed that the E-divisive algorithm outperformed the PELT algorithm in the majority of data sets. Even when anomalies were either smoothed by taking the rolling median of their neighbors, or by removing them altogether, the E-divisive algorithm still outperformed PELT. Of the variants of the EDM algorithm (EDM head, EDM tail, and EDM-exact), the EDM-tail variant (i.e. the one using the most recent observations) was also quickest to execute. However, due to fewer assumptions about the nature of the underlying generating distributions, the various E-divisive algorithms take longer to execute than the PELT algorithm, with its stronger assumptions, but worse general performance. To summarize, the EDM algorithms outperform PELT in the presence of anomalies, and generally speaking, the EDM-tail variant seems to work best when considering computational running time as well.

The next section dealt with the history and applications of change-point/breakout detection algorithms, in fields such as finance, medical applications, and signal processing. As finance is of a particular interest, the paper acknowledges the ARCH and various flavors of GARCH models, along with the work of James and Matteson in devising a trading strategy based on change-point detection. Applications in genomics to detect cancer exist as well. In any case, the paper cites many sources showing the extension and applications of change-point/breakout detection algorithms, of which finance is one area, especially through work done by Matteson. This will be covered further in the literature review.

To conclude, the paper proposes a new algorithm called the E-divisive with medians, complete with a new statistical permutation test using advanced distance statistics to determine whether or not a time series has had a change in its median. This method makes fewer assumptions about the nature of the underlying distribution than a competitive algorithm, and is robust in the face of anomalies, such as those found in twitter’s production cloud data. This algorithm outperforms a competing algorithm which possessed stronger assumptions about the underlying distribution, detecting a breakout sooner in a time series, even if it took longer to run. The applications of such work range from finance to medical devices, and further beyond. As change-point detection is a technique around which trading strategies can be constructed, it has particular relevance to trading applications.

Statement of Hypothesis

Breakouts can occur in data which does not conform to any known regular distribution, thus rendering techniques that assume a certain distribution less effective. Using the E-divisive with medians algorithm, the paper attempts to predict the presence of breakouts using time series with innovations from no regular distribution as inputs, and if effective, will outperform an existing algorithm that possesses stronger assumptions about distributions. To validate or refute a more general form of this hypothesis, which is the ability of the algorithm to detect breakouts in a timely fashion, this summary test it on the cumulative squared returns of the S&P 500, and compare the analysis created by the breakpoints to the analysis performed by Dr. Robert J. Frey of Keplerian Finance, a former managing director at Renaissance Technologies.

Literature Review

Motivation

A good portion of the practical/applied motivation of this paper stems from the explosion of growth in mobile internet applications, A/B testing, and other web-specific reasons to detect breakouts. For instance, longer loading time on a mobile web page necessarily results in lower revenues. To give another example, machines in the cloud regularly fail.

However, the more salient literature regarding the topic is the literature dealing with the foundations of the mathematical ideas behind the paper.

Key References

Paper 1:

David S. Matteson and Nicholas A. James. A nonparametric approach for multiple change point analysis of multivariate data. Journal of the American Statistical Association, 109(505):334–345, 2013.

Thesis of work: this paper is the original paper for the e-divisive and e-agglomerative algorithms, which are offline, nonparametric methods of detecting change points in time series. Unlike Paper 3, this paper lays out the mathematical assumptions, lemmas, and proofs for a formal and mathematical presentation of the algorithms. Also, it documents performance against the PELT algorithm, presented in Paper 6 and technically documented in Paper 5. This performance compares favorably. The source paper being replicated builds on the exact mathematics presented in this paper, and the subject of this report uses the ecp R package that is the actual implementation/replication of this work to form a comparison for its own innovations.

Paper 2:

M. L. Rizzo and G. J. Sz´ekely. DISCO analysis: A nonparametric extension of analysis of variance. The Annals of Applied Statistics, 4(2):1034–1055, 2010

Thesis of work: this paper generalizes the ANOVA using distance statistics. This technique aims to find differences among distributions outside their sample means. Through the use of distance statistics, the techniques aim to more generally answer queries about the nature of distributions (EG identical means, but different distributions as a result of different factors). Its applicability to the source paper is that it forms the basis of the ideas for the paper’s divergence measure, as detailed in its second section.

Paper 3:

Nicholas A. James and David S. Matteson. ecp: An R package for nonparametric multiple change point analysis of multivariate data. Technical report, Cornell University, 2013.

Thesis of work: the paper introduces the ecp package which contains the e-agglomerative and e-divisive algorithms for detecting change points in time series in the R statistical programming language (in use on at least one elite trading desk). The e-divisive method recursively partitions a time series and uses a permutation test to determine change points, but it is computationally intensive. The e-agglomerative algorithm allows for inputs from the user for initial time-series segmentation and is a computationally faster algorithm. Unlike most academic papers, this paper also includes examples of data and code in order to facilitate the use of these algorithms. Furthermore, the paper includes applications to real data, such as the companies found in the Dow Jones Industrial Index, further proving the effectiveness of these methods. This paper is important to the topic in question as the E-divisive algorithm created by James and Matteson form the base changepoint detection process for which the paper builds its own innovations for, and visually compares against; furthermore, the source paper restates many of the techniques found in this paper.

Paper 4:

Owen Vallis, Jordan Hochenbaum, and Arun Kejariwal. A novel technique for long-term anomaly detection in the cloud. In 6th USENIX Workshop on Hot Topics in Cloud Computing (HotCloud 14), June 2014.

Thesis of work: the paper proposes the use of piecewise median and median absolute deviation statistics to detect anomalies in time series. The technique builds upon the ESD (Extreme Studentized Deviate) technique and uses piecewise medians to approximate a long-term trend, before extracting seasonality effects from periods shorter than two weeks. The piecewise median method of anomaly detection has a greater F-measure of detecting anomalies than does the standard STL (seasonality trend loess decomposition) or quantile regression techniques. Furthermore, piecewise median executes more than three times faster. The relevance of this paper to the source paper is that it forms the idea of using robust statistics and building the techniques in the paper upon the median as opposed to the mean.

Paper 5:

Rebecca Killick and Kaylea Haynes. changepoint: An R package for changepoint analysis

Thesis of work: manual for the implementation of the PELT algorithm written by Rebecca Killick and Kaylea Haynes. This package is a competing change-point detection package, mainly focused around the Pruned Extraction Linear Time algorithm, although containing other worse algorithms, such as the segment neighborhoods algorithm. Essentially, it is a computational implementation of the work in Paper 2. Its application toward the source paper is that the paper at hand compares its own methodology against PELT, and often outperforms it.

Paper 6:

Rebecca Killick, Paul Fearnhead, and IA Eckley. Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association, 107(500):1590–1598, 2012

Thesis of work: the paper proposes an algorithm (PELT) that scales linearly in running time with the size of the input time series to detect exact locations of change points. The paper aims to replace both an approximate binary partitioning algorithm, and an optimal segmentation algorithm that doesn’t involve a pruning mechanism to speed up the running time. The paper uses an MLE algorithm at the heart of its dynamic partitioning in order to locate change points. The relevance to the source paper is that through the use of the non-robust MLE procedure, this algorithm is vulnerable to poor performance due to the presence of anomalies/outliers in the data, and thus underperforms the new twitter change point detection methodology which employs robust statistics.

Paper 7:

Wassily Hoeffding. The strong law of large numbers for u-statistics. Institute of Statistics mimeo series, 302, 1961.

Thesis of work: this paper establishes a convergence of the mean of tuples of many random variables to the mean of said random variables, given enough such observations. This paper is a theoretical primer on establishing the above thesis. The mathematics involve use of measure theory and other highly advanced and theoretical manipulations. Its relevance to the source paper is in its use to establish a convergence of an estimated characteristic function.

Similar Work

In terms of financial applications, the papers covering direct applications of change points to financial time series are listed above. Particularly, David Matteson presented his ecp algorithms at R/Finance several years ago, and his work is already in use on at least one professional trading desk. Beyond this, the paper cites works on technical analysis and the classic ARCH and GARCH papers as similar work. However, as this change point algorithm is created to be a batch process, direct comparison with other trend-following (that is, breakout) methods would seem to be a case of apples and oranges, as indicators such as MACD, Donchian channels, and so on, are online methods (meaning they do not have access to the full data set like the e-divisive and the e-divisive with medians algorithms do). However, they are parameterized in terms of their lookback period, and are thus prone to error in terms of inaccurate parameterization resulting from a static lookback value.

In his book Cycle Analytics for Traders, Dr. John Ehlers details an algorithm for computing the dominant cycle of a security—that is, a way to dynamically parameterize the lookback parameter, and if this were to be successfully implemented in R, it may very well allow for improved breakout detection methods than the classic parameterized indicators popularized in the last century.

References With Implementation Hints

Reference 1: Breakout Detection In The Wild

This blog post is a reference contains the actual example included in the R package for the model, written by one of the authors of the source paper. As the data used in the source paper is proprietary twitter production data, and the model is already implemented in the package discussed in this blog post, this makes the package and the included data the go-to source for starting to work with the results presented in the source paper.

Reference 2: Twitter BreakoutDetection R package evaluation

This blog post is that of a blogger altering the default parameters in the model. His analysis of traffic to his blog contains valuable information as to greater flexibility in the use of the R package that is the implementation of the source paper.

Data

The data contained in the source paper comes from proprietary twitter cloud production data. Thus, it is not realistic to obtain a copy of that particular data set. However, one of the source paper’s co-authors, Arun Kejariwal, was so kind as to provide a tutorial, complete with code and sample data, for users to replicate at their convenience. It is this data that we will use for replication.

Building The Model

Stemming from the above, we are fortunate that the results of the source paper have already been implemented in twitter’s released R package, BreakoutDetection. This package has been written by Nicholas A. James, a PhD candidate at Cornell University studying under Dr. David S. Matteson. His page is located here.

In short, all that needs to be done on this end is to apply the model to the aforementioned data.

Validate the Results

To validate the results—that is, to obtain the same results as one of the source paper’s authors, we will execute the code on the data that he posted on his blog post (see Reference 1).

require(devtools)
install_github(repo="BreakoutDetection", username="twitter")
require(BreakoutDetection)

data(Scribe)
res = breakout(Scribe, min.size=24, method='multi', beta=.001, degree=1, plot=TRUE)
res$plot

This is the resulting image, identical from the blog post.

Validation of the Hypothesis

This validation was inspired by the following post:

The Relevance of History

The post was written by Dr. Robert J. Frey, professor of Applied Math and Statistics at Stony Brook University, the head of its Quantitative Finance program, and former managing director at Renaissance Technologies (yes, the Renaissance Technologies founded by Dr. Jim Simons). While the blog is inactive at the moment, I sincerely hope it will become more active again.

Essentially, it uses mathematica to detect changes in the slope of cumulative squared returns, and the final result is a map of spikes, mountains, and plains, the x-axis being time, and the y-axis the annualized standard deviation. Using the more formalized e-divisive and e-divisive with medians algorithms, this analysis will attempt to detect change points, and use the PerformanceAnalytics library to compute the annualized standard deviation from the data of the GSPC returns itself, and output a similarly-formatted plot.

Here’s the code:

require(quantmod)
require(PerformanceAnalytics)

getSymbols("^GSPC", from = "1984-12-25", to = "2013-05-31")
monthlyEp <- endpoints(GSPC, on = "months")
GSPCmoCl <- Cl(GSPC)[monthlyEp,]
GSPCmoRets <- Return.calculate(GSPCmoCl)
GSPCsqRets <- GSPCmoRets*GSPCmoRets
GSPCsqRets <- GSPCsqRets[-1,] #remove first NA as a result of return computation
GSPCcumSqRets <- cumsum(GSPCsqRets)
plot(GSPCcumSqRets)

This results in the following image:

So far, so good. Let’s now try to find the amount of changepoints that Dr. Frey’s graph alludes to.

t1 <- Sys.time()
ECPmonthRes <- e.divisive(X = GSPCsqRets, min.size = 2)
t2 <- Sys.time()
print(t2 - t1)

t1 <- Sys.time()
BDmonthRes <- breakout(Z = GSPCsqRets, min.size = 2, beta=0, degree=1)
t2 <- Sys.time()
print(t2 - t1)

ECPmonthRes$estimates
BDres$loc

With the following results:

> ECPmonthRes$estimates
[1]   1 285 293 342
> BDres$loc
[1] 47 87

In short, two changepoints for each. Far from the 20 or so regimes present in Dr. Frey’s analysis. So, not close to anything that was expected. My intuition tells me that the main reason for this is that these algorithms are data-hungry, and there is too little data for them to do much more than what they have done thus far. So let’s go the other way and use daily data.

dailySqRets <- Return.calculate(Cl(GSPC))*Return.calculate(Cl(GSPC))
dailySqRets <- dailySqRets["1985::"]

plot(cumsum(dailySqRets))

And here’s the new plot:

First, let’s try the e-divisive algorithm from the ecp package to find our changepoints, with a minimum size of 20 days between regimes. (Blog note: this is a process that takes an exceptionally long time. For me, it took more than 2 hours.)

t1 <- Sys.time()
ECPres <- e.divisive(X = dailySqRets, min.size=20)
t2 <- Sys.time()
print(t2 - t1)
Time difference of 2.214813 hours

With the following results:

index(dailySqRets)[ECPres$estimates]
 [1] "1985-01-02" "1987-10-14" "1987-11-11" "1998-07-21" "2002-07-01" "2003-07-28" "2008-09-15" "2008-12-09"
 [9] "2009-06-02" NA   

The first and last are merely the endpoints of the data. So essentially, it encapsulates Black Monday and the crisis, among other things. Let’s look at how the algorithm split the volatility regimes. For this, we will use the xtsExtra package for its plotting functionality (thanks to Ross Bennett for the work he did in implementing it).

require(xtsExtra)
plot(cumsum(dailySqRets))
xtsExtra::addLines(index(dailySqRets)[ECPres$estimates[-c(1, length(ECPres$estimates))]], on = 1, col = "blue", lwd = 2)

With the resulting plot:

In this case, the e-divisive algorithm from the ecp package does a pretty great job segmenting the various volatility regimes, as can be thought of roughly as the slope of the cumulative squared returns. The algorithm’s ability to accurately cluster the Black Monday events, along with the financial crisis, shows its industrial-strength applicability. How does this look on the price graph?

plot(Cl(GSPC))
xtsExtra::addLines(index(dailySqRets)[ECPres$estimates[-c(1, length(ECPres$estimates))]], on = 1, col = "blue", lwd = 2)

In this case, Black Monday is clearly visible, along with the end of the Clinton bull run through the dot-com bust, the consolidation, the run-up to the crisis, the crisis itself, the consolidation, and the new bull market.

Note that the presence of a new volatility regime may not necessarily signify a market top or bottom, but the volatility regime detection seems to have worked very well in this case.

For comparison, let’s examine the e-divisive with medians algorithm.

t1 <- Sys.time()
BDres <- breakout(Z = dailySqRets, min.size = 20, beta=0, degree=1)
t2 <- Sys.time()
print(t2-t1)

BDres$loc
index(dailySqRets)[BDres$loc]

With the following result:

Time difference of 2.900167 secs
> BDres$loc
[1] 5978
> BDres$loc
[1] 5978
> index(dailySqRets)[BDres$loc]
[1] "2008-09-12"

So while the algorithm is a lot faster, its volatility regime detection, it only sees the crisis as the one major change point. Beyond that, to my understanding, the e-divisive with medians algorithm may be “too robust” (even without any penalization) against anomalies (after all, the median is robust to changes in 50% of the data). In short, I think that while it clearly has applications, such as twitter cloud production data, it doesn’t seem to obtain a result that’s in the ballpark of two other separate procedures.

Lastly, let’s try and create a plot similar to Dr. Frey’s, with spikes, mountains, and plains.

require(PerformanceAnalytics)
GSPCrets <- Return.calculate(Cl(GSPC))
GSPCrets <- GSPCrets["1985::"]
GSPCrets$regime <- ECPres$cluster
GSPCrets$annVol <- NA

for(i in unique(ECPres$cluster)) {
  regime <- GSPCrets[GSPCrets$regime==i,]
  annVol <- StdDev.annualized(regime[,1])
  GSPCrets$annVol[GSPCrets$regime==i,] <- annVol
}

plot(GSPCrets$annVol, ylim=c(0, max(GSPCrets$annVol)), main="GSPC volatility regimes, 1985 to 2013-05")

With the corresponding image, inspired by Dr. Robert Frey:

This concludes the research replication.

********************************

Whew. Done. While I gained some understanding of what change points are useful for, I won’t profess to be an expert on them (some of the math involved uses PhD-level mathematics such as characteristic functions that I never learned). However, it was definitely interesting pulling together several different ideas and uniting them under a rigorous process.

Special thanks for this blog post:

Brian Peterson, for the process paper and putting a formal structure to the research replication process (and requesting this post).
Robert J. Frey, for the “volatility landscape” idea that I could objectively point to as an objective benchmark to validate the hypothesis of the paper.
David S. Matteson, for the ecp package.
Nicholas A. James, for the work done in the BreakoutDetection package (and clarifying some of its functionality for me).
Arun Kejariwal, for the tutorial on using the BreakoutDetection package.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

An Update On EAA and a Volatility Strategy

Again, before starting this post, I’d like to inform readers that the book Quantitative Trading With R, written by Harry Georgakopoulos, with contributions from myself, is now available for order on Amazon. Already, it has garnered a pair of five-star reviews, and it deals not only with quantstrat, but with aspects such as spread trading, high frequency data, and options. I highly recommend it.

So, first things first, I want to inform everyone that EAA (that is, Elastic Asset Allocation, the new algorithm recently released by Dr. Wouter Keller a couple of weeks ago) is now in my IKTrading package. I made some modifications to deal with incongruous security starting dates (that is, handled NA momentum, and so on, similarly to the process in FAA). Again, no particular guarantees, but at this point, I think the algorithm won’t regularly break (but I may be missing some edge case, so feedback is always appreciated). Also, after thinking about it a bit more, I don’t foresee EAA as it stands being able to make use of a conditional correlation algorithm, since rather than using correlation simply for security selection, it uses correlations to make weighting decisions, which raises the question of what the correlation value of the first security would be. 0? -1? Ideas on how to address this are always welcome, since applying conditional correlation outside of a ranking context is a topic now of interest to me.

Furthermore, TrendXplorer has recently posted his own post on EAA after seeing mine on his blog. It is *very* comprehensive, and for those that are more inclined towards AmiBroker, you’ll be in Nirvana. It can be found here. Also, it seems he has done some work with another SeekingAlpha contributor named Cliff Smith (and seems to have worked hand in hand with him), and thus, had a far more positive experience than I did going solo replicating Harry Long’s strategies (or, if some of you may like, marketing materials). TrendXplorer did some work with a strategy called QTS, which I hope I’ll be able to cover in the near future. That can all be found here. So, I’d like to formally extend thanks to TrendXplorer for the work he has done with both EAA, and also pointing me towards yet another viable asset allocation strategy.

In terms of my own updated EAA, to test it out, I added Tesla Motors to the original seven securities. So let’s look at the as-of-now-current EAA.

"EAA" <- function(monthlyPrices, wR=1, wV=0, wC=.5, wS=2, errorJitter=1e-6, 
                cashAsset=NULL, bestN=1+ceiling(sqrt(ncol(monthlyPrices))),
                enableCrashProtection = TRUE, returnWeights=FALSE, monthlyRiskFree=NULL) {
  returns <- Return.calculate(monthlyPrices)
  returns <- returns[-1,] #return calculation uses one observation
  if(!is.null(monthlyRiskFree)) {
    returnsRF <- Return.calculate(monthlyRiskFree)
    returnsRF <- returnsRF[-1,]
  }
  
  if(is.null(cashAsset)) {
    returns$zeroes <- 0
    cashAsset <- "zeroes"
    warning("No cash security specified. Recommended to use one of: quandClean('CHRIS/CME_US'), SHY, or VFISX. 
            Using vector of zeroes instead.")
  }
  
  cashCol <- grep(cashAsset, colnames(returns))
  
  weights <- list()
  for(i in 1:(nrow(returns)-11)) {
    returnsData <- returns[i:(i+11),] #each chunk will be 12 months of returns data
    #per-month mean of cumulative returns of 1, 3, 6, and 12 month periods
    periodReturn <- ((returnsData[12,] + Return.cumulative(returnsData[10:12,]) + 
                      Return.cumulative(returnsData[7:12,]) + Return.cumulative(returnsData)))/22
    
    if(!is.null(monthlyRiskFree)) {
      rfData <- returnsRF[i:(i+11),]
      rfReturn <- ((rfData[12,] + Return.cumulative(rfData[10:12,]) + 
                    Return.cumulative(rfData[7:12,]) + Return.cumulative(rfData)))/22
      periodReturn <- periodReturn - as.numeric(rfReturn)
    }
    
    vols <- StdDev.annualized(returnsData) 
    mktIndex <- xts(rowMeans(returnsData, na.rm=TRUE), order.by=index(returnsData)) #equal weight returns of universe
    cors <- cor(returnsData, mktIndex) #correlations to market index
    
    weightedRets <- periodReturn ^ wR
    weightedCors <- (1 - as.numeric(cors)) ^ wC
    weightedVols <- (vols + errorJitter) ^ wV
    wS <- wS + errorJitter
    
    z <- (weightedRets * weightedCors / weightedVols) ^ wS #compute z_i and zero out negative returns
    z[periodReturn < 0] <- 0
    crashProtection <- sum(z==0, na.rm=TRUE)/sum(!is.na(z)) #compute crash protection cash cushion
    
    orderedZ <- sort(as.numeric(z), decreasing=TRUE)
    selectedSecurities <- z >= orderedZ[bestN]
    preNormalizedWeights <- z*selectedSecurities #select top N securities, keeping z_i scores
    periodWeights <- preNormalizedWeights/sum(preNormalizedWeights, na.rm=TRUE) #normalize
    if (enableCrashProtection) {
      periodWeights <- periodWeights * (1-crashProtection) #CP rule
    }
    periodWeights[is.na(periodWeights)] <- 0
    weights[[i]] <- periodWeights
  }
  
  weights <- do.call(rbind, weights)
  weights[, cashCol] <- weights[, cashCol] + 1-rowSums(weights) #add to risk-free asset all non-invested weight
  strategyReturns <- Return.rebalancing(R = returns, weights = weights) #compute strategy returns
  if(returnWeights) {
    return(list(weights, strategyReturns))
  } else {
    return(strategyReturns)
  }
}

Essentially, little changed aside from some lines dealing with NAs (AKA securities that were not yet around at the time whose prices are given as NAs).

To test out whether the algorithm worked, I added TSLA to see if it didn’t break the code. Here is the new test code.

require(quantmod)
require(PerformanceAnalytics)

symbols <- c("VTSMX", "FDIVX", "VEIEX", "VBMFX", "VFISX", "VGSIX", "QRAAX", "TSLA")

getSymbols(symbols, from="1990-01-01")
prices <- list()
for(i in 1:length(symbols)) {
  prices[[i]] <- Ad(get(symbols[i]))  
}
prices <- do.call(cbind, prices)
colnames(prices) <- gsub("\\.[A-z]*", "", colnames(prices))
ep <- endpoints(prices, "months")
prices <- prices[ep,]
prices <- prices["1997-03::"]

getSymbols("^IRX", from="1990-01-01")
dailyYield <- (1+(Cl(IRX)/100))^(1/252) - 1
threeMoPrice <- cumprod(1+dailyYield)
threeMoPrice <- threeMoPrice["1997-03::"]
threeMoPrice <- threeMoPrice[endpoints(threeMoPrice, "months"),]

offensive <- EAA(prices, cashAsset="VBMFX", bestN=3)
defensive <- EAA(prices, cashAsset="VBMFX", bestN=3, wS=.5, wC=1)
offRF <- EAA(prices, cashAsset="VBMFX", bestN=3, monthlyRiskFree = threeMoPrice)
defRF <- EAA(prices, cashAsset="VBMFX", bestN=3, wS=.5, wC=1, monthlyRiskFree = threeMoPrice)
compare <- cbind(offensive, defensive, offRF, defRF)
colnames(compare) <- c("Offensive", "Defensive", "OffRF", "DefRF")
stats <- rbind(Return.annualized(compare)*100, StdDev.annualized(compare)*100, maxDrawdown(compare)*100, SharpeRatio.annualized(compare))
rownames(stats)[3] <- "Worst Drawdown"
charts.PerformanceSummary(compare)
stats

With the following statistics table and equity curve:

> stats
                                 Offensive Defensive      OffRF     DefRF
Annualized Return               17.6174693 13.805683 16.7376777 13.709368
Annualized Standard Deviation   22.7328695 13.765444 22.3854966 13.504313
Worst Drawdown                  25.3534015 12.135310 25.3559118 12.146654
Annualized Sharpe Ratio (Rf=0%)  0.7749778  1.002923  0.7477019  1.015184

Essentially, TSLA — a high momentum, high-volatility stock causes some consternation in the offensive variant of the algorithm. Let’s look at the weight statistics of TSLA when it was in the portfolio.

test <- EAA(prices, cashAsset = "VBMFX", bestN=3, returnWeights=TRUE)
weights <- test[[1]]
summary(weights$TSLA[weights$TSLA > 0])

With the results:

    Index                 TSLA        
 Min.   :2011-07-29   Min.   :0.01614  
 1st Qu.:2012-09-14   1st Qu.:0.32345  
 Median :2013-07-31   Median :0.48542  
 Mean   :2013-06-20   Mean   :0.51415  
 3rd Qu.:2014-04-15   3rd Qu.:0.75631  
 Max.   :2014-12-31   Max.   :0.95793  

Also, to be clear, R’s summary function was not created with xts type objects in mind, so the Index statistics are just pure nonsense (R is trying to do summary statistics on the underlying numerical values of the date index — they have no relation to the TSLA weights), so if you ever call summary on anything in an xts, be aware that it isn’t actually providing you the dates of the corresponding weights (if they exist at all — E.G. the mean of the weights isn’t an actual weight at any point in time).

In any case, it seems that the offensive variant of the algorithm is susceptible to creating portfolios that are very poorly diversified, since the offensive variant doesn’t place any weight on security volatility–simply correlation. So if there was a very volatile instrument that was on a roaring trend, EAA would tell you to just place your entire portfolio in that one instrument–which of course, can be the correct thing to do if you know for certain that said trend will continue, until, of course, it doesn’t.

I’m sure there are still some methods to account for instruments of wildly different risk/return profiles, even without the need of additional code, by varying the parameters. I just wanted to demonstrate the need to be aware of this phenomenon, which I happened upon simply by testing the portfolio for incongruous starting dates and just so happened to pick a “hot topic” stock.

Last (for this post), I’d like to make readers aware that the blogger Volatility Made Simple has created a variant of a strategy I had written about earlier (again, thanks to Mr. Helmuth Vollmeier for providing the initial foundation), in which he mixed signals from the three variants I had found to be in stable regions, and I’m really happy he has done so, as he’s one of the first people who have explicitly extended my work.

Unfortunately, said strategy is currently in drawdown. However, looking at its drawdown curve against that of XIV itself, it seems that volatility has been doing crazy things lately, and the drawdown has been worse in the past. I am concerned, however, that it may be a strategy prone to overfitting, and it’s a constant reminder that there is still more to learn, and more techniques to use to convince oneself that a backtest isn’t just an overfit, data-mined, sample-dependent illusion with good marketing that will break down immediately upon looking at a larger sample. However, as I did not originate the strategy myself, I’d at least like to hope that whoever was the first person who came up with the VXV/VXMT ratio idea had some good rationale for the strategy to begin with.

In the immediate future, I’ll be looking into change point analysis and twitter’s new breakout detection package.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

Adding a Risk-Free Rate To Your Analyses

First off, before beginning this post, I’d like to make my readers aware of the release of a book that I contributed almost an entire chapter for.

Quantitative Trading With R is a primer on quantitative trading in R written by Harry Georgakopoulos, one of Chicago’s better quants. I contributed almost the entire chapter on quantstrat. If you’ve been able to follow and understand the code I write on this blog, then said chapter will mostly be review and a basic nuts and bolts reference. But for those of my readers who gloss over the code and wait for the punchline, I highly recommend it. In addition, there are chapters on high frequency trading, options, spreads, and other things that I do not believe are available in any other book that actually teaches readers the details of implementation. Now, onto the post.

As part of my continuation of Elastic Asset Allocation, I wanted to cover how to implement a measure of a risk-free rate in your analyses. In this post, I’ll analyze two slight variations of EAA from the last post.

Essentially, the idea that rather than look at absolute return, we should look at *excess* return over some sort of risk-free rate, such as the 13-week treasury bill.

Luckily, Yahoo actually *has* a way of getting the returns of the risk-free asset, namely, IRX. But first, let’s get the similarities to the last post out of the way.

require(quantmod)
require(PerformanceAnalytics)

symbols <- c("VTSMX", "FDIVX", "VEIEX", "VBMFX", "VFISX", "VGSIX", "QRAAX")

getSymbols(symbols, from="1990-01-01")
prices <- list()
for(i in 1:length(symbols)) {
  prices[[i]] <- Ad(get(symbols[i]))  
}
prices <- do.call(cbind, prices)
colnames(prices) <- gsub("\\.[A-z]*", "", colnames(prices))
ep <- endpoints(prices, "months")
prices <- prices[ep,]
prices <- prices["1997-03::"]

Okay, everything fine so far, same as before. Now here’s the new innovation, brought to my attention by TrendXplorer. It turns out that the IRX index is actually the annualized yield for the short-term (three month) treasuries. So by adding 1, raising it to the 252nd root, and taking the cumulative product, we can actually get the “price” of the risk-free rate, and from that, compute daily returns (this is most likely redundant, but I want all my returns computed the same way).

getSymbols("^IRX", from="1990-01-01")
dailyYield <- (1+(Cl(IRX)/100))^(1/252) - 1
threeMoPrice <- cumprod(1+dailyYield)
threeMoPrice <- threeMoPrice["1997-03::"]
threeMoPrice <- threeMoPrice[endpoints(threeMoPrice, "months"),]

So how does this fit into EAA? Well, simply, I added a new argument called monthlyRiskFree, which will let a user pass in the monthly price series of the risk-free asset, in this case the derived IRX price series. That information is then used to compute a risk-free return, which is subtracted from the returns of all assets, and rather than taking the absolute return of the assets in the universe, instead the algorithm computes the return in excess of the risk-free asset.

Here’s the modified function:

EAA <- function(monthlyPrices, wR=1, wV=0, wC=.5, wS=2, errorJitter=1e-6, 
                cashAsset=NULL, bestN=1+ceiling(sqrt(ncol(monthlyPrices))),
                enableCrashProtection = TRUE, returnWeights=FALSE, monthlyRiskFree=NULL) {
  returns <- Return.calculate(monthlyPrices)
  returns <- returns[-1,] #return calculation uses one observation
  if(!is.null(monthlyRiskFree)) {
    returnsRF <- Return.calculate(monthlyRiskFree)
    returnsRF <- returnsRF[-1,]
  }
  
  if(is.null(cashAsset)) {
    returns$zeroes <- 0
    cashAsset <- "zeroes"
    warning("No cash security specified. Recommended to use one of: quandClean('CHRIS/CME_US'), SHY, or VFISX. 
            Using vector of zeroes instead.")
  }
  
  cashCol <- grep(cashAsset, colnames(returns))
  
  weights <- list()
  for(i in 1:(nrow(returns)-11)) {
    returnsData <- returns[i:(i+11),] #each chunk will be 12 months of returns data
    #per-month mean of cumulative returns of 1, 3, 6, and 12 month periods
    periodReturn <- ((returnsData[12,] + Return.cumulative(returnsData[10:12,]) + 
                      Return.cumulative(returnsData[7:12,]) + Return.cumulative(returnsData)))/22
    
    if(!is.null(monthlyRiskFree)) {
      rfData <- returnsRF[i:(i+11),]
      rfReturn <- ((rfData[12,] + Return.cumulative(rfData[10:12,]) + 
                    Return.cumulative(rfData[7:12,]) + Return.cumulative(rfData)))/22
      periodReturn <- periodReturn - as.numeric(rfReturn)
    }
    
    vols <- StdDev.annualized(returnsData) 
    mktIndex <- xts(rowMeans(returnsData), order.by=index(returnsData)) #equal weight returns of universe
    cors <- cor(returnsData, mktIndex) #correlations to market index
    
    weightedRets <- periodReturn ^ wR
    weightedCors <- (1 - as.numeric(cors)) ^ wC
    weightedVols <- (vols + errorJitter) ^ wV
    wS <- wS + errorJitter
    
    z <- (weightedRets * weightedCors / weightedVols) ^ wS #compute z_i and zero out negative returns
    z[periodReturn < 0] <- 0
    crashProtection <- sum(z==0)/length(z) #compute crash protection cash cushion
    
    orderedZ <- sort(as.numeric(z), decreasing=TRUE)
    selectedSecurities <- z >= orderedZ[bestN]
    preNormalizedWeights <- z*selectedSecurities #select top N securities, keeping z_i scores
    periodWeights <- preNormalizedWeights/sum(preNormalizedWeights) #normalize
    if (enableCrashProtection) {
      periodWeights <- periodWeights * (1-crashProtection) #CP rule
    }
    weights[[i]] <- periodWeights
  }
  
  weights <- do.call(rbind, weights)
  weights[, cashCol] <- weights[, cashCol] + 1-rowSums(weights) #add to risk-free asset all non-invested weight
  strategyReturns <- Return.rebalancing(R = returns, weights = weights) #compute strategy returns
  if(returnWeights) {
    return(list(weights, strategyReturns))
  } else {
    return(strategyReturns)
  }
}

Essentially, the key new block of code is this:

    if(!is.null(monthlyRiskFree)) {
      rfData <- returnsRF[i:(i+11),]
      rfReturn <- ((rfData[12,] + Return.cumulative(rfData[10:12,]) + 
                    Return.cumulative(rfData[7:12,]) + Return.cumulative(rfData)))/22
      periodReturn <- periodReturn - as.numeric(rfReturn)
    }

Which does exactly as I mentioned above — computes the EAA variant of the returns for the period for the risk-free asset and subtracts it from the other returns.

So how does using this new innovation compare to simply looking at absolute returns?

Let’s find out.

offensive <- EAA(prices, cashAsset="VBMFX", bestN=3)
defensive <- EAA(prices, cashAsset="VBMFX", bestN=3, wS=.5, wC=1)
offRF <- EAA(prices, cashAsset="VBMFX", bestN=3, monthlyRiskFree = threeMoPrice)
defRF <- EAA(prices, cashAsset="VBMFX", bestN=3, wS=.5, wC=1, monthlyRiskFree = threeMoPrice)
compare <- cbind(offensive, defensive, offRF, defRF)
colnames(compare) <- c("Offensive", "Defensive", "OffRF", "DefRF")
stats <- rbind(Return.annualized(compare)*100, StdDev.annualized(compare)*100, maxDrawdown(compare)*100, SharpeRatio.annualized(compare))
rownames(stats)[3] <- "Worst Drawdown"
charts.PerformanceSummary(compare)

Here’s the table of statistics:

                                Offensive Defensive     OffRF     DefRF
Annualized Return               12.206133 10.262766 11.415583 10.269146
Annualized Standard Deviation   11.352728  8.615134 10.551722  8.129250
Worst Drawdown                  12.629251  8.134785 14.351895  9.376533
Annualized Sharpe Ratio (Rf=0%)  1.075172  1.191249  1.081869  1.263234

And the corresponding chart:

In short, nope. No dice there. On the defensive portfolio, the change is negligible, and on the offensive side, it seems to encourage more risk than necessary, with…nothing to show for it, really.

That stated, just because this method didn’t pan out doesn’t mean that the actual mechanics of obtaining risk-free data are without value.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

For A New Year, A New Asset Allocation System Just Published in SSRN

Happy New Year!

So, this is something I’ve been working on before its official publication (so this is the first place on the entire internet that you’ll see it outside SSRN, and certainly one of the few places that will extend it), directly in contact with the original paper author, Dr. Wouter Keller, who published the Flexible Asset Allocation algorithm that I improved on earlier. It’s called Elastic Asset Allocation, and seems to be a simpler yet more general algorithm. Here’s the link to the SSRN

Essentially, the algorithm can be explained simply:

Use monthly data.

For a 12 month period, we want the monthly average of the 1, 3, 6, and 12 month cumulative returns (that is, sum all those up, divide by 22), the volatility of the individual monthly returns, and the correlation of the universe of returns to the monthly average of the returns. Then arrange those values in the following expression:

z_i = (r_i ^ wR * (1-c_i) ^ wC / (v_i + error) ^ wV) ^ (wS + error) if r_i > 0, 0 otherwise, where:

r_i are the average returns
c_i are the correlations
v_i are the volatilities

wR, wC, and wV are the respective weights for returns, correlations, and volatilities.

Next, select the top N of P assets to include in the portfolio.

Then, the weights for each security can be expressed as the normalized values of the sum of selected z_i’s.

If a crash protection rule is enabled, compute CP, which is the fraction of the universe of securities (not selected securities, but the entire universe) below zero, divided by the size of the universe, and multiply all weights by 1-CP. Reinvest the remainder in the cash asset (something like VBMFX, VFISX, SHY, etc.). In FAA, I called this the risk-free asset, but in this case, it’s simply the cash asset.

The error term is 1e-6, or some other small value to avoid divide by zero errors.

Finally, wS is an aggression parameter. Setting this value to 0 essentially forces an equal weight portfolio, and infinity will simply select the best asset each month.

Anyhow, let’s look at a prototype of the code (will bug with NA returns, and still doesn’t include excess returns over some treasury security), using the original FAA securities. The weights are wR=1, wC=.5, and wV=0, with wS = 2 for an offensive scheme or wS = .5 and wC = 1 for a more defensive scheme. Crash protection is enabled.

require(quantmod)
require(PerformanceAnalytics)

symbols <- c("VTSMX", "FDIVX", "VEIEX", "VBMFX", "VFISX", "VGSIX", "QRAAX")

getSymbols(symbols, from="1990-01-01")
prices <- list()
for(i in 1:length(symbols)) {
  prices[[i]] <- Ad(get(symbols[i]))  
}
prices <- do.call(cbind, prices)
colnames(prices) <- gsub("\\.[A-z]*", "", colnames(prices))
ep <- endpoints(prices, "months")
prices <- prices[ep,]
prices <- prices["1997-03::"]

EAA <- function(monthlyPrices, wR=1, wV=0, wC=.5, wS=2, errorJitter=1e-6, 
                cashAsset=NULL, bestN=1+ceiling(sqrt(ncol(monthlyPrices))),
                enableCrashProtection = TRUE, returnWeights=FALSE) {
  returns <- Return.calculate(monthlyPrices)
  returns <- returns[-1,] #return calculation uses one observation
  
  if(is.null(cashAsset)) {
    returns$zeroes <- 0
    cashAsset <- "zeroes"
    warning("No cash security specified. Recommended to use one of: quandClean('CHRIS/CME_US'), SHY, or VFISX. 
            Using vector of zeroes instead.")
  }
  
  cashCol <- grep(cashAsset, colnames(returns))
  
  weights <- list()
  for(i in 1:(nrow(returns)-11)) {
    returnsData <- returns[i:(i+11),] #each chunk will be 12 months of returns data
    #per-month mean of cumulative returns of 1, 3, 6, and 12 month periods
    periodReturn <- ((returnsData[12,] + Return.cumulative(returnsData[10:12,]) + 
      Return.cumulative(returnsData[7:12,]) + Return.cumulative(returnsData)))/22
    vols <- StdDev.annualized(returnsData) 
    mktIndex <- xts(rowMeans(returnsData), order.by=index(returnsData)) #equal weight returns of universe
    cors <- cor(returnsData, mktIndex) #correlations to market index
    
    weightedRets <- periodReturn ^ wR
    weightedCors <- (1 - as.numeric(cors)) ^ wC
    weightedVols <- (vols + errorJitter) ^ wV
    wS <- wS + errorJitter
    
    z <- (weightedRets * weightedCors / weightedVols) ^ wS #compute z_i and zero out negative returns
    z[periodReturn < 0] <- 0
    crashProtection <- sum(z==0)/length(z) #compute crash protection cash cushion
    
    orderedZ <- sort(as.numeric(z), decreasing=TRUE)
    selectedSecurities <- z >= orderedZ[bestN]
    preNormalizedWeights <- z*selectedSecurities #select top N securities, keeping z_i scores
    periodWeights <- preNormalizedWeights/sum(preNormalizedWeights) #normalize
    if (enableCrashProtection) {
      periodWeights <- periodWeights * (1-crashProtection) #CP rule
    }
    weights[[i]] <- periodWeights
  }
  
  weights <- do.call(rbind, weights)
  weights[, cashCol] <- weights[, cashCol] + 1-rowSums(weights) #add to risk-free asset all non-invested weight
  strategyReturns <- Return.rebalancing(R = returns, weights = weights) #compute strategy returns
  if(returnWeights) {
    return(list(weights, strategyReturns))
  } else {
    return(strategyReturns)
  }
}

offensive <- EAA(prices, cashAsset="VBMFX", bestN=3)
defensive <- EAA(prices, cashAsset="VBMFX", bestN=3, wS=.5, wC=1)
compare <- cbind(offensive, defensive)
colnames(compare) <- c("Offensive", "Defensive")
stats <- rbind(Return.annualized(compare)*100, StdDev.annualized(compare)*100, maxDrawdown(compare)*100, SharpeRatio.annualized(compare))
rownames(stats)[3] <- "Worst Drawdown"
charts.PerformanceSummary(compare)

And here are the results:

                                Offensive Defensive
Annualized Return               12.085183 10.197450
Annualized Standard Deviation   11.372610  8.633327
Worst Drawdown                  12.629251  8.134785
Annualized Sharpe Ratio (Rf=0%)  1.062657  1.181173

And the resultant equity curves:

Risk and return on full display here. The defensive variant possessing a higher CAGR than max drawdown make it look attractive at first glance, but I’m fairly certain it’s due to investing heavily in bond funds during QE.

Not a bad algorithm, though the fact that there are only 7 securities here leave it open to some idiosyncratic risk. The low-hanging fruit, of course, is that the correlation uses a single-pass variant, meaning that it is quite possible to select correlated assets that are not correlated to non-selected assets, but are indeed correlated to each other. However, since the actual correlations are used here as opposed to correlation rank, I am thinking that a stepwise correlation selection process would have to be specifically written for this particular algorithm. I’ll no doubt run this by David Varadi and see how to properly set up a stepwise correlation algorithm when working with the interplay of correlations with other values (returns, volatility).

In any case, what I like about this algorithm is that not only is it a security *selection* algorithm, which was what FAA was, but it also includes another aspect, which is the actual weighting aspect, rather than simply leaving all assets at equal weight.

At the moment, this is still in its prototypical phase, and I am aware of the bugs of assets that don’t start at the same time, which will be fixed before I commit this to my IKTrading package. But I wanted to put this out there for those that wished to experiment with it and provide feedback if they came across something they believe is worth sharing.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

Why Backtesting On Individual Legs In A Spread Is A BAD Idea

So after reading the last post, the author of quantstrat had mostly critical feedback, mostly of the philosophy that prompted its writing in the first place. Basically, the reason I wrote it, as I stated before, is that I’ve seen many retail users of quantstrat constantly ask “how do I model individual spread instruments”, and otherwise try to look like they’re sophisticated by trading spreads.

The truth is that real professionals use industrial-strength tools to determine their intraday hedge ratios (such a tool is called a spreader). The purpose of quantstrat is not to be an execution modeling system, but to be a *strategy* modeling system. Basically, the purpose of your backtest isn’t to look at individual instruments, since in the last post, the aggregate trade statistics told us absolutely nothing about how our actual spread trading strategy performed. The backtest was a mess as far as the analytics were concerned, and thus rendering it more or less useless. So this post, by request of the author of quantstrat, is about how to do the analysis better, and looking at what matters more–the actual performance of the strategy on the actual relationship being traded–namely, the *spread*, rather than the two components.

So, without further ado, let’s look at the revised code:

require(quantmod)
require(quantstrat)
require(IKTrading)

getSymbols("UNG", from="1990-01-01")
getSymbols("DGAZ", from="1990-01-01")
getSymbols("UGAZ", from="1990-01-01")
UNG <- UNG["2012-02-22::"]
UGAZ <- UGAZ["2012-02-22::"]

spread <- 3*OHLC(UNG) - OHLC(UGAZ)

initDate='1990-01-01'
currency('USD')
Sys.setenv(TZ="UTC")
symbols <- c("spread")
stock(symbols, currency="USD", multiplier=1)

strategy.st <- portfolio.st <- account.st <-"spread_strategy_done_better"
rm.strat(portfolio.st)
rm.strat(strategy.st)
initPortf(portfolio.st, symbols=symbols, initDate=initDate, currency='USD')
initAcct(account.st, portfolios=portfolio.st, initDate=initDate, currency='USD')
initOrders(portfolio.st, initDate=initDate)
strategy(strategy.st, store=TRUE)

#### paramters

nEMA = 20

### indicator

add.indicator(strategy.st, name="EMA",
              arguments=list(x=quote(Cl(mktdata)), n=nEMA),
              label="ema")

### signals

add.signal(strategy.st, name="sigCrossover",
           arguments=list(columns=c("Close", "EMA.ema"), relationship="gt"),
           label="longEntry")

add.signal(strategy.st, name="sigCrossover",
           arguments=list(columns=c("Close", "EMA.ema"), relationship="lt"),
           label="longExit")

### rules

add.rule(strategy.st, name="ruleSignal", 
         arguments=list(sigcol="longEntry", sigval=TRUE, ordertype="market", 
                        orderside="long", replace=FALSE, prefer="Open", orderqty=1), 
         type="enter", path.dep=TRUE)

add.rule(strategy.st, name="ruleSignal", 
         arguments=list(sigcol="longExit", sigval=TRUE, orderqty="all", ordertype="market", 
                        orderside="long", replace=FALSE, prefer="Open"), 
         type="exit", path.dep=TRUE)

#apply strategy
t1 <- Sys.time()
out <- applyStrategy(strategy=strategy.st,portfolios=portfolio.st)
t2 <- Sys.time()
print(t2-t1)

In this case, things are a LOT simpler. Rather than jumping through the hoops of pre-computing an indicator, along with the shenanigans of separate rules for both the long and the short end, we simply have a spread as it’s theoretically supposed to work–three of an unleveraged ETF against the 3x leveraged ETF, and we can go long the spread, or short the spread. In this case, the dynamic seems to be on the up, and we want to capture that.

So how did we do?

#set up analytics
updatePortf(portfolio.st)
dateRange <- time(getPortfolio(portfolio.st)$summary)[-1]
updateAcct(portfolio.st,dateRange)
updateEndEq(account.st)

#trade statistics
tStats <- tradeStats(Portfolios = portfolio.st, use="trades", inclZeroDays=FALSE)
tStats[,4:ncol(tStats)] <- round(tStats[,4:ncol(tStats)], 2)
print(data.frame(t(tStats[,-c(1,2)])))
(aggPF <- sum(tStats$Gross.Profits)/-sum(tStats$Gross.Losses))
(aggCorrect <- mean(tStats$Percent.Positive))
(numTrades <- sum(tStats$Num.Trades))
(meanAvgWLR <- mean(tStats$Avg.WinLoss.Ratio[tStats$Avg.WinLoss.Ratio < Inf], na.rm=TRUE))

And here’s the output:

> print(data.frame(t(tStats[,-c(1,2)])))
                   spread
Num.Txns            76.00
Num.Trades          38.00
Net.Trading.PL       9.87
Avg.Trade.PL         0.26
Med.Trade.PL        -0.10
Largest.Winner       7.76
Largest.Loser       -1.06
Gross.Profits       21.16
Gross.Losses       -11.29
Std.Dev.Trade.PL     1.68
Percent.Positive    39.47
Percent.Negative    60.53
Profit.Factor        1.87
Avg.Win.Trade        1.41
Med.Win.Trade        0.36
Avg.Losing.Trade    -0.49
Med.Losing.Trade    -0.46
Avg.Daily.PL         0.26
Med.Daily.PL        -0.10
Std.Dev.Daily.PL     1.68
Ann.Sharpe           2.45
Max.Drawdown        -4.02
Profit.To.Max.Draw   2.46
Avg.WinLoss.Ratio    2.87
Med.WinLoss.Ratio    0.78
Max.Equity          13.47
Min.Equity          -1.96
End.Equity           9.87
> (aggPF <- sum(tStats$Gross.Profits)/-sum(tStats$Gross.Losses))
[1] 1.874225
> (aggCorrect <- mean(tStats$Percent.Positive))
[1] 39.47
> (numTrades <- sum(tStats$Num.Trades))
[1] 38
> (meanAvgWLR <- mean(tStats$Avg.WinLoss.Ratio[tStats$Avg.WinLoss.Ratio < Inf], na.rm=TRUE))
[1] 2.87

In other words, the typical profile for a trend follower, rather than the uninformative analytics from the last post. Furthermore, the position sizing and equity curve chart actually make sense now. Here they are.

To conclude, while it’s possible to model spreads using individual legs, it makes far more sense in terms of analytics to actually examine the performance of the strategy on the actual relationship being traded, which is the spread itself. Furthermore, after constructing the spread as a synthetic instrument, it can be treated like any other regular instrument in the context of analysis in quantstrat.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

A Way To Model Execution On Individual Legs Of A Spread In Quantstrat

In this post, I’ll attempt to address a question I’ve seen tossed around time and again regarding quantstrat.

“How do I model executions on individual underlying instruments in spread trading?”

First off, a disclaimer: this method is a bit of a kludge, and in using it, you’ll lose out on quantstrat’s inbuilt optimization functionality. Essentially, it builds upon the pre-computed signal methodology I described in a previous post.

Essentially, by appending a column with the same name but with different values to two separate instruments, I can “trick” quantstrat into providing me desired behavior by modeling trading on two underlying instruments.

SO here’s the strategy:

Go long 3 shares of the UNG (natural gas) ETF against 1 share of UGAZ (3x bull) when the spread crosses above its 20-day exponential moving average, otherwise, do nothing. Here’s the reasoning as to why:

require(quantmod)
require(quantstrat)
require(IKTrading)

getSymbols("UNG", from="1990-01-01")
getSymbols("DGAZ", from="1990-01-01")
getSymbols("UGAZ", from="1990-01-01")
UNG <- UNG["2012-02-22::"]
UGAZ <- UGAZ["2012-02-22::"]

spread <- 3*OHLC(UNG) - OHLC(UGAZ)

nEMA=20

chart_Series(spread)
add_TA(EMA(Cl(spread), n=nEMA), on=1, col="blue", lwd=1.5)
legend(x=5, y=50, legend=c("EMA 20"),
       fill=c("blue"), bty="n")

With the corresponding plot:

So, as you can see, we have a spread that drifts upward (something to do with the nature of the leveraged ETF)? So, let’s try and capture that with a strategy.

The way I’m going to do that is to precompute a signal–whether or not the spread’s close is above its EMA20, and append that signal to UNG, with the negative of said signal appended to UGAZ, and then encapsulate it in a quantstrat strategy. In this case, there’s no ATR order sizing function or initial equity–just a simple 3 UNG to 1 UGAZ trade.

signal <- Cl(spread) > EMA(Cl(spread), n=nEMA)
UNG$precomputedSig <- signal
UGAZ$precomputedSig <- signal*-1

initDate='1990-01-01'
currency('USD')
Sys.setenv(TZ="UTC")
symbols <- c("UNG", "UGAZ")
stock(symbols, currency="USD", multiplier=1)

strategy.st <- portfolio.st <- account.st <-"spread_strategy"

rm.strat(portfolio.st)
rm.strat(strategy.st)
initPortf(portfolio.st, symbols=symbols, initDate=initDate, currency='USD')
initAcct(account.st, portfolios=portfolio.st, initDate=initDate, currency='USD')
initOrders(portfolio.st, initDate=initDate)
strategy(strategy.st, store=TRUE)

#long rules
add.signal(strategy.st, name="sigThreshold",
           arguments=list(column="precomputedSig", threshold=.5, 
                          relationship="gt", cross=TRUE),
           label="longEntry")

add.signal(strategy.st, name="sigThreshold",
           arguments=list(column="precomputedSig", threshold=.5, 
                          relationship="lt", cross=TRUE),
           label="longExit")

#short rules
add.signal(strategy.st, name="sigThreshold",
           arguments=list(column="precomputedSig", threshold=-.5, 
                          relationship="lt", cross=TRUE),
           label="shortEntry")

add.signal(strategy.st, name="sigThreshold",
           arguments=list(column="precomputedSig", threshold=-.5, 
                          relationship="gt", cross=TRUE),
           label="shortExit")

#buy 3
add.rule(strategy.st, name="ruleSignal", 
         arguments=list(sigcol="longEntry", sigval=TRUE, ordertype="market", 
                        orderside="long", replace=FALSE, prefer="Open", orderqty=3), 
         type="enter", path.dep=TRUE)

add.rule(strategy.st, name="ruleSignal", 
         arguments=list(sigcol="longExit", sigval=TRUE, orderqty="all", ordertype="market", 
                        orderside="long", replace=FALSE, prefer="Open"), 
         type="exit", path.dep=TRUE)

#short 1
add.rule(strategy.st, name="ruleSignal", 
         arguments=list(sigcol="shortEntry", sigval=TRUE, ordertype="market", 
                        orderside="short", replace=FALSE, prefer="Open", orderqty=-1), 
         type="enter", path.dep=TRUE)

add.rule(strategy.st, name="ruleSignal", 
         arguments=list(sigcol="shortExit", sigval=TRUE, orderqty="all", ordertype="market", 
                        orderside="short", replace=FALSE, prefer="Open"), 
         type="exit", path.dep=TRUE)

t1 <- Sys.time()
out <- applyStrategy(strategy=strategy.st,portfolios=portfolio.st)
t2 <- Sys.time()
print(t2-t1)

#set up analytics
updatePortf(portfolio.st)
dateRange <- time(getPortfolio(portfolio.st)$summary)[-1]
updateAcct(portfolio.st,dateRange)
updateEndEq(account.st)

So, did our spread trade work?

#trade statistics
tStats <- tradeStats(Portfolios = portfolio.st, use="trades", inclZeroDays=FALSE)
tStats[,4:ncol(tStats)] <- round(tStats[,4:ncol(tStats)], 2)
print(data.frame(t(tStats[,-c(1,2)])))
(aggPF <- sum(tStats$Gross.Profits)/-sum(tStats$Gross.Losses))
(aggCorrect <- mean(tStats$Percent.Positive))
(numTrades <- sum(tStats$Num.Trades))
(meanAvgWLR <- mean(tStats$Avg.WinLoss.Ratio[tStats$Avg.WinLoss.Ratio < Inf], na.rm=TRUE))

> (aggPF <- sum(tStats$Gross.Profits)/-sum(tStats$Gross.Losses))
[1] 1.101303
> (aggCorrect <- mean(tStats$Percent.Positive))
[1] 51.315
> (numTrades <- sum(tStats$Num.Trades))
[1] 76
> (meanAvgWLR <- mean(tStats$Avg.WinLoss.Ratio[tStats$Avg.WinLoss.Ratio < Inf], na.rm=TRUE))
[1] 1.065

Sort of. However, when you think about it–looking at the statistics on a per-instrument basis in a spread trade is a bit of a red herring. After all, outside of a small spread, what one instrument makes, another will lose, so the aggregate numbers should be only slightly north of 1 or 50% in most cases, which is what we see here.

A better way of looking at whether or not the strategy performs is to look at the cumulative sum of the daily P&L.

#portfolio cash PL
portString <- paste0("portfolio.", portfolio.st)
portPL <- .blotter[[portString]]$summary$Net.Trading.PL
portPL <- portPL[-1,] #remove initialization date
plot(cumsum(portPL))

With the following equity curve:

Is this the greatest equity curve? Probably not. In fact, after playing around with the strategy a little bit, it’s better to actually get in at the close of the next day than the open (apparently there’s some intraday mean-reversion).

Furthermore, one thing to be careful of is that in this backtest, I made sure that for UNG, my precomputedSig would only take values 1 and 0, and vice versa for the UGAZ variant, such that I could write the rules I did. If it took the values 1, 0, and -1, or 1 and -1, the results would not make sense.

In conclusion, the method I showed was essentially a method building on a previous technique of pre-computing signals. Doing this will disallow users to use quantstrat’s built-in optimization functionality, but will allow users to backtest individual leg execution.

To answer one last question, if one wanted to short the spread as well, the thing to do using this methodology would be to pre-compute a second column called, say, precomputedSig2, that behaved the opposite way.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

The ZOMMA Warthog Index

Harry Long posted another article on SeekingAlpha. As usual, it’s another strategy to be taken with a grain of salt, but very much worth discussing. So here’s the strategy:

Rebalance annually:

50% XLP, 15% GLD, 35% TLT — aka a variation of the stocks/bonds portfolio, just with some gold added in. The twist is the idea that XLP is less risky than SPY. The original article, written by Harry Long, can be found here.

This strategy is about as simple as they come, and is tested just as easily in R.

Here’s the original replication:

require(PerformanceAnalytics)
require(quantmod)
getSymbols("XLP", from="1990-01-01")
getSymbols("GLD", from="1990-01-01")
getSymbols("TLT", from="1990-01-01")
prices <- cbind(Ad(XLP), Ad(GLD), Ad(TLT))
prices <- prices[!is.na(prices[, 2]),]
rets <- Return.calculate(prices)
warthogRets <- Return.portfolio(rets, weights=c(.5, .15, .35), rebalance_on = "years")
getSymbols("SPY", from="1990-01-01")
SPYrets <- Return.calculate(Ad(SPY))
comparison <- merge(warthogRets, SPYrets, join='inner')
charts.PerformanceSummary(comparison)

And the resulting replicated (approximately) equity curve:

So far, so good.

Now, let’s put the claim of outperforming over an entire cycle to the test.

years <- apply.yearly(comparison, Return.cumulative)
years <- years[-1,] #remove 2004
years
sum(years[,1] > years[,2])/nrow(years)
sapply(years, mean)

And here’s the output:

> years
           portfolio.returns SPY.Adjusted
2005-12-30        0.07117683   0.04834811
2006-12-29        0.10846819   0.15843582
2007-12-31        0.14529234   0.05142241
2008-12-31        0.05105269  -0.36791039
2009-12-31        0.03126667   0.26344690
2010-12-31        0.14424559   0.15053339
2011-12-30        0.20384853   0.01897321
2012-12-31        0.07186807   0.15991238
2013-12-31        0.04234810   0.32309145
2014-12-10        0.16070863   0.11534450
> sum(years[,1] > years[,2])/nrow(years)
[1] 0.5
> sapply(years, mean)
portfolio.returns      SPY.Adjusted 
       0.10302756        0.09215978 

So, even with SPY’s 2008 in mind, this strategy seems to break even on a year-to-year basis, and the returns are slightly better, which I suppose should be par for the course for a backtest against SPY through the recession.

Now, let’s just remove that one year and see how things stack up:

years <- years[-4,] #remove 2008
sapply(years, mean)

And we obtain the following results:

> sapply(years, mean)
portfolio.returns      SPY.Adjusted 
        0.1088025         0.1432787 

In short, the so-called outperformance is largely sample-dependent on the S&P having a terrible year. Beyond that, things don’t look so stellar for the strategy.

Lastly, anytime anyone makes a claim about a strategy outperforming a benchmark, one very, very easy test to do is to simply look at the equity curve of a “market neutral” strategy that shorts the benchmark against the strategy, and see if that looks like an equity curve one would be comfortable holding. Here are the results:

diff <- comparison[,1] - comparison[,2]
charts.PerformanceSummary(diff, main="short SPY against portfolio")

And the resulting equity curve of the outperformance:

Basically, ever since the crisis passed, the strategy has been flat to underperforming against SPY, which is fair enough considering that the market has been on a roar since then.

To hammer the point, let’s look at the equity curves from 2009 onwards.

#strategy against SPY post-crisis
charts.PerformanceSummary(comparison["2009::"])

With the result:

To its credit, the strategy’s equity curve is certainly a smoother ride. At the same time, the underperformance is clear as well. I suppose this makes sense as a portfolio invested 35% in bonds should expect to underperform a portfolio fully invested in equities on an absolute returns basis.

Now, as I did before, in It’s Amazing How Well Dumb Things Get Marketed, I’m going to backtest this strategy as far as I can using proxies from Quandl for gold and bonds, by going directly to the futures. Refer to that post to see that the proxies are relatively decent for the ETFs. In this case, I switched from adjusted to close prices, but the concept should still hold.

require(IKTrading)
goldFutures <- quandClean(stemCode = "CHRIS/CME_GC", verbose = TRUE)
thirtyBond <- quandClean(stemCode="CHRIS/CME_US", verbose = TRUE)
longPrices <- cbind(Cl(XLP), Cl(goldFutures), Cl(thirtyBond))
longRets <- Return.calculate(longPrices)
longRets <- longRets[!is.na(longRets[,1]),]
names(longRets) <- c("XLP", "Gold", "Bonds")
longWarthog <- Return.portfolio(longRets, weights=c(.5, .15, .35), rebalance_on = "years")
longComparison <- merge(longWarthog, SPYrets, join='inner')
charts.PerformanceSummary(longComparison)

And here’s the resultant equity curve comparison:

So, the idea that the strategy outperforms SPY on an absolute returns basis? It seems to be due to SPY’s horrendous performance in 2008, which is a once-in-several-decades event. As you can see, if you started from XLP’s inception to approximately 2005, the strategy was actually in drawdown for all those years, at least if you started in 1998.

Let’s look at the risk/reward metrics for this timeframe:

stats <- rbind(Return.annualized(longComparison), 
      maxDrawdown(longComparison),
      SharpeRatio.annualized(longComparison),
      Return.annualized(longComparison)/maxDrawdown(longComparison)
)
rownames(stats)[4] <- "Return Over MaxDD"
stats

And the results:

> stats
                                portfolio.returns SPY.Adjusted
Annualized Return                       0.0409267   0.05344120
Worst Drawdown                          0.2439091   0.55186722
Annualized Sharpe Ratio (Rf=0%)         0.4846184   0.26295594
Return Over MaxDD                       0.1677949   0.09683706

To its credit, the strategy has a significantly better return to risk profile than SPY does, which I hope would be the case for any backtest going up against SPY with 2008 in mind. Overall, the strategy has a bit of merit as a defensive play, and for a simple static allocation based strategy that plays on the defensive side, it delivers what you might expect from a risk-averse strategy, which is reflected in a much better return to drawdown profile than the market. Nevertheless, one should be careful not to misinterpret a strategy whose focus is to control risk for a strategy that aims to swing for the fences. Furthermore, this is a strategy that Harry Long has been so kind as to give away, so take it for what it’s worth–a simple-to-implement defensive strategy which delivers on its defensive promises.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.

An Update to the Robustness Heuristic and a Variation of a Volatility Strategy

So, before revealing a slight wrinkle on the last strategy I wrote about, I’d like to clear up a bit of confusion regarding Jaekle and Tomasini’s idea of a stable region.

Essentially, the entire idea *is* that similar parameter configurations behave in very similar ways, and so, are supposed to be highly correlated. It does not mean the strategy may not be overfit in other ways, but that incremental changes to a parameter should mean incremental changes to performance, rather than seeing some sort of lucky spike in a sea of poor performance.

In any case, the one change to the strategy from last week is that rather than get in at the current close (aka observe close, execute at close), to get in at the next day’s close.

Again, here’s the strategy script:

require(downloader)
require(quantmod)
require(PerformanceAnalytics)
require(TTR)

download("http://www.cboe.com/publish/scheduledtask/mktdata/datahouse/vxvdailyprices.csv", 
         destfile="vxvData.csv")
download("http://www.cboe.com/publish/ScheduledTask/MktData/datahouse/vxmtdailyprices.csv", 
         destfile="vxmtData.csv")

vxv <- xts(read.zoo("vxvData.csv", header=TRUE, sep=",", format="%m/%d/%Y", skip=2))
vxmt <- xts(read.zoo("vxmtData.csv", header=TRUE, sep=",", format="%m/%d/%Y", skip=2))
ratio <- Cl(vxv)/Cl(vxmt)


download("https://dl.dropboxusercontent.com/s/jk6der1s5lxtcfy/XIVlong.TXT",
         destfile="longXIV.txt")

download("https://dl.dropboxusercontent.com/s/950x55x7jtm9x2q/VXXlong.TXT", 
         destfile="longVXX.txt") #requires downloader package

xiv <- xts(read.zoo("longXIV.txt", format="%Y-%m-%d", sep=",", header=TRUE))
vxx <- xts(read.zoo("longVXX.txt", format="%Y-%m-%d", sep=",", header=TRUE))

xiv <- merge(xiv, ratio, join='inner')
vxx <- merge(vxx, ratio, join='inner')
colnames(xiv)[5] <- colnames(vxx)[5] <- "ratio"

xivRets <- Return.calculate(Cl(xiv))
vxxRets <- Return.calculate(Cl(vxx))

retsList <- list()
count <- 1
for(i in 10:200) {
  ratioSMA <- SMA(ratio, n=i)
  vxxSig <- lag(ratio > 1 & ratio > ratioSMA, 2)
  xivSig <- lag(ratio < 1 & ratio < ratioSMA, 2)
  rets <- vxxSig*vxxRets + xivSig*xivRets
  colnames(rets) <- i
  retsList[[i]]  <- rets
  count <- count+1  
}
retsList <- do.call(cbind, retsList)
colnames(retsList) <- gsub("X", "", colnames(retsList))
charts.PerformanceSummary(retsList)
retsList <- retsList[!is.na(retsList[,191]),]
retsList <- retsList[-1,]

The one change I made is that rather than go with the default lag value, I went with a lag of 2. A lag of zero induces look-ahead bias. In any case, let’s run through the process again of analyzing for robustness.

rankComparison <- function(rets, perfAfun="Return.cumulative") {
  fun <- match.fun(perfAfun)
  monthlyFun <- apply.monthly(rets, fun)
  monthlyRank <- t(apply(monthlyFun, MARGIN=1, FUN=rank))
  meanMonthlyRank <- apply(monthlyRank, MARGIN=2, FUN=mean)
  rankMMR <- rank(meanMonthlyRank)
  
  aggFun <- fun(rets)
  aggFunRank <- rank(aggFun)
  
  bothRanks <- data.frame(cbind(aggFunRank, rankMMR, names(rankMMR)), stringsAsFactors=FALSE)
  names(bothRanks) <- c("aggregateRank", "averageMonthlyRank", "configName")
  bothRanks$aggregateRank <- as.numeric(bothRanks$aggregateRank)
  bothRanks$averageMonthlyRank <- as.numeric(bothRanks$averageMonthlyRank)
  bothRanks$sum <- bothRanks[,1] + bothRanks[,2]
  bothRanks <- bothRanks[order(bothRanks$sum, decreasing=TRUE),]
  
  plot(aggFunRank~rankMMR, main=perfAfun)
  print(cor(aggFunRank, meanMonthlyRank))
  return(bothRanks)
}

retRank <- rankComparison(retsList)
sharpeRank <- rankComparison(retsList, perfAfun="SharpeRatio.annualized")

In this case, I added some functionality to not only do the plotting and correlation, but to spit out a table comparing both the aggregate metric along with the rank of the average monthly rank (again, dual ranking layer), and ordered the table by the sum of both the aggregate and the monthly metric, starting with the highest.

For instance, here’s the output from the returns comparison:

> retRank <- rankComparison(retsList)
[1] 0.736377

> head(retRank, 20)
    aggregateRank averageMonthlyRank configName sum
62            190                191         62 381
63            189                187         63 376
60            185                189         60 374
66            191                182         66 373
65            187                183         65 370
59            184                185         59 369
56            181                186         56 367
64            188                179         64 367
152           174                190        152 364
61            183                178         61 361
67            186                167         67 353
151           165                188        151 353
57            179                173         57 352
153           167                184        153 351
58            182                164         58 346
154           170                175        154 345
53            164                180         53 344
155           166                176        155 342
158           163                177        158 340
150           157                181        150 338

So, for this configuration, the correlation went down from above .8 to around .74…which is still strong and credence that the strategy configurations have validity outside some lucky months. The new feature I added was the data frame of the two ranks side by side, along with their configuration name (in this case, my names were simply the SMA parameter, but the names could be anything such as say, SMA_60_lag_2), and the sum of the two rankings, which orders the configurations. As there were 191 configurations (SMA ranging from 10 to 200), the best score that could be achieved was 382. Furthermore, note that although there seems to be a strong region from SMA 53 to SMA 67, there also seems to be another region, at least when it comes to absolute return, of an SMA parameter at SMA 150+.

Here’s the same table for annualized Sharpe (this variation takes a bit longer to compute due to the monthly annualized Sharpes).

> sharpeRank <- rankComparison(retsList, perfAfun="SharpeRatio.annualized")
[1] 0.5590881
> head(sharpeRank, 20)
    aggregateRank averageMonthlyRank configName   sum
62            190              191.0         62 381.0
59            185              190.0         59 375.0
61            183              186.5         61 369.5
60            186              181.0         60 367.0
63            189              175.0         63 364.0
66            191              164.0         66 355.0
152           166              173.0        152 339.0
58            182              155.0         58 337.0
56            181              151.0         56 332.0
53            174              153.0         53 327.0
57            179              148.0         57 327.0
151           159              162.0        151 321.0
76            177              143.0         76 320.0
150           152              163.0        150 315.0
54            173              140.0         54 313.0
77            178              131.0         77 309.0
65            187              119.0         65 306.0
143           146              156.0        143 302.0
74            167              132.0         74 299.0
153           161              138.0        153 299.0

So, largely the same sort of results as we see with the annualized returns. A correlation of .5 gives some cause for concern, which will hopefully show up in the line plot of the rank of the four metrics (returns, Sharpe, drawdowns, and return to drawdown), which will reveal the regions with strong performance, and not-so-strong performances.

Here’s the ranking line plot.

aggReturns <- Return.annualized(retsList)
aggSharpe <- SharpeRatio.annualized(retsList)
aggMAR <- Return.annualized(retsList)/maxDrawdown(retsList)
aggDD <- maxDrawdown(retsList)

plot(rank(aggReturns)~as.numeric(colnames(aggReturns)), type="l", ylab="annualized returns rank", xlab="SMA", 
     main="Risk and return rank comparison")
lines(rank(aggSharpe)~as.numeric(colnames(aggSharpe)), type="l", ylab="annualized Sharpe rank", xlab="SMA", col="blue")
lines(rank(aggMAR)~as.numeric(colnames(aggMAR)), type="l", ylab="Max return over max drawdown", xlab="SMA", col="red")
lines(rank(-aggDD)~as.numeric(colnames(aggDD)), type="l", ylab="max DD", xlab="SMA", col="green")
legend("bottomright", c("Return rank", "Sharpe rank", "MAR rank", "Drawdown rank"), pch=0, col=c("black", "blue", "red", "green"))

And the resulting plot:

There are several regions that show similar, strong metrics for similar parameter choices for the value of SMA when we use a “delayed” entry. Namely, the regions around the 60 day SMA, the 150 day SMA, and the 125 day SMA.

Let’s look at those configurations.

truncRets <- retsList[,c(51, 116, 141)]
stats <- data.frame(cbind(t(Return.annualized(truncRets)),
                 t(SharpeRatio.annualized(truncRets)),
                 t(maxDrawdown(truncRets))))
colnames(stats) <- c("A.Return", "A.Sharpe", "Worst_Drawdown")
stats$MAR <- stats[,1]/stats[,3]
stats <- round(stats, 3)

And the results:

> stats
    A.Return A.Sharpe Worst_Drawdown   MAR
60     1.103    2.490          0.330 3.342
125    0.988    2.220          0.368 2.683
150    0.983    2.189          0.404 2.435

And the resulting performance, on both a regular, and log scale:

charts.PerformanceSummary(truncRets)

logRets <- log(cumprod(1+truncRets))
chart.TimeSeries(logRets)


Perfect strategies? There’s probably room for improvement. As good if not better than the volatility strategies posted elsewhere on the internet? Probably. Is there more investigation that can be done regarding the differences in signal delay? Yes.

So, in conclusion for this post, I’m hoping that the rank comparison heuristic and its new output gives people another tool to consider, along with another vol strategy to consider as well.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.