Dividend Maximization in the Cramer-Lundberg Model using Homotopy Analysis Method

Problem statement: We used the Homotopy Analysis Method (HAM) to numerically compute the value function of the dividend payment in the basic insurance process. Approach: The process is a constant income stream from premiums which is subtracted a claim process of the Poisson type. Further, an allowance for payment of dividends to share holders was incorporated. Results: The case when the claims are exponential has an analytical solution. The HAM was then applied to the resulting Hamilton-Jacobi-Bellman equation and the numerical results obtained were compared to the theoretical results in order to check the validity of the method. Conclusion: The HAM was then applied to the model to check for other claim size distributions. The results obtained are very encouraging.


INTRODUCTION
The introduction of the classical collective risk model in by Lundberg (1903) to represent the surplus process of an insurance company sparked off lots of research on the probability of ruin of such a company for a long time and this approach was used to assess the stability of the company. It took not less than 50 years for De Finetti (1957) to criticize this approach and thereby laying the ground work for the study of dividend payouts since it was evident that a trajectory of the surplus process that did not lead to ruin in this model exceeded every finite level, which was typically unrealistic in practice.
De Finetti (1957) modeled the surplus of the insurance company using a discrete time surplus process with increments of ±1. He also assumed that the premium income per unit time was 1. He suggested that it was possible for a company to maximize the expected present value of all dividend payouts before ruin. He found an explicit solution for the value function and the optimal dividend strategy. He showed that the optimal strategy for his model was a barrier strategy.
The idea of dividend maximization later attracted much attention among researchers e.g., Gerber (1969;1974;1979). However the problem of dividend maximization in a classical continuous time model (Cramer-Lundberg model) was first studied by Buhlmann (1996). In section 6.4 of this text book, he discussed this problem for exponential claims and obtained explicit results. He showed that the optimal strategy was a barrier strategy.
Since then many research groups have tried to address this optimality question under more general and more realistic model assumptions and until nowadays this turns out to be a rich and challenging field of research that needs the combination of tools from analysis, probability, stochastic control and numerical methods. However, not many of these address the dividend optimization problem using numerical methods. Kasozi and Paulsen (2005a) used the block-byblock numerical method to study the flow of dividends under a constant force of interest. Their study culminated into a linear Volterra integro equation of the second kind. The block-by-block numerical method is fully developed in Paulsen et al. (2005). More pertinent literature is available in e.g., Paulsen (2003) and Paulsen and Gjessing (1997). Kasozi and Paulsen (2005b) also solved the problem of ultimate ruin probabilities in the classical model. They used an order four Block-by-block method in conjunction with Simpson's rule to solve the resulting integral equation.
A very candid review of the history and studies about dividend payouts was documented by Avanzi (2009).
In this study, we apply a numerical method called Homotopy Analysis Method to maximize dividend payouts prior to ruin. The method was proposed by Liao (1995) for his Ph.D dissertation. The method appears in his later works e.g., Liao (Liao, 1995;2004;.
In this method, the solution is considered as the summation of an infinite series, which usually converges rapidly to the exact solution. The HAM is based on homotopy, a fundamental concept in topology and differential geometry. The method relies on constructing a continuous mapping of an initial guess approximation to the exact solution of considered equations.
An auxiliary linear operator is chosen to construct such kind of continuous mapping and an auxiliary parameter is used to ensure the convergence of the series. The method enjoys great freedom in choosing initial approximations and auxiliary linear operators. The approximations obtained by the HAM are uniformly valid not only for small parameters, but also for very large parameters. Until recently, the application of the homotopy analysis method in nonlinear problems has attracted the attention of scientists and engineers.
This method provides a series solution of the linear integro differential equation. By use of an initial approximation, it is possible to construct a series solution of the exact solution of the IDE in question. This is achieved by choosing "well" an auxiliary parameter, an auxiliary function and an auxiliary Linear operator. This method can be used to solve both linear and non linear problems in science. It has been used by very many researchers to solve differential equations, partial differential equations, Integro differential equations and other types of equations that arise in scientific research. This method was used by Adawi et al. (2009) to numerically solve linear Volterra integral equations of the first and second kind and the results were very encouraging. In their study, they proposed theorems to show existence of the solution. They used MATLAB 7 to carry out the computations.
El-Nahhas (2007) used the HAM to numerically solve Volterra population equations. The results in that study were also highly acceptable. Like Adawi et al. (2009) andLundberg (1903), he also stated a theorem for existence of the series solution and he provided a detailed proof for it. The work in this study is well detailed and even shows how to get the auxiliary parameter. Jaradat et al. (2008), used the same method to provide a general way of numerically solving first order linear integro-differential equations. The validity of the method was checked using some examples in that study and the results were found to be very satisfactory.
A careful comparison of HAM with homotopy peturbation method was done by He (2004).
Most of the papers on HAM are handling first order Integro Differential Equations (IDEs) whose initial value is assumed to be known. However, in dividend maximization problems, the initial value of the resulting IDEs in the classical model is unknown. This poses a great challenge and it is another concern of this study.
The model: To give a mathematical formulation of the model, we assume that all the random variables and processes are defined on the stochastic basis (Ω,F,{ℱ} t ∈R + , P) here, Ω is a sample space. It is the set of all possible scenarios that can occur. The elements of this sample space will be denoted by ω, F is a σ-algebra on Ω. It is the collection of all events A⊆Ω{ℱ} t ∈R + is a filtration. That is an increasing and right continuous family of sub σ algebras of F The probability measure defined on F is P.
In risk theory, the Cramer-Lundberg model or the classical risk model is: Where: Y t = The surplus of the insurance company at time t y =Y 0 = The initial surplus (capital) p = The rate of premium income {N} t ∈R + = A Poisson process with rate λ (which is the counting process for the claims), the {X i } i ∈R + ∈N are independent and identically distributed (iid) random variables representing the claim sizes, with distribution function F. F has finite expectation and finite variance. N and X are independent.
Using Itˆo's formula, the infinitesimal generator for Y is given by the following integro-differential operator: A natural generalization of (1) is to incorporate de Finetti's ideas by allowing payment of dividends to the share holders. We predetermine a surplus level b in such a way that when the surplus level hits b, premium income is paid to shareholders as dividends. If the surplus is b and a claim occurs, the claim is paid and no dividends are paid until the surplus hits b again. Eventually and with certainty, ruin will occur at some stage in the future. When ruin occurs, the process stops. Question at hand is what b ensures that the insurer pays maximum dividends to shareholders before ruin? When dividends D(t) are distributed according to such a strategy and if D(t) is a non-decreasing process that represents the sum of dividends over a time interval (0,t], then (1) becomes: The dividend value function, which is the Expected Present Value (EPP) of the dividend payouts before time of ruin, is given by: The dividend value function then satisfies the so called Hamilton-Jacobi-Bellman (HJB) equation below: It should be noted that HJB equations are in general somehow hard to solve analytically. However attempts have been made to solve the above HJB equation when claims are exponential. That is F(s) = 1e -αs .
The IDE thus becomes: Converting the above IDE into a differential equation gives: where, m 1 <0 and m 2 >0 are the roots of the equation below: We now sort for a solution that n general works for all claim size distributions. This will be compared against results given by Eq. 9 to test its validity and there after applied for other claim size distributions.

MATERIALS AND METHOD
Homotopy involves defining a path by mapping points in the interval from 0-1 to points in the region in a continuous pattern; that is, so that the neighboring points on the interval correspond to neighboring points on the path.
Let us consider the problem of solving a general first order IDE below: Where: a = A constant; p(y), g(y) and K(y, t) and are given functions h(y) = The function to be determined We shall let N be an operator such that N[h(y)] = 0, h 0 (y) denote an initial approximation of the exact solution h(y); h ≠ 0 be an auxiliary parameter; H(y) ≠ 0 be an auxiliary function; and L be an auxiliary Linear operator such that L[h(y)] = 0 when h(x) = 0.
If we use q ∈ [0, 1] as an embedding parameter, we can construct that homotopy as; If we make the homotopy (12 ) equal to zero, we get the zero-order deformation: Since N[h(y)] = 0 and h ≠ 0, H(y) ≠ 0, we have: Notice that as q increases from 0-1, φ(y,q) varies continuously from the initial approximation h 0 (y) to the exact solution h(y). This is called deformation in homotopy. Using Taylor's theorem: To obtain the m th -order deformation equation, we differentiate the zero-order deformation (13) m times with respect to q and then divide by m! and then set q = 0. We make use of Leibnitz's rule for derivatives of products to get:

L[h (y) h (y)] hH(y)[h ' (y) p(y)
and the homotopy series solution is given by (18). To complete the iterative formula, we need to choose the Linear operator L, the auxiliary parameter h, the auxiliary function H(y) and the initial approximation h 0 (y).
To implement HAM to the IDE (11), we suppose the solution to Eq. 6 is of the form: Now (27) and (11) give: To find the value function, we use: The optimal barrier level b * , we solve the equation h"(b * ) = 0.
All calculations were done on a Pentium M computer using Maple 7. It should be noted that MATLAB or Mathematica could as well do the job very well.
The numerical solution for exponential claims are compared to the analytical solution in order to test the effectiveness of HAM. Figure 1 shows the h-curve for the problem and clearly the range of the auxiliary parameter h is (-2,-1). Let b N be the numerical optimal barrier level and b A be the exact optimal barrier.
Solving the iterative scheme (35) On working further to the more iterations, we notice that on solving for b there is convergence realized both in the numerically computed barrier b N and the value function V N . At the fourth iteration, the barrier approximates to 10.19323790. The value functions are shown in Table 2 and Fig. 2.

DISCUSSION
The HAM has been used to calculate for both the value function and the corresponding optimal barrier level for general claim size distributions. The method provides a series solution which converges to the exact solution after a few iterations. The method also requires less computational time. We have also labored to discuss the results as they fell by. We do recommend that this method can still be extended to second order IDE problems that arise when the basic insurance process is compounded by investments of the Black and Scholes type and general stochastic return on investments process.

CONCLUSION
We have derived value functions for both light tailed and heavy tailed distributions representing small and large claims using the HAM. This method can be used to find value functions for the diffusion model. The method gives results that are encouraging and it assures quick convergence to the exact solution. However, it should be noted that this method involves repeated analytical integration and thus only works with symbolic software.

ACKNOWLEDGEMENT
The researchers extend their gratitude to NOMA for the financial support and anonymous referees for the valuable comments.