# session  options
options(repr.plot.width = 14, repr.plot.height=7, warn = -1)

library(tidyverse)
library(tidyfit)
library(glue)
library(CVXR)
library(tibbletime)
library(roll)
library(patchwork)
pacman::p_load_current_gh("Robot-Wealth/rsims", dependencies = TRUE) 

# chart options
theme_set(theme_bw())
theme_update(text = element_text(size = 20))

# data
perps <- read_csv("https://github.com/Robot-Wealth/r-quant-recipes/raw/master/quantifying-combining-alphas/binance_perp_daily.csv")
head(perps)

# get same universe as before - top 30 by rolling 30-day dollar volume, no stables

# remove stablecoins
# list of stablecoins from defi llama
url <- "https://stablecoins.llama.fi/stablecoins?includePrices=true"
response <- httr::GET(url)

stables <- response %>%
  httr::content(as = "text", encoding = "UTF-8") %>%
  jsonlite::fromJSON(flatten = TRUE) %>%
  pluck("peggedAssets") %>%
  pull(symbol)

# sort(stables)

perps <- perps %>% 
  filter(!ticker %in% glue::glue("{stables}USDT")) 

# get the top 30 by trailing 30-day volume
# exclude from universe until we have enough data to do a covariance estimate (see later)
trading_universe_size <- 30
ew_cov_init_wdw <- 30

universe <- perps %>%
  group_by(ticker) %>% 
  # also calculate returns for later
  mutate(
    total_return_simple = funding_rate + (close - lag(close, 1))/lag(close, 1),
    total_return_log = log(1 + total_return_simple),
    total_fwd_return_simple = dplyr::lead(funding_rate, 1) + (dplyr::lead(close, 1) - close)/close,
    total_fwd_return_log = log(1 + total_fwd_return_simple),
    trail_volume = roll_mean(dollar_volume, 30)
  ) %>% 
  mutate(days_since_listing = row_number()) %>%   # the number of days we have data for each ticker
  na.omit() %>%
  ungroup() %>% 
  # exclude from volume ranking - won't matter if we use a rolling volume calculated with more than the covariance initialisation window as will drop these rows using na.omit
  mutate(to_exclude = days_since_listing < ew_cov_init_wdw) %>%   
  group_by(date) %>%
  mutate(
    volume_rank = if_else(to_exclude, NA_integer_, row_number(-trail_volume)),
    is_universe = volume_rank <= trading_universe_size,
  ) 

universe %>%
  group_by(date, is_universe) %>%
  summarize(count = n(), .groups = "drop") %>%
  ggplot(aes(x=date, y=count, color = is_universe)) + 
  geom_line() + 
  labs(
    title = 'Universe size'
)

# calculate features as before
rolling_days_since_high_20 <- purrr::possibly(
  tibbletime::rollify(
    function(x) {
      idx_of_high <- which.max(x)
      days_since_high <- length(x) - idx_of_high
      days_since_high
    }, 
    window = 20, na_value = NA), 
  otherwise = NA
)

features <- universe %>%
  group_by(ticker) %>%
  arrange(date) %>%
  mutate(
    breakout = 9.5 - rolling_days_since_high_20(close),  # puts this feature on a scale -9.5 to +9.5
    momo = close - lag(close, 10)/close,
    carry = funding_rate
  ) %>%
  ungroup() %>% 
  na.omit()

# create a model df on our universe with momo and carry features scaled into deciles
model_df <- features %>%
  filter(is_universe) %>% 
  group_by(date) %>%
  mutate(
    carry_decile = ntile(carry, 10),
    momo_decile = ntile(momo, 10),
    # also calculate demeaned return for everything in our universe each day for later
    demeaned_return = total_return_simple - mean(total_return_simple, na.rm = TRUE),
    demeaned_fwd_return = total_fwd_return_simple - mean(total_fwd_return_simple, na.rm = TRUE)
  ) %>% 
  ungroup()

  # start simulation from date we first have n tickers in the universe
start_date <- features %>%
  group_by(date, is_universe) %>%
  summarize(count = n(), .groups = "drop") %>%
  filter(count >= trading_universe_size) %>%
  head(1) %>%
  pull(date)

# use a 90-day window and refit every 10 days
is_days <- 90
step_size <- trading_universe_size*10  

# rolling model for cross-sectional features
roll_xs_coeffs_df <- model_df %>% 
  filter(date >= start_date) %>% 
  regress(
    demeaned_fwd_return ~ carry_decile + momo_decile,
    m("lm", vcov. = "HAC"),  
    .cv = "sliding_index", 
    .cv_args = list(lookback = days(is_days), step = step_size, index = "date"),  
    .force_cv = TRUE, 
    .return_slices = TRUE
  )

# rolling model for time series features 
breakout_cutoff <- 5.5  # below this level, we set our expected return to zero
roll_ts_coeffs_df <- model_df %>% 
  filter(date >= start_date) %>% 
  # setting regression weights to zero when breakout < breakout_cutoff will give these data points zero weight in estimating coefficients
  mutate(regression_weights = case_when(breakout < breakout_cutoff ~ 0, TRUE ~ 1)) %>% 
  regress(
    total_fwd_return_simple ~ breakout,
    m("lm", vcov. = "HAC"), 
    .weights = "regression_weights", 
    .cv = "sliding_index", 
    .cv_args = list(lookback = days(is_days), step = step_size, index = "date"),  
    .force_cv = TRUE, 
    .return_slices = TRUE
  )

