R/VaR.R

In zhangkd5/riskCalc: Risk Calculation System -- MATH G4082 Project

#' Calculate VaR with Various Methods
#' 
#' This generic function computes VaR and creates an S3 class object that contains data and VaR results.
#' 
#' @param data an object that contains structured data of stocks, options, or a portfolio.
#' @param TimeLen the time horizon in days that the VaR is defined over.
#' @param conf.level the confidence level of VaR
#' @param method the method used to calculate VaR. Must be one of \code{"parametric"} (default), 
#' \code{"historical"}, or \code{"Monte-Carlo"}, can be abbreviated.
#' @param lookback.years the time horizon in years that historical simulation is based on.
#' @param iteration the number of iteration Monte-Carlo simulation to be performed.
#' @param ... arguments to be passed to or from methods.
#' @param V0 the total value which VaR is measured upon.
#'  
#' @return 
#' An object of class \code{"VaR"}.
#' 
#' @seealso
#' \code{\link{backtest}} to perform VaR backtesting.
#' 
#' @details
#' This function is able to compute portfolio VaR for portfolio data, and is also able to 
#' compute VaR for each individual stocks and/or options given stocks or options data. 
#' 
#' In parametric method and Monte-Carlo method, the portfolio components should be firstly calibrated
#' to a GBM model. In parametric method, the portfolio is assumed normally distributed and the options are
#' linearly approximated by underlying stocks. In historical and Monte-Carlo method, 
#' the option values are fully evaluated.
#' 
#' See model documentation for details about VaR calculation methods.
#' 
#' @export
VaR <- function(data, ...)
  UseMethod("VaR")

