On non-negative solutions of SDDEs with an application to CARMA processes

This note provides a simple sufficient condition ensuring that solutions of stochastic delay differential equations (SDDEs) driven by subordinators are non-negative. While, to the best of our knowledge, no simple non-negativity conditions are available in the context of SDDEs, we compare our result to the literature within the subclass of invertible continuous-time ARMA (CARMA) processes. In particular, we analyze why our condition cannot be necessary for CARMA($p,q$) processes when $p=2$, and we show that there are various situations where our condition applies while existing results do not as soon as $p\geq 3$. Finally, we extend the result to a multidimensional setting.


Introduction
Many quantities, such as wind speeds or (local) volatility of assets, are non-negative and behave in a stationary manner, and thus any reasonable model for these phenomena should comply with such constraints. Barndorff-Nielsen and Shephard [2] advocated the use of the stationary Ornstein-Uhlenbeck process driven by a subordinator (or, equivalently, a non-negative Lévy process) (L t ) t∈R , that is, the unique stationary solution to for some λ ∈ (0, ∞). Since the solution of (1.1) is explicitly given by a convolution between a non-negative kernel t → e −λt and a non-negative random measure dL, it is automatically non-negative. However, due to the simplicity of the Ornstein-Uhlenbeck process, its autocorrelation function is Corr(X 0 , X h ) = e −λ|h| with the only parameter being the rate of decay λ; in particular, the autocorrelation function must be monotonically decreasing. Consequently, there has been a need for working with more flexible modeling classes. A particular popular one consists of the continuous-time ARMA (CARMA) processes, which (as the name suggests) is the natural continuous-time analogue to the discrete-time ARMA processes. A CARMA process (X t ) t∈R can be characterized as a moving average of the form where the Fourier transform of g : [0, ∞) → R is rational. Consequently, a CARMA process is made up of a Lévy process (L t ) t∈R , a denominator (autoregressive) polynomial P , and a numerator (moving average) polynomial Q. Here it is required that deg(P ) > deg(Q) and that the zeroes of P belong to {z ∈ C : Re(z) < 0}. For more details on CARMA processes, see [7,8,9,12] or Section 3. Another important class, which also extends the Ornstein-Uhlenbeck process, is the one formed by solutions of stochastic delay differential equations (SDDEs) of the form where φ is a signed measure. Such equations have, for instance, been studied in [3,11,13]. Under suitable conditions, which will be stated in Section 2, there exists a unique stationary solution to (1.4) and it is, as well, a moving average of the form (1.3) with g being characterized through its Fourier transform. While both CARMA processes and solutions of SDDEs give rise to increased flexibility in the kernel g, and hence in the autocorrelation function, it is no longer guaranteed that it stays nonnegative on [0, ∞). This means that (X t ) t∈R is not necessarily non-negative although (L t ) t∈R is so, and one must instead place additional restrictions on (P, Q) and φ. A necessary and sufficient, but unfortunately rather implicit, condition ensuring nonnegativity of g is given by the famous Bernstein theorem on completely monotone functions [6]. In the context of CARMA processes, more explicit sufficient conditions were provided by [1,19] as they argued that a CARMA process is non-negative if the zeroes of the associated polynomials P and Q are real and negative, and if they respect a certain ordering (see (3.4)). To the best of our knowledge, this is the only available easy-to-check condition. In addition to the fact that this condition is not able to identify all non-negative CARMA processes, the Fourier transform of the driving kernel g of the solution of an SDDE (1.4) is most often not rational, meaning that the condition is not applicable in such setting. In this paper we show that, in order for the unique solution (X t ) t∈R of (1.4) to be non-negative, it is sufficient that φ is a non-negative measure when restricted to (0, ∞), that is, This result is presented in Section 2, in which we also give various examples and illustrate through simulations that simple violations of (1.5) results in solutions which can indeed go negative. Furthermore, in Section 3 we exploit the relation between SDDEs and invertible CARMA processes (that is, CARMA processes whose associated moving average polynomial Q only has zeroes in {z ∈ C : Re(z) < 0}) to establish conditions ensuring that a CARMA process is non-negative. Specifically, we observe that it is sufficient to show complete monotonicity of the rational function R/Q instead of Q/P on [0, ∞), where R is the (negative) remainder polynomial obtained from division of P with Q. This result is compared to the findings of [1,19] and it is shown that the two approaches identify different, but overlapping, regions of parameter values for which the CARMA(3, 2) process is non-negative. Finally, in Section 4 we extend the theory to the multivariate SDDEs (as introduced in [3]). In particular, when leaving the univariate setting we find that, in addition to a condition similar to (1.5) for a matrix-valued signed measure φ, one needs to require that λ := −φ({0}) is a so-called M-matrix. Section 5 contains proofs of the stated results and a couple of auxiliary lemmas. Before turning to the above-mentioned sections, we devote a paragraph to introduce relevant notations as well as essential background knowledge.
Preliminaries The models that we will consider are built on subordinators or, equivalently, non-negative Lévy processes. Recall that a real-valued stochastic process (L t ) t≥0 , L 0 ≡ 0, is called a one-sided subordinator if it has càdlàg sample paths and its increments are stationary, independent, and non-negative. These properties imply that the law of (L t ) t≥0 (that is, of all its finite-dimensional marginals) is completely determined by that of L 1 , which is infinitely divisible and, thus, by the Lévy-Khintchine formula and [10, Proposition 3.10]. Here γ ∈ [0, ∞) and ν is a σ-finite measure on (0, ∞) with ∞ 0 (1∧x) ν(dx) < ∞ (ν is a Lévy measure). For any subordinator (L t ) t≥0 the induced pair (γ, ν) is unique and, conversely, given such pair, there exists subordinator (which is unique in law) inducing this pair. We define a two-sided subordinator (L t ) t∈R from two independent one-sided subordinators (L 1 t ) t≥0 and (L 2 t ) t≥0 with the same law by Examples of subordinators include the gamma Lévy process, the inverse Gaussian Lévy process, and any compound Poisson process constructed from a sequence of non-negative i.i.d. random variables. We may also refer to multivariate subordinators, which are simply Lévy processes with values in R d , d ≥ 2, and entrywise non-negative increments (see [16] for details). We will say that a real-valued set function µ, defined for any Borel set of [0, ∞), is a signed measure if µ = µ + − µ − for two singular finite Borel measures µ + and µ − on [0, ∞). Note that the variation |µ| := µ + + µ − of µ is a measure on [0, ∞) and that integration with respect to µ can be defined in an obvious manner for any function f : [0, ∞) → R which is integrable with respect to |µ|.
In the last part of the paper we will consider a multivariate setting, and to distinguish this from the univariate one, we shall denote a matrix A with bold font and, unless stated otherwise, refer to its (j, k)-th entry by A jk . Finally, integration of matrix-valued functions against matrix-valued signed measures (that is, matrices whose entries are signed measures) can be defined in an obvious manner by means of the usual rules for matrix multiplication.

