R/StatAndUtil.R

In alejandro-sotolongo/InvTools:

#' @title Stanardize CME IDs
#' @param x = vector of CME IDs
#' @details Removes '_Ex', '_M', and other idiosyncratic CME Id issues, also
#' changes to upper case to remove case sensitivity
#' @export
standardCmeId <- function(x) {
  x <- x %>%
    toupper() %>%
    gsub(pattern = "_EX", replacement = "") %>%
    gsub(pattern = "_M", replacement = "") %>%
    gsub(pattern = "DJUBSTR INDEX", replacement = "BCOMTR INDEX") %>%
    gsub(pattern = "CSLABGS  INDEX", replacement = "CSLABGS INDEX") %>%
    gsub(pattern = "CAUS_BOGE INDEX", replacement = "CA_EXUS_BOGE INDEX") %>%
    gsub(pattern = "CAUS_PE INDEX", replacement = "CA_EXUS_PE INDEX")
  return(x)
}

#' @title Return to Price
#' @param ret return time-series with date in first column
#' @export
retToPrice <- function(ret) {
  ret <- checkTimeSeries(ret)
  addOne <- function(x) x + 1
  price <- ret %>%
    replace(is.na(.), 0) %>%
    mutate_at(vars(-date), addOne) %>%
    mutate_at(vars(-date), cumprod)
  return(price)
}

#' @title Price to Return
#' @param price price time-series with date in first column
#' @export
priceToRet <- function(price) {
  price <- checkTimeSeries (price)
  calcDelta <- function(x) x / lag(x, 1) - 1
  ret <- price %>%
    tidyr::fill(-date) %>% # replace missing values with last value
    mutate_at(vars(-date), calcDelta) %>%
    slice(2:nrow(price))
  return(ret)
}

#' @title Check time-series
#' @param time_series data.frame or tibble of time-series data with dates
#' in the first column
#' @param check_date option to arrange the time-series by date
#' (oldest to newest)
#' @export
checkTimeSeries <- function(time_series, arrange_date = TRUE) {
  if (!'Date' %in% class(time_series[[1]])) {
    warning('First column in time-series needs to be date observations.')
    try(dt <- as.Date(time_series[[1]]))
    if (class(dt) == 'try-error') {
      stop('could not convert first column of time-series to date')
    } else {
      time_series[[1]] <- dt
    }
  }
  if (colnames(time_series)[1] != 'date') {
    warning('renaming time_series first column \'date\'')
    colnames(time_series)[1] <- 'date'
  }
  if (arrange_date) {
    time_series <- time_series %>%
      arrange(date)
  }
  return(time_series)
}

#' @title Check frequency
#' @param freq = string to specify time-series frequency
#' @export
checkFreq <- function(freq) {
  freq <- toupper(freq)
  if (!freq %in% c('D', 'W', 'M', 'Q', 'A')) {
    stop('freq misspecified, must be \'D\', \'W\', \'M\', \'Q\', or \'A\'')
  }
  return(freq)
}

#' @title Convert asset objects to a return matrix
#' @param asset list of Assets
#' @param freq optional parameter to convert all assets to a frequency
#' @param trunc_incept optional parameter to truncate the time-series at
#'    the maximum (most recent) common inception
#' @export
assetToRet <- function(asset, freq = NULL, trunc_incept = TRUE) {
  if (!is.null(freq)) {
    asset <- changeAssetFreq(asset, freq)
  }
  retl <- lapply(asset, '[[', 'time_series')
  ret <- Reduce(function(x, y) full_join(x, y, by = 'date'), retl)
  ret <- ret %>% arrange(date)
  if (trunc_incept) {
    comm_incept <- sapply(retl, getIncept) %>%
      as.Date(origin = '1970-01-01') %>%
      max()
    ret <- ret %>% filter(date >= comm_incept)
  }
  return(ret)
}