roll_xs_coeffs_df %>% head
roll_ts_coeffs_df %>% select(-settings) %>% head

# for this to work, need to install.packages("sandwich", "lmtest")
xs_coefs <- roll_xs_coeffs_df %>% 
  coef()

xs_coefs_df <- xs_coefs %>% 
  ungroup() %>% 
  select(term, estimate, slice_id) %>% 
  pivot_wider(id_cols = slice_id, names_from = term, values_from = estimate) %>% 
  mutate(slice_id = as_date(slice_id)) %>% 
  # need to lag slice id to make it oos
  # slice_id_oos is the date we start using the parameters
  mutate(slice_id_oos = lead(slice_id)) %>% 
  rename("xs_intercept" = `(Intercept)`)

ts_coefs <- roll_ts_coeffs_df %>% 
  coef()

ts_coefs_df <- ts_coefs %>% 
  ungroup() %>% 
  select(term, estimate, slice_id) %>% 
  pivot_wider(id_cols = slice_id, names_from = term, values_from = estimate) %>% 
  mutate(slice_id = as_date(slice_id)) %>% 
  # need to lag slice id to make it oos
  # slice_id_oos is the date we start using the parameters
  mutate(slice_id_oos = lead(slice_id)) %>% 
  rename("ts_intercept" = `(Intercept)`)

xs_coefs_df %>% head
# ts_coefs_df %>% head

# plot cross-sectional estimates
xs_coefs_df %>% 
  select(-slice_id) %>% 
  pivot_longer(cols = -slice_id_oos, names_to = "coefficient", values_to = "estimate") %>% 
  ggplot(aes(x = slice_id_oos, y = estimate)) +
  geom_line() +
  facet_wrap(~coefficient, ncol = 1, scales = "free_y") +
  labs(
    title = "Cross-sectional features model parameters",
    x = "Date",
    y = "Estimate"
  )

# plot time-series estimates
# ts_coefs_df %>% 
#   select(-slice_id) %>% 
#   pivot_longer(cols = -slice_id_oos, names_to = "coefficient", values_to = "estimate") %>% 
#   ggplot(aes(x = slice_id_oos, y = estimate)) +
#   geom_line() +
#   facet_wrap(~coefficient, ncol = 1, scales = "free_y") +
#   labs(
#     title = "Time-series features model parameters",
#     x = "Date",
#     y = "Estimate"
#   )

# join and fill using slice_id to designate when the model goes oos
exp_return_df <- model_df %>% 
  left_join(
    xs_coefs_df %>% left_join(ts_coefs_df, by = c("slice_id", "slice_id_oos")), 
    by = join_by(closest(date > slice_id_oos)), suffix = c("_factor", "_coef")
  ) %>% 
  na.omit() %>% 
  # forecast cross-sectional expected return as
  mutate(expected_xs_return = carry_decile_factor*carry_decile_coef + momo_decile_factor*momo_decile_coef + xs_intercept ) %>% 
  # mean expected xs return each day is zero
  # let total expected return be xs return + ts return - allows time series expected return to tilt weights
  mutate(expected_ts_return = case_when(breakout_factor >= 5.5 ~ 0.002, TRUE ~ 0)) %>% 
  ungroup() 

# long-short the xs expected return
# layer ts expected return on top
# position by expected return

# 1 in the numerator lets it get max 100% long due to breakout
max_ts_pos <- 0.5/trading_universe_size

strategy_df <- exp_return_df %>% 
  filter(date >= start_date) %>% 
  group_by(date) %>% 
  mutate(xs_position = expected_xs_return - mean(expected_xs_return, na.rm = TRUE)) %>% 
  # scale positions so that leverage is 1
  group_by(date) %>% 
  mutate(xs_position = if_else(xs_position == 0, 0, xs_position/sum(abs(xs_position)))) %>%
  # layer ts expected return prediction
  ungroup() %>% 
  mutate(ts_position = sign(expected_ts_return)) %>% 
  # constrain maximum delta added by time series prediction
  mutate(ts_position = if_else(ts_position >= 0, pmin(ts_position, max_ts_pos), pmax(ts_position, -max_ts_pos))) %>% 
  mutate(position = xs_position + ts_position) %>% 
  # strategy return
  mutate(strat_return = position*total_fwd_return_simple) %>% 
  # scale back to leverage 1
  group_by(date) %>% 
  mutate(position = if_else(position == 0, 0, position/sum(abs(position)))) 
  
returns_plot <- strategy_df %>%
  group_by(date) %>% 
  summarise(total_ret = sum(strat_return)) %>% 
  ggplot(aes(x = date, y = cumsum(log(1+total_ret)))) +
  geom_line() +
  labs(
    title = "Cumulative strategy return",
    y = "Cumulative return"
  )

weights_plot <- strategy_df %>%
  summarise(net_pos = sum(position)) %>% 
  ggplot(aes(x = date, y = net_pos)) +
  geom_line() +
  labs(
    x = "Date",
    y = "Net Weight"
  )

returns_plot / weights_plot + plot_layout(heights = c(2,1))

# tickers that were ever in the universe
universe_tickers <- features %>%
  filter(is_universe) %>%
  pull(ticker) %>%
  unique()

# universe_tickers

