This note provides a simple sufficient condition ensuring that solutions of stochastic delay differential equations (SDDEs) driven by subordinators are nonnegative. While, to the best of our knowledge, no simple nonnegativity 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≥3. Finally, we extend the result to a multidimensional setting.

CARMA processescomplete monotonocitynonnegative stationary processesstochastic delay differential equationssubordinators60G1060G1760H0560H10This work was supported by the Danish Council for Independent Research (grants 4002-00003 and 9056-00011B).Introduction

Many quantities, such as wind speeds or (local) volatility of assets, are nonnegative 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 nonnegative Lévy process) (Lt)t∈R, that is, the unique stationary solution to
dXt=−λXtdt+dLt,t∈R,
for some λ∈(0,∞). Since the solution of (1) is explicitly given by
Xt=∫−∞te−λ(t−s)dLs,t∈R,
a convolution between a nonnegative kernel t↦e−λt and a nonnegative random measure dL, it is automatically nonnegative. However, due to the simplicity of the Ornstein–Uhlenbeck process, its autocorrelation function is Corr(X0,Xh)=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 (Xt)t∈R can be characterized as a moving average of the form
Xt=∫−∞tg(t−s)dLs,t∈R,
where the Fourier transform of g:[0,∞)→R is rational. Consequently, a CARMA process is made up of a Lévy process (Lt)t∈R, a denominator (autoregressive) polynomial P, and a numerator (moving average) polynomial Q. Here it is required that deg(P)>deg(Q)\deg (Q)$]]> and that the zeroes of P belong to {z∈C:Re(z)<0}. For more details on CARMA processes, see [7–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
dXt=∫[0,∞)Xt−sϕ(ds)dt+dLt,t∈R,
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 (4) and it is, as well, a moving average of the form (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 (Xt)t∈R is not necessarily nonnegative although (Lt)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, 21] as they argued that a CARMA process is nonnegative if the zeroes of the associated polynomials P and Q are real and negative, and if they respect a certain ordering (see (22)). 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 nonnegative CARMA processes, the Fourier transform of the driving kernel g of the solution of an SDDE (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 (Xt)t∈R of (4) to be nonnegative, it is sufficient that ϕ be a nonnegative measure when restricted to (0,∞), that is,
ϕ(B∩(0,∞))≥0for all Borel measurable setsB.
This result is presented in Section 2, in which we also give various examples and illustrate through simulations that simple violations of (5) result 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 nonnegative. 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, 21] and it is shown that the two approaches identify different, but overlapping, regions of parameter values for which the CARMA(3,2) process is nonnegative. 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 (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, nonnegative Lévy processes. Recall that a real-valued stochastic process (Lt)t≥0, L0≡0, is called a one-sided subordinator if it has càdlàg sample paths and its increments are stationary, independent, and nonnegative. These properties imply that the law of (Lt)t≥0 (that is, all its finite-dimensional marginal distributions) is completely determined by that of L1, which is infinitely divisible and, thus,
logE[eiθL1]=iθγ+∫0∞(eiθx−1)ν(dx),θ∈R,
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)<∞ (in particular, ν is a Lévy measure). For any subordinator (Lt)t≥0 the induced pair (γ,ν) is unique and, conversely, given such pair, there exists a subordinator (which is unique in law) inducing this pair. We define a two-sided subordinator (Lt)t∈R from two independent one-sided subordinators (Lt1)t≥0 and (Lt2)t≥0 with the same law by
Lt=Lt11[0,∞)(t)−L(−t)−21(−∞,0)(t),t∈R.
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 nonnegative i.i.d. random variables. We may also refer to multivariate subordinators, which are simply Lévy processes with values in Rd, d≥2, and entrywise nonnegative increments (see [17] 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 mutually 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 Ajk. 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 (Lt)t∈R be a two-sided subordinator with a nonzero Lévy measure ν satisfying ∫1∞xν(dx)<∞ or, equivalently, E[L1]<∞. Moreover, let ϕ be a signed measure on [0,∞) with
∫0∞t2|ϕ|(dt)<∞.
A stochastic process (Xt)t∈R is said to be a solution of the corresponding SDDE if it is stationary, has finite first moments (that is, E[|X0|]<∞), and
Xt−Xs=∫st∫[0,∞)Xu−vϕ(dv)du+Lt−Ls,s<t.
By (7), we mean that the equality holds almost surely for each fixed pair (s,t)∈R2 with s<t. We will often write the equation in differential form as in (4).

The condition that a stationary solution (Xt)t∈R must satisfy E[|X0|]<∞ is imposed in order to ensure that the integral on the right-hand side of (7) can be defined in the usual way. Indeed, in this case
∫st∫[0,∞)|Xu−v|ϕ(dv)du<∞
almost surely since its mean is finite. Note also that, if (7) has any such solution (Xt)t∈R, then E[L1]=−ϕ([0,∞))E[X0]. Consequently, our assumption that E[L1]<∞ is innocent and only imposed to avoid trivial cases. If ϕ has bounded support, (7) makes perfect sense without the condition E[|X0|]<∞, and existence and uniqueness of solutions were established in [11, 13] under milder conditions on (Lt)t∈R. While it should be doable to study (7) under milder conditions when ϕ is unbounded as well, we leave this for future research.

Set C+:={z∈C:Re(z)≥0}. Properties such as existence and uniqueness of solutions of (7) are closely related to the zeroes of the function hϕ:C+→C given by
hϕ(z):=z−∫[0,∞)e−ztϕ(dt),z∈C+.
Note that hϕ is analytic on {z∈C:Re(z)>0}0\}$]]> and real-valued on {z∈C+:Im(z)=0}. 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 existence and uniqueness result of [3]:

(Basse-O’Connor et al. [<xref ref-type="bibr" rid="j_vmsta177_ref_003">3</xref>]).

Lethϕbe given as in (8) and assume thathϕ(z)≠0for allz∈C+. Then the unique solution of (7) is given byXt=∫−∞tgϕ(t−u)dLu,t∈R,wheregϕ:[0,∞)→Ris characterized by∫0∞e−itygϕ(t)dt=1hϕ(iy),y∈R.

The assumption that hϕ(z)≠0 for all z∈C+ as well as the form of the solution (9) are, in fact, rather natural. Indeed, heuristically speaking the equation (7) reads (AϕX)t=(DL)t with Aϕ being the linear operator given by (Aϕf)(t)=(Df)(t)+∫[0,∞)f(t−s)ϕ(ds) and D=ddt the derivative with respect to t. Since
∫Re−ity(Aϕf)(t)dt=hϕ(iy)∫Re−ityf(t)dt,y∈R,
the assumption on hϕ can be interpreted as an invertibility assumption on Aϕ, and its inverse Aϕ−1 takes the form (Aϕ−1f)(t)=∫−∞tgϕ(t−s)f(s)ds. While it might sometimes be challenging to show that hϕ(z)≠0 when z∈C+, as required by Theorem 1, a necessary and easy-to-check condition is that ϕ([0,∞))<0. Indeed, this follows from the facts that hϕ(0)=−ϕ([0,∞)), hϕ is continuous, and hϕ(x)→∞ as x→∞ (where x∈[0,∞)).

To get an idea of the stationary processes that can be generated from (7), we provide a couple of examples where the condition of Theorem 1 can be checked.

(CARMA processes).

Let λ∈R and, for a given p∈N, suppose b∈Rp and A∈Rp×p with a spectrum σ(A) contained in {z∈C:Re(z)>0}0\}$]]>. Consider the delay measure
ϕ(dt)=−λδ0(dt)+b⊤e−Ate1dt,
where δ0 is the Dirac measure at 0 and e1 is the first canonical basis vector of Rp. Note that e1 is merely used as a normalization: the effect of replacing e1 by an arbitrary vector c∈Rp 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, (6) is satisfied). The function hϕ takes the form hϕ(z)=z+λ−b⊤(A+zIp)−1e1, where Ip is the p×p identity matrix. By the fraction decomposition it follows that b⊤(A+zIp)−1e1=R(z)/Q(z) for z∈C+, where (Q,R) 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)\text{deg}(R)$]]>, and (iii) Q and R have no common zeroes. Consequently, it follows from Theorem 1 that, as long as
P(z):=(z+λ)Q(z)−R(z)≠0for allz∈C+,
there exists a unique stationary solution (Xt)t∈R to (7), and it is given by (9) with gϕ satisfying
∫0∞e−itygϕ(t)dt=Q(iy)P(iy),y∈R.
In other words, (Xt)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
Xt=∫−∞tg(t−s)dLs,t∈R,
with g:[0,∞)→R characterized by (12) 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 a process constitutes a stationary solution of the SDDE (7) with a delay measure ϕ of the form (11). Indeed, this is due the fact that a function f:[0,∞)→R with a rational Laplace transform R/Q can always be represented as f(t)=b⊤e−Ate1 for suitable b∈Rdeg(Q) and A∈Rdeg(Q)×deg(Q) with σ(A)⊆{z∈C:Re(z)>0}0\}$]]>. For more on the relation between solutions of SDDEs and CARMA processes, see [3, Section 4.3].

(Discrete delay).