#' @title Change frequency of a list of Assets
#' @param asset_in list of Assets
#' @param freq string to represent frequency
#' @export
changeAssetFreq <- function(asset_in, freq) {
  asset_out <- list()
  for (i_asset in 1:length(asset_in)) {
    asset_out[[i_asset]] <- asset_in[[i_asset]]$clone(deep = FALSE)
    asset_out[[i_asset]]$changeFreq(freq)
  }
  return(asset_out)
}

#' @title Utility to get inception from time-seires
#' @param x time-series with dates in first column named 'date'
#' @export
getIncept <- function(x) {
  x$date[1]
}

#' @title Change frequency of a time-series data.frame
#' @param time_series data.frame with dates in first column
#' @param freq = desired frequency
#' @param data_type = either RETURN for % change or LEVEL
#' @export
changeTimeSeriesFreq <- function(time_series, freq, data_type = 'RETURN') {
  time_series <- checkTimeSeries(time_series)
  freq <- checkFreq(freq)
  if (freq == 'D') {
    return(time_series)
  }
  if (data_type == 'RETURN') {
    time_series <- retToPrice(time_series)
  }
  date <- time_series$date
  new_date <- switch(
    freq,
    D = paste0(lubridate::day(date), '-', lubridate::year(date)),
    W = paste0(lubridate::isoweek(date), '-', lubridate::year(date)),
    M = paste0(lubridate::month(date), '-', lubridate::year(date)),
    Q = paste0(lubridate::quarter(date), '-', lubridate::year(date)),
    A = paste0(lubridate::year(date), '-', lubridate::year(date)),
  )
  time_series <- time_series %>%
    mutate(new_date = new_date) %>%
    group_by(new_date) %>%
    filter(date == max(date)) %>%
    ungroup() %>%
    select(-new_date)
  if (data_type == 'RETURN') {
    time_series <- priceToRet(time_series)
  }
  if (freq %in% c('M', 'Q', 'A')) {
    time_series <- time_series %>%
      mutate(date = lubridate::ceiling_date(date, unit = 'months') - 1)
  }
  return(time_series)
}

#' @title Truncate time-series
#' @param time_series = data.frame with dates in first column
#' @param date_start = date to start truncation, leave NULL to not truncate the start
#' @param date_end = date to end truncation, leave NULL to not truncate the end
#' @export
trunTimeSeries <- function(time_series, date_start = NULL, date_end = NULL) {
  time_series <- checkTimeSeries(time_series)
  if (is.null(date_end)) {
    date_end <- Sys.Date()
  }
  if (is.null(date_start)) {
    date_start <- as.Date('0000-01-01')
  }
  time_series <- time_series %>%
    filter(date >= date_start) %>%
    filter(date <= date_end)
  return(time_series)
}

#' @title Convert numeric data to string % format
#' @export
fPercent <- function(x, digits = 2) {
  x_fmt <- formatC(x * 100, digits = digits, format = 'f')
  x_fmt_abs <- formatC(abs(x) * 100, digits = digits, format = 'f')
  x_per <- paste0(x_fmt, '%')
  less_0 <- x < 0
  less_0[is.na(less_0)] <- FALSE
  x_per[less_0] <- paste0('(', x_fmt_abs[less_0], '%)')
  x_per[x_per == ' NA%'] <- '-'
  return(x_per)
}

#' @title Convert numeric data to 0.00 character format
#' @export
fNum <- function(x, digits = 2) {
  x <- formatC(x, digits = 2, format = 'f')
  x[x == ' NA'] <- '-'
  return(x)
}


#' @title Convert character frequency to corresponding numeric scale
#' @export
freqToScale <- function(x) {
  switch(
    toupper(x),
    D = 252,
    W = 52,
    M = 12,
    Q = 4,
    S = 2,
    A = 1
  )
}

