R/moment.functions.R

In R-Finance/PortfolioAnalytics: Portfolio Analysis, including Numerical Methods for Optimization of Portfolios

###############################################################################
# R (http://r-project.org/) Numeric Methods for Optimization of Portfolios
#
# Copyright (c) 2004-2014 Brian G. Peterson, Peter Carl, Ross Bennett, Kris Boudt
#
# This library is distributed under the terms of the GNU Public License (GPL)
# for full details see the file COPYING
#
# $Id$
#
###############################################################################

#' compute comoments for use by lower level optimization functions when the conditional covariance matrix is a CCC GARCH model
#' 
#' it first estimates the conditional GARCH variances, then filters out the 
#' time-varying volatility and estimates the higher order comoments on the innovations 
#' rescaled such that their unconditional covariance matrix is the conditional covariance matrix forecast
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
#' @param momentargs list containing arguments to be passed down to lower level functions, default NULL
#' @param \dots any other passthru parameters
#' @export
CCCgarch.MM = function(R, momentargs = NULL , ... )
{
    stopifnot("package:fGarch" %in% search() || require("fGarch",quietly=TRUE))
    if (!hasArg(momentargs) | is.null(momentargs)) 
        momentargs <- list()
    cAssets = ncol(R)
    T = nrow(R)
    if (!hasArg(mu)){ 
        mu = apply(R, 2, "mean")
    }else{ mu = match.call(expand.dots = TRUE)$mu }
    R = R - matrix( rep(mu,T) , nrow = T , byrow = TRUE )
    momentargs$mu = mu
    S = nextS = c();
    for( i in 1:cAssets ){
       gout =  garchFit(formula ~ garch(1,1), data = R[,i],include.mean = F, cond.dist="QMLE", trace = FALSE )
       if( as.vector(gout@fit$coef["alpha1"]) < 0.01 ){
               sigmat = rep( sd( as.vector(R[,i])), length(R[,i]) ); nextSt = sd( as.vector(R[,i]))
        }else{
               sigmat = gout@sigma.t; nextSt = predict(gout)[1,3]
        }
        S = cbind( S , sigmat); nextS = c(nextS,nextSt)
    }
    U = R/S; #filtered out time-varying volatility
    if (!hasArg(clean)){ 
        clean = match.call(expand.dots = TRUE)$clean
    }else{ clean = NULL }
    if(!is.null(clean)){ 
        cleanU <- try(Return.clean(U, method = clean))
        if (!inherits(cleanU, "try-error")) { U = cleanU }
    }
    Rcor = cor(U)
    D = diag( nextS ,ncol=cAssets )
    momentargs$sigma = D%*%Rcor%*%D
    # set volatility of all U to last observation, such that cov(rescaled U)=sigma 
    uncS = sqrt(diag( cov(U) ))
    U = U*matrix( rep(nextS/uncS,T  ) , ncol = cAssets , byrow = T )
    momentargs$m3 = PerformanceAnalytics:::M3.MM(U)
    momentargs$m4 = PerformanceAnalytics:::M4.MM(U)
    return(momentargs)
} 

#' set portfolio moments for use by lower level optimization functions
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
#' @param constraints an object of type "constraints" specifying the constraints for the optimization, see \code{\link{constraint}}
#' @param momentargs list containing arguments to be passed down to lower level functions, default NULL
#' @param \dots any other passthru parameters
#' @export
set.portfolio.moments_v1 <- function(R, constraints, momentargs=NULL,...){

    if(!hasArg(momentargs) | is.null(momentargs)) momentargs<-list()
    if(is.null(constraints$objectives)) {
        warning("no objectives specified in constraints")
        next()
    } else {

        lcl <- grep('garch', constraints)
        if (!identical(lcl, integer(0))) {
            for (objective in constraints[lcl]) {
                objective = unlist(objective)
                if( is.null( objective$garch ) ) next
                if (objective$garch){
                   if (is.null(momentargs$mu)|is.null(momentargs$sigma)|is.null(momentargs$m3)|is.null(momentargs$m4))
                   {
                        momentargs =  CCCgarch.MM(R,clean=objective$arguments.clean,...)
                   }
               }
           }
        }


        lcl<-grep('clean',constraints)
        if(!identical(lcl,integer(0))) {
            for (objective in constraints[lcl]){
                objective = unlist(objective)
                #if(!is.null(objective$arguments$clean)) {
                if (!is.null(objective$arguments.clean)){
                   	if (is.null(momentargs$mu)|is.null(momentargs$sigma)|is.null(momentargs$m3)|is.null(momentargs$m4))
                   	{
                       	# cleanR<-try(Return.clean(R,method=objective$arguments$clean))
                       	cleanR <- try(Return.clean(R, method = objective$arguments.clean,...))
                    	if(!inherits(cleanR,"try-error")) {
                        	momentargs$mu = matrix( as.vector(apply(cleanR,2,'mean')),ncol=1);
                        	momentargs$sigma = cov(cleanR);
                        	momentargs$m3 = PerformanceAnalytics:::M3.MM(cleanR)
                        	momentargs$m4 = PerformanceAnalytics:::M4.MM(cleanR)
                        	#' FIXME NOTE: this isn't perfect as it overwrites the moments for all objectives, not just one with clean='boudt'
                    	}
                	}
            	}    
        	}
        }
        for (objective in constraints$objectives){
            switch(objective$name,
                    sd =,
                    StdDev = { 
                        if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean', na.rm=TRUE)),ncol=1);
                        if(is.null(momentargs$sigma)) momentargs$sigma = cov(R, use='pairwise.complete.obs')
                    },
                    var =,
                    mVaR =,
                    VaR = {
                        if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean')),ncol=1);
                        if(is.null(momentargs$sigma)) momentargs$sigma = cov(R)
                        if(is.null(momentargs$m3)) momentargs$m3 = PerformanceAnalytics:::M3.MM(R)
                        if(is.null(momentargs$m4)) momentargs$m4 = PerformanceAnalytics:::M4.MM(R)
                    },
                    es =,
                    mES =,
                    CVaR =,
                    cVaR =,
                    ES = {
                        if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean')),ncol=1);
                        if(is.null(momentargs$sigma)) momentargs$sigma = cov(R)
                        if(is.null(momentargs$m3)) momentargs$m3 = PerformanceAnalytics:::M3.MM(R)
                        if(is.null(momentargs$m4)) momentargs$m4 = PerformanceAnalytics:::M4.MM(R)
                    }
            ) # end switch on objectives    
        }    
    }    
    return(momentargs)
}

