# ST303 Stochastic Simumlation
Class 10

Philipp Sterzinger

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

11 December 2025

# Recap

## Brownian motion

- Starting point: \(x_0\)
- Independent Increments: \(X_t - X_s \perp \mathcal{F}_s\) (\(\mathcal{F}_s \approx\) info at time \(s \leq t\))
- Distribution of increments: \(X_t - X_s \sim \mathcal{N}(\mu (t - s), \sigma^2(t - s))\)

## Brownian motion

### Simulation

- Discretise the time line on which we want to simulate the Brownian motion
- Simulate the increments at the time points

Code

    brownian_path <- function(grid_step = 0.01, start = 0, stop = 2, seed = 123){
        set.seed(seed)
        time <- seq(from = start, to = stop, by = grid_step)
        num_steps <- length(time)
        rvs <- sqrt(grid_step) * rnorm(num_steps)
        cumsum(rvs)
    }

    N <- 10
    time <- seq(0, 2, 0.01)
    brownian_paths <- vapply(1:N, function(seed){brownian_path(seed = seed)}, numeric(length(time)))
    matplot(time, brownian_paths, type = "l")

## Homogenous Poisson process

- Counting process \(\{N(t), t \leq 0\}\) which counts the number of events up to time \(t\)
- Number of events in interval \((a, b]\) is an independent Poisson random variable with parameter \((b - a) \lambda\) (not useful for simulation)
- Waiting times \(W_1,\ldots, W_k\) are i.i.d. exponential random variables with parameter \(\lambda\) (wasteful simulation)

### Conditional distribution of event times

\(T_1,\ldots,T_k \mid N(t) = k\) follows the distribution of an ordered sample \(t_1 \leq t_2 \leq \ldots \leq t_k\) with \(t_i \sim \mathcal{U}(0, t)\)

Simulating \(N(t)\):

1. Draw a Poisson random variable \(N(t)\) with parameter \(t \lambda\)
2. Draw \(k = N(t)\) uniform random variables on \((0,t)\)
3. Order the sample, these are our event times

# Questions

## Question 1

> Simulate a standard Brownian motion from time \(0\) to time \(2\) using a grid with steps of size \(0.01\). We are interested in the median of the path. It is known that its distribution is the same as the distribution of \(|X| − |Y |\), where \(X\) and \(Y\) are two independent standard normal random variables.

> 1. Generate a sample of size \(10000\) for the median of the Brownian motion and a sample of the same size from the distribution of \(|X| − |Y|\).

> 2. Compare the two samples with various plots and calculating various quantities. The reason for doing this is the following. We noticed that there is a bias when generating the maximum of a Brownian motion and we want to see if other path related quantities such as the median are also affected.

## Question 1.a

Code

    N <- 10^5
    time <- seq(0, 2, 0.01)
    brownian_paths <- vapply(1:N, function(seed){brownian_path(seed = seed)}, numeric(length(time)))

    medians <- apply(brownian_paths, 2, median)
    med_br <- mean(medians)
    var_br <- var(medians) / N

    br <- c(med_br, var_br)

    set.seed(123)
    med_controls <- abs(rnorm(N)) - abs(rnorm(N))
    med_control <- mean(med_controls)
    var_control <- var(med_controls) / N

    cont <- c(med_control, var_control)

    df <- data.frame(cbind(br, cont))
    colnames(df) <- c("BR_med", "Control_med")
    rownames(df) <- c("mean", "var")

    print(df)

               BR_med   Control_med
    mean 4.592620e-04 -2.098764e-03
    var  7.324927e-06  7.332610e-06

Code

    plot(density(medians), col = "red", xlab = "median", main = "Probabiltiy Density")
    lines(density(med_controls), col = "blue")
    legend("topright", legend = c("Brownian motions median", "|X| - |Y|"),
           col = c("red", "blue"), lty = rep(1, 2), cex = 0.8)

## Question 2

> The loss to a general insurer is \[ S(t) = \sum_{i = 1}^{N(t)} Y_i \,, \] where \(N (t)\) is a Poisson process and \(\{Y_i\}_{i = 1,2, \ldots, N(t)}\) is a sequence of i.i.d random variables that is independent of \(N (t)\). The premium income is linear with rate \(c\). Let \(T\) be a fixed time.

> 1. Explain why we would like to generate samples with the same distribution as \[ Z = \max_{0 < t < T} \left( \sum_{i = 1}^{N(t)} Y_i - ct \right) \,. \]

> 2. Generate a sample from this distribution when \(T = 3\), the Poisson rate is \(10\), \(c = 11\) and \(Y_i\) has a Gamma distribution with parameter \(2\) and \(2\).

## Question 2.a

- Premium income at time \(t\): \(ct\)
- Loss at time \(t\): \(\sum_{i = 1}^{N(t)} Y_i\) (At each event up to time \(t\), incur loss \(Y_i\))

### Ruin time

First time \(t \in (0,T]\) where loss is greater than income

\[ t : \max_{0 < t < T} \left( \sum_{i = 1}^{N(t)} Y_i - ct \right) > 0 \,. \]

Random variable \(Z = \max_{0 < t < T} \left( \sum_{i = 1}^{N(t)} Y_i - ct \right)\) can be used to simulate the ruin time

## Question 2.b

- Number of claims: \(N(t) \sim \textrm{Pois}(10t)\)
- Time horizon: \((0, 3]\), so \(N(T) \sim \textrm{Pois}(30)\)
- Arrival times of claims: \(T_1,\ldots, T_{k} \mid N(t) = k\) ordered \(k\)-sample of \(\mathcal{U}(0,t)\) random variables
- Losses at event times: \(Y_i \sim \Gamma(2,2)\)
- Premium income up to event time \(T_i\): \(c T_i\)

Code

    sim_Z <- function(T, lambda = 10, gamma_shape = 2, gamma_rate = 2, ins_premium = 11){
      N_T <- rpois(1, T * lambda)

      T_s <- sort(T * runif(N_T))

      Y_s <- rgamma(N_T, shape = gamma_shape, rate = gamma_rate)

      claims <- cumsum(Y_s)
      income <- ins_premium * T_s

      max(claims - income) ## max can only be attained at jumps
    }

    sim_ruin <- function(T, lambda = 10, gamma_shape = 2, gamma_rate = 2, ins_premium = 11){
      N_T <- rpois(1, T * lambda)

      T_s <- sort(T * runif(N_T))

      Y_s <- rgamma(N_T, shape = gamma_shape, rate = gamma_rate)

      claims <- cumsum(Y_s)
      income <- ins_premium * T_s

      t_ruin <- which((claims - income) > 0)[1]

      if(length(t_ruin) == 1){
        T_s[t_ruin]
      }else{
        NaN
      }
    }

    set.seed(123)
    N <- 10^5
    Z_s <- vapply(1:N, function(x){sim_Z(3)}, 0.0)
    cat("Probability of ruin in (0, 3] with initial capital = 10: ", sum(Z_s > 10) / N)

    Probability of ruin in (0, 3] with initial capital = 10:  0.07786

Code

    t_ruin_s <- vapply(1:N, function(x){sim_ruin(3)}, 0.0)
    hist(t_ruin_s, freq = FALSE)
    lines(density(t_ruin_s[!is.na(t_ruin_s)]))