# 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
#' @importFrom methods setOldClass
#' @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:
#' if (require(xts)) {
#' data(stock_returns)
#' IBM <- stock_returns[,'IBM']
#' asr <- as.sr(IBM,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)
#' }
#'
#' # using real data.
#' if (require(xts)) {
#' data(stock_returns)
#' asro <- as.sropt(stock_returns)
#' }
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)
#'
#' # using real data.
#' if (require(xts)) {
#' data(stock_returns)
#' # hedge out SPY
#' G <- diag(dim(stock_returns)[2])[3,]
#' asro <- as.del_sropt(stock_returns,G=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.
#' @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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.