ST303 Stochastic Simumlation

Class 7

Philipp Sterzinger

p.sterzinger@lse.ac.uk

Office hours: Tuesday, COL 5.03

19 November 2025

Recap

Variance reduction

Monte Carlo integral approximation

\[ \begin{aligned} I(g) &= \int_0^1 g(x) dx \\ &= \textrm{E}_{\mathcal{U}(0,1)}\left[g(U) \right] \\ &= \lim\limits_{n \to \infty} \frac{1}{n} \sum_{i = 1}^n g(U_i) \quad (\textrm{almost surely}) \\ & \approx \frac{1}{n} \sum_{i = 1}^n g(U_i) \\ & \equiv \hat{I}(g) \end{aligned} \]

Approximation error

\(\textrm{var}\left(\hat{I}(g) - I(g) \right) = \mathcal{O}(n^{-1})\) but may be substantial in small samples

Variance reduction

Suppose we have two random variables \(X,Y\) with the same mean and variance \(\sigma^2\). Then

\[ \textrm{var} \left(\frac{X + Y}{2} \right) = \frac{1}{2} \sigma^2 (1 + \rho) < \sigma^2 \iff \rho < 1 \,, \] where \(\rho \in (-1, 1)\) is the correlation coefficient.

Variance reduction

BUT: We can always double the sample size in our simulation

Hence, we should compare \((X+Y) / 2\) to the mean of two independent samples of \(X\), say \(X_1, X_2\).

To acheive effective variance reduction, we thus need

\[ \textrm{var} \left(\frac{X + Y}{2} \right) = \frac{1}{2} \sigma^2 (1 + \rho) < \textrm{var} \left(\frac{X_1 + X_2}{2} \right) = \frac{1}{2} \sigma^2 \iff \rho < 0 \,, \]

Variance reduction

Antithetic variates

Idea: minimize variance of MC estimator \(\hat{I}(g)\) by cooking up second estimator \(\tilde{I}(g)\) with same mean and variance but negative correlation.

Theorem 1 If \(g(x)\) is a monotone function, and \(U \sim \mathcal{U}(0,1)\) a uniform random variable, then \[ \textrm{cov}(g(U), g(1-U)) \leq 0 \]

Thus, letting \(\tilde{I}(g) = \frac{1}{2n} \sum_{i = 1}^n g(U_i) + \frac{1}{2n} \sum_{i = 1}^n g(1 - U_i)\) we have that

\[ \textrm{var}\left( \tilde{I}(g) \right) < \textrm{var}\left( \hat{I}(g) \right) \]

Questions

Question 1

Improve your estimates for all integrals in Exercise \(3\) of Seminar \(1\) and Exercise \(2\) of Seminar \(2\) by the method of Antithetic Variables. Approximate the new standard errors.

Question 1.a

\[ I = \int_0^1 \exp\{-x^{1.5}\} dx = \textrm{E}_{\mathcal{U}(0,1)} \left[ g(U) \right] \] for \(g(x) = \exp\{-x^{1.5}\}\).

Standard MC

\[ \hat{I} = \frac{1}{n}\sum_{i = 1}^n g(U_i), \quad U_i \overset{\textrm{i.i.d.}}{\sim} \mathcal{U}(0, 1) \]

Question 1.a

Antithetic variates

\[ \tilde{I} = \frac{1}{2n} \sum_{i = 1}^n g(U_i) + \frac{1}{2n} \sum_{i = 1}^n g(1 - U_i), \quad U_i \overset{\textrm{i.i.d.}}{\sim} \mathcal{U}(0, 1) \]

Monotonicity:

Code

f <- function(x){exp(-1.5 * x)} 
plot.new() 
curve(f, from = 0, to = 1, lwd = 2, col = "red", ylab = "g(x)")

Question 1.a

Antithetic variates

Code

g <- function(x){exp(-x^(1.5))} 

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

U <- runif(N) 