#' set portfolio moments for use by lower level optimization functions
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
#' @param portfolio an object of type "portfolio" specifying the constraints and objectives for the optimization, see \code{\link{portfolio.spec}}
#' @param momentargs list containing arguments to be passed down to lower level functions, default NULL
#' @param \dots any other passthru parameters
#' @aliases set.portfolio.moments
#' @rdname set.portfolio.moments
#' @export
set.portfolio.moments_v2 <- function(R, portfolio, momentargs=NULL,...){
  
  if(!hasArg(momentargs) | is.null(momentargs)) momentargs<-list()
  if(is.null(portfolio$objectives)) {
    warning("no objectives specified in portfolio")
    next()
  } else {
    
    # How would this be specified in the new portfolio.spec? As a constraint or in the portfolio part?
    # 
    lcl <- grep('garch', portfolio)
    if (!identical(lcl, integer(0))) {
      for (objective in portfolio[lcl]) {
        objective = unlist(objective)
        if( is.null( objective$garch ) ) next
        if (objective$garch){
          if (is.null(momentargs$mu)|is.null(momentargs$sigma)|is.null(momentargs$m3)|is.null(momentargs$m4))
          {
            momentargs =  CCCgarch.MM(R,clean=objective$arguments.clean,...)
          }
        }
      }
    }
    
    
    lcl<-grep('clean',portfolio)
    if(!identical(lcl,integer(0))) {
      for (objective in portfolio[lcl]){
        objective = unlist(objective)
        #if(!is.null(objective$arguments$clean)) {
        if (!is.null(objective$arguments.clean)){
          if (is.null(momentargs$mu)|is.null(momentargs$sigma)|is.null(momentargs$m3)|is.null(momentargs$m4))
          {
            # cleanR<-try(Return.clean(R,method=objective$arguments$clean))
            cleanR <- try(Return.clean(R, method = objective$arguments.clean,...))
            if(!inherits(cleanR,"try-error")) {
              momentargs$mu = matrix( as.vector(apply(cleanR,2,'mean')),ncol=1);
              momentargs$sigma = cov(cleanR);
              momentargs$m3 = PerformanceAnalytics:::M3.MM(cleanR)
              momentargs$m4 = PerformanceAnalytics:::M4.MM(cleanR)
              #' FIXME NOTE: this isn't perfect as it overwrites the moments for all objectives, not just one with clean='boudt'
            }
          }
        }    
      }
    }
    for (objective in portfolio$objectives){
      switch(objective$name,
             mean = {
               if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean', na.rm=TRUE)),ncol=1)
               },
             var =,
             sd =,
             StdDev = { 
               if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean', na.rm=TRUE)),ncol=1);
               if(is.null(momentargs$sigma)) momentargs$sigma = cov(R, use='pairwise.complete.obs')
             },
             mVaR =,
             VaR = {
               if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean')),ncol=1);
               if(is.null(momentargs$sigma)) momentargs$sigma = cov(R)
               if(is.null(momentargs$m3)) momentargs$m3 = PerformanceAnalytics:::M3.MM(R)
               if(is.null(momentargs$m4)) momentargs$m4 = PerformanceAnalytics:::M4.MM(R)
             },
             es =,
             mES =,
             CVaR =,
             cVaR =,
             ETL=,
             mETL=,
             ES = {
               # We don't want to calculate these moments if we have an ES 
               # objective and are solving as an LP problem.
               if(hasArg(ROI)) ROI=match.call(expand.dots=TRUE)$ROI else ROI=FALSE
               if(!ROI){
                 if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean')),ncol=1);
                 if(is.null(momentargs$sigma)) momentargs$sigma = cov(R)
                 if(is.null(momentargs$m3)) momentargs$m3 = PerformanceAnalytics:::M3.MM(R)
                 if(is.null(momentargs$m4)) momentargs$m4 = PerformanceAnalytics:::M4.MM(R)
               }
             }
      ) # end switch on objectives    
    }    
  }    
  return(momentargs)
}

