%\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{Using the SharpeR Package} %\VignetteKeyword{Finance} %\VignetteKeyword{Sharpe} %\VignettePackage{SharpeR} \documentclass[10pt,a4paper,english]{article} % front matter%FOLDUP \usepackage{url} \usepackage{amsmath} \usepackage{amsfonts} \usepackage{hyperref} \usepackage[square,numbers]{natbib} %\usepackage[authoryear]{natbib} %\usepackage[iso]{datetime} %\usepackage{datetime} \makeatletter \makeatother %\input{sr_defs.tex} \usepackage{SharpeR} % knitr setup%FOLDUP <<'preamble', include=FALSE, warning=FALSE, message=FALSE>>= library(knitr) # set the knitr options ... for everyone! # if you unset this, then vignette build bonks. oh, joy. #opts_knit$set(progress=TRUE) opts_knit$set(eval.after='fig.cap') # for a package vignette, you do want to echo. # opts_chunk$set(echo=FALSE,warning=FALSE,message=FALSE) opts_chunk$set(warning=FALSE,message=FALSE) #opts_chunk$set(results="asis") opts_chunk$set(cache=TRUE,cache.path="cache/SharpeR") #opts_chunk$set(fig.path="figure/",dev=c("pdf","cairo_ps")) opts_chunk$set(fig.path="figure/SharpeR",dev=c("pdf")) opts_chunk$set(fig.width=5,fig.height=4,dpi=64) # doing this means that png files are made of figures; # the savings is small, and it looks like shit: #opts_chunk$set(fig.path="figure/",dev=c("png","pdf","cairo_ps")) #opts_chunk$set(fig.width=4,fig.height=4) # for figures? this is sweave-specific? #opts_knit$set(eps=TRUE) # this would be for figures: #opts_chunk$set(out.width='.8\\textwidth') # for text wrapping: options(width=64,digits=2) opts_chunk$set(size="small") opts_chunk$set(tidy=TRUE,tidy.opts=list(width.cutoff=50,keep.blank.line=TRUE)) compile.time <- Sys.time() # from the environment # only recompute if FORCE_RECOMPUTE=True w/out case match. FORCE_RECOMPUTE <- (toupper(Sys.getenv('FORCE_RECOMPUTE',unset='False')) == "TRUE") # compiler flags! # not used yet LONG.FORM <- FALSE library(quantmod) options("getSymbols.warning4.0"=FALSE) @ %UNFOLD %UNFOLD % document incantations%FOLDUP \begin{document} \title{Using SharpeR} \author{Steven E. Pav % \thanks{\email{shabbychef@gmail.com}}} %\date{\today, \currenttime} \maketitle %UNFOLD \begin{abstract}%FOLDUP The SharpeR package provides basic functionality for testing significance of the \txtSR of a series of returns, and of the Markowitz portfolio on a number of possibly correlated assets.\cite{Sharpe:1966} The goal of the package is to make it simple to estimate profitability (in terms of risk-adjusted returns) of strategies or asset streams. \end{abstract}%UNFOLD \section{The \txtSR and Optimal \txtSR}%FOLDUP Sharpe defined the `reward to variability ratio', now known as the `\txtSR', as the sample statistic $$ \ssr = \frac{\smu}{\ssig}, $$ where \smu is the sample mean, and \ssig is the sample standard deviation. \cite{Sharpe:1966} The \txtSR was later redefined to include a `risk-free' or `disastrous rate of return': $\ssr = \wrapParens{\smu - \rfr}/\ssig.$ It is little appreciated in quantitative finance that the \txtSR is identical to the sample statistic proposed by Gosset in 1908 to test for zero mean when the variance is unknown. \cite{student08ttest} The `\tstat-test' we know today, which includes an adjustment for sample size, was formulated later by Fisher. \cite{Fisher1925} Knowing that the \txtSR is related to the \tstat-statistic provides a `natural arbitrage,' since the latter has been extensively studied. Many of the interesting properties of the \tstat-statistic can be translated to properties about the \txtSR. Also little appreciated is that the multivariate analogue of the \tstat-statistic, Hotelling's \Tstat, is related to the Markowitz portfolio. Consider the following portfolio optimization problem: \begin{equation} \max_{\sportw : \qform{\svsig}{\sportw} \le R^2} \frac{\trAB{\sportw}{\svmu} - \rfr}{\sqrt{\qform{\svsig}{\sportw}}}, \label{eqn:port_prob} \end{equation} where \svmu, \svsig are the sample mean vector and covariance matrix, \rfr is the risk-free rate, and $R$ is a cap on portfolio `risk' as estimated by \svsig. (Note this differs from the traditional definition of the problem which imposes a `self-financing constraint' which does not actually bound portfolio weights.) The solution to this problem is $$ \sportwopt \defeq \frac{R}{\sqrt{\qform{\minv{\svsig}}{\svmu}}} \minv{\svsig}\svmu. $$ The \txtSR of this portfolio is \begin{equation} \ssropt \defeq \frac{\trAB{\sportwopt}{\svmu} - \rfr}{\sqrt{\qform{\svsig}{\sportwopt}}} = \sqrt{\qform{\minv{\svsig}}{\svmu}} - \frac{\rfr}{R} = \sqrt{\Tstat / \ssiz} - \frac{\rfr}{R}, \label{eqn:ssr_uncons} \end{equation} where \Tstat is Hotelling's statistic, and \ssiz is the number of independent observations (\eg `days') used to construct \svmu. The term $\rfr / R$ is a deterministic `drag' term that merely shifts the location of \ssropt, and so we can (mostly) ignore it when testing significance of \ssropt. Under the (typically indefensible) assumptions that the returns are generated \iid from a normal distribution (multivariate normal in the case of the portfolio problem), the distributions of \ssr and \ssrsqopt are known, and depend on the sample size and the population analogues, \psr and \psnrsqopt. In particular, they are distributed as rescaled non-central \tlaw{} and \flaw{} distributions. Under these assumptions on the generating processes, we can perform inference on the population analogues using the sample statistics. The importance of each of these assumptions (\viz homoskedasticity, independence, normality, \etc) can and should be checked. \cite{Lumley:2002,Opdyke2007} The reader must be warned that this package is distributed without any warranty of any kind, and in no way should any analysis performed with this package be interpreted as implicit investment advice by the author(s). The units of \smu are `returns per time,' while those of \ssig are `returns per square root time.' Consequently, the units of \ssr are `per square root time.' Typically the \txtSR is quoted in `annualized' terms, \ie \yrto{-1/2}, but the units are omitted. I believe that units should be included as it avoids ambiguity, and simplifies conversions. There is no clear standard whether arithmetic or geometric returns should be used in the computation of the \txtSR. Since arithmetic returns are always greater than the equivalent geometric returns, one would suspect that arithmetic returns are \emph{always} used when advertising products. However, I suspect that geometric returns are more frequently used in the analysis of strategies. Geometric returns have the attractive property of being `additive', meaning that the geometric return of a period is the sum of those of subperiods, and thus the sign of the arithmetic mean of some geometric returns indicates whether the final value of a strategy is greater than the initial value. Oddly, the arithmetic mean of arithmetic returns does not share this property. On the other hand, arithmetic returns are indeed additive \emph{contemporaneously}: if \vreti is the vector of arithmetic returns of several stocks, and \sportw is the dollar proportional allocation into those stocks at the start of the period, then \trAB{\vreti}{\sportw} is the arithmetic return of the portfolio over that period. This holds even when the portfolio holds some stocks `short.' Often this portfolio accounting is misapplied to geometric returns without even an appeal to Taylor's theorem. For more details on the \txtSR and portfolio optimization, see the vignette, ``Notes on the \txtSR'' distributed with this package. %UNFOLD \section{Using the \Robject{sr} Class}%FOLDUP An \Robject{sr} object encapsulates one or more \txtSR statistics, along with the degrees of freedom, the rescaling to a \tstat statistic, and the annualization and units information. One can simply stuff this information into an \Robject{sr} object, but it is more straightforward to allow \Rfunction{as.sr} to compute the \txtSR for you. <<'babysteps'>>= library(SharpeR) # suppose you computed the Sharpe of your strategy to # be 1.3 / sqrt(yr), based on 1200 daily observations. # store them as follows: my.sr <- sr(sr=1.3,df=1200-1,ope=252,epoch="yr") print(my.sr) # multiple strategies can be tracked as well. # one can attach names to them. srstats <- c(0.5,1.2,0.6) dim(srstats) <- c(3,1) rownames(srstats) <- c("strat. A","strat. B","benchmark") my.sr <- sr(srstats,df=1200-1,ope=252,epoch="yr") print(my.sr) @ Throughout, \Rcode{ope} stands for `Observations Per Epoch', and is the (average) number of returns observed per the annualization period, called the \Rcode{epoch}. At the moment there is not much hand holding regarding these parameters: no checking is performed for sane values. The \Rfunction{as.sr} method will compute the \txtSR for you, from numeric, \Robject{data.frame}, \Robject{xts} or \Robject{lm} objects. In the latter case, it is assumed one is performing an attribution model, and the statistic of interest is the fit of the \Rcode{(Intercept)} term divided by the residual standard deviation. Here are some examples: <<'showoff'>>= set.seed(as.integer(charToRaw("set the seed"))) # Sharpe's 'model': just given a bunch of returns. returns <- rnorm(253*8,mean=3e-4,sd=1e-2) asr <- as.sr(returns,ope=253,epoch="yr") print(asr) # a data.frame with a single strategy asr <- as.sr(data.frame(my.strategy=returns),ope=253,epoch="yr") print(asr) @ When a \Robject{data.frame} with multiple columns is given, the \txtSR of each is computed, and they are all stored: <<'more_data_frame'>>= # a data.frame with multiple strategies asr <- as.sr(data.frame(strat1=rnorm(253*8),strat2=rnorm(253*8,mean=4e-4,sd=1e-2)), ope=253,epoch="yr") print(asr) @ Here is an example using \Robject{xts} objects. In this case, if the \Rcode{ope} is not given, it is inferred from the time marks of the input object. % MOCK it up. <<'stock_loading',eval=FALSE,echo=TRUE>>= require(quantmod) # get price data, compute log returns on adjusted closes get.ret <- function(sym,warnings=FALSE,...) { # getSymbols.yahoo will barf sometimes; do a trycatch trynum <- 0 while (!exists("OHCLV") && (trynum < 7)) { trynum <- trynum + 1 try(OHLCV <- getSymbols(sym,auto.assign=FALSE,warnings=warnings,...),silent=TRUE) } adj.names <- paste(c(sym,"Adjusted"),collapse=".",sep="") if (adj.names %in% colnames(OHLCV)) { adj.close <- OHLCV[,adj.names] } else { # for DJIA from FRED, say. adj.close <- OHLCV[,sym] } rm(OHLCV) # rename it colnames(adj.close) <- c(sym) adj.close <- adj.close[!is.na(adj.close)] lrets <- diff(log(adj.close)) #chop first lrets[-1,] } get.rets <- function(syms,...) { some.rets <- do.call("cbind",lapply(syms,get.ret,...)) } @ <<'stock_loading_sneaky',eval=TRUE,echo=FALSE>>= # sleight of hand to load precomputed data instead. get.rets <- function(syms,from='2003-01-01',to='2013-01-01',...) { fname <- system.file('extdata','ret_data.rda',package='SharpeR') if (fname == "") { fname <- 'ret_data.rda' } # poofs all.ret here load(fname) sub.data <- all.ret[paste(from,to,sep="::"),colnames(all.ret) %in% syms] return(sub.data) } @ <<'helper_function'>>= require(quantmod) # quantmod::periodReturn does not deal properly with multiple # columns, and the straightforward apply(mtms,2,periodReturn) barfs my.periodReturn <- function(mtms,...) { per.rets <- do.call(cbind,lapply(mtms, function(x) { retv <- periodReturn(x,...) colnames(retv) <- colnames(x) return(retv) })) } # convert log return to mtm, ignoring NA lr2mtm <- function(x,...) { x[is.na(x)] = 0 exp(cumsum(x)) } @ <<'some_stocks'>>= some.rets <- get.rets(c("IBM","AAPL","XOM"), from="2007-01-01",to="2013-01-01") print(as.sr(some.rets)) @ The annualization of an \Robject{sr} object can be changed with the \Rfunction{reannualize} method. The name of the epoch and the observation rate can both be changed. Changing the annualization will not change statistical significance, it merely changes the units. <<'reannualize'>>= yearly <- as.sr(some.rets[,"XOM"]) monthly <- reannualize(yearly,new.ope=21,new.epoch="mo.") print(yearly) # significance should be the same, but units changed. print(monthly) @ \subsection{Attribution Models}%FOLDUP When an object of class \Robject{lm} is given to \Rfunction{as.sr}, the fit \Rcode{(Intercept)} term is divided by the residual volatility to compute something like the \txtSR. In terms of Null Hypothesis Significance Testing, nothing is gained by summarizing the \Robject{sr} object instead of the \Robject{lm} object. However, confidence intervals on the \txtSR are quoted in the more natural units of reward to variability, and in annualized terms (or whatever the epoch is.) As an example, here I perform a CAPM attribution to the monthly returns of \StockTicker{AAPL}. Note that the statistical significance here is certainly tainted by selection bias, a topic beyond the scope of this note. <<'AAPL'>>= # get the returns (see above for the function) aapl.rets <- get.rets(c("AAPL","SPY"),from="2003-01-01",to="2013-01-01") # make them monthly: mo.rets <- my.periodReturn(lr2mtm(aapl.rets),period='monthly',type='arithmetic') rm(aapl.rets) # cleanup # look at both of them together: both.sr <- as.sr(mo.rets) print(both.sr) # confindence intervals on the Sharpe: print(confint(both.sr)) # perform a CAPM attribution, using SPY as 'the market' linmod <- lm(AAPL ~ SPY,data=mo.rets) # convert attribution model to Sharpe CAPM.sr <- as.sr(linmod,ope=both.sr$ope,epoch="yr") # statistical significance does not change (though note the sr summary # prints a 1-sided p-value) print(summary(linmod)) print(CAPM.sr) # the confidence intervals tell the same story, but in different units: print(confint(linmod,'(Intercept)')) print(confint(CAPM.sr)) @ %UNFOLD \subsection{Testing Sharpe and Power}%FOLDUP The function \Rfunction{sr\_test} performs one- and two-sample tests for significance of \txtSR. Paired tests for equality of \txtSR can be performed via the \Rfunction{sr\_equality\_test}, which applies the tests of Leung \etal or of Wright \etal \cite{Leung2008,Wright2012} <<'SPDRcheck'>>= # get the sector 'spiders' secto.rets <- get.rets(c("XLY","XLE","XLP","XLF","XLV","XLI","XLB","XLK","XLU"), from="2003-01-01",to="2013-01-01") # make them monthly: mo.rets <- my.periodReturn(lr2mtm(secto.rets),period='monthly',type='arithmetic') # one-sample test on utilities: XLU.monthly <- mo.rets[,"XLU"] print(sr_test(XLU.monthly),alternative="two.sided") # test for equality of Sharpe among the different spiders print(sr_equality_test(as.matrix(secto.rets))) # perform a paired two-sample test via sr_test: XLF.monthly <- mo.rets[,"XLF"] print(sr_test(x=XLU.monthly,y=XLF.monthly,ope=12,paired=TRUE)) @ %UNFOLD %UNFOLD \section{Using the \Robject{sropt} Class}%FOLDUP The class \Robject{sropt} stores the `optimal' \txtSR, which is that of the optimal (`Markowitz') portfolio, as defined in \eqnref{port_prob}, as well as the relevant degrees of freedom, and the annualization parameters. Again, the constructor can be used directly, but the helper function is preferred: <<'sropt_basics'>>= set.seed(as.integer(charToRaw("7bf4b86a-1834-4b58-9eff-6c7dec724fec"))) # from a matrix object: ope <- 253 n.stok <- 7 n.yr <- 8 # somewhat unrealistic: independent returns. rand.rets <- matrix(rnorm(n.yr * ope * n.stok),ncol=n.stok) asro <- as.sropt(rand.rets,ope=ope) rm(rand.rets) print(asro) # under the alternative, when the mean is nonzero rand.rets <- matrix(rnorm(n.yr * ope * n.stok,mean=6e-4,sd=1e-2),ncol=n.stok) asro <- as.sropt(rand.rets,ope=ope) rm(rand.rets) print(asro) # from an xts object some.rets <- get.rets(c("IBM","AAPL","XOM"), from="2007-01-01",to="2013-01-01") asro <- as.sropt(some.rets) print(asro) @ One can compute confidence intervals for the population parameter $\psnropt \defeq \sqrt{\qform{\minv{\pvsig}}{\pvmu}}$, called the `optimal signal-noise ratio', based on inverting the non-central \flaw{}-distribution. Estimates of \psnropt can also be computed, via Maximum Likelihood Estimation, or a `shrinkage' estimate. \cite{kubokawa1993estimation,MC1986216} <<'sropt_estim'>>= # confidence intervals: print(confint(asro,level.lo=0.05,level.hi=1)) # estimation print(inference(asro,type="KRS")) print(inference(asro,type="MLE")) @ A nice rule of thumb is that, to a first order approximation, the MLE of \psnropt is zero exactly when $\ssrsqopt \le \nlatf/\ssiz,$ where \nlatf is the number of assets. \cite{kubokawa1993estimation,MC1986216} Inspection of this inequality confirms that \ssropt and \ssiz can be expressed `in the same units', meaning that if \ssropt is in \yrto{-1/2}, then \ssiz should be the number of \emph{years}. For example, if the Markowitz portfolio on 8 assets over 7 years has a \txtSR of 1\yrto{-1/2}, the MLE will be zero. This can be confirmed empirically as below. <<'MLE_rule'>>= ope <- 253 zeta.s <- 0.8 n.check <- 1000 df1 <- 10 df2 <- 6 * ope rvs <- rsropt(n.check,df1,df2,zeta.s,ope,drag=0) roll.own <- sropt(z.s=rvs,df1,df2,drag=0,ope=ope,epoch="yr") MLEs <- inference(roll.own,type="MLE") zerMLE <- MLEs <= 0 crit.value <- 0.5 * (max(rvs[zerMLE])^2 + min(rvs[!zerMLE])^2) aspect.ratio <- df1 / (df2 / ope) cat(sprintf("empirical cutoff for zero MLE is %2.2f yr^{-1}\n", crit.value)) cat(sprintf("the aspect ratio is %2.2f yr^{-1}\n",aspect.ratio)) @ \subsection{Approximating Overfit}%FOLDUP A high-level sketch of quant work is as follows: construct a trading system with some free parameters, \stratrc, backtest the strategy for $\stratrc[1],\stratrc[2],\ldots,\stratrc[m]$, then pick the \stratrc[i] that maximizes the \txtSR in backtesting. `Overfit bias' (variously known as `datamining bias,' `garbatrage,' `backtest arb,' \etc) is the upward bias in one's estimate of the true signal-to-noise of the strategy parametrized by \stratrc[i^*] due to one using the same data to select the strategy and estimate it's performance. \cite{Aronson2007}[Chapter 6] As an example, consider a basic Moving Average Crossover strategy. Here \stratrc is a vector of 2 window lengths. One longs the market exactly when one moving average exceeds the other, and shorts otherwise. One performs a brute force search of all allowable window sizes. Before deciding to deploy money into the MAC strategy parametrized by \stratrc[i^*], one has to estimate its profitability. There is a way to roughly estimate overfit bias by viewing the problem as a portfolio problem, and performing inference on the optimal \txtSR. To do this, suppose that \vreti[i] is the \ssiz-vector of backtested returns associated with \stratrc[i]. (I will assume that all the backtests are over the same period of time.) Then approximately embed the backtested returns vectors in the subset of a \nlatf-dimensional subspace. That is, by a process like PCA, make the approximation: %A(?): Make PCA-like linear approximation of returns vectors: %$$\wrapNeBraces{\vreti[1],\ldots,\vreti[m]} \approx \mathcal{L} \subset \setwo{\sum_{1\le j \le \nlatf} k_j \vretj[j]}{k_j \in \reals}$$ $$\wrapNeBraces{\vreti[1],\ldots,\vreti[m]} \approx \mathcal{K} \subset \mathcal{L} \defeq \setwo{\mretj \sportw}{\sportw \in \reals{\nlatf}}$$ Abusing notation, let $\funcit{\ssr}{\stratrc}$ be the sample \txtSR associated with parameters \stratrc, and also let $\funcit{\ssr}{\vreti}$ be the \txtSR associated with the vector of returns \vreti. Then make the approximation $$ \funcit{\ssr}{\stratrc[*]} \defeq \max_{\stratrc[1],\ldots,\stratrc[m]} \funcit{\ssr}{\stratrc[i]} \approx \max_{\vreti \in \mathcal{K}} \funcit{\ssr}{\vreti} \le \max_{\vreti \in \mathcal{L}} \funcit{\ssr}{\vreti} = \ssropt. $$ This is a conservative approximation: the true maximum over $\mathcal{L}$ is presumably much larger than \funcit{\ssr}{\stratrc[*]}. One can then use \funcit{\ssr}{\stratrc[*]} as \ssropt over a set of \nlatf assets, perform inference on \psnropt, which, by a series of approximations as above, is an approximate upper bound on \funcit{\psnr}{\stratrc[*]}. This approximate attack on overfitting will work better when one has a good estimate of \nlatf, when $m$ is relatively large and \nlatf relatively small, and when the linear approximation to the set of backtested returns is good. Moreover, the definition of $\mathcal{L}$ explicitly allows shorting, whereas the backtested returns vectors $\vreti[i]$ may lack the symmetry about zero to make this a good approximation. By way of illustration, consider the case where the trading system is set up such that different \stratrc produce minor variants on a clearly losing strategy: in this case we might have $\funcit{\ssr}{\stratrc[*]} < 0$, which cannot hold for \ssropt. One can estimate \nlatf via Monte Carlo simulations, by actually performing PCA, or via the `SWAG' method. Surprisingly, often one's intuitive estimate of the true `degrees of freedom' in a trading system is reasonably good. <<'estimate_overfit',fig.cap=paste("Q-Q plot of",n.sim,"achieved optimal \\txtSR values from brute force search over both windows of a Moving Average Crossover under the null of driftless log returns with zero autocorrelation versus the approximation by a 2-parameter optimal \\txtSR distribution is shown.")>>= require(TTR) # brute force search two window MAC brute.force <- function(lrets,rrets=exp(lrets)-1,win1,win2=win1) { mtms <- c(1,exp(cumsum(lrets))) # prepend a 1. # do all the SMAs; SMA1 <- sapply(win1,function(n) { SMA(mtms,n=n) }) symmetric <- missing(win2) if (!symmetric) SMA2 <- sapply(win2,function(n) { SMA(mtms,n=n) }) mwin <- max(c(win1,win2)) zeds <- matrix(NaN,nrow=length(win1),ncol=length(win2)) upb <- if (symmetric) length(win1) - 1 else length(win1) # 2FIX: vectorize this! for (iidx in 1:upb) { SM1 <- SMA1[,iidx] lob <- if (symmetric) iidx + 1 else 1 for (jidx in lob:length(win2)) { SM2 <- if (symmetric) SMA1[,jidx] else SMA2[,jidx] trades <- sign(SM1 - SM2) dum.bt <- trades[mwin:(length(trades)-1)] * rrets[mwin:length(rrets)] # braindead backtest. mysr <- as.sr(dum.bt) zeds[iidx,jidx] <- mysr$sr # abuse symmetry of arithmetic returns if (symmetric) zeds[jidx,iidx] <- - zeds[iidx,jidx] } } retv <- max(zeds,na.rm=TRUE) return(retv) } # simulate one. sim.one <- function(nbt,win1,...) { lrets <- rnorm(nbt+max(win1),sd=0.01) retv <- brute.force(lrets,win1=win1,...) return(retv) } # set everything up set.seed(as.integer(charToRaw("e23769f4-94f8-4c36-bca1-28c48c49b4fb"))) ope <- 253 n.yr <- 4 n.obs <- ceiling(ope * n.yr) LONG.FORM <- FALSE n.sim <- if (LONG.FORM) 2048 else 1024 win1 <- if (LONG.FORM) c(2,4,8,16,32,64,128,256) else c(4,16,64,256) # run them system.time(max.zeds <- replicate(n.sim,sim.one(n.obs,win1))) # qqplot; qqplot(qsropt(ppoints(length(max.zeds)),df1=2,df2=n.obs),max.zeds, xlab = "Theoretical Approximate Quantiles", ylab = "Sample Quantiles") qqline(max.zeds,datax=FALSE,distribution = function(p) { qsropt(p,df1=2,df2=n.obs) }, col=2) @ Here I illustrate the quality of the approximation for the two-window simple MAC strategy. I generate log returns which are homoskedastic, driftless, and have zero autocorrelation. In this case, \emph{every} MAC strategy has zero expected return (ignoring trading costs). In spite of this deficiency in the market, I find the best combination of window sizes by looking at \Sexpr{n.yr} years of daily data. By selecting the combination of windows with the highest \txtSR, then using that maximal value as an estimate of the selected model's true signal-noise-ratio, I have subjected myself to overfit bias. I repeat this experiment \Sexpr{n.sim} times, then Q-Q plot the maximal \txtSR values over those experiments versus an optimal \txtSR distribution assuming $\nlatf=2$ in \figref{estimate_overfit}. The fit is reasonable\footnote{And much better when the overfitting is more aggressive, which takes more processing time to simulate.} except in the case where the maximal in-sample \txtSR is very low (recall that it can be negative for this brute-force search, whereas the optimal \txtSR distribution does not produce negative values). This case is unlikely to lead to a trading catastrophe, however. It behooves me to replicate the above experiment `under the alternative,' \eg when the market has autocorrelated returns, to see if the approximation holds up when $\psnropt > 0$. I leave this for future iterations. Instead, I apply the $\nlatf=2$ approximation to the brute-force MAC overfit on \StockTicker{SPY}. <<'now_on_spy'>>= # is MAC on SPY significant? SPY.lret <- get.rets(c('SPY'),from="2003-01-01",to="2013-01-01") # oops! there might be NAs in there! mysr <- as.sr(SPY.lret,na.rm=TRUE) # just to get the ope print(mysr) # get rid of NAs SPY.lret[is.na(SPY.lret)] <- 0 # try a whole lot of windows: win1 <- seq(4,204,by=10) zeds <- brute.force(SPY.lret,win1=win1) SPY.MAC.asro <- sropt(z.s=zeds,df1=2,df2=length(SPY.lret) - max(win1),ope=mysr$ope) print(SPY.MAC.asro) print(inference(SPY.MAC.asro,type="KRS")) @ The \txtSR for \StockTicker{SPY} over this period is \Sexpr{mysr$sr}\yrto{-1/2}. The optimal \txtSR for the tested MAC strategies is \Sexpr{SPY.MAC.asro$sropt}\yrto{-1/2}. The KRS estimate (using the $\nlatf=2$ approximation) for the upper bound on signal-to-noise of the optimal MAC strategy is only \Sexpr{inference(SPY.MAC.asro,type="KRS")}\yrto{-1/2}. \cite{kubokawa1993estimation} This leaves little room for excitement about MAC strategies on \StockTicker{SPY}. %UNFOLD %UNFOLD \section{Using the \Robject{del\_sropt} Class}%FOLDUP The class \Robject{del\_sropt} stores the `optimal' \txtSR of the optimal hedge-constrained portfolio. See the `Notes on \txtSR' vignette for more details. There is an object constructor, but likely the \Rcode{as.del\_sropt} function will prove more useful. <<'del_sropt_basics'>>= set.seed(as.integer(charToRaw("364e72ab-1570-43bf-a1c6-ee7481e1c631"))) # from a matrix object: ope <- 253 n.stok <- 7 n.yr <- 8 # somewhat unrealistic: independent returns, under the null rand.rets <- matrix(rnorm(n.yr * ope * n.stok),ncol=n.stok) # the hedge constraint: hedge out the first stock. G <- diag(n.stok)[1,] asro <- as.del_sropt(rand.rets,G,ope=ope) print(asro) # hedge out the first two G <- diag(n.stok)[1:2,] asro <- as.del_sropt(rand.rets,G,ope=ope) print(asro) # under the alternative, when the mean is nonzero rand.rets <- matrix(rnorm(n.yr * ope * n.stok,mean=6e-4,sd=1e-2),ncol=n.stok) G <- diag(n.stok)[1,] asro <- as.del_sropt(rand.rets,G,ope=ope) print(asro) @ Here I present an example of hedging out \StockTicker{SPY} from a portfolio holding \StockTicker{IBM}, \StockTicker{AAPL}, and \StockTicker{XOM}. The 95\% confidence interval on the optimal hedged signal-noise ratio contains zero: <<'del_sropt_hedging'>>= # from an xts object some.rets <- get.rets(c("SPY","IBM","AAPL","XOM"), from="2007-01-01",to="2013-01-01") # without the hedge, allowing SPY position asro <- as.sropt(some.rets) print(asro) # hedge out SPY! G <- diag(dim(some.rets)[2])[1,] asro.hej <- as.del_sropt(some.rets,G) print(asro.hej) @ One can compute confidence intervals for, and perform inference on, the population parameter, $\sqrt{\psnrsqoptG{\eye} - \psnrsqoptG{\hejG}}$. This is the optimal signal-noise ratio of a hedged portfolio. These are all reported here in the same `annualized' units of \txtSR, and may be controlled by the \Rcode{ope} parameter. <<'del_sropt_estim'>>= # confidence intervals: print(confint(asro,level.lo=0.05,level.hi=1)) # estimation print(inference(asro,type="KRS")) print(inference(asro,type="MLE")) @ %UNFOLD \section{Hypothesis Tests}%FOLDUP The function \Rfunction{sr\_equality\_test}, implements the tests of Leung \etal and of Wright \etal for testing the equality of \txtSR. \cite{Leung2008,Wright2012} I have found that these tests can reject the null `for the wrong reason' when the problem is ill-scaled. Here I confirm that the tests give approximately uniform p-values under the null in three situations: when stock returns are independent Gaussians, when they are independent and \tlaw{}-distributed, and when they are Gaussian and correlated. Visual confirmation is via \figref{sr_eq_test_normal} through \figref{sr_eq_test_correlated}, showing Q-Q plots against uniformity of the putative p-values in Monte Carlo simulations with data drawn from the null. <<'sr_eq_test_normal',fig.cap=paste("P-values from \\Rfunction{sr\\_equality\\_test} on",n.sim,"multiple realizations of independent, normally distributed returns series are Q-Q plotted against uniformity.")>>= # function for Q-Q plot against uniformity qqunif <- function(x,xlab="Theoretical Quantiles under Uniformity", ylab=NULL,...) { if (is.null(ylab)) ylab=paste("Sample Quantiles (",deparse(substitute(x)),")",sep="") qqplot(qunif(ppoints(length(x))),x,xlab=xlab,ylab=ylab,...) abline(0,1,col='red') } # under normality. LONG.FORM <- FALSE n.ro <- if (LONG.FORM) 253*4 else 253*2 n.co <- if (LONG.FORM) 20 else 4 n.sim <- if (LONG.FORM) 1024 else 512 set.seed(as.integer(charToRaw("e3709e11-37e0-449b-bcdf-9271fb1666e5"))) afoo <- replicate(n.sim,sr_equality_test(matrix(rnorm(n.ro*n.co),ncol=n.co),type="F")) sr_eq.pvals <- unlist(afoo["p.value",]) qqunif(sr_eq.pvals) @ <<'sr_eq_test_hetero',fig.cap=paste("P-values from \\Rfunction{sr\\_equality\\_test} on",n.sim,"multiple realizations of independent, \\tlaw{}-distributed returns series are Q-Q plotted against uniformity.")>>= # try heteroskedasticity? set.seed(as.integer(charToRaw("81c97c5e-7b21-4672-8140-bd01d98d1d2e"))) afoo <- replicate(n.sim,sr_equality_test(matrix(rt(n.ro*n.co,df=4),ncol=n.co),type="F")) sr_eq.pvals <- unlist(afoo["p.value",]) qqunif(sr_eq.pvals) @ <<'sr_eq_test_correlated',fig.cap=paste("P-values from \\Rfunction{sr\\_equality\\_test} on",n.sim,"multiple realizations of correlated, normally-distributed returns series are Q-Q plotted against uniformity.")>>= # try correlated returns n.fact <- max(2,n.co - 5) gen.ret <- function(n1,n2,f=max(2,n2-2),fuzz=0.1) { A <- matrix(rnorm(n1*f),nrow=n1) B <- matrix(rnorm(f*n2),nrow=f) C <- sqrt(1-fuzz^2) * A %*% B + fuzz * matrix(rnorm(n1*n2),nrow=n1) } set.seed(as.integer(charToRaw("e4d61c2c-efb3-4cba-9a6e-5f5276ce2ded"))) afoo <- replicate(n.sim,sr_equality_test(gen.ret(n.ro,n.co,n.fact),type="F")) sr_eq.pvals <- unlist(afoo["p.value",]) qqunif(sr_eq.pvals) @ <<'sr_eq_test_mtime_pre',echo=FALSE>>= n.co <- if (LONG.FORM) 100 else 50 n.sim <- if (LONG.FORM) 1024 else 128 @ As a followup, I consider the case where the returns are those of several randomly generated market-timing strategies which are always fully invested long or short in the market, with equal probability. I look at \Sexpr{n.sim} simulations of \Sexpr{n.co} random market timing strategies rebalancing weekly on a 10 year history of \StockTicker{SPY}. The 'p-values' from \Rfunction{sr\_equality\_test} applied to the returns are Q-Q plotted against uniformity in \figref{sr_eq_test_mtime}. The type I rate of this test is far larger than the nominal rate: it is rejecting for `uninteresting reasons'. I suspect this test breaks down because of small sample size or small 'aspect ratio' (ratio of weeks to strategies in this example). In summary, one should exercise caution when interpreting the results of the Sharpe equality test: it would be worthwhile to compare the results against a `Monte Carlo null' as illustrated here. <<'sr_eq_test_mtime',fig.cap=paste("P-values from \\Rfunction{sr\\_equality\\_test} on",n.sim,"multiple realizations of",n.co,"market timing strategies' returns series are Q-Q plotted against uniformity. Nominal coverage is not maintained: the test is far too liberal in this case.")>>= n.co <- if (LONG.FORM) 100 else 50 n.sim <- if (LONG.FORM) 1024 else 128 SPY.lret <- get.rets(c('SPY'),from="2003-01-01",to="2013-01-01") SPY.wk.rret <- my.periodReturn(lr2mtm(SPY.lret),period='weekly', type='arithmetic') gen.tim <- function(n2) { mkt.timing.signal <- sign(rnorm(n2*length(SPY.wk.rret))) mkt.ret <- matrix(rep(SPY.wk.rret,n2) * mkt.timing.signal,ncol=n2) } set.seed(as.integer(charToRaw("447cfe85-b612-4b14-bd01-404e6e99aca4"))) system.time(afoo <- replicate(n.sim,sr_equality_test(gen.tim(n.co),type="F"))) sr_eq.pvals <- unlist(afoo["p.value",]) qqunif(sr_eq.pvals) @ \subsection{Equality Test with Fancy Covariance} The \Rcode{sr\_equality\_test} function now accepts an optional function to perform the inner variance-covariance estimation. I do not think this will correct the problems noted previously for the ill-scaled case. It does, however, allow one to take into account \eg heteroskedasticity and autocorrelation, as follows. Note that the use of a `fancy' covariance estimator does not change the conclusions of the Sharpe equality test in this case. <<'sr_fancy_eq'>>= # get returns some.rets <- as.matrix(get.rets(c("IBM","AAPL","XOM"), from="2007-01-01",to="2013-01-01")) # using the default vcov test.vanilla <- sr_equality_test(some.rets,type="F") print(test.vanilla) if (require(sandwich)) { # and a fancy one: test.HAC <- sr_equality_test(some.rets,type="F",vcov.func=vcovHAC) print(test.HAC) } @ %UNFOLD \section{Asymptotics}%FOLDUP \subsection{Variance Covariance of the \txtSR}%FOLDUP Given \ssiz observations of the returns of \nlatf assets, the covariance of the sample \txtSR can be estimated via the Delta method. This operates by stacking the \nlatf vector of returns on top of the \nlatf vector of the returns squared element-wise. One then estimates the covariance of this $2\nlatf$ vector, using one's favorite covariance estimator. This estimate is then translated back to an estimate of the covariance of the \nlatf vector of \txtSR values. This process was described by Lo, and Ledoit and Wolf for the $\nlatf=1$ case, and is used in the Sharpe equality tests of Leung \etal and of Wright \etal. \cite{lo2002,Ledoit2008850,Leung2008,Wright2012} This basic function is available via the \Rcode{sr\_vcov} function, which allows one to pass in a function to perform the covariance estimate of the $2\nlatf$ vector. The passed in function should operate on an \Rcode{lm} object. The default is the \Rcode{vcov} function; other sensible choices include \Rcode{sandwich:vcovHC} and \Rcode{sandwich:vcovHAC} to deal with heteroskedasticity and autocorrelation. <<'sr_vcov'>>= # get returns some.rets <- as.matrix(get.rets(c("IBM","AAPL","XOM"), from="2007-01-01",to="2013-01-01")) ope <- 252 vanilla.Sig <- sr_vcov(some.rets,ope=ope) print(vanilla.Sig) if (require(sandwich)) { HC.Sig <- sr_vcov(some.rets,vcov=vcovHC,ope=ope) print(HC.Sig$Ohat) } if (require(sandwich)) { HAC.Sig <- sr_vcov(some.rets,vcov=vcovHAC,ope=ope) print(HAC.Sig$Ohat) } @ %UNFOLD \subsection{Variance Covariance of the Markowitz Portfolio}%FOLDUP Given a vector of returns, \vreti, prepending a one and taking the uncentered moment gives a `unified' parameter: $$ \pvsm \defeq \E{\ogram{\avreti}} = \twobytwo{1}{\tr{\pvmu}}{\pvmu}{\pvsig + \ogram{\pvmu}}, $$ where $\avreti = \asvec{1,\tr{\vreti}}$. The inverse of \pvsm contains the (negative) Markowitz portfolio: $$ \minv{\pvsm} = \twobytwo{1 + \qiform{\pvsig}{\pvmu}}{-\tr{\pvmu}\minv{\pvsig}}{-\minv{\pvsig}\pvmu}{\minv{\pvsig}} = \twobytwo{1 + \psnrsqopt}{-\tr{\pportwopt}}{-\pportwopt}{\minv{\pvsig}}. $$ These computations hold for the sample estimator $\svsm \defeq \oneby{\ssiz}\gram{\amreti}$, where \amreti is the matrix whose rows are the \ssiz vectors \tr{\avreti[i]}. By the multivariate Central Limit Theorem, \cite{wasserman2004all} $$ \sqrt{\ssiz}\wrapParens{\fvech{{\svsm}} - \fvech{{\pvsm}}} \rightsquigarrow \normlaw{0,\pvvar}, $$ for some \pvvar which can be estimated from the data. Via the delta method the asymptotic distribution of \fvech{\minv{\svsm}}, which contains $-\minv{\svsig}\svmu$, can be found. See the other vignette for more details. This functionality is available via the \Rcode{ism\_vcov} function. This function returns the sample estimate and estimated asymptotic variance of \fvech{\minv{\svsm}} (with the leading one removed). This function also allows one to pass in a callback to perform the covariance estimation. The default is the \Rcode{vcov} function; other sensible choices include \Rcode{sandwich:vcovHC} and \Rcode{sandwich:vcovHAC}. <<'marko_vcov'>>= # get returns some.rets <- get.rets(c("IBM","AAPL","XOM"), from="2007-01-01",to="2013-01-01") ism.wald <- function(X,vcov.func=vcov) { # negating returns is idiomatic to get + Markowitz ism <- ism_vcov(- as.matrix(X),vcov.func=vcov.func) ism.mu <- ism$mu[1:ism$p] ism.Sg <- ism$Ohat[1:ism$p,1:ism$p] retval <- ism.mu / sqrt(diag(ism.Sg)) dim(retval) <- c(ism$p,1) rownames(retval) <- rownames(ism$mu)[1:ism$p] return(retval) } wald.stats <- ism.wald(some.rets) print(t(wald.stats)) if (require(sandwich)) { wald.stats <- ism.wald(some.rets,vcov.func=sandwich::vcovHAC) print(t(wald.stats)) } @ %UNFOLD %UNFOLD \section{Miscellanea}%FOLDUP \subsection{Distribution Functions}%FOLDUP There are \Rcode{dpqr} functions for the \txtSR distribution, as well as the `optimal \txtSR' distribution. These are merely rescaled non-central \tstat and \Fstat distributions, provided for convenience, and for testing purposes. See the help for \Rfunction{dsr} and \Rfunction{dsropt} for more details. %UNFOLD %UNFOLD % bibliography%FOLDUP \nocite{CambridgeJournals:4493808,lo2002,Lecoutre2007,Johnson:1940,Ledoit2008850,sref2011} \nocite{1307.0450} \nocite{pav_ssc} %\bibliographystyle{jss} %\bibliographystyle{ieeetr} \bibliographystyle{plainnat} %\bibliographystyle{acm} \bibliography{SharpeR} %UNFOLD \end{document} %for vim modeline: (do not edit) % vim:fdm=marker:fmr=FOLDUP,UNFOLD:cms=%%s:syn=rnoweb:ft=rnoweb:cole=0