This example investigates the performance of R in comparison to Julia language. Additionally shows how to easily call Julia inside R code. With that being said, we will load JuliaCall library that enables us to do so. Alternatively, there is also XRJulia library available.

library(JuliaCall)

It is necessery to tell R where is Julia.exe stored, so the loaded library can communicate with it accordingly.

julia_setup(JULIA_HOME = "path to the julia.exe")
# sometimes JuliaCall package requires to also set there a working directory 
# setwd("path to the folder with julia.exe")

Next we will define function that price European options based on Monte Carlo simulation. Namely, simulates random path of prices many times for a given set of parameters and calculates value of option based on expected value of simulated payoffs.

The model of motion of prices will be a following commonly used stochastic differential equation, called geometric Brownian motion:

\[dS_t = \mu S_t dt + \sigma S_t dW_t\] Where \(W_t\), \(\mu\) and \(\sigma\) are respectively a Wiener process, drift and volatility.

Price dynamics model as well as chosen derivative is not very complex. The point of herein document is to compare R and Julia performance on some well-known example. As show below, it is also easy to price such an option with the same price dynamics using a well known workhorse in finance, namely Black-Scholes formula

MC_option_pricing <- function(K, S, r, mu, sig, t, n_sim, call_put){
  
  # defining function that generate geometric brownian motion
  rand_walk <- function(mu, sig, n_sim, r){ cumprod((rnorm(n_sim, mu, sig) + 1)/(r + 1)) }
  
  # initializing empty matrix
  simulations <- matrix(0, ncol = n_sim, nrow = t)
  
  # simulating random walk for a given number of times
  for (i in 1:n_sim) {  simulations[,i] <- rand_walk(mu, sig, t, r) }
  
  # different payoff for put and call option
  if (call_put == 0) { payoff <- simulations[t,] * K - S
    
  } else if (call_put == 1) {
    
    payoff <- S - simulations[t,] * K}
  
  # calculating expected value of option
  mean(pmax(payoff, 0))}

Now we can check if the function is working properly. To do so, we will define Black-Scholes formula with a following equation for pricing call options:

\[ C(S_T, t) = N(d_1)S_t - N(d_2)PV(K)\] \[d_1 = \frac{1}{\sigma \sqrt{T} }\bigg[ln \bigg(\frac{S_t}{K}\bigg) + \bigg(r + \frac{\sigma^2}{2} \bigg)t \bigg]\] \[d_2 = d_1 - \sigma \sqrt{t}\] \[PV(K) = Ke^{-rt}\] Where,

BlackScholes <- function(K, S, r, sig, T, call_put){
  
  d1 <- (log(S/K) + (r + sig^2/2)*T) / (sig*sqrt(T))
  d2 <- d1 - sig*sqrt(T)
  
  if(call_put == 1){
    value <- S*pnorm(d1) - K*exp(-r*T)*pnorm(d2)
    return(value)}
  
  if(call_put == 0){
    value <-  (K*exp(-r*T)*pnorm(-d2) - S*pnorm(-d1))
    return(value)}}

S <- 40
K <- 50
mu <- 0
sig <- 0.03
T <- 100
r <- 0.001

BlackScholes(K, S, r, sig, T, 0); MC_option_pricing(K, S, r, mu, sig, T, 10e3, call_put = 0)
## [1] 8.121264
## [1] 7.967556

Prices are almost the same. The difference is insignificant and is only related to the number of simulations executed by MC_option_pricing function. Results from this function will converge to the result of analytical function as number of simulations increases.

We will use julia_command() function to define very similar function but in Julia, which can be used in R.

julia_command("
function MC_option_pricing(K, S, mu, r, sig, t, n_sim, call_put)
              
function rand_walk(mu, sig, t, r)
    dist = Normal(mu, sig)
    rets = rand(dist, t)
    rets = cumprod([x + 1 for x in rets] .* (1 + r)) 
    return rets
end

simulations = zeros(t, n_sim)

for i in 1:n_sim
  simulations[1:t, i] = rand_walk(mu, sig, t, r)
end

if call_put == 0
  payoff = simulations[t,:] .* K .- S
else
  payoff = S .- simulations[t,:] .* K
end

payoff[findall(payoff .<= 0)] .= 0
return mean(payoff)
end")
## MC_option_pricing (generic function with 1 method)

Alternatively, we can import whole script to R.

#julia_source("path/monte_carlo_option_pricing.jl")

Next, we can define variables for both of the functions to compute. The most important will be the number of simulations that a particular function have to perform. This will ultimately determine how computationally expensive function will be.

mu <- 0 # drift term
sig <- 0.02 # volatility (standard deviation)
t <- 100 # time to expiration
n_sim <- 10^6 # number of simulations
K <- 40 # strike price at expiration
S <- 60 # stock price at time tt = 0
call_put <- 0 # 1 for call option, 0 for put
r <- 0.002 # daily risk free rate

Script for Julia have dependency, so we should also load a Distributions library.

julia_library("Distributions")

Now we will use both functions to price put option based on 1000000 simulations, each with 100 rows/observations.

Note that Julia is sensitive to object types and sometimes need to have objects specified in explicit way. (use of as.integer() function).

r_time <- system.time( r_price <- MC_option_pricing(K, S, mu, r, sig, t, n_sim, call_put = 1) )

julia_time <- system.time(julia_price <- julia_call("MC_option_pricing", K, S, mu, r, sig, as.integer(t), as.integer(n_sim), 1)  )

As before, Monte Carlo method have stochastic properities so the price vary a little. That is why, in order to obtain very precise value, it is neccessery to execute many simulations. According to R the price is 11.992 and Julia said it costs 12.006, a difference of 0.014.

The computation took 15.87 seconds for R and 6.42 seconds for Julia, a difference of 9.45 seconds. Julia was 2.472 times faster.

It might be also interesting how long it takes to execute these functions many times but with fewer simulations. Package microbenchmark enables to do it easily.

library(microbenchmark)
set.seed(123)

test <- microbenchmark(R = MC_option_pricing(K, S, mu, r, sig, t, 10^5, call_put = 1),
               Julia = julia_call("MC_option_pricing", K, S, mu, r, sig, as.integer(t), as.integer(10^5), 1))

summary(test)
##    expr       min        lq      mean    median        uq       max neval
## 1     R 1577.6810 1600.7104 1634.3757 1619.5581 1642.2906 1978.7081   100
## 2 Julia  462.8935  557.4521  562.7707  563.0522  570.5963  677.0913   100

As we see in the summary, Julia function is on average 2.904 faster than R (by a median it’s x 2.876). Over 100 repeats, none of R executions was faster than the slowest one of Julia.

library(ggplot2)

autoplot(test)+
  labs(title = "R vs Julia", 
       subtitle = "Average time of calculating put option value with monte carlo simulation (10000 simulations)")