# ST303 Stochastic Simumlation
Class 8

Philipp Sterzinger

p.sterzinger@lse.ac.uk
Office hours: Thursday, COL 8.05

04 December 2025

# Recap

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

# Questions

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

Code

    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

Code

    # 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

Code

    # 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

Code

    ## 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

Code

    ## 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 \]

Code

    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