Let λ,ξ∈R and τ∈(0,∞), and assume that |ξ|≤τ−1. Consider the following SDDE written in differential form:
dXt=(−λXt+ξXt−τ)dt+dLt,t∈R.
To show existence of a unique stationary solution using Theorem 1, we must argue that z+λ−ξe−zτ≠0 whenever z∈C+ or, equivalently, that the two equations
x+λ−ξe−τxcos(τy)=0andy+ξe−τxsin(τy)=0
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 (14) implies that y=0. If this is the case we have that
x+λ−ξe−τxcos(τy)≥λ−ξ,
since u↦u−ξe−τu is increasing on [0,∞). In view of the last part of Remark 2, we conclude that hϕ(z)≠0 for all z∈C+ if and only if ξ<λ.

Turning to the question of whether a given solution (Xt)t∈R to (7) is nonnegative, we observe initially that, since it will necessarily take the form (9), there exist a number of equivalent statements for nonnegativity; these are presented in Theorem 2 below. Before formulating the result, recall that a function f:(0,∞)→R is called completely monotone if it is infinitely differentiable and
(−1)ndndtnf(t)≥0for alln∈N0andt>0.0\text{.}\]]]>
If f is nonnegative and infinitely differentiable, and its derivative Df is completely monotone, it is called a Bernstein function. A convenient property of such function is that f∘g is completely monotone as long as g is so (in fact, this property characterizes the class of Bernstein functions). For further details on completely monotone functions, Bernstein functions, and their relations, see [20].

Lethϕbe given as in (8) and assume thathϕ(z)≠0for allz∈C+. Furthermore, let(Xt)t∈Randgϕbe defined through (9) and (10), respectively. The following statements are equivalent:

Xt≥0almost surely for somet∈R.

Xt≥0almost surely for allt∈R.

gϕis nonnegative almost everywhere.

1/hϕis completely monotone on(0,∞).

Let the setting be as in Theorem 2. The notion of almost sure nonnegativity is the best possible given that the process (Xt)t∈R itself is only defined for each fixed t∈R up to a set of probability zero. However, by (7), (Xt)t∈R can always be chosen to have càdlàg sample paths, in which case the property Xt≥0 will hold across all t∈R outside a set of probability zero.

While Theorem 2 tells that it is sufficient to show nonnegativity 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 (13) with λ>00$]]> and ξ=0 is the Ornstein–Uhlenbeck process, so gϕ(t)=e−λt (by (10) as hϕ(z)=z+λ) and, thus, proves that the solution is nonnegative. 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, provides a sufficient and simple condition on the delay ϕ which ensures that the solution is nonnegative.

Lethϕbe defined as in (8) and assume thathϕ(z)≠0for allz∈C+. Suppose further that ϕ is a nonnegative measure when restricted to(0,∞), that is,ϕ(B∩(0,∞))≥0for all Borel measurable setsB.Then (i)–(iv) of Theorem2are satisfied.

In light of Theorem 3, when one is trying to model nonnegative processes, it is natural to write (7) as
dXt=−λXtdt+∫0∞Xt−sη(ds)dt+dLt,t∈R,
and then search for nonnegative measures η on (0,∞) satisfying η((0,∞))<λ (by the last part of Remark 2 this inequality must hold if hϕ(z)≠0 for all z∈C+). Here we use the convention ∫0∞:=∫(0,∞).

As will, for instance, appear from Example 3, the condition (16) is not necessary for (Xt)t∈R to be nonnegative. According to Theorem 2 it is necessary and sufficient that (15) is satisfied for 1/hϕ. The case n=0 is always true when hϕ(z)≠0 for all z∈C+, and by relying on Faà di Bruno’s formula (see, e.g., [16, Theorem 2]) it can be checked that
(−1)ndndxn1hϕ(x)=∑α∈N0n:∑j=1njαj=nn!|α|!(1−Lη(1)(x))α1α1!(1!)α1⋯αn!(n!)αnhϕ(x)|α|+1∏j=2n((−1)jLη(j)(x))αj
when n≥1, where Lη(j) denotes the jth derivative of the Laplace transform Lη(x):=∫0∞e−xtη(dt) of η:=ϕ(·∩(0,∞)). Moreover, αj refers to the jth entry of α∈N0n and |α|=∑j=1nαj. From this expression it is easy to see that (15) is satisfied if η is nonnegative, since then Lη is completely monotone and, hence, each term of the sum in (18) will be nonnegative. It does, however, also show that η cannot be “too negative” if the solution (Xt)t∈R is required to be nonnegative; for instance, the restriction for n=1 implies in particular that
∫0∞tη(dt)≥−1.
Unfortunately, the complete set of restrictions implied by (15) and (18) seem to be very difficult to analyze.