strategy_df <- exp_return_df %>% 
  filter(date >= start_date) %>% 
  group_by(date) %>% 
  mutate(expected_xs_return = expected_xs_return - mean(expected_xs_return, na.rm = TRUE)) %>% 
  ungroup() %>% 
  mutate(expected_return = expected_xs_return + expected_ts_return) %>% 
  select(date, ticker, expected_return) %>% 
  # join back onto df of prices for all tickers that were ever in the universe
  # so that we have prices before and after a ticker comes into or out of the universe
  # for backtesting purposes
  right_join(
    features %>% 
      filter(ticker %in% universe_tickers) %>% 
      select(date, ticker, close, total_fwd_return_simple, funding_rate), 
    by = c("date", "ticker")
  ) %>% 
  arrange(date, ticker) %>%
  filter(date >= start_date)

# get same length exp returns vector for each date
exp_returns_wide <- strategy_df %>% 
  select(ticker, date, expected_return) %>% 
  # will give NA where a ticker didn't exist
  pivot_wider(id_cols = date, names_from = ticker, values_from = expected_return)

# EWMA covariance estimate
# note definition of lambda in line with Risk Metrics
# ie higher values of lambda put less weight on the most recent returns and more weight on historical returns.
ewma_cov <- function(x, y, lambda, initialisation_wdw = 100) {
  # check that x and y are the same length and greater than initialisation_wdw
  stopifnot("x and y must be of equal length" = length(x) == length(y))
  if(length(x) <= initialisation_wdw) {
    ewma_cov <- rep(NA, length(x))
    return(ewma_cov)
  }

  # create initialisation window and estimation window
  init_x = x[1:initialisation_wdw]
  init_y = y[1:initialisation_wdw]

  num_obs <- length(x)

  # initial covariance and mean return estimates
  old_cov <- cov(init_x, init_y)
  old_x <- mean(init_x)
  old_y <- mean(init_y)

  # preallocate output vector
  ewma_cov <- vector(mode = "numeric", length = num_obs)

  # pad with NA for initialisation window
  ewma_cov[1:initialisation_wdw] <- NA

  # covariance estimate
  for(i in c((initialisation_wdw+1):num_obs)) {
    ewma_cov[i] <- lambda*old_cov + (1 - lambda)*(old_x * old_y)
    old_cov <- ewma_cov[i]
    old_x <- x[i]
    old_y <- y[i]
  }
  ewma_cov
}

cov_lambda <- 0.995
wdw <- ew_cov_init_wdw

returns <- features %>% 
  filter(ticker %in% universe_tickers) %>% 
  select(ticker, date, total_return_simple)

# long dataframe of pairwise EW covariances
ewma_covs <- returns %>%
  full_join(returns, by = "date") %>%
  na.omit() %>%
  ungroup() %>%
  # get all combinations (tickers) and remove duplicate combos (eg BTC-ETH, ETH-BTC)
  mutate(tickers = ifelse(ticker.x < ticker.y, glue("{ticker.x}, {ticker.y}"), glue("{ticker.y}, {ticker.x}"))) %>%
  distinct(date, tickers, .keep_all = TRUE) %>%
  # calculate rolling pairwise ewma correlations
  group_by(tickers) %>%
  arrange(date, .by_group = TRUE) %>%
  mutate(ewma_cov = ewma_cov(total_return_simple.x, total_return_simple.y, lambda = cov_lambda, initialisation_wdw = ew_cov_init_wdw)) %>%
  select(date, tickers, ewma_cov, ticker.x, ticker.y) %>%
  na.omit() %>% 
  ungroup()

tail(ewma_covs)

# plot pairwise rolling covariances
ewma_covs %>%
  filter(str_starts(tickers, "BTCUSDT, S")) %>%
  ggplot(aes(x = date, y = ewma_cov, colour = tickers)) +
    geom_line(size=0.8) +
    labs(
      title = "Rolling Window Covariance Estimates",
      colour = "Tickers"
    )

# functions for wrangling covariance matrix
library(Matrix) 
library(matrixcalc)

get_today_covmat <- function(ewma_covs_long, today_date) {
    today_cov_df <- ewma_covs_long %>%
        ungroup() %>% 
        filter(date == today_date) %>% 
        select(ticker.x, ticker.y, ewma_cov) %>%
        pivot_wider(names_from = ticker.y, values_from = ewma_cov) %>%  # , names_prefix = "cov_"
        column_to_rownames(var = "ticker.x")
    
    # check rownames and colnames match
    stopifnot("error making today's cov matrix" = all.equal(rownames(today_cov_df), colnames(today_cov_df)))
    today_covmat <- as.matrix(today_cov_df)
    today_covmat[lower.tri(today_covmat)] <- t(today_covmat)[lower.tri(today_covmat)]
    
    # return nearest positive semi-definite matrix
    today_covmat <- nearPD(today_covmat)
    return(today_covmat$mat)
}

#' shrinks covs towards zero, reflecting uncertainty in cov estimates. may increases variances.
#' @param shrinkage_intensity: a value between 0 and 1. 0 means no shrinkage (use the original covmat), 1 means full shrinkage (use only the target_matrix)
shrink_covmat <- function(covmat, shrinkage_intensity = 0.9) {
    # average variance from the diagonal 
    average_variance <- mean(diag(covmat))
    # shrinkage target matrix (identity matrix scaled by average variance)
    target_matrix <- diag(average_variance, nrow = ncol(covmat))
    # shrunk covariance matrix
    covmat_shrunk <- (1 - shrinkage_intensity) * covmat + shrinkage_intensity * target_matrix

    return(covmat_shrunk)
}

