Ruin Probability via Stochastic Process Simulation [Julia] vs [Python]

Problem Introduction and the Julia Language

I have recently been using the Julia programming language for simulation problems or problems that require heavy computational processing. I even used it in my undergraduate thesis, which was a great way to boost my learning of this new tool.

I chose it because of its ability to solve the two-language problem, by allowing (i) rapid prototyping of functional code and (ii) achieving excellent computational performance, all without the need to rewrite the algorithm in a compiled language like C/C++.

So I want to share a case where using Julia was computationally advantageous and that many actuarial science students will likely encounter during their studies: calculating the ruin probability of an insurance company through Cramér-Lundberg simulations via Monte Carlo simulation.

If you want to check all the implemented code without going through the entire step-by-step below, here is the repository on github.

The algorithm simulates an insurance company's net worth over a discrete time period of 10 years (120 months), aiming to obtain the company's ruin probability by calculating the ratio between simulated cases that resulted in ruin and the total number of simulated cases. Here I will cover the classic case, but I may bring more advanced variations that incorporate additional relevant variables in the future.

I will also present the same steps in Julia and Python so that the reader can reproduce and compare the processing time of both on their personal machine.

Problem Assumptions

• 120 months (10 years) of the company's net worth will be evaluated
• A constant premium of 500 reais per client is assumed
• New client acquisition follows a \(Poisson(\lambda)\) distribution with mean 50. (\(n_p\))
• A 5% churn rate is incorporated, considered before new acquisitions
• For the Aggregate Claim calculation, \(S_t\):
  • The number of claims \(n_t\) follows a \(Poisson(\lambda)\) distribution with mean equal to 10% of the client base that month
  • Each claim, \(X_t\), follows an \(Exp(\alpha)\) distribution with mean 5000.

Number of policyholders with 5% churn + new clients

Considering an initial number of clients (\(n_0\)) and the assumptions of 5% churn of current policyholders plus the addition of new policyholders (\(n_p\)) following a \(Poisson(\lambda)\) distribution, we can derive the following steps for time t=1 and t=2.

\[ n_1 = 0,95 \cdot n_0 + n_p \\ n_2 = 0,95 \cdot n_1 + n_p \\ n_2 = 0,95(0,95 \cdot n_0 + n_p) + n_p \\ n_2 = 0,95^2 \cdot n_0 + 0,95 \cdot n_p + n_p \\ n_2 = 0,95^2 \cdot n_0 + (0,95 + 1)n_p \]

When we generalize this function in t, we get:

\[n_t = 0,95^tn_0 + (0,95^{t-1} + ... + 0,95 + 1)n_p\]

Since this summation is a geometric series, we can obtain the following equivalences \((0,95^{t-1} + ... + 0,95 + 1) = \frac{1-0,95^t}{1-0,95} = \frac{1-0,95^t}{0,05} = 20 * (1 - 0,95^t)\), thus resulting in the final equation applied in the algorithm.

\[n_t = 0,95^tn_0 + 20(1-0,95^t)n_p\]

Methodology

The classic model of Cramer-Lundberg is a stochastic process that simulates the behavior of net worth through the following discrete formula:

\[U_t = U_{t-1} + P_t - S_t, \;\;\; t \ge 1\]

where \(t = 1,...,n\)

We assume an initial net worth value (\(U_0\)) at the start of time \(t\), which will be a scenario study parameter. Next, we assume \(P_t = c \cdot n_t\), where the premium value for the period is the premium constant times the number of clients in that period, according to the adopted assumption.

Finally, we move to the calculation of the Aggregate Claim or Convoluted Claim (\(S_t\)), which is basically the convolution of the following two random variables:

• \(N \sim Poisson(\alpha = 0.1*n_t)\), which simulates the number of claims that occurred.
• \(X \sim Exp(\lambda = 1/5000)\), which simulates the severity of the claims.

Thus, the Aggregate Claim is given by the following formula:

