5.1 Poisson process

In this chapter, we have only looked at independent random variables so far. Thus, the quantities we tried to compute were only dependent on the random variables’ marginal distributions. Obviously, not everything can be modeled via independent random variables which often complicates things analytically.

For instance, let us look at a stochastic process \(X = \{X(t), t \geq 0 \}\) which is nothing but a family of random variables indexed by a set \(T\) (here \(T = [0, \infty)\)). More specifically, let us look at a Poisson process \(N = \{N(t), t \geq 0 \}\) with rate \(\lambda > 0\) which is defined as a stochastic process that fulfills \[\begin{align*} N(t) = \max\bigg\{n \in \mathbb{N}\, \bigg\vert\, \sum_{k = 1}^n Y_k \leq t \bigg\} \end{align*}\] where \(Y_k\), \(k \in \mathbb{N}\), are iid. \(\text{Exp}(\lambda)\) random variables. Consequently, \(X(t_1)\) and \(X(t_2)\) are not independent as they both depend on the distribution of (multiple) \(Y_k\). A trajectory or path, i.e. a realization, of such a Poisson process may look something like this

Figure 5.1: Realization of a Poisson process with rate \(\lambda = 0.5\)

Notice that simulation of such a path is straight-forward. Simply simulate independent \(\text{Exp}(\lambda)\) random variables, compute the cumulative sum of these realizations and increase \(N\) by one whenever the cumulative sum increases.

With the Monte Carlo method it is easy to compute quantities such as \(\mathbb{E}\inf\{t > 0 \vert N(t) \geq 10\}\). Since the process always jumps by at most one, we simply simulate 10 exponential random variables and compute their sum. If we do this often enough and keep denoting the sum of the realizations, we can average these times to get an approximation of \(\mathbb{E}\inf\{t > 0 \vert N(t) \geq 10\}\).

set.seed(123)
threshold <- 10    # Threshold to reach
N <- 100000        # Number of paths
lambda <- 1        # Intensity of exponential rvs

time_to_threshold <- numeric(N)
for (i in seq_along(time_to_threshold)) {
  time_to_threshold[i] <- sum(rexp(threshold, lambda))
}
# Estimate
mean(time_to_threshold)
#> [1] 9.998061

Admittedly, we could have computed this ourselves by realizing that the sum of \(n\) iid. exponential random variables follows a gamma distribution whose expected value can be easily computed to be \(n\). Nevertheless, it is good to have seen this approach as it generalizes well.

Now, imagine that instead of making unit jumps, the Poisson process can make arbitrary jumps. This kind of thinking leads us to the notion of the compound Poisson process.

Suppose that \((\xi_n)_{n \in \mathbb{N}}\) is a sequence of iid. random variables. These will represent the jumps of the compound Poisson process. Thus, a compound Poisson process is a stochastic process \(Y = \{Y(t), t \geq 0 \}\) where \(Y(t)\) can be characterized by \[\begin{align*} Y(t) = \sum_{k = 1}^{N(t)} \xi_k \end{align*}\] where \(N = \{ N(t), t \geq 0 \}\) is a Poisson process.

If \(\xi_1 \sim \text{U}(-2, 2)\), then a path of a compound Poisson process with intensity \(\lambda = 0.5\) may look like this

Theoretically, any probabilities and expectations involving the compound Poisson process come with the additional complexity of having to know the distribution of the \(n\)-fold convolution of \(\xi_k\), \(k = 1, \ldots, n\) where \(n\) can be any arbitrary positive integer. As for simulations, the only real threshold we need to cross is that we have to be able simulate realizations from the distribution of \(\xi_1\).

In ruin theory, a part of actuarial science focusing on the possibility of an insurer’s financial ruin, compound Poisson processes are a popular tool. For instance, you may imagine that an insurer is vulnerable to claims it has to pay as part of its insurance contracts. These claims can be random in numbers and height. In return, the insurer gets a steady stream of premiums \(c > 0\) from its customers. Additionally, it may currently have a amount of cash \(c_0\).

Now, the so-called classical Cramér-Lundberg model describes the assets of an insurer over time via a stochastic process \(X = \{ X(t), t \geq 0 \}\) defined via \[\begin{align*} X(t) = c_0 + ct - \sum_{k = 1}^{N(t)} \xi_k \end{align*}\] where \(N = \{ N(t), t \geq 0 \}\) is a Poisson process with rate \(\lambda > 0\).

Naturally, one may ask what the probability is that the insurer will experience financial ruin in the next, say, 10 years, i.e. we are interested in the quantity \(\mathbb{P}(\inf\{ t \geq 0\, \vert\, X(t) < 0 \} \leq 10)\). Let us compute this quantity via the Monte Carlo method assuming that \(\xi_1 \sim \text{Exp}(\lambda_{\text{claims}})\), with \(\lambda_{\text{claims}} > 0\).