mc_mean <- mean(g(U)) 
mc_sd <- sd(g(U)) 

av_vals <- 0.5 * g(U) + 0.5 * g(1 - U)
av_mean <- mean(av_vals) 
av_sd <- sd(av_vals) 

reps <- 1000
Us <- matrix(runif(reps * N), nrow = N, ncol = reps) 

av_means <- vapply(1:reps, function(col){U <- Us[, col]; mean(0.5 * g(U) + 0.5 * g(1 - U))}, 0.0)
mc_means <- vapply(1:reps, function(col){U <- Us[, col]; mean(g(U))}, 0.0)

av <- cbind(av_mean, av_sd) 
mc <- cbind(mc_mean, mc_sd) 

comparison <- rbind(mc, av) 
colnames(comparison) <- c("Mean", "Standard deviation") 
rownames(comparison) <- c("Monte Carlo", "Antithetic variates") 

par(mfrow = c(1, 2))

hist(g(U), col = "#5b5be0b0", freq = FALSE, breaks = 100, xlab = "MC/AV summands", ylab = "Distribution of Summands", main = "Histograms of MC summands") 
hist(av_vals, freq=FALSE, breaks=100, col="#e05b5bb0", border="#e05b5bb0", add = TRUE) 
abline(v = 0.699792, col = "green", lty = 2, lwd = 2) 

legend("topright", 
       legend = c("MC", "AV", "truth"), 
       fill = c("#5b5be0b0", "#e05b5bb0", NA), 
       border = NA, 
       lty = c(NA, NA, 2), 
       lwd = c(NA, NA, 2), 
       pch = c(NA, NA, NA), 
       col = c(NA, NA, "green"))

hist(mc_means, col = "#5b5be0b0",  freq = FALSE,  ylab = "Distribution of MC estimates", xlab = "MC means", main = "Histogram of MC estimates") 
hist(av_means, col = "#e05b5bb0",  freq = FALSE,  ylab = "Distribution of MC/AV estimates", xlab = "MC/AV means", main = "Histograms of MC/AV estimates", add = TRUE, border = "#e05b5bb0") 
abline(v = 0.699792, col = "green", lty = 2, lwd = 2) 
legend("topright", 
       legend = c("MC", "AV", "truth"), 
       fill = c("#5b5be0b0", "#e05b5bb0", NA), 
       border = NA, 
       lty = c(NA, NA, 2), 
       lwd = c(NA, NA, 2), 
       pch = c(NA, NA, NA), 
       col = c(NA, NA, "green"))

Code

comparison

                         Mean Standard deviation
Monte Carlo         0.7003113        0.196876091
Antithetic variates 0.6998029        0.004407909

Question 1.b

\[ \begin{aligned} I &= \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 \\ &= \textrm{E}_{\mathcal{U}(0,1)} \left[ h(U) \right] \end{aligned} \] for \[ h(y) = \frac{1}{y^2}\exp\left\{-\left(\frac{1}{y}-1\right)^{1.5}\right\} \]

Question 1.b

Monotonicity of transformation

Code

h <- function(x){exp(- (1 / x - 1)^(3 / 2)) / x^2}  
plot.new() 
curve(h, from = 0, to = 1, lwd = 2, col = "red", ylab = "h(x)")

Question 1.b

Code

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

U <- runif(N) 

cat("cov(h(u), h(1-u)): ", cov(h(U), h(1-U)))

cov(h(u), h(1-u)):  -0.1782384

Code

mc_mean <- mean(h(U)) 
mc_sd <- sd(h(U)) 

av_vals <- 0.5 * h(U) + 0.5 * h(1 - U)
av_mean <- mean(av_vals) 
av_sd <- sd(av_vals) 


av <- cbind(av_mean, av_sd) 
mc <- cbind(mc_mean, mc_sd) 

comparison <- rbind(mc, av) 
colnames(comparison) <- c("Mean", "Standard deviation") 
rownames(comparison) <- c("Monte Carlo", "Antithetic variates") 

