functional-simulations.R

# This example simulates a stock price 1 time and 10k times.
# It uses a functional approach and a loop approach.
# A speed test is done at the end.

# This example for GBM ignores the fact that 
# 1) There is an explicit solution to GBM where simulation isn't needed
# 2) Along the same lines, there is a vectorized form where pure matrix mult can be used.
# These are ignored for the sake of the example, because there are often cases
# where this isn't the case and you HAVE to do the discrete loop.

library(purrr)

n_steps <- 100
r <- 0.02
sigma <- .25
dt <- .01
normal_increments <- rnorm(n_steps)
s_0 <- 50

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Geometric brownian motion - functional
# 1 simulation

# The accepted answer of:
# https://quant.stackexchange.com/questions/7125/simulation-of-gbm

# Focusing on the RHS of the simulation equation
# No indexing needed
gbm_increment <- function(s_t, Z_t, r, sigma, dt) {
  s_t + r * s_t * dt + sigma * s_t * sqrt(dt) * Z_t
}

sims <- accumulate(
  .x = normal_increments, 
  .f = gbm_increment, 
  r = r, sigma = sigma, dt = dt, 
  .init = s_0
)

head(sims)
#> [1] 50.00000 50.78881 51.12444 51.55642 49.20616 51.56984

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Geometric brownian motion - loop
# 1 simulation

# Notes:
# Have to create a "holder" s_t vector
# Loop requires keeping track of indexing, potentially error prone
# Might be clearer for beginners to understand with the explicit indexing

s_t <- vector("numeric", length = 101L)
s_t[1] <- s_0

for(i in seq_len(n_steps)) {
  s_t[i+1] <- s_t[i] + r * s_t[i] * dt + sigma * s_t[i] * sqrt(dt) * normal_increments[i]
}

head(s_t)
#> [1] 50.00000 50.78881 51.12444 51.55642 49.20616 51.56984

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Geometric brownian motion - functional
# 10k simulations

# What if I wanted to do 10k sims?
n_sims <- 10000

# A 100 element list where each element is 10000 rnorm numbers
# Each element represents a step, so that we what we accumulate over
normal_increments_sims <- replicate(n_steps, rnorm(n_sims), simplify = FALSE)

# This code doesnt change!
sims_list <- accumulate(
  .x = normal_increments_sims, 
  .f = gbm_increment, 
  r = r, sigma = sigma, dt = dt, 
  .init = s_0
)

# Bind all the steps into 1 matrix
sims <- do.call(cbind, sims_list)

sims[1:10, 1:10]
#>       [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
#>  [1,]   50 48.82488 49.28438 49.84095 50.29673 48.55607 48.87270 46.80843
#>  [2,]   50 49.77376 51.37026 54.94362 53.63327 56.37143 57.15921 56.57089
#>  [3,]   50 50.85879 51.10557 51.84187 53.47533 55.29510 55.47696 54.98492
#>  [4,]   50 52.50587 53.79625 54.42480 55.75047 55.51358 57.65000 56.42236
#>  [5,]   50 47.54502 47.80591 48.42428 47.98959 49.15185 48.50118 48.47622
#>  [6,]   50 50.00669 50.75497 50.84408 51.38521 50.44118 49.02976 46.60753
#>  [7,]   50 46.71363 46.38224 46.90848 44.99075 42.96552 41.52559 40.48895
#>  [8,]   50 46.49652 45.73792 45.63420 45.90111 47.64764 47.18049 44.56368
#>  [9,]   50 49.45114 49.10725 49.45577 51.40958 50.84070 51.07022 51.10069
#> [10,]   50 50.18387 51.45386 52.06531 53.59077 53.72025 53.14685 56.07085
#>           [,9]    [,10]
#>  [1,] 47.56217 46.72485
#>  [2,] 56.20900 57.19786
#>  [3,] 57.12318 56.79404
#>  [4,] 55.44803 54.14621
#>  [5,] 48.51532 49.02127
#>  [6,] 46.64433 46.26615
#>  [7,] 39.37709 41.23779
#>  [8,] 45.30542 44.71432
#>  [9,] 52.00825 52.67660
#> [10,] 54.95495 57.11199

