{'fraction_non_alphanumeric': 0.042439510547588724, 'fraction_numerical': 0.013011461159404242, 'mean_word_length': 4.441030275643922, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Backward Simulation of Multivariate Mixed Poisson Processes

The simulation of correlated multivariate Poisson processes with negative correlation between their components has many important applications in Finance, Insurance, Geophysics, and many other areas of applied probability. Introduced in our earlier work, the Backward Simulation (BS) approach to the simulation of correlated multivariate Poisson processes is able to capture a wide range of correlation, including extreme positive and extreme negative correlation, that is not possible with other approaches such as the forward simulation approach. Moreover, the BS approach enables simple and efficient generation of sample paths of correlated multivariate Poisson processes. In this work, we extend the BS approach to multivariate mixed Poisson processes.


Introduction
The simulation of dependent Poisson processes is an important problem having many applications in Insurance, Finance, Geophysics and many other areas of applied probability-see [1,2,4,5,6,7,14,27,28,32] and references therein. For example, in Operational risk, estimating the losses resulting from operational events requires the simulation of multivariate Poisson processes that must be calibrated to historical correlation matrices of operational events; see [13]. The Poisson and Negative Binomial processes are some of the most popular underlying models amongst practitioners for describing the operational losses of the business units of a financial organization and the moments of claim arrivals in the insurance industry [21,25]. Dependence between Poisson processes can be achieved by various operations applied to independent processes. One of the most popular approaches, often considered in actuarial modeling, is the Common Shock Model (CSM) [21,29,32], where a third Poisson process is used to couple two independent processes. For example, let (ν (1) t , ν (2) t , ν (3) t ) be three independent Poisson processes with intensities (λ 1 , λ 2 , λ 3 ), each defined as for j = 1, 2, 3 where the T (j) i , exponentially distributed random variables, are the arrival moments corresponding to the j th process. Through superposition, we can obtain two correlated Poisson processes N The correlation coefficient between the Poisson processes N (1) t and N (2) t , having intensities λ = λ 1 + λ 2 and µ = λ 2 + λ 3 , in the CSM satisfy The latter relation immediately implies 0 ≤ ρ ≤ min(λ, µ) max(λ, µ) .
It is clear that, in such a model, negative correlations cannot be obtained and that correlations are constant in time. The extreme correlation problem was considered in [18] and [24] where an optimization problem for the joint distribution was solved numerically. The problem was reduced to that of random vectors having specified marginal distributions in [20], where the Extreme Joint Distribution (EJD) method, a pure probabilistic, efficient, and rather simple algorithm to find the joint distributions with extreme correlations applicable to any discrete marginal probability distribution was proposed. Connections with some classical results obtained in [16], [19], and [33] were also discussed. The Backward Simulation (BS) method, in conjunction with the EJD method, was developed in order to address the restrictions in the correlation structure of multivariate Poisson processes constructed using classical approaches such as the CSM. The Backward Simulation approach, considered in [13], allows for a wider range of correlations, both positive and negative, to be attained in comparison to the CSM. Moreover, BS allows for a dynamic correlation structure versus static correlation in the CSM; specifically, it is a linear function, in time, of the terminal correlation, ρ(T ), at the end of the simulation interval [0, T ] [20]. There are two general approaches to the simulation of multivariate Poisson processes-Forward and Backward simulation. The Forward approach consists of generating exponentially distributed inter-arrival times until the simulation time is at or past the simulation interval [0, T ]. Our Backward approach is based on exploiting the conditional uniformity of Poisson processes-we first construct a joint distribution satisfying the marginal distributions with the desired correlation structure at the terminal simulation time T and then generate the corresponding number of arrival moments using the conditional uniformity of the arrival times 1 . This is one of the major advantages of the Backward approach-only the ability to sample from a suitable joint distribution at the terminal time is required; the arrival moments of the multivariate Poisson process can be generated uniformly in a coordinate-wise manner. The BS approach was extended in [20] to the class of bivariate processes containing both Poisson and Wiener components. It also led to the introduction of the Forward Continuation (FC) of Backward Simulation, a method for extending the process simulated by BS to subsequent intervals [nT, (n+ 1)T ], where n is some integer, that preserves the joint distribution at various grid points nT [9].
The EJD method enables the construction of joint distributions that exhibit extreme dependence between the components; in other words, the EJD method constructs extreme joint distributions that extremize ρ(T ). Extreme joint distributions are used to generate extreme admissible correlations, from which all correlations within the admissible range can be obtained. It was extended to the multivariate setting in [9].
In this paper, we extend the Backward Simulation approach for Poisson processes to the class of Mixed Poisson processes (MPPs), which are a natural generalization of the class of Poisson processes that can be represented as a Poisson process with a random intensity [17]. Our contribution is a method of simulation for multivariate mixed Poisson processes such that 1. any desired correlation that is admissible at the terminal simulation time can be matched, 2. the correlation structure is a function of time.
Moreover, we describe the time structure of correlations for this class of processes and analyze the Forward Continuation of the Backward Simulation in the finite interval. This approach allows us to extend the model to arbitrary times.
The plan of this paper is as follows. In Section 2, we briefly discuss the basics of MPPs. Section 3 extends the BS method for Poisson processes to multivariate MPPs, allowing us to fill in our process back to time 0. The extension to MPPs is first discussed in the bivariate setting. Section 4 briefly reviews the EJD method, necessary for the construction of joint distributions needed at the terminal time. We also discuss in Section 4 how to sample from Extreme Joint Distributions. In Section 5, we extend the FC approach for Poisson processes to MPPs. This allows us to propagate the process forward in time to some possibly infinite horizon. Finally, we make some concluding remarks in Section 6.

Mixed Poisson Process
We begin by reviewing some properties of MPPs. The main results of the theory of MPPs can be found in [17]. Recent results on the characterization of the multivariate MPP are in [34]. We consider a counting process with arrival moments 0 < T 1 < · · · < T i < · · · , where 1(·) is the indicator function. Also, for convenience, we let T 0 = 0. The classical Poisson process is defined as a process with independent increments such that the inter-arrival times between the events ∆T i := T i − T i−1 form a sequence of independent identically distributed random variables having an exponential distribution with parameter λ. It is well known that the number of events of X t in the interval [0, t] has the Poisson distribution with parameter λt: A natural generalization of the Poisson distribution is to randomize the intensity parameter λ, leading to the Mixed Poisson Distribution (MPD). [17]). A discrete random variable X is said to be mixed Poisson distributed, MP(U), with structure distribution U, if

