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

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

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

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

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\)

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

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)]))