1 Introduction↩︎

In actuarial sciences, the ultimate ruin probability is a key measure of the risk associated with the surplus process of an insurance company, and it has been extensively studied in various risk models. However, even for the classical risk model, the evaluation of the ruin probability is not an easy task (e.g. Rolski et al. [1], Asmussen and Albrecher [2]).

The paper concerns the classical risk model with surplus process \[\label{classical-risk-model} U(t) = u + ct-\sum_{i=1}^{N(t)}X_i, \quad t\geq 0,\tag{1}\] where \(u\) is the initial surplus, \(c>0\) is the rate of premium per unit time and \(S(t)=\sum_{i=1}^{N(t)}X_i\) denotes the aggregate claims up to time t with \(S(t)=0\) if \(N(t)=0\). More explicitly, \(\{X_1, X_2,...\}\) is a sequence of independent and identically distributed (i.i.d.) non negative random variables, which represents individual claim sizes, with distribution function \(P(x)=1-\overline{P}(x)=Pr(X\leq x)\) and mean \(\mathbb{E}[X]=\int_{0}^{\infty}zdP(z)\). We assume that the sizes of these claims arrive at an insurer according to a Poisson process \(\{N(t):t\geq 0\}\) with intensity \(\lambda\). We further assume that the claims are independent of the claim-arrival process.

Let \(T=inf\{t: U(t)<0\}\) denote the time of ruin, then the probability of ruin can be expressed as \[\psi(u)=Pr(T<\infty|U(0)=u), \quad u\geq 0.\] Moreover, the premium loading factor is given by \(c= (1+\theta)\lambda EX\), where \(\theta > 0\) is the security loading, ensuring that the model 1 is stable.

Define also \[\label{non-ruin32probability} \phi(u)=1-\psi(u)\tag{2}\] to be the probability that ruin never occurs starting from initial surplus \(u\), a probability also known as the probability of non-ruin. In general, explicit expressions for the (non-)ruin probability are known only in some cases, for instance, when claim sizes have exponential or phase-type distribution (see Asmussen and Albrecher [2], pg 14-15).

In recent years, a considerable number of contributions have focused on numerical methods for the computation of the ruin probability in the classical risk model. These papers roughly follow two main directions, namely attempting to solve the relevant integral or integro-differential equation numerically, and applying Monte Carlo estimation techniques. As for the second direction, several studies have concentrated on improving the Monte Carlo simulation techniques of ruin probability computation. Mnatsakanov et al. [3] proposed an empirical-type estimator of the Laplace transform of the ruin probability and subsequently recovered by a regularized Laplace inversion method. Within the same framework, Zhang et al. [4] studied the estimation of ruin probability in the classical ruin model, where estimator was constructed via Fourier inversion. Moreover, an estimation for ruin probability by calculating convolutions in the formula with quasi-Monte Carlo method proposed in [5]. For the estimation of ruin probability, many parametric, semi-paremetric and non-parametric estimators have been proposed in the literature; for references, see for example [6]–[8] among others.

On the other hand, a variety of numerical and semi-analytical techniques have been proposed in the ruin-theory literature. Considering the ruin probability evaluation in a classical risk model, a product-integration method to compute ruin probabilities was proposed by Ramsay and Usabel [9], while a Laplace transform based representation involving exponential integrals for Pareto claim sizes was derived by Ramsay [10]. Albrecher et al. [11] developed an efficient quadrature along with a rational-approximation procedure, motivated by the Trefethen–Weideman–Schmelzer method [12]. Moment-based Padé approximations of ultimate ruin probabilities via exponential mixtures are studied by Avram et. al. [13]. In the work of Goffard et al. [14], an orthogonal-polynomial expansion method was proposed to approximate the ultimate ruin probability in the compound Poisson ruin model. Constantinescu et al. [15] present three formulations for the ruin probabilities with gamma distributed claims. The results are obtained as inverse Laplace transforms of integro-differential equations with two involving Mittag–Leffler functions and the third based on moments of the claim size distribution. Fourier-cosine methods are introduced in [16] to analyze the sensitivity of ruin probabilities to perturbations in the moments of the claim-size and claim-arrival distributions. Gzyl et al. [17] propose a maximum entropy approach based on fractional exponential moments is employed to determine the probability of ruin. For further numerical approaches to ruin probability approximation in the classical risk model, see also [18]–[21].

In the present work, we investigate numerical approximations for ruin probabilities through the associated Volterra integro-differential equation (VIDE). The methods developed here are closely related to the recent literature on numerical Runge–Kutta (RK) schemes for VIDEs (see [22], [23]). In particular, we investigate high-order Runge–Kutta methods for the numerical computation of ruin probabilities in the classical risk model through two complementary approaches. The first combines the classical fourth-order one-step RK method with Simpson’s \(1/3\) quadrature rule, together with an equivalent ordinary differential equation formulation whenever possible. The second employs two-step Runge–Kutta (TSRK) methods coupled with Gaussian quadrature formulas for the integral term of the VIDE, while in the Pareto case we use a Gauss–Jacobi quadrature adapted to the associated heavy-tailed kernel structure.

The remainder of this paper is organized as follows. In Section 2, we present the framework of the one-step and two-step Runge–Kutta methods together with quadrature approximation of the convolution term. Section 3 develops the RK4 and TSRK schemes for Gamma claim-size distributions. In section 4, we consider the Pareto case using Gauss–Jacobi quadrature formulas. Section 5 provides numerical results that indicate the accuracy and efficiency of the proposed methods. Finally, Section 6 contains concluding remarks.

2 Runge–Kutta-based Methods for the integro-differential equation for ruin probability↩︎

This section describes the implementation of Runge–Kutta numerical methods for the approximation of the ruin probability in the classical risk model. Specifically, we consider combinations of Runge–Kutta schemes and numerical quadrature rules in order to solve the governing integro-differential equation.

