Monte Carlo simulation of the moments of a copula-dependent risk process with Weibull interwaiting time

Alhabshi, Sharifah Farah Syed Yusoff; Zamzuri, Zamira Hasanah; Ramli, Siti Norafidah Mohd Article Monte Carlo simulation of the moments of a copuladependent risk process with Weibull interwaiting time Risks Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Alhabshi, Sharifah Farah Syed Yusoff; Zamzuri, Zamira Hasanah; Ramli, Siti Norafidah Mohd (2021) : Monte Carlo simulation of the moments of a copula-dependent risk process with Weibull interwaiting time, Risks, ISSN 2227-9091, MDPI, Basel, Vol. 9, Iss. 6, pp. 1-13, https://doi.org/10.3390/risks9060109 This Version is available at: https://hdl.handle.net/10419/258197 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ risks Article Monte Carlo Simulation of the Moments of a Copula-Dependent Risk Process with Weibull Interwaiting Time Sharifah Farah Syed Yusoff Alhabshi , Zamira Hasanah Zamzuri and Siti Norafidah Mohd Ramli *   Citation: Syed Yusoff Alhabshi, Sharifah Farah, Zamira Hasanah Zamzuri, and Siti Norafidah Mohd Ramli. 2021. Monte Carlo Simulation of the Moments of a CopulaDependent Risk Process with Weibull Interwaiting Time. Risks 9: 109. https://doi.org/10.3390/risks9060109 Academic Editor: Mogens Steffensen Received: 15 March 2021 Accepted: 5 May 2021 Published: 3 June 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Department of Mathematical Sciences, Universiti Kebangsaan Malaysia, Bandar Baru Bangi 43600, Selangor, Malaysia; [email protected] (S.F.S.Y.A.); [email protected] (Z.H.Z.) *Correspondence: [email protected] Abstract: The widely used Poisson count process in insurance claims modeling is no longer valid if the claims occurrences exhibit dispersion. In this paper, we consider the aggregate discounted claims of an insurance risk portfolio under Weibull counting process to allow for dispersed datasets. A copula is used to define the dependence structure between the interwaiting time and its subsequent claims amount. We use a Monte Carlo simulation to compute the higher-order moments of the risk portfolio, the premiums and the value-at-risk based on the New Zealand catastrophe historical data. The simulation outcomes under the negative dependence parameter θ , shows the highest value of moments when claims experience exhibit overdispersion. Conversely, the underdispersed scenario yields the highest value of moments when θ is positive. These results lead to higher premiums being charged and more capital requirements to be set aside to cope with unfavorable events borne by the insurers. Keywords: aggregate discounted claims; overdispersed counting process; value-at-risk; solvency capital requirement; premium principle 1. Introduction The primary objective of an insurance company is to ensure that it can pay its promised obligations and remain solvent in the business. This can be achieved by managing surplus processes effectively and charging the appropriate premium amount, which will then guarantee an adequate reserve and capital requirement to cope with any unfavorable events. The ability to estimate accurately the aggregate discounted claims would help insurers manage their future liabilities and estimate appropriate premiums to be charged for an insurance contract under a particular line of business. The aggregate discounted claims, or the present value of the total claims amount paid out by an insurance company, are based on the number of claims occurrences up to time t , claims arrival time and claims amount, as well as the discount factor. Past studies related to the first two moments of the aggregate discounted claims were seen in Léveillé and Garrido (2001a) using renewal theory argument while, per Jang (2004), obtaining the Laplace transformation of the distribution using a different martingale approach. Léveillé and Garrido (2001b) applied the renewal theory arguments and conditioning on the first claim arrival to derive the m-th recursive aggregate moment, assuming independence of claim arrival time from its severity. The independence assumption between claims size and inter-claim arrival time in a classical risk model, as seen in Jang (2004), Waters (1983) and Yang and Zhang (2001) may no longer be appropriate for insurance risk portfolio modeling. With the increased frequency of catastrophic events (MunichRE 2018), the independence assumption will underestimate the risks faced by insurers, especially the values related to pricing, reserving and ruin measures of a risk portfolio. Past literature has relaxed the independence assumption used in the classical risk theory to include dependency elements between Risks 2021,9, 109. https://doi.org/10.3390/risks9060109 https://www.mdpi.com/journal/risks Risks 2021,9, 109 2 of 21 claims amounts and claim occurrences that follows a Poisson process (Barges et al. 2011; Mohd Ramli and Jang 2014;Sun et al. 2020;Woo and Cheung 2013). Barges et al. (2011) and Woo and Cheung (2013) were among the first to derive the m-th moment of the aggregate discounted claims using the Laplace transformation and relax the independence assumptions through a copula. Mohd Ramli and Jang (2014) extended the works by explicitly finding the Neumann series expressions and numerically solving the first two moments, where the dependence structure is captured by Farlie–Gumbel–Morgenstern (FGM), Gumbel, and Gaussian copulae. Mohd Ramli et al. (2018), on the other hand solved the recursive moments which take the form of a Volterra integral equation (VIE) using the Laplace transform, assuming an exponential claims size under the FGM, Frank and the HRT copulae. Recent study by Sun et al. (2020) derived the explicit expressions for the first three moments using a Laplace transformation by adopting the Spearman copula to cover a wide range of positive dependency, and developed the convex approximation of the copula. The studies mentioned above illustrated the use of the Poisson counting process to represent event occurrence. Owing to the practical limitation in the Poisson process which requires an equidispersed set of data, researchers has also adopted the Poisson-gamma or also known as negative binomial count process that only allows for an overdispersed dataset (Lora and Singer 2011). The Weibull count process, on the other hand, is a good alternative to allow for both underdispersed and overdispersed datasets, as commonly seen in applied data analysis with heterogeneous populations (Kreer et al. 2015;McShane et al. 2008; Winkelmann 1995). McShane et al. (2008) derived a closed form of the renewal process based on Weibull inter-arrival times using a Taylor series expansion. Additionally, as the constant hazard rate assumption is not valid in real practice, the duration dependence or time-varying function under a renewal process replaces the constant hazard function. The non-constant hazard rates that vary according to the duration of the interwaiting time (IWT) would be useful in modeling the unexpected random shocks that lead to the breakdown of many engineering systems, as seen in the work of Liu (2019). Compared with a simple Poisson count model with independent marginals, the Weibull count model with a copula is found to be more useful for predicting the number of goals scored in a football match (Boshnakov et al. 2017). The Weibull process has also been widely used in insurance and financial studies. Kreer et al. (2015) and Hasumi et al. (2009) mentioned that the probability distribution of the inter-occurrence time of the earthquake events and the automotive claim sizes distribution can be described by the Weibull distribution. It has also been applied in ruin theory see (Albrecher et al. 2011) to establish explicit formulas for ruin probabilities under a Sparre–Andersen risk process. Past researchers worked with explicit solutions of the moments derived from the Laplace transform (Barges et al. 2011;Li and Lu 2018;Mohd Ramli et al. 2018) or Neumann series (Mohd Ramli and Jang 2014) , which may require heavy computational capacity resulting from lengthy expressions. The variance expression derived from the Neumann series (Mohd Ramli and Jang 2014) has up to sextuple numerical integration even with a simple Poisson counting process. Thus, Mohd Ramli et al. (2019) work with numerical inversion of the Laplace transform via Gaver–Stehfest algorithm to compute the first two moments of Weibull risk process. However, the slow convergence of the numerical solutions under copulae other than FGM may be improved with better computational techniques. This study seeks to overcome the mathematical complexity under a Weibull risk process by using the Monte Carlo simulation to allow for faster computation and more flexibility to compute the higher-order moments through a wider choice of copula. We organize the remainder of this article as follows: In Section 2, we introduce the aggregate discounted claims model under Weibull counting process together with copulae to represent the dependency between the IWT and the subsequent claim sizes. In Section 3, we first show the comparison of the simulated moments under Monte Carlo simulation for a special case when λ= 1, with results obtained in previous studies to ensure consistency. Risks 2021,9, 109 3 of 21 We then fit the distributions and estimate the associated parameters to the New Zealand catastrophe occurrences and losses amount data to illustrate the dispersion in the dataset as allowed by a Weibull process. We also find the best copula to accommodate the dependency in the dataset. We then perform a scenario analysis to examine the characteristics of the resulting risk portfolio through the moments, the VaR and the premium computation. Section 4concludes the article. 2. Materials and Methods 2.1. Aggregate Risk Model We define the aggregate discounted claims as in Léveillé and Garrido (2001b), whereby Z={Z(t)}t≥0with: Z(t)=(∑N(t) i=1e−δTiXiif N(t)>0 0 if N(t)=0, (1) where Xi , i= 1, 2, . . . is a continuous, non-negative, independent and identically distributed (i.i.d.) random variable that represents claims sizes occurring at random times Ti , i= 1, 2, . . . , N(t) . The instantaneous rate of net interest δ is assumed to be deterministic. 2.1.1. Weibull Counting Process In this study, we let the counting process N={N(t)}t≥0 follow a basic Weibull counting process as in Mohd Ramli et al. (2019), whereby the corresponding continuous random variable of the inter-claim arrival time, Wi, is defined as: Wi=(Tifor i=1 Ti−Ti−1for i=2, 3, ..., (2) which follows a Weibull distribution. The time elapsed between two successive claims arrival is called the interwaiting time (henceforth IWT) of the counting process N(t) . Each pair of the joint variables forms a sequence of i.i.d. random vectors {(Xi,Wi)}i∈N . The mean and variance of a Weibull count model with shape parameter λ and scale parameter β, derived by McShane et al. (2008) are given by: E[N(t)] = ∞ ∑ n=1 ∞ ∑ j=n n(−1)j+ntλ/βjτn j Γ(λj+1), (3) Var[N(t)] =EhN2(t)i−E2[N(t)] = ∞ ∑ n=2 ∞ ∑ j=n n2(−1)j+ntλ/βjτn j Γ(λj+1)−  ∞ ∑ n=1 ∞ ∑ j=n n(−1)j+ntλ/βjτn j Γ(λj+1)  2 , (4) where τn j=     Γ(λj+1) Γ(j+1), for n=0, j=0, 1, 2, ... j−1 ∑ m=n−1 τn−1 j Γ(λj−λm+1) Γ(j−m+1), for n=1, 2, 3, ..., j=n,n+1, n+2, ..., in which the parameter ndenotes the number of events that occur in the interval [0, t]. Although the memoryless property of the exponentially distributed IWT allows for mathematical convenience, the Poisson counting process is only adequate if the data satisfy the restrictive assumption of equidispersion in which the variance of the data is equal to their mean. Hence, the Weibull count model is a better alternative for its ability to accommodate both overdispersed and underdispersed datasets. It also nests commonly used count models, including the Poisson and the negative binomial distributions. Additionally, Risks 2021,9, 109 4 of 21 the Weibull count model also allows for non-constant hazard rates that vary according to the duration of the IWT (McShane et al. 2008). The hazard function characterizes the IWT distribution and relates them to the type of dispersion observed in the corresponding count data. Past studies (Jose and Abraham 2011;McShane et al. 2008;Winkelmann 1995) have verified that the underlying IWT displays negative (positive) duration dependence under a decreasing (increasing) hazard for shape parameter of 0 <λ< 1 ( λ> 1) which causes overdispersion (underdispersion). The lack of duration dependence for λ = 1 leads to the Poisson distribution with a constant hazard function, as shown in Figure 1. We also note that there are two types of dependency applied when modeling the IWT under the Weibull counting process. First dependency lying between the duration of the IWT and the frequency of the claims arrival represented by non-constant hazard rates. As mentioned in McShane et al. (2008) and Winkelmann (1995), the negative (positive) duration dependence indicates a higher (lower) probability of claims occurring immediately after the occurrence of previous claims and decreases (increases) steadily as the IWT increases. This implies that the overdispersed (underdispersed) cases produce higher claims frequencies under short (long) IWT and vice versa. Figure 1. Weibull hazard function with different shape parameters to indicate different dispersions. 2.1.2. Copula The second dependency is between the IWT and its subsequent claims amount, which is captured by a copula. A copula provides a more convenient way to model any joint distribution of two or more random variables. We can model the marginal distributions of each variable by itself, and the copula can link these into a joint distribution. As we allow dependency between the claims size Xi , and the inter-arrival time Wi , we have a dependent Sparre–Andersen risk process with a dependence structure defined by a copula or other joint probability functions. Please note that when the marginals are negatively (positively) correlated, a large claim amount Xi will occur following a short (long) inter-arrival claim Risks 2021,9, 109 5 of 21 time Wi , and vice versa. In this study, the dependence structure between the two marginals of X and W are described by a few copulae including the Clayton, Frank and FGM copulae. Both the Clayton and Frank copulae belong to the same Archimedean family, in which the cumulative distribution function (CDF) of each copula is given by: CClayton θ(FX(x),FW(w)) =maxnFX(x)−θ+FW(w)−θ−1, 0o−1 θ, (5) where θ∈[−1, ∞), CFrank θ(FX(x),FW(w)) =−1 θln 1+(e−θFX(x)−1)(e−θFW(w)−1) (e−θ−1)!, (6) where θ∈R. The CDF of an FGM copula, which is from a different class of copula family is given by: CFGM θ(FX(x),FW(w)) =FX(x)FW(w)+θ(1−FX(x))(1−FW(s)), (7) where θ∈[−1, 1] represents the dependence parameter (Klugman and Parsa 1999;Ly et al. 2019). The copulae return to an independence copula when θ= 0. The FGM copula is a popular choice in extreme value analysis and risk management (Mao and Yang 2015) due to its simplicity and analytical tractability, despite its ability to capture only moderate dependency. The moderate dependence structure for FGM copula can be seen in Figure 2k–o , in which θ∈[−1, 1] is restricted to Kendall’s tau, τ∈[−2/9, 2/9] . In contrast, the Clayton copula (Figure 2a,e) and Frank copula (Figure 2f,j) allow for a wider range of dependency with τ∈[−1, 1] , as shown by the almost linear scatter plots of strong negative and positive dependency. However, the Frank copula could neither capture lower nor upper tail dependence, whereas the Clayton copula allows for lower tail dependence, or dependency on small values (Ly et al. 2019). The lower tail dependency can be seen clearly in its ’ice cream cone’ shapes as in Figure 2b,d. 2.2. Recursive Moment Expressions The general form of the m-th moment of aggregate discounted claims Z(t) , as given here, µ(m) z(T)=E[Zm(T)] =ZT 0fw(s)e−mδsE(Xm|W=s)ds +ZT 0fw(s)e−mδsµ(m) z(T−s)ds + m−1 ∑ j=1m jZT 0fw(s)e−mδsEXj|W=sµ(m−j) z(T−s)ds, (8) has a recursive nature (Barges et al. 2011;Mohd Ramli and Jang 2014). Working with explicit solutions of Equation (8) require heavy computational capacity due to lengthy expressions as seen in Mohd Ramli and Jang (2014) and Mohd Ramli et al. (2018). We therefore use the Monte Carlo simulation to compute the higher order of moments in a much shorter period with a common computational capacity. Hence, the mathematical properties of the insurance risk portfolio can be examined in a more efficient and viable manner despite the complicated expressions, which would be welcomed by the industrial community (Driels and Shin 2004). Risks 2021,9, 109 6 of 21 Figure 2. The simulated scatter plots of Clayton, Frank and FGM copulae under different dependence parameters. ( a ) Clayton with θ=− 1; ( b ) Clayton with θ=− 0.7; ( c ) Clayton with θ= 0; ( d ) Clayton with θ= 2; ( e ) Clayton with θ= 100; ( f ) Frank with θ=− 100; ( g ) Frank with θ=− 6; ( h ) Frank with θ= 0; ( i ) Frank with θ= 6; ( j ) Frank with θ= 100; (k) FGM with θ=−1; (l) FGM with θ=−0.5; (m) FGM with θ=0; (n) FGM with θ=0.5; (o) FGM with θ=1. 2.3. Monte Carlo Simulation Monte Carlo simulation is a popular and powerful quantitative tool often used in risk analysis such as in evaluation of insurers’ capital requirements (Casarano et al. 2017), estimation of claim size distribution (Bar-Lev and Ridder 2019), mortality projection (Zamzuri and Hui 2020) and bootstrapping approach in Structural Equation Modeling (Razak et al. 2019). We apply the Monte Carlo simulations to incorporate and understand the impact of risk and uncertainty in our model by computing the probability of different outcomes in a risk process (Zamzuri and Hui 2020). In this study, we simulate the aggregate discounted claims process using the estimated parameters and distributions from a real insurance dataset. In simulating the copula function, we adopt the same approach as Kelly (2007), whereby the first stage involves modeling the unidimensional marginal distributions while the second stage involves modeling the dependence structure. The steps taken in the simulation process are the following: Risks 2021,9, 109 7 of 21 1. Generate pairs of dependent random variates {(Xi,Wi)} from multivariate distributions constructed from the chosen copula. The multivariate distributions are based on the best-fit distribution obtained from the insurance dataset. 2. Compute the random time Ti, from the accumulated Wias in Equation (2). 3. Compute the aggregate discounted claims Z(T), by assuming deterministic δ. 4. Stop running the iterations once the Ti is above a pre-determined term of a policy contract. 5. Repeat the process from step 1 to 4 for nsimulations. 6. Determine the moments, premium and VaR from the simulated risk process Z(T). 3. Results and Discussion In this section, we organize the results and discussion by first showing the consistency of results returned by the Monte Carlo simulation method with results obtained in previous studies (Section 3.1). Next, we fit the aggregate discounted claims model to the New Zealand catastrophe data in Section 3.2. Using the estimated parameters of the claims size distribution, the IWT distribution and the dependency between these marginals, we then examine the risk portfolio of the New Zealand catastrophe dataset. This includes computing the value of moments, premium and VaR in Section 3.3. Finally, we perform scenario analysis in Section 3.4 by varying the parameters of the Weibull process and the copulae to examine the impact on the risk portfolio. 3.1. Results Verification We first consider the special case of a Weibull IWT with the shape parameter λ= 1 to ensure consistency of the Monte Carlo simulations with the numerical values obtained through Laplace transform as in Barges et al. (2011), Mohd Ramli and Jang (2014) and Mohd Ramli et al. (2018). We simulate pairs of the joint variables of random vectors (Xi,Wi)(i∈N) with an FGM copula defining the dependency between the variables. The values were computed using the same parameters in which X∼ Exp (α= 0.01 ) , W∼Exp(β=1) , δ= 0.04 and T= 5 for different values of θ . The average number of iterations for each simulation is 5, which is consistent with the average of the exponentially distributed IWT. The process is then repeated for n=100, 000 simulations. The mean and variance of the aggregate discounted claims obtained via Monte Carlo simulation are presented in Table 1and we compare the results with the explicit solution under the Laplace transform as in Mohd Ramli et al. (2018). The simulation results show only slight deviation from the values computed using the Laplace transform, of less than 1% differences for the mean and 2% for the variance. Table 1. Verification of the Monte Carlo Simulation. Moments θMonte Carlo Laplace Transform Relative Deviation −0.999 479.23 477.66 0.330% −0.9 474.83 475.23 0.084% −0.5 469.16 465.43 0.803% Mean 0 452.79 453.17 0.084% 0.5 440.59 440.92 0.074% 0.9 433.43 431.12 0.537% 0.999 430.23 428.69 0.359% −0.999 106,554.00 106,351.84 0.190% −0.9 105,554.20 103,929.50 1.563% −0.5 94,099.89 94,253.78 0.163% Variance 0 80,630.85 82,420.23 2.171% 0.5 69,601.11 70,874.44 1.797% 0.9 61,217.14 61,845.86 1.017% 0.999 59.182.13 59,638.74 0.766% Risks 2021,9, 109 8 of 21 3.2. Fitting Distribution and Parameter Estimation of Insurance Datasets In this section, we fit the aggregate discounted claim models on the New Zealand catastrophe occurrences and losses amount data from 1968 to 2014 which was retrieved from the CASdatasets package in R Statistical Software (Dutang and Charpentier 2020). We estimate the parameters of the fitted distributions of the claims size and the IWT, as well as the dependency between the marginals, using the maximum likelihood estimation method. 3.2.1. The Claim Sizes Distribution The Cullen and Frey (Cullen et al. 1999) graph in Figure 3illustrates the skewnesskurtosis plot of the empirical distribution of the New Zealand catastrophe loss claims. The values are computed on 1000 bootstrap samples to allow for the uncertainty of the estimated values of kurtosis and skewness from data. This graph exhibits some common distributions on the plot to help with the choice of distributions to fit to data. The normal, uniform, and logistic distributions are not good fits for the data since the skewness-kurtosis points of the distribution are far from the empirical value. The descriptive statistic in Table 2 shows a positive skewness and a heavy-tailed kurtosis. Thus, the fit of three common rightskewed distributions could be considered, which are the Weibull, Gamma and Lognormal distributions, together with other heavy-tailed distributions such as Pareto and Burr. The results in Table 3and Figure 4show that the lognormal distribution, with a meanlog of 1.35969 and an sdlog of 1.52856, provides the best fit for the claims size distribution (values in million USD). The fitted lognormal distribution has the lowest negative log likelihood, Akaike information criterion (AIC) and Bayesian information criterion (BIC) values, with significant p -values under the Kolmogorov–Smirnov (KS) and Anderson-Darling (AD) test. Table 2. Descriptive statistics for New Zealand catastrophe claims sizes. Min Max Median Mean Estimated Std Dev Estimated Skewness Estimated Kurtosis 0.01 112 3.6 10.67252 17.14509 2.997792 14.29496 Figure 3. Skewness-kurtosis plot for the New Zealand catastrophe claims sizes. Risks 2021,9, 109 15 of 21 Figure 7. Premium amount and VaR under a Clayton copula across different values for θ. Figure 8. Premium amount and VaR under a Frank copula across different values for θ. Risks 2021,9, 109 16 of 21 Table 12. Simulated VaR and premium ( κ= 0.1) under various dispersion of claims occurrences across different θunder a Frank copula when T= 5. Overdispersed WeibullC(λ=0.5, β=0.945) θVaR ($m) Premium Amount ($m) 95% 99.50% Mean Principle Std Dev Principle −100 225.816 566.751 73.363 76.605 −50 227.023 554.466 73.356 76.737 −10 222.041 544.253 72.451 75.522 0 179.770 475.843 59.577 62.559 10 74.381 223.736 31.722 32.858 50 34.979 45.038 21.940 21.018 100 33.045 38.072 21.427 20.435 Equidispersed WeibullC(λ=1, β=1.17835) θVaR ($m) Premium Amount ($m) 95% 99.50% Mean Principle Std Dev Principle −100 203.215 536.925 67.420 71.469 −50 204.435 532.854 67.100 70.079 −10 200.478 504.214 66.215 69.807 0 161.197 429.608 55.876 58.620 10 94.393 300.319 39.512 41.085 50 58.900 177.422 31.990 32.082 100 52.942 124.696 30.383 29.798 Underdispersed WeibullC(λ=2, β=1.5) θVaR ($m) Premium Amount ($m) 95% 99.50% Mean Principle Std Dev Principle −100 162.916 459.141 52.387 55.664 −50 161.985 445.691 52.355 55.655 −10 153.385 432.808 50.560 53.662 0 128.944 368.501 44.232 47.049 10 92.540 297.590 35.405 37.388 50 81.145 275.849 33.254 35.059 100 77.572 252.967 32.445 33.838 Figure 9. Premium amount and VaR under a FGM copula across different values for θ. Risks 2021,9, 109 17 of 21 Table 13. Simulated VaR and premium ( κ= 0.1) under various dispersion of claims occurrences across different θunder a FGM copula when T= 5. Overdispersed WeibullC(λ=0.5, β=0.945) θVaR ($m) Premium Amount ($m) 95% 99.50% Mean Principle Std Dev Principle −0.999 204.390 515.814 66.564 69.693 −0.5 190.822 509.914 62.891 66.048 0 178.699 480.092 59.106 62.356 0.5 165.488 448.837 55.451 58.282 0.999 154.182 426.935 51.789 54.361 Equidispersed WeibullC(λ=1, β=1.17835) θVaR ($m) Premium Amount ($m) 95% 99.50% Mean Principle Std Dev Principle −0.999 178.922 492.664 60.988 64.017 −0.5 172.109 458.163 58.608 61.381 0 161.823 436.838 56.028 58.640 0.5 151.357 423.809 53.247 56.144 0.999 140.298 411.168 50.556 53.115 Underdispersed WeibullC(λ=2, β=1.5) θVaR ($m) Premium Amount ($m) 95% 99.50% Mean Principle Std Dev Principle −0.999 140.496 403.169 46.940 49.749 −0.5 133.912 392.826 45.231 47.897 0 126.499 372.115 43.538 46.174 0.5 120.662 367.651 42.261 44.904 0.999 115.166 349.602 40.982 43.725 Figure 10. Cont. Risks 2021,9, 109 18 of 21 Figure 10. Premium amount under all copulae across different values for τ. Figure 11. VaR under all copulae across different values for τ. Risks 2021,9, 109 19 of 21 4. Conclusions The ability of an insurance firm to manage surplus effectively is important to ensure capital adequacy, as required by the Solvency II standard. To allow for dispersed datasets, we propose the Weibull counting process to represent the claims arrival process of an insurance risk portfolio and an arbitrary continuous claims size. Under this model, two types of dependency were accommodated which were captured by the hazard function and the copula. We then fit the model to the historical data of New Zealand catastrophe events occurrences and losses, which was best represented by an overdispersed Weibull process with claim amounts following a lognormal distribution and a Clayton copula with a weak dependency. The moments estimation via Monte Carlo simulations shows the highest mean and variance of the aggregate discounted claims were obtained under an overdispersed claim arrivals with a strong negative correlation. The positive skewness and the high value of kurtosis of the aggregate discounted claims with Weibull IWT also implied an asymmetric and a heavy-tailed risk distribution. The scenario analysis conducted indicated that insurers will need to charge higher premiums and set aside more capital requirements as the claims experience exhibits overdispersed (underdispersed) IWT with negative (positive) dependency captured by copulae, in comparison to claims experience with equidispersed IWT. Other than computing the respective ruin probability under the Weibull risk process and applying commonly used copulae in finance and insurance, we may also examine the applicability of the model on a random annuity plan with random payments being paid at random times. Author Contributions: Conceptualization, S.N.M.R.; methodology, S.N.M.R. and Z.H.Z.; software, Z.H.Z.; validation, S.N.M.R. and Z.H.Z.; formal analysis, S.F.S.Y.A.; investigation, S.F.S.Y.A.; resources, S.N.M.R.; data curation, Z.H.Z. and S.F.S.Y.A.; writing–original draft preparation, S.F.S.Y.A.; writing– review and editing, S.N.M.R. and Z.H.Z.; visualization, S.F.S.Y.A.; supervision, S.N.M.R. and Z.H.Z.; project administration, S.N.M.R. and Z.H.Z.; funding acquisition, S.N.M.R. and Z.H.Z. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Universiti Kebangsaan Malaysia Grant Number GUP-2020-033 . Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: http://cas.uqam.ca/,http://dutangc.perso.math.cnrs.fr/RRepository/ (accessed on 15 March 2021). Conflicts of Interest: The authors declare no conflict of interest. The funder had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Abbreviations The following abbreviations are used in this manuscript: AD Anderson-Darling AIC Akaike information criterion BIC Bayesian information criterion CDF Cumulative distribution function FGM Farlie-Gumbel-Mogenstern GOF Goodness-of-fit IWT Interwaiting time KS Kolmogorov–Smirnov MDPI Multidisciplinary Digital Publishing Institute SCR Solvency capital requirement Std Dev Standard deviation VaR Value-at-Risk Risks 2021,9, 109 20 of 21 References Albrecher, Hansjörg, Corina Constantinescu, and Stéphane Loisel. 2011. Explicit ruin formulas for models with dependence among risks. Insurance: Mathematics and Economics 48: 265–70. [CrossRef] Bar-Lev, Shaul, and Ad Ridder. 2019. Monte Carlo methods for insurance risk computation. International Journal of Statistics and Probability 8: 54–74. [CrossRef] Barges, Mathieu, Hélene Cossette, Stéphane Loisel, and Etienne Marceau. 2011. On the moments of aggregate discounted claims with dependence introduced by a FGM copula. ASTIN Bulletin: The Journal of the IAA 41: 215–38. Boshnakov, Georgi, Tarak Kharrat, and Ian G. McHale. 2017. A bivariate Weibull count model for forecasting association football scores. International Journal of Forecasting 33: 458–66. [CrossRef] Casarano, Giuseppe, Gilberto Castellani, Luca Passalacqua, Francesca Perla, and Paolo Zanetti. 2017. Relevant applications of Monte Carlo simulation in Solvency II. Soft Computing 21: 1181–92. [CrossRef] Christiansen, Marcus C., and Andreas Niemeyer. 2014. Fundamental definition of the solvency capital requirement in Solvency II. Astin Bulletin 44: 501–33. [CrossRef] Cullen, Alison C., H. Christopher Frey, and Christopher H. Frey. 1999. Probabilistic Techniques in Exposure Assessment: A Handbook for Dealing with Variability and Uncertainty in Models and Inputs. Berlin: Springer Science & Business Media. Driels, Morris R., and Young S. Shin. 2004. Determining the Number of Iterations for Monte Carlo Simulations of Weapon Effectiveness. Technical Report. Monterey: Naval Postgraduate School Monterey. Dutang, Christophe, and Arthur Charpentier. 2020. CASdatasets: Insurance Datasets. R Package Version 1.0-11. Vienna: R Core Team. Hasumi, Tomohiro, Takuma Akimoto, and Yoji Aizawa. 2009. The Weibull–log Weibull distribution for inter-occurrence times of earthquakes. Physica A: Statistical Mechanics and Its Applications 388: 491–98. [CrossRef] Jang, Jiwook. 2004. Martingale approach for moments of discounted aggregate claims. Journal of Risk and Insurance 71: 201–11. [CrossRef] Jose, Kanichukattu K., and Bindu Abraham. 2011. A count model based on Mittag-Leffler interarrival times. Statistica 71: 501–14. Kelly, Dana L. 2007. Using copulas to model dependence in simulation risk assessment. Paper presented at the ASME 2007 International Mechanical Engineering Congress and Exposition, Seattle, WA, USA, November 11–15; vol. 43084, pp. 81–89. Klugman, Stuart A., and Rahul Parsa. 1999. Fitting bivariate loss distributions with copulas. Insurance: Mathematics and Economics 24: 139–48. [CrossRef] Kreer, Markus, Ay¸se Kızılersü, Anthony W. Thomas, and Alfredo D Egídio dos Reis. 2015. Goodness-of-fit tests and applications for left-truncated Weibull distributions to non-life insurance. European Actuarial Journal 5: 139–63. [CrossRef] Léveillé, Ghislain, and Jose Garrido. 2001a. Moments of compound renewal sums with discounted claims. Insurance: Mathematics and Economics 28: 217–31. Léveillé, Ghislain, and José Garrido. 2001b. Recursive moments of compound renewal sums with discounted claims. Scandinavian Actuarial Journal 2001: 98–110. [CrossRef] Li, Shuanming, and Yi Lu. 2018. On the moments and the distribution of aggregate discounted claims in a Markovian environment. Risks 6: 59. [CrossRef] Liu, Hanlin. 2019. Reliability and maintenance modeling for competing risk processes with Weibull inter-arrival shocks. Applied Mathematical Modelling 71: 194–207. [CrossRef] Lora, Mayra Ivanoff, and Julio M. Singer. 2011. Beta-binomial/gamma-Poisson regression models for repeated counts with random parameters. Brazilian Journal of Probability and Statistics 25: 218–35. [CrossRef] Ly, Sel, Kim-Hung Pho, Sal Ly, and Wing-Keung Wong. 2019. Determining distribution for the product of random variables by using copulas. Risks 7: 23. [CrossRef] Mao, Tiantian, and Fan Yang. 2015. Risk concentration based on expectiles for extreme risks under FGM copula. Insurance: Mathematics and Economics 64: 429–39. [CrossRef] McShane, Blake, Moshe Adrian, Eric T. Bradlow, and Peter S. Fader. 2008. Count models based on Weibull interarrival times. Journal of Business & Economic Statistics 26: 369–78. Mohd Ramli, Siti Norafidah, and Jiwook Jang. 2014. Neumann series on the recursive moments of copula-dependent aggregate discounted claims. Risks 2: 195–210. [CrossRef] Mohd Ramli, Siti Norafidah, Nur Atikah Mohamed Rozali, Sharifah Farah Syed Yusoff Alhabshi, and Ishak Hashim. 2018. Laplace transform on the recursive moments of copula-dependent aggregate discounted claims. In AIP Conference Proceedings. Melville: AIP Publishing LLC, vol. 1974, p. 020110. Mohd Ramli, Siti Norafidah, Nur Atikah Mohamed Rozali, Sharifah Farah Syed Yusoff Alhabshi, and Ishak Hashim. 2019. Laplace transform on the recursive moments of aggregate discounted claims with Weibull interwaiting time. In AIP Conference Proceedings. Melville: AIP Publishing LLC, vol. 2184, p. 050017. MunichRE. 2018. NatCatSERVICE. Available online: http://natcatservice.munichre.com (accessed on 14 February 2018). Razak, Nor Iza Anuar, Zamira Hasanah Zamzuri, and Nur Riza Mohd Suradi. 2019. The Implementation Of Double Bootstrap Method In Structural Equation Modeling. Asm Science Journal 12: 8–14 Sun, Weiwei, Xiang Hu, and Lianzeng Zhang. 2020. Moments of discounted aggregate claims with dependence based on Spearman copula. Journal of Computational and Applied Mathematics 377: 112889. [CrossRef] Risks 2021,9, 109 21 of 21 Waters, Howard R. 1983. Probability of ruin for a risk process with claims cost inflation. Scandinavian Actuarial Journal 1983: 148–64. [CrossRef] Winkelmann, Rainer. 1995. Duration dependence and dispersion in count-data models. Journal of Business & Economic Statistics 13: 467–74. Woo, Jae-Kyung, and Eric C. K. Cheung. 2013. A note on discounted compound renewal sums under dependency. Insurance: Mathematics and Economics 52: 170–79. [CrossRef] Yang, Hailiang, and Lihong Zhang. 2001. On the distribution of surplus immediately after ruin under interest force. Insurance: Mathematics and Economics 29: 247–55. [CrossRef] Zamzuri, Zamira Hasanah, and Gwee Jia Hui. 2020. Comparing and forecasting using stochastic mortality models: A Monte Carlo simulation. Sains Malaysiana 49: 2013–22. [CrossRef]