Stochastic delay differential equations and nonnegative solutions
Let (L t ) t∈R be a two-sided subordinator with a non-zero Lévy measure ν satisfying A stochastic process (X t ) t∈R is said to be a solution of the corresponding SDDE if it is stationary, has finite first moments (that is, E[|X 0 |] < ∞), and By (2.2), we mean that the equality holds almost surely for each fixed pair (s, t) ∈ R 2 with s < t. We will often write the SDDE in differential form as Set C + := {z ∈ C : Re(z) ≥ 0}. Properties such as existence and uniqueness of solutions of (2.2) are closely related to the zeroes of the function h φ : C + → C given by In case φ has bounded support, existence and uniqueness of solutions of (2.2) were established in [11,13] (under even milder conditions on (L t ) t∈R ). As we will later translate our findings into the framework of CARMA processes, which correspond to a particular class of delay measures φ with unbounded support, we will rely on the following result of [3]: Connor et al. [3]). Let h φ be given as in (2.3) and assume that h φ (z) = 0 for all z ∈ C + . Then the unique solution of (2.2) is given by Remark 2.2. Although it might sometimes be challenging to show that h φ (z) = 0 when z ∈ C + , as required by Theorem 2.1, a necessary and easy-to-check condition is that φ([0, ∞)) < 0. Indeed, this follows from the facts that To get an idea of which stationary processes that can be generated from (2.2), we provide a couple of examples where the condition of Theorem 2.1 can be checked.
where δ 0 is the Dirac measure at 0 and e 1 is the first canonical basisvector of R p . Note that e 1 is merely used as a normalization: the effect of replacing e 1 by an arbitrary vector c ∈ R p can be incorporated in b and A. By the assumption on σ(A), all the entries of e −At are exponentially decaying as t → ∞ and, thus, |φ| is a finite measure with moments of any order (in particular, is the unique pair of real polynomials Q, R : C → C such that (i) Q is monic and has no zeroes on C + , (ii) deg(Q) > deg(R), and (iii) Q and R have no common zeroes. Consequently, it follows from Theorem 2.1 that, as long as there exists a unique stationary solution (X t ) t∈R to (2.2), and it is given by In other words, (X t ) t∈R is a causal and invertible CARMA process with autoregressive polynomial P and moving average polynomial Q (see Section 3 or [8, Remark 4]). Conversely, given a moving average with g : [0, ∞) → R characterized by (2.7) for some polynomials P and Q having no zeroes in C + , and which satisfy deg(P ) = deg(Q) + 1, one can choose a unique constant λ ∈ R such that the polynomial R(z) := (z +λ)Q(z)−P (z) meets deg(R) < deg(Q). Thus, it follows that such process constitutes a stationary solution of the SDDE (2.2) with a delay measure φ of the form (2.6). Indeed, this is due the fact that a function f : [0, ∞) → R with a rational Laplace transform R/Q can always be represented as Example 2.4 (Discrete delay). Let λ, ξ ∈ R and τ ∈ (0, ∞), and assume that |ξ| ≤ τ −1 . Consider the following SDDE written in differential form: To show existence of a unique stationary solution using Theorem 2.1, we must argue that z + λ − ξe −zτ = 0 whenever z ∈ C + or, equivalently, that the two equations x + λ − ξe −τ x cos(τ y) = 0 and y + ξ sin(τ y) = 0 (2.9) cannot hold simultaneously if x ≥ 0 and y ∈ R. Since the only real solution of an equation of the form u = sin(αu) is u = 0 when |α| ≤ 1, the second equation in (2.9) implies that y = 0. If this is the case we have that In view of Remark 2.2, we conclude that h φ (z) = 0 for all z ∈ C + if and only if ξ < λ.
Turning to the question of whether a given solution (X t ) t∈R to (2.2) is nonnegative, we observe initially that, since it will necessarily take the form (2.4), there exist a number of equivalent statements for non-negativity: Let h φ be given as in (2.3) and assume that h φ (z) = 0 for all z ∈ C + . Furthermore, let (X t ) t∈R and g φ be defined through (2.4) and (2.5), respectively. The following statements are equivalent: (ii) X t ≥ 0 almost surely for all t ∈ R.
(iii) g φ is non-negative almost everywhere.
Remark 2.6. Let the setting be as in Theorem 2.5. The notion of almost sure nonnegativity is the best possible given that the process (X t ) t∈R itself is only defined for each fixed t ∈ R up to a set of probability zero. However, by (2.2), (X t ) t∈R can always be chosen to have càdlàg sample paths, in which case the property X t ≥ 0 will hold across all t ∈ R outside a set of probability zero.
While Theorem 2.5 tells that it is sufficient to show non-negativity of g φ , the kernel is often not tractable-not even when φ is rather simple. To give an example where this strategy does indeed work out, note that the solution of (2.8) with λ > 0 and ξ = 0 is the Ornstein-Uhlenbeck process, so g φ (t) = e −λt (by (2.5) as h φ (z) = z + λ) and, thus, proves that the solution is non-negative. In contrast, as soon as ξ = 0, g φ cannot be explicitly determined since its structure depends on the infinitely many solutions of the equation z + λ − ξe −τ z = 0 for z ∈ C (see [11,Lemma 2.1] for details). The following result, which relies on part (iv) of Theorem 2.5, provides a sufficient and simple condition on the delay φ which ensures that the solution is non-negative.
In light of Theorem 2.7, when one is trying to model non-negative processes, it is natural to write (2.2) as and then search for non-negative measures η on (0, ∞) satisfying η((0, ∞)) < λ (by Remark 2.2 this inequality must hold if h φ (z) = 0 for all z ∈ C + ). Here we use the convention ∞ 0 := (0,∞) . Remark 2.8. As will, for instance, appear from Example 3.3, the condition (2.11) is not necessary for (X t ) t∈R to be non-negative. According to Theorem 2.5 it is necessary and sufficient that (2.10) is satisfied. The case n = 0 is always true when h φ (z) = 0 for all z ∈ C + , and by relying on Faà di Bruno's formula it can be checked that . Moreover, α j refers to the jth entry of α ∈ N n 0 and |α| = n j=1 α j . From this expression it is easy to see that is satisfied if η is non-negative, since then L η is completely monotone and, hence, each term of the sum in (2.13) will be non-negative. It does, however, also show that η cannot be "too negative" if the solution (X t ) t∈R is required to be non-negative; for instance, the restriction for n = 1 implies in particular that Unfortunately, the complete set of restrictions implied by (2.10) and (2.13) seem to be very difficult to analyze.
In Figure 1 we simulate the stationary solution of (2.12) when η = ξδ τ , δ τ being the Dirac measure at τ , with the specific values λ = τ = 1 and ξ = 0.2. For comparison we rerun the simulation with ξ = −0.8. Note that existence and uniqueness of the stationary solution is ensured by Example 2.4 in both cases. As should be expected, it appears from the simulations that the solution stays nonnegative when ξ = 0.2, but when ξ = −0.8 (where non-negativity is not guaranteed by Theorem 2.7), the solution eventually becomes negative. The latter observation can, for instance, be proved theoretically by checking that

Define the companion matrix
Alternatively, the kernel g(t) = b ⊤ e At e p can be characterized in the frequency domain by the relation Due to the well-known form X t = t −∞ e A(t−s) e p dL s of the stationary Ornstein-Uhlenbeck process with drift parameter A and driven by (e p L t ) t∈R (see, for example, [17,18]), (3.1) shows immediately that (Y t ) t∈R admits the following state-space representation: The intuition behind any of the above (equivalent) definitions of the CARMA process is that (Y t ) t∈R should be the solution of the formal differential equation where D = d dt denotes derivative with respect to t. (For instance, by heuristically computing the Fourier transform of (Y t ) t∈R from (3.3), one can deduce that Y t should indeed take the form (3.1).) In general, a wide range of theoretical and applied aspects of CARMA processes have been studied in the literature: for further details, see the survey of Brockwell [8] and references therein.
Concerning non-negativity of CARMA processes driven by subordinators, the results are less conclusive; we briefly review existing results here. The simplest CARMA process, obtained with P (z) = z + λ for λ > 0 and Q(z) = 1, is the Ornstein-Uhlenbeck process, which (as noted in previous sections) is always nonnegative when the driving noise (L t ) t∈R is a subordinator. A particularly appealing property of CARMA processes, in contrast to general solutions of SDDEs, is that the Fourier transform (3.2) of the kernel g is rational. To be specific, with α 1 , . . . , α p and β 1 , . . . , β q being the zeroes of P and Q, respectively, Ball [1] showed that g is non-negative if α 1 , . . . , α p , β 1 , . . . , β q ∈ {z ∈ C : Re(z) < 0, Im(z) = 0} where the latter condition assumes that the zeroes are ordered such that α 1 ≥ · · · ≥ α p and β 1 ≥ · · · ≥ β q . This result was later complemented by [19], which gives two conditions (one necessary and one sufficient) for CARMA(p, 0) or, simply, CAR(p) processes to be non-negative. Furthermore, it is argued that (3.4) is both necessary and sufficient to ensure non-negativity of CARMA(2, 1) processes. For many purposes, but not all, CARMA(p, p − 1) processes are the most useful in practice; for instance, CARMA(p, q) processes have differentiable sample paths when q < p − 1 (see [12,Proposition 3.32]). When q = p − 1 and Q has no zeroes on C + (invertibility), we saw in Example 2.3 that the corresponding causal and invertible CARMA process (Y t ) t∈R can be identified as the unique solution of an SDDE of the form (2.12) with η(dt) = f (t) dt and f : [0, ∞) → R being characterized by Here R(z) := (z + λ)Q(z) − P (z) and λ ∈ R is uniquely determined by the condition deg(R) < q; it is explicitly given by Note that, in case the zeroes β 1 , . . . , β q of Q are distinct, Cauchy's residue theorem implies that

6)
Q ′ being the derivative of Q (see also [8,Remark 5]). In view of Theorem 2.7, it follows that (Y t ) t∈R is non-negative if f is so. To put it differently, while it is necessary and sufficient that Q/P is completely monotone on [0, ∞) for (Y t ) t∈R to be non-negative, it is sufficient to verify complete monotonicity of R/Q, a rational function where both the numerator and denominator are of lower order than the original one. We state this finding in the following result.
Theorem 3.1. Let (Y t ) t∈R be a causal and invertible CARMA(p, p − 1) process associated to the pair (P, Q), and let R be given as above. Then (Y t ) t∈R is nonnegative if R/Q is completely monotone on [0, ∞) or, equivalently, if the function f characterized by (3.5) is non-negative.
Remark 3.2. In situations where q < p − 1, Theorem 3.1 can still be useful for obtaining non-negative CARMA(p, q) processes as it can be applied in conjunction with the results of [19]. To be specific, suppose that (i) Theorem 3.1 applies to the CARMA(q + 1, q) process associated with a given the pair (P 1 , Q), and (ii) one of the sufficient conditions of [19,Theorem 1] applies to the CAR(p−q −1) process associated with a given autoregressive polynomial P 2 .
Then the CARMA(p, q) process associated with the pair (P 1 ×P 2 , Q) is non-negative (× denoting multiplication). This follows immediately from the fact that its kernel is the convolution between the kernels of the two moving averages associated to (i) and (ii).
While the first restriction in (3.8) is necessary, the second is not. In general, we see that (3.7) and (3.8) differ the most when |α − β| is small. In Figure 2 we have marked the feasible area for γ implied by (3.7) and (3.8), respectively, for a given polynomial P with zeroes α, β ∈ (−∞, 0). While the condition of Theorem 3.1 falls short in the context of CARMA(2, 1) processes, Corollary 3.4 below shows that it extends the class of CARMA(p, p − 1) processes which is known to be non-negative as soon as p ≥ 3. Specifically, this result gives sufficient conditions ensuring that CARMA (3,2) processes are non-negative in situations where those of [1,19] do not (and the latter are, to the best our knowledge, the only available easy-to-check conditions in the literature). One could aim for a similar result for a general p ≥ 3, but this would involve more complicated expressions.
In particular, if either (i) or (ii) is satisfied, then (Y t ) t∈R is non-negative.
(3.9) (Here, as in (3.4), the zeroes are ordered such that α 1 ≥ α 2 ≥ α 3 and β 1 ≥ β 2 .) Comparing (3.9) to the conditions obtained in Corollary 3.4, it follows that they are quite different. In particular, Corollary 3.4 does not require that α 1 , α 2 , and α 3 are real numbers, and even if this is the case, the conditions of this result are sometimes met while those of (3.9) are not, and vice versa. To give a simple example of the former case, suppose that α 1 > α 2 = α 3 and β 1 = β 2 . If β 1 is larger than α 2 , but smaller than the local extremum x * of P located in the interval (α 2 , α 1 ), it follows that both P (β 1 ) ≤ 0 and P ′ (β 1 ) ≤ 0. It is easy to check that x * = (2α 1 + α 2 )/3, so we conclude that the conditions of Corollary 3.4 are satisfied while those of (3.9) are not whenever On the other hand, a simple example where (3.9) applies, but Corollary 3.4 does not, is when α 3 < β 1 < α 2 , since this implies P (β 1 ) > 0.

