Random variable generation

Exercise

Consider the quadratic random variable \(X\) with pdf: \[ \begin{aligned} f_X(x) &= \left\{ \begin{array}{ll} \frac{3}{64} (x-2)^2 &, x \in [-2,2] \\ 0 &, x \notin [-2,2] \end{array} \right. \\ \end{aligned} \] Propose an inverse transform algorithm and acceptance-rejection algorithm using \(Y \sim N(0,1)\) to generate \(X\)
(modified from RMSC4001 2014S final Q3)

Solution

Note that for \(x \in [-2,2]\), \[ \begin{aligned} F_X(x) &= \int_{-2}^x \frac{3}{64} (t-2)^2 dt \\ &= \left[ \frac{1}{64} (t-2)^3 \right]_{-2}^x \\ &= \frac{1}{64} (x-2)^3 +1 \\ U = \frac{1}{64} (X-2)^3 +1 \implies (X-2)^3 &= 64U-64 \\ X &= (64U-64)^{\frac{1}{3}} +2 \end{aligned} \] Hence the inverse transform algorithm is as follows:
1) Generate \(U \sim U(0,1)\)
2) Set \(X = (64U-64)^{\frac{1}{3}} +2\)

We verify the algorithm with the following R implementation:

n = 10000
U = runif(n)
X = -(64-64*U)^(1/3)+2 #R cannot handle cube root of negative
hist(X, freq=F)
px = seq(-2,2,by=0.01)
py = 3/64*(px-2)^2
lines(px, py, col="blue") #add target pdf of X

For the acceptance-rejection algorithm, note that \[ \begin{aligned} f_Y(y) &= \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2} y^2} \\ \frac{f_X(y)}{f_Y(y)} &= \left\{ \begin{array}{ll} \frac{3\sqrt{2\pi}}{64} (y-2)^2 e^{\frac{1}{2} y^2} &, y \in [-2,2] \\ 0 &, y \notin [-2,2] \end{array} \right. \\ c = \max_y \frac{f_X(y)}{f_Y(y)} &= \frac{f_X(-2)}{f_Y(-2)} = \frac{3 e^2 \sqrt{2\pi}}{4} \end{aligned} \] since \((y-2)^2\) is decreasing for \(y \in [-2,2]\) while \(e^{\frac{1}{2} y^2}\) is maximized at \(y=-2\) and \(y=2\). Hence the algorithm is as follows:
1) Generate \(Y \sim N(0,1)\)
2) If \(Y \notin [-2,2]\), return to step 1 (due to different support)
3) Generate \(U \sim U(0,1)\)
4) If \(U \leq \frac{1}{c} \cdot \frac{f_X(Y)}{f_Y(Y)} = \frac{1}{16} (Y-2)^2 e^{\frac{1}{2} Y^2-2}\), set \(X=Y\). Otherwise return to step 1

We verify the algorithm with the following R implementation:

# Use result of the above
X = vector(length=n)
for (i in 1:n)
{
  repeat
  {
    repeat
    {
      Y = rnorm(1)
      if (Y >= -2 && Y <= 2)
        break
    }
    U = runif(1)
    if (U <= 1/16 *(Y-2)^2 *exp(0.5*Y^2-2))
    {
      X[i] = Y
      break
    }
  }
}
hist(X, freq=F)
lines(px, py, col="blue")

Antithetic variables

Exercise

Consider an arithmetic floating strike Asian put option which payoff is \(\max(A_T-S_T, 0)\), where \(S_T\) is the terminal price of the underlying asset and, for \(t_{j+1}-t_j = \Delta t, j=0,1,...,n\) and \(t_n=T\), \[ A_{t_m} = \frac{1}{m} \sum_{j=1}^m S_{t_j},\quad \forall1 \leq m \leq n \] Suppose that \(T=1, n=12, r=0.02, \sigma=0.4, S_0=100\) and the underlying asset price follows GBM. Propose an algorithm with antithetic variables to value this Asian option
(modified from RMSC4001 2015S final Q5)

Solution

Note that \(\Delta t = \frac{1}{12}\). The algorithm is as follows:
1) Generate \(Z_j \sim N(0,1)\) for \(j=1,..,n\)
2) Set \(S_{t_j}^{(i)} = S_{t_{j-1}}^{(i)} \exp \left[ (r-\frac{1}{2} \sigma^2) \Delta t +\sigma \sqrt{\Delta t} Z_j \right]\) for \(j=1,..,n\)
3) Set \(S_{t_j}^{(m+i)} = S_{t_{j-1}}^{(m+i)} \exp \left[ (r-\frac{1}{2} \sigma^2) \Delta t -\sigma \sqrt{\Delta t} Z_j \right]\) for \(j=1,..,n\)
4) Repeat step 1 to 3 for \(i=1,...,m\)
5) Compute \(A_T^{(i)} = \frac{1}{n} \sum_{j=1}^n S_{t_j}^{(i)}\) for \(i=1,...,2m\)
6) Compute \(\hat\theta_{ant} = \frac{e^{-rT}}{2m} \sum_{i=1}^{2m} \max(A_T^{(i)}-S_T^{(i)}, 0)\)