set.seed(123)
# Initial parameters
lambda_N <- 2
lambda_claims <- 3
c_0 <- 0.25
c <- 1
N <- 10000
time_horizon <- 10

# Monte Carlo simulation
ruin_occurred <- logical(N)
for (i in seq_along(ruin_occurred)) {
  t <- 0
  X <- c_0
  ruin <- FALSE
  while (t <= time_horizon) {
    waiting_time <- rexp(1, lambda_N)
    xi <- rexp(1, lambda_claims)
    t <- t + waiting_time
    X <- X + c * waiting_time - xi
    if (X < 0) {
      ruin <- TRUE
      break  # break out of while-loop
    }
  }
  ruin_occurred[i] <- ruin
}

# Estimate probability
mean(ruin_occurred)
#> [1] 0.5107

Thus, with the initial parameters that were chosen here, the insurer has approximately a 50-50 chance to go bankrupt. Of course, the initial parameters were chosen arbitrarily. In practice, an insurer will have to use its historic records to figure out what a good choice of the initial parameters is.

5.2 Brownian Motion / Wiener Process

In financial mathematics, there is an infamous market model called Black-Scholes model. In its simplest form, this model assumes that the market consists of

• one risk-free asset that delivers a constant rate of return \(r\) and
• one risky stock whose price is described by a so-called geometric Brownian motion with constant drift \(\mu\) and volatility \(\sigma\).

Let us untangle the latter asset step by step. A Brownian motion or Wiener process is a stochastic process \(W = \{W(t), t \geq 0\}\) that fulfills

1. \(W(0) = 0\)
2. \(W\) has continuous trajectories
3. \(W\) has independent increments which implies that \(W(t + s) - W(t)\), \(s \geq 0\), is independent of the past values \(W(x)\), \(x \leq t\)
4. \(W(t + s) - W(t) \sim \mathcal{N}(0, s)\) for all \(s, t \geq 0\)

The Poisson process we looked at earlier has condition 1 and 3 in common with the Brownian motion. Of course, there is a lot more to know about the Brownian motion as it is one of the best studied stochastic processes. But as we are only trying to use these type of processes for simulation purposes, let us stick to just that.

To do so, it is worth noting that from the above conditions one can deduce that for all \(0 < s < t\) it holds that \(W(t) = W(s) + \sqrt{t - s} Z\) where \(Z\) is a standard normal random variable. Combined with \(W(0) = 0\) this allows us to iteratively simulate a trajectory of a Brownian motion on a grid \(t = 0, \text{dt}, 2\text{dt}, \ldots, T\) where \(\text{dt}\) decribes the finenness of the grid. This may look like this:

Now, a geometric Brownian motion \(S = \{ S(t), t \geq 0 \}\) with drift \(\mu\) and volatility \(\sigma\) can be computed via \[\begin{align*} S(t) = S(0) \exp \bigg\{ \bigg( \mu - \frac{\sigma^2}{2} \bigg)t + \sigma B(t) \bigg\} \end{align*}\] where \(B = \{ B(t), t \geq 0 \}\) is a Brownian motion and \(S(0)\) is some initial value of \(S\). Consequently, the above realization of a Brownian motion can be used to give the following realization of a geometric Brownian motion with drift \(\mu = 0.06\), volatility \(\sigma = 0.15\) and initial value \(S(0) = 1\).

mu <- 0.06
sigma <- 0.15
S_0 <- 1

# t is time scale of previous Brownian motion
# B is previous realization of Brownian motion
S <- S_0 * exp((mu - sigma^2 / 2) * t + sigma * B) 
tibble(
  t = t,
  S = S
) %>% 
  ggplot() +
  geom_line(aes(t, S))

In the Black-Scholes model, this now may be considered as the stock price evolution over time. One reason why this model is so popular is that it enables us to compute the value of a so-called European call option on the stock \(S\).

A European call option is a financial contract between two parties that gives the buyer of said option the opportunity to buy a stock share at a so-called strike price \(K\) at a predetermined time \(T\) in the future.

On one hand, the buyer of a call option can make a profit if the strike price \(K\) is lower than the stock price at the time of execution, i.e. if \(K < S(T)\). On the other hand, if \(K > S(T)\) then the buyer will simply not use the option (because why would anyone pay more than necessary at the regular market). Thus, the payout of a call option at time \(T\) can be summarized as \(C = \max\{S(T) - K, 0\}\).

Of course, this begs the question: What is the value of such an option, i.e. what would be a good price to pay for such an option? In the Black-Scholes model there is a theoretical answer to this. At any time \(t \in [0, T]\) this answer depends on the

