A MAP(2) is fully described by the two transition rate matrices (D₀, D₁) given in terms of six rate parameters (λ₁₁, λ₁₂, λ₂₁, λ₂₂, σ₁, σ₂), i.e.

where (λ₁, λ₂) = (λ₁₁+λ₁₂, λ₂₁+λ₂₂). Each off - diagonal rate parameter of D₀ represents the transition rate out of a given phase into another phase but without invoking an arrival whereas each rate parameter of D₁ represents the transition rate invoking an arrival with or without changing phase. Let Q = D₀ + D₁ and P = (-D₀)^-1D₁. Then, Q is an infinitesimal generator for the continuous time Markov chain whereas P is the transition probability matrix for the embedded discrete time Markov chain. Let p be the stationary probability vector for P, i.e. pP = p and pe = 1 where e is a vector of ones.

Let T_i be the i-th stationary interval of a MAP. Then, the marginal moments and the lag-1 joint moments are obtained as follows

respectively. In order to simplify the notation, we use reduced moments defined as

A MAP(2) can be completely described by three marginal moments (r₁, r₂, r₃) and one joint moment r₁₁. In fact, a MAP(n) is completely described by (2n-1) marginal moments and (n-1)² lag-1 joint moments; see Casale et al. (2010) and Telek et al. (2007) for a minimal set of moments for MAP(n)s.

It was shown in Kim (2016) that the lag-1 joint Laplace transform (LT) of stationary intervals of MAP(n)s can be written in terms of n² parameters. Especially for MAP(2)s, the marginal and lag-1 joint LTs can be written in terms of (a₀, a₁, b₁, c₁₁) as follows

where (a₀, a₁, b₁, c₁₁) = -D0, Trace D0,pD1e,pD12e. A MAP(2) can be completely described by Eq. (2.2) which is given in terms of four coefficients (a₀, a₁, b₁, c₁₁).

Another minimal representation used in our study is the Jordan representation given in two matrices (E, R); see Telek et al. (2007) and Kim (2020) for details. A MAP(n) can be represented by two n × nmatrices (E, R) given in terms of n² parameters. The diagonal entries of the matrix E are eigenvalues of (-D₀)^-1. The matrix R accounts for lag-1 correlation and satisfies Re = e; see Telek et al. (2007) for details.

For a MAP(2), let (ν₁, ν₂) be the eigenvalues of (-D₀)^-1 with ν₁ ≥ ν₂. For the case of distinct eigenvalues, (E, R) can be written in terms of (ν₁, ν₂, ν₃, ν₄), i.e.

If ν₁ = ν₂, the E is written as

The marginal moments and the lag-1 joint moments are obtained as

where v = (ν₄/(ν₃+ν₄), ν₃/(ν₃+ν₄)). Again, by matching moments (ν₁, ν₂, ν₃, ν₄) can be determined by (r₁, r₂, r₃, r₁₁).

If (-D₀)^-1 has distinct eigenvalues, then (-1/ν₁,-1/ν₂) can be determined as roots of the characteristic polynomial equation given as s²+a₁s+a₀ = 0. Assuming, WLOG, that ν₁ > ν₂, we have

and

by moment matching based on Eq. (2.3). Again, by moment matching, the following equations can be obtained to write the marginal and joint LTs in (2.1) and (2.2) in terms of (ν₁, ν₂, ν₃, ν₄)

By inversion of the LT in (2.1) and (2.2) given in terms of (ν₁, ν₂, ν₃, ν₄), we get the following marginal and lag-1 joint density functions

where

Note that (ϕ₁, ϕ₂) = (ϕ₁₁+ϕ₁₂, ϕ₂₁+ϕ₂₂). Depending on the sign of ν₃, the marginal distribution of a stationary interval is either hyper-exponential or mixed generalized Erlang (MGE). In order to simplify the notation, let u₁ (x) and u₂ (x) be the exponential density function each with mean ν₁ and ν₂ respectively, i.e.

Also, let u₁₂ (x) be the density function of the hypo-exponential distribution that is a sum of two independent exponential random variables each with mean ν₁ and ν₂, i.e.

Below, we denote this hypo-exponential distribution by Hypoexp(1/ν₁, 1/ν₂).

If (-D₀)^-1 has identical eigenvalues, i.e. ν₁ = ν₂, then we have a12=4a0 in Eq. (2.1) and we have ν₁ = a₁/2 and

along with

Let u₁₁ (x) be the density function of the Erlang(2, 1/ν₁) distribution that is a sum of two independent and identical exponential random variables with mean ν₁ i.e.

Then, by inverse LT, the marginal LT in Eq. (2.1) is converted into the following marginal density function

where (ψ₁, ψ₁₁) = (1-ϕ₁/ν₁, ϕ₁/ν₁). That is,

The joint LT given in Eq. (2.2), by inverse LT, is converted into the following joint density function.

where

By the definition of conditional probability, the conditional density function is given as

by which the conditional distribution of Y given X is either exponential or Erlang. That is,

or if X ~ Erlang(2,1/ν₁), then the conditional distribution of Y is

Consider a MAP(2) with the following set of moments (r₁, r₂, r₃, r₁₁) = (3/7, 11/63, 13/189, 103/567) for which we have (a₀, a₁, b₁, c₁₁) = (9, 6, 15/7, 31/7), (ν₁, ν₂) = (1/3, 1/3), and (ν₃, ν₄) = (-1/2, 4/3). We have (ϕ₁, ϕ₂) = (19/27, 2/27) and the marginal distribution is MGE(2) which can be generated as follows

by Eq. (3.9). Since (ψ_1,1, ψ_1,12, ψ_12,1, ψ_12,12) = (31/63, 14/63, 14/63, 4/63), the next interval Y can be generated conditional to the distribution of X, i.e.

or

by Eqs. (3.10) and (3.11). The next interval is generated conditional to the distribution of Y.