%\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(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 <- 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 <- 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}

%\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
