Control variables

Want to estimate some quantity \(\theta = \textrm{E}(X)\), where \(X\) is a simulation outcome

For sample quantity \(Y\) such that \(\textrm{E}(Y) = 0\), then \(X + cY\) is also an unbiased estimate of \(\theta\)

Minimising variance

\[ \min_{c} \textrm{var}(X + cY) = \min_{c} \textrm{var}(X) + c^2 \textrm{var}(Y) + 2 c\textrm{cov}(X, Y) \]

\[ 2 c \textrm{var}(Y) + 2 \textrm{cov}(X, Y) = 0 \iff c^* = \frac{\textrm{cov}(X, Y)}{\textrm{var}(Y)} \]

Then we get that

\[ \textrm{var}(X + c^*Y) = \textrm{var}(X) - \frac{(\textrm{cov}(X, Y))^2}{\textrm{var}(Y)} \]

Control variables

Sample quantities

\[ \hat{c} = \frac{\frac{1}{n - 1} \sum_{i = 1}^n (X_i - \bar{X})(Y_i - \bar{Y})}{\textrm{var}(Y)} \quad (\textrm{var}(Y) \textrm{ known}) \]

\[ \hat{c} = \frac{ \sum_{i = 1}^n (X_i - \bar{X})(Y_i - \bar{Y})}{\sum_{i = 1}^n (Y_i - \bar{Y})^2} \quad (\textrm{var}(Y) \textrm{ unknown}) \]

Control variables

Choosing \(Y\)

Find approximation \(\tilde{X}\) of \(X\) with \(\textrm{E}(\tilde{X}) = \mu\), \(\mu\) known

Then hopefully we have a high correlation between \(X\) and \(\tilde{X}\)

Set \(Y = \tilde{X} - \mu\)

Possible choice:

If \(X = \frac{1}{n} \sum_{i = 1}^n h(U_i)\), we can approximate \(h(U)\) using a first or second order Taylor series \(T(U)\) and let \(Y = T(U) - \textrm{E}(T(U))\)

Conditioning

\[ \textrm{var}(X) = \textrm{var}(\textrm{E}[X \mid Y]) + \textrm{E}(\textrm{var}(X \mid Y)) \]

so that

\[ \textrm{var}(\textrm{E}[X \mid Y]) \leq \textrm{var}(X) \]

Useful if we deal with compositions of random variables, e.g. \(\mathbb{1}\{X + Y > a\}\), where we can condition on one and then simplify

Stratified sampling

\[ F_X(x) = \sum_{i = 1}^J \alpha_j F_j(x) \]

• Interested in \(\mathrm E(X) =\mathrm E_J(\mathrm E(X\mid J = j)) = \sum_{i = 1}^J \alpha_j \mu_j\)

• Convential MC: \(\mathrm E(X) \approx \frac{1}{n} \sum_{i = 1}^n \mathbb{1}\{J_i = j\} X_i^{(j)}\)

• Conditioning: \(\bar{X} = \frac{1}{n} \sum_{j = 1}^J \sum_{i = 1}^{\alpha_j n} X_i^{(j)}\)

• \(\textrm{var}(\bar X) = \frac{1}{n}\sum_{j = 1}^J \alpha_j \textrm{var}(X \mid J = j) \leq \frac{1}{n} \textrm{var}(X)\)

Stratified sampling

Idea: Break a sum or integral into distinct regions and sample from within these. Then average over all regions

\[ \int_0^1 h(x) dx \]

• Break \((0,1)\) into regions \([0, 1/n], (1/n, 2/n], \ldots, ((n-1)/n, n]\)
• Within each interval randomly sample from \(\mathcal{U}((j-1)/n, j / n)\), compute \(h(U)\)
• Average all samples

The variance of this estimator is \[\textrm{E}_{J}(\textrm{var}(h(U) \mid J = j)) \leq \textrm{var}(h(U))\,,\] where the conditioning variable \(J\) identifies the interval from which we are sampling.

Question 1

In section 5.2, we tried to estimate \[ \int_0^1 \frac{1}{1 + x} dx = \log(2) \approx 0.6931472\,, \] using the method of control variables. We used a uniform \((0, 1)\) random variable \(U\) and estimated the expected value of \(\frac{1}{1 + U}\) and \(\frac{1}{2} - U\) as control variable. A suggestion was made at the end of the example, that maybe \(Y = \frac{1}{3} - U + \frac{U^2}{2}\) is a better control variable. Try it.

N <- 10^5 
set.seed(123) 

U <- runif(N) 

X <- 1 / (1 + U) 

## Old option 
Y_1 <- 0.5 - U 
c_1 <- - cov(X, Y_1) / var(Y_1) 