#' extend covmat to cover all tickers in the tradeable universe
#' put average variance on the diagonal, zeroes elsewhere
#' intent is that optimisation will explicitly constrain weights for assets not yet in the universe to 0
extend_covmat <- function(covmat, all_tickers) {
    current_tickers <- colnames(covmat) 
    
    full_covmat <- matrix(0, nrow = length(all_tickers), ncol = length(all_tickers),
         dimnames = list(all_tickers, all_tickers))
    
    current_tickers_ordered <- intersect(all_tickers, rownames(covmat))

    # ensure covmats are ordered correctly
    covmat_reordered <- covmat[current_tickers_ordered, current_tickers_ordered] %>% 
      # ensure is Matrix type for compatibility
      as.matrix()
    full_covmat <- full_covmat[all_tickers, all_tickers]

    # replace in full_covmat
    full_covmat[current_tickers_ordered, current_tickers_ordered] <- covmat_reordered

    # replace zeros on diagonal with average variance
    zero_diag_indices <- which(diag(full_covmat) == 0)
    average_variance <- mean(diag(covmat))
    diag(full_covmat)[zero_diag_indices] <- average_variance
    
    return(full_covmat)
}

wrangle_covmat <- function(ewma_covs_long, date, tickers, shrinkage_intensity = 0.9) {
    today_covmat <- get_today_covmat(ewma_covs_long, date) %>% 
      shrink_covmat(shrinkage_intensity = shrinkage_intensity) %>% 
      extend_covmat(all_tickers = tickers)

    # check
    if(is.positive.semi.definite(today_covmat, tol=1e-8)){
        if(isSymmetric(today_covmat)) {
            glue("covmat crated for {date} not symmetric")
        }
    } else {
        glue("covmat crated for {date} not PSD")
    }    
    
    return(today_covmat)
}

# precompute covariances matrixes 
simulation_start_date <- start_date + days(100)  # allow for first in-sample period - otherwise we run lots of simulation loops for zero weights
exp_returns_sim_df <- exp_returns_wide %>% filter(date >= simulation_start_date)

dates <- exp_returns_sim_df$date
covmat_list <- map(dates, ~wrangle_covmat(ewma_covs, .x, tickers = universe_tickers, shrinkage_intensity = 0.9))

# define our mvo function
mvo_with_costs <- function(expected_returns, current_weights, na_idxs = c(), costs, covmat, lambda = 1, tau=1) {
  # define our weights vector as a CVXR::Variable
  weights <- Variable(length(expected_returns))
  # define an alpha term as the weighted sum of expected returns
    # again, express using linear algebra
  alpha_term <- (t(weights) %*% expected_returns)
  # define a costs term. depends on:
    # cost of trading - needs to be expressed such that it scales with expected returns
      # calculate as elementwise cost * absolute value of weights - current_weights
      # use CVXR::multiply and CVXR::abs
    # absolute distance of current_weights to our weights variable
      # the more our target weights differ from current weights, the more it costs to trade
      # this is a decent representation of fixed percentage costs, but doesn't capture minimum commissions
    # sum_entries is a CVXR function for summing the elements of a vector
  costs_term <- sum_entries(multiply(costs, abs(weights - current_weights))) # elementwise abs, multiply
  # define a risk term as w*Sigma*w
    # quad_form is a CVXR function for doing w*Sigma*w
  risk_term <- quad_form(weights, covmat)
  # define our objective
    # maximise our alpha less our risk term multiplied by some factor, lambda, less our costs term multiplied by tau
  objective <- Maximize(alpha_term - lambda*risk_term - tau*costs_term)
  # apply our no leverage constraint
  constraints <- list(cvxr_norm(weights, 1) <= 1, weights[na_idxs] == 0)
  # specify the problem
  problem <- Problem(objective, constraints)
  # solve
  result <- solve(problem)

  # return the values of the variable we solved for
  result$getValue(weights)
}

library(doParallel)
library(foreach)

# register parallel backend
num_cores <- detectCores() - 1  # leave one core free
cl <- makeCluster(num_cores)
registerDoParallel(cl)

# use clusterEvalQ to define functions and objects available to each worker
clusterExport(
  cl, 
  varlist = c("get_today_covmat", "shrink_covmat", "extend_covmat", "wrangle_covmat", "mvo_with_costs", "universe_tickers", "exp_returns_sim_df", "covmat_list")
)

# load libraries in each worker
clusterEvalQ(cl, {
  library(CVXR); library(Matrix); library(matrixcalc); library(dplyr); library(tidyr); library(glue)
})

num_days <- nrow(exp_returns_wide)
lambdas <- c(0.1, 0.3, 1, 3)
taus <- c(0.003, 0.1, 0.3, 1, 3)
date_idxs_to_do <- c(1:500)
costs <- 0.15/100

# custom combine function for flattening results of foreach to one big list
# ensures that each result from the inner loop is appended to the main list, rather than nesting it further 
# results in a flat structure where each element is a list corresponding to an iteration over lambda and tau without deep nesting
# requres .init arg of foreach to be list() - initialises an empty list for storing the results
flatten_lists <- function(x, y) {
  c(x, y)
}

