# R/sr.r In shabbychef/SharpeR: Statistical Significance of the Sharpe Ratio

# Copyright 2012-2014 Steven E. Pav. All Rights Reserved.
# Author: Steven E. Pav

# This file is part of SharpeR.
#
# SharpeR is free software: you can redistribute it and/or modify
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# SharpeR is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with SharpeR.  If not, see <http://www.gnu.org/licenses/>.

# env var:
# nb:
# todo:
# changelog:
#
# Created: 2013.04.16
# Copyright: Steven E. Pav, 2012-2013
# Author: Steven E. Pav

#' @include utils.r
#' @include distributions.r
#' @include estimation.r

setOldClass(c('sr','sropt','del_sropt','summary.sr','summary.sropt'))

########################################################################
# Sharpe Ratio#FOLDUP

# spawn a "SR" object.
#' @title Create an 'sr' object.
#'
#' @description
#'
#' Spawns an object of class \code{sr}.
#'
#' @details
#'
#' The \code{sr} class contains information about a rescaled t-statistic.
#' The following are list attributes of the object:
#' \describe{
#' \item{sr}{The Sharpe ratio statistic.}
#' \item{df}{The d.f. of the equivalent t-statistic.}
#' \item{c0}{The drag 'risk free rate' used.}
#' \item{ope}{The 'observations per epoch'.}
#' \item{rescal}{The rescaling parameter.}
#' \item{epoch}{The string name of the 'epoch'.}
#' }
#'
#' The stored Sharpe statistic, \code{sr} is equal to the t-statistic
#' times \eqn{rescal * sqrt{ope}}{rescal * sqrt(ope)}.
#'
#' For the most part, this constructor should \emph{not} be called directly,
#' Sharpe ratio.
#'
#' @param sr a Sharpe ratio statistic.
#' @param df the degrees of freedom of the equivalent t-statistic.
#' @param c0 the 'risk-free' or 'disastrous' rate of return. this is
#'        assumed to be given in the same units as x, \emph{not}
#'        in 'annualized' terms.
#' @template param-ope
#' @param rescal the rescaling parameter.
#' @param epoch the string representation of the 'epoch', defaulting
#'        to 'yr'.
#' @template param-cumulants
#' @keywords univar
#' @return a list cast to class \code{sr}.
#' @rdname sr
#' @export sr
#' @template etc
#' @template sr
#'
#' @note
#' 2FIX: allow rownames?
#'
#' @examples
#' ope <- 253
#' zeta <- 1.0
#' n <- 3 * ope
#' rvs <- rsr(1,n,zeta,ope=ope)
#' roll.own <- sr(sr=rvs,df=n-1,ope=ope,rescal=sqrt(1/n))
#' # put a bunch in. naming becomes a problem.
#' rvs <- rsr(5,n,zeta,ope=ope)
#' roll.own <- sr(sr=rvs,df=n-1,ope=ope,rescal=sqrt(1/n))
#'
sr <- function(sr,df,c0=0,ope=1,rescal=sqrt(1/(df+1)),epoch="yr",cumulants=NULL) {
retval <- list(sr = sr,df = df,c0 = c0,
ope = ope,rescal = rescal,epoch = epoch,
cumulants=cumulants)
class(retval) <- "sr"
return(retval)
}
# compute SR in only one place. I hope.
.compute_sr <- function(mu,c0,sigma,ope) {
sr <- (mu - c0) / sigma
if (!missing(ope))
sr <- sr * sqrt(ope)
return(sr)
}
.get_strat_names <- function(xstat) {
retval <- if (is.null(rownames(xstat))) rep(c("Sharpe"),length(xstat)) else rownames(xstat)
return(retval)
}
#' @title Compute the Sharpe ratio.
#'
#' @description
#'
#' Computes the Sharpe ratio of some observed returns.
#'
#' @details
#'
#' Suppose \eqn{x_i}{xi} are \eqn{n}{n} independent returns of some
#' asset.
#' Let \eqn{\bar{x}}{xbar} be the sample mean, and \eqn{s}{s} be
#' the sample standard deviation (using Bessel's correction). Let \eqn{c_0}{c0}
#' be the 'risk free rate'.  Then
#' \deqn{z = \frac{\bar{x} - c_0}{s}}{z = (xbar - c0)/s}
#' is the (sample) Sharpe ratio.
#'
#' The units of \eqn{z}{z} are \eqn{\mbox{time}^{-1/2}}{per root time}.
#' Typically the Sharpe ratio is \emph{annualized} by multiplying by
#' \eqn{\sqrt{\mbox{ope}}}{sqrt(ope)}, where \eqn{\mbox{ope}}{ope}
#' is the number of observations
#' per year (or whatever the target annualization epoch.)
#'
#' Note that if \code{ope} is not given, the converter from \code{xts}
#' attempts to infer the observations per year, without regard to
#' the name of the \code{epoch} given.
#'
#' @param x vector of returns, or object of class \code{data.frame}, \code{xts},
#'        or \code{lm}.
#' @param c0 the 'risk-free' or 'disastrous' rate of return. this is
#'        assumed to be given in the same units as x, \emph{not}
#'        in 'annualized' terms.
#' @template param-ope
#' @param na.rm logical.  Should missing values be removed?
#' @param epoch the string representation of the 'epoch', defaulting
#'        to 'yr'.
#' @param higher_order  a Boolean. If true, we compute
#'        cumulants of the returns to leverage higher order accuracy formulae
#'        when possible.
#' @keywords univar
#' @return a list containing the following components:
#' \describe{
#' \item{sr}{the annualized Sharpe ratio.}
#' \item{df}{the t-stat degrees of freedom.}
#' \item{c0}{the risk free term.}
#' \item{ope}{the annualization factor.}
#' \item{rescal}{the rescaling factor.}
#' \item{epoch}{the string epoch.}
#' }
#' cast to class \code{sr}.
#' @rdname as.sr
#' @export as.sr
#' @template etc
#' @template sr
#' @template ref-Lo
#'
#' @examples
#' # Sharpe's 'model': just given a bunch of returns.
#' asr <- as.sr(rnorm(253*3),ope=253)
#' # or a matrix, with a name
#' my.returns <- matrix(rnorm(253*3),ncol=1)
#' colnames(my.returns) <- c("my strategy")
#' asr <- as.sr(my.returns)
#' # given an xts object:
#' \dontrun{
#' if (require(quantmod)) {
#'   IBM <- getSymbols('IBM',auto.assign=FALSE)
#'   asr <- as.sr(lrets,na.rm=TRUE)
#' }
#' }
#' # on a linear model, find the 'Sharpe' of the residual term
#' nfac <- 5
#' nyr <- 10
#' ope <- 253
#' set.seed(as.integer(charToRaw("determinstic")))
#' Factors <- matrix(rnorm(ope*nyr*nfac,mean=0,sd=0.0125),ncol=nfac)
#' Betas <- exp(0.1 * rnorm(dim(Factors)[2]))
#' Returns <- (Factors %*% Betas) + rnorm(dim(Factors)[1],mean=0.0005,sd=0.012)
#' APT_mod <- lm(Returns ~ Factors)
#' asr <- as.sr(APT_mod,ope=ope)
#' # try again, but make the Returns independent of the Factors.
#' Returns <- rnorm(dim(Factors)[1],mean=0.0005,sd=0.012)
#' APT_mod <- lm(Returns ~ Factors)
#' asr <- as.sr(APT_mod,ope=ope)
#'
#' # compute the Sharpe of a bunch of strategies:
#' my.returns <- matrix(rnorm(253*3*4),ncol=4)
#' asr <- as.sr(my.returns)  # without sensible colnames?
#' colnames(my.returns) <- c("strat a","strat b","strat c","strat d")
#' asr <- as.sr(my.returns)
#'
as.sr <- function(x,c0=0,ope=1,na.rm=FALSE,epoch="yr",higher_order=FALSE) {
UseMethod("as.sr", x)
}
.as.sr.unified <- function(x,mu,sigma,c0,ope,na.rm,epoch,cumulants=NULL) {
z <- .compute_sr(mu,c0,sigma,ope)
dim(z) <- c(length(mu),1)

# what bother; try to find the rownames
if (length(dimnames(x)) == 1)
rnz <- unlist(dimnames(x))
else
rnz <- unlist(dimnames(x)[2])
if (is.null(rnz))
rnz <- deparse(substitute(x))

if (length(rnz) == dim(z)[1]) {
rownames(z) <- rnz
} else {
rstub <- deparse(substitute(x))
rownames(z) <- sapply(1:dim(z)[1],
function(n) { sprintf('%s_%d',rstub,n) })
}

# make it a col vector if not otherwise defined.
if (is.null(dim(x)))
dim(x) <- c(length(x),1)
if (na.rm)
df <- apply(!is.na(x),2,sum)
else
df <- dim(x)[1]

# n.b. sr$df is n-1 retval <- sr(z,df=df-1,c0=c0,ope=ope, rescal=1/sqrt(df),epoch=epoch,cumulants=cumulants) return(retval) } #' @rdname as.sr #' @method as.sr default #' @export as.sr.default <- function(x,c0=0,ope=1,na.rm=FALSE,epoch="yr",higher_order=FALSE) { if (na.rm) { x <- na.omit(x) } mu <- mean(x) sigma <- sd(x) if (higher_order) { # find this in sr_bias.r cumulants <- matrix(.smplgamma(x,muy=mu),ncol=1) # make it a matrix. rownames(cumulants) <- paste0('gamma_',1:4) } else { cumulants <- NULL } retval <- .as.sr.unified(x=x,mu=mu,sigma=sigma,c0=c0,ope=ope, na.rm=na.rm,epoch=epoch,cumulants=cumulants) return(retval) } #' @rdname as.sr #' @method as.sr matrix #' @export as.sr.matrix <- function(x,c0=0,ope=1,na.rm=FALSE,epoch="yr",higher_order=FALSE) { mu <- apply(x,2,mean,na.rm=na.rm) sigma <- apply(x,2,sd,na.rm=na.rm) # remember that the cumulants in general will be a matrix... if (higher_order) { # find this in sr_bias.r cumulants <- apply(x,2,.smplgamma,na.rm=na.rm) rownames(cumulants) <- paste0('gamma_',1:4) } else { cumulants <- NULL } retval <- .as.sr.unified(x=x,mu=mu,sigma=sigma,c0=c0,ope=ope, na.rm=na.rm,epoch=epoch,cumulants=cumulants) return(retval) } #' @rdname as.sr #' @method as.sr data.frame #' @export as.sr.data.frame <- function(x,c0=0,ope=1,na.rm=FALSE,epoch="yr",higher_order=FALSE) { retval <- as.sr.matrix(x,c0=c0,ope=ope,na.rm=na.rm,epoch=epoch,higher_order=higher_order) return(retval) } #' @rdname as.sr #' @method as.sr lm #' @export as.sr.lm <- function(x,c0=0,ope=1,na.rm=FALSE,epoch="yr",higher_order=FALSE) { mu <- x$coefficients["(Intercept)"]
sigma <- sqrt(deviance(x) / x$df.residual) if (higher_order) { warning('cannot compute higher order cumulants on lm object; ignoring the flag') } z <- .compute_sr(mu,c0,sigma,ope) dim(z) <- c(1,1) rownames(z) <- deparse(substitute(x)) XXinv <- vcov(x) / sigma^2 rescal <- sqrt(XXinv["(Intercept)","(Intercept)"]) retval <- sr(z,df=x$df.residual,c0=c0,ope=ope,
rescal=rescal,epoch=epoch,cumulants=NULL)
return(retval)
}
#' @rdname as.sr
#' @method as.sr xts
#' @export
as.sr.xts <- function(x,c0=0,ope=1,na.rm=FALSE,epoch="yr",higher_order=FALSE) {
if (missing(ope) && missing(epoch)) {
ope <- .infer_ope_xts(x)
epoch <- "yr"
}
retval <- as.sr.data.frame(as.data.frame(x),c0=c0,ope=ope,na.rm=na.rm,epoch=epoch,higher_order=higher_order)
return(retval)
}
#' @rdname as.sr
#' @method as.sr timeSeries
#' @export
as.sr.timeSeries <- function(x,c0=0,ope=1,na.rm=FALSE,epoch="yr",higher_order=FALSE) {
# you want to do this, but requires xts package. oops.
#retval <- as.sr.xts(as.xts(x),c0=c0,ope=ope,na.rm=na.rm,epoch=epoch)
if (missing(ope) && missing(epoch)) {
ope <- .infer_ope_xts(x)
epoch <- "yr"
}
retval <- as.sr.data.frame(as.data.frame(x),c0=c0,ope=ope,na.rm=na.rm,epoch=epoch,higher_order=higher_order)
return(retval)
}
#' @title Is this in the "sr" class?
#'
#' @description
#'
#' Checks if an object is in the class \code{'sr'}
#'
#' @details
#'
#' To satisfy the minimum requirements of an S3 class.
#'
#' @param x an object of some kind.
#' @return a boolean.
#' @seealso sr
#' @template etc
#' @family sr
#' @export
#'
#' @examples
#' rvs <- as.sr(rnorm(253*8),ope=253)
#' is.sr(rvs)
is.sr <- function(x) inherits(x,"sr")