In Figure 1 we simulate the stationary solution of (17) 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 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 nonnegativity is not guaranteed by Theorem 3), the solution eventually becomes negative. The latter observation can, for instance, be proved theoretically by checking that
hϕ(x)3d2dx21hϕ(x)|x=0=ξ2+5ξ+2
when λ=τ=1, and hence 1/hϕ is not completely monotone on [0,∞) if ξ=−0.8.

Simulations of X1,…,X200 from the model (17) with η=ξδ1 when ξ=0.2 (left) and ξ=−0.8 (right). In both cases, λ=1

Nonnegative CARMA processes

Let P(z)=zp+a1zp−1+⋯+ap and Q(z)=b0+b1z+⋯+bq−1zq−1+zq be two real monic polynomials with p>qq$]]>, and assume that P has no zeroes on C+ (causality). Define the companion matrix
A=010⋯0001⋯0⋮⋮⋱⋱⋮00⋯01−ap−ap−1⋯−a2−a1,
and let b=(b0,b1,…,bq−1,1,0,…,0)⊤ and ep=(0,0,…,0,1)⊤ (both being elements of Rp). The causal CARMA(p,q) process (Yt)t∈R driven by the subordinator (Lt)t∈R, and associated to the autoregressive polynomial P and moving average polynomial Q, is defined by
Yt=∫−∞tb⊤eA(t−s)epdLs,t∈R.
Alternatively, the kernel g(t)=b⊤eAtep can be characterized in the frequency domain by the relation
∫0∞e−ityg(t)dt=Q(iy)P(iy),y∈R.
Due to the well-known form Xt=∫−∞teA(t−s)epdLs of the stationary Ornstein–Uhlenbeck process with drift parameter A and driven by (epLt)t∈R (see, for example, [19, 18]), (19) shows immediately that (Yt)t∈R admits the following state-space representation:
dXt=AXtdt+epdLt,Yt=b⊤Xt.
The intuition behind any of the above (equivalent) definitions of the CARMA process is that (Yt)t∈R should be the solution of the formal differential equation
P(D)Yt=Q(D)(DL)t,t∈R.
(For instance, by heuristically computing the Fourier transform of (Yt)t∈R from (21), one can deduce that Yt should indeed take the form (19).) 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 nonnegativity 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 λ>00$]]> and Q(z)=1, is the Ornstein–Uhlenbeck process, which (as noted in previous sections) is always nonnegative when the driving noise (Lt)t∈R is a subordinator. A particularly appealing property of CARMA processes, in contrast to general solutions of SDDEs, is that the Fourier transform (20) 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 nonnegative if
α1,…,αp,β1,…,βq∈{z∈C:Re(z)<0,Im(z)=0}and∑j=1kαj≥∑j=1kβjfork=1,…,q,
where the latter condition assumes that the zeroes are ordered so that α1≥⋯≥αp and β1≥⋯≥βq. This result was later complemented by [21], which gives two conditions (one necessary and one sufficient) for CARMA(p,0) or, simply, CAR(p) processes to be nonnegative. Furthermore, it is argued that (22) is both necessary and sufficient to ensure nonnegativity 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 1 that the corresponding causal and invertible CARMA process (Yt)t∈R can be identified as the unique solution of an SDDE of the form (17) with η(dt)=f(t)dt and f:[0,∞)→R being characterized by
∫0∞e−ityf(t)dt=R(iy)Q(iy),y∈R.
Here R(z):=(z+λ)Q(z)−P(z) and λ∈R is uniquely determined by the condition deg(R)<q; it is explicitly given by
λ=a1−bp−2=−αp+∑j=1p−1(βj−αj).
Note that, in case the zeroes β1,…,βq of Q are distinct, Cauchy’s residue theorem implies that
f(t)=−∑j=1qP(βj)Q′(βj)eβjt,t≥0,Q′ being the derivative of Q (see also [8, Remark 5]). In view of Theorem 3, it follows that (Yt)t∈R is nonnegative if f is so. To put it differently, while it is necessary and sufficient that Q/P is completely monotone on [0,∞) for (Yt)t∈R to be nonnegative, 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.

Let(Yt)t∈Rbe a causal and invertible CARMA(p,p−1) process associated to the pair(P,Q), and let R be given as above. Then(Yt)t∈Ris nonnegative ifR/Qis completely monotone on[0,∞)or, equivalently, if the function f characterized by (23) is nonnegative.

In situations where q<p−1, Theorem 4 can still be useful for obtaining nonnegative CARMA(p,q) processes as it can be applied in conjunction with the results of [21]. To be specific, suppose that

Theorem 4 applies to the CARMA(q+1,q) process associated with a given the pair (P1,Q), and

one of the sufficient conditions of [21, Theorem 1] applies to the CAR(p−q−1) process associated with a given autoregressive polynomial P2.

Then the CARMA(p,q) process associated with the pair (P1×P2,Q) is nonnegative (× 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).