Multivariate extensions
In this section we extend the results of Section 2 to multivariate SDDEs or, in short, MSDDEs which were studied in [3]. To be specific, for a given d ≥ 2 we let A stochastic process X t = (X 1 t , . . . , X d t ) ⊤ , t ∈ R, is a solution of the corresponding MSDDE if it is stationary, has finite first moments (that is, E[ X 0 ] < ∞), and satisfies for j = 1, . . . , d. In line with (2.7) we use the decomposition φ = −λδ 0 + η, where λ := −φ({0}) and η := φ( · ∩ (0, ∞)). If we define the convolution (µ * f )(t) = (f * µ)(t) := f (t − s) µ(ds) between a signed measure and a real-valued function f (assuming that the integral exists) and extend the definition to matrix-valued quantities by the usual rules of matrix multiplication, (4.1) can be written compactly as dX t = −λX t dt + (η * X)(t) dt + dL t , t ∈ R.
Among stationary processes which can be identified as solutions of (4.1) are the multivariate Ornstein-Uhlenbeck process (η ≡ 0) and, more generally, multivariate CARMA processes [12] (η(dt) = f (t) dt with f being of exponential type as in (2.6)). For further details, we refer to [3]. In general, solutions of MSDDEs are related to a function which is similar to (2.3), namely