results <- foreach(lambda = lambdas, .combine = 'flatten_lists', .init = list()) %:%   # 
  foreach(tau = taus, .combine = 'flatten_lists', .init = list()) %dopar% {
   
   weights <- list()
   errors <- list()
   w0 <- rep(0, length(universe_tickers))
   for( i in date_idxs_to_do ) {
      # today's date
      d = exp_returns_sim_df$date[i]

      # cov estimate
      today_covmat <- covmat_list[[i]]

      # check cov matrix and exp returns vector are ticker-aligned
      # all.equal(exp_returns_sim_df %>% filter(date == d) %>% pivot_longer(-date, names_to = "ticker", values_to = "expected_return")%>% arrange(match(ticker, universe_tickers)) %>% pull(ticker), colnames(today_covmat))
      
      # get row of expected returns as a vector in the same order as the columns of the covariance matrix
      # contains NA at this point
      exp_rets <- exp_returns_sim_df %>% filter(date == d) %>% pivot_longer(-date, names_to = "ticker", values_to = "expected_return")%>% arrange(match(ticker, universe_tickers)) %>% pull(expected_return)
      
      # get na indexes
      # we'll explicitly constrain these to get zero weight
      na_idxs <- which(is.na(exp_rets))
      # convert NA to zero
      exp_rets[is.na(exp_rets)] <- 0
      
      # initialise w in case we get a solver error (can retain w0 in these cases)
      w <- w0
      
      # solve for next period's weights
      tryCatch({
         w <- mvo_with_costs(expected_returns = exp_rets, current_weights = w0, na_idxs = na_idxs, costs = costs, covmat = today_covmat, lambda = lambda, tau = tau)
      }, error = function(e) {
            # in case of solver error, log and retain existing weights
            errors[[i]] <- e
         }
      )
      
      # store weights
      # w is a one-column matrix
      weights[[i]] <- w
      w0 <- w 
    }
    # return a list with weights, parameters, and any errors
    list(list(weights = weights, lambda = lambda, tau = tau, errors = errors))
}
  
stopCluster(cl)

# save to disk
saveRDS(results, file = "lambda_tau_search_results_1_500.rds")

# read results from disk
# results <- readRDS("lambda_tau_search_results_1_500.rds")

# function for converting a day's weights vector to a dataframe
weights_mat_to_df <- function(x, d, lambda, tau, tickers) {
  x %>% 
    as_tibble(.name_repair = ~paste0("lambda_",lambda,"_tau_",tau)) %>% 
    mutate(date = d) %>%
    mutate(ticker = tickers)
}

# make one dataframe of weights for each lambda-tau combination
weights_dfs <- list()
for(i in seq_along(results)) {
  this_entry <- results[[i]]
  these_weights <- this_entry$weights
  this_lambda <- this_entry$lambda
  this_tau <- this_entry$tau

  weights_df <- purrr::map2(
    these_weights, 
    exp_returns_wide$date[date_idxs_to_do], 
    ~weights_mat_to_df(.x, .y, lambda = this_lambda, tau = this_tau, tickers = universe_tickers)
  ) %>% 
    bind_rows()
  
  weights_dfs[[i]] <- weights_df
}

# combine dataframes for each lmbda-tau combination into one big dataframe
all_weights <- weights_dfs %>% purrr::reduce(left_join, by = c("date", "ticker")) %>% 
  relocate(c(date, ticker))

head(all_weights)

all_weights %>% 
  pivot_longer(cols = c(-date, -ticker), names_to = "parameterisation", values_to = "weight") %>% 
  filter(ticker == "BTCUSDT") %>% 
  ggplot(aes(x = date, y = weight, colour = parameterisation)) +
  geom_line() +
  facet_wrap(~parameterisation) +
  theme(legend.position = "none") +
  labs(
    title = "BTC weights for various paramterisations"
  )

all_weights %>% 
  pivot_longer(cols = c(-date, -ticker), names_to = "parameterisation", values_to = "weight") %>% 
  group_by(parameterisation, date) %>% 
  summarise(total_weight = sum(abs(weight)), net_weight = (sum(weight)), .groups = "drop") %>%   
  pivot_longer(cols = c(total_weight, net_weight), names_to = "total_net", values_to = "weight")  %>% 
  ggplot(aes(x = date, y = weight, colour = total_net)) + 
  geom_line() + 
  facet_wrap(~parameterisation) +
  labs(
    title = "Total and net portfolio weight for various paramterisations"
  )

backtest_df <- all_weights %>% 
  pivot_longer(cols = c(-date, -ticker), names_to = "parameterisation", values_to = "weight") %>%  
  left_join(strategy_df, by = c("ticker", "date"))

head(backtest_df)

params <- unique(backtest_df$parameterisation)
params

# simulation for each parameterisation

