R/sr.r

Defines functions summary.sropt summary.sr print.del_sropt is.del_sropt as.del_sropt.xts as.del_sropt.default as.del_sropt del_sropt.asF del_sropt sropt.pval T2sropt sropt2T print.sropt is.sropt as.sropt.xts as.sropt.default as.sropt sropt as.markowitz.default as.markowitz markowitz sub.markowitz reannualize.sropt reannualize.sr reannualize psr t2sr sr2t print.sr format.sr is.sr as.sr.timeSeries as.sr.xts as.sr.lm as.sr.data.frame as.sr.matrix as.sr.default as.sr.unified as.sr get_strat_names compute_sr sr

Documented in as.del_sropt as.del_sropt.default as.del_sropt.xts as.sr as.sr.data.frame as.sr.default as.sr.lm as.sr.matrix as.sropt as.sropt.default as.sropt.xts as.sr.timeSeries as.sr.xts del_sropt is.del_sropt is.sr is.sropt print.del_sropt print.sr print.sropt psr reannualize reannualize.sr reannualize.sropt sr sropt summary.sr summary.sropt

# 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
# it under the terms of the GNU Lesser General Public License as published by
# 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: 
# see also:
# todo:
# changelog: 
#
# Created: 2013.04.16
# Copyright: Steven E. Pav, 2012-2013
# Author: Steven E. Pav
# Comments: 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,
#' rather \code{\link{as.sr}} should be called instead to compute the
#' 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}.
#' @seealso \code{\link{as.sr}}
#' @rdname sr
#' @export sr
#' @template etc
#' @template sr
#'
#' @note
#' 2FIX: allow rownames? 
#'
#' @examples 
#' # roll your own.
#' 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}.
#' @seealso sr-distribution functions, \code{\link{dsr}, \link{psr}, \link{qsr}, \link{rsr}}
#' @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)
#'   lrets <- diff(log(IBM[,"IBM.Adjusted"]))
#'   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)
#' # roll your own.
#' 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)
}
# @hadley's suggested form
# 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,
#' rather \code{\link{as.sropt}} should be called instead to compute the
#' 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}
#'        term given to \code{\link{sr}}.
#' @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.}
#' }
#' @seealso \code{\link{as.sropt}}
#' @rdname sropt
#' @export 
#' @template etc
#' @template sropt
#'
#' @note
#' 2FIX: allow rownames?
#'
#' @examples 
#' # roll your own.
#' 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}.
#' @seealso \code{\link{sropt}}, \code{\link{sr}}, sropt-distribution functions, 
#' \code{\link{dsropt}, \link{psropt}, \link{qsropt}, \link{rsropt}}
#' @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,...)
#'     lrets <- diff(log(OHLCV[,paste(c(sym,"Adjusted"),collapse=".",sep="")]))
#'     # 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,
#' rather \code{\link{as.del_sropt}} should be called instead to compute the
#' 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.}
#' }
#'
#' @seealso \code{\link{as.del_sropt}}
#' @rdname del_sropt
#' @export 
#' @template etc
#' @template del_sropt
#' @template warning
#'
#' @note
#' 2FIX: allow rownames?
#'
#' @examples 
#' # roll your own.
#' 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}.
#' @seealso \code{\link{del_sropt}}, \code{\link{sropt}}, 
#' \code{\link{sr}}
#' @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,...)
#'     lrets <- diff(log(OHLCV[,paste(c(sym,"Adjusted"),collapse=".",sep="")]))
#'     # 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

Try the SharpeR package in your browser

Any scripts or data that you put into this service are public.

SharpeR documentation built on Oct. 8, 2018, 1:05 a.m.