2)
I d being the d × d identity matrix. In particular, it was shown in [3] that if det(h φ (z)) = 0 for all z ∈ C + , the unique solution (X t ) t∈R is given by In Theorem 4.1, which is a generalization of Theorem 2.7, we give sufficient conditions for the solution (X t ) t∈R to have non-negative entries. Before formulating this result we recall the definition of the so-called M-matrices: A ∈ R d×d is called an M-matrix if it can be expressed as A = αI d − B for some α ∈ [0, ∞) and B ∈ R d×d with non-negative entries and spectral radius ρ(B) ≤ α. The main point in introducing these is their relation to entrywise non-negativity of the matrix exponential e −At for all t ≥ 0 (see Lemma 5.2). For further details and other characterizations of M-matrices, we refer to [5,14] and references therein.
Theorem 4.1. Let h φ be defined as in (4.2) and assume that det(h φ (z)) = 0 for all z ∈ C + . Suppose further that all the entries of η are non-negative measures and that λ is an M-matrix. Then the entries of the unique solution (X t ) t∈R to (2.2), which is given by (2.4), are non-negative in the sense that X j t ≥ 0 for j = 1, . . . , d almost surely for each fixed t ∈ R.
Remark 4.2. When comparing Theorem 4.1 to Theorem 2.7, it seems that we have to impose an additional assumption of λ being an M-matrix, which was not needed in the univariate case. However, requiring that the one-dimensional quantity λ = −φ({0}) is an M-matrix would naturally be interpreted as a non-negativity constraint, but this is automatically satisfied under the stated assumptions on h φ , cf. Remark 2.6.
In the same way as we relied on the findings of Section 2 to obtain conditions for a CARMA process to be non-negative in Section 3, Theorem 4.1 can be used in conjunction with the relation between MSDDEs and multivariate CARMA processes outlined in [3] to obtain similar conditions in the multivariate setting.

