Monte Carlo simulations: Anthitetic variates
In this post I will discuss different variance reduction methods applied to option pricing theory with the help of Monte Carlo simulations.
There are different methods to improve efficiency of Monte Carlo simulations. One way is to increase the number of experiments, n. It is obvious that increasing n will lead to convergence of error size to 0, however, it comes with a tradeoff in computational time.
Let’s denote confidence level alpha as a constant in 0.05.
From the first computation of error estimate, we can conclude that as n goes to inifinty, error estimate will converge into 0. Increasing n has possible drawbacks, for example:
- Monte Carlo experiments are random in nature. Using law of large numbers approach will cost us low accuracy.
- Exponential time complexity slows down computations, which is a high cost especially in dynamic stock markets.
What if we don’t want to change n? Then, we should focus on numerator of this equation, C.
If we look closely to C, we can see that one of the main multipliers is variance. And today’s topic is how to decrease it. Before going into it, I would like to also take a look at the final equation above, number of experiments. Decrease in variance, a multiplier of C will also decrease n. We can conclude that computational time is also proportional to the variance.
For changing variance, we should transform our random variable Y in a way that it’s expected value remains same, but variance decreases.
This kind of transformations are available through several methods:
- Anthitetic variates
- Control Variates
- Importance Sampling
- Stratified sampling and etc.
Antithetic Variates.
As it can be concluded by the name, it has something to do with opposite sign. And that is what it does. We find another random variable Y_new which is negatively correlated(anthitetic) to our originial variable Y_old. They have same expected value, as they have same distribution. Then we create a new random variable Z and assign it to average of Y_new and Y_old.
From the equations we can see step by step how variance of Z becomes smaller that original Y_old.
Smart reader can ask: What about computational time? If we create 2 random variables, doesn’t it slow down our computer 2 times more?
That’s right! If you look closely to correlation coefficient, you can find the answer how to find the spot for cutoff rate. If we have p<0, variance will decrease more than 2 times. As decrease in variance is proportional to time complexity, we should always choose p smaller that 0 to gain a speedup.
How to generate values of anthitetic variate?
- Look at the originial distribution and find the mean.
- Subtract original variable from 2*mean. That’s your Y_new’s values.
Implementation.
- Create a generator which returns a matrix M with coulumns Y_old and Y_new
- Define the function of interest. For ex, g(Y)=5*y+4, etc.
- Define a new function for returning average of anthitetic variates: g2(Y_old, Y_new)=(g(Y_old)+g(Y_new))/2
- Use Monte Carlo simulation with generator and g2.
Simple Example.
MC2 <- function(g,Xgen,error,alpha){
n0 <- 1000 #how many variables to generate at each iteration
sum_Y <- 0
sum_Y2 <- 0
n <- 0
error_estimate <- error+1
while(error_estimate>error){
X <- Xgen(n0)
n <- n+n0
Y <- g(X)
sum_Y <- sum_Y+sum(Y)
sum_Y2 <- sum_Y2+sum(Y**2)
sigma_y <- sqrt(abs(sum_Y2-sum_Y**2/n)/(n-1))
error_estimate <- -qnorm(alpha/2)*sigma_y/sqrt(n)
}
return(c(sum_Y/n,n))
}gen=function(n){
return(runif(n, min=0,max=2))
}
g=function(x){
return((x+x^3)/(1-0.1*sqrt(x)))
}e=0.02
alpha=0.05answer=MC2(g,gen,e,alpha)gen_2=function(n){
y_old=runif(n, min=0,max=2)
y_new=2-runif(n, min=0,max=2)
return(cbind(y_old,y_new))
}g_2=function(x){
return((g(x[,1])+g(x[,2]))/2)
}answer2=MC2(g_2,gen_2,e,alpha)
Pricing stock option
We want to compute expected value of an option price with payoff function given. If the payoff function is monotone, it is highly likely out variance reduction will work better:
As a generator, I will use euler’s function (which I will discuss in another post). Main takeaway is: by transforming eurler’s function to return a matrix of 2 variables, and calculating with a new g function I will get similar answer with less number of computations and more precision
MC2 <- function(g,genX,error,alpha){
n0 <- 1000 #how many variables to generate at each iteration
sum_Y <- 0
sum_Y2 <- 0
n <- 0
error_estimate <- error+1
while(error_estimate>error){
X <- genX(n0)
n <- n+n0
Y <- g(X)
sum_Y <- sum_Y+sum(Y)
sum_Y2 <- sum_Y2+sum(Y**2)
sigma_y <- sqrt(abs(sum_Y2-sum_Y**2/n)/(n-1))
error_estimate <- -qnorm(alpha/2)*sigma_y/sqrt(n)
}
return(c(sum_Y/n,n))
}S_euler <- function(n,S0,m,T,mu,sigma){
dt <- T/m
t <- seq(0,T,length.out=m+1)
S <- matrix(S0,ncol=n,nrow=m+1)
for(i in 1:m){
dB <-sqrt(dt)*rnorm(n)
S[i+1,]<- S[i,]*(1+mu(t=t[i],s=S[i,])*dt+sigma(t=t[i],s=S[i,])*dB)
}
return(S[m+1,])
}E <-43
r <- 0.02
D <- 0.01
T <-0.5
S0 <- 49
sigma <- function(t,s){
return(0.3+04/(1+0.01*(s-50)^2))
}
mu <- function(t,s){
return(r-D)
}
g <- function(s){
return(exp(-r*T)*pmax(s-E,0))
}
m <- 20
alpha <- 0.05
error <- 0.05
gen <- function(n){
return(S_euler(n,S0,m,T,mu,sigma))
}answer1 <- MC2(g,gen,error,alpha)
#use antithetic variatesS_euler_antithetic <- function(n,S0,m,T,mu,sigma){
dt <- T/m
t <- seq(0,T,length.out=m+1)
S <- matrix(S0,ncol=n,nrow=m+1)
Stilde <- matrix(S0,ncol=n,nrow=m+1) #antithetic
for(i in 1:m){
dB <-sqrt(dt)*rnorm(n)
S[i+1,]<- S[i,]*(1+mu(t=t[i],s=S[i,])*dt+sigma(t=t[i],s=S[i,])*dB)
Stilde[i+1,]<- Stilde[i,]*(1+mu(t=t[i],s=Stilde[i,])*dt+sigma(t=t[i],s=Stilde[i,])*(-dB))
}
return(cbind(S[m+1,],Stilde[m+1,]))
}gen_antithetic <- function(n){
return(S_euler_antithetic(n,S0,m,T,mu,sigma))
}
gen_antithetic(10)
g_antithetic <- function(x){
return((g(x[,1])+g(x[,2]))/2)
}
answer2 <-MC2(g_antithetic,gen_antithetic,error,alpha)
answer2[1]-answer1[1]
#compute speedup factor
answer1[2]/(2*answer2[2])
