Plotting ruin in R

Question

I'm trying to recreate something similar to an image in modern actuarial risk theory using R: https://www.academia.edu/37238799/Modern_Actuarial_Risk_Theory (page 89)

Click here for image

In my case, the drops are of size based on an exponential distribution with parameter 1/2000 and they are spaced apart with Poisson inter arrival times which means they are distributed exponentially with a rate parameter of 0.25 (in my model)

The value of U is given by an initial surplus plus a premium income (c) per unit time (for an amount of time determined by the inter arrival distribution) minus a claim amount which would be random from the exponential distribution mentioned above.

I have a feeling a loop will need to be used and this is what I have so far:

lambda <- 0.25
EX <- 2000
theta <- 0.5

c <- lambda*EX*(1+theta)

x <- rexp(1, 1/2000)
s <- function(t1){for(t1 in 1:10){v <- c(rep(rexp(t1,1/2000)))
print(sum(v))}}


u <- function(t){10000+c*t}
plot(u, xlab = "t", xlim = c(-1,10), ylim = c(0,20000))
abline(v=0)


for(t1 in 1:10){v <- c(rep(rexp(t1,1/2000))) 
print(sum(v))}

The end goal is to run this simulation say 10,000 times over a 10 year span and use it as a visible representation as the rate of ruin for an insurance company.

Any help appreciated.

Answer 1

I think you're looking for something like this, all wrapped up in a neat function which by default draws the plot, but if wanted simply returns "ruin" or "safe" so you can run it in simulation:

simulate_ruin <- function(lambda = 0.25, EX = 2000,
                          theta  = 0.5, initial_amount = 10000,
                          max_time = 10, draw = TRUE) {

  income_per_year <- lambda * EX * (1 + theta)
  
  # Simulate a Poisson process. Include the initial time 0,
  # and replicate every other time point so we have values "before" and
  # "after" each drop
  times <- c(0, rep(cumsum(rexp(1000, lambda)), each = 2))
  times <- c(times[times < max_time], max_time)
  
  # This would be our income if there were no drops (a straight line)
  total_without_drops <- initial_amount + (income_per_year * times)
    
  # Now simulate some drops. 
  drop_size <- rexp((length(times) - 1) / 2, 1/2000)
  
  # Starting from times[3], we apply our cumulative drops every second index:
  payout_total <- rep(c(0, cumsum(drop_size)), each = 2)
  
  total <- total_without_drops - payout_total
  
  if(draw) {
    plot(times, total, type = "l", ylim = c(-1000, 20000))
    abline(h = 0, lty = 2)
  } else {
    if(any(total < 0)) 
      return("ruin") 
    else 
      return("safe")
  }
}

So we can call it once for a simulation:

simulate_ruin()

And again for a different simulation

simulate_ruin()

And table the results of 10,000 simulations to find the rate of ruin, which turns out to be around 3%

table(sapply(1:10000, function(x) simulate_ruin(draw = FALSE)))
#> 
#> ruin safe 
#>  305 9695

^{Created on 2022-04-06 by the reprex package (v2.0.1)}