#' @title Net a constant fee from a time-series
#'
#' @param ret return time-series with dates in first column and values in the
#' second column
#' @param annual_fee annual_fee to deduct, e.g., 2% = 0.02
#' @param ret_freq return frequency
#'
#' @details The annual fee is scaled by the frequency and subtracted from each
#' observation in the time-series, i.e., (1 + annual_fee)^(1 / freq_scale) - 1
#'
#' @export
netFee <- function(ret, annual_fee, ret_freq = 'm') {
  mult <- freqToScale(ret_freq)
  fee_vec <- rep((1 + annual_fee)^(1 / mult) - 1, nrow(ret))
  ret[, 2] <- ret[, 2] - fee_vec
  return(ret)
}

#' @title Calendar time helper
#'
#' @param xwin calendar window
#' @param asof as of date
#'
#' @return beginning date of window
#'
#' @export
calTime <- function(xwin = c('dtd', 'wtd', 'mtd', 'qtd', 'ytd', 'ttm',
                             '3 yr', '5 yr', '10 yr', '20 yr'),
                    asof = Sys.Date()) {
  xwin <- xwin[1] %>% tolower()
  xyear <- lubridate::year(asof)
  xmon <- lubridate::month(asof)
  xday <- lubridate::day(asof)
  switch (xwin,
          dtd = asof,
          wtd = lubridate::floor_date(asof, 'week') + 1,
          mtd = lubridate::floor_date(asof, 'month'),
          qtd = lubridate::floor_date(asof, 'quarter'),
          ytd = lubridate::floor_date(asof, 'year'),
          ttm = paste0(xyear - 1, '-', xmon, '-', xday) %>% as.Date(),
          '3 yr' = paste0(xyear - 3, '-', xmon, '-', xday) %>% as.Date(),
          '5 yr' = paste0(xyear - 5, '-', xmon, '-', xday) %>% as.Date(),
          '10 yr' = paste0(xyear - 10, '-', xmon, '-', xday) %>% as.Date(),
          '20 yr' = paste0(xyear - 20, '-', xmon, '-', xday) %>% as.Date()
  )
}

#' @title Frequency abbreviation to name
#'
#' @param x frequency abbreviation
#'
#' @export
freqToStr <- function(x) {
  switch(toupper(x),
         D = 'Daily',
         W = 'Weekly',
         M = 'Monthly',
         Q = 'Quarterly',
         A = 'Annual')
}

#' @title Omega ratio
#' @param ret data.frame or tibble organized by columns with dates in first
#' column
#' @param mar minimum acceptable return
#'
#' @details Calculates based on emperical density. CDF of returns are subsetted
#' by CDF of return above MAR (UpCDF) and below MAR (DownCDF). Then omega =
#' (UpCDF \* Returns above MAR) / (DownCDF \* Returns below MAR).
#'
#' @export
omegaRatio <- function(ret, mar = 0) {
  ret <- ret %>%
    replace(is.na(.), 0)
  den <- apply(ret[, 2:ncol(ret)], 2, density)
  sapply(den, omegaCalc, mar = mar)
}

#' @title Utility to calculate omega ratio
#'
#' @param den density result from base density() function
#' @param mar minimum acceptable return
#'
#' @note see omegaRatio
#'
#' @export
omegaCalc <- function(den, mar = 0) {
  abs(sum(den$y[den$x >= mar] * den$x[den$x >= mar])) /
    abs(sum(den$y[den$x < mar] * den$x[den$x < mar]))
}

#' @title Sortino ratio
#'
#' @param ret data.frame or tibble organized by columns with dates in first
#' column
#' @param rf risk-free time-series similar to \code{ret}
#' @param freq time-series frequency
#'
#' @export
sortinoRatio <- function(ret, rf, freq) {
  ann_ret <- annRet(ret, freq)
  ann_rf <- annRet(rf, freq)
  ann_down_sd <- apply(ret[, 2:ncol(ret)], 2, downSd) *
    sqrt(freqToScale(freq))
  (ann_ret - ann_rf) / ann_down_sd
}