• time to maturity \(T - t\),
• current price of the risky asset \(S(t)\),
• strike price \(K\),
• volatility \(\sigma\) of the risky asset and
• the risk-free return \(r\) from the risk-free asset in the market.

Given these quantities, the Black-Scholes formula delivers the value of an European call option \(C_t\) at time \(t\) by \[\begin{align*} C_t = \Phi(d_1)S_t - \Phi(d_2) K e^{-r(T - t)}, \end{align*}\] where \(\Phi\) is the distribution function of the standard normal distribution and \[\begin{align*} d_1 &= \frac{1}{\sigma \sqrt{T - t}} \bigg[ \log \frac{S_t}{K} + \bigg( r + \frac{\sigma^2}{2} \bigg) (T - t) \bigg]\\ d_2 &= d_1 - \sigma\sqrt{T-t} \end{align*}\]

Let us use these (maybe complicated looking) formulas as a base line to determine the value of a European call option Monte Carlo-style. This can be accomplished as follows:

1. Simulate a lot of paths of the risky asset (under the so-called risk-neutral measure which here simply means that \(\mu = r\))
2. Compute the payoff of the option for each path
3. Average the payoffs
4. Since the option pays in the future we need to compute their value at time \(t\) via discounting the average by multiplying it by \(e^{-r(T-t)}\).

In order for step 1 to work, I have implemented a function generate_geom_BM(mu, sigma, S_0, time_horizon, dt = 0.025) that simulates a path from 0 to time_horizon in steps of length dt = 0.025 with parameters mu, sigma and S_0 (you will get to implement that in the exercises). The output of said function is as follows

mu <- 0.06
sigma <- 0.15
S_0 <- 1
time_horizon <- 10
generate_geom_BM(mu, sigma, S_0, time_horizon)
#> # A tibble: 401 x 2
#>       t     S
#>   <dbl> <dbl>
#> 1 0     1    
#> 2 0.025 0.989
#> 3 0.05  0.996
#> 4 0.075 0.984
#> 5 0.1   1.01 
#> 6 0.125 1.02 
#> # ... with 395 more rows

Now we can follow the above steps and compute the call option price at time \(t = 0\) and compare it to the theoretical value.

set.seed(123)
sigma <- 0.15
S_0 <- 1
time_horizon <- 10
r <- 0.02
K <- 1
N <- 10000

# Theoretical price
d_1 <- 1 / (sigma * sqrt(time_horizon)) * (
  log(S_0 / K) + (r + sigma^2 / 2) * time_horizon
)
d_2 <- d_1 - sigma * sqrt(time_horizon)
val_theo <- pnorm(d_1) * S_0 - pnorm(d_2) * K * exp(-r * time_horizon)

# Simulated price
payoff <- numeric(N)
for (i in seq_along(payoff)) {
  # Choose mu = r here for "risk-neutral measure"
  S <- generate_geom_BM(r, sigma, S_0, time_horizon)$S
  payoff[i] <- max(S[length(S)] - K, 0)
}
val_simu <- exp(-r * time_horizon) *  mean(payoff)

# Output
glue::glue("Theoretical value: {val_theo}")
#> Theoretical value: 0.275713492487472
glue::glue("Monte-Carlo value: {val_simu}")
#> Monte-Carlo value: 0.271898095322233

As you can see, the Monte-Carlo approach delivers a result which is very close to the theoretical result. Again, the Monte-Carlo approach is not really “necessary” here since there already is a formula to price the call option. However, Monto-Carlo option pricing is hugely popular among practitioners as it allows to price even more exotic financial derivatives when there is no analytic solution.

5.3 Exercises

5.3.1 The Poisson Process

Write a function simulate_PP(lambda, time, grid) which simulates a Poisson process with rate lambda from time zero to time on a grid with fineness grid and returns a tibble consisting of the columns t and X. Here, t contains all grid points from zero to time and X denotes the value of the Poisson process at time t. For instance, the result may look like this.

set.seed(1234)
simulate_PP(lambda = 0.5, time = 5, grid = 0.5)
#> # A tibble: 11 x 2
#>       t     X
#>   <dbl> <dbl>
#> 1   0       1
#> 2   0.5     1
#> 3   1       1
#> 4   1.5     1
#> 5   2       1
#> 6   2.5     1
#> # ... with 5 more rows

Use it to compute multiple box plots summarizing the distributions of \(X(n)\), \(n \in \{1, \ldots, 10\}\), where \(X\) is a Poisson process with rate \(\lambda = 1\). Also, indicate the theoretical expected value \(\mathbb{E}X(n) = n\) via a red dot. Your plots should be based on 1000 realizations of the Poisson process and should look something like this.

Hint: You may collect the simulated processes in a tibble via iteratively extending a tibble in a for loop. Even though this is not efficient, we have not covered lists in great detail yet which are a convenient format for such purposes. So, the slow method is fine, too.