hist(g(U), 
  col = "#5b5be0b0", 
  freq = FALSE, breaks = 100, 
  xlab = "Estimates", 
  ylab = "Distribution of Estimates", 
  main = "Histograms of estimates", 
  xlim = c(.3, max(c(h(U), av_vals)))) 
hist(av_vals, freq=FALSE, breaks=100, col="#e05b5bb0", add = TRUE) 
abline(v = 0.902745, col = "green", lty = 2, lwd = 2) 

legend("topright", 
       legend = c("MC", "AV", "truth"), 
       fill = c("#5b5be0b0", "#e05b5bb0", NA), 
       border = NA, 
       lty = c(NA, NA, 2), 
       lwd = c(NA, NA, 2), 
       pch = c(NA, NA, NA), 
       col = c(NA, NA, "green"))

Code

comparison

                         Mean Standard deviation
Monte Carlo         0.9015331          0.6264433
Antithetic variates 0.9032674          0.3275354

Question 2

The loss to insurance portfolio \(A\) is of the form \[ S_A = \sum_{i = 1}^{n_A} X_i \,, \] where \(N_A\) is a Poisson random variable with parameter \(10\) and \(X_i\) are i.i.d random variables that are log-normal with parameters \(2\) and \(8\) (\(8\) is the variance of the underlying normal, not the standard deviation). The loss to another (independent of \(A\)) insurance portfolio \(B\) is of the form \[ S_B = \sum_{i = 1}^{n_B} Y_i \,, \] where \(N_B\) is a Poisson random variable with parameter \(10\) and \(Y_i\) are i.i.d random variables that are log-normal with parameter \(0\) and \(12\).

Calculate \(\Pr(S_A > S_B)\) and give a \(95\%\) confidence interval for your estimate.

Question 2

Simulating \(S_A\)

Draw \(n_A\) from the Poisson distribution
Given \(n_A = n\), draw \(n\) independent samples \(X_i\)
Compute \(\tilde{S}_A = \sum_{i = 1}^n X_i\)

\[ \begin{aligned} \Pr \left(\tilde{S}_A \leq t \right) &= \sum_{k = 0}^\infty \Pr\left(\sum_{i = 1}^k X_i \leq t {\Bigg |} n_A = k \right) \Pr(n_A = k) \\ &= \sum_{k = 0}^\infty \Pr\left(\sum_{i = 1}^k X_i \leq t \cap n_A = k \right) \\ &= \Pr\left(\bigcup_{k = 0}^\infty \left\{\sum_{i = 1}^{k} X_i \leq t \cap n_A = k \right\} \right) \\ &= \Pr\left(\sum_{i = 1}^{n_A} X_i \leq t \right) \\ &= \Pr(S_A \leq t) \end{aligned} \]

Question 2

\[ \begin{aligned} \Pr(S_A > S_B) &= \mathrm{E} \left[\mathbb{1} \left\{S_A > S_B \right\} \right] \\ & \overset{a.s.}{=}\lim\limits_{n \to \infty} \frac{1}{n} \sum_{i = 1}^n \mathbb{1} \left\{S_A > S_B \right\} \\ & \approx \frac{1}{n} \sum_{i = 1}^n \mathbb{1} \left\{S_A > S_B \right\} \end{aligned} \]

Algorithm

For \(n\) reps, do:

Draw lengths \(n_A\), \(n_B\) from respective Poisson distributions
Draw \(n_A (n_B)\) samples of \(X_i (Y_i)\)
Compute \(S_A\), \(S_B\)

\(\Pr(S_A > S_B) \approx \frac{1}{n} \sum_{i = 1}^n \mathbb{1} \left\{S_A > S_B \right\}\)

Question 2

Code

