Ruin Analysis on a New Risk Model with Stochastic Premiums and Dependence Based on Time Series for Count Random Variables

In this paper, we propose a new discrete-time risk model of an insurance portfolio with stochastic premiums, in which the temporal dependence among the premium numbers of consecutive periods is fitted by the first-order integer-valued autoregressive (INAR(1)) process and the temporal dependence among the claim numbers of consecutive periods is described by the integer-valued moving average (INMA(1)) process. To measure the risk of the model quantitatively, we study the explicit expression for a function whose solution is defined as the Lundberg adjustment coefficient and give the Lundberg approximation formula for the infinite-time ruin probability. In the case of heavy-tailed claim sizes, we establish the asymptotic formula for the finite-time ruin probability via the large deviations of the aggregate claims. Two numerical examples are provided in order to illustrate our theoretical findings.


Introduction
As an absolutely necessary part of the modern financial system, insurance is one of the most effective ways for people to manage risks, such that it plays a significant role in our daily life. A very important task of insurance companies is to quantitatively analyze future claims. Consequently, risk theory has become an active research field of actuarial science. For the classical mathematical risk model, the so-called Lundberg-Cramér surplus process has the following form: in which u ≥ 0 is the initial capital of an insurance portfolio, c > 0 is the constant rate of premium income, {N 0 t , t ≥ 0} is a homogeneous Poisson process with intensity λ, the total claim numbers are denoted up to time t, and Y i describes the size of the ith claim. In the literature, Asmussen and Albrecher [1] presented excellent reviews about this well-known and important model.
In model (1), independent structures are usually assumed. For example, the claim amount {Y i , i ≥ 1} is a sequence of non-negative independent and identically distributed (i.i.d.) random variables, and the claim numbers of different periods are assumed to be a sequence of i.i.d. random variables. However, these are not always true in practice because of the increasing complexity of individual risks. To avoid this restriction, a growing number of actuaries have been paying attention to the model with dependent risks. As stated in Yang and Zhang [2], there are mainly two kinds of correlation in insurance: one is the correlation among lines of businesses, and the other is temporal dependence, such as the correlation between the current claim and the previous claims. For recent works about the first type of correlation, Refs. [3,4] studied the dependence among individual risks, Refs. [5,6] discussed the two-dimensional risk models with dependent surplus processes, and [7,8] examined the risk models that have multiple classes of insurance business with thinning dependence structure. The relevant results have been used in a variety of actuarial areas, including, among others, value at risk, dividend strategies, reinsurance, capital allocation, etc.
In this paper, we focus on the second type. To deal with this problem, the use of a time series is a critical method. Gerber [9] considered the calculation of ruin probabilities in a Gaussian linear risk model; Gourieroux and Jasiak [10] applied the integer-valued time series model to update the premiums in vehicle insurance; and many researchers have extensively revisited the relevant results afterwards. Considering that the compound distributions are the cornerstones of a great number of risk models in risk theory, Cossette et al. [11] proposed some new discrete-time risk models, where the first-order integer-valued moving average (INMA(1)) and first-order integer-valued autoregressive (INAR(1)) processes are used to describe the dependence structures among the number of claims for each period. The authors derived expressions for the functions that allow people to find the Lundberg adjustment coefficients and discussed the Lundberg approximation formulas for infinite-time ruin probabilities. Along the same line, Cossette et al. [12] determined the distributions of aggregate claim amount and provided an effective way to measure some related risk quantities, including VaR, TVaR, and the stop-loss premium. Shi and Wang [13] gave an approximation method for the risk model with the Poisson INAR(1) claim number process in order to obtain the upper bound of the infinite-time ruin probability. Zhang et al. [14] solved the problem of optimal reinsurance strategy for the risk model with the INMA(1) claim number process. Afterwards, Hu et al. [15] and Chen and Hu [16] further generalized this kind of model by replacing the Poisson innovations with compound Poisson innovations in the INAR(1) and INMA(1) claim number processes, respectively. Guan and Hu [17] even utilized an INAR(1) process with an arbitrary innovations' distribution to specify the temporal dependence among the claim numbers.
In the papers mentioned above, it should be noted that the incomes of all the risk models are linear functions of time t, because the premiums are collected continuously with positive deterministic constant rate c, providing great convenience for risk analysis. However, this assumption is obviously lacking in terms of describing the real situation of insurance portfolios; for example, it cannot capture the uncertainty of the customers' arrivals. As an alternative to a fixed premium rate, Boikov [18] supposesed that the premium income also follows a compound Poisson process and calculates the ruin probability. From then on, the risk models with stochastic premiums have been extensively improved by many actuaries. Wang et al. [19] studied the investment problem of such models. Labbé and Sendova [20] discussed the Gerber-Shiu function. Zhao and Yin [21] proposed a renewal risk model with stochastic incomes. Recently, Su et al. [22] provided a statistical method for estimating the Gerber-Shiu function; Ragulina [23] investigated the De Vylder approximation for the ruin probability and a constant dividend strategy in the risk model with stochastic premiums; and Dibu and Jacob [24] focused on a double barrier hybrid dividend strategy. Wang et al. [25] quantitatively assessed the impact of the stochastic income process on some ruin quantities in detail.
Similar to the classical risk model, the premium numbers of different periods are commonly set to be a sequence of i.i.d. random variables in the aforementioned papers. To better characterize the uncertainty and capture the variability of an insurer's income process, Guan and Wang [26] proposed modeling the temporal dependence among the premium numbers of each period by a Poisson INAR(1) process. In this paper, we follow this trend of research. We also aim to study a new dependent risk model with stochastic premiums based on time series for count random variables, in which the INAR(1) process and INMA (1) process are applied to fit the temporal dependence among the premium numbers and the temporal dependence among the claim numbers of consecutive periods, respectively. Our goal is to approximate the infinite-time ruin probability of the proposed surplus process by the Lundberg adjustment coefficient and discuss the asymptotic formula for the finite-time ruin probability when the claim sizes follow distributions with heavy tails.
Our model generalizes the classical discrete-time surplus process of an insurance portfolio with stochastic premiums to a new dependent risk model, and our results extend what has been studied in the existing literature. The contributions of our paper mainly include the following two aspects: • In contrast to the assumption that either claim numbers or premium numbers have a temporally dependent structure, we propose a new risk model of an insurance portfolio with both claim numbers and premium numbers being dependent within the integer-valued time series framework, which is more flexible in insurance practice. • In addition to studying the distribution of the aggregate claims, the Lundberg adjustment coefficient, and the Lundberg approximation formula for the infinite-time ruin probability in the case of light-tailed claim sizes, we also explore the large deviations of the aggregate claims and the asymptotic formula for the finite-time ruin probability when the claim sizes are heavy-tailed, which enlarges the applicability of the risk model.
The remainder of the paper is organized as follows: Section 2 introduces our concerned risk model and considers some probabilistic properties of the proposed model. Section 3 defines the Lundberg adjustment coefficient via the solution of an explicit equation. Section 4 establishes an exponential asymptotic estimation for the infinite-time ruin probability. Section 5 studies the large deviations of the aggregate claims when the claim sizes follow a class of heavy-tailed distributions and presents an asymptotic formula for the finite-time ruin probability. Section 6 illustrates the main results by numerical simulations. Section 7 finally concludes this paper.

