Exponentially weighted covariance in an Equal Risk Contribution portfolio optimisation problem

library(FRAPO)
library(tidyverse)
library(roll)

# Set chart options
options(repr.plot.width = 14, repr.plot.height=7, warn=-1)
theme_set(theme_bw())
theme_update(text = element_text(size = 20))

# source helper functions
source("erc_helpers.R")

# get data
tickers <- c("TLT", "SPY", "GLD", "VWO", "EMB")
prices <- yahoo_prices(
    tickers,
    from_date = "2008-01-01",
    to_date = "2023-11-16"
  ) %>% 
  arrange(Date, Ticker)
 
head(prices)

# we'll need simple and log returns
returns <- prices %>%
  select(Ticker, Date, Adj.Close) %>%
  arrange(Date) %>%
  group_by(Ticker) %>%
  mutate(
    simple_return = (Adj.Close - dplyr::lag(Adj.Close))/dplyr::lag(Adj.Close),
    fwd_simple_return = dplyr::lead(simple_return),
    log_return = log(Adj.Close/dplyr::lag(Adj.Close)),
    fwd_log_return = dplyr::lead(log_return)
  )

# make wide returns dataframe
wide_rets <- returns %>%
  select(Ticker, Date, log_return) %>%
  pivot_wider(Date, names_from = Ticker, values_from = log_return) %>%
  na.omit()

# make pairwise returns dataframe 
pairwise_returns <- returns %>%
  full_join(returns, by = "Date") %>%
  ungroup() %>% 
  # get all combinations (tickers) and remove duplicate combos (eg A-AA, AA-A)
  mutate(tickers = ifelse(Ticker.x < Ticker.y, glue("{Ticker.x}, {Ticker.y}"), glue("{Ticker.y}, {Ticker.x}"))) %>%
  distinct(Date, tickers, .keep_all = TRUE) %>%
  # arrange for rolling pairwise ewma correlations
  arrange(Date, .by_group = TRUE) %>%
  na.omit() 

head(pairwise_returns)

# wrangle ewma covariances
lambda <- 0.99
wdw <- 200  # a reasonable approximation of an "equivalent" window length is 2/(1-lambda)
ewma_covs <- pairwise_returns %>%
  group_by(tickers) %>%
  arrange(Date, .by_group = TRUE) %>%
  mutate(cov = ewma_cov(log_return.x, log_return.y, lambda = lambda, initialisation_wdw = wdw)) %>%
  select(Date, tickers, cov) %>%
  na.omit()

ewma_covs_long <- ewma_covs %>%
  separate(tickers, into = c("ticker.x", "ticker.y")) 

# daily ERC weights from exponentially weighted covariance matrix
days <- ewma_covs_long %>% distinct(Date) %>% pull()
tickers <- ewma_covs_long %>% distinct(ticker.x) %>% pull()
num_assets <- length(tickers)

erc_wgts <- vector(mode = "list", length = length(days))
non_psd <- vector(mode = "numeric", length = length(days))
i <- 1
for(day in days) {
  
  cov_mat <- ewma_covs_long %>%
    dplyr::filter(Date == day) %>%
    recover_covmat(tickers, num_assets)

  # check that covmat is PSD, if not convert to nearest
  if(any(round(eigen(cov_mat)$values, 10) < 0)) {
    capture.output(psd <- Matrix::nearPD(cov_mat, keepDiag = TRUE, ensureSymmetry = TRUE, base.matrix = TRUE, maxit = 100), file = "NULL", type = "message")
    cov_mat <- psd$mat
    # if still non-psd or non symmetrical (eg due to non-convergence of nearPD), use hacky method
    if(any(round(eigen(cov_mat)$values, 10) < 0) || !isSymmetric(cov_mat)) {
      cov_mat <- get_near_psd(cov_mat)
      # sometimes this will still fail as we force symmetry - just note this for now
      if(any(round(eigen(cov_mat)$values, 10) < 0)) {
        non_psd[[i]] <- day
      }
    }
  }

  # ERC weights based on covariance matrix 
  capture.output(these_wgts <- FRAPO::PERC(Sigma = cov_mat), file = "NULL")  # too much output...
  erc_wgts[[i]] <- Weights(these_wgts) 
  i <- i + 1
}

erc_wgts <- wrangle_weights(erc_wgts, dates = ewma_covs_long %>% distinct(Date) %>% arrange(Date), tickers)

options(repr.plot.width = 14, repr.plot.height=10)
show_performance(erc_wgts, title = "ERC Portfolio Performance with Exponentially Weighted Covariance Matrix")

# calculate rolling pairwise covariances over period as a long dataframe
wdw = 200 # a reasonable approximation of an "equivalent" window length is 2/(1-lambda)
rolling_covs <- pairwise_returns %>%
  group_by(tickers) %>%
  arrange(Date, .by_group = TRUE) %>%
  mutate(cov = roll::roll_cov(  # replace with cov here, in previous function, and in previous use
    x = log_return.x, 
    y = log_return.y, 
    width = wdw
  )) %>% 
  select(Date, tickers, cov) %>%
  na.omit()

rolling_covs_long <- rolling_covs %>%
  separate(tickers, into = c("ticker.x", "ticker.y")) 

# daily ERC weights from exponentially weighted covariance matrix
days <- rolling_covs_long %>% distinct(Date) %>% pull()
tickers <- rolling_covs_long %>% distinct(ticker.x) %>% pull()
num_assets <- length(tickers)

erc_wgts <- vector(mode = "list", length = length(days))
non_psd <- vector(mode = "numeric", length = length(days))
i <- 1
for(day in days) {
  
  cov_mat <- rolling_covs_long %>%
    dplyr::filter(Date == day) %>%
    recover_covmat(tickers, num_assets)

  # check that covmat is PSD, if not convert to nearest
  if(any(round(eigen(cov_mat)$values, 10) < 0)) {
    capture.output(psd <- Matrix::nearPD(cov_mat, keepDiag = TRUE, ensureSymmetry = TRUE, base.matrix = TRUE, maxit = 100), file = "NULL", type = "message")
    cov_mat <- psd$mat
    # if still non-psd or non symmetrical (eg due to non-convergence of nearPD), use hacky method
    if(any(round(eigen(cov_mat)$values, 10) < 0) || !isSymmetric(cov_mat)) {
      cov_mat <- get_near_psd(cov_mat)
      # sometimes this will still fail as we force symmetry - just note this for now
      if(any(round(eigen(cov_mat)$values, 10) < 0)) {
        non_psd[[i]] <- day
      }
    }
  }

  # ERC weights based on covariance matrix 
  capture.output(these_wgts <- FRAPO::PERC(Sigma = cov_mat), file = "NULL")  # too much output...
  erc_wgts[[i]] <- Weights(these_wgts) 
  i <- i + 1
}

erc_wgts <- wrangle_weights(erc_wgts, dates = rolling_covs_long %>% distinct(Date) %>% arrange(Date), tickers)

show_performance(erc_wgts, title = "ERC Portfolio Performance with Rolling Window Covariance Matrix")

eq_wgts <- returns %>% 
  select(Date, Ticker, simple_return, fwd_simple_return) %>% 
  mutate(weight = 100/3)

show_performance(eq_wgts, title = "Equally Weighted Portfolio Performance")