Suppose that (Yt)t∈R is a causal and invertible CARMA(2,1) process, so that its polynomials P and Q are representable as P(z)=(z−α)(z−β) and Q(z)=z−γ for suitable α,β∈{z∈C:Re(z)<0} and γ∈(−∞,0). It follows from (22) that for (Yt)t∈R to be nonnegative it is necessary and sufficient that
Im(α)=Im(β)=0andγ≤max{α,β}.
On the other hand, from (24) we have f(t)=−P(γ)eγt, and so one would need to require that P(γ)=(γ−α)(γ−β)≤0 in order to apply Theorem 4. Imposing this condition is equivalent to assuming that
Im(α)=Im(β)=0andmin{α,β}≤γ≤max{α,β}.
While the first restriction in (26) is necessary, the second is not. In general, we see that (25) and (26) differ the most when |α−β| is small. In Figure 2 we have marked the feasible area for γ implied by (25) and (26), respectively, for a given polynomial P with zeroes α,β∈(−∞,0).

An autoregressive polynomial P with zeroes α,β∈(−∞,0) (gray) together with the feasible values for the zero γ∈(−0,∞) of the corresponding moving average polynomial Q implied by (25) (red: solid and dashed) and (26) (red: solid)

While the condition of Theorem 4 falls short in the context of CARMA(2,1) processes, Corollary 1 below shows that it extends the class of CARMA(p,p−1) processes which is known to be nonnegative as soon as p≥3. Specifically, this result gives sufficient conditions ensuring that CARMA(3,2) processes are nonnegative in situations where those of [1, 21] 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.

Let(Yt)t∈Rbe a causal and invertible CARMA(3,2) process with polynomials P and Q, and letβ1andβ2,Re(β1)≥Re(β2), denote the zeroes of Q. Then the function f characterized by (23) is nonnegative if and only ifIm(β1)=Im(β2)=0and one of the following conditions holds:

β1>β2{\beta _{2}}$]]>andP(β1)≤min{P(β2),0}.

β1=β2andmax{P(β1),P′(β1)}≤0.

In particular, if either (i) or (ii) is satisfied, then(Yt)t∈Ris nonnegative.

Let (Yt)t∈R be the causal and invertible CARMA(3,2) associated to real polynomials P and Q representable as
P(z)=(z−α1)(z−α2)(z−α3)andQ(z)=(z−β1)(z−β2),z∈C,
for some α1,α2,α3∈{z∈C:Re(z)<0} and β1,β2∈(−∞,0). If one seeks to apply the result of [1, 21] stated in (22), it must be required that
Im(αj)=0,β1≤α1,andβ1+β2≤α1+α2.
(Here, as in (22), the zeroes are ordered such that α1≥α2≥α3 and β1≥β2.) Comparing (27) to the conditions obtained in Corollary 1, it follows that they are quite different. In particular, Corollary 1 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 (27) are not, and vice versa. To give a simple example of the former case, suppose that α1>α2=α3{\alpha _{2}}={\alpha _{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 1 are satisfied while those of (27) are not whenever
α1+α22<β1≤α1+α22+α1−α26.
On the other hand, a simple example where (27) applies, but Corollary 1 does not, is when α3<β1<α2, since this implies P(β1)>00$]]>.

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 Lt=(Lt1,…,Ltd)⊤, t∈R, be a d-dimensional subordinator with E[‖L1‖]<∞ and ϕ=(ϕjk) a d×d matrix-valued signed measure on [0,∞) with
∫0∞t2|ϕjk|(dt)<∞,j,k=1,…,d.
A stochastic process Xt=(Xt1,…,Xtd)⊤, t∈R, is a solution of the corresponding MSDDE if it is stationary, has finite first moments (that is, E[‖X0‖]<∞), and satisfies
Xtj−Xsj=∑k=1d∫st∫[0,∞)Xu−vkϕjk(dv)du+Ltj−Lsj,s<t,
for j=1,…,d. In line with (3) 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, (28) can be written compactly as
dXt=−λXtdt+(η∗X)(t)dt+dLt,t∈R.
Among stationary processes which can be identified as solutions of (28) 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 (11)). For further details, we refer to [3]. In general, solutions of MSDDEs are related to a function which is similar to (8), namely
hϕ(z)=zId−∫[0,∞)e−ztϕ(dt),z∈C+,Id 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 (Xt)t∈R is given by
Xt=∫−∞tgϕ(t−s)dLs,t∈R,
where gϕ:[0,∞)→Rd×d is characterized by
∫0∞e−itygϕ(t)dt=hϕ(iy)−1,y∈R.
In Theorem 5, which is a generalization of Theorem 3, we give sufficient conditions for the solution (Xt)t∈R to have nonnegative entries. Before formulating this result we recall the definition of the so-called M-matrices: A∈Rd×d is called an M-matrix if it can be expressed as A=αId−B for some α∈[0,∞) and B∈Rd×d with nonnegative entries and spectral radius ρ(B)≤α. The main point in introducing these is their relation to entrywise nonnegativity of the matrix exponential e−At for all t≥0 (see Lemma 2). For further details and other characterizations of M-matrices, we refer to [5, 14] and references therein.

