Algorithm to Estimate the Probability That aCasuality Insurance Company’s Capital Will Remain Positive Throughout an Accounting Period

We want to do a simulation run to estimate the probability that a casuality insurance company’s capital will remain positive throughout the accounting period [0, T]. We will assume that the initial capital is a0, and that the company initially has n0 policyholders. We assume that policyholders generate claims (experience losses) according to independent Poisson processes with common rate λ. Individual claim amounts are continuous r.v.’s having a severity distribution with c.d.f. FX. Suppose also that new customers buy a policy according to a Poisson process with rate ν. Suppose that each policyholder remains with the company for an exponentially distributed time with rate µ. Also, suppose that each policyholder pays the firm at a fixed rate c per unit time.

Starting with n0 customers at time t = 0, and with an initial capital a0, we want to simulate the process as a way of estimating the probability that the company will be able to meet its obligations (claims payments) throughout the time period [0, T].

Possible events that could happen are: a) a new policyholder is acquired, b) an existing policyholder is lost, or c) a policyholder submits a claim (experiences a loss).

The system state variables are: n = number of policyholders at time t, a = current capital at time t. We could also consider S = total amount of claims paid by time t.

If (n, a) is the state of the system at time t, then the next event will occur at time t + X, where X ~ Exponential(ν + nµ + nλ). Regardless of when the event occurs, it will be

i)A new policyholder with probability , producing state (n + 1, a) at time

t+X, or

ii)A lost policyholder with probability , producing state (n - 1, a) at time

t+X, or

iii)A claim of (random) amount Y with probability , producing state (n, a-Y) at time t + X.

The indicator of the state of the system at the end of the period will be: if the company’s capital is positive at time T, or if not.

There are many possible loss distributions that have been used in actuarial applications (lognormal, Weibull, gamma, Pareto, etc.) In this example, for simplicity, it will be assumed that the loss distribution is an exponential distribution having a mean of β. Since we want to estimate the probability that Indic = 1, we will do the simulation a large number, K, of times, and calculate the fraction of times that Indic = 1 at the end of a simulation.

The algorithm is as follows:

Step 1: Declare all variables

Step 2: Input initial number of policyholders, initial amount of capital, rate of payment per unit time by each policyholder, Poisson rate of generation of claims, Poisson rate of acquisition of new policyholders, Poisson rate of loss of existing policyholders, parameter(s) of the severity (loss) distribution, sample size K, and initial seed value.

Step 3: Initialize Prop = 0.

Step 4: Do loop generating sample of runs of process, for index value L = 1 to K.

A)Call function RINSURE to generate the Lth value of INDIC.

B)Increment Prop = Prop + (INDIC(L))/K.

Step 5: Print: ‘The probability that the firm’s assets will still be positive at time T is approximately’, Proportion.

Step 6: Append RINSURE function, Exponential function, (Loss function), and RAN2 function.

The algorithm for the RINSURE function is as follows:

Step 1: Set Indic = 1, t = 0, , .

Step 2: Call Exponential function with mean ν + nµ + nλ to generate a new event time tE= X.

Step 3: If tE > T, Return.

Step 4: If t E ≤ T, then

A)a = a + nc(tE – t)

B)t = tE

C)Call RAN2 function to generate a unit uniform random deviate U.

i) If , then n = n + 1

ii)Else if, then n = n - 1

iii)Else .

i)Call Loss distribution function to generate claim amount Y.

ii)If Y > a, then set Indic = 0 and Return.

Else set a = a – Y.

Step 5: Call Exponential function with mean ν + nµ + nλ to generate a new event time X, set

and go to Step 3.

The FORTRAN program to implement the simulation is as follows:

PROGRAM INSURANCE

INTEGER N0, ISEED, K

REAL MU, NU, LAMBDA, TMAX, A0, C, BETA

REAL RK, INDIC(20000), PROP

PRINT*, 'Input initial capital.'

READ*, A0

PRINT*, 'Input length of accounting period.'

READ*, TMAX

PRINT*, 'Input initial pool size.'

READ*, N0

PRINT*, 'Input policy income (premium) rate (per customer).'

READ*, C

PRINT*, 'Input customer acquisition rate.'

READ*, NU

PRINT*, 'Input customer loss rate.'

READ*, MU

PRINT*, 'Input claim rate.'

READ*, LAMBDA

PRINT*, 'Input Loss distribution parameter(s).'

READ*, BETA

PRINT*, 'Input initial seed value.'

READ*, ISEED

PRINT*, 'Input sample size.'

READ*, K

RK = K

DO L = 1, K

INDIC(L) = 1.

ENDDO

PROP = 0.

DO 20 L = 1, K

INDIC(L) = RINSURE(A0,N0,C,NU,MU,LAMBDA,BETA,TMAX,ISEED)

PROP = PROP + (INDIC(L)/RK)

20 CONTINUE

PRINT*, 'The estimated probability is ', PROP, '.'

END

REAL FUNCTION RINSURE(A0,N0,C,NU,MU,LAMBDA,BETA,TMAX,ISEED)

REAL A, EXPPAR, P1, P2, X, TE, T, Y, U, LAMBDA, MU, NU

INTEGER N, ISEED

RINSURE = 1.

T = 0.

A = A0

N = N0

EXPPAR = NU + (N*MU) + (N*LAMBDA)

P1 = NU/EXPPAR

P2 = (NU + (N*MU))/EXPPAR

X = EXPON(ISEED, EXPPAR)

TE = X

10 IF (TE .GT. TMAX) THEN

RETURN

ELSE IF (TE .LE. TMAX) THEN

A = A + (N*C*(TE - T))

T = TE

U = RAN2(ISEED)

IF (U .LT. NU/EXPPAR) THEN

N = N + 1

ELSE IF (P1 .LE. U .AND. U .LT. P2) THEN

N = N - 1

ELSE

Y = EXPON(ISEED,1./BETA)

IF (Y .GT. A) THEN

RINSURE = 0.

RETURN

ELSE IF (Y .LE. A) THEN

A = A - Y

END IF

END IF

EXPPAR = NU + (N*MU) + (N*LAMBDA)

P1 = NU/EXPPAR

P2 = (NU + (N*MU))/EXPPAR

X = EXPON(ISEED,EXPPAR)

TE = T + X

GO TO 10

END IF

RETURN

END

REAL FUNCTION EXPON(ISEED,RATE)

U = RAN2(ISEED)

EXPON = -LOG(U)/RATE

RETURN

END

REAL FUNCTION RAN2(IDUM)

PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6)

DATA IFF /0/

DIMENSION IR(97)

IY = 0

IF(IDUM .LT. 0. .OR. IFF .EQ. 0)THEN

IFF=1

IDUM=MOD(IC-IDUM,M)

DO 11 J=1,97

IDUM=MOD(IA*IDUM+IC,M)

IR(J)=IDUM

11 CONTINUE

IDUM=MOD(IA*IDUM+IC,M)

IY=IDUM

ENDIF

J=1+(97*IY)/M

IF(J .GT. 97 .OR. J .LT. 1)PAUSE

IY=IR(J)

RAN2=IY*RM

IDUM=MOD(IA*IDUM+IC,M)

IR(J)=IDUM

RETURN

END