\[S_t = \sum_{i=1}^{n_t} X_i\] where the number of claims (\(n_t\)) is first simulated through a Poisson distribution, the claim values (\(X_i\)) are drawn from an Exponential distribution, and finally all values are summed, resulting in the period's Aggregate Claim (\(S_t\)).

With this, we can simulate, for example, 5000 trajectories of the net worth (\(U_t\)) and check how many cases resulted in company ruin (\(U_t < 0\)). Thus, the proportion of ruin cases to total simulated cases equals the insurance company's ruin probability, given the adopted assumptions \(\psi(U_0) = \frac{N_{ruins}}{N_{simulations}}\).

Coding the Algorithm

First, we import the necessary packages for both languages.

Julia

using Distributions

Python

import numpy as np
from numpy.random import exponential as expon
from numpy.random import poisson as poiss

function novos_clientes(n₀, poiss, t)
    C = n₀*0.95^t + 20*(1 - 0.95^t)*rand(poiss)
    round(Int, C)
end

def novos_clientes(N, lbds, i):
    clientes = N*0.95**(i-1) + 20*(1 -(0.95)**(i-1))*poiss(lbds)
    return int(clientes)

function sin_convol(expn, cliente)
    N = rand(Poisson(0.1 * cliente))
    sum(rand(expn, N))
end

def sin_convoluto(expn, cliente):
    N = poiss(0.1 * cliente)
    S = expon(expn, N)
    return np.sum(S)

function testa_ruina(k, N₀, λ, α, U₀, c)

    # vetor do patrimônio líquido
    U = Vector{Float64}(undef, k)

    # gera os tipos de distribuições do modelo
    poiss = Poisson(λ)
    expn  = Exponential(α)

    # executa o primeiro mês para considerar o argumento U₀
    ncliente = novos_clientes(N₀, poiss, 1)
    P₁ = c * ncliente
    S₁ = sin_convol(expn, ncliente)
    U[1] = U₀ + P₁ - S₁

    # executa os demais caminhos para verificar se houve ruína
    R = 0
    for i in 2:k
        ncliente = novos_clientes(N₀, poiss, i)
        Pᵢ = c * ncliente
        Sᵢ = sin_convol(expn, ncliente)
        U[i] = U[i-1] + Pᵢ - Sᵢ
        if U[i] < 0
            R = 1
            break
        end
    end

    return R
end

def testa_ruina(K, N0, lbds, alpha, U0, C):
    U = np.array([])

    ncliente = novos_clientes(N0, lbds, 1)
    P1 = C * ncliente
    S1 = sin_convoluto(alpha, ncliente)
    U = np.append(U, U0)
    U[0] = U[0] + P1 + S1

    R = 0
    for i in range(1, K):
        ncliente = novos_clientes(N0, lbds, i)
        P = C * ncliente
        S = sin_convoluto(alpha, ncliente)
        new_U = U[i-1] + P - S
        U = np.append(U, new_U)
        if U[i] < 0:
            R = 1
            break

    return R

function prob_ruina(n_sim; k=120, n₀=50, λ=50.0, α=5000.0, U₀=10000.0, c=500.0)
    R = Vector{Float64}(undef, n_sim)
    for i in 1:n_sim
        R[i] = testa_ruina(k, n₀, λ, α, U₀, c)
    end
    return sum(R)/n_sim
end

def prob_ruina(NSIM, K=120, N0=50, lbds=50, alpha=5000, U0=10000, C=500):
    ruinas = np.array([])
    for _ in range(NSIM):
        ruinas = np.append(ruinas, testa_ruina(K, N0, lbds, alpha, U0, C))
    return np.sum(ruinas) / NSIM

using BenchmarkTools

tempo_julia = @belapsed prob_ruina(5000)

> 0.0551965

from timeit import default_timer as timer
from datetime import timedelta
import statistics

t = []
for _ in range(30):
    t1 = timer()
    _ = prob_ruina(5000)
    t2 = timer()
    t.append(t2-t1)
    
tempo_python = statistics.mean(t)

print(tempo_python, "seconds")

> 5.638724613333332 seconds

# razao do tempo de execução (python/julia)
razao_tempo = py$tempo_python / julia_eval("tempo_julia")
razao_tempo