Lethϕbe defined as in (30) and assume thatdet(hϕ(z))≠0for allz∈C+. Suppose further that all the entries ofηare nonnegative measures and thatλis an M-matrix. Then the entries of the unique solution(Xt)t∈Rto (28), which is given by (9), are nonnegative in the sense thatXtj≥0forj=1,…,dalmost surely for each fixedt∈R.

When comparing Theorem 5 to Theorem 3, 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 nonnegativity constraint, but this is automatically satisfied under the stated assumptions on hϕ, cf. Remark 3.

In the same way as we relied on the findings of Section 2 to obtain conditions for a CARMA process to be nonnegative in Section 3, Theorem 5 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. For brevity, however, we will only address nonnegativity of the multivariate Ornstein–Uhlenbeck process and the multivariate CARMA(2,1) process in the following.

Let λ∈Rp×p with a spectrum σ(λ) contained in {z∈C:Re(z)>0}0\}$]]>, and consider the Ornstein–Uhlenbeck equation:
dXt=−λXtdt+dLt,t∈R.
The stationary solution (Xt)t∈R takes the form Xt=∫−∞te−λ(t−s)dLs and thus, to ensure that Xt is nonnegative (and without further restrictions on the driving Lévy process (Lt)t∈R), it must be required that all entries of e−λt are nonnegative for every t≥0. According to Lemma 2 this is the case if and only if λ is an M-matrix, which is exactly the condition of Theorem 5 for nonnegativity.

Let A1,A2,B0∈Rd×d be given such that both
z↦det(Idz2+A1z+A2)andz↦det(Idz+B0)
have no roots on C+ (causality and invertibility). Due to the same reasoning as in Section 3, the corresponding d-dimensional CARMA(2,1) process is naturally characterized by
Yt=∫−∞tg(t−s)dLs,t∈R,
where the Fourier transform of g:[0,∞)→Rd×d is
∫0∞e−ityg(t)dt=(−Idy2+A1iy+A2)−1(Idiy+B0),y∈R.
(See also [12, Theorem 3.22].) By [3, Theorem 4.8], (Yt)t∈R can also be identified as the unique solution to (29) with λ=A1−B0 and
η(dt)=e−B0tdt[B0(A1−B0)−A2].
In view of Theorem 5 and Lemma 2 we conclude that the entries of (Yt)t∈R are nonnegative if both B0 and A1−B0 are M-matrices, and all entries of the matrix B0(A1−B0)−A2 are nonnegative.

Proofs

The equivalence between (i) and (ii) is an immediate consequence of stationarity of (Xt)t∈R. The fact that (iii) implies (ii) is rather obvious—this is also readily seen by representing the subordinator (Lt)t∈R as a cumulative sum of its positive jumps plus a nonnegative drift. To see that (ii) implies (iii), note by [15] that Xt is infinitely divisible with a Lévy measure νX given by
νX(B):=(Leb×ν)({(s,x)∈[0,∞)×R:xgϕ(s)∈B∖{0}}),
where Leb(·) refers to the Lebesgue measure. Since Xt≥0 almost surely, νX is concentrated on (0,∞), and hence we deduce that
0=νX((−∞,0))≥Leb({s∈[0,∞):gϕ(s)<0})ν((0,∞)).
Since (Lt)t∈R is a subordinator and ν is nonzero, ν((0,∞))>00$]]>. By combining this observation with (33) we conclude that gϕ is nonnegative almost everywhere. In view of Bernstein’s theorem on monotone functions [6], which states that gϕ is nonnegative 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. □

This follows immediately from the discussion in Remark 4 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]). □

Since f takes the form (23) and deg(Q)=2, it follows from [21, Remark 2] that, necessarily, Im(β1)=Im(β2)=0 (as we also discussed in relation to (22)). If β1>β2{\beta _{2}}$]]>, we may use (24) and represent f explicitly as
f(t)=1β1−β2(P(β2)eβ2t−P(β1)eβ1t),t≥0.
Here we have also used the fact that Q′(β1)=−Q′(β2)=β1−β2. By considering t=0 and t→∞ we conclude that f(t)≥0 for all t if and only if both P(β1)≤P(β2) and P(β1)≤0, and this thus establishes (i). If β1=β2=:β, we observe initially that
R(z)=(z+a1+2β)(z−β)2−P(z)=2β3+a1β2−a3−(3β2+2βa1−a2)z=P′(β)β−P(β)−P′(β)z.
Since the function having Fourier transform y↦(iy−β)−2 is t↦teβt, it follows from (23) and the above expression for R that
f(t)=(P′(β)β−P(β))teβt−P′(β)(βteβt+eβt)=−(P(β)t+P′(β))eβt.
By considering t=0 and t→∞ once again we deduce that f is nonnegative on [0,∞) if and only if P(β)≤0 and P′(β)≤0, which concludes the proof. □

