Monte Carlo option pricing - comparison of R and Julia languages
Published:
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,
- $N(\cdot)$ is a C.D.F. of a standard normal distribution
- $S_t$ - Current price of underlying asset
- $K$ - Option strike price
- $t$ - time to option expiry
- $r$ - Interest rate of risk-free asset
- $\sigma$ - Returns volatility of underlying asset price
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] 8.097749Prices 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.993 and Julia said it costs 12.011, a difference of 0.018.
The computation took 18.56 seconds for R and 4.97 seconds for Julia, a difference of 13.59 seconds. Julia was 3.734 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 1471.9590 1539.2106 1700.6300 1569.5385 1736.1726 2794.069 100
## 2 Julia 514.8352 590.9109 641.1392 607.7939 635.8842 1272.219 100As we see in the summary, Julia function is on average 2.653 faster than R (by a median it’s x 2.582). 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)")