# no trade buffer (trades controlled by tau parameter and are embedded in weights)
equity_curves <- vector("list", length = length(params))
results_dfs <- vector("list", length = length(params))
i <- 1
for( p in params ) {
    # get weights as a wide matrix 
    # note that date column will get converted to unix timestamp
    backtest_weights <- backtest_df %>%
        dplyr::filter(parameterisation == p) %>% 
        pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, weight)) %>%  # pivot wider guarantees prices and theo_weight are date aligned
        select(date, starts_with("weight")) %>%
        mutate(date = as.numeric(date)) %>% 
        data.matrix()

    # NA weights should be zero
    backtest_weights[is.na(backtest_weights)] <- 0

    # get prices as a wide matrix
    # note that date column will get converted to unix timestamp
    backtest_prices <- backtest_df %>% 
        dplyr::filter(parameterisation == p) %>% 
        pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, weight)) %>%  # pivot wider guarantees prices and theo_weight are date aligned
        select(date, starts_with("close_")) %>%
        mutate(date = as.numeric(date)) %>% 
        data.matrix()

    # get funding as a wide matrix
    # note that date column will get converted to unix timestamp
    backtest_funding <- backtest_df %>%
        dplyr::filter(parameterisation == p) %>% 
        pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, funding_rate)) %>%  # pivot wider guarantees prices and funding_returns_simple are date aligned
        select(date, starts_with("funding_rate_")) %>%
        mutate(date = as.numeric(date)) %>% 
        data.matrix()

    # simulation
    results_df <- fixed_commission_backtest_with_funding(
        prices = backtest_prices,
        target_weights = backtest_weights,
        funding_rates = backtest_funding,
        trade_buffer = 0.,
        initial_cash = 10000,
        margin = 0.05,
        commission_pct = 0.0015,
        capitalise_profits = FALSE
        ) %>%
        mutate(ticker = str_remove(ticker, "close_")) %>%
        # remove coins we don't trade from results
        drop_na(Value)

    margin <- results_df %>%
        group_by(Date) %>%
        summarise(Margin = sum(Margin, na.rm = TRUE))

    cash_balance <- results_df %>%
        filter(ticker == "Cash") %>%
        select(Date, Value) %>%
        rename("Cash" = Value)

    equity <- cash_balance %>%
        left_join(margin, by = "Date") %>% 
        mutate(Equity = Cash + Margin) %>% 
        select(Date, Equity) %>% 
        rename(!!p := Equity)
    
    equity_curves[[i]] <- equity
    results_dfs[[i]] <- results_df

    i <- i + 1
}

equity_curves %>% 
  purrr::discard(is.null) %>% 
  reduce(left_join, by = "Date") %>% 
  pivot_longer(-Date, names_to = "param", values_to = "equity") %>%
  ggplot(aes(x = Date, y = equity, colour = param)) +
  geom_line() +
  facet_wrap(~param) +
  theme(legend.position = "none") +
  labs(
    title = "After-costs equity curves for 500-day simulations"
  )

initial_cash <- 10000
i <- 1
dfs <- list()
for (r in results_dfs) {
  # get lambda and tau
  these_params <- str_extract_all(params[i],"[0-9]+\\.*[0-9]*", simplify = TRUE)
  this_lambda <- these_params[1, 1]
  this_tau <- these_params[1, 2]

  # average daily turnover
  average_turnover <- r %>%
    filter(ticker != "Cash") %>% 
    group_by(Date) %>% 
    summarise(turnover = 100*sum(abs(TradeValue))/initial_cash, .groups = "drop") %>% 
    summarise(ave_turnover = mean(turnover)) %>% 
    pull(ave_turnover)

  # total return
  equity <- equity_curves[[i]][, 2, drop = TRUE]

  fin_eq <- equity %>%
    tail(1) 

  init_eq <- equity%>%
    head(1) 

  total_return <- (fin_eq/init_eq - 1) * 100
  days <- length(equity)
  ann_return <- total_return * 365/days
  
  # sharpe
  returns <- na.omit(equity/lag(equity) - 1)
  sharpe <- sqrt(365)*mean(returns)/sd(returns)

  # store results
  df <- data.frame(
    lambda = this_lambda,
    tau = this_tau,
    sharpe = sharpe,
    total_return = total_return,
    average_turnover = average_turnover
  )
  dfs[[i]] <- df

  i <- i + 1
}

metrics <- dfs %>% bind_rows()
metrics %>% arrange(desc(sharpe))

lambda <- 1
tau <- 1
costs <- 0.15/100
   
weights <- list()
errors <- list()
w0 <- rep(0, length(universe_tickers))
for( i in c(1:length(exp_returns_sim_df$date)) ) {
    # today's date
    d = exp_returns_sim_df$date[i]

    # cov estimate
    today_covmat <- covmat_list[[i]]

    # check cov matrix and exp returns vector are ticker-aligned
    # all.equal(exp_returns_sim_df %>% filter(date == d) %>% pivot_longer(-date, names_to = "ticker", values_to = "expected_return")%>% arrange(match(ticker, universe_tickers)) %>% pull(ticker), colnames(today_covmat))
    
    # get row of expected returns as a vector in the same order as the columns of the covariance matrix
    # contains NA at this point
    exp_rets <- exp_returns_sim_df %>% filter(date == d) %>% pivot_longer(-date, names_to = "ticker", values_to = "expected_return")%>% arrange(match(ticker, universe_tickers)) %>% pull(expected_return)
    # get na indexes
    # we'll explicitly constrain these to get zero weight
    na_idxs <- which(is.na(exp_rets))
    # convert NA to zero
    exp_rets[is.na(exp_rets)] <- 0
    # initialise w in case we get a solver error (can retain w0 in these cases)
    w <- w0
    tryCatch({
        w <- mvo_with_costs(expected_returns = exp_rets, current_weights = w0, na_idxs = na_idxs, costs = costs, covmat = today_covmat, lambda = lambda, tau = tau)
    }, error = function(e) {
        # in case of solver error, log and retain existing weights
        errors[[i]] <- e
        }
    )
    
    # w is a one-column matrix
    weights[[i]] <- w
    w0 <- w 
}

results <- list(weights = weights, lambda = lambda, tau = tau, errors = errors)

saveRDS(results, file = "lambda_1_tau_1.rds")

# wrangle weights into a dataframe
weights_df <- purrr::map2(
  weights, 
  exp_returns_sim_df$date, 
  ~weights_mat_to_df(.x, .y, lambda = lambda, tau = tau, tickers = universe_tickers)
) %>% 
  bind_rows() %>% 
  rename("weight" = lambda_1_tau_1)