dim(sims)
#> [1] 10000   101

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Geometric brownian motion - loop
# 10k simulations

normal_increments_mat <- matrix(rnorm(n_steps * n_sims), nrow = n_sims)

# Have to worry about setting up the correct dim sizes...
s_t <- matrix(0, nrow = n_sims, ncol = n_steps+1)
s_t[,1] <- s_0

# This changed a little bit from before. Now we do columns indexing
for(i in seq_len(n_steps)) {
  s_t[,i+1] <- s_t[,i] + r * s_t[,i] * dt + sigma * s_t[,i] * sqrt(dt) * normal_increments_mat[,i]
}

s_t[1:10, 1:10]
#>       [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
#>  [1,]   50 50.50125 49.89321 49.67687 49.48158 50.44275 49.87563 50.68925
#>  [2,]   50 50.48235 51.33004 49.47028 49.95901 50.55156 53.60091 55.46067
#>  [3,]   50 50.08999 47.91631 47.22120 44.74997 48.30522 48.82862 47.20558
#>  [4,]   50 51.48490 51.05863 52.40778 52.72697 52.88700 51.80046 51.37699
#>  [5,]   50 49.13463 50.50223 49.50573 51.47593 51.83414 51.53101 53.38273
#>  [6,]   50 48.40819 47.71540 46.88044 46.35158 43.55721 42.37722 42.39574
#>  [7,]   50 49.83166 48.55261 49.85788 48.29407 46.09949 46.67961 47.18883
#>  [8,]   50 48.16406 49.03634 47.73413 50.27268 49.79826 49.50948 51.22986
#>  [9,]   50 50.14662 48.58000 49.43702 47.11591 49.37879 51.72572 51.73143
#> [10,]   50 47.80738 47.41565 47.45949 48.58250 48.73871 50.46086 50.72911
#>           [,9]    [,10]
#>  [1,] 50.05278 49.84897
#>  [2,] 57.30995 58.35765
#>  [3,] 46.52706 47.48566
#>  [4,] 50.33741 50.56660
#>  [5,] 55.91655 56.16715
#>  [6,] 43.28351 42.30380
#>  [7,] 47.91071 46.37111
#>  [8,] 50.21492 48.24507
#>  [9,] 53.55376 53.44589
#> [10,] 49.36376 49.00284

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Speed test for the 10k sims

microbenchmark::microbenchmark(
  # Accumulate
  {
    sims_list <- accumulate(
      .x = normal_increments_sims, 
      .f = gbm_increment, 
      r = r, sigma = sigma, dt = dt, 
      .init = s_0
    )
    sims <- do.call(cbind, sims_list)
  },
  
  # Loop
  for(i in seq_len(n_steps)) {
    s_t[,i+1] <- s_t[,i] + r * s_t[,i] * dt + sigma * s_t[,i] * sqrt(dt) * normal_increments_mat[,i]
  },
  
  times = 50
)
#> Unit: milliseconds
#>                                                                                                                                                                         expr
#>  {     sims_list <- accumulate(.x = normal_increments_sims, .f = gbm_increment,          r = r, sigma = sigma, dt = dt, .init = s_0)     sims <- do.call(cbind, sims_list) }
#>                         for (i in seq_len(n_steps)) {     s_t[, i + 1] <- s_t[, i] + r * s_t[, i] * dt + sigma * s_t[,          i] * sqrt(dt) * normal_increments_mat[, i] }
#>        min        lq     mean   median       uq       max neval
#>   7.344971  9.206522 18.20160 10.71141 14.43719  63.02775    50
#>  37.378124 42.675036 58.35127 46.91390 74.80985 131.69133    50

#' Created on 2018-02-23 by the [reprex package](http://reprex.tidyverse.org) (v0.1.1.9000).