# ' @export
# ' @export
format.sr <- function(x,...) {
# oh! ugly! ugly!
retval <- capture.output(print.sr(x,...))
return(retval)
}
#' @title Print values.
#'
#' @description
#'
#' Displays an object, returning it \emph{invisibly},
#' (via \code{invisible(x)}.)
#'
#' @param x an object of class \code{sr} or \code{sropt}.
#' @template param-ellipsis
#'
#' @return the object, wrapped in \code{invisible}.
#' @rdname print
#' @method print sr
#' @export
#' @export
#' @template etc
#' @template sr
#' @examples
#' # compute a 'daily' Sharpe
#' mysr <- as.sr(rnorm(253*8),ope=1,epoch="day")
#' print(mysr)
#' ope <- 253
#' zeta <- 1.0
#' n <- 6 * ope
#' rvs <- rsr(1,n,zeta,ope=ope)
#' roll.own <- sr(sr=rvs,df=n-1,ope=ope,rescal=sqrt(1/n))
#' print(roll.own)
#' # put a bunch in. naming becomes a problem.
#' rvs <- rsr(5,n,zeta,ope=ope)
#' roll.own <- sr(sr=rvs,df=n-1,ope=ope,rescal=sqrt(1/n))
#' print(roll.own)
#' # for sropt objects:
#' nfac <- 5
#' nyr <- 10
#' ope <- 253
#' # simulations with no covariance structure.
#' # under the null:
#' set.seed(as.integer(charToRaw("be determinstic")))
#' Returns <- matrix(rnorm(ope*nyr*nfac,mean=0,sd=0.0125),ncol=nfac)
#' asro <- as.sropt(Returns,drag=0,ope=ope)
#' print(asro)
print.sr <- function(x,...) {
tval <- .sr2t(x)
pval <- pt(tval,x$df,lower.tail=FALSE) serr <- se(x,type="t") coefs <- cbind(x$sr,serr,tval,pval)
#colnames(coefs) <- c("stat","t.stat","p.value")
colnames(coefs) <- c(paste(c("SR/sqrt(",x$epoch,")"),sep="",collapse=""), "Std. Error","t value","Pr(>t)") rownames(coefs) <- .get_strat_names(x$sr)
printCoefmat(coefs,P.values=TRUE,has.Pvalue=TRUE,
digits=max(2, getOption("digits") - 3),
cs.ind=c(1,2),tst.ind=c(3),dig.tst=2)
# invisible return
invisible(x)
}
# print.sr <- function(x,...) cat(format(x,...), "\n")