Before turning to the proof of Theorem 5 we will need a couple of auxiliary lemmas. In relation to the formulation of the first result we recall the notation
(μ∗f)jk=∑l∫flk(·−t)μjl(dt)
for a matrix-valued function f=(fjk) and a matrix-valued signed measure μ=(μjk) of suitable dimensions and such that the involved integrals are well defined.

Suppose thatdet(hϕ(z))≠0for allz∈C+and letgϕ:[0,∞)→Rd×dbe characterized by (32). Then it holds thatgϕ(t)=gϕ(t−s)gϕ(s)+∫stgϕ(t−u)(ϕ∗(gϕ1[0,s]))(u)dufor arbitrarys,t∈[0,∞)withs<t.

To lighten notation, and in line with Remark 4, we denote the Laplace transform by L. Specifically, for a given signed measure μ and an integrable function f:[0,∞)→R, set L[μ](z)=L[μ(dt)](z):=∫[0,∞)e−ztμ(dt) for z∈C+ and L[f]:=L[f(t)dt]. From [3, Proposition 5.1] it follows that
gϕ(t)=gϕ(s)+∫st(gϕ∗ϕ)(u)du,t>s≥0,s\ge 0,\]]]>
where (gϕ∗ϕ)jk:=∑l=1d∫[0,∞)gϕ,jl(·−t)ϕlk(dt). Since L[gϕ](z)=hϕ(z)−1 by (32), it is easy to check that L[gϕ] and L[ϕ] commute, and hence
gϕ∗ϕ=ϕ∗gϕ.
For a fixed s∈[0,∞), it thus follows from (34) and (35) that
L[gϕ1(s,∞)](z)=L[1(s,∞)](z)gϕ(s)+L[∫0·(ϕ∗gϕ)(u)1(s,∞)(u)du](z)=1z(L[δs](z)gϕ(s)+L[(ϕ∗(gϕ1[0,s]))1(s,∞)](z)+L[ϕ](z)L[gϕ1(s,∞)](z))
for z∈C with Re(z)>00$]]>. After rearranging terms we obtain
hϕ(z)L[gϕ1(s,∞)](z)=L[δs](z)gϕ(s)+L[(ϕ∗(gϕ1[0,s]))1(s,∞)](z).
By multiplying L[gϕ](z) from the left on both sides of (36) and using uniqueness of the Laplace transform, this shows that
gϕ(t)=(gϕ∗δs)(t)gϕ(s)+(gϕ∗(ϕ∗(gϕ1[0,s])))(t)=gϕ(t−s)gϕ(s)+∫stgϕ(t−u)(ϕ∗(gϕ1[0,s]))(u)du
for t>ss$]]>, and this completes the proof. □

The following lemma, which will be used in the proof of Theorem 5 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.

LetA∈Rd×d. IfAis an M-matrix, thene−Athas nonnegative entries for allt≥0. Conversely, ife−Athas nonnegative entries for allt≥0and the spectrumσ(A)ofAbelongs to{z∈C:Re(z)>0}0\}$]]>, thenAis an invertible M-matrix.

If A is an M-matrix, we have that e−At=e−αteBt for some α∈[0,∞) and B∈Rd×d with nonnegative entries. Since eBt has nonnegative entries, we conclude that the same holds for e−At. If instead the entries of e−At are nonnegative and σ(A)⊆{z∈C:Re(z)>0}0\}$]]>, it follows that the entries of
∫0∞e−xte−Atdt=(Idx+A)−1
are nonnegative for all x∈[0,∞) as well. By [14, Theorem 2] this implies that A is an invertible M-matrix. □

