R/simulations.R

In shill1729/KellyCriterion: Log-optimal stock allocation

#' Simulate Kelly strategy for generalized coin toss
#'
#' @param bankroll initial bankroll to bet with
#' @param p probability of winning each trial
#' @param a payout on winning toss on top of wager
#' @param b wager lost on each trial on losing toss
#' @param trials number of trials to simulate
#'
#' @description {Simulating the generalized coin tossing gambling game.}
#' @details {Start off with a given value of wealth and bet the Kelly fraction each round.}
#' @return vector
#' @export sim_binary
sim_binary <- function(bankroll, p, a, b, trials = 100)
{
  # Initial wealth process
  z <- matrix(data = 0, nrow = trials+1)
  z[1] <- bankroll
  y <- kelly_binary(p, a, b)

  # Simulate the Bernoulli outcomes 0,1 (loss, win)
  outcomes <- rbinom(n = trials, size = 1, prob = p)
  # Convert to bet outcomes:
  epsilon <- a*outcomes-b*(1-outcomes)
  for(i in 1:(trials))
  {
    z[i+1] <- z[i]+epsilon[i]*z[i]*y
  }
  return(z)
}

#' Simulate a series of IID trials of moneyline under log-optimal growth
#'
#' @param bankroll initial bankroll to bet
#' @param p true chance of outcome
#' @param wagers the vector of wagers, see details
#' @param trials the number of trials to simulate
#'
#' @description {Simulate a series of moneyline trials under log-optimal growth.}
#' @details {The argument \code{wagers} must contain four wagers, the first two the odds
#' for the favorite, the latter two the odds for the underdog.}
#' @return vector
#' @export sim_moneyline
sim_moneyline <- function(bankroll, p, wagers, trials = 100)
{
  a <- wagers[1]
  b <- wagers[2]
  u <- wagers[3]
  v <- wagers[4]

  # Initial wealth process
  z <- matrix(data = 0, nrow = trials+1)
  z[1] <- bankroll
  # Optimal strategy
  x <- kelly_moneyline(p, a, b, u, v)
  y <- 1-x

  # Simulate the Bernoulli outcomes 0,1 (loss, win)
  outcomes <- rbinom(n = trials, size = 1, prob = p)
  # Convert to bet outcomes:
  epsilon <- a*outcomes-b*(1-outcomes) # Favorite wins
  delta <- u*(1-outcomes)-v*outcomes # Underdog wins
  for(i in 1:(trials))
  {
    z[i+1] <- z[i]+epsilon[i]*z[i]*x+delta[i]*z[i]*y
  }
  return(z)
}


#' Simulate log-optimal strategy on mixture diffusions
#'
#' @param bankroll initial bankroll to invest
#' @param t time horizon to trade over
#' @param spot the initial stock price
#' @param rate the return of the bond
#' @param parameters matrix of parameters of mixture, see details
#'
#' @description {Simulate log-optimal strategy for mixture diffusion prices}
#' @details {The matrix must contain a row of probabilities, a row of means, and
#' a row of standard deviations.}
#' @return data.frame of solution
#' @export sim_mixdiff
sim_mixdiff <- function(bankroll, t, spot, rate, parameters)
{
  probs <- parameters[1, ]
  mus <- parameters[2, ]
  sigmas <- parameters[3, ]

  # Drift and volatility coefficients for mixture diffusion
  mu <- function(t, s) sde::drift_lvm(s, t, probs, mus, sigmas, spot)
  volat <- function(t, s) sde::volat_lvm(s, t, probs, mus, sigmas, spot)
  IC <- list(s = spot, x = bankroll)
  control <- function(t, s) kelly_mixdiff(t, s, rate, parameters, spot)

  f <- list(function(t, s, x) mu(t, s)*s,
            function(t, s, x) (rate+(mu(t, s)-rate)*control(t, s))*x
  )
  g <- list(function(t, s, x) volat(t, s)*s,
            function(t, s, x) control(t, s)*volat(t, s)*x
  )
  z <- sde::sde_system(f, g, IC, NULL, t0 = 0, tn = t, n = 1000)
  sol <- z
  sol$s <- log(sol$s)-log(IC$s)
  sol$x <- log(sol$x)-log(IC$x)
  graphics::par(mfrow = c(1, 2))
  odeSolveR::plot_trajectories(sol, legend_names = c("Stock", "Portfolio"),
                               legend_size = 0.4, legend_loc = "topleft")
  odeSolveR::plot_phase_portrait(sol)
  return(z)
}

#' Simulate log-optimal strategy for DTFM
#'
#' @param n number of days to simulate
#' @param spot initial stock price
#' @param bankroll initial bankroll
#' @param distr distribution of daily arithmetic returns
#' @param param parameters of the distribution
#' @param rate risk-neutral rate
#' @param plotG boolean for plotting on call
#'
#' @description {Simulates daily stock prices under the model and implements the strategy.
#' Assuming no trading costs, etc. Returns mean and total growth plus pnl.}
#' @return list
#' @importFrom numerics rgmm rstable
#' @export sim_dtfm
sim_dtfm <- function(n, spot, bankroll, distr, param, rate = 0, plotG = TRUE)
{
  rdistr <- paste("r", distr, sep = "")
  rdistr <- get(rdistr)
  x <- kelly_dtfm(distr, param, rate)
  if(distr == "norm" || distr == "unif")
  {
    args <- c(n, as.list(param))
  } else if(distr == "gmm")
  {
    args <- list(n, param[1, ], param[2, ], param[3, ])

  } else if(distr == "stable")
  {
    args <- list(n, param)
  }
  r <- do.call(rdistr, args)
  s <- spot*c(1, cumprod(1+r))
  p <- bankroll*c(1, cumprod(1+rate+x*(r-rate)))
  sx <- log(s/spot)
  px <- log(p/bankroll)
  bounds <- c(min(px, sx), max(px, sx))

  if(plotG)
  {
    plot(px, ylim = bounds, type = "l", ylab = c("Total-log return"), xlab = "Day")
    graphics::lines(sx, col = "red")
  }
  result <- list(total_growth = px[n+1], profit = p[n+1]-p[1])

  return(result)
}