# SR methods#FOLDUP
# 2FIX: make the x$rescal * sqrt(x$ope) occur in one place only ...
# get the t-stat associated with an SR object.
.sr2t <- function(x) {
tval <- x$sr / (x$rescal * sqrt(x$ope)) return(tval) } # and the reverse .t2sr <- function(x,tval) { srval <- tval * (x$rescal * sqrt(x$ope)) return(srval) } .psr <- function(q,zeta,...) { retv <- prt(q$sr,df=q$df,K=(q$rescal * sqrt(q$ope)),rho=zeta,...) return(retv) } # no longer needed? #.dsr <- function(q,zeta,...) { #retv <- drt(q$sr,df=q$df,K=(q$rescal * sqrt(q$ope)),rho=zeta,...) #return(retv) #} #' @title Change the annualization of a Sharpe ratio. #' #' @description #' #' Changes the annualization factor of a Sharpe ratio statistic, or the rate at #' which observations are made. #' #' @param object an object of class \code{sr} or \code{sropt}. #' @param new.ope the new observations per epoch. If none given, it is #' not updated. #' @param new.epoch a string representation of the epoch. If none given, it is not #' updated. #' @return the input object with the annualization and/or epoch updated. #' @seealso sr #' @family sr #' @seealso sropt #' @family sropt #' @template etc #' @export #' @rdname reannualize #' #' @examples #' # compute a 'daily' Sharpe #' mysr <- as.sr(rnorm(253*8),ope=1,epoch="day") #' # turn into annual #' mysr2 <- reannualize(mysr,new.ope=253,new.epoch="yr") #' #' # for sropt #' ope <- 253 #' zeta.s <- 1.0 #' df1 <- 10 #' df2 <- 6 * ope #' rvs <- rsropt(1,df1,df2,zeta.s,ope,drag=0) #' roll.own <- sropt(z.s=rvs,df1,df2,drag=0,ope=ope,epoch="yr") #' # make 'monthly' #' roll.monthly <- reannualize(roll.own,new.ope=21,new.epoch="mo.") #' # make 'daily' #' roll.daily <- reannualize(roll.own,new.ope=1,new.epoch="day") reannualize <- function(object,new.ope=NULL,new.epoch=NULL) { UseMethod("reannualize", object) } #' @rdname reannualize #' @method reannualize sr #' @export reannualize.sr <- function(object,new.ope=NULL,new.epoch=NULL) { if (!is.sr(object)) stop("must give sr object") if (!missing(new.ope)) { object$sr <- object$sr * sqrt(new.ope / object$ope)
object$ope <- new.ope } if (!missing(new.epoch)) object$epoch <- new.epoch
return(object)
}
#' @rdname reannualize
#' @method reannualize sropt
#' @export
reannualize.sropt <- function(object,new.ope=NULL,new.epoch=NULL) {
if (!is.sropt(object)) stop("must give sropt object")
if (!missing(new.ope)) {
object$sropt <- object$sropt * sqrt(new.ope / object$ope) object$ope <- new.ope
}
if (!missing(new.epoch)) object$epoch <- new.epoch return(object) } #UNFOLD #UNFOLD ######################################################################## # Markowitz #FOLDUP # markowitz, possibly on a subset, or basket, of the # assets. .sub.markowitz <- function(mu,Sigma,w=NULL,H=NULL) { if (!is.null(H)) { mu <- H %*% mu Sigma <- .qoform(H,Sigma) } if (is.null(w)) w <- solve(Sigma,mu) zeta.sq <- t(mu) %*% w if (!is.null(H)) w <- t(H) %*% mu retval <- list(w=w,zeta.sq=zeta.sq) return(retval) } markowitz <- function(mu,Sigma,df2,w=NULL,H=NULL) { # also computes Hotelling's T2 retval <- .sub.markowitz(mu,Sigma,w,H) df1 <- ifelse(is.null(H),length(mu),dim(H)[1]) retval <- c(retval,list(mu=mu,Sigma=Sigma, df1=df1,df2=df2, T2=df2*retval$zeta.sq))
class(retval) <- c("markowitz")
return(retval)
}