#' @rdname VaR
#' @export
VaR.portfolioData <- function(data, TimeLen=5, conf.level=0.99, 
                              method=c("parametric", "historical", "Monte-Carlo"), 
                              lookback.years=5, iteration=10000, ...)
{
  method <- match.arg(method)
  if(method!="historical" && !inherits(data, "GBM"))
    stop(gettextf("method %s requires to firstly fit GBM model to data", deparse(method)))
  
  w.raw <- matrix(data$compWeight, nrow=data$N, ncol=data$DimPf, byrow=TRUE)
  if(attr(data$compWeight, "unit")=="share") {
    for(i in seq_along(data$compNames)) {
      if(data$compType[i]=="stock")
        w.raw[,i] <- w.raw[,i] * data$price[,data$compNames[i]]
      else
        w.raw[,i] <- w.raw[,i] * data$premium[,data$compNames[i]]
    }
  }
  w.raw <- w.raw / rowSums(abs(w.raw)) * data$V0
  
  if(method=="parametric") 
  {
    opNames <- data$compNames[data$compType=="option"]
    rawIdx <- match(data$stockNames, data$compNames)    
    w <- w.raw[,rawIdx,drop=FALSE]
    w[,is.na(rawIdx)] <- 0
    for(Name in opNames) {
      idx <- match(Name, data$optionNames)
      und <- data$underlying[idx]
      delta <- BSDelta(data$price[,und], data$strike[,Name], data$rate[,2], 
                       data$impVol[,Name], data$maturity[idx], data$callput[idx], data$divRate[,und])
      i.w <- match(und, data$stockNames)
      i.raw <- match(Name, data$compNames)
      w[,i.w] <- w[,i.w] + delta * data$price[,und] * w.raw[,i.raw] / data$premium[,Name]
    }
    VaR <- as.numeric(ComputeVaR.PF(TimeLen/252, conf.level, data.matrix(data$mu[,-1,drop=FALSE]), 
                                    data.matrix(data$vol[,-1,drop=FALSE]), corr=data$corr, w))
  }
  else if(method=="historical") {
    lookback <- lookback.years * 252
    if(lookback > data$N)
      stop("cannot compute historical VaR because of insufficient data")
    
    stkRet <- matrix(nrow=data$N, ncol=data$Dim, dimnames=list(NULL, data$stockNames))
    stkRet[seq(data$N-TimeLen),] <- 
      -diff(log(data.matrix(data$price[,-1,drop=FALSE])), lag=TimeLen)
    
    VaR <- c(as.numeric(sapply(seq(data$N-lookback+1), function(nr) {
      Vhist <- numeric(lookback)
      for(i in seq_along(data$compNames)) {
        if(data$compType[i]=="stock") {
          Vhist <- Vhist + w.raw[nr,i] * (exp(stkRet[nr:(nr+lookback-1),data$compNames[i]]))
        } else { ## option
          idx <- match(data$compNames[i], data$optionNames)
          und <- data$underlying[idx]
          undPrice <- data$price[nr,und] * exp(stkRet[nr:(nr+lookback-1),und])
          newPremium <- BSPrice(undPrice, data$strike[nr,idx+1], data$rate[nr,2], data$impVol[nr,idx+1],
                                data$maturity[idx]-TimeLen/252, data$callput[idx], data$divRate[nr,und])
          Vhist <- Vhist + w.raw[nr,i] * (newPremium / data$premium[nr,idx+1])
        }
      }
      sum(w.raw[nr,]) - quantile(Vhist, prob=1-conf.level, type=4, na.rm=TRUE)
    })), rep(NA, lookback-1))
  }
  else { # Monte-Carlo
    VaR <- numeric(data$N)
    for(nr in seq(data$N)) {
      stkRet <- structure(SimulateRet(TimeLen/252, data$mu[nr,-1], data$vol[nr,-1], 
                                      data$corr[[nr]], iteration), dimnames=list(NULL, data$stockNames))
      VSim <- numeric(iteration)
      for(i in seq_along(data$compNames)) {
        if(data$compType[i]=="stock") {
          VSim <- VSim + w.raw[nr,i] * exp(stkRet[,data$compNames[i]])
        } else {
          idx <- match(data$compNames[i], data$optionNames)
          und <- data$underlying[idx]
          undPrice <- data$price[nr,und] * exp(stkRet[,und])
          newPremium <- BSPrice(undPrice, data$strike[nr,idx+1], data$rate[nr,2], data$impVol[nr,idx+1],
                                data$maturity[idx]-TimeLen/252, data$callput[idx], data$divRate[nr,und])
          VSim <- VSim + w.raw[nr,i] * (newPremium / data$premium[nr,idx+1])
        }
      }
      VaR[nr] <- sum(w.raw[nr,]) - quantile(VSim, prob=1-conf.level, type=4, na.rm=TRUE)
    }
  }
  VaR <- data.frame(Date=data$price[,1], VaR=VaR)
  r <- structure(list(data=data, VaR=VaR, TimeLen=TimeLen, conf.level=conf.level, method=method), class="VaR")
  if(method=="historical") r$lookback.years <- lookback.years
  if(method=="Monte-Carlo") r$iteration <- iteration
  return(r)
}

#' @rdname VaR
#' @export
VaR.optionsData <- function(data, TimeLen=5, conf.level=0.99, 
                            method=c("parametric", "historical", "Monte-Carlo"), V0=10000, 
                            lookback.years=5, iteration=10000, ...)
{
  r <- NextMethod()
  for(Name in data$optionNames) {
    pf <- portfolio(data, Name, 1, V0)
    r$VaR[,Name] <- VaR(pf, TimeLen, conf.level, method, lookback.years, iteration)$VaR[,2]
  }
  return(r)
}