Risk Model and Basic Properties
In this section, we first describe the new dependent risk model, and then, provide some moment results of the premiums and claims. Let U t be the surplus of an insurance portfolio at the end of period t, and we define the surplus process by the dynamic equation where U 0 = u ≥ 0 is the initial surplus level; P t = M t ∑ k=1 X t,k aggregates the premiums during period t, in which M t counts the number of individual income and X t,k represents the amount of the kth premium income for the insurance portfolio during period t; is the aggregate claims during period t, in which N t denotes the number of claims and Y t,j is the size of the jth payment to the insured in period t. For mathematical tractability, the following assumptions are made: (1) Both {X t,k , t = 1, 2, · · · , k = 1, 2, · · · } and {Y t,j , j = 1, 2, · · · , k = 1, 2, · · · } are arrays of i.i.d. random variables, which have the same distributions as non-negative X and Y, respectively.
The dependence structures of the model are constructed in the following ways: (i) {M t , t = 1, 2, · · · } constitutes a Poisson INAR(1) process that satisfies where the so-called binomial thinning operator "•" is given by in which the following statements are true: • The thinning parameter α ∈ [0, 1).

Remark 1.
Time series analysis is one of the most important methods for dealing with dependent data and has attracted a lot of interest during the last decades. However, the classical real-valued time series models with continuous ranges can not account for discreteness, so they are of limited use for fitting the premium numbers and the claim numbers, which are typical count random variables fairly common in practice. Their poor performances in modeling this class of data mainly include: (1) the data generating mechanism can not be explained; (2) the approximate errors are big; and (3) the forecast results are not integer-valued. Therefore, models and methods for integer-valued time series have been covered by a large number of papers in recent years. Refs. [27][28][29][30] present comprehensive surveys on this fascinating research area. As two core models of integer-valued time series, INAR(1) process and INMA(1) process have been extensively applied in the literature of actuarial science, and the relevant results have been widely used in a variety of risk management.