as.markowitz <- function(X,...) {
UseMethod("as.markowitz", X)
}

# compute the (constrained) markowitz portfolio
as.markowitz.default <- function(X,mu=NULL,Sigma=NULL,...) {
X <- na.omit(X)
if (is.null(mu))
mu <- colMeans(X)
if (is.null(Sigma))
Sigma <- cov(X)
df2 <- dim(X)[1]
retval <- markowitz(mu,Sigma,df2,...)
return(retval)
}

#UNFOLD

########################################################################
# Optimal Sharpe ratio#

# spawn a "SROPT" object.
#' @title Create an 'sropt' object.
#'
#' @description
#'
#' Spawns an object of class \code{sropt}.
#'
#' @details
#'
#' The \code{sropt} class contains information about a rescaled T^2-statistic.
#' The following are list attributes of the object:
#' \describe{
#' \item{sropt}{The (optimal) Sharpe ratio statistic.}
#' \item{df1}{The number of assets.}
#' \item{df2}{The number of observations.}
#' \item{drag}{The drag term, which is the 'risk free rate' divided by
#' the maximum risk.}
#' \item{ope}{The 'observations per epoch'.}
#' \item{epoch}{The string name of the 'epoch'.}
#' }
#'
#' For the most part, this constructor should \emph{not} be called directly,
#' needed statistics.
#'
#' @param z.s an optimum Sharpe ratio statistic.
#' @inheritParams dsropt
#' @param drag the 'drag' term, \eqn{c_0/R}{c0/R}. defaults to 0. It is assumed
#'        that \code{drag} has been annualized, \emph{i.e.} has been multiplied
#'        by \eqn{\sqrt{ope}}{sqrt(ope)}. This is in contrast to the \code{c0}
#' @param T2 the Hotelling \eqn{T^2}{T^2} statistic. If not given, it is
#'        computed from the given information.
#' @template param-ope
#' @param epoch the string representation of the 'epoch', defaulting
#'        to 'yr'.
#' @keywords univar
#' @return a list cast to class \code{sropt}, with the following attributes:
#' \describe{
#' \item{sropt}{the optimal Sharpe statistic.}
#' \item{df1}{the number of assets.}
#' \item{df2}{the number of observed vectors.}
#' \item{drag}{the input \code{drag} term.}
#' \item{ope}{the input \code{ope} term.}
#' \item{epoch}{the input \code{epoch} term.}
#' \item{T2}{the Hotelling \eqn{T^2} statistic.}
#' }
#' @rdname sropt
#' @export
#' @template etc
#' @template sropt
#'
#' @note
#' 2FIX: allow rownames?
#'
#' @examples
#' ope <- 253
#' zeta.s <- 1.0
#' df1 <- 10
#' df2 <- 6 * ope
#' set.seed(as.integer(charToRaw("fix seed")))
#' rvs <- rsropt(1,df1,df2,zeta.s,ope,drag=0)
#' roll.own <- sropt(z.s=rvs,df1,df2,drag=0,ope=ope)
#' print(roll.own)
#' # put a bunch in. naming becomes a problem.
#' rvs <- rsropt(5,df1,df2,zeta.s,ope,drag=0)
#' roll.own <- sropt(z.s=rvs,df1,df2,drag=0,ope=ope)
#' print(roll.own)
sropt <- function(z.s,df1,df2,drag=0,ope=1,epoch="yr",T2=NULL) {
retval <- list(sropt = z.s,df1 = df1,df2 = df2,
drag = drag,ope = ope,epoch = epoch)
if (is.null(T2))
T2 <- .sropt2T(retval)
retval$T2 <- T2 # not sure I should do this ... retval$pval <- pT2(T2,retval$df1,retval$df2,lower.tail=FALSE)
class(retval) <- "sropt"
return(retval)
}
#' @title Compute the Sharpe ratio of the Markowitz portfolio.
#'
#' @description
#'
#' Computes the Sharpe ratio of the Markowitz portfolio of some observed returns.
#'
#' @details
#'
#' Suppose \eqn{x_i}{xi} are \eqn{n}{n} independent draws of a \eqn{q}{q}-variate
#' normal random variable with mean \eqn{\mu}{mu} and covariance matrix
#' \eqn{\Sigma}{Sigma}. Let \eqn{\bar{x}}{xbar} be the (vector) sample mean, and
#' \eqn{S}{S} be the sample covariance matrix (using Bessel's correction). Let
#' \deqn{\zeta(w) = \frac{w^{\top}\bar{x} - c_0}{\sqrt{w^{\top}S w}}}{zeta(w) = (w'xbar - c0)/sqrt(w'Sw)}
#' be the (sample) Sharpe ratio of the portfolio \eqn{w}{w}, subject to
#' risk free rate \eqn{c_0}{c0}.
#'
#' Let \eqn{w_*}{w*} be the solution to the portfolio optimization problem:
#' \deqn{\max_{w: 0 < w^{\top}S w \le R^2} \zeta(w),}{max {zeta(w) | 0 < w'Sw <= R^2},}
#' with maximum value \eqn{z_* = \zeta\left(w_*\right)}{z* = zeta(w*)}.
#' Then
#' \deqn{w_* = R \frac{S^{-1}\bar{x}}{\sqrt{\bar{x}^{\top}S^{-1}\bar{x}}}}{%
#' w* = R S^-1 xbar / sqrt(xbar' S^-1 xbar)}
#' and
#' \deqn{z_* = \sqrt{\bar{x}^{\top} S^{-1} \bar{x}} - \frac{c_0}{R}}{%
#' z* = sqrt(xbar' S^-1 xbar) - c0/R}
#'
#' The units of \eqn{z_*}{z*} are \eqn{\mbox{time}^{-1/2}}{per root time}.
#' Typically the Sharpe ratio is \emph{annualized} by multiplying by
#' \eqn{\sqrt{\mbox{ope}}}{sqrt(ope)}, where \eqn{\mbox{ope}}{ope}
#' is the number of observations
#' per year (or whatever the target annualization epoch.)
#'
#' Note that if \code{ope} and \code{epoch} are not given, the
#' converter from \code{xts} attempts to infer the observations per epoch,
#' assuming yearly epoch.
#'
#' @param X matrix of returns, or \code{xts} object.
#' @inheritParams sropt
#' @keywords univar
#' @return An object of class \code{sropt}.
#' @rdname as.sropt
#' @export
#' @template etc
#' @template sropt
#' @examples
#' nfac <- 5
#' nyr <- 10
#' ope <- 253
#' # simulations with no covariance structure.
#' # under the null:
#' set.seed(as.integer(charToRaw("be determinstic")))
#' Returns <- matrix(rnorm(ope*nyr*nfac,mean=0,sd=0.0125),ncol=nfac)
#' asro <- as.sropt(Returns,drag=0,ope=ope)
#' # under the alternative:
#' Returns <- matrix(rnorm(ope*nyr*nfac,mean=0.0005,sd=0.0125),ncol=nfac)
#' asro <- as.sropt(Returns,drag=0,ope=ope)
#' # generating correlated multivariate normal data in a more sane way
#' if (require(MASS)) {
#'   nstok <- 10
#'   nfac <- 3
#'   nyr <- 10
#'   ope <- 253
#'   X.like <- 0.01 * matrix(rnorm(500*nfac),ncol=nfac) %*%
#'     matrix(runif(nfac*nstok),ncol=nstok)
#'   Sigma <- cov(X.like) + diag(0.003,nstok)
#'   # under the null:
#'   Returns <- mvrnorm(ceiling(ope*nyr),mu=matrix(0,ncol=nstok),Sigma=Sigma)
#'   asro <- as.sropt(Returns,ope=ope)
#'   # under the alternative
#'   Returns <- mvrnorm(ceiling(ope*nyr),mu=matrix(0.001,ncol=nstok),Sigma=Sigma)
#'   asro <- as.sropt(Returns,ope=ope)
#' }
#' \dontrun{
#' # using real data.
#' if (require(quantmod)) {
#'   get.ret <- function(sym,...) {
#'     OHLCV <- getSymbols(sym,auto.assign=FALSE,...)
#'     # chomp first NA!
#'     lrets[-1,]
#'   }
#'   get.rets <- function(syms,...) {
#'		some.rets <- do.call("cbind",lapply(syms,get.ret,...))
#'	 }
#'   some.rets <- get.rets(c("IBM","AAPL","A","C","SPY","XOM"))
#'   asro <- as.sropt(some.rets)
#' }
#' }
as.sropt <- function(X,drag=0,ope=1,epoch="yr") {
UseMethod("as.sropt", X)
}
#' @rdname as.sropt
#' @method as.sropt default
#' @export
as.sropt.default <- function(X,drag=0,ope=1,epoch="yr") {
# somehow call sropt!
hotval <- as.markowitz(X)
hotval <- c(hotval,list(ope=ope,drag=drag))
# get the zeta.s
z.s <- .T2sropt(hotval)
dim(z.s) <- c(1,1)

# this stinks
retv <- sropt(z.s=z.s,df1=hotval$df1,df2=hotval$df2,
drag=drag,ope=ope,epoch=epoch,T2=hotval$T2) # 2FIX: merge markowitz in? return(retv) } #' @rdname as.sropt #' @method as.sropt xts #' @export as.sropt.xts <- function(X,drag=0,ope=1,epoch="yr") { if (missing(ope)) { ope <- .infer_ope_xts(X) } retval <- as.sropt.default(X,drag=drag,ope=ope,epoch=epoch) return(retval) } #' @title Is this in the "sropt" class? #' #' @description #' #' Checks if an object is in the class \code{'sropt'} #' #' @details #' #' To satisfy the minimum requirements of an S3 class. #' #' @param x an object of some kind. #' @return a boolean. #' @seealso sropt #' @template etc #' @family sropt #' @export #' #' @examples #' nfac <- 5 #' nyr <- 10 #' ope <- 253 #' # simulations with no covariance structure. #' # under the null: #' set.seed(as.integer(charToRaw("be determinstic"))) #' Returns <- matrix(rnorm(ope*nyr*nfac,mean=0,sd=0.0125),ncol=nfac) #' asro <- as.sropt(Returns,drag=0,ope=ope) #' is.sropt(asro) is.sropt <- function(x) inherits(x,"sropt") #' @rdname print #' @method print sropt #' @export print.sropt <- function(x,...) { Tval <- x$T2
pval <- .sropt.pval(x)
ci <- confint(x,level=0.95)
coefs <- cbind(x$sropt,sric(x),ci,Tval,pval) colnames(coefs) <- c(paste(c("SR/sqrt(",x$epoch,")"),sep="",collapse=""),
paste(c("SRIC/sqrt(",x$epoch,")"),sep="",collapse=""), colnames(ci)[1],colnames(ci)[2], "T^2 value","Pr(>T^2)") rownames(coefs) <- .get_strat_names(x$sropt)
printCoefmat(coefs,P.values=TRUE,has.Pvalue=TRUE,
digits=max(2, getOption("digits") - 3),
cs.ind=c(1,2,3,4),tst.ind=c(5),dig.tst=2)
}