We verify the algorithm with the following R implementation:

rm(list=ls()) #clear variables from last question
nEstimate = 1000
nSample = 1000
estimate = matrix(nrow=nEstimate, ncol=2)
colnames(estimate) = c("standard", "antithetic")
n = 12
dt = 1/n
r = 0.02
sigma = 0.4
thetaH.std = function(m){
  S = matrix(100, nrow=m, ncol=n+1) #set S0=100
  for (i in 1:n)
  {
    Z = rnorm(m)
    S[,i+1] = S[,i]*exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*Z)
  }
  A = rowMeans(S[,-1]) #remove S0 in rowMeans
  return(exp(-r)*mean(pmax(A-S[,n+1],0)))
}
thetaH.ant = function(m){
  S = matrix(100, nrow=m, ncol=n+1) #set S0=100
  for (i in 1:n)
  {
    Z = rnorm(m/2) #the other half are -Z for antithetic variables
    S[1:(m/2),i+1] = S[1:(m/2),i]*exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*Z)
    S[(m/2+1):m,i+1] = S[(m/2+1):m,i]*exp((r-0.5*sigma^2)*dt-sigma*sqrt(dt)*Z)
  }
  A = rowMeans(S[,-1]) #remove S0 in rowMeans
  return(exp(-r)*mean(pmax(A-S[,n+1],0)))
}
for (i in 1:nEstimate)
{
  estimate[i,1] = thetaH.std(nSample)
  estimate[i,2] = thetaH.ant(nSample)
}
apply(estimate,2,mean)
##   standard antithetic 
##   8.099422   8.113494
apply(estimate,2,var)
##   standard antithetic 
## 0.09172547 0.03439894

We can see that the price is similar but variance is reduced

Control variate

Exercise

With the same setting from last question, i.e. arithmetic floating strike Asian put option, propose an algorithm using \(A_{t_9}\) as control variate to value this Asian option
(modified from RMSC4001 2015S final Q5)

Solution

Since calculating \(E(A_{t_9})\) can be difficult, we propose to estimate it via pilot simulation:
1) Generate \(Z_j \sim N(0,1)\) for \(j=1,..,n\)
2) Set \(S_{t_j}^{(i)} = S_{t_{j-1}}^{(i)} \exp \left[ (r-\frac{1}{2} \sigma^2) \Delta t +\sigma \sqrt{\Delta t} Z_j \right]\) for \(j=1,..,n\)
3) Repeat step 1 to 2 for \(i=1,...,m\) (sample size of pilot)
4) Compute \(A_{t_9}^{(i)} = \frac{1}{9} \sum_{j=1}^9 S_{t_j}^{(i)}\) for \(i=1,...,m\)
5) Compute \(A_T^{(i)} = \frac{1}{n} \sum_{j=1}^n S_{t_j}^{(i)}\) for \(i=1,...,m\)
5) Compute pilot sample estimates \(\hat\mu_{A_{t_9}}, \hat\sigma^2_{A_{t_9}}, \hat\sigma_{A_{t_9},\max(A_T-S_T,0)}\)

The main algorithm is as follows:
1) Generate \(Z_j \sim N(0,1)\) for \(j=1,..,n\)
2) Set \(S_{t_j}^{(i)} = S_{t_{j-1}}^{(i)} \exp \left[ (r-\frac{1}{2} \sigma^2) \Delta t +\sigma \sqrt{\Delta t} Z_j \right]\) for \(j=1,..,n\)
3) Repeat step 1 to 2 for \(i=1,...,m\)
4) Compute \(A_{t_9}^{(i)} = \frac{1}{9} \sum_{j=1}^9 S_{t_j}^{(i)}\) for \(i=1,...,m\)
5) Compute \(A_T^{(i)} = \frac{1}{n} \sum_{j=1}^n S_{t_j}^{(i)}\) for \(i=1,...,m\)
6) Compute \(\hat{\pi}_T^{(i)} = \max(A_T^{(i)}-S_T^{(i)},0) -\frac{\hat\sigma_{A_{t_9},\max(A_T-S_T,0)}} {\hat\sigma^2_{A_{t_9}}} (A_{t_9}^{(i)} -\hat\mu_{A_{t_9}})\)
7) Compute \(\hat\theta_{ctr} = \frac{e^{-rT}}{m} \sum_{i=1}^{m} \hat{\pi}_T^{(i)}\)

We verify the algorithm with the following R implementation:

# Use result of the above
colnames(estimate) = c("standard", "control")
# pilot simulation
mPilot = 1000
S = matrix(100, nrow=mPilot, ncol=n+1) #set S0=100
for (i in 1:n)
{
  Z = rnorm(mPilot)
  S[,i+1] = S[,i]*exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*Z)
}
A9 = rowMeans(S[,2:10])
AT = rowMeans(S[,-1]) #remove S0 in rowMeans
muH = mean(A9)
cH = cov(A9,pmax(AT-S[,n+1],0))
vH = var(A9)
cor(A9,pmax(AT-S[,n+1],0)) #small correlation
## [1] -0.1189858
thetaH.ctr = function(m){
  # main simulation
  S = matrix(100, nrow=m, ncol=n+1) #set S0=100
  for (i in 1:n)
  {
    Z = rnorm(m)
    S[,i+1] = S[,i]*exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*Z)
  }
  A9 = rowMeans(S[,2:5])
  AT = rowMeans(S[,-1]) #remove S0 in rowMeans
  piH = pmax(AT-S[,n+1],0) -cH/vH *(A9-muH)
  return(exp(-r)*mean(piH)) 
}
for (i in 1:nEstimate)
{
  estimate[i,1] = thetaH.std(nSample)
  estimate[i,2] = thetaH.ctr(nSample)
}
apply(estimate,2,mean)
## standard  control 
## 8.102662 8.111057
apply(estimate,2,var)
##   standard    control 
## 0.10538334 0.09952177

We can see that the price is similar but variance reduction effect is very small (or even none). This is actually expected from pilot simulation as correlation between \(A_{t_9}\) and payoff is small. Therefore, it is important to choose a good control variate in practice
(the original question has a better control variate but its derivation is not covered in this course)

Importance sampling

Exercise

Consider the simulation of a European put option on a stock with current price 100, strike price 80, interest rate 3%, volatility 40% and maturity 6 months. Assume the Black-Scholes model for the underlying stock. Provide an importance sampling algorithm based on the Laplace random variable \(Y\) with pdf \[f_Y(y) = \frac{1}{2} e^{-|y|}\] (modified from RMSC4001 2018S final Q1)

Solution

Let \(X \sim N(0,1)\). The likelihood ratio is \[ \begin{aligned} \frac{f_X(y)}{f_Y(y)} &= \frac{2}{\sqrt{2\pi}} e^{-\frac{1}{2} y^2+|y|} \end{aligned} \] The target function is \[ h(X) = \max(K-S_T,0) = \max(K-S_0 e^{(r-\frac{1}{2} \sigma^2)T +\sigma \sqrt{T} X}, 0) \] where \(K=80, S_0=100, r=0.03, \sigma=0.4, T=\frac{1}{2}\). To generate \(Y \sim Laplace(0,1)\), note that \[ \begin{aligned} F_Y(y) &= \frac{1}{2} \int_{-\infty}^y e^{-|t|} dt \\ &= \left\{ \begin{array}{ll} 1-\frac{1}{2} e^{-y} &, y \geq 0 \\ \frac{1}{2} e^{y} &, y < 0 \end{array} \right. \\ U = F_Y(Y) \implies Y &= \left\{ \begin{array}{ll} -\ln(2-2U) &, U \geq 0.5 \\ \ln(2U) &, U < 0.5 \end{array} \right. \\ \end{aligned} \] Hence the algorithm is as follows:
1) Generate \(U_i \sim U(0,1)\)
2) If \(U_i \geq 0.5\), set \(Y_i=-\ln(2-2U_i)\). Otheriwse, set \(Y_i=\ln(2U_i)\)
3) Compute \(L_i = \frac{2}{\sqrt{2\pi}} e^{-\frac{1}{2} Y_i^2+|Y_i|}\)
4) Compute \(h_i = \max(K-S_0 e^{(r-\frac{1}{2} \sigma^2)T +\sigma \sqrt{T} Y_i}, 0)\)
5) Repeat step 1 to 4 for \(i=1,...,n\)
6) Compute \(\hat\theta_{imp} = \frac{e^{-rT}}{n} \sum_{i=1}^{n} h_i \cdot L_i\)

We verify the algorithm with the following R implementation:

rm(list=ls()) #clear variables from last question
nEstimate = 1000
nSample = 1000
estimate = matrix(nrow=nEstimate, ncol=2)
colnames(estimate) = c("standard", "importance")
K = 80
S0 = 100
T = 1/2
r = 0.03
sigma = 0.4
thetaH.std = function(n){
  X = rnorm(n)
  ST = S0*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*X)
  return(exp(-r*T)*mean(pmax(K-ST,0)))
}
thetaH.imp = function(n){
  U = runif(n)
  Y = ifelse(U<0.5, log(2*U), -log(2-2*U))
  L = 2/sqrt(2*pi) *exp(-0.5*Y^2+abs(Y))
  h = pmax(K-S0*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*Y),0)
  return(exp(-r*T)*mean(h*L))
}
for (i in 1:nEstimate)
{
  estimate[i,1] = thetaH.std(nSample)
  estimate[i,2] = thetaH.imp(nSample)
}
apply(estimate,2,mean)
##   standard importance 
##   2.781887   2.786455
apply(estimate,2,var)
##   standard importance 
## 0.04397127 0.03246934

We can see that the price is similar but variance is reduced as Laplace r.v. has more extreme value to “escape” from the out-of-the-money situation