Remark 2.
The INAR(1) process (3) shows that the premium number in period t is composed of two parts: ε t denotes the new incomes arriving between period t − 1 and t, and α • M t−1 presents a random proportion of the premium number in the previous period. This can be reasonably explained for the insurance practice that states: every insured entity could continue to pay a premium with probability α; or withdraw from the contract with probability 1 − α in the next period. When α = 0, (3) becomes M t = ε t , meaning that the premium number in period t could be totally determined by ε t , and our model (2) will reduce to the classical discrete-time risk model with stochastic premiums, where the premium numbers of different periods are independent (please see Appendix A for details).

Remark 3.
The INMA(1) process (5) reveals that the claim number in period t also consists of two parts: η t is the new claim during period t, and β • η t−1 indicates the claims of period t − 1 that could produce another accident with probability β in period t. Instead of (3), we use the INMA(1) process (5) to fit the temporal dependence among the claim numbers for each period, considering that the insured parties cannot receive benefits every year for some insurance products. Taking unemployment insurance as an example, every time the claimant is out of work, they could receive the benefits for up to 2 years, if the premiums for at least 1 year have been paid. Another appropriate example might be some medical insurance contracts, which state that no matter how long the patient stays in the hospital, the insurer would pay the benefits for at most (for instance) 2 months. Similarly, if β = 0, our proposed model will reduce to the classical case, where the claim numbers of different periods are independent.
As stated in Al-Osh and Alzaid [31], under the condition of 0 ≤ α < 1, it follows that the process of premium numbers {M t , t = 1, 2, · · · } is a stationary and ergodic Markov chain. Furthermore, if we assume ε t ∼ P(λ 1 ), then M t is also Poisson distributed with mean λ 1 1−α . Hence, by the law of iterated expectation and the assumption that {X t,k , k = 1, 2, · · · } and M t are independent, it is easy to find that Meanwhile, by the law of total variance, we can obtain Furthermore, Al-Osh and Alzaid [31] show that Similarly, for the process of claim numbers {N t , t = 1, 2, · · · }, under the condition of 0 ≤ β < 1, its marginal distribution is uniquely determined by the law of {η t , t = 0, 1, · · · }. Therefore, the assumption of η t ∼ P(λ 2 ) will result in N t being Poisson distributed with a mean of (1 + β)λ 2 . Consequently, by the same method to drive (7)-(9), we have These results are consistent with those in [11].