# SROPT methods#
# get the T2-stat associated with an SROPT object.
.sropt2T <- function(x,Tval=x$T2) { if (is.null(Tval)) { tval <- .deannualize(x$sropt + x$drag,x$ope)
Tval <- x$df2 * (tval^2) } return(Tval) } # and the reverse; with a sticky zero. ouch. .T2sropt <- function(x,Tval=x$T2) {
z.star <- sqrt(pmax(Tval,0) / x$df2) z.star <- .annualize(z.star,x$ope)
z.star <- z.star - x$drag return(z.star) } .sropt.pval <- function(x,...) { Tval <- .sropt2T(x,...) retv <- x$pval
if (is.null(retv))
retv <- pT2(Tval,x$df1,x$df2,lower.tail=FALSE)
return(retv)
}
#
#

########################################################################
# Delta Squared Optimal Sharpe ratio#FOLDUP

# spawn a "DEL_SROPT" object.
#' @title Create an 'del_sropt' object.
#'
#' @description
#'
#' Spawns an object of class \code{del_sropt}.
#'
#' @details
#'
#' The \code{del_sropt} class contains information about the difference
#' between two rescaled T^2-statistics, useful for spanning
#' tests, and inference on hedged portfolios.
#' The following are list attributes of the object:
#' \describe{
#' \item{sropt}{The (optimal) Sharpe ratio statistic of
#' the 'full' set of assets.}
#' \item{sropt_sub}{The (optimal) Sharpe ratio statistic on
#' some subset, or linear subspace, of the assets.}
#' \item{df1}{The number of assets.}
#' \item{df2}{The number of observations.}
#' \item{df1.sub}{The number of degrees of freedom in the
#' hedge constraint.}
#' \item{drag}{The drag term, which is the 'risk free rate' divided by
#' the maximum risk.}
#' \item{ope}{The 'observations per epoch'.}
#' \item{epoch}{The string name of the 'epoch'.}
#' }
#'
#' For the most part, this constructor should \emph{not} be called directly,
#' needed statistics.
#'
#' @param z.s an optimum Sharpe ratio statistic, on some set of assets.
#' @param z.sub an optimum Sharpe ratio statistic, on a linear subspace
#' of the assets.  If larger than \code{z.s} an error is thrown.
#' @param df1.sub the rank of the linear subspace of the hedge
#' constraint.
#' by restricting attention to the subspace.
#' @inheritParams sropt
#' @template param-ope
#' @keywords univar
#' @return a list cast to class \code{del_sropt}, with attributes
#' \describe{
#' \item{sropt}{the optimal Sharpe statistic.}
#' \item{sropt.sub}{the optimal Sharpe statistic on the subspace.}
#' \item{df1}{the number of assets.}
#' \item{df2}{the number of observed vectors.}
#' \item{df1.sub}{the input \code{df1.sub} term.}
#' \item{drag}{the input \code{drag} term.}
#' \item{ope}{the input \code{ope} term.}
#' \item{T2}{the Hotelling \eqn{T^2} statistic.}
#' \item{T2.sub}{the Hotelling \eqn{T^2} statistic on the subspace.}
#' }
#'
#' @rdname del_sropt
#' @export
#' @template etc
#' @template del_sropt
#' @template warning
#'
#' @note
#' 2FIX: allow rownames?
#'
#' @examples
#' ope <- 253
#'
#' set.seed(as.integer(charToRaw("be determinstic")))
#' n.stock <- 10
#' X <- matrix(rnorm(1000*n.stock),nrow=1000)
#' Sigma <- cov(X)
#' mu <- colMeans(X)
#' w <- solve(Sigma,mu)
#' z <- t(mu) %*% w
#' n.sub <- 6
#' w.sub <- solve(Sigma[1:n.sub,1:n.sub],mu[1:n.sub])
#' z.sub <- t(mu[1:n.sub]) %*% w.sub
#' df1.sub <- n.stock - n.sub
#'
#' roll.own <- del_sropt(z.s=z,z.sub=z.sub,df1=10,df2=1000,
#'  df1.sub=df1.sub,ope=ope)
#' print(roll.own)
#'
del_sropt <- function(z.s,z.sub,df1,df2,df1.sub,drag=0,ope=1,epoch="yr") {
retval <- list(sropt = z.s,sropt.sub = z.sub,
sropt.del = sqrt((z.s)^2 - (z.sub)^2),
df1 = df1,df2 = df2,df1.sub = df1.sub,
drag = drag,ope = ope,epoch = epoch)

retval$T2.sub <- .sropt2T(list(sropt=z.sub,drag=drag,ope=ope,df2=df2)) #retval$T2 <- .sropt2T(retval)
retval$T2 <- .sropt2T(list(sropt=z.s,drag=drag,ope=ope,df2=df2)) retval$T2.del <- retval$T2 - retval$T2.sub

class(retval) <- "del_sropt"
return(retval)
}
# the statistical guts.
.del_sropt.asF <- function(x) {
# see Giri eqn (1.9) and section 3.
# make Z ~ B(bdf1,bdf2) under the null
nmin <- x$df2 - 1 Z.num <- (1 + x$T2.sub/nmin)
Z.denom <- (1 + x$T2/nmin) R1 <- 1 - (1 / Z.num) Z <- Z.num / Z.denom bdf1 <- (x$df2 - x$df1) / 2 bdf2 <- (x$df1 - x$df1.sub) / 2 # transform beta to F; # define the F so that it is non-central under the alternative # (on the top) df1 <- 2*bdf2 df2 <- 2*bdf1; Fval <- (df2 * (1-Z)) / (df1 * Z) pval <- pf(Fval,df1,df2,lower.tail=FALSE) # Giri's R1 retval <- list(Fval=Fval,df1=df1,df2=df2,N=x$df2,pval=pval,R1=R1)
return(retval)
}
#' @title Compute the Sharpe ratio of a hedged Markowitz portfolio.
#'
#' @description
#'
#' Computes the Sharpe ratio of the hedged Markowitz portfolio of some observed returns.
#'
#' @details
#'
#' Suppose \eqn{x_i}{xi} are \eqn{n}{n} independent draws of a \eqn{q}{q}-variate
#' normal random variable with mean \eqn{\mu}{mu} and covariance matrix
#' \eqn{\Sigma}{Sigma}. Let \eqn{G} be a \eqn{g \times q}{g x q} matrix
#' of rank \eqn{g}.
#' Let \eqn{\bar{x}}{xbar} be the (vector) sample mean, and
#' \eqn{S}{S} be the sample covariance matrix (using Bessel's correction).
#' Let
#' \deqn{\zeta(w) = \frac{w^{\top}\bar{x} - c_0}{\sqrt{w^{\top}S w}}}{zeta(w) = (w'xbar - c0)/sqrt(w'Sw)}
#' be the (sample) Sharpe ratio of the portfolio \eqn{w}{w}, subject to
#' risk free rate \eqn{c_0}{c0}.
#'
#' Let \eqn{w_*}{w*} be the solution to the portfolio optimization
#' problem:
#' \deqn{\max_{w: 0 < w^{\top}S w \le R^2,\,G S w = 0} \zeta(w),}{max {zeta(w) | 0 < w'Sw <= R^2, G S w = 0},}
#' with maximum value \eqn{z_* = \zeta\left(w_*\right)}{z* = zeta(w*)}.
#'
#' Note that if \code{ope} and \code{epoch} are not given, the
#' converter from \code{xts} attempts to infer the observations per epoch,
#' assuming yearly epoch.
#'
#' @param X matrix of returns, or \code{xts} object.
#' @param G an \eqn{g \times q}{g x q} matrix of hedge constraints. A
#' garden variety application would have \code{G} be one row of the
#' identity matrix, with a one in the column of the instrument to be
#' 'hedged out'.
#' @inheritParams sropt
#' @keywords univar
#' @return An object of class \code{del_sropt}.
#' @rdname as.del_sropt
#' @export
#' @template etc
#' @template del_sropt
#' @examples
#' nfac <- 5
#' nyr <- 10
#' ope <- 253
#' # simulations with no covariance structure.
#' # under the null:
#' set.seed(as.integer(charToRaw("be determinstic")))
#' Returns <- matrix(rnorm(ope*nyr*nfac,mean=0,sd=0.0125),ncol=nfac)
#' # hedge out the first one:
#' G <- matrix(diag(nfac)[1,],nrow=1)
#' asro <- as.del_sropt(Returns,G,drag=0,ope=ope)
#' print(asro)
#' G <- diag(nfac)[c(1:3),]
#' asro <- as.del_sropt(Returns,G,drag=0,ope=ope)
#' # compare to sropt on the remaining assets
#' # they should be close, but not exact.
#' asro.alt <- as.sropt(Returns[,4:nfac],drag=0,ope=ope)
#' \dontrun{
#' # using real data.
#' if (require(quantmod)) {
#'   get.ret <- function(sym,...) {
#'     OHLCV <- getSymbols(sym,auto.assign=FALSE,...)
#'     # chomp first NA!
#'     lrets[-1,]
#'   }
#'   get.rets <- function(syms,...) {
#'		some.rets <- do.call("cbind",lapply(syms,get.ret,...))
#'	 }
#'   some.rets <- get.rets(c("IBM","AAPL","A","C","SPY","XOM"))
#'   # hedge out SPY
#'   G <- diag(dim(some.rets)[2])[5,]
#'   asro <- as.del_sropt(some.rets,G)
#' }
#' }
as.del_sropt <- function(X,G,drag=0,ope=1,epoch="yr") {
UseMethod("as.del_sropt",X)
}
#' @rdname as.del_sropt
#' @method as.del_sropt default
#' @export
as.del_sropt.default <- function(X,G,drag=0,ope=1,epoch="yr") {
# somehow call sropt!
hotval <- as.markowitz(X)
hotval <- c(hotval,list(ope=ope,drag=drag))
# get the zeta.s
z.s <- .T2sropt(hotval)
dim(z.s) <- c(1,1)
if (!is.matrix(G)) G <- matrix(G,nrow=1)

if (length(hotval$w) != dim(G)[2]) stop("wrong size hedge constraint?") hotval.sub <- markowitz(hotval$mu,hotval$Sigma,hotval$df2,H=G)
hotval.sub <- c(hotval.sub,list(ope=ope,drag=drag))
# get the zeta.sub
z.sub <- .T2sropt(hotval.sub)
dim(z.sub) <- c(1,1)

# this stinks
retv <- del_sropt(z.s=z.s,z.sub=z.sub,df1=hotval$df1,df2=hotval$df2,
df1.sub=dim(G)[1],drag=drag,ope=ope,epoch=epoch)

return(retv)
}
#' @rdname as.del_sropt
#' @method as.del_sropt xts
#' @export
as.del_sropt.xts <- function(X,G,drag=0,ope=1,epoch="yr") {
if (missing(ope)) {
ope <- .infer_ope_xts(X)
}
retval <- as.del_sropt.default(X,G,drag=drag,ope=ope,epoch=epoch)
return(retval)
}
#' @title Is this in the "del_sropt" class?
#'
#' @description
#'
#' Checks if an object is in the class \code{'del_sropt'}
#'
#' @details
#'
#' To satisfy the minimum requirements of an S3 class.
#'
#' @param x an object of some kind.
#' @return a boolean.
#' @seealso del_sropt
#' @template etc
#' @family del_sropt
#' @export
#'
#' @examples
#' roll.own <- del_sropt(z.s=2,z.sub=1,df1=10,df2=1000,df1.sub=3,ope=1,epoch="yr")
#' is.sropt(roll.own)
is.del_sropt <- function(x) inherits(x,"del_sropt")