## Method 1 
insurance_loss <- function(N, lambdapois = 1, meanlog = 0, sdlog = 1){
  Ns <- rpois(N, lambdapois) 
  vapply(Ns, function(x){sum(exp(rnorm(x, mean = meanlog, sd = sdlog)))}, 0.0) 
}

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

ptm <- proc.time()
S_A <- insurance_loss(N, 10, 2, sqrt(8)) 
S_B <- insurance_loss(N, 10, 0, sqrt(12))   
sim_time_method_1 <- (proc.time() - ptm)[3]

pr_method_1 <- mean(S_A > S_B)

## Method 2 
CP_A <- function(x){
  n <- rpois(1, 10)
  x <- exp(rnorm(n, 2, sqrt(8)))
  s <- sum(x)
  return(s)
}

CP_B <- function(x){
  n <- rpois(1, 10)
  x <- exp(rnorm(n, 0, sqrt(12)))
  s <- sum(x)
  return(s)
}

set.seed(123)
ptm <- proc.time()
x <- 1:10^5 
CP_SA <- sapply(x, CP_A)
CP_SB <- sapply(x, CP_B)
sim_time_method_2 <- (proc.time() - ptm)[3]

pr_method_2 <- mean(CP_SA > CP_SB)


## Method 3 
set.seed(123) 
ptm <- proc.time()
CP_A <- rep(0, N)
CP_B <- rep(0, N)
N_A <- rpois(N, 10)
N_B <- rpois(N, 10)
for(i in 1 : N){
  CP_A[i] <- sum(exp(rnorm(N_A[i], 2, sqrt(8))))
  CP_B[i] <- sum(exp(rnorm(N_B[i], 0, sqrt(12))))
}

sim_time_method_3 <- (proc.time() - ptm)[3]
pr_method_3 <- mean(CP_A > CP_B)

cbind(pr_method_1, pr_method_2, pr_method_3)

     pr_method_1 pr_method_2 pr_method_3
[1,]     0.67875     0.67804     0.67954

Code

cbind(sim_time_method_1, sim_time_method_2, sim_time_method_3)

        sim_time_method_1 sim_time_method_2 sim_time_method_3
elapsed             0.245             0.321             0.178

Question 3

Same as Question 2 but now \(n_A = n_B \sim \textrm{Pois}(10)\).

Code

## Method 1 
insurance_loss <- function(N, lambdapois = 10, meanlog = c(2, 0), sdlog = c(sqrt(8), sqrt(12))){
  Ns <- rpois(N, lambdapois) 
  vapply(Ns, function(x){
     sum(exp(rnorm(x, meanlog[1], sdlog[1])) 
      - exp(rnorm(x, meanlog[2], sdlog[2])))}, 
    0.0) 
}

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

ptm <- proc.time()
ins_loss <- insurance_loss(N) 
sim_time_method_1 <- (proc.time() - ptm)[3]

pr_method_1 <- mean(ins_loss > 0)

## Method 2 
CP_AB <- function(x){
  n <- rpois(1, 10)
  x <- exp(rnorm(n, 2, sqrt(8)))
  y <- exp(rnorm(n, 0, sqrt(12)))
  s <- sum(x - y)
  return(s)
}

set.seed(123)
ptm <- proc.time()
x <- 1:10^5 
CP_SAB <- sapply(x, CP_AB)
sim_time_method_2 <- (proc.time() - ptm)[3]
pr_method_2 <- mean(CP_SAB > 0)


## Method 3 
set.seed(123) 

insurance_loss_loop <- function(N, lambdapois = 10, meanlog = c(2, 0), sdlog = c(sqrt(8), sqrt(12))){
  CP_SAB_loop <- rep(0, N)
  N_AB <- rpois(N, lambdapois)
  for(i in 1 : N){
    CP_SAB_loop[i] <- sum(exp(rnorm(N_AB[i], meanlog[1], sdlog[1])) 
      - exp(rnorm(N_AB[i], meanlog[2], sdlog[2])))
  }
  CP_SAB_loop 
}