#' @title Downside standard deviation
#'
#' @param ret data.frame or tibble organized by columns with dates in the first
#' column
#'
#' @export
downSd <- function(ret) {
  ret_down <- ret[ret < 0]
  sd(ret_down, na.rm = TRUE)
}

#' @title Annualzed (geomtric) return
#'
#' @param ret data.frame or tibble organized by columns with dates in the first
#' column
#' @param freq return frequency to scale (annualize) return(s)
#'
#' @export
annRet <- function(ret, freq) {
  mult <- freqToScale(freq)
  gross_prod <- apply(1 + ret[, 2:ncol(ret), drop = FALSE], 2, prod, na.rm = TRUE)
  gross_prod^(mult / nrow(ret)) - 1
}

#' @title Sharpe ratio
#'
#' @param ret data.frame or tibble organized by columns with dates in first
#' column
#' @param rf risk-free time-series similar to \code{ret}
#' @param freq time-series frequency
#'
#' @export
sharpeRatio <- function(ret, rf, freq) {
  ann_ret <- annRet(ret, freq)
  ann_rf <- annRet(rf, freq)
  ann_sd <- apply(ret[, 2:ncol(ret)], 2, sd, na.rm = TRUE) *
    sqrt(freqToScale(freq))
  (ann_ret - ann_rf) / ann_sd
}

#' @title Up and down capture ratios
#'
#' @param fund data.frame or tibble of fund return time-series (1 fund only)
#' with dates in the first column and returns in the second
#' @param bench data.frame or tibble of benchmark time-series (1 benchmark)
#' with dates in the first column and returns in the second
#'
#' @export
captureRatio <- function(fund, bench) {
  fund_up <- fund[bench[, 2] >= 0, ]
  bench_up <- bench[bench[, 2] >= 0, ]
  fund_down <- fund[bench[, 2] < 0, ]
  bench_down <- bench[bench[, 2] < 0,  ]
  res <- c(NA, NA)
  names(res) <- c('up', 'down')
  res[1] <- (prod(1 + fund_up[, 2], na.rm = TRUE) - 1) /
    (prod(1 + bench_up[, 2], na.rm = TRUE) - 1)
  res[2] <- (prod(1 + fund_down[, 2], na.rm = TRUE) - 1) /
    (prod(1 + bench_down[, 2], na.rm = TRUE) - 1)
  return(res)
}

#' @title Count of positive or negative observations in a vector
#'
#' @param x vector of observations
#' @param direct up or down: return positive or negative count
#'
#' @export
nUpDown <- function(x, direct = c('up', 'down')) {
  if (direct[1] == 'up') {
    return(sum(x >= 0, na.rm = TRUE))
  } else {
    return(sum(x < 0, na.rm = TRUE))
  }
}

#' @title Tracking error
#'
#' @export
trackingError <- function(fund, bench, freq) {
  er <- excessRet(fund, bench)
  sd(er[[2]], na.rm = TRUE) * sqrt(freqToScale(freq))
}

#' @export
infoRatio <- function(fund, bench, freq) {
  er <- excessRet(fund, bench)
  te <- sd(er[[2]], na.rm = TRUE) * sqrt(freqToScale(freq))
  annRet(er, freq) / te
}

#' @export
excessRet <- function(fund, rf) {
  dat <- left_join(ret, rf, by = 'date')
  comm_start <- max(ret$date[1], rf$date[1])
  dat <- dat %>%
    filter(date >= comm_start)
  dat[is.na(dat[, 3]), 3] <- 0
  tibble(date = dat$date, value = dat[[2]] - dat[[3]])
}

#' @export
maxDrawdown <- function(ret) {
  dd <- drawdown(ret)
  apply(dd[, 2:ncol(dd)], 2, min)
}