# Alias for set.portfolio.moments
#' @export
set.portfolio.moments <- set.portfolio.moments_v2

garch.mm <- function(R,mu_ts, covlist,momentargs=list(),...) {
    #momentargs<-list()
    #momentargs$mu<-mu_ts[last(index(R)),]
    momentargs$mu<-mu_ts[last(index(R)),]
    
    momentargs$sigma<-covlist[as.character(last(index(R)))]
    if(is.null(momentargs$m3)) momentargs$m3 = PerformanceAnalytics:::M3.MM(R)
    if(is.null(momentargs$m4)) momentargs$m4 = PerformanceAnalytics:::M4.MM(R)
    return(momentargs)
}

#' Portfolio Moments
#' 
#' Set portfolio moments for use by lower level optimization functions using
#' a statistical factor model based on the work of Kris Boudt.
#' 
#' @note If any of the objectives in the \code{portfolio} object have 
#' \code{clean} as an argument, the cleaned returns are used to fit the model. 
#' 
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of 
#' asset returns
#' @param portfolio an object of type \code{portfolio} specifying the 
#' constraints and objectives for the optimization, see 
#' \code{\link{portfolio.spec}}
#' @param momentargs list containing arguments to be passed down to lower level 
#' functions, default NULL
#' @param k number of factors used for fitting statistical factor model
#' @param \dots any other passthru parameters
#' @export
portfolio.moments.boudt <- function(R, portfolio, momentargs=NULL, k=1, ...){
  
  # Fit the statistical factor model
  # If any of the objectives have clean as an argument, we fit the factor
  # model with cleaned returns. Is this the desired behavior we want?
  clean <- unlist(lapply(portfolio$objectives, function(x) x$arguments$clean))
  if(!is.null(clean)){
    if(length(unique(clean)) > 1){
      warning(paste("Multiple methods detected for cleaning returns, default to use clean =", tmp[1]))
    }
    # This sets R as the cleaned returns for the rest of the function
    # This is proably fine since the only other place R is used is for the 
    # mu estimate
    R <- Return.clean(R, method=clean[1])
  }
  fit <- statistical.factor.model(R=R, k=k)
  
  if(!hasArg(momentargs) | is.null(momentargs)) momentargs<-list()
  if(is.null(portfolio$objectives)) {
    warning("no objectives specified in portfolio")
    next()
  } else {
    for (objective in portfolio$objectives){
      switch(objective$name,
             mean = {
               if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean', na.rm=TRUE)),ncol=1)
             },
             var =,
             sd =,
             StdDev = { 
               if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean', na.rm=TRUE)),ncol=1)
               if(is.null(momentargs$sigma)) momentargs$sigma = extractCovariance(fit)
             },
             mVaR =,
             VaR = {
               if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean')),ncol=1)
               if(is.null(momentargs$sigma)) momentargs$sigma = extractCovariance(fit)
               if(is.null(momentargs$m3)) momentargs$m3 = extractCoskewness(fit)
               if(is.null(momentargs$m4)) momentargs$m4 = extractCokurtosis(fit)
             },
             es =,
             mES =,
             CVaR =,
             cVaR =,
             ETL=,
             mETL=,
             ES = {
               # We don't want to calculate these moments if we have an ES 
               # objective and are solving as an LP problem.
               if(hasArg(ROI)) ROI=match.call(expand.dots=TRUE)$ROI else ROI=FALSE
               if(!ROI){
                 if(is.null(momentargs$mu)) momentargs$mu = matrix( as.vector(apply(R,2,'mean')),ncol=1)
                 if(is.null(momentargs$sigma)) momentargs$sigma = extractCovariance(fit)
                 if(is.null(momentargs$m3)) momentargs$m3 = extractCoskewness(fit)
                 if(is.null(momentargs$m4)) momentargs$m4 = extractCokurtosis(fit)
               }
             }
      ) # end switch on objectives    
    }    
  }    
  return(momentargs)
}

###############################################################################
# $Id$
###############################################################################