CV_1 <- mean(X + c_1 * Y_1 )

## New option 
Y_2 <- 1 / 3 - U + U^2 / 2 
c_2 <- - cov(X, Y_2) / var(Y_2) 

CV_2 <- mean(X + c_2 * Y_2)

ln2 <- log(2) 

cbind(CV_1, CV_2, ln2) 

          CV_1      CV_2       ln2
[1,] 0.6930806 0.6931596 0.6931472

Question 1

# Assessing variance 
cv_comparison <- function(seed, N = 10^4){
    set.seed(seed) 
    U <- runif(N) 
    X <- 1 / (1 + U)
    Y_1 <- 0.5 - U 
    c_1 <- - cov(X, Y_1) / var(Y_1) 
    CV_1 <- X + c_1 * Y_1
    Y_2 <- 1 / 3 - U + U^2 / 2 
    c_2 <- - cov(X, Y_2) / var(Y_2) 
    CV_2 <- X + c_2 * Y_2
    cbind(CV_1, CV_2) 
}

cvs <- cv_comparison(seed = 123, N = 10^4)

mean_row <- c(colMeans(cvs), log(2)) 
var_row <- c(apply(cvs, 2, var), 0.0)

df <- data.frame(
  CV_1 = c(mean_row[1], var_row[1]),
  CV_2 = c(mean_row[2], var_row[2]),
  ln_2 = c(mean_row[3], var_row[3])
)
rownames(df) <- c("mean", "var")
print(df)

             CV_1         CV_2      ln_2
mean 0.6928039693 0.6933639547 0.6931472
var  0.0006139236 0.0001312156 0.0000000

Question 2

Consider the integral \[ \int_0^1 \exp \left\{-x^{1.5} \right\}dx \,. \] Try to estimate by estimating the expected value of \(\exp \left\{-U^{1.5} \right\}\) with \(Y = \frac{2}{5}− U^{1.5}\) as control variable. \(U\) is a uniform \((0, 1)\) random variable.

Question 2

How we got here

\[ \exp{x} \approx 1 + x \] for small \(x\) (Taylor series of degree one)

Hence \[ \exp \left\{-U^{1.5} \right\} \approx 1 - U^{1.5} \]

\[ \textrm{E}(U^{1.5}) = \int_0^1 x^{1.5} dx = \left[\frac{2}{5} x^{5 / 2} \right]_0^1 = \frac{2}{5} \]

\[ Y = 1 - U^{1.5} - \textrm{E}\left[1 - U^{1.5}\right] = 1 - U^{1.5} - 1 + \frac{2}{5} = \frac{2}{5} - U^{1.5} \]

Question 2

# Assessing variance 
cv <- function(seed, N = 10^5){
    set.seed(seed) 
    U <- runif(N) 
    X <- exp(-U^1.5)
    Y <- 2 / 5 - U^1.5
    c <- - cov(X, Y) / var(Y) 
    cbind(X, X + c * Y)
}


cvs <- cv(seed = 123, N = 10^4)

mean_row <- c(colMeans(cvs), 0.699792327761494) 
var_row <- c(apply(cvs, 2, var), 0.0)

df <- data.frame(
  CV_1 = c(mean_row[1], var_row[1]),
  CV_2 = c(mean_row[2], var_row[2]),
  truth = c(mean_row[3], var_row[3])
)
rownames(df) <- c("mean", "var")
print(df)

           CV_1         CV_2     truth
mean 0.70145244 0.6994077390 0.6997923
var  0.03829033 0.0006029402 0.0000000

Question 3

Improve your estimate of Question 5, Exercise 2 using the conditioning method in section 5.3.1.

Let \(Y_t\) be the yield of a fund at \(t\) year and \(Y_t = \exp(X_t)\), where \(X_1, X_2, \ldots , X_t\) is an independent identically distributed sequence of normal random variables with mean \(0.09\) and variance \(0.1\).

An investor invests an amount \(1\) at times \(t = 0, 1, 2\) and wants to know what is the probability that the accumulation at \(t = 3\) will be higher than \(4\) for a vector of \(10^6\) possible returns.

\[ \textrm{Accumulation at } t_3 \, = \exp\{X_1\} + \exp\{X_1 + X_2\} + \exp\{X_1 + X_2 + X_3\} \]

Question 3

\[ \begin{aligned} A_3 &= \exp\{X_1\} + \exp\{X_1 + X_2\} + \exp\{X_1 + X_2 + X_3\} \\ &= \exp\{0.09 +\sqrt{0.1}Z_1\} \\ &+ \exp\{0.18 + \sqrt{0.1}(Z_1 + Z_2)\} \\ &+ \exp\{0.27 + \sqrt{0.1}(Z_1 + Z_2 + Z_3)\} \\ &= \exp\{\sqrt{0.1}Z_1\} \left(\exp\{0.09\} + \exp\{0.18 + \sqrt{0.1}Z_2\} + \exp\{0.27 + \sqrt{0.1}(Z_2 + Z_3)\} \right) \end{aligned} \]