Proofs
Proof of Theorem 2.5. The equivalence between (i) and (ii) is an immediate consequence of stationarity of (X t ) t∈R . The fact that (iii) implies (ii) is rather obviousthis is also readily seen by representing the subordinator (L t ) t∈R as a cumulative sum of its positive jumps plus a non-negative drift. To see that (ii) implies (iii), note by [15] that X t is infinitely divisible with a Lévy measure ν X given by where Leb( · ) refers to the Lebesgue measure. Since X t ≥ 0 almost surely, ν X is concentrated on (0, ∞), and hence we deduce that Since (L t ) t∈R is a subordinator and ν is non-zero, ν((0, ∞)) > 0. By combining this observation with (5.1) we conclude that g φ is non-negative almost everywhere.
In view of Bernstein's theorem on monotone functions [6], which states that g φ is non-negative if and only if its Laplace transform is completely monotone on [0, ∞), we have that (iii) holds if and only if (iv) does, and this completes the proof.
Proof of Theorem 2.7. This follows immediately from the discussion in Remark 2.8 or by observing that x → h φ (x) is a Bernstein function, and hence its reciprocal x → 1/h φ (x) must be completely monotone (see [4,Theorem 5.4]).
By considering t = 0 and t → ∞ once again we deduce that f is non-negative on [0, ∞) if and only if P (β) ≤ 0 and P ′ (β) ≤ 0, which concludes the proof.
Before turning to the proof of Theorem 4.1 we will need a couple of auxiliary lemmas. In relation to the formulation of the first result we recall the notation for a matrix-valued function f = (f jk ) and a matrix-valued signed measure µ = (µ jk ) of suitable dimensions and such that the involved integrals are well-defined.
Proof. To lighten notation, and in line with Remark 2.8, we denote the Laplace transform by L. Specifically, for a given signed measure µ and an integrable function , it is easy to check that L[g φ ] and L[φ] commute, and hence For a fixed s ∈ [0, ∞), it thus follows from (5.2) and (5.3) that for z ∈ C with Re(z) > 0. After rearranging terms we obtain By multiplying L[g φ ](z) from the left on both sides of (5.4) and using uniqueness of the Laplace transform, this shows that for t > s, and this completes the proof.
The following lemma, which will be used in the proof of Theorem 4.1 as well, presents a key property of M-matrices. While the result is well-known, we have not been able to find a proper reference, and hence we include a small proof of the statement here.
Lemma 5.2. Let A ∈ R d×d . If A is an M-matrix, then e −At has non-negative entries for all t ≥ 0. Conversely, if e −At has non-negative entries for all t ≥ 0 and the spectrum σ(A) of A belongs to {z ∈ C : Re(z) > 0}, then A is an invertible M-matrix.
Proof. If A is an M-matrix, we have that e −At = e −αt e Bt for some α ∈ [0, ∞) and B ∈ R d×d with non-negative entries. Since e Bt has non-negative entries, we conclude that the same holds for e −At . If instead the entries of e −At are non-negative and σ(A) ⊆ {z ∈ C : Re(z) > 0}, it follows that the entries of Proof of Theorem 4.1. For ε ∈ (0, ∞) define the measure φ ε := −λδ 0 +η(· ∩(ε, ∞)).
We start by observing that, when ε is sufficiently small, det(h φε (z)) = 0 for all z ∈ C + . (5.5) To see this suppose, for the sake contradiction, there exist sequences (ε n ) n≥1 ⊆ (0, ∞) and (z n ) n≥1 ⊆ C + such that det(h φε n (z n )) = 0 for all n ≥ 1 and ε n → 0 as n → ∞. Note that the length of any each entry of L[φ ε ](z) is bounded by the constant max j,k=1,...,d |φ jk |([0, ∞)), for all ε and z, and hence inf ε | det(h φε (z))| ∼ |z| d as |z| → ∞ by the Leibniz formula (the notation ∼ means that the ratio tends to one). This shows in particular that the sequence (z n ) n≥1 must be bounded, and so it has a subsequence (z n k ) k≥1 which converges to a point z * ∈ C + . However, this would imply that det(h φ (z * )) = lim k→∞ det(h φε n k (z n k )) = 0, which contradicts the original assumption that det(h φ (z)) = 0 for all z ∈ C + . Consequently, the condition (5.5) is satisfied as long as ε is smaller than a certain threshold, say, ε * . Thus, for ε ≤ ε * we can define the corresponding function g φε : [0, ∞) → R d×d through (4.4) (with φ replaced by φ ε ). By for t ∈ [0, ε]. By uniqueness of solutions of such differential equation it follows that for g φε (t) = e −λt for t ∈ [0, ε], in which case we can rely on Lemma 5.2 to deduce that its entries are non-negative. Now, given that g φε (t) has non-negative entries for t ∈ [0, kε] for some positive integer k ≥ 1, we can apply Lemma 5.1 with s = kε to deduce that this remains true when t ∈ (kε, (k + 1)ε]. Consequently, it follows by induction that g φε (t) has non-negative entries for all t ∈ [0, ∞).