head(weights_df)

# join weights onto prices
backtest_df <- weights_df %>% 
  left_join(strategy_df, by = c("ticker", "date"))

head(backtest_df)

# simulate

# get weights as a wide matrix 
backtest_weights <- backtest_df %>%
  pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, weight)) %>%  # pivot wider guarantees prices and theo_weight are date aligned
  select(date, starts_with("weight")) %>%
  mutate(date = as.numeric(date)) %>% 
  data.matrix()

# NA weights should be zero
backtest_weights[is.na(backtest_weights)] <- 0

# get prices as a wide matrix
backtest_prices <- backtest_df %>% 
  pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, weight)) %>%  # pivot wider guarantees prices and theo_weight are date aligned
  select(date, starts_with("close_")) %>%
  mutate(date = as.numeric(date)) %>% 
  data.matrix()

# get funding as a wide matrix
backtest_funding <- backtest_df %>%
  pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, funding_rate)) %>%  # pivot wider guarantees prices and funding_returns_simple are date aligned
  select(date, starts_with("funding_rate_")) %>%
  mutate(date = as.numeric(date)) %>% 
  data.matrix()

# simulation
results_df <- fixed_commission_backtest_with_funding(
  prices = backtest_prices,
  target_weights = backtest_weights,
  funding_rates = backtest_funding,
  trade_buffer = 0.,
  initial_cash = 10000,
  margin = 0.05,
  commission_pct = 0.0015,
  capitalise_profits = FALSE
) %>%
mutate(ticker = str_remove(ticker, "close_")) %>%
# remove coins we don't trade from results
drop_na(Value)

margin <- results_df %>%
  group_by(Date) %>%
  summarise(Margin = sum(Margin, na.rm = TRUE))

cash_balance <- results_df %>%
  filter(ticker == "Cash") %>%
  select(Date, Value) %>%
  rename("Cash" = Value)

equity <- cash_balance %>%
  left_join(margin, by = "Date") %>% 
  mutate(Equity = Cash + Margin) %>% 
  select(Date, Equity) 

# plot results
equity_plot <- equity %>% 
  ggplot(aes(x = Date, y = Equity)) + 
  geom_line() +
  labs(
    title = "Simulation of carry, momentum & breakout factors with lambda 1, tau 1",
    y = "Equity, $"
  )

weights_plot <- backtest_df %>%
  group_by(date) %>% 
  summarise(
    total_weight = sum(abs(weight)),
    net_weight = sum(weight)
  ) %>% 
  pivot_longer(-date, names_to = "type", values_to = "weight") %>% 
  ggplot(aes(x = date, y = weight, colour = type)) +
  geom_line() +
  labs(
    title = "Total and Net Weight",
    y = "Weight",
    colour = ""
  ) +
  theme(legend.position = "bottom")

turnover_plot <- results_df  %>% 
  filter(ticker != "Cash") %>% 
  group_by(Date) %>% 
  summarise(Turnover = 100*sum(abs(TradeValue))/initial_cash) %>% 
  ggplot(aes(x = Date, y = Turnover)) +
    geom_line() +
    geom_point() +
    labs(
      title = "Turnover as % of trading capital",
      y = "Turnover, %"
    )

options(repr.plot.width = 14, repr.plot.height=12)
equity_plot / weights_plot / turnover_plot + plot_layout(heights = c(2, 1, 1))
options(repr.plot.width = 14, repr.plot.height=7)

# define our mvo function
mvo_with_costs_constrained <- function(expected_returns, current_weights, na_idxs = c(), costs, covmat, lambda = 1, tau = 1) {
  # define our weights vector as a CVXR::Variable
  weights <- Variable(length(expected_returns))
  # define an alpha term as the weighted sum of expected returns
    # again, express using linear algebra
  alpha_term <- (t(weights) %*% expected_returns)
  # define a costs term. depends on:
    # cost of trading - needs to be expressed such that it scales with expected returns
      # calculate as elementwise cost * absolute value of weights - current_weights
      # use CVXR::multiply and CVXR::abs
    # absolute distance of current_weights to our weights variable
      # the more our target weights differ from current weights, the more it costs to trade
      # this is a decent representation of fixed percentage costs, but doesn't capture minimum commissions
    # sum_entries is a CVXR function for summing the elements of a vector
  costs_term <- sum_entries(multiply(costs, abs(weights - current_weights))) # elementwise abs, multiply
  # define a risk term as w*Sigma*w
    # quad_form is a CVXR function for doing w*Sigma*w
  risk_term <- quad_form(weights, covmat)
  # define our objective
    # maximise our alpha less our risk term multiplied by some factor, lambda, less our costs term multiplied by tau
  objective <- Maximize(alpha_term - lambda*risk_term - tau*costs_term)
  # constrain to a maximum net weight +/- 0.5, maximum leverage 1
  constraints <- list(cvxr_norm(weights, 1) <= 1, abs(sum_entries(weights)) <= 0.5, weights[na_idxs] == 0)
  # specify the problem
  problem <- Problem(objective, constraints)
  # solve
  result <- solve(problem)

  # return the values of the variable we solved for
  result$getValue(weights)
}