#' @rdname VaR
#' @export
VaR.stocksData <- function(data, TimeLen=5, conf.level=0.99, 
                           method=c("parametric", "historical", "Monte-Carlo"), V0=10000, 
                           lookback.years=5, iteration=10000, ...)
{
  method <- match.arg(method)
  VaR <- data$price[,1,drop=FALSE]
  for(Name in data$stockNames) {
    pf <- portfolio(data, Name, 1, V0)
    VaR[,Name] <- VaR(pf, TimeLen, conf.level, method, lookback.years, iteration)$VaR[,2]
  }
  data$V0 <- V0
  r <- structure(list(data=data, VaR=VaR, TimeLen=TimeLen, conf.level=conf.level, method=method), class="VaR")
  if(method=="historical") r$lookback.years <- lookback.years
  if(method=="Monte-Carlo") r$iteration <- iteration
  return(r)
}


# 
# ComputeVaR.GBM <- function(TimeLen, p, mu, vol, V0) {
#   stopifnot(p>0, p<1, TimeLen>0, na.omit(vol>=0))
#   n = max(length(TimeLen), length(p), length(mu), length(vol), length(V0))
#   TimeLen = rep(TimeLen, length.out = n)
#   p = rep(p, length.out = n)
#   mu = rep(mu, length.out = n)
#   vol = rep(vol, length.out = n)
#   V0 = rep(V0, length.out = n)
#   
#   return(V0 - V0 * exp(vol*sqrt(TimeLen)*qnorm(1-p) + (mu-0.5*vol^2)*TimeLen))
# }

# SimulateVaR.GBM.Option <- function(TimeLen, p, mu, vol, 
#                                    price, strike, rate, maturity, ImpVol, callput=1, divRate=0,
#                                    V0=1, w.stock=0, w.option=1, iter=10000)
# {
#   stopifnot(TimeLen>0, p>0, p<1, all(vol>0|is.na(vol)), all(ImpVol>0|is.na(ImpVol)))
#   n = length(mu)
#   nS = V0 * w.stock / price
#   P0 = BSPrice(price, strike, rate, ImpVol, maturity, callput, divRate)
#   nP = V0 * w.option / P0
#   q = numeric(n)
#   for(i in seq(n)) {
#     St = price[i] * exp((mu[i]-vol[i]^2/2)*TimeLen + vol[i]*sqrt(TimeLen)*rnorm(iter))
#     Pt = BSPrice(St, strike[i], rate[i], ImpVol[i], maturity-TimeLen, callput, divRate[i])
#     Vt.P = nS[i] * St + nP[i] * Pt
#     q[i] = V0 - quantile(Vt.P, prob=1-p, type=4, na.rm=TRUE)
#   }
#   return(q)
# }

ComputeVaR.PF <- function(TimeLen, p, mu, vol, corr, w)
{
  n <- nrow(mu)
  if(is.vector(w, "numeric") && length(w)==ncol(mu))
    w <- t(w)
  w <- as.matrix(w)
  w <- w[rep(seq(nrow(w)),length.out=n),,drop=FALSE]
  #  w <- w / rowSums(abs(w)) * V0
  if(!is.list(corr)) corr <- rep(list(corr),length.out=n)
  V0 <- rowSums(w)
  ESt <- exp(mu*TimeLen)
  EVt <- rowSums(ESt * w)
  VarVt <- sapply(seq(n), function(nr) {
    w[nr,,drop=FALSE] %*% ((t(ESt[nr,,drop=FALSE]) %*% ESt[nr,,drop=FALSE]) * 
                             (exp((t(vol[nr,,drop=FALSE]) %*% vol[nr,,drop=FALSE])*corr[[nr]]*TimeLen)-1)) %*%
      t(w[nr,,drop=FALSE])
  })
  V0 - (EVt - sqrt(VarVt) * qnorm(p))
}

#' @import mvtnorm
SimulateRet <- function(TimeLen, mu, vol, corr, iter)
{
  mu <- as.numeric(mu)
  vol <- as.numeric(vol)
  if(any(is.na(corr)))
    return(matrix(NA_real_, nrow=iter, ncol=length(mu)))
  Z <- rmvnorm(iter, sigma = corr)
  t((mu-vol^2)*TimeLen + vol*sqrt(TimeLen)*t(Z))
}