How Much Damage Can I Do Turbo-Punting Shitcoins?

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

# load and install packages
install.packages("tibbletime")  # Rcpp is a dependency
install.packages("patchwork")
install.packages("tidyverse")

library(tidyverse)
library(tibbletime)
library(patchwork)
library(Rcpp)

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

# hourly GBM simulator with autocorrelation and jump diffusion

# C++ function for doing the autogregressive and jump diffusion bit
cppFunction("
NumericMatrix generateAR1WithJumps(int nsim, int t, double phi, double lambda, double mu_j, double sigma_j) {
  NumericMatrix epsilon(t-1, nsim);
  NumericVector eta(nsim);
  NumericVector jumpSize(nsim);

  // Initial values from normal distribution
  for(int j = 0; j < nsim; ++j) {
    epsilon(0, j) = R::rnorm(0, 1);
  }

  // AR(1) process with jumps
  for(int i = 1; i < t-1; ++i) {
    for(int j = 0; j < nsim; ++j) {
      eta[j] = R::rnorm(0, 1);
      epsilon(i, j) = phi * epsilon(i-1, j) + eta[j];
      
      // Generate jumps
      if(R::runif(0, 1) < lambda) {
        jumpSize[j] = R::rnorm(mu_j, sigma_j);
        epsilon(i, j) += jumpSize[j];
      }
    }
  }
  
  return epsilon;
}
")

# gbm generator function
#' @param phi: autocorrelation
#' @param lambda: jump intensity
#' @param mu_j @param sigma_j: define the distribution of the jump sizes 
#' other params as per standard GBM model
#' NOTE simulates hourly prices using annualised volatility and returns
gbm_autocor_jumps <- function(nsim = 100, t = 25, mu = 0, sigma = 0.1, S0 = 100, dt = 1./(365*24), phi = 0.5, lambda = 0.01, mu_j = 0, sigma_j = 0.1) {
  
  # generate epsilon with AR(1) and jumps
  epsilon <- generateAR1WithJumps(nsim, t, phi, lambda, mu_j, sigma_j)
  
  # get GBM and convert to price paths
  gbm <- exp((mu - sigma * sigma / 2) * dt + sigma * epsilon * sqrt(dt))
  gbm <- apply(rbind(rep(S0, nsim), gbm), 2, cumprod)

  return(gbm)
}

# make some random price series with shitcoin-like parameters
nsim <- 100
t <- 100*24  # 100 days
mu <- -1
sigma <- 1.5
S0 <- 100
phi <- 0.1

set.seed(503)  # so we can reproduce results
gbm <- gbm_autocor_jumps(nsim, t, mu, sigma, S0, phi = phi, lambda = 0.2, mu_j = 0, sigma_j = 1)

gbm_df <- as.data.frame(gbm) %>%
  mutate(ix = 1:nrow(gbm)) %>%
  pivot_longer(-ix, names_to = 'sim', values_to = 'price')

gbm_df %>%
  ggplot(aes(x=ix, y=price, color=sim)) +
  geom_line() +
  theme(legend.position = 'none') +
  labs(x = "hour")

# ape new highs strategy

# calculate positions
hold_period <- 5*24  # in hours

positions_df <- gbm_df %>%
  group_by(sim) %>%
  mutate(all_time_high = cummax(price)) %>%
  mutate(idx_all_time_high = match(all_time_high, unique(price))) %>%
  mutate(lagged_idx_all_time_high = dplyr::lag(idx_all_time_high)) %>%
  mutate(periods_since_high = ix - lagged_idx_all_time_high) %>%
  mutate(position = case_when(
    between(periods_since_high, 1, hold_period) & lagged_idx_all_time_high > 1 ~ 1, TRUE ~ 0
  ))

tail(positions_df)

# example plot of when we're in a position for a given series
positions_df %>%
  mutate(position_plot = case_when(position == 1 ~ price, TRUE ~ NA_real_)) %>%
  filter(sim == "V2") %>%
  ggplot(aes(x = ix, y = price)) +
  geom_line() +
  geom_point(aes(x = ix, y = position_plot), colour = "blue") +
  labs(title = "Positions in a single coin", x = "hour")

# simple returns
returns_df <- positions_df %>%
  group_by(sim) %>%
  # assume we enter a position at the hourly price at which we make a new high 
  mutate(simple_return = (price - dplyr::lag(price))/dplyr::lag(price)) %>%
  mutate(strategy_simple_return = position*simple_return) %>%
  na.omit() %>%
  ungroup()

# will need to do in a for loop since actual leverage will be dependent on drawdown
# calculate portfolio pnl for each day 
# adjust leverage for next day if needs be (if we're at rebal frequency or positions change)

# need a function for defining how we modify our leverage
set_leverage <- function(max_leverage, current_drawdown, max_acceptable_drawdown=0.9) {
  # decrease leverage linearly with drawdown
  # in practive, you might want to decrease it faster - maybe some non-linear function
  leverage <- max(min(max_leverage - max_leverage/max_acceptable_drawdown*current_drawdown, max_leverage), 0)
  return(leverage)
}

# this bit is going to be slow
# could optimise by replacing dataframes with matrixes if we want to do lots
starting_leverage <- 5
rebalance_frequency <- 24  # in hours
starting_leverage <- 5
all_time_equity_high <- 1 # starting ATH
portfolio_results <- list()

for(i in min(returns_df$ix):max(returns_df$ix)) {
  if(i == min(returns_df$ix)) {
    # start by sizing all positions using starting leverage
    leverage <- starting_leverage
    num_positions <- returns_df %>%
      filter(ix == i) %>%
      # calculate current positions for target leverage
      summarise(num_positions = sum(position)) %>%
      pull(num_positions)
  
    size_per_position <- starting_leverage/num_positions

    # adjust positions for all future rows (these will be updated in subsequent steps)  
    returns_df <- returns_df %>%
      mutate(leveraged_position = position*size_per_position) %>%
      mutate(strategy_simple_return = leveraged_position*simple_return)
  }

  # calculate current positions for target leverage
  current_num_positions <- returns_df %>%
      filter(ix == i) %>%
      summarise(num_positions = sum(position)) %>%
      pull(num_positions)
  
  # assume we rebalance at the start of the period if we change positions or hit our rebalance frequency    
  if(i %% rebalance_frequency == 0 | current_num_positions != num_positions) {
    leverage <- set_leverage(starting_leverage, current_drawdown)

    if(current_num_positions > 0)
      size_per_position <- leverage/current_num_positions

    # adjust positions for all future rows (these will be updated in subsequent steps)  
    returns_df <- returns_df %>%
      mutate(leveraged_position = position*size_per_position) %>%
      mutate(strategy_simple_return = leveraged_position*simple_return)
  }

  num_positions <- current_num_positions
    
  # calculate leveraged portfolio returns
  leveraged_port_returns <- returns_df %>%
    filter(ix == i) %>%
    # leveraged portfolio returns will be the equal-weighted strategy returns scaled to the target leverage
    group_by(ix) %>%
    summarise(leveraged_returns = sum(strategy_simple_return, na.rm = TRUE)) %>%
    mutate(leverage = leverage)

  # store results
  portfolio_results[[i-1]] <- leveraged_port_returns
    
  # keep a running portfolio equity and drawdown
  current_port_returns <- portfolio_results %>%
    bind_rows() %>%
    # just use simple returns
    mutate(cum_port_returns = cumprod(1 + leveraged_returns)) %>%
    tail(1) %>%
    pull(cum_port_returns)
  
  all_time_equity_high <- max(current_port_returns, all_time_equity_high)
  current_drawdown <- all_time_equity_high - current_port_returns
}

# plot results
results_df <- bind_rows(portfolio_results) %>%
  mutate(cum_port_returns = cumprod(1 + leveraged_returns)) 

# get all-time-high to plot drawdown target
drawdown_target <- max(results_df$cum_port_returns)*0.1

returns_plot <- results_df %>%
  ggplot(aes(x = ix, y = cum_port_returns)) +
  geom_line() +
  geom_hline(yintercept = drawdown_target, colour = "red", size = 1, linetype = "dashed") +
  labs(
    x = "Hour",
    y = "Portfolio return",
    title = glue::glue("Portfolio returns and leverage, rebal frequency {rebalance_frequency} hours"),
    subtitle = "Drawdown target shown in red (90% off ATH)"
  )

leverage_plot <- results_df %>%
  ggplot(aes(x = ix, y = leverage)) +
  geom_line() +
  labs(
    x = "Hour",
    y = "Leverage"
  )

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