Definition 2.1 (Mixed Poisson Distribution
where Λ is a random variable distributed according to U . The structure distribution U can be viewed as a prior distribution, which allows us to view (2) as a conditional distribution, given a realization of the intensity parameter Λ = λ and (3) as an unconditional distribution. Another interpretation of (3) is that it is a mixture of Poisson distributions.
The MPP is a Poisson process with a non-negative random intensity.
Lundberg [22] also showed that there exists a MPP for each structure distribution U and that the process is uniquely defined.
In what follows, we denote by MPP(U ), the class of MPPs with structure distribution U . It is not difficult to see that if X t ∈ MPP(U ), then the probability generating function takes the form and

Conditional distribution of arrival moments
It is well known that the intervals ∆T i = T i − T i−1 of a Poisson process with intensity λ form a sequence of independent, exponentially distributed random variables: This forms the basis of the forward approach to the simulation of Poisson processes. Given n events to be generated and a positive intensity λ, one can sequentially generate exponentially distributed intervals, ∆T i , and determine the arrival moments, However, there is an alternative approach based on the fundamental property of the conditional distribution of the arrival moments [10]. Let T = {T 1 , T 2 , . . . , T n } be a sequence of n independent random variables having a uniform distribution in the interval [0, T ]: Denote by τ k the kth order statistic of T, (k = 1, 2, . . . , n): Theorem 2.3. The distribution of the arrival moments of a Poisson process, X t , with finite intensity in the interval [0, T ] conditional on the number of arrivals, X T = n, coincides with the distribution of the order statistics: The converse statement was proved in [20]: If a process, X t , is represented as a random sum where {T k } N k=1 are independent, identically distributed random variables having a uniform conditional distribution, P(T k ≤ t | N ) = tT −1 , k = 1, 2, . . . , N in the interval [0, T ] and the random variable N ∼ Pois(λT ), then X t is a Poisson process with intensity λ in the interval [0, T ].
This result leads to the BS algorithm for the multivariate Poisson processes considered in [13] and [20]. In Section 3, we generalize Proposition 2.4 for the class of MPPs, which can be obtained by a similar construction using the random variable N ∼ MP(U ).

Negative Binomial Process
Let us now consider the Negative Binomial (NB) process, which is a MPP with the structure distribution U being the gamma distribution. The Negative Binomial process is widely used for count data that exhibit overdispersion because, unlike the Poisson process, it does not have the restriction that its mean must equal its variance. For this reason, we use the Negative Binomial distribution in numerical experiments in Section 3.2 and Section 5.1 to compare the processes generated by BS in the mixed Poisson versus the Poisson case.
The generating function of the Negative Binomial process is given in the following Lemma, the proof of which can be found in standard texts [15].
where b corresponds to the probability of success and r corresponds to the number of failures until the process is stopped. Remark 1. Notice that our process is not a Lévy process-the inter-arrival times are not independent but only conditionally independent.

Backward Simulation of Mixed Poisson Processes
Backward Simulation of Poisson processes, studied in [13] and [20], relies on the conditional uniformity of the arrival moments. BS requires sampling the corresponding joint distribution, at terminal time, to obtain a vector of the number of events for each coordinate in the simulation interval [0, T ]. The dependency structure manifests itself in the joint distribution (Section 4 discusses how correlated joint distributions can be obtained), i.e., in the sampled vector of terminal events. Each coordinate is simulated independently by drawing the corresponding number of uniform variates, which are then ordered to give the arrival moments of events.

Fundamental Results
Let us now show that the distribution of the arrival moments, conditional on the number of events, is also uniform for MPPs. Moreover, we show that the process generated by BS remains a MPP. First, we introduce the following two lemmas, the proofs of which can be found in [20]. To this end, we introduce some useful notation. Let X t be a MPP and consider the points with non-negative integer coordinates, k j ≥ 0, we denote the norm of the vector by For any d-dimensional vector, x = (x 1 , x 2 , . . . , x d ), with non-negative coordinates, and k ∈ Z d + , we denote The conditional probability of the number of events in each non-overlapping interval given the number of terminal events takes the form [17] with the multinomial coefficient Consider a discrete random variable, ξ, taking non-negative integer values with probabilities, p k = P(ξ = k), k = 0, 1, 2, . . . , and denote its generating function byp Then, for any fixed x ∈ [0, 1], the sequence {q k (x)} is a probability distribution and its generating function,q(z; Lemma 3.2. Consider a discrete random variable, ξ, taking non-negative integer values with probabilities, p k = P(ξ = k), k = 0, 1, 2, . . . , and denote its generating function byp With Lemma 3.2, the vector analogue of Lemma 3.1, we can prove the main theorem of this section.
Theorem 3.3. Let the process X t be represented as a random sum where the number of random events N ∼ MP(U ) and {T k } N k=1 are independent, identically distributed random variables having a uniform conditional distribution in the interval Proof. We prove the following two statements.

At any time
The increments of the process X t over disjoint intervals are conditionally independent random variables.
The theorem follows immediately from the two results above. Let us prove the first statement. As noted in (4), the generating function of X T is In our case,p(z) = ∞ 0 e xT (z−1) dU (x). Thus, the probabilities q k = p k (t) := P (X t = k | X T ). Taking x = tT −1 in (9), we obtain from (12) and Lemma 3.1 that as was to be proved.
The second statement is proved using Lemma 3.2. Let x = (x 1 , x 2 , . . . , x d ), satisfying the conditions listed in the statement of the lemma.
Then from (8), we have Now consider the generating function The latter relation implies that the increments of X t are conditionally independent, as was to be proved. Theorem 3.3 is intuitively appealing. Indeed, X t is a Poisson process with random intensity, λ, which is determined at time t = 0 and, therefore, measurable with respect to the filtration {F t } t≥0 generated by the process X t . The conditional distribution of the arrival moments is uniform in the interval [0, T ] and does not depend on the parameters of the process.
Algorithm 1 describes the Backward Simulation of multivariate MPPs (MMPPs) in detail. A correlated multivariate mixed Poisson process, N t , has as its marginals MPPs. Since the marginals are correlated and the joint distribution does not factorize, a joint distribution that has the desired correlation structure with the marginalized distributions satisfying the constraints of the given marginals is needed (Step 1 of Algorithm 1). This is discussed in Section 4. Given a vector of the counts of the number of events from the joint distribution, Theorem 3.3 applies to each marginal distribution independently.

Time Structure of Correlations
Let us now analyze the time structure of the correlations of multivariate MPPs generated by BS. Since correlations are inherently pairwise in nature, the analysis carried out in the bivariate setting corresponds to pairs of variables in the multivariate setting.
Theorem 3.4 (Time Structure of the Correlation Coefficient). Consider a bivariate process (X t , Y t ) such that X t and Y t possess the conditional uniformity property. The sample paths of the processes X t and Y t are generated by BS in the interval [0, T ]. Let the correlation coefficient where and σ 2 (X t ) denotes the variance of X t .
Proof. First, we show that the generating function of the process, To this end, note that for 0 ≤ m ≤ k and 0 ≤ n ≤ l, since at the end of the simulation interval T there are k events in total for X T and l events in total for Y T , the probability of the number of events m and n by a certain time t can be viewed as a Bernoulli trial with probability of success (t/T ). Taking expectation, we obtain Denote P kl = P (X T = k, Y T = l). Then we find Equation (15) is derived. Differentiatingĝ(t, z, w) twice, we find from which we obtain Theorem 3.4 is thus proved.
In the Poisson case, the auxiliary function Z(T )/Z(t) in (14) reduces to tT −1 . Thus, the correlation structure is linear in time in the simulation interval [0, T ]. This is not true in general for MPPs. For example, for the Negative Binomial processes, the auxiliary functions take the form The graph of the correlation function is presented in Figure 1, where the good agreement of the theoretical and empirical results can be seen.

Forward vs Backward Simulation and their Correlation Boundaries
Given a simulation interval [0, T ], stochastic processes are usually simulated forwards in time. This is due to the fact that it is conceptually natural and technically simpler to do so. However, it is not always the most suitable choice. This can be seen in the Forward Simulation (FS) of a correlated bivariate Poisson process (N k } k≥0 denotes the sequence of inter-arrival times for process i ∈ {1, 2}. Forward Simulation of counting processes like Poisson processes consists of repeatedly simulating the inter-arrival times {∆T The sequence of inter-arrival times represents a sample path of the counting process. In the Poisson case, the inter-arrival times are exponentially distributed, We must rely on the Fréchet-Hoeffding Theorem 3 in [16,19] to induce dependence between the marginal distributions of the inter-arrival times in the case of FS, which gives us the relations to obtain extremal positive dependence between the distributions of the inter-arrival times and to obtain extremal negative dependence. We claim that the relations (18) and (19) lead to extreme correlations of the process under FS. In the case of extremal positive dependence, (18) implies that Define κ = µ 1 /µ 2 . Obviously, 0 < κ < ∞. We show that for all t > 0, κt .
Suppose that N

This implies that T
(2) m ≤ κt < T (2) m+1 which in turns implies that N t ). The latter relation implies that Similar reasoning in the case 0 < κ < 1 leads to This allows us to compare the correlations obtained from the FS case to the correlations obtained in the BS case. From the last two equations above, we can see that the notion of extreme dependence obtained via the Fréchet-Hoeffding theorem results in a correlation coefficient that is a function of the intensities. This is very restrictive and precludes the possibility of calibrating to data. In contrast, the BS approach allows for the construction of processes with any desired correlation that is within the range of admissible correlations. Further comparisons of the Forward vs the Backward approaches can be found in [8,20].

Backward Simulation and Extreme Joint Distributions
We showed in Section 3 that the conditional uniformity property holds for the class of Mixed Poisson processes and that the process resulting from Backward Simulation with the number of events at terminal simulation time T generated by a MPD is indeed a MPP. Backward Simulation for the class of Mixed Poisson processes relies on the knowledge of the joint MPD at terminal time T , but how do we construct a multivariate MPD with some desired dependency structure in the first place? In this section, we briefly review the work in [9] and [20] in order for this paper to be self-contained. Moreover, some details are explained more clearly here than in [9] and [20]. We address the general problem of constructing multivariate joint distributions from given marginal distributions such that the linear correlation coefficient between the marginals are equal to some desired correlations. We also discuss how to sample from such multivariate joint distributions.

The Bivariate Case
We begin by describing the main ideas in 2-dimensions to build some intuition before presenting the general d-dimensional case. To that end, suppose we have a discrete bivariate distribution P with marginals Q (1) and Q (2) . Clearly, the admissible linear correlation coefficient C between the marginals is bounded by some maximum attainable correlationĈ (1) and some minimum attainable correlationĈ (2) . Moreover, every admissible correlation C, can be represented as a convex combination of the extreme correlations for some w ∈ [0, 1]. The extreme correlations are clearly extreme points. Extreme Measures in the bivariate case are defined as follows subject to (2) j = 1. Extremize denotes either max or min and the objective function is where P ij = P(X 1 = i, X 2 = j).
For completeness, we mention that the infinite dimensional LP 4 (23) is a Monge Kantorovich Problem (MKP). This aspect of the problem is not immediately relevant to us; we refer to standard references such as [31] for more details.
The solution to (23) is an Extreme Joint Distribution that determines the Extreme Measureŝ P (1) andP (2) which have a one-to-one relationship to the extreme correlationsĈ (1) andĈ (2) [8]. The Extreme Measures (23) lead to extreme correlations since extremizing the bivariate expectation extremizes the linear correlation coefficient as can be seen in (24). Moreover, let where w is the solution of (22) andP (1) andP (2) are the Extreme Measures having extreme correlationsĈ (1) andĈ (2) , respectively. Then, it is not hard to show that P is a discrete bivariate probability distribution with marginals Q (1) and Q (2) and correlation coefficient C. This insight allows us to reduce the problem of calibration to a simpler problem of solving a linear equation. Thus, if extreme joint distributions can be computed, they can be used to generate extreme correlations (extreme points) to calibrate to the given correlation. If the calibration fails-there is no solution to the linear equation (22) with w ∈ [0, 1]-it implies that no bivariate process with the marginal distributions Q (1) and Q (2) and correlation C exists. In such an event, the assumptions of the parameter values (including the inference procedures to obtain them) and any raw data should be checked.

The General Case
The 2-dimensional case described in the previous subsection generalizes to the d-dimensional case (d > 2) described below. However, instead of dealing with a single correlation coefficient, in the general case, we consider a d × d correlation matrix C where C ij represents the linear correlation coefficient between marginal distributions Q (i) and Q (j) . Clearly, we now have to consider more general notions of extremal dependency. One concept of extremal dependency consistent with observations of correlations matrices is to consider only pairwise extremal dependence. That is, we consider pairwise monotonicity, which represents the strongest type of association between random variables and implies maximally positive (comonotonicity) and negative values (antimonotonicity) for the linear correlation coefficient; see [30] for more details on extremal dependence concepts in multivariate settings and [20] for details on monotonicity as it relates to distributions. In contrast to the bivariate case, there are n = 2 d−1 extreme correlation matriceŝ C (j) , which are also extreme points [9], each described by a monotonicity structure.
k , then marginal distributions Q (i) and Q (k) have a comonotone dependency relationship and an antimonotone dependency relationship otherwise.
Remark 2. Note that whether e (j) 1 is initially set to 0 or 1 does not matter and is an arbitrary choice.
Similar to the bivariate case, for each j = 1, 2, . . . , n, each extreme correlation matrixĈ (j) is associated with an Extreme MeasureP (j) , described by a monotonicity structure e (j) , as described below.
represents the k-th given marginal distribution and each objective function takes the form where, similarly, Remark 3. There are m = d(d − 1)/2 objective functions, where each h k,l (p) extremizes the dependency between a pair of coordinates.
The multi-objective program (27) is, in fact, a multi-objective multi-marginal MKP, the solutions of which determine Extreme Measures. Potential solutions of (27) are multivariate probability measures P which are tensors. Thus, the multi-objective problem is not only tedious to program but practically prohibitively expensive to compute (in terms of both time and storage) for moderate d [8]. One approach that leads to a computable solution to these infinite dimensional problems (23) and (27) is given by the Extreme Joint Distribution (EJD) Theorem, which gives a semi-analytic form describing completely the extreme joint distribution.
An accompanying algorithm, the EJD algorithm, provides an efficient numerical method to solve (27) by computing the extreme joint distributions in (29) and their corresponding supports; see [20] and [9] for more details. Note thatP (j) is very sparse in that mostP

Sampling from Multivariate Extreme Measures
There are two attractive features of the EJD approach which make sampling from multivariate Extreme Measures simple. The first is that Extreme MeasuresP (k) are monotone distributions [20]. Consequently, their support remains a graph in higher dimensions. This is very convenient for sampling as this means that Extreme Measures can be sampled from via the inverse CDF method. Second, any discrete multivariate probability measure with specified marginals and some desired dependency structure can be represented as a convex combination of Extreme Measures. That is, the one-to-one relationship between (22) and (25) extends to the multidimensional case as follows. We first find coefficients (w 1 , . . . , w n ) that satisfy w j ≥ 0 for j = 1, 2, . . . , n, n j=1 w j = 1 and C = w 1Ĉ (1) + · · · + w nĈ (n) (31) and then set P = w 1P (1) + · · · + w nP (n) (32) whereP (j) is the extreme measure satisfying the LP (27) and having the extreme correlation matrixĈ (j) . Since w i ≥ 0 for i = 1, 2, ..., n and n i=1 w i = 1, it follows immediately that P is a probability measure. Moreover, it follows from the linearity of sums that P has the correlation matrix C given on the left side of (31). In addition, since eachP (j) has marginal distributions Q (1) , . . . , Q (d) , it follows that P also has marginal distributions Q (1) , . . . , Q (d) .
Note that (31) can be converted to a constrained system of linear equations by flattening each extreme correlation matrixĈ (j) into a column vector A j ∈ R m where m = d(d − 1)/2. Since C and allĈ (j) are symmetric with 1s on their diagonal, this can be done by taking each row in the strictly upper triangular part of eachĈ (j) , appending them into a row vector and taking the transpose to be A j to obtain A = [A 1 , . . . , A n ] ∈ R m×n , representing the extreme points of our problem in correlation space. Similarly, we can flatten the correlation matrix C on the left side of (31) to a vector b ∈ R m . Then (31) and the constraints w j ≥ 0 for j = 1, 2, ..., n and b j=1 w j = 1 are equivalent to the constrained system of linear equations There are many possible solutions to the constrained system of equations (33). One approach is to choose a suitable objective function 5 and then use (33) as the constraints for an optimization problem with that objective function. However, if the goal is just to find any solution to (33), then a simpler approach is to reformulate (33) aŝ whereÂ is A with the row 1 T appended to the bottom of it andb is b with a 1 appended to the bottom of it. Then note that (34) has the form of the standard constraints for a Linear Programming Problem (LPP). Moreover, the first stage of many LPP codes finds a solution to (34). As explained in Section 13.5 of [26], one standard approach to finding a solution to (34) is to solve the LPP where z ∈ R m+1 and E is a (m + 1) × (m + 1) diagonal matrix such that E ii = +1 ifb i ≥ 0 and E ii = −1 ifb i < 0. Clearly, w = 0 and z = | b | satisfies the constraints (35b) and (35c). So, we can use w = 0 and z = | b | as a staring point for the simplex method to solve (35). It's clear from the constraint z ≥ 0 that the solution satisfies 1 T z ≥ 0. Moreover, if 1 T z = 0 then z = 0. Hence, (35) has a solution 1 T z = 0 if and only ifÂw =b, w ≥ 0 has a solution. Hence, the simplex method applied to (35) will find a solution to (33), if a solution exists. After obtaining {w j } n j=1 through calibration, as described above, the decomposition of the desired measure as a convex combination of Extreme Measures (32) provides an easy method for simulation. Since w j ≥ 0 for j = 1, 2, . . . , n and n j=1 w j = 1 we can view w j as the probability that a draw of P comes from the Extreme MeasureP (j) . Moreover, as noted at the beginning of this subsection, we can easily sample fromP (j) . Therefore, we can utilize methods of discrete random variate simulation (see [12] for a comprehensive exposition) to generate a random variable that has the distribution P .
Remark 4. Despite the fact that the problem size grows exponentially in d, due to the structure of the problem (33) and the fact that the simplex method needs to explicitly access m+1 columns ofÂ at a time (assuming you have some clever way to decide which new vector to bring into the active set at each step of the simplex method without explicitly accessing all the columns of A that are in the inactive set) the LP (35) can be solved for a surprisingly large d, e.g., d = 51, which corresponds to n = 2 50 ≈ 10 15 ; see [8] and [23].

Applications
The novelty in our Backward Simulation approach to the construction of multivariate mixed Poisson processes is that we are able to calibrate to a given correlation structure in the form of a correlation matrix and attain extreme correlations. One application of our BS approach is to Operational Risk [13,20] where it is necessary to compute regulatory capital that is allocated in the event of an operational loss. Generally, this is accomplished using Monte Carlo scenario generation, which requires calibration to historical data. The BS approach captures ranges of correlation that have been empirically observed (such as negative correlations ) but could not be captured before. The BS approach has also found application in the modeling of credit losses of loans under the default model with correlated defaults. The paper [PenH] emphasizes the importance of modelling correlated defaults (in particular, when the correlations are negative) in credit risk modelling.

Forward Continuation of the Backward Simulation for Mixed Poisson Processes
Up to this point we have discussed the construction of MPPs within some interval [0, T ]. A natural question to ask is, what if we want to simulate the process forwards in time, past the original simulation interval. One solution to this was introduced in the Poisson setting in [9], known as the Forward Continuation (FC) to Backward Simulation (BS), where the joint distribution was preserved at various future time points, with the interval in between filled in by BS. The main idea of the FC method is to retain the independent increments property. Note that since the notion of the linear correlation coefficient and our choice of multivariate extremal dependency concept is pairwise in nature, the discussion in the bivariate setting generalizes immediately to the multivariate case.
In the more general setting of MPPs, the conditional independence of the increments poses a challenge. We show that with the right construction, the arguments in the Poisson setting extend naturally to the MPP setting. Consider a sequence of time intervals [0, T ), [T, 2T ), . . . , [mT, (m + 1)T ). Suppose that a bivariate MPP (X t , Y t ) has been simulated in [0, T ) using BS and that we wish to continue forward the process (X T +τ , Y T +τ ) = (X T + ∆X τ , Y T + ∆Y τ ) in [T, 2T ). At the end of the interval [T, 2T ), draw a new independent version of (X T , Y T ) and add it to the original (X T , Y T ) to obtain (X 2T , Y 2T ) at time 2T , i.e, take (X 2T , Y 2T ) d = (X T , Y T ). For 0 ≤ τ < T , the process (X T +τ , Y T +τ ) can be constructed by conditional uniformity given the number of events in [T, 2T ), i.e by Backward Simulation. This retains the independence of the marginal increments since X T is independent of ∆X τ , and similarly for Y T , yet preserves the joint distribution of the bivariate process (X t , Y t ) even though, in the MPP setting, the marginal increments are only conditionally independent. Note that X T +τ remains a MPP due to our construction and the superposition property of MPPs [17].

Forward Time Structure of Correlations
Now that we have a method to extend a bivariate MPP defined initially within the interval [0, T ) to an interval [mT, (m+1)T ), it is natural to analyze the behavior of the correlation coefficient on each interval. It turns out that we can prove asymptotic stationarity of the correlation coefficient under Forward Continuation. Indeed, the covariance of the processes X t and Y t at time T + τ can be written as Cov since X T is independent of ∆Y τ and vice versa. From (17), we have: Dividing each side by σ(X T +τ ) σ(Y T +τ ), we obtain the correlation coefficient .
Passing to the limit as m → ∞ in the standard manner, we obtain that lim m→∞ ρ(mT + τ ) = ρ(T ) as was to be proved.
A graphical illustration of Theorem 5.1 can be found in Figure 2, which shows the good agreement between the analytic and numerical results. Note that while the correlations at the calibrated integer grid points nT are exact, the correlation structure in between grid points, generated via the FC method, requires a few time periods in order to settle to the asymptotic value of ρ(T ).
The first few time periods can be used as "burn in" periods to achieve a more stable and accurate desired value for the correlation between processes for the in-between time periods as filled in by BS. This can be also used to construct bivariate MPPs that display a constant time structure of correlation at some given correlation value instead of exhibiting the linear behavior in the interval [0, T ]. Note that, depending on the choice of simulation, the process obtained by simulating purely through BS on [0,2T) will have a different correlation structure than a process simulated using BS on [0,T) and FC on [T,2T). This is a user-defined choice that is dependent on the needs of the application.

Concluding remarks
In this paper, we extended the Backward Simulation (BS) method and the Forward Continuation of the BS method from the class of Poisson processes to the more general class of multivariate Mixed Poisson Processes. The advantages of the Backward approach over the Forward approach for generating sample paths of multivariate Mixed Poisson processes in some simulation interval [0, T ] are numerous: simple and efficient simulation in d-dimensions; specification of a dependency structure in the form of a given correlation matrix C at terminal simulation time T ; a wider range of possible correlations between the marginal distributions.
The Backward Simulation approach is applicable to any process that exhibits the order statistic property. For example, Backward Simulation is applicable to the class of Negative Binomial Lévy Processes [3]. In fact, it is applicable to a more general class of processes known as Compound Poisson Processes which also posses a linear correlation structure under BS. It is also applicable to the inhomogeneous Poisson Processes. Extending the BS approach to Compound Poisson Processes and inhomogeneous Poisson Processes will be the subject of our forthcoming work.