#' @export
drawdown <- function(ret) {

  # converts returns to drawdowns
  # ret = data.frame organized by columns, with dates in the first column

  row_zero <- data.frame(
    date = ret$date[1] - 1,
    matrix(0, nrow = 1, ncol(ret) - 1)
  )
  colnames(row_zero) <- colnames(ret)
  dd <- ret %>%
    rbind(row_zero) %>%
    arrange(date) %>%
    mutate_at(vars(-date), vecDrawdown)
  return(dd)
}

#' @export
vecDrawdown <- function(x) {

  # drawdown utility
  # x = return vector

  cum_ret <- cumprod(x + 1)
  cum_ret / cummax(cum_ret) - 1
}


#' @title Calendar returns
#'
#' @param asset_list list of Asset or Portfolio Objects
#' @param xwin window for calendar return
#' @param asof as of date
#'
#' @details periods longer than one year are annualized
#'
#' @export
calRet <- function(asset_list, xwin = c('dtd', 'wtd', 'mtd', 'qtd', 'ytd',
                                        'ttm', '3 yr', '5 yr', '10 yr', '20 yr'),
                   asof = Sys.Date()) {
  date_start <- calTime(xwin, asof)
  asset_ret <- assetToRet(asset_list, trunc_incept = FALSE) %>%
    filter(date >= date_start) %>%
    filter(date <= asof)
  cum_ret <- apply(asset_ret[, -1] + 1, 2, prod, na.rm = TRUE)
  a <- 1
  switch (xwin,
          '3 yr' = a <- 3,
          '5 yr' = a <- 5,
          '10 yr' = a <- 10,
          '20 yr' = a <- 20
  )
  cum_ret <- cum_ret^(1 / a) - 1
  incept <- sapply(asset_list, getTimeSeriesStart) %>%
    as.Date(origin = '1970-01-01')
  cum_ret[incept > date_start] <- NA
  return(cum_ret)
}

#' @export
calSd <- function(asset_list,
                  freq = 'w',
                  xwin = c('mtd', 'qtd', 'ytd', 'ttm', '3 yr', '5 yr', '10 yr',
                           '20 yr'),
                  asof = Sys.Date()) {
  freq <- toupper(freq)
  date_start <- calTime(xwin, asof)
  asset_ret <- assetToRet(asset_list) %>%
    filter(date >= date_start) %>%
    filter(date <= asof)
  asset_ret <- changeTimeSeriesFreq(asset_ret, freq)
  period_sd <- apply(asset_ret[, -1], 2, sd, na.rm = TRUE)
  period_sd <- period_sd * sqrt(freqToScale(freq))
  incept <- sapply(asset_list, getTimeSeriesStart) %>%
    as.Date(origin = '1970-01-01')
  period_sd[incept > date_start] <- NA
  return(period_sd)
}

#' @export
getTimeSeriesStart <- function(asset) {
  asset$time_series$date[1]
}

#' @title Calculate Wealth Index
#'
#' @param ret data.frame or tibble of returns with dates in the first column
#' labeled date
#' @param init_val starting value for wealth index
#'
#' @export
wealthIndex <- function(ret, init_val = 1) {
  if (colnames(ret)[1] != 'date') {
    warning('ret first column is not named \'date\', trying to force')
    colnames(ret)[1] <- 'date'
  }
  ret[is.na(ret)] <- 0
  wi <- ret
  wi[, -1] <- apply(ret[, -1, drop = FALSE] + 1, 2, cumprod) * init_val
  init_val_row <- matrix(init_val, nrow = 1, ncol = ncol(ret) - 1)
  first_row <- as_tibble(init_val_row) %>%
    add_column(date = ret$date[1] - 1, .before = 1)
  colnames(first_row) <- colnames(ret)
  wi <- rbind(first_row, wi)
  return(wi)
}

#' @title Utility for beta calc of two time-series
#'
#' @param x = data.frame, tibble, matrix, or matrix-like structure containing
#' two columns with time-series for beta calc
#'
#' @details uses max-pairwise
#'
#' @export
calcBeta <- function(x) {
  xcov <- cov(x, use = 'pairwise.complete.obs')
  xcov[2, 1] / xcov[2, 2]
}