lambda <- 1
tau <- 1
weights <- list()
errors <- list()
w0 <- rep(0, length(universe_tickers))
for( i in c(1:length(exp_returns_sim_df$date)) ) {
    # today's date
    d = exp_returns_sim_df$date[i]

    # cov estimate
    today_covmat <- covmat_list[[i]]

    # check cov matrix and exp returns vector are ticker-aligned
    # all.equal(exp_returns_sim_df %>% filter(date == d) %>% pivot_longer(-date, names_to = "ticker", values_to = "expected_return")%>% arrange(match(ticker, universe_tickers)) %>% pull(ticker), colnames(today_covmat))
    
    # get row of expected returns as a vector in the same order as the columns of the covariance matrix
    # contains NA at this point
    exp_rets <- exp_returns_sim_df %>% filter(date == d) %>% pivot_longer(-date, names_to = "ticker", values_to = "expected_return")%>% arrange(match(ticker, universe_tickers)) %>% pull(expected_return)
    # get na indexes
    # we'll explicitly constrain these to get zero weight
    na_idxs <- which(is.na(exp_rets))
    # convert NA to zero
    exp_rets[is.na(exp_rets)] <- 0
    # initialise w in case we get a solver error (can retain w0 in these cases)
    w <- w0
    tryCatch({
        w <- mvo_with_costs_constrained(expected_returns = exp_rets, current_weights = w0, na_idxs = na_idxs, costs = costs, covmat = today_covmat, lambda = lambda, tau = tau)
    }, error = function(e) {
        # in case of solver error, log and retain existing weights
        errors[[i]] <- e
        }
    )
    
    # w is a one-column matrix
    weights[[i]] <- w
    w0 <- w 
}

saveRDS(results, file = "lambda_1_tau_1_net_weight_constrained.rds")

# wrangle weights into a dataframe
weights_df <- purrr::map2(
  weights, 
  exp_returns_sim_df$date, 
  ~weights_mat_to_df(.x, .y, lambda = lambda, tau = tau, tickers = universe_tickers)
) %>% 
  bind_rows() %>% 
  rename("weight" = lambda_1_tau_1)

# join weights onto prices
backtest_df <- weights_df %>% 
  left_join(strategy_df, by = c("ticker", "date"))

head(backtest_df)

# simulate

# get weights as a wide matrix 
backtest_weights <- backtest_df %>%
  pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, weight)) %>%  # pivot wider guarantees prices and theo_weight are date aligned
  select(date, starts_with("weight")) %>%
  mutate(date = as.numeric(date)) %>% 
  data.matrix()

# NA weights should be zero
backtest_weights[is.na(backtest_weights)] <- 0

# get prices as a wide matrix
backtest_prices <- backtest_df %>% 
  pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, weight)) %>%  # pivot wider guarantees prices and theo_weight are date aligned
  select(date, starts_with("close_")) %>%
  mutate(date = as.numeric(date)) %>% 
  data.matrix()

# get funding as a wide matrix
backtest_funding <- backtest_df %>%
  pivot_wider(id_cols = date, names_from = ticker, values_from = c(close, funding_rate)) %>%  # pivot wider guarantees prices and funding_returns_simple are date aligned
  select(date, starts_with("funding_rate_")) %>%
  mutate(date = as.numeric(date)) %>% 
  data.matrix()

# simulation
results_df <- fixed_commission_backtest_with_funding(
  prices = backtest_prices,
  target_weights = backtest_weights,
  funding_rates = backtest_funding,
  trade_buffer = 0.,
  initial_cash = 10000,
  margin = 0.05,
  commission_pct = 0.0015,
  capitalise_profits = FALSE
) %>%
mutate(ticker = str_remove(ticker, "close_")) %>%
# remove coins we don't trade from results
drop_na(Value)

margin <- results_df %>%
  group_by(Date) %>%
  summarise(Margin = sum(Margin, na.rm = TRUE))

cash_balance <- results_df %>%
  filter(ticker == "Cash") %>%
  select(Date, Value) %>%
  rename("Cash" = Value)

equity <- cash_balance %>%
  left_join(margin, by = "Date") %>% 
  mutate(Equity = Cash + Margin) %>% 
  select(Date, Equity) 

# plot results
equity_plot <- equity %>% 
  ggplot(aes(x = Date, y = Equity)) + 
  geom_line() +
  labs(
    title = "Simulation of carry, momentum & breakout factors with lambda 1, tau 1",
    subtitle = "Portfolio net position constrained to +/- 0.5",
    y = "Equity, $"
  )

weights_plot <- backtest_df %>%
  group_by(date) %>% 
  summarise(
    total_weight = sum(abs(weight)),
    net_weight = sum(weight)
  ) %>% 
  pivot_longer(-date, names_to = "type", values_to = "weight") %>% 
  ggplot(aes(x = date, y = weight, colour = type)) +
  geom_line() +
  labs(
    title = "Total and Net Weight",
    y = "Weight",
    colour = ""
  ) +
  theme(legend.position = "bottom")

turnover_plot <- results_df  %>% 
  filter(ticker != "Cash") %>% 
  group_by(Date) %>% 
  summarise(Turnover = 100*sum(abs(TradeValue))/initial_cash) %>% 
  ggplot(aes(x = Date, y = Turnover)) +
    geom_line() +
    geom_point() +
    labs(
      title = "Turnover as % of trading capital",
      y = "Turnover, %"
    )

options(repr.plot.width = 14, repr.plot.height=12)
equity_plot / weights_plot / turnover_plot + plot_layout(heights = c(2, 1, 1))
options(repr.plot.width = 14, repr.plot.height=7)