Definition of the Lundberg Adjustment Coefficient
In this section, we first consider how to calculate the moment generating functions (m.g.f.) of the aggregate premiums and aggregate claims up to period t, and then, define the Lundberg adjustment coefficient of the proposed dependent risk model with stochastic premiums based on time series for count random variables by means of a equation.
After recursive calculation, we can rewrite model (2) as Y i,j represent the aggregate premium incomes and aggregate claim payments up to time t, respectively. As for the m.g.f. of W t and S t , by the definition, we have that where M X (·) denotes the m.g.f. of X and P M 1 +···+M t (·) presents the probability generating function (p.g.f.) of the total premium number up to period t of the proposed model (2). Similarly, it holds that where M Y (·) denotes the m.g.f. of Y, and P N 1 +···+N t (·) presents the p.m.f. of the total claim number up to period t of the proposed model (2). In order to compute M W t (r) and M S t (r), we find the explicit expressions for P M 1 +···+M t (·) and P N 1 +···+N t (·) in the following two lemmas, respectively. Lemma 1. For t = 1, 2, · · · , when 0 ≤ s ≤ 1, the p.g.f. of M 1 + · · · + M t is given by Proof. Since M 1 ∼ P( λ 1 1−α ), it is obvious that When t ≥ 2, we denote By the property of the binomial thinning operator (see Scotto et al. [28] for example), we can rewrite M 1 + · · · + M t as For the p.g.f. calculation, we have in which h 1 (s) = s and h t (s) = s(αh t−1 (s) Similarly, we can obtain (17) Combining (15)- (17), it follows that Moreover, from the definition h t (s) = s(αh t−1 (s) Then, recursive calculation results in Finally, inserting (19) into (18), we can obtain This completes the proof.
To further analyze the insurance portfolio, we write and let Then, the positive solution to the equation c(r) = 0 can be defined as the Lundberg adjustment coefficient, which is denoted by R and can be used to approximate the infinitetime ruin probability of the proposed model (2). The following result gives the explicit expression for c(r).
Proof. Due to the non-negativity of r and X, it follows that Then, by Lemma 1, we have On the other hand, from (13) and (20), we obtain Then, combining (24) and (25) with (21) and (22) yields This completes the proof.

Remark 4.
When α = 0, the proposed model (2) degenerates to the discrete-time risk model based on the Poisson INMA(1) process studied by [11,14], where only the temporal dependence among the claim numbers of consecutive periods is considered. Consequently, (23) becomes which corresponds to (7) in [11,14].

Remark 5.
When β = 0, the proposed model (2) reduces to the discrete-time risk model with stochastic premiums and dependence based on the Poisson INAR(1) process studied by [26], where only the temporal dependence among the premium numbers of consecutive periods is considered. As a result, (23) becomes which corresponds to (3.10) in [26].

Lundberg Approximation Formula for the Infinite-Time Ruin Probability
Let the ruin time of our proposed surplus process (2) be T = inf U t goes below 0 at least once; otherwise, take T = +∞. As a consequence, the infinite-time ruin probability ψ(u) is defined by Ruin probability ψ(u) is well-known as one of the most common and important quantities used to measure the riskiness of an insurance portfolio in the risk-theoretic context. However, as can be seen from the expression (11), our proposed model releases the condition that P t and L t are independent of U t−1 , which is a key but defective assumption in the classical risk model with stochastic premiums and allows for the temporal dependence among the premium numbers and claim numbers. Therefore, P t and L t are correlated with U t−1 , and {U t , t = 1, 2, · · · } is no longer a Lévy process with stationary independent increments in our model. Consequently, it is not easy to derive the upper bounds and explicit expression for the infinite-time ruin probability such as those in some classical models. As an efficient alternative, the following result gives an asymptotic estimation for ψ(u) of our proposed model (2).
we can obtain the Lundberg approximation formula for the infinite-time ruin probability ψ(u), which has the following expression: where u and R are the initial capital and the Lundberg adjustment coefficient, respectively.
Proof. According to Theorem 2.1 in Müller and Pflug [32], it is sufficient for us to prove that the equation c(r) = 0 has a unique positive solution, which can be defined as the Lundberg adjustment coefficient R. To this end, we derive the following four properties of the function c(r). Firstly, noting that M X (0) = M Y (0) = 1, we have Secondly, it is easy to calculate that Together with the fact M X (0) = E(X) and M Y (0) = E(Y), we obtain Thirdly, it is easy to verify the convexity of c(r), which results from the fact that c t (r) is convex and the definition of c(r) = lim t→+∞ c t (r).
Finally, when the m.g.f. of Y exists, i.e., there exists some quantity r 0 , 0 < r 0 ≤ +∞, such that M Y (r) is finite for all r < r 0 with Therefore, it can be concluded that there exists a unique positive solution to the equation c(r) = 0, and then, (27) follows immediately.

Remark 6.
In risk and ruin theory, the assumption (26) is the so-called relative safety loading condition, which implies that the expected premium incomes should be more than the expected claim expenses to guarantee that the insurance company can operate normally and profitably.

Remark 7.
As a result of the approximation formula (27), we can asymptotically estimate the infinite-time ruin probability ψ(u) by if the initial surplus u becomes large enough.
From (9) and (10), it can be seen that the thinning parameters α and β could quantitatively measure the degree of the dependence in the risk model (2); hence, it is necessary for us to discuss their impacts on the adjustment coefficient and further on the risk of the insurance portfolio.

Proposition 1.
As a function of the two thinning parameters, the Lundberg adjustment coefficient R of our proposed risk model (2) increases with respect to α and decreases with respect to β.
Proof. For convenience, we now rewrite c(r) as c(α, β, r); the Lundberg adjustment coefficient R is determined by c(α, β, R) = 0 and can be taken as a function of α and β. By the properties derived in the proof of Theorem 2, we know that Meanwhile, with R > 0 in mind, it follows that 0 ≤ M X (−R) < 1. Thus, we take the partial derivative of c(α, β, R) with respect to variable α and then have ∂c(α, β, R) ∂α As a result, using implicit function theorem, it holds that implying that R increases with respect to α. Similarly, because M Y (R) > 1 for R > 0, taking the partial derivative of c(α, β, R) with respect to variable β yields ∂c(α, β, R) ∂β from which we apply implicit function theorem again and obtain ∂R ∂β meaning that R decreases with respect to β.

Remark 8.
As shown in Proposition 1, the degree of riskiness can be measured and quantified by the Lundberg adjustment coefficient R, in the sense that it decreases with the thinning parameter α, while it increases with the thinning parameter β. In insurance practice, it can be naturally explained that when α increases, the insured parties would like to renew their insurance contracts with a higher probability in the next period, which would lower the risk of the portfolio. On the contrary, when β increases, a reported claim becomes more likely to produce another insurance accident in the next period, which could make the portfolio much riskier.

Asymptotic Formula for the Finite-Time Ruin Probability
In this section, we turn our focus to the case of heavy-tailed claim sizes, which are frequently used in insurance practice for catastrophe risks, such as earthquakes, hurricanes, floods, financial crises, agricultural disasters, and so on. In these instances, the Lundberg adjustment coefficient and Lundberg approximation estimation for infinite-time ruin probability can no longer be applied because M Y (r) (the m.g.f. of Y) does not exist for r > 0. Therefore, increasing numbers of researchers have increasingly paid close attention to the precise large deviations in the aggregate of claims, as well as the asymptotic formulas for infinite-time and finite-time ruin probabilities. The relevant study was initiated by Klüppelberg and Mikosch [33] and then has been revisited by many researchers afterwards. We refer to Chen et al. [34] and Fu et al. [35] for some recent contributions on this topic. Cheng and Wang [36], Yang et al. [37], and Jing et al. [38] considered the asymptotic ruin probabilities in risk models with dependence among the claim sizes. Xun et al. [39] obtained the uniformly asymptotic result of ruin probability in a general risk model with stochastic premiums. Yu [40] derived the precise large deviations of the aggregate amount of claims for a risk model with the Poisson ARCH claim number process. Along the same line, in this section, we investigate our proposed model (2) when the distribution of claim sizes belongs to a heavy-tailed class.
First, we give some brief notations. Let a(x) and b(x) be two positive functions. We de- We denote the common distribution functions of premium amount X and claim size Y with F X (x) and F Y (y), respectively.
Then, we recall a class of heavy-tailed distributions and one of its important properties. More detailed discussions can be found in Embrechts et al. [41], Asmussen and Albrecher [1], etc.
A distribution function F on [0, ∞] is said to have a consistently varying tail, denoted where F(x) is the tail probability with F(x) = 1 − F(x). The class C is a wide class of distributions commonly used in actuarial science, including the well-known Pareto, Burr, and loggamma distributions. Ng et al. [42] established a very useful result for the distributions of class C, which is given in the following lemma.

Lemma 3.
Suppose that {Y j , j = 1, 2, · · · } is a sequence of i.i.d. non-negative random variables with common distribution function F Y (y) ∈ C and E(Y) < +∞. Taking Q t = t ∑ j=1 Y j , for any fixed γ > 0, it holds uniformly for all y > γt that in which the uniformity is understood in the following sense: Analogous to the infinite-time ruin probability ψ(u), for any fixed t = 1, 2, · · · , we define the finite-time ruin probability ψ(u, t) of the discrete-time risk model (2) as ψ(u, t) = P(T ≤ t|U 0 = u).
In order to further study the asymptotic formula of ψ(u, t), which is also a core actuarial quantity, we revise Lemma 3 as follows.

Lemma 4.
Suppose that {Y j , j = 1, 2, · · · } is a sequence of i.i.d. non-negative random variables with the common distribution function F Y (y) ∈ C and E(Y) < +∞. Define Q t = t ∑ j=1 Y j ; then, for any fixed γ > 0 and δ > 0, it holds uniformly for all y > γt 1+δ that P(Q t > y) ∼ tF Y (y), t → +∞. (34) Proof. By the definition of class C, it follows for any fixed θ > 0 and sufficiently large y that from which we can obtain Hence, it holds that Furthermore, by Lemma 3 and (35), it follows that The proof is then completed. Now, we give the precise large deviations of the aggregate claims, S t , which is described in model (11). (2), let F Y (y) and E(Y) be the common distribution function and expectation of the claim sizes, respectively. Assuming F Y (y) ∈ C and E(Y) < +∞, then for any fixed γ > 0 and δ > 0, it holds uniformly for all y > γt 1+δ that

Theorem 3. For our proposed model
Proof. Let {Y j , j = 1, 2, · · · } be a sequence of i.i.d. non-negative random variables, with their common distribution function denoted by F Y (y). Suppose that ϕ S t (r) is the With the same method to derive (12) and (13), we can obtain where ϕ Y (r) is the characteristic function of Y.
Combining (40)- (42) gives On the other hand, it holds that in which because of the fact that Generally, letting η ↓ 0 in (43) and (46) and keeping (39) in mind, we finally conclude that Then, the proof is completed.
With the help of the above conclusion, we can manifest the asymptotic formula for the finite-time ruin probability in the following theorem.
Proof. From the definition of finite-time ruin probability, it is clear that X i,k and keeping (9) in mind, for any η > 0, we have Then, for any θ > 0, if t is sufficiently large such that u is sufficiently large, it holds that Furthermore, by Theorem 3 and let θ ↓ 0, we have Plugging (49) into (48) gives On the other hand, for any fixed γ > 0 and δ > 0, we have uniformly for all u > γt 1+δ that Therefore, we complete the proof by combining (50) and (51).

Remark 9.
Applying Lemma 3 instead of Lemma 4 in Theorem 3, it is not difficult to see that the precise large deviation (36) also holds uniformly for all y > γt. In this paper, we restrict ourselves to the interval y > γt 1+δ in order to provide convenience for investigating the finite-time ruin probability ψ(u, t). Moreover, we can prove that the asymptotic formula (47) in Theorem 4 holds uniformly for all u ∈ Ω = {u; t = o(u)}, which includes u > γt 1+δ as a special case. In practice, when t is large enough, we can asymptotically estimate ψ(u, t) by λ 2 (1 + β)tF Y (u), as the size of claims belong to the distributions of class C and the insurer's initial surplus is adequate in the sense of u > γt 1+δ .

Numerical Examples
In this section, we aim to perform some numerical simulations to demonstrate and assess the Lundberg adjustment coefficient and the Lundberg approximation results for the infinite-time ruin probability ψ(u), as well as the asymptotic formula for the finite-time ruin probability ψ(u, t), of our proposed model.

Example 1.
We suppose that the gain amount X and the claim size Y follow exponential distributions that have means 1/µ 1 and 1/µ 2 , respectively. Therefore, we have the moment generating functions of X and Y as follows: Then, from Theorem 2, c(r) = 0 is equivalent to The unique positive solution to Equation (53) can be found but appears tedious. In what follows, we give some numerical results to show the properties and performance of R and e −Ru . Without loss of generality, we set λ 1 = 1, λ 2 = 0.4, µ 1 = 1, and µ 2 = 0.5, and then, calculate and discuss the Lundberg adjustment coefficient R and the approximated ruin probability e −Ru for different values of α and β. When we consider the impacts of α and β on the main results, it should be noted that the relative safety loading condition (26) has to be satisfied, i.e., which, in our parameter scenario, implies Table 1 gives the computed values of Lundberg adjustment coefficients corresponding to different values of α and β. We also illustrate these results in Figure 1, from which it can be clearly seen that R increases as α increases, implying that the insurance portfolio would become less and less dangerous because the approximated infinite-time ruin probability e −Ru decreases. In the same sense, when β increases, R will decrease, meaning that there could be higher risks in the insurance portfolio. (In Table 1, the notation " − " means that the values of α and β do not satisfy the relative safety loading condition, and the Lundberg adjustment coefficients are not considered for these situations.)  In order to evaluate the performance of the approximated infinite-time ruin probability e −Ru , we fix α = β = 0.5 in the proposed risk model and compute the true ruin probabilities corresponding to different values of u by the Monte Carlo method used in Albrecher and Kantor [44]. For this purpose, we randomly draw sample paths according to the Poisson INAR(1) process and the Poisson INMA(1) process for the premium arrivals {M t , t = 1, 2, · · · } and the claim numbers {N t , t = 1, 2, · · · }, respectively. Afterwards, we simulate the surplus process (2) starting at U 0 = u, with the premium amounts and claim sizes following the given exponential distributions. These simulations are replicated n = 3000 times; then, the trajectories with negative values (i.e., ruin event occurs) are counted, and we denote this number by n 1 . Hence, the infinite-time ruin probability ψ(u) can be estimated byψ In addition, because of the fact that U t → +∞ with probability one as t → +∞ when the relative safety loading condition holds, we know that U t will never become negative when t is large enough. Therefore, it is necessary for us to choose a suitable T st at which we should stop the simulated surplus process for each sample path if the ruin event does not occur before this time. As a consequence, (56) is actually the estimate of the finite-time ruin probability ψ(u, T st ) = P(T ≤ T st |U 0 = u). In this paper, we set T st = 1000. In practice, we can choose larger values for T st so that the bias of the estimator for ψ(u) is less significant. In Table 2 and Figure 2, we compare the simulated ruin probability with the approximated ruin probability. As can be seen, when u grows, both of the ruin probabilities approach zero. However, as alternatives to the true ruin probabilities, the approximations do not work well when the values of u are small. We can explain these results with the following three reasons. Firstly, as the limit of ψ(u) as u → +∞, e −Ru may be very different than ψ(u) at the beginning. Secondly, the simulated infinite-time ruin probabilities are indeed the estimated values for the finite-time ruin probability ψ(u, T st ), which are smaller than the true values of ψ(u). Thirdly, the total number of simulated trajectories n and the chosen time T st affect the simulated results. We could increase n and T st to improve the performance, but a longer run time is needed.
On the other hand, the values of simulated ruin probability and the approximated ruin probability become closer and closer with the increase in u, implying that the approximation method could work better as u grows. To strengthen this statement, we define γ(u) =ψ (u) e −Ru , and then, calculate the values of γ(u) with respect to different values of u; the results are listed in the last column of Table 2. It can be seen that γ(u) approaches 1 asymptotically, which indicates that it is valid to take e −Ru as the approximated result forψ(u) and, furthermore, for ψ(u) when u is large enough. Figure 3 also illustrates this conclusion visually. In practice, an insurer is always required to hold a huge number of initial surplus to guarantee its solvency under certain regulatory frameworks; therefore, the approximation method is of importance and is applicable in the risk management of insurance.

Example 2.
We suppose that the gain amount X is distributed by the exponential distribution with mean of 1/µ 1 , and the claim size Y follows the Pareto distribution, which has shape parameter τ 1 and scale parameter τ 2 , i.e., the distribution function F Y (y) is given by To perform the calculations, we set λ 1 = 1, µ 1 = 1, λ 2 = 0.1, τ 1 = 3, τ 2 = 16, and α = β = 0.5. It is not difficult to check that these values satisfy the relative safety loading condition (26). Our goal is to compare the asymptotic result λ 2 (1 + β)tF Y (y) (AS for simplification) with the simulated results of the finite-time ruin probabilities obtained using the Monte Carlo method (MC for simplification). As can be seen from Table 3, the ratio of MC to AS becomes closer and closer to one as t increases for different u, indicating that the asymptotic formula stated in Theorem 4 is valid and applicable in practice.

Conclusions
In this paper, we examine a generalization of the classical discrete-time risk model of an insurance portfolio with stochastic premiums, using a Poisson INAR(1) process and a Poisson INMA(1) process to fit the temporal dependence among the premium numbers and the temporal dependence among the claim numbers, respectively. We give the explicit expression for the function satisfied by the Lundberg adjustment coefficient and find the Lundberg approximation formula for the infinite-time ruin probability. Furthermore, we discuss and analyze the impact of the two thinning parameters and manifest that the dependence structure in the model has a significant influence on the risk of the surplus process in an insurance company. When the claim sizes follow a class of heavy-tailed distributions, we establish the large deviations of the aggregate claims and investigate the asymptotic formula for the finite-time ruin probability. In the numerical examples, we use MATLAB to randomly draw the sample paths of the proposed surplus process and compute estimates of the true ruin probabilities corresponding to different values of u using the Monte Carlo method. From the simulated results, it can be seen that the approximation formula and asymptotic formula we obtained are effective. Furthermore, these two formulas are much simpler to use for calculating and estimating the ruin probabilities than the Monte Carlo method.
As for future work, we could implement the same methodology by applying the time series for count data with other distributed innovations or an arbitrary innovations' distribution. Generally, using the same approach as that in Lemma 1 and Lemma 2, we can extend (14) and (20) to respectively. Therefore, if we could derive the explicit expression of c(r), the properties of the solution to the equation c(r) = 0 can be discussed, and the adjustment coefficient can be obtained to measure the risk. Additionally, we could adopt some higher-order processes to make the insurance risk model much more practical and flexible. In this situation, it becomes more challenging to find the expressions of P M 1 +···+M t (s) and P N 1 +···+N t (s). As a consequence, there might be some difficulties in deriving c(r) and defining the adjustment coefficient for an insurance portfolio.
On the other hand, instead of fixing the distributions and the parameters to illustrate the results by simulation, we can use the real dataset to fit the distributions and obtain the statistical estimates of the parameters, so that the ruin problems of the risk model could be analyzed in a more scientific way.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this appendix, we explicate the motivation and rationale of our proposed model (2) by the following descriptions, in order to make this paper more accessible to the readers. Let us begin with the classical discrete-time Lundberg-Cramér risk model where U t corresponds to the surplus of an insurance portfolio at time t, with U 0 = u being the initial surplus; c being the constant premium income per period, and L t representing the aggregate claim amount in period t that is defined as in which N t denotes the number of claims and Y t,j is the size of the jth payment to the insured in period t. After recursively calculating, it is easy to see that the risk model (A1) can be rewritten as which is equivalent to the model (1) by denoting N 0 For simplicity, it is assumed that the claim numbers of different periods are independent in the Lundberg-Cramér risk model, i.e., the claim number process {N t , t = 1, 2, · · · } is a sequence of i.i.d. random variables, which is certainly not realistic. As a consequence, ref. [11] proposes some new discrete-time risk models, where the Poisson INMA(1) process and Poisson INAR(1) process are used to describe the dependence structures among the numbers of claims. That is to say, the claim number process {N t , t = 1, 2, · · · } satisfies N t = α • N t−1 + ε t , t = 2, 3, · · · , or N t = β • η t−1 + η t , t = 1, 2, · · · .
On the other hand, both of the above two types of risk models suppose that the premiums are collected with positive deterministic constant rate c, which is also lacks the ability of describing the real situation of insurance portfolio. As an alternative to this case, ref. [18] proposes the risk model with stochastic premiums that can be expressed as where L t is defined by (A2), and P t aggregates the premiums in period t that is defined as in which M t counts the number of individual income, and X t,k represents the amount of the kth premium income for the insurance portfolio in period t. In the risk model (A3), it should be noted that both the premium number process {M t , t = 1, 2, · · · } and the claim number process {N t , t = 1, 2, · · · } are supposed to be sequences of i.i.d. random variables. Therefore, the goal of this paper is to introduce the idea of [11] into the risk model with stochastic premiums by using time series for count random variables to fit the temporal dependence among {M t , t = 1, 2, · · · } and {N t , t = 1, 2, · · · }, respectively. Furthermore, considering some insurance practices (please see Remarks 2 and 3), we assume that {M t , t = 1, 2, · · · } constitutes a Poisson INAR(1) process that satisfies M t = α • M t−1 + ε t , t = 2, 3, · · · , and {N t , t = 1, 2, · · · } constitutes a Poisson INMA(1) process that satisfies N t = β • η t−1 + η t , t = 1, 2, · · · . Our proposed risk model can also be generalized in several aspects to make itself more flexible and applicable, as discussed in the Conclusions.