Consider the equation satisfied by the survival probability \(\phi\) (see [1], Section 5.3): \[\phi'(u) = \frac{\lambda}{c}\phi(u) - \frac{\lambda}{c}\int_0^u \phi(u-z)\,dP(z), \label{eq:survival}\tag{3}\] where \(\phi'(u)\) denotes the derivative of \(\phi\) with respect to \(u\). Using the relation for \(\psi\) in \(\eqref{non-ruin32probability}\) and substituting into 3 , it follows that \[\psi'(u) = \frac{\lambda}{c}\psi(u) - \frac{\lambda}{c}I(u) - \frac{\lambda}{c}\overline{P}(u), \label{eq:ruin95compact}\tag{4}\] where \[I(u) = \int_0^u \psi(u-z)\,dP(z) = \int_0^u \psi(z)\,p(u-z)\,dz\] denotes the convolution term and \(p(x)\) is the density corresponding to \(P(x)\). The density formulation is used to avoid the reflected Stieltjes differential notation.

Hence, with the notation introduced above, equation \(\eqref{eq:ruin95compact}\) may be expressed as \[\mathsf{f}\bigl(u,\psi(u),I(u)\bigr) = \frac{\lambda}{c}\psi(u) - \frac{\lambda}{c}I(u) - \frac{\lambda}{c}\overline{P}(u),\] which places \(\psi'(u)=\mathsf{f}(u, \psi(u), I(u))\) within the general framework of Volterra integro-differential initial value problems: \[\left\{ \begin{align} y'(x) &= \mathsf{f}\left(x, y(x), \int_{x_0}^{x} K(x,s,y(s))\,ds\right), \quad x\in [x_0, T],\\ y(x_0) &= y_0. \end{align} \right. \label{eq:standard95vide}\tag{5}\] where \[\mathsf{f}:[x_0, T] \times \mathbb{R}\times \mathbb{R}, \quad K:\{(x,s):x_0\leq s \leq x \leq T\} \times \mathbb{R} \rightarrow \mathbb{R}\] are supposed to be sufficiently smooth to guarantee the existence and uniqueness of the solution (see Gripenberg et al. [24]).
Since the initial condition at \(u=0\) is explicitly available in the classical risk model, the associated VIDE of the form \(\eqref{eq:standard95vide}\) naturally suggests the application of Runge–Kutta methods together with suitable quadrature approximations for the integral term. In particular, fourth-order one-step and two-step Runge–Kutta discretizations are combined with Simpson’s \(1/3\) rule and Gaussian quadrature approximations, respectively.

2.1 Fourth-order Runge–Kutta (RK4) method with Simpson’s 1/3 rule↩︎

It is well known that \(\phi\) is a (proper) probability distribution with a mass of size \(1-\psi(0)\), where \[\label{initial32condition32psi40041} \psi(0)=\frac{1}{1+\theta}, \quad \theta>0.\tag{6}\] By equation \(\eqref{initial32condition32psi40041}\), we set the initial condition \(\psi(0)=\psi_0\). Let \[u_n=nh, \qquad \psi_n\approx \psi(u_n).\] For a stage displacement \(\delta\geq 0\), define \(g(z)=\psi(z)p(u_n+\delta-z)\). Thus, the history part of the convolution at the stage point \(u_n+\delta\) is \[\label{Hd} H_n(\delta) = \int_0^{u_n} g_{n,\delta}(z)\,dz,\qquad \delta \in \left\{0,\frac{h}{2},h\right\}.\tag{7}\] Now, using Simpson’s \(1/3\) rule, we approximate \[H_n(\delta) \approx \mathcal{S}^{(n)}_{1/3}[g_{n,\delta}].\] For \(n=0\), we set \(H_0(\delta)=0\). If \(n\) is even, then \[\mathcal{S}^{(n)}_{1/3}[\phi] = \frac{h}{3}\left[\phi(u_0)+\phi(u_n)+4\sum_{\substack{1\leq j\leq n-1\\ j\;\mathrm{odd}}}\phi(u_j)+2\sum_{\substack{2\leq j\leq n-2\\ j\;\mathrm{even}}} \phi(u_j)\right].\] If \(n\) is odd, Simpson’s \(1/3\) rule is applied on \([0,u_{n-1}]\) and the final interval is treated by the trapezoidal rule: \[\mathcal{S}^{(n)}_{1/3}[\phi]=\frac{h}{3} \left[\phi(u_0)+\phi(u_{n-1})+4\sum_{\substack{1\leq j\leq n-2\\ j\;\mathrm{odd}}}\phi(u_j)+2\sum_{\substack{2\leq j\leq n-3\\ j\;\mathrm{even}}}\phi(u_j)\right] +\frac{h}{2} \left[\phi(u_{n-1})+\phi(u_n)\right]. \label{eq:general-S13-odd}\tag{8}\] The first-stage convolution is \(I_1=H_n(0)\), and hence \[\label{k1-one-step-S1473} k_1=\frac{\lambda}{c}[\psi_n-I_1-\overline{P}(u_n)].\tag{9}\] For the second stage, setting \(Y_2=\psi_n+\frac{h}{2}k_1\), then the convolution at \(u_n+h/2\) is approximated by \[I_2=H_n\left(\frac{h}{2}\right)+\frac{h}{4} \left[\psi_n p\left(\frac{h}{2}\right)+Y_2 p(0)\right],\] and \[\label{k2-one-step-S1473} k_2=\frac{\lambda}{c}\left[Y_2-I_2-\overline{P}\left(u_n+\frac{h}{2}\right)\right].\tag{10}\] For the third stage, setting \(Y_3=\psi_n+\frac{h}{2}k_2\) and the corresponding convolution approximation is \[I_3=H_n\left(\frac{h}{2}\right)+\frac{h}{4} \left[\psi_n p\left(\frac{h}{2}\right)+Y_3 p(0) \right].\] Therefore, \[\label{k3-one-step-S1473} k_3=\frac{\lambda}{c} \left[Y_3-I_3-\overline{P}\left(u_n+\frac{h}{2}\right) \right].\tag{11}\] For the fourth stage, setting \(Y_4=\psi_n+h k_3\), then the convolution at \(u_n+h\) is approximated by \[I_4=H_n(h)+\frac{h}{6} \left[ \psi_n p(h)+4Y_3 p\left(\frac{h}{2}\right)+Y_4 p(0) \right].\] Hence, \[\label{k4-one-step-S1473} k_4=\frac{\lambda}{c} \left[Y_4-I_4-\overline{P}(u_n+h) \right].\tag{12}\] Finally, using the stage values defined in 9 – 12 , the solution is advanced by \[\label{psi95n431} \psi_{n+1}=\psi_n+\frac{h}{6}[k_1+2k_2+2k_3+k_4].\tag{13}\]

2.2 Two-Step Runge–Kutta (TSRK) methods with Gauss quadrature↩︎

A general \(\mathsf{m}\)-stage two-step Runge–Kutta discretization for 4 may be written in the form \[\psi_{n+1} = \theta_1\psi_n+\theta_2\psi_{n-1} + h\sum_{i=1}^{\mathsf{m}} v_i k_i^{[n-1]} + h\sum_{i=1}^{\mathsf{m}} w_i k_i^{[n]}, \quad n\geq 1, \label{eq:tsrk95ruin95update}\tag{14}\] where \[k_i^{[n]} = \mathsf{f}\bigl(u_n+c_i h,\Psi_i^{[n]},I_i^{[n]}\bigr), \quad i=1,\dots,\mathsf{m},\] and \[\Psi_i^{[n]} = \delta_{i1}\psi_n+\delta_{i2}\psi_{n-1} + h\sum_{j=1}^{\mathsf{m}} a_{ij}k_j^{[n-1]} + h\sum_{j=1}^{\mathsf{m}} b_{ij}k_j^{[n]}, \quad i=1,\dots,\mathsf{m}. \label{eq:tsrk95ruin95internal}\tag{15}\] where \(\Psi_i^{[n]} \approx \psi(u_n + c_i h).\) At the stage point \(u_n+c_i h\), the convolution term is decomposed into a history contribution over \([0,u_n]\) and a local contribution over \([u_n,u_n+c_i h]\). Accordingly, we define \[\begin{align} I_i^{[n]} &= I_n^{\mathrm{hist}}(u_n+c_i h) + h\sum_{j=1}^{\mu_0} w_{ij}\, p\bigl((c_i-d_{ij})h\bigr) \Bigl( \eta_{ij1}\psi_n+\eta_{ij2}\psi_{n-1} \notag\\ & +\, h\sum_{l=1}^{\mathsf{m}}\alpha_{ijl}k_l^{[n-1]} + h\sum_{l=1}^{\mathsf{m}}\beta_{ijl}k_l^{[n]} \Bigr), \quad i=1,\dots,\mathsf{m}, \label{eq:tsrk95ruin95integral95stage} \end{align}\tag{16}\] where \(I_n^{\mathrm{hist}}(u)\) denotes a quadrature approximation to the history integral \[I_n^{\mathrm{hist}}(u)\approx \int_0^{u_n}\psi(z)\,p(u-z)\,dz. \label{eq:tsrk95ruin95history}\tag{17}\] In direct analogy with the general TSRK construction, one may write \[\begin{align} I_n^{\mathrm{hist}}(u_n+c_i h) &= h\sum_{j=1}^{\mu_1}\omega_j\, p\bigl((n+c_i-\xi_j)h\bigr) \Bigl( \psi_0+\frac{h}{2}\sum_{l=1}^{2\mathsf{m}+1}\widetilde{\gamma}_{jl}k_l^{[0]} \Bigr) \notag\\ &+ h\sum_{\nu=1}^{n-1}\sum_{j=1}^{\mu_1}\omega_j\, p\bigl((n-\nu+c_i-\xi_j)h\bigr) \notag\\ & \times \Bigl( \zeta_{j1}\psi_\nu+\zeta_{j2}\psi_{\nu-1} + h\sum_{l=1}^{\mathsf{m}}\rho_{jl}k_l^{[\nu-1]} + h\sum_{l=1}^{\mathsf{m}}\gamma_{jl}k_l^{[\nu]} \Bigr). \label{eq:tsrk95ruin95history95formula} \end{align}\tag{18}\] Thus \(I_n^{\mathrm{hist}}(u_n+c_i h)\) approximates the accumulated contribution from the previously computed interval \([0,u_n]\), whereas the second term in 16 represents a local quadrature correction over the current interval \([u_n,u_n+c_i h]\).

Hence the resulting method advances the ruin probability by combining a two-step Runge–Kutta discretization of the differential term with a quadrature approximation of the convolution integral \(I(u)\).

Since the TSRK scheme is a two-step method, the values \(\psi_0\) and \(\psi_1\) must be available before the recurrence can be started; in our implementation, \(\psi_1\) is computed from the prescribed initial value \(\psi_0\) by one step of the classical fourth-order Runge–Kutta method.

In the rest of this subsection we implement a Fourth-order TSRK on the ruin probability equation. In particular, we apply the TSRK method with order 4 and \(\mathsf{m}=1\) Take \(\mu_0=\mu_1=2\), then the one-stage two-step Runge–Kutta method may be written as \[\begin{align} \psi_{n+1} &= \theta_1\psi_n+\theta_2\psi_{n-1} +h v_1 k_1^{[n-1]}+h w_1 k_1^{[n]}, \quad n\ge 1, \\ k_1^{[n]} &= \frac{\lambda}{c} \Bigl( \delta_{11}\psi_n+\delta_{12}\psi_{n-1} +h a_{11}k_1^{[n-1]}+h b_{11}k_1^{[n]} \Bigr) -\frac{\lambda}{c}\,I_1^{[n]} -\frac{\lambda}{c}\,\overline{P}(u_n+c_1h), \end{align}\] where \(I_1^{[n]}\), \(I_n^{\mathrm{hist}}(u_n+c_1h)\) are given by 16 and 18 for the specific choices of \(p\), \(\mathsf{m}\), \(\mu_0\), \(\mu_1\). The local order conditions corresponding to \(d_n^{(1)}=O(h^5)\) are \[\begin{align} 1-\theta_1-\theta_2 &= 0,\\ \frac{1}{\tau!} -\theta_2\frac{(-1)^\tau}{\tau!} -v_1\frac{(c_1-1)^{\tau-1}}{(\tau-1)!} -w_1\frac{c_1^{\,\tau-1}}{(\tau-1)!} &=0, \quad \tau=1,\dots,4. \end{align}\] The stage order conditions corresponding to \(d_{n,1}^{(2)}=O(h^4)\) are \[\begin{align} 1-\delta_{11}-\delta_{12} &= 0,\\ \frac{c_1^{\,\tau}}{\tau!} -\delta_{12}\frac{(-1)^\tau}{\tau!} -a_{11}\frac{(c_1-1)^{\tau-1}}{(\tau-1)!} -b_{11}\frac{c_1^{\,\tau-1}}{(\tau-1)!} &=0, \end{align}\]
where \(\tau=1,2,3 \, denotes the condition index in the order and interpolation conditions below\).
The quadrature conditions are \[\begin{align} \sum_{j=1}^{2}\omega_j\xi_j^k &= \frac{1}{k+1}, \quad k=0,1,2,3,\\ \sum_{j=1}^{2}w_{1j}d_{1j}^k &= \frac{c_1^{\,k+1}}{k+1}, \quad k=0,1,2,3. \end{align}\] The consistency conditions for the starting and history approximations are \[\begin{align} \frac{\xi_j^{\,\tau}}{\tau!} -\sum_{l=1}^{3}\widetilde{\gamma}_{jl}\frac{\widetilde{c}_l^{\,\tau-1}}{(\tau-1)!} &=0, \quad \tau=1,2,3, \\ 1-\zeta_{j1}-\zeta_{j2} &= 0, \\ \frac{\xi_j^{\,\tau}}{\tau!} -\frac{(-1)^\tau}{\tau!}\zeta_{j2} -\rho_{j1}\frac{(c_1-1)^{\tau-1}}{(\tau-1)!} -\gamma_{j1}\frac{c_1^{\,\tau-1}}{(\tau-1)!} &=0, \quad \tau=1,2,3, \\ 1-\eta_{1j1}-\eta_{1j2} &= 0, \\ \frac{\tau_{1j}^{\,\tau}}{\tau!} -\frac{(-1)^\tau}{\tau!}\eta_{1j2} -\alpha_{1j1}\frac{(c_1-1)^{\tau-1}}{(\tau-1)!} -\beta_{1j1}\frac{c_1^{\,\tau-1}}{(\tau-1)!} &=0, \quad \tau=1,2,3. \end{align}\]

Choosing \[\omega_1=\omega_2=\frac{1}{2}, \quad \xi_1=\frac{1}{2}-\frac{1}{2\sqrt3}, \quad \xi_2=\frac{1}{2}+\frac{1}{2\sqrt3},\] and \[w_{11}=\frac{c_1}{2}, \quad w_{12}=\frac{c_1}{2}, \quad d_{11}=\frac{c_1}{2}-\frac{c_1}{2\sqrt3}, \quad d_{12}=\frac{c_1}{2}+\frac{c_1}{2\sqrt3},\] we get the remaining coefficients by solving the above algebraic system.

3 Implementation of Runge–Kutta methods for Gamma claims↩︎

Gamma claim-size distributions have been considered in several ruin theory studies. For example, Constantinescu et al. [15] proposed three numerical procedures based on Mittag–Leffler functions, moments of the claim-size distribution and Laplace transform techniques for solving integro-differential equations. Since Gamma distributions with integer-valued shape parameter belong to phase-type distributions, explicit formulas for ruin probabilities are available, allowing us to assess the accuracy of our approximations (see Chapter 9 of [2]). In the following, the integro-differential equation of the ruin probability with Gamma claim-sizes is approximated through the RK4 and TSRK methods developed for VIDEs (see Wen et al. [22]).

3.1 Implementation of RK with Simpson’s 1/3 rule↩︎

Consider Gamma-distributed claim sizes with shape parameter \(2\) and rate parameter \(\beta>0\) with right tail \[\label{Gamma-tail} \overline{P}(u)=e^{-\beta u}(1+\beta u), \quad \beta>0.\tag{19}\] Since the mean claim size is \(2/\beta\), the safety loading relation implies \[\frac{\lambda}{c}=\frac{\beta}{2(1+\theta)}.\] In this case, the stage-dependent convolution also requires the auxiliary quantity \[J_n \approx \int_0^{u_n} q_n(z)\,dz, \label{eq:gamma-q}\tag{20}\] where \(q_n(z)=\psi(z)\beta^2 e^{-\beta(u_n-z)}\). Since, for an intermediate Runge–Kutta stage point \(u_n+\delta\), where \(\delta=h/2\) or \(\delta=h\), the history part of the convolution is \[H_n(\delta)=\int_0^{u_n} \psi(z)\beta^2(u_n+\delta-z) e^{-\beta(u_n+\delta-z)} \,dz. \label{RK4-Gamma-1}\tag{21}\] Also, equation \(\eqref{RK4-Gamma-1}\) can be rewritten as \[\begin{align} H_n(\delta) &= \int_0^{u_n} \psi(z)\beta^2\bigl((u_n-z)+\delta\bigr) e^{-\beta(u_n-z)}e^{-\beta\delta} \,dz \notag\\ &= e^{-\beta\delta} \left[ \int_0^{u_n} \psi(z)\beta^2(u_n-z)e^{-\beta(u_n-z)} \,dz + \delta \int_0^{u_n} \psi(z)\beta^2 e^{-\beta(u_n-z)} \,dz \right]. \end{align}\] The first integral is precisely the history convolution \(I_n\). The second integral motivates the definition of the auxiliary quantity \(J_n\) (see Equation \(\eqref{eq:gamma-q}\)).

Consequently, the shifted history contribution can be written as \[H_n(\delta) = e^{-\beta\delta} \left( I_n+\delta J_n \right).\] For \(n=0\), we set \(I_0=J_0=0\). For \(n\geq 1\), both \(I_n\) and \(J_n\) are computed by Simpson’s \(1/3\) history quadrature. Namely, \[I_n=\mathcal{S}^{(n)}_{1/3}[g_n], \qquad J_n=\mathcal{S}^{(n)}_{1/3}[q_n],\] where, if \(n\) is even, \[\mathcal{S}^{(n)}_{1/3}[\phi]=\frac{h}{3} \left[\phi(u_0)+\phi(u_n)+4\sum_{\substack{1\leq j\leq n-1\\ j\;\mathrm{odd}}}\phi(u_j)+ 2\sum_{\substack{2\leq j\leq n-2\\ j\;\mathrm{even}}} \phi(u_j) \right].\] If \(n\) is odd, Simpson’s \(1/3\) rule is applied on \([0,u_{n-1}]\) and the final interval is treated by the trapezoidal rule: \[\mathcal{S}^{(n)}_{1/3}[\phi]=\frac{h}{3} \left[\phi(u_0)+\phi(u_{n-1})+ 4\sum_{\substack{1\leq j\leq n-2\\ j\;\mathrm{odd}}}\phi(u_j)+ 2\sum_{\substack{2\leq j\leq n-3\\ j\;\mathrm{even}}} \phi(u_j)\right] +\frac{h}{2} [\phi(u_{n-1})+\phi(u_n)].\] The Runge–Kutta stages are then evaluated as follows. First, we define \[\label{k1-one-step-Gamma} k_1=\frac{\lambda}{c} \left[ \psi_n-I_n-e^{-\beta u_n}(1+\beta u_n) \right].\tag{22}\] Now, setting \(Y_2=\psi_n+\frac{h}{2}k_1\), then the convolution at \(u_n+h/2\) is approximated by \[I_2=e^{-\beta h/2}\left(I_n+\frac{h}{2}J_n \right)+\frac{\beta^2h^2}{8}\psi_n e^{-\beta h/2}.\] Hence, \[\label{k2-one-step-Gamma} k_2=\frac{\lambda}{c} \left[Y_2-I_2-e^{-\beta(u_n+h/2)} \left(1+\beta(u_n+h/2)\right) \right].\tag{23}\] For the third stage, we set \[Y_3=\psi_n+\frac{h}{2}k_2,\] and \[I_3=e^{-\beta h/2} \left(I_n+\frac{h}{2}J_n \right)+ \frac{\beta^2h^2}{8} \psi_n e^{-\beta h/2}.\] Thus, \[\label{k3-one-step-Gamma} k_3=\frac{\lambda}{c} \left[Y_3-I_3-e^{-\beta(u_n+h/2)} \left(1+\beta(u_n+h/2)\right) \right].\tag{24}\] Finally, setting \(Y_4=\psi_n+h k_3\), then the convolution at \(u_n+h\) is approximated by \[I_4=e^{-\beta h}\left( I_n+hJ_n \right)+\frac{h}{6} \left[\psi_n\beta^2 h e^{-\beta h}+4Y_3\beta^2\frac{h}{2}e^{-\beta h/2}\right],\] and the fourth stage is \[\label{k4-one-step-Gamma} k_4=\frac{\lambda}{c}[Y_4-I_4-e^{-\beta(u_n+h)} \left(1+\beta(u_n+h)\right) ].\tag{25}\] The numerical solution is advanced by the classical fourth-order Runge–Kutta formula 13 , using the stage values defined in 22 –25 .

3.2 Implementation of TSRK method for Gamma\((2,\beta)\) Claims↩︎

Now, we employ the TSRK method with order 4, \(\mathsf{m}=1\) according to Remark 2 for Gamma-distributed claims with right tail given by \(\eqref{Gamma-tail}\). In this case, the convolution term can be represented by introducing the auxiliary variable \[K_1(u)= \int_0^u \psi(z)e^{-\beta(u-z)}\,dz.\] Consequently, the scalar integro-differential equation is transformed into the following finite-dimensional system of ordinary differential equations: \[\begin{cases} \psi'(u) = \frac{\lambda}{c}(\psi(u)-I(u)- e^{-\beta u}(1+\beta u)), \\[0.4em] K_1'(u) = \psi(u)-\beta K_1(u), \\[0.4em] I'(u) = \beta^2 K_1(u)-\beta I(u). \end{cases}\] Equivalently, defining \[Y(u) = \begin{pmatrix} \psi(u)\\ K_1(u)\\ I(u) \end{pmatrix},\] the system can be written in the linear non-homogeneous form \[Y'(u)=MY(u)+g(u),\] where \[M= \begin{pmatrix} \frac{\lambda}{c} & 0 & -\frac{\lambda}{c}\\ 1 & -\beta & 0\\ 0 & \beta^2 & -\beta \end{pmatrix}, \qquad g(u) = \begin{pmatrix} -\frac{\lambda}{c} e^{-\beta u}(1+\beta u)\\ 0\\ 0 \end{pmatrix}.\] The initial condition is \[Y(0) = \begin{pmatrix} \psi(0)\\ K_1(0)\\ I(0) \end{pmatrix} = \begin{pmatrix} \dfrac{1}{1+\theta}\\ 0\\ 0 \end{pmatrix}.\] The numerical method used in the implementation applies the two-step Runge–Kutta method directly to this ODE system. For a uniform mesh \(u_n=nh\), let \(Y_n\approx Y(u_n)\). A general \(\mathsf{m}\)-stage TSRK method has the form \[Y_{n+1} = \theta_1Y_n+\theta_2Y_{n-1} +h\sum_{i=1}^\mathsf{m} v_i k_i^{[n-1]} +h\sum_{i=1}^\mathsf{m} w_i k_i^{[n]},\] where \[k_i^{[n]} = F(u_n+c_i h,Y_i^{[n]}), \qquad F(u,Y)=MY+g(u),\] and the internal stages satisfy \[Y_i^{[n]} = \delta_{i1}Y_n+\delta_{i2}Y_{n-1} +h\sum_{j=1}^\mathsf{m} a_{ij}k_j^{[n-1]} +h\sum_{j=1}^\mathsf{m} b_{ij}k_j^{[n]}.\] Since \(F(u,Y)=MY+g(u)\) is linear in \(Y\), the stage equations are linear algebraic systems. Therefore, at each step, the current stage values can be obtained by solving a finite-dimensional linear system.

For the fourth-order one-stage method (\(\mathsf{m}=1\)), the stage equation is \[Y_1^{[n]} = \delta_{11}Y_n+\delta_{12}Y_{n-1} +h a_{11}k_1^{[n-1]} +h b_{11}k_1^{[n]},\] with \[k_1^{[n]}=M Y_1^{[n]}+g(u_n+c_1h).\] Substituting the expression for \(k_1^{[n]}\) into the stage equation yields \[\bigl(\mathbb{I}_3-hb_{11}M\bigr)Y_1^{[n]} = \delta_{11}Y_n+\delta_{12}Y_{n-1} +h a_{11}k_1^{[n-1]} +h b_{11}g(u_n+c_1h),\] where \(\mathbb{I}_3\) is the \(3\times3\) identity matrix. Once \(Y_1^{[n]}\) is found, we compute \[k_1^{[n]}=M Y_1^{[n]}+g(u_n+c_1h),\] and update \[Y_{n+1} = \theta_1Y_n+\theta_2Y_{n-1} +h v_1 k_1^{[n-1]} +h w_1 k_1^{[n]}.\]

For the sixth-order two-stage method \((\mathsf{m=2})\), the two current stages satisfy \[Y_i^{[n]} = \delta_{i1}Y_n+\delta_{i2}Y_{n-1} +h\sum_{j=1}^2 a_{ij}k_j^{[n-1]} +h\sum_{j=1}^2 b_{ij}k_j^{[n]}, \qquad i=1,2,\] with \[\label{k95i-TSRK} k_i^{[n]}=M Y_i^{[n]}+g(u_n+c_i h), \quad i=1,2.\tag{26}\]

This leads to a block linear system for \(Y_1^{[n]}\) and \(Y_2^{[n]}\): \[\begin{pmatrix} \mathbb{I}_3-hb_{11}M & -hb_{12}M\\ -hb_{21}M & I_3-hb_{22}M \end{pmatrix} \begin{pmatrix} Y_1^{[n]}\\ Y_2^{[n]} \end{pmatrix} = \begin{pmatrix} R_1^{[n]}\\ R_2^{[n]} \end{pmatrix},\] where \[R_i^{[n]} = \delta_{i1}Y_n+\delta_{i2}Y_{n-1} +h\sum_{j=1}^2 a_{ij}k_j^{[n-1]} +h\sum_{j=1}^2 b_{ij}g(u_n+c_jh), \qquad i=1,2.\] After solving this block system, the stage derivatives are evaluated as \(\eqref{k95i-TSRK}\) and the new approximation is computed from \[Y_{n+1}=\theta_1Y_n+\theta_2Y_{n-1}+h\sum_{i=1}^2 v_i k_i^{[n-1]}+h\sum_{i=1}^2 w_i k_i^{[n]}.\] In the present ODE system formulation, the special Erlang structure of the Gamma\((2,\beta)\) density makes it possible to represent the convolution term exactly by the auxiliary variables \(K_1(u)\) and \(I(u)\). Thus, the quadrature-based convolution approximation is replaced by an equivalent finite-dimensional ODE system.

4 Implementation of Runge–Kutta methods for Pareto claims↩︎

Among heavy-tailed claim-size distributions, the Pareto distribution is of particular interest in actuarial science and ruin theory. Practical computation algorithms have been proposed in the literature for ruin probabilities under different versions of the Pareto distribution (see [9], [10], [25]). In this subsection, we implement two numerical approaches based on the RK4 method combined with Simpson’s \(1/3\) quadrature rule and a two-step Runge–Kutta (TSRK) method combining the Gauss–Jacobi quadrature for the computation of ruin probabilities with Pareto claim sizes.

4.1 Implementation of RK with Simpson’s 1/3 rule↩︎

A single Pareto claim-size distribution with support on the positive real axis and mean \(\mathbb{E}[X]=1\) is considered. Since Pareto distributions with integer parameter are commonly used in risk theory, the right tail is given by \[\overline{P}(u)=\left(\frac{m}{u+m}\right)^{m+1}, \quad u\geq 0 \quad \text{and} \quad m=1,2, \ldots \label{pareto-tail}\tag{27}\] with corresponding density \[\label{Pareto-density} p(u)=\frac{m+1}{m}\left( \frac{m}{u+m} \right)^{m+2}, \quad u\geq 0.\tag{28}\] Recall 7 8 , then for the Pareto–Lomax density the first-stage convolution is \(I_1=H_n(0)\) and \[\label{k1-Pareto} k_1=\frac{\lambda}{c} \left[ \psi_n-I_1-\left(\frac{m}{u_n+m}\right)^{m+1} \right].\tag{29}\] For the second stage, setting \(Y_2=\psi_n+\frac{h}{2}k_1\), the convolution at \(u_n+h/2\) is approximated by \[I_2=H_n\left(\frac{h}{2}\right)+\frac{h}{4} \left[ \psi_n p\left(\frac{h}{2}\right)+Y_2 p(0) \right],\] where \(p(0)\) and \(p(h/2)\) are given by \(\eqref{Pareto-density}\) at \(u=0\) and \(u=h/2\), respectively. Hence, \[k_2=\frac{\lambda}{c} \left[Y_2-I_2-\left( \frac{m}{u_n+h/2+m}\right)^{m+1} \right].\] For the third stage, setting \(Y_3=\psi_n+\frac{h}{2}k_2\), the corresponding convolution approximation is \[I_3=H_n\left(\frac{h}{2}\right)+\frac{h}{4} \left[ \psi_n p\left(\frac{h}{2}\right)+Y_3 p(0) \right].\] Thus, \[k_3=\frac{\lambda}{c}\left[Y_3-I_3-\left(\frac{m}{u_n+h/2+m}\right)^{m+1} \right].\] For the fourth stage, setting \(Y_4=\psi_n+h k_3\), then the convolution at \(u_n+h\) is approximated by \[I_4=H_n(h)+\frac{h}{6} \left[ \psi_n p(h)+4Y_3 p\left(\frac{h}{2}\right)+Y_4 p(0) \right],\] where \(p(h)\) is obtained from \(\eqref{Pareto-density}\) for \(u=h\). Therefore, \[\label{k4-Pareto} k_4=\frac{\lambda}{c} \left[ Y_4-I_4-\left( \frac{m}{u_n+h+m}\right)^{m+1} \right].\tag{30}\] Finally, from the stage values 29 –30 , the solution is obtained by the classical RK4 formula 13 .

4.2 Implementation of TSRK method with order 4 and one internal stage↩︎

We now specialize the TSRK discretization to Pareto claims. The TSRK update 14 and internal stage 15 are unchanged; only the convolution term is evaluated by a Pareto-adapted history–local quadrature.

At the stage point \[s_1^{[n]}=u_n+c_1 h,\] the stage derivative is computed as \[k_1^{[n]} = \frac{\lambda}{c}\Psi_1^{[n]} - \frac{\lambda}{c}I_1^{[n]} - \frac{\lambda}{c}\overline{P}\left(s_1^{[n]}\right), \label{eq:pareto-tsrk-stage-history-local}\tag{31}\] where \[I_1^{[n]} = \int_0^{s_1^{[n]}} \psi(z)p(s_1^{[n]}-z)\,dz .\] With the lag variable \(x=s_1^{[n]}-z\) and \(\delta_1=c_1 h\), this convolution is split as \[I_1^{[n]} = \int_{\delta_1}^{s_1^{[n]}} \psi\left(s_1^{[n]}-x\right)p(x)\,dx + \int_0^{\delta_1} \psi\left(s_1^{[n]}-x\right)p(x)\,dx . \label{eq:pareto-history-local-split}\tag{32}\] The first integral contains only previously computed solution values, whereas the second is local to the current TSRK step.

For Pareto claims we use the change of variables \(y=m/(x+m)\) and \(x=m(1-y)/y\). Thus, for \(0\leq a<b\), it follows that \[\int_a^b f(x)p(x)\,dx=(m+1)\int_{Y_b}^{Y_a} f\!\left(m\frac{1-y}{y}\right)y^m\,dy, \label{eq:pareto-truncated-weight-integral}\tag{33}\] where \[Y_a=\frac{m}{a+m},\quad Y_b=\frac{m}{b+ m}.\] We approximate the last integral by a \(q\)-point Gaussian rule for the truncated weight \(y^m\) on \([Y_b,Y_a]\): \[\int_a^b f(x)p(x)\,dx \simeq (m+1) \sum_{r=1}^q \widetilde{\omega}_r^{[a,b]} f\!\left(x_r^{[a,b]}\right), \label{eq:pareto-truncated-gaussian-rule}\tag{34}\] where \[x_r^{[a,b]}=m\frac{1-y_r^{[a,b]}}{y_r^{[a,b]}} .\]

The nodes and weights are obtained from the moments \[M_k^{[a,b]}=\int_{Y_b}^{Y_a} y^{k+m}\,dy =\frac{Y_a^{m+k+1}-Y_b^{m+k+1}}{m+k+1}, \qquad k=0,\ldots,2q-1 .\] Let us define \[H_0^{[a,b]}= \left(M_{j+k}^{[a,b]}\right)_{j,k=0}^{q-1}, \qquad H_1^{[a,b]}= \left(M_{j+k+1}^{[a,b]}\right)_{j,k=0}^{q-1}.\] Then \(y_r^{[a,b]}\) are the generalized eigenvalues of \[H_1^{[a,b]}v=y\,H_0^{[a,b]}v,\] and the weights satisfy \[\sum_{r=1}^q \widetilde{\omega}_r^{[a,b]} \left(y_r^{[a,b]}\right)^k=M_k^{[a,b]}, \qquad k=0,\ldots,q-1 .\]

Applying this rule to the local part of 32 , with \(a=0\) and \(b=\delta_1\), gives \[I_{1,\mathrm{loc}}^{[n]}\simeq(m+1) \sum_{r=1}^q \widetilde{\omega}_r^{[0,\delta_1]} \psi\!\left(s_1^{[n]}-x_r^{[0,\delta_1]}\right). \label{eq:pareto-tsrk-local-part}\tag{35}\]

For the history part, split \([0,u_n]\) into \([u_\ell,u_{\ell+1}]\), \(\ell=0,\ldots,n-1\). The corresponding lag interval is \[a_{\ell 1}^{[n]}=s_1^{[n]}-u_{\ell+1}, \qquad b_{\ell 1}^{[n]}=s_1^{[n]}-u_\ell .\] Hence, \[I_{1,\mathrm{hist}}^{[n]}\simeq(m+1) \sum_{\ell=0}^{n-1} \sum_{r=1}^q \widetilde{\omega}_{\ell 1 r}^{[n]} \psi\!\left(s_1^{[n]}-x_{\ell 1 r}^{[n]}\right), \label{eq:pareto-tsrk-history-part}\tag{36}\] where \[x_{\ell 1 r}^{[n]}=m\frac{1-y_{\ell 1 r}^{[n]}}{y_{\ell 1 r}^{[n]}}\] and \(y_{\ell 1 r}^{[n]}\), \(\widetilde{\omega}_{\ell 1 r}^{[n]}\) are the nodes and weights of the same Gaussian rule on \[\left[ \frac{m}{b_{\ell 1}^{[n]}+m}, \frac{m}{a_{\ell 1}^{[n]}+m} \right].\] Thus, \[I_1^{[n]} \simeq I_{1,\mathrm{hist}}^{[n]} + I_{1,\mathrm{loc}}^{[n]}. \label{eq:pareto-tsrk-total-convolution}\tag{37}\] The Pareto tail entering 31 is \[\overline{P}(s_1^{[n]}) = \left( \frac{m}{s_1^{[n]}+m} \right)^{m+1}.\] The values of \(\psi\) at off-grid points are evaluated by the cubic interpolation used in the TSRK construction: local points by the TSRK stage interpolant and history points from the already computed numerical solution. The resulting Pareto-weighted history–local scheme is denoted by TSRK4–PHL.

5 Numerical results↩︎

This section summarizes the findings of comparing different orders of Runge–Kutta (RK) methods combined with various numerical integration schemes to calculate the probabilities of ruin under both Gamma and Pareto claim size distributions.

We consider Gamma-distributed claim sizes with different parameter choices, summarized in Tables 1 and 2. In particular, Table 1 presents numerical approximations of the survival probability when claim sizes follow Gamma\((2, 2.4)\) with safety loading \(\theta=0.2\). The numerical approximations obtained using the fourth-order Runge–Kutta method combined with Simpson’s \(1/3\) quadrature rule (denoted by RK4–S\(1/3\)), together with the two-step Runge–Kutta method of order four, indicated by TSRK4–G, respectively. We compare our results with the three infinite-series formulas for ruin probabilities proposed in [15]. We also include the results obtained by Martire in [21], where the ruin probability is computed through a numerical solution of the associated Volterra integral equation based on a Lipschitz-continuity approach. The notations M-1, M-2, and M-3 correspond to Methods 1–3 in [15], associated with different discretization steps \(\Delta\), while AM denotes the results obtained using the method in [21].

Table 2 presents the corresponding numerical approximations of the ruin probability for a claim size follows Gamma\((2,1)\) with \(\theta=1.5\). The same numerical methods are employed on the interval \(u\in [0, 14.892]\) using the same step size \(h=0.0016\). In this case, the results are compared with approximated values obtained through the scaled Laplace transform inversion method of Mnatsakanov et al. in [18], as reported in [14]. The scaled Laplace transform technique has been presented in [26] and applied to the approximation of ruin probabilities in [18].

Table 3 presents numerical approximations of the ruin probability for Pareto claim sizes with parameter \(m=1\) and for several values of the safety loading parameter \(\theta\). The approximations are obtained using the classical fourth-order Runge–Kutta method combined with Simpson’s \(1/3\) quadrature rule (RK–S\(1/3\)), as described in Subsection \(\ref{sec:section-RK4-S1473-Pareto}\), together with the two-step Runge–Kutta method combined with Gauss–Jacobi quadrature (TSRK4-PHL) presented in Subsection \(\ref{sec:section-tsrk-pareto-history-local}\). The computations are performed with step-size \(h=0.01\). The numerical results are compared with exact known ruin probabilities for the Pareto case given in [10] (denoted by RAM), as well as with the approximations obtained by Martire in [21] (denoted by AM).

8pt

6 Conclusions↩︎

In this paper, one-step and two-step Runge–Kutta methods have been developed for approximating ruin probabilities through the associated integro-differential Volterra equation. The proposed framework provides a natural setting for the application of Runge–Kutta discretizations combined with suitable quadrature approximations for the convolution term.

For the special Gamma\((2, \beta)\) case, the convolution term admits an equivalent finite-dimensional ODE representation, allowing the direct implementation of TSRK methods on the resulting system of ordinary equations. For Pareto claims, the convolution term was treated numerically through adapted quadrature procedures based on Gauss–Jacobi formulas.

Numerical results indicate that the proposed methods provide accurate approximations for both light-tailed and heavy-tailed claim size distributions, remaining in close agreement with analytical formulas and benchmark numerical results in the literature. The comparisons presented in Tables 1–3 show that the proposed RK and TSRK frameworks perform effectively across the considered examples. Moreover, the numerical results for Gamma and Pareto claim-size distributions suggest that the TSRK methods yield consistently improved approximations over the one-step schemes based on Simpson’s rules.

Possible extensions of the present work include the study of multistep Runge–Kutta methods in this context and applications to the Gerber–Shiu function and Sparre Andersen-type risk models with diffusion perturbations, where the corresponding equations can be formulated as Volterra integral equations of the second kind and related integro-differential equations.

Data availability↩︎

No data was used for the research described in the article.

References↩︎