ptm <- proc.time()
ins_loss_loop <- insurance_loss_loop(N) 
sim_time_method_3 <- (proc.time() - ptm)[3]
pr_method_3 <- mean(ins_loss_loop > 0)

cbind(pr_method_1, pr_method_2, pr_method_3)

     pr_method_1 pr_method_2 pr_method_3
[1,]     0.69181     0.69373     0.69181

Code

cbind(sim_time_method_1, sim_time_method_2, sim_time_method_3)

        sim_time_method_1 sim_time_method_2 sim_time_method_3
elapsed             0.222             0.232             0.181

Question 4

Improve your estimate in Question 2 and 3 by the method of antithetic variables. Do at least one of them and approximate the new standard errors.

Question 4

Recall that \[ X_i = e^{N_i}, \, N_i \sim \mathrm{N}(2,8) \] and note that \(4 - N_i \sim \mathrm{N}(2,8)\).

Now \[ 4 - N_i \sim \mathrm{N}(2,8) \]

and

\[ \textrm{cov}(N_i, 4 - N_i) = - \textrm{var}(N_i) \implies \rho = -1 \]

Futher \(\exp\{x\}\) is an increasing monotonic transformation, so that

\[ \textrm{cov}(\exp\{4 - N_i\}, \exp\{N_i\}) < 0 \]

Question 4

Hence, we can use \[ \bar{S}_A = \sum_{i = 1}^n \exp\{4 - N_i\} \] as the antithetic variate to \(S_A = \sum_{i = 1}^n \exp\{N_i\}\)

Similarly, for \(S_B\) we use the antithetic variate \[ \bar{S}_B = \sum_{i = 1}^n \exp\{ - M_i\} \] for \(S_B =\sum_{i = 1}^n \exp\{M_i\}\), \(M_i \sim \mathcal{N}(0, 12)\)

Question 4

Estimate \(\Pr(S_A > S_B)\) as

\[ \Pr(S_A > S_B) = \textrm{E} [\mathbb{1}\{S_A > S_B\}] \approx \frac{1}{N} \sum_{i = 1}^N \frac{\mathbb{1}\left\{S_A^i > S_B^i \right\} + \mathbb{1}\left\{\bar{S}_A^i > \bar{S}_B^i \right\}}{2} \]

Code

## Method 1 
pr_insurance_loss <- function(N, lambdapois = 10, meanlog = c(2, 0), sdlog = c(sqrt(8), sqrt(12))){
  Ns <- rpois(N, lambdapois) 
  loss <- vapply(Ns, function(x){
      X <- rnorm(x, meanlog[1], sdlog[1])
      Y <-  rnorm(x, meanlog[2], sdlog[2])
      loss <- sum(exp(X)) - sum(exp(Y))  
      antithetic_loss <-   sum(exp(4 - X)) - sum(exp(-Y))   
      ((loss > 0) + (antithetic_loss > 0)) / 2
  }, 
  0.0) 
  c(mean(loss), sd(loss))
}

set.seed(123)
av_est <- pr_insurance_loss(10^5)

ins_loss <- insurance_loss(10^5) > 0
mc_est <- c(mean(ins_loss), sd(ins_loss)) 

rbind(av_est, mc_est)

          [,1]      [,2]
av_est 0.69204 0.3140726
mc_est 0.69342 0.4610757

ST303 Stochastic Simumlation Class 7

Recap

Variance reduction

Monte Carlo integral approximation

Approximation error

Variance reduction

Variance reduction

Variance reduction

Antithetic variates

Questions

Question 1

Question 1.a

Standard MC

Question 1.a

Antithetic variates

Question 1.a

Antithetic variates

Question 1.b

Question 1.b

Monotonicity of transformation

Question 1.b

Question 2

Question 2

Simulating \(S_A\)

Question 2

Algorithm

Question 2

Question 3

Question 4

Question 4

Question 4

Question 4

ST303 Stochastic Simumlation

Class 7