For ε∈(0,∞) define the measure ϕε:=−λδ0+η(·∩(ε,∞)). We start by observing that, when ε is sufficiently small,
det(hϕε(z))≠0for allz∈C+.
To see this suppose, for the sake of contradiction, there exist sequences (εn)n≥1⊆(0,∞) and (zn)n≥1⊆C+ such that det(hϕεn(zn))=0 for all n≥1 and εn→0 as n→∞. Note that the length of any entry of L[ϕε](z) is bounded by the constant
maxj,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 (zn)n≥1 must be bounded, and so it has a subsequence (znk)k≥1 which converges to a point z∗∈C+. However, this would imply that
det(hϕ(z∗))=limk→∞det(hϕεnk(znk))=0,
which contradicts the original assumption that det(hϕ(z))≠0 for all z∈C+. Consequently, the condition (37) is satisfied as long as ε is smaller than a certain threshold, say, ε∗. Thus, for ε≤ε∗ we can define the corresponding function gϕε:[0,∞)→Rd×d through (32) (with ϕ replaced by ϕε). By (34) and (35),
gϕε(t)=Id+∫0t(ϕε∗gϕε)(s)ds=Id−λ∫0tgϕε(s)ds
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 2 to deduce that its entries are nonnegative. Now, given that gϕε(t) has nonnegative entries for t∈[0,kε] for some positive integer k≥1, we can apply Lemma 1 with s=kε to deduce that this remains true when t∈(kε,(k+1)ε]. Consequently, it follows by induction that gϕε(t) has nonnegative entries for all t∈[0,∞).

Next, for ε∈[0,ε∗] define the Fourier transform
Fϕε,jk(y):=(hϕε(iy)−1)jk,y∈R,
of the (j,k)-th entry gϕε,jk of gϕε. By Cramer’s rule it follows that Fϕε,jk(y)=det(hϕε(iy)ˆ)/det(hϕε(iy)), where hϕε(iy)ˆ is the matrix formed by replacing the jth column of hϕε(iy) by the kth canonical basis vector. By the same type of arguments as above (in particular, by relying on the Leibniz formula),
infε≤ε∗|det(hϕε(iy))|≥c1(1∨|y|d)andsupε≤ε∗|det(hϕε(iy)ˆ)|≤c2(1∨|y|d−1)
for suitably chosen constants c1,c2∈(0,∞). Consequently, we conclude that
sup0≤ε≤ε∗|Fϕε,jk(y)|≤c1−1c2(1∧|y|−1),y∈R.
Since Fϕε,jk converges pointwise to Fϕ0,jk (the Fourier transform of the (j,k)-th entry gϕ,jk of gϕ) as ε→0, it follows by the Plancherel theorem and (38) that
∫0∞(gϕε,jk(t)−gϕ,jk(t))2dt=∫R|Fϕε,jk(y)−Fϕ0,jk(y)|2dy→0,ε→0.
From this we establish that gϕεn,jk→gϕ,jk almost everywhere for a suitable sequence (εn)n≥1⊆(0,ε∗] with εn→0 as n→∞. This shows that gϕ,jk(t)≥0 for almost all t∈[0,∞) and, thus, completes the proof. □

Acknowledgement

We wish to thank the referees for their constructive comments which led to an improvement of the manuscript.

ReferencesBall, K.: Completely monotonic rational functions and Hall’s marriage theorem. Barndorff-Nielsen, O.E., Shephard, N.: Non-Gaussian Ornstein–Uhlenbeck-based models and some of their uses in financial economics. Basse-O’Connor, A., Nielsen, M.S., Pedersen, J., Rohde, V.: Multivariate stochastic delay differential equations and CAR representations of CARMA processes. Berg, C.: Stieltjes–Pick–Bernstein–Schoenberg and their connection to complete monotonicity. Positive Definite Functions: From Schoenberg to Space-Time Challenges, 15–45 (2008)Berman, A., Plemmons, R.J.: Bernstein, S.: Sur les fonctions absolument monotones. Brockwell, P.J.: Lévy-driven CARMA processes. In: Brockwell, P.J.: Recent results in the theory and applications of CARMA processes. Brockwell, P.J., Davis, R.A., Yang, Y.: Estimation for non-negative Lévy-driven CARMA processes. Cont, R., Tankov, P.: Gushchin, A.A., Küchler, U.: On stationary solutions of delay differential equations driven by a Lévy process. Marquardt, T., Stelzer, R.: Multivariate CARMA processes. Mohammed, S.E.A., Scheutzow, M.K.R.: Lyapunov exponents and stationary solutions for affine stochastic delay equations. Plemmons, R.J.: M-matrix characterizations. I. Nonsingular M-matrices. Rajput, B.S., Rosiński, J.: Spectral representations of infinitely divisible processes. Roman, S.: The formula of Faà di Bruno. Sato, K.: Sato, K., Yamazato, M.: Stationary processes of Ornstein-Uhlenbeck type. In: Sato, K., Watanabe, T., Yamazato, M.: Recurrence conditions for multidimensional processes of Ornstein-Uhlenbeck type. Schilling, R.L., Song, R., Vondraček, Z.: Tsai, H., Chan, K.S.: A note on non-negative continuous time processes.