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 \mathrm{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. \end{aligned} \] Set \(U = (X-2)^3/64 +1\), we have \[ (X-2)^3 = 64U-64 \implies X = (64U-64)^{\frac{1}{3}} +2. \] Hence the inverse transform algorithm is as follows:
1) Generate \(U \sim \mathrm{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 number
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. \end{aligned} \] Since \((y-2)^2\) is decreasing for \(y \in [-2,2]\) while \(e^{y^2/2}\) is maximized at \(y=-2\) and \(y=2\), we have \[ 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}. \] The algorithm is as follows:
1) Generate \(Y \sim \mathrm{N}(0,1)\).
2) If \(Y \notin [-2,2]\), return to step 1 (due to different support).
3) Generate \(U \sim \mathrm{U}(0,1)\).
4) If \(U \leq \frac{1}{c} \cdot \frac{f_X(Y)}{f_Y(Y)} = (Y-2)^2 \exp(0.5Y^2-2)/16\), 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 <= (Y-2)^2 *exp(0.5*Y^2-2) /16) {
      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 = 1/12\). The algorithm is as follows:
1) Generate \(Z_j \sim \mathrm{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)} = \sum_{j=1}^n S_{t_j}^{(i)}/n\) for \(i=1,...,2m\).
6) Compute \(\hat\theta_{\mathrm{ant}} = e^{-rT} \sum_{i=1}^{2m} \max(A_T^{(i)}-S_T^{(i)}, 0)/(2m)\).

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.102948   8.111333
apply(estimate,2,var)
##   standard antithetic 
## 0.09079188 0.03474802

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 the control variable to value this Asian option.
(modified from RMSC4001 2015S final Q5)

Solution

Since calculating \(\mathbb{E}(A_{t_9})\) can be difficult, we propose to estimate it via pilot simulation:
1) Generate \(Z_j \sim \mathrm{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,...,k\), where \(k\) is the sample size of pilot simulation.
4) Compute \(A_{t_9}^{(i)} = \sum_{j=1}^9 S_{t_j}^{(i)}/9\) for \(i=1,...,k\).
5) Compute \(A_T^{(i)} = \sum_{j=1}^n S_{t_j}^{(i)}/n\) for \(i=1,...,k\).
6) 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 \mathrm{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)} = \sum_{j=1}^9 S_{t_j}^{(i)}/9\) for \(i=1,...,m\).
5) Compute \(A_T^{(i)} = \sum_{j=1}^n S_{t_j}^{(i)}/n\) for \(i=1,...,m\).
6) Compute \(\hat{\pi}_T^{(i)} = \max(A_T^{(i)}-S_T^{(i)},0) -\hat\sigma_{A_{t_9},\max(A_T-S_T,0)} (A_{t_9}^{(i)} -\hat\mu_{A_{t_9}}) /\hat\sigma^2_{A_{t_9}}\).
7) Compute \(\hat\theta_{\mathrm{ctr}} = e^{-rT} \sum_{i=1}^{m} \hat{\pi}_T^{(i)}/m\).

We verify the algorithm with the following R implementation:

# Use result of the above
colnames(estimate) = c("standard","control")
# pilot simulation
kPilot = 1000
S = matrix(100, nrow=kPilot, ncol=n+1) #set S0=100
for (i in 1:n) {
  Z = rnorm(kPilot)
  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.1438001
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:10])
  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.090873 8.117448
apply(estimate,2,var)
##  standard   control 
## 0.1070240 0.1040946

We can see that the price is similar but the variance reduction effect is very small (or even none). This is actually expected from the pilot simulation as correlation between \(A_{t_9}\) and the payoff is small. Therefore, it is important to choose a good control variable in practice.
(the original question has a better control variable 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 \mathrm{N}(0,1)\). The likelihood ratio is \[ \frac{f_X(y)}{f_Y(y)} = \frac{2}{\sqrt{2\pi}} e^{-\frac{1}{2} y^2+|y|}. \] 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=1/2\). To generate \(Y\), 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. \end{aligned} \] Set \(U = F_Y(Y)\), we have \[ Y = \left\{ \begin{array}{ll} -\ln(2-2U), & U \geq 0.5; \\ \ln(2U), & U < 0.5. \end{array} \right. \] Hence the algorithm is as follows:
1) Generate \(U_i \sim \mathrm{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_{\mathrm{imp}} = e^{-rT} \sum_{i=1}^{n} h_i L_i/n\).

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.784974   2.790092
apply(estimate,2,var)
##   standard importance 
## 0.04057495 0.03645575

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