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) ${({L_{t}})_{t\in \mathbb{R}}}$, that is, the unique stationary solution to for some $\lambda \in (0,\infty )$. Since the solution of (1) is explicitly given by a convolution between a nonnegative kernel $t\mapsto {e^{-\lambda t}}$ and a nonnegative random measure $\operatorname{d}L$, it is automatically nonnegative. However, due to the simplicity of the Ornstein–Uhlenbeck process, its autocorrelation function is $\text{Corr}({X_{0}},{X_{h}})={e^{-\lambda |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\in \mathbb{R}}}$ can be characterized as a moving average of the form where the Fourier transform of $g:[0,\infty )\to \mathbb{R}$ is rational. Consequently, a CARMA process is made up of a Lévy process ${({L_{t}})_{t\in \mathbb{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\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{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 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,\infty )$. This means that ${({X_{t}})_{t\in \mathbb{R}}}$ is not necessarily nonnegative although ${({L_{t}})_{t\in \mathbb{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 ${({X_{t}})_{t\in \mathbb{R}}}$ of (4) to be nonnegative, it is sufficient that ϕ be a nonnegative measure when restricted to $(0,\infty )$, that is, 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\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{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,\infty )$, 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 $\boldsymbol{\phi }$, one needs to require that $\boldsymbol{\lambda }:=-\boldsymbol{\phi }(\{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 ${({L_{t}})_{t\ge 0}}$, ${L_{0}}\equiv 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 ${({L_{t}})_{t\ge 0}}$ (that is, all its finite-dimensional marginal distributions) 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 $\gamma \in [0,\infty )$ and ν is a σ-finite measure on $(0,\infty )$ with ${\textstyle\int _{0}^{\infty }}(1\wedge x)\hspace{0.1667em}\nu (\operatorname{d}x)<\infty $ (in particular, ν is a Lévy measure). For any subordinator ${({L_{t}})_{t\ge 0}}$ the induced pair $(\gamma ,\nu )$ 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 ${({L_{t}})_{t\in \mathbb{R}}}$ from two independent one-sided subordinators ${({L_{t}^{1}})_{t\ge 0}}$ and ${({L_{t}^{2}})_{t\ge 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 nonnegative i.i.d. random variables. We may also refer to multivariate subordinators, which are simply Lévy processes with values in ${\mathbb{R}^{d}}$, $d\ge 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,\infty )$, is a signed measure if $\mu ={\mu _{+}}-{\mu _{-}}$ for two mutually singular finite Borel measures ${\mu _{+}}$ and ${\mu _{-}}$ on $[0,\infty )$. Note that the variation $|\mu |:={\mu _{+}}+{\mu _{-}}$ of μ is a measure on $[0,\infty )$ and that integration with respect to μ can be defined in an obvious manner for any function $f:[0,\infty )\to \mathbb{R}$ which is integrable with respect to $|\mu |$.

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 $\boldsymbol{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.

Let ${({L_{t}})_{t\in \mathbb{R}}}$ be a two-sided subordinator with a nonzero Lévy measure ν satisfying ${\textstyle\int _{1}^{\infty }}x\hspace{0.1667em}\nu (\operatorname{d}x)<\infty $ or, equivalently, $\mathbb{E}[{L_{1}}]<\infty $. Moreover, let ϕ be a signed measure on $[0,\infty )$ with A stochastic process ${({X_{t}})_{t\in \mathbb{R}}}$ is said to be a solution of the corresponding SDDE if it is stationary, has finite first moments (that is, $\mathbb{E}[|{X_{0}}|]<\infty $), and By (7), we mean that the equality holds almost surely for each fixed pair $(s,t)\in {\mathbb{R}^{2}}$ with $s<t$. We will often write the equation in differential form as in (4).

Set ${\mathbb{C}_{+}}:=\{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)\ge 0\}$. Properties such as existence and uniqueness of solutions of (7) are closely related to the zeroes of the function ${h_{\phi }}:{\mathbb{C}_{+}}\to \mathbb{C}$ given by Note that ${h_{\phi }}$ is analytic on $\{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\text{Re}(z)>0\}$ and real-valued on $\{z\in {\mathbb{C}_{+}}\hspace{0.1667em}:\hspace{0.1667em}\text{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]:

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.

Turning to the question of whether a given solution ${({X_{t}})_{t\in \mathbb{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,\infty )\to \mathbb{R}$ is called completely monotone if it is infinitely differentiable and 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\circ 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].

While Theorem 2 tells that it is sufficient to show nonnegativity of ${g_{\phi }}$, 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 $\lambda >0$ and $\xi =0$ is the Ornstein–Uhlenbeck process, so ${g_{\phi }}(t)={e^{-\lambda t}}$ (by (10) as ${h_{\phi }}(z)=z+\lambda $) and, thus, proves that the solution is nonnegative. In contrast, as soon as $\xi \ne 0$, ${g_{\phi }}$ cannot be explicitly determined since its structure depends on the infinitely many solutions of the equation $z+\lambda -\xi {e^{-\tau z}}=0$ for $z\in \mathbb{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.

In light of Theorem 3, when one is trying to model nonnegative processes, it is natural to write (7) as and then search for nonnegative measures η on $(0,\infty )$ satisfying $\eta ((0,\infty ))<\lambda $ (by the last part of Remark 2 this inequality must hold if ${h_{\phi }}(z)\ne 0$ for all $z\in {\mathbb{C}_{+}}$). Here we use the convention ${\textstyle\int _{0}^{\infty }}:={\textstyle\int _{(0,\infty )}}$.

In Figure 1 we simulate the stationary solution of (17) when $\eta =\xi {\delta _{\tau }}$, ${\delta _{\tau }}$ being the Dirac measure at τ, with the specific values $\lambda =\tau =1$ and $\xi =0.2$. For comparison we rerun the simulation with $\xi =-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 $\xi =0.2$, but when $\xi =-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 when $\lambda =\tau =1$, and hence $1/{h_{\phi }}$ is not completely monotone on $[0,\infty )$ if $\xi =-0.8$.

Let $P(z)={z^{p}}+{a_{1}}{z^{p-1}}+\cdots +{a_{p}}$ and $Q(z)={b_{0}}+{b_{1}}z+\cdots +{b_{q-1}}{z^{q-1}}+{z^{q}}$ be two real monic polynomials with $p>q$, and assume that P has no zeroes on ${\mathbb{C}_{+}}$ (causality). Define the companion matrix and let $\boldsymbol{b}={({b_{0}},{b_{1}},\dots ,{b_{q-1}},1,0,\dots ,0)^{\top }}$ and ${\boldsymbol{e}_{p}}={(0,0,\dots ,0,1)^{\top }}$ (both being elements of ${\mathbb{R}^{p}}$). The causal CARMA($p,q$) process ${({Y_{t}})_{t\in \mathbb{R}}}$ driven by the subordinator ${({L_{t}})_{t\in \mathbb{R}}}$, and associated to the autoregressive polynomial P and moving average polynomial Q, is defined by Alternatively, the kernel $g(t)={\boldsymbol{b}^{\top }}{e^{\boldsymbol{A}t}}{\boldsymbol{e}_{p}}$ can be characterized in the frequency domain by the relation Due to the well-known form ${\boldsymbol{X}_{t}}={\textstyle\int _{-\infty }^{t}}{e^{\boldsymbol{A}(t-s)}}{\boldsymbol{e}_{p}}\hspace{0.1667em}\operatorname{d}{L_{s}}$ of the stationary Ornstein–Uhlenbeck process with drift parameter $\boldsymbol{A}$ and driven by ${({\boldsymbol{e}_{p}}{L_{t}})_{t\in \mathbb{R}}}$ (see, for example, [19, 18]), (19) shows immediately that ${({Y_{t}})_{t\in \mathbb{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\in \mathbb{R}}}$ should be the solution of the formal differential equation (For instance, by heuristically computing the Fourier transform of ${({Y_{t}})_{t\in \mathbb{R}}}$ from (21), one can deduce that ${Y_{t}}$ 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+\lambda $ for $\lambda >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\in \mathbb{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 ${\alpha _{1}},\dots ,{\alpha _{p}}$ and ${\beta _{1}},\dots ,{\beta _{q}}$ being the zeroes of P and Q, respectively, Ball [1] showed that g is nonnegative if where the latter condition assumes that the zeroes are ordered so that ${\alpha _{1}}\ge \cdots \ge {\alpha _{p}}$ and ${\beta _{1}}\ge \cdots \ge {\beta _{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 ${\mathbb{C}_{+}}$ (invertibility), we saw in Example 1 that the corresponding causal and invertible CARMA process ${({Y_{t}})_{t\in \mathbb{R}}}$ can be identified as the unique solution of an SDDE of the form (17) with $\eta (\operatorname{d}t)=f(t)\hspace{0.1667em}\operatorname{d}t$ and $f:[0,\infty )\to \mathbb{R}$ being characterized by Here $R(z):=(z+\lambda )Q(z)-P(z)$ and $\lambda \in \mathbb{R}$ is uniquely determined by the condition $\deg (R)<q$; it is explicitly given by Note that, in case the zeroes ${\beta _{1}},\dots ,{\beta _{q}}$ of Q are distinct, Cauchy’s residue theorem implies that ${Q^{\prime }}$ being the derivative of Q (see also [8, Remark 5]). In view of Theorem 3, it follows that ${({Y_{t}})_{t\in \mathbb{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,\infty )$ for ${({Y_{t}})_{t\in \mathbb{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.

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\ge 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\ge 3$, but this would involve more complicated expressions.

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\ge 2$ we let ${\boldsymbol{L}_{t}}={({L_{t}^{1}},\dots ,{L_{t}^{d}})^{\top }}$, $t\in \mathbb{R}$, be a d-dimensional subordinator with $\mathbb{E}[\| {\boldsymbol{L}_{1}}\| ]<\infty $ and $\boldsymbol{\phi }=({\phi _{jk}})$ a $d\times d$ matrix-valued signed measure on $[0,\infty )$ with A stochastic process ${\boldsymbol{X}_{t}}={({X_{t}^{1}},\dots ,{X_{t}^{d}})^{\top }}$, $t\in \mathbb{R}$, is a solution of the corresponding MSDDE if it is stationary, has finite first moments (that is, $\mathbb{E}[\| {\boldsymbol{X}_{0}}\| ]<\infty $), and satisfies for $j=1,\dots ,d$. In line with (3) we use the decomposition $\boldsymbol{\phi }=-\boldsymbol{\lambda }{\delta _{0}}+\boldsymbol{\eta }$, where $\boldsymbol{\lambda }:=-\boldsymbol{\phi }(\{0\})$ and $\boldsymbol{\eta }:=\boldsymbol{\phi }(\hspace{0.2222em}\cdot \hspace{0.2222em}\cap (0,\infty ))$. If we define the convolution $(\mu \ast f)(t)=(f\ast \mu )(t):=\textstyle\int f(t-s)\hspace{0.1667em}\mu (\operatorname{d}s)$ 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 Among stationary processes which can be identified as solutions of (28) are the multivariate Ornstein–Uhlenbeck process ($\boldsymbol{\eta }\equiv \mathbf{0}$) and, more generally, multivariate CARMA processes [12] ($\boldsymbol{\eta }(\operatorname{d}t)=\boldsymbol{f}(t)\hspace{0.1667em}\operatorname{d}t$ with $\boldsymbol{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 ${\boldsymbol{I}_{d}}$ being the $d\times d$ identity matrix. In particular, it was shown in [3] that if $\det ({\boldsymbol{h}_{\phi }}(z))\ne 0$ for all $z\in {\mathbb{C}_{+}}$, the unique solution ${({\boldsymbol{X}_{t}})_{t\in \mathbb{R}}}$ is given by where ${\boldsymbol{g}_{\phi }}:[0,\infty )\to {\mathbb{R}^{d\times d}}$ is characterized by In Theorem 5, which is a generalization of Theorem 3, we give sufficient conditions for the solution ${({\boldsymbol{X}_{t}})_{t\in \mathbb{R}}}$ to have nonnegative entries. Before formulating this result we recall the definition of the so-called M-matrices: $\boldsymbol{A}\in {\mathbb{R}^{d\times d}}$ is called an M-matrix if it can be expressed as $\boldsymbol{A}=\alpha {\boldsymbol{I}_{d}}-\boldsymbol{B}$ for some $\alpha \in [0,\infty )$ and $\boldsymbol{B}\in {\mathbb{R}^{d\times d}}$ with nonnegative entries and spectral radius $\rho (\boldsymbol{B})\le \alpha $. The main point in introducing these is their relation to entrywise nonnegativity of the matrix exponential ${e^{-\boldsymbol{A}t}}$ for all $t\ge 0$ (see Lemma 2). For further details and other characterizations of M-matrices, we refer to [5, 14] and references therein.

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.

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 for a matrix-valued function $\boldsymbol{f}=({f_{jk}})$ and a matrix-valued signed measure $\boldsymbol{\mu }=({\mu _{jk}})$ of suitable dimensions and such that the involved integrals are well defined.

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.