> [1] 102.1573

Evaluating Different Scenarios

Aiming to understand the impact of initial capital and the initial number of policyholders on the ruin probability, I generated 24 different scenarios considering values of 50, 80, and 100 policyholders, and for initial capital ranging from 0 to 1 million. Below are all 24 scenarios, which are all combinations of the two variables.

capital = Float64[0, 10000, 30000, 50000, 70000, 100000, 500000, 1000000];
nseg = [50, 80, 100];
combin_prod = Base.product(nseg, capital) |> collect;
combin = [combin_prod...]

24-element Vector{Tuple{Int64, Float64}}:
 (50, 0.0)
 (80, 0.0)
 (100, 0.0)
 (50, 10000.0)
 (80, 10000.0)
 (100, 10000.0)
 (50, 30000.0)
 (80, 30000.0)
 (100, 30000.0)
 (50, 50000.0)
 ⋮
 (50, 100000.0)
 (80, 100000.0)
 (100, 100000.0)
 (50, 500000.0)
 (80, 500000.0)
 (100, 500000.0)
 (50, 1.0e6)
 (80, 1.0e6)
 (100, 1.0e6)

To run the 5000 simulations across the 24 different scenarios, I pass the scenario vector to the map() function to apply the simulation function to each scenario.

k = 5000;

probs_ruina = map(combin) do x
    prob_ruina(k, n₀ = x[1], U₀ = x[2]);
end;

results =[vcat(collect.(combin)'...) probs_ruina];

Scenario Analysis

With the results, we can generate the base table with the scenario outcomes.

using DataFrames, PrettyTables

df_ruina = DataFrame(nseg = results[:,1],
                    cap_ini = results[:,2],
                    prob_ruina = 100 * results[:,3]);

pretty_table(df_ruina,
            nosubheader=true,
            title = "\nResultado das 5 mil simulacoes";
            formatters = (ft_printf("%6.0f", 2), 
                          ft_printf("%5.2f", 3)) )

Resultado das 5 mil simulacoes
┌──────┬─────────┬────────────┐
│ nseg │ cap_ini │ prob_ruina │
├──────┼─────────┼────────────┤
│   50 │       0 │      97.00 │
│   80 │       0 │      96.62 │
│  100 │       0 │      95.64 │
│   50 │   10000 │      95.46 │
│   80 │   10000 │      95.58 │
│  100 │   10000 │      95.42 │
│   50 │   30000 │      93.98 │
│   80 │   30000 │      93.36 │
│  100 │   30000 │      93.78 │
│   50 │   50000 │      91.32 │
│   80 │   50000 │      90.90 │
│  100 │   50000 │      91.40 │
│   50 │   70000 │      88.10 │
│   80 │   70000 │      89.02 │
│  100 │   70000 │      89.44 │
│   50 │  100000 │      84.34 │
│   80 │  100000 │      85.58 │
│  100 │  100000 │      84.64 │
│   50 │  500000 │      44.30 │
│   80 │  500000 │      44.36 │
│  100 │  500000 │      44.94 │
│   50 │ 1000000 │      14.18 │
│   80 │ 1000000 │      14.74 │
│  100 │ 1000000 │      14.78 │
└──────┴─────────┴────────────┘

Increasing the Number of Simulations

Since Julia allows simulating a larger number of scenarios than Python, I tested running 5000, 50000, and 100000 simulations to evaluate the convergence of results and generated a plot with R that summarizes the scenario analysis and the number of simulations well.

The chart clearly shows the inverse relationship between initial capital and number of policyholders with the insurance company's ruin probability. The higher the initial capital or the initial number of policyholders, the lower the chance of the company going bankrupt. Furthermore, with 100,000 simulations, good convergence is observed with little variance in the results, thus eliminating the need to increase the number of simulations to obtain satisfactory results.

Final Julia Code

Below is an image with all the Julia code implemented in a summary image created on the Carbon platform.

In the following repository on github you can find more detailed code, including the R script I used to generate the final plot.