#' @rdname print
#' @method print del_sropt
#' @export
print.del_sropt <- function(x,...) {
Fandp <- .del_sropt.asF(x)
ci <- confint(x,level=0.95)
coefs <- cbind(x$sropt.del,ci,Fandp$Fval,Fandp$pval) colnames(coefs) <- c(paste(c("SR/sqrt(",x$epoch,")"),sep="",collapse=""),
colnames(ci)[1],colnames(ci)[2],
"F value","Pr(>F)")
rownames(coefs) <- .get_strat_names(x$sropt) printCoefmat(coefs,P.values=TRUE,has.Pvalue=TRUE, digits=max(2, getOption("digits") - 3), cs.ind=c(1),tst.ind=c(2),dig.tst=2) } #UNFOLD ## test them: #X <- matrix(rnorm(1000*10,0.1),ncol=10) #G <- diag(dim(X)[2])[1,] #foo <- as.del_sropt(X,G) #print(foo) #X <- matrix(rnorm(1000*10,0.00001),ncol=10) #foo <- as.del_sropt(X,G) #print(foo) # create a summary#FOLDUP #' @title Summarize a Sharpe, or (delta) optimal Sharpe object. #' #' @description #' #' Computes a \sQuote{summary} of an object, adding in some statistics. #' #' @details #' #' Enhances an object of class \code{sr}, \code{sropt} or \code{del_sropt} to also #' include t- or T-statistics, p-values, and so on. #' #' @param object an object of class \code{sr}, \code{sropt} or \code{del_sropt}. #' @param ... additional arguments affecting the summary produced, though #' ignored here. #' @inheritParams summary #' @return When an \code{sr} object is input, the object cast to class \code{summary.sr} with some #' additional fields: #' \describe{ #' \item{tval}{the equivalent t-statistic.} #' \item{pval}{the p-value under the null.} #' \item{serr}{the standard error of the Sharpe ratio.} #' } #' When an \code{sropt} object is input, the object cast to class \code{summary.sropt} with some #' additional fields: #' \describe{ #' \item{pval}{the p-value under the null.} #' \item{SRIC}{the SRIC value, see \code{\link{sric}}.} #' } #' #' @seealso \code{\link{print.sr}}. #' @rdname summary #' @template etc #' @template sr #' #' @examples #' # Sharpe's 'model': just given a bunch of returns. #' set.seed(1234) #' asr <- as.sr(rnorm(253*3),ope=253) #' summary(asr) #' @method summary sr #' @export summary.sr <- function(object,...) { object$tval <- .sr2t(object)
object$pval <- pt(object$tval,object$df,lower.tail=FALSE) object$serr <- se(object,type="t")
class(object) <- "summary.sr"
object
}
#' @rdname summary
#' @method summary sropt
#' @export
summary.sropt <- function(object,...) {
object$pval <- .sropt.pval(object) object$SRIC <- sric(object)
#... anything else? KRS?
class(object) <- "summary.sropt"
object
}

#UNFOLD

#for vim modeline: (do not edit)
# vim:fdm=marker:fmr=FOLDUP,UNFOLD:cms=#%s:syn=r:ft=r

shabbychef/SharpeR documentation built on May 29, 2019, 8:05 p.m.