#' @title Calculate Rolling Beta of 2 time-series
#'
#' @param ret data.frame or tibble of two time-series with dates in the first
#' column
#' @param xwin rolling window of observations for beta calculation
#'
#' @export
rollBeta <- function(ret, xwin = 63) {
  if (ncol(ret) > 3) {
    warning('ret is greater than 3 columns, taking beta of column 2 w/r/t 3')
  }
  n_obs <- nrow(ret)
  xbeta <- rep(NA, n_obs)
  for (i in xwin:n_obs) {
    xbeta[i] <- calcBeta(ret[(i - xwin + 1):i, 2:3])
  }
  tibble(date = ret[[1]], beta = xbeta)
}

#' @export
rollCor <- function(ret, xwin = 63) {
  if (ncol(ret) > 3) {
    warning('ret is greater than 3 columns, taking corr of column 2 w/r/t 3')
  }
  n_obs <- nrow(ret)
  xcor <- rep(NA, n_obs)
  for (i in xwin:n_obs) {
    xcor[i] <- cor(ret[(i - xwin + 1):i, 2], ret[(i - xwin + 1):i, 3],
                   use = 'pairwise.complete.obs')
  }
  tibble(date = ret[[1]], cor = xcor)
}

#' @export
riskWgt <- function(cap_wgt, ecov) {
  cap_wgt <- matrix(cap_wgt, nrow = length(cap_wgt), ncol = 1)
  (cap_wgt * (ecov %*% cap_wgt)) / (t(cap_wgt) %*% ecov %*% cap_wgt)[1]
}

#' @title Generates FFLA Factors
#'
#' @param ts TickStore object
#'
#' @export
FFLAFactorModel <- function(ts) {
  tick <- c('MXWD Index', 'FFA-MOM-D', 'FFA-VALUE-D', 'FFA-QUALITY-D',
            'FFA-CARRY-D', 'FFA-MACRO-D', 'FFA_US_TERM_D', 'FFA_US_DEF_D',
            'BCOMTR Index', 'DXY Curncy')
  fact_list <- lapply(tick, ts$genAsset)
  names(fact_list) <- tick
  fact_list$`DXY Curncy`$time_series <- priceToRet(fact_list$`DXY Curncy`$time_series)
  return(fact_list)
}

#' @export
factRiskWgt <- function(xbeta, cov_fact, y_variance) {

  # calculates factor contribution to risk
  # xbeta = factor betas
  # cov_fact = factor covariance matrix
  # y_variance = variance of y variable

  xbeta <- matrix(xbeta, nrow = length(xbeta), ncol = 1)
  risk_wgt <- (xbeta * (cov_fact %*% xbeta)) / y_variance[1]
  return(risk_wgt)
}

#' @export
covEWMA <- function(ret, freq) {
  ret <- checkTimeSeries(ret)
  retn <- ret[, 2:ncol(ret)] %>%
    replace(is.na(.), 0)
  freq <- checkFreq(freq)
  mult <- freqToScale(freq)
  n_obs <- nrow(retn)
  n_assets <- ncol(retn)
  mu <- colMeans(retn)
  mu_mat <- matrix(mu, nrow = n_obs, ncol = n_assets, byrow = TRUE)
  retc <- retn - mu_mat
  lambda <- 1 - 2 / (n_obs + 1)
  exp_cov <- matrix(0, nrow = n_assets, ncol = n_assets)
  for (t in 1:n_obs) {
    roll_ret <- retc[t, ]
    eps <- t(roll_ret) %*% as.numeric(roll_ret) * mult
    exp_cov <- lambda * exp_cov + (1 - lambda) * eps
  }
  return(exp_cov)
}

#' @export
vecGeoRet <- function(x) {
  prod((x + 1), na.rm = TRUE)^(1 / length(x)) - 1
}