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.

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] 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.

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.

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.