\(Z_i \sim \mathrm{N}(0,1)\)

Conditioning

\(h(Z_2, Z_3) = \exp\{0.09\} + \exp\{0.18 + \sqrt{0.1}Z_2\} + \exp\{0.27 + \sqrt{0.1}(Z_2 + Z_3)\}\)

Question 3

\[ \begin{aligned} \Pr \left(A_3 > a \right) &= \textrm{E} \left(\mathbb{1} \{A_3 > a \} \right) \\ &= \textrm{E} \left(\mathbb{1} \{ \exp\{\sqrt{0.1} Z_1 \}h(Z_2,Z_3) > a \} \right) \\ &= \textrm{E} \left(\textrm{E} \left[ \mathbb{1} \left\{ \exp\{\sqrt{0.1} Z_1 \}h(Z_2,Z_3) > a \right\} \mid Z_2, Z_3 \right] \right) \\ &= \textrm{E} \left(\textrm{E} \left[ \mathbb{1} \left\{ \exp\{\sqrt{0.1} Z_1 \} > \frac{a}{h(Z_2,Z_3)}\right\}\mid Z_2, Z_3 \right] \right) \\ &= \textrm{E} \left(\textrm{E} \left[ \mathbb{1} \left\{ \sqrt{0.1} Z_1 > \log \left( \frac{a}{h(Z_2,Z_3)} \right) \right\}\mid Z_2, Z_3 \right] \right) \\ &= \textrm{E} \left(\textrm{E} \left[ \mathbb{1} \left\{ Z_1 > \frac{\log(a) - \log(h(Z_2, Z_3))}{\sqrt{0.1}} \right\} \mid Z_2, Z_3 \right] \right) \\ &= \textrm{E} \left( \Pr \left( Z_1 > \frac{\log(a) - \log(h(Z_2, Z_3))}{\sqrt{0.1}} \mid Z_2, Z_3 \right) \right) \\ &= \textrm{E} \left( \Psi \left(- \frac{\log(a) - \log(h(Z_2, Z_3))}{\sqrt{0.1}} \right) \right) \\ &\approx \frac{1}{n} \sum_{i=1}^n \Psi \left(- \frac{\log(a) - \log(h(Z_{2,i}, Z_{3,i}))}{\sqrt{0.1}} \right) \end{aligned} \]

Question 3

Algorithm

1. Sample \(n\)-vectors \(Z_2, Z_3\) from standard normal distribution
2. Compute \(n\)-vector \(h(Z_2,Z_3)\)
3. Compute \(n\)-vector \(\Psi \left(- \frac{\log(a) - \log(h(Z_{2}, Z_{3}))}{\sqrt{0.1}} \right)\)
4. Take average

## Old method for comparison 
prb_old <- function(threshold, N = 10^4, seed = 123){
  set.seed(seed)
  X1 <- rnorm(N, 0.09, sqrt(0.1)) 
  X2 <- rnorm(N, 0.09, sqrt(0.1))
  X3 <- rnorm(N, 0.09, sqrt(0.1)) 
  X <- exp(X1 + X2 + X3) + exp(X2 + X3) + exp(X3)
  X > threshold
}

## Conditioning method 
prb_new <- function(threshold, N = 10^4, seed = 123){
  set.seed(seed)
  Z2 <- rnorm(N)
  Z3 <- rnorm(N) 
  h <- exp(0.09) + exp(0.18 + sqrt(0.1) * Z2) + exp(0.27 + sqrt(0.1) * (Z2 + Z3)) 
  p <- pnorm(- (log(threshold) - log(h)) / sqrt(0.1)) 
  p
}



old_probs <- prb_old(threshold = 4, seed = 4)
new_probs <-  prb_new(threshold = 4, seed = 4)

mean_row <- colMeans(cbind(old_probs, new_probs))
var_row <-  apply(cbind(old_probs, new_probs), 2, var)

df <- data.frame(
  MC = c(mean_row[1], var_row[1]),
  conditioning = c(mean_row[2], var_row[2])
)
rownames(df) <- c("mean", "var")
print(df)

            MC conditioning
mean 0.4136000   0.41009092
var  0.2425593   0.06044932

Question 4

In section 5.3.2 in stratified sampling, we estimated \[ \int_0^1 h(x)dx \] with a uniform \((0, 1)\) sample of size \(n\) and calculate \[ \frac{1}{n} \sum_{j = 1}^n h\left( \frac{U_j + j - 1}{n} \right) \,. \] An improvement is \[ \frac{1}{2n} \sum_{j = 1}^n h\left( \frac{U_j + j - 1}{n} \right) + \frac{1}{2n} \sum_{j = 1}^n h\left( \frac{j - U_j}{n} \right) \,. \]

Question 4

Try them for the integrals \[ \int_0^1 \exp \{-x^{1.5}\} dx, \quad \int_0^\infty \exp \{-x^{1.5}\} dx \,. \]

Question 4

First integral is straightforward

## MC integral for comparison 
MC_int <- function(seed = 123, N = 10^2){
  set.seed(seed)
  U <- runif(N) 
  exp(-U^(1.5))
}

## Stratified sampling 
strat_int <- function(seed = 123, N = 10^2){
  set.seed(seed)
  U_strat <- (runif(N) + 1:N - 1) / N 
  exp(-U_strat^(1.5))
}

## Stratified sampling + antithetic variates 
strat_atv_int <- function(seed = 123, N = 10^2){
  set.seed(seed)
  U <- runif(N) 
  U_strat_1 <- (U + 1:N - 1) / N 
  U_strat_2 <- (1:N - U) / N 
  0.5 * exp(-U_strat_1^(1.5)) + 0.5 * exp(-U_strat_2^(1.5))
}

## AV 
atv_int <-function(seed = 123, N = 10^2){
  set.seed(seed)
  U <- runif(N) 
  0.5 * exp(-U^1.5) + 0.5 * exp(-(1 - U)^1.5)
}


mc_ints <- MC_int(seed = 123, N = 10^2)
strat_ints <-  strat_int(seed = 123, N = 10^2)
strat_atv_ints <-  strat_atv_int(seed = 123, N = 10^2)
av_ints <- atv_int(seed = 123, N = 10^2)

mean_row <- apply(cbind(mc_ints, strat_ints, strat_atv_ints, av_ints), 2, mean)
var_row <-  apply(cbind(mc_ints, strat_ints, strat_atv_ints, av_ints), 2, var)

df <- data.frame(
  MC = c(mean_row[1], var_row[1]),
  stratified = c(mean_row[2], var_row[2]), 
  stratified_atv = c(mean_row[3], var_row[3]), 
  int = c(0.699792, 0.0) 
)
rownames(df) <- c("mean", "var")
print(df)

             MC stratified stratified_atv      int
mean 0.70135346 0.69985020     0.69979218 0.699792
var  0.03821751 0.03913765     0.03920729 0.000000

Question 4

Estimating actual variances

              MC   stratified stratified_atv
var 0.0004043397 3.572664e-08   4.159339e-13

Question 4

For the second integral, recall from Class 1 that, using \(y = \frac{1}{1 + x}\), we have

\[ \int_0^\infty \exp \{-x^{1.5}\} dx = \int_0^1 \frac{1}{y^2} \exp\left\{-\left(\frac{1}{y} - 1\right)^{1.5}\right\} dy \]

integrand <- function(x){1 / x^2 * exp(-(1 / x - 1)^(1.5))} 

## MC integral for comparison 
MC_int <- function(seed = 123, N = 10^2){
  set.seed(seed)
  U <- runif(N) 
  integrand(U)
}

## Stratified sampling 
strat_int <- function(seed = 123, N = 10^2){
  set.seed(seed)
  U_strat <- (runif(N) + 1:N - 1) / N 
  integrand(U_strat)
}

## Stratified sampling + antithetic variates 
strat_atv_int <- function(seed = 123, N = 10^2){
  set.seed(seed)
  U <- runif(N) 
  U_strat_1 <- (U + 1:N - 1) / N 
  U_strat_2 <- (1:N - U) / N 
  0.5 * integrand(U_strat_1) + 0.5 * integrand(U_strat_2)
}

seed <- 123
N <- 10^4

mc_ints <- MC_int(seed = seed, N = N)
strat_ints <-  strat_int(seed = seed, N = N)
strat_atv_ints <-  strat_atv_int(seed = seed, N = N)

mean_row <- apply(cbind(mc_ints, strat_ints, strat_atv_ints), 2, mean)
var_row <-  apply(cbind(mc_ints, strat_ints, strat_atv_ints), 2, var)

df <- data.frame(
  MC = c(mean_row[1], var_row[1]),
  stratified = c(mean_row[2], var_row[2]), 
  stratified_atv = c(mean_row[3], var_row[3]), 
  int = c(0.902745, 0.0) 
)
rownames(df) <- c("mean", "var")
print(df)

            MC stratified stratified_atv      int
mean 0.9044889  0.9027455      0.9027453 0.902745
var  0.3934238  0.3929561      0.3929563 0.000000