## 1 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) ${({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 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 where

##### (1)

\[ \operatorname{d}{X_{t}}=-\lambda {X_{t}}\operatorname{d}t+\operatorname{d}{L_{t}},\hspace{2em}t\in \mathbb{R},\]##### (2)

\[ {X_{t}}={\int _{-\infty }^{t}}{e^{-\lambda (t-s)}}\hspace{0.1667em}\operatorname{d}{L_{s}},\hspace{2em}t\in \mathbb{R},\]*λ*; 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##### (3)

\[ {X_{t}}={\int _{-\infty }^{t}}g(t-s)\hspace{0.1667em}\operatorname{d}{L_{s}},\hspace{2em}t\in \mathbb{R},\]*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##### (4)

\[ \operatorname{d}{X_{t}}={\int _{[0,\infty )}}{X_{t-s}}\hspace{0.1667em}\phi (\operatorname{d}s)\hspace{0.1667em}\operatorname{d}t+\operatorname{d}{L_{t}},\hspace{2em}t\in \mathbb{R},\]*ϕ*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 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

*ϕ*be a nonnegative measure when restricted to $(0,\infty )$, that is,##### (5)

\[ \phi (B\cap (0,\infty ))\ge 0\hspace{2em}\text{for all Borel measurable sets}\hspace{2.5pt}B\text{.}\]*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,

\[ \log \mathbb{E}[{e^{i\theta {L_{1}}}}]=i\theta \gamma +{\int _{0}^{\infty }}({e^{i\theta x}}-1)\hspace{0.1667em}\nu (\operatorname{d}x),\hspace{2em}\theta \in \mathbb{R},\]

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

\[ {L_{t}}={L_{t}^{1}}{\mathbb{1}_{[0,\infty )}}(t)-{L_{(-t)-}^{2}}{\mathbb{1}_{(-\infty ,0)}}(t),\hspace{2em}t\in \mathbb{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 ${\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.

## 2 Stochastic delay differential equations and nonnegative solutions

Let ${({L_{t}})_{t\in \mathbb{R}}}$ be a two-sided subordinator with a nonzero Lévy measure By (7), we mean that the equality holds almost surely for each

*ν*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##### (7)

\[ {X_{t}}-{X_{s}}={\int _{s}^{t}}{\int _{[0,\infty )}}{X_{u-v}}\hspace{0.1667em}\phi (\operatorname{d}v)\hspace{0.1667em}\operatorname{d}u+{L_{t}}-{L_{s}},\hspace{2em}s<t.\]*fixed*pair $(s,t)\in {\mathbb{R}^{2}}$ with $s<t$. We will often write the equation in differential form as in (4).##### Remark 1.

The condition that a stationary solution ${({X_{t}})_{t\in \mathbb{R}}}$ must satisfy $\mathbb{E}[|{X_{0}}|]<\infty $ 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

\[ {\int _{s}^{t}}{\int _{[0,\infty )}}|{X_{u-v}}|\hspace{0.1667em}\phi (\operatorname{d}v)\hspace{0.1667em}\operatorname{d}u<\infty \]

almost surely since its mean is finite. Note also that, if (7) has any such solution ${({X_{t}})_{t\in \mathbb{R}}}$, then $\mathbb{E}[{L_{1}}]=-\phi ([0,\infty ))\mathbb{E}[{X_{0}}]$. Consequently, our assumption that $\mathbb{E}[{L_{1}}]<\infty $ is innocent and only imposed to avoid trivial cases. If *ϕ*has bounded support, (7) makes perfect sense without the condition $\mathbb{E}[|{X_{0}}|]<\infty $, and existence and uniqueness of solutions were established in [11, 13] under milder conditions on ${({L_{t}})_{t\in \mathbb{R}}}$. While it should be doable to study (7) under milder conditions when*ϕ*is unbounded as well, we leave this for future research.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

##### (8)

\[ {h_{\phi }}(z):=z-{\int _{[0,\infty )}}{e^{-zt}}\hspace{0.1667em}\phi (\operatorname{d}t),\hspace{2em}z\in {\mathbb{C}_{+}}.\]*ϕ*with unbounded support, we will rely on the following existence and uniqueness result of [3]:##### Theorem 1 (Basse-O’Connor et al. [3]).

*Let*${h_{\phi }}$

*be given as in*(8)

*and assume that*${h_{\phi }}(z)\ne 0$

*for all*$z\in {\mathbb{C}_{+}}$

*. Then the unique solution of*(7)

*is given by*

##### (9)

\[ {X_{t}}={\int _{-\infty }^{t}}{g_{\phi }}(t-u)\hspace{0.1667em}\operatorname{d}{L_{u}},\hspace{2em}t\in \mathbb{R},\]*where*${g_{\phi }}:[0,\infty )\to \mathbb{R}$

*is characterized by*

##### Remark 2.

The assumption that ${h_{\phi }}(z)\ne 0$ for all $z\in {\mathbb{C}_{+}}$ as well as the form of the solution (9) are, in fact, rather natural. Indeed, heuristically speaking the equation (7) reads ${({\mathcal{A}_{\phi }}X)_{t}}={(DL)_{t}}$ with ${\mathcal{A}_{\phi }}$ being the linear operator given by $({\mathcal{A}_{\phi }}f)(t)=(Df)(t)+{\textstyle\int _{[0,\infty )}}f(t-s)\hspace{0.1667em}\phi (\operatorname{d}s)$ and $D=\frac{\operatorname{d}}{\operatorname{d}t}$ the derivative with respect to

*t*. Since \[ {\int _{\mathbb{R}}}{e^{-ity}}({\mathcal{A}_{\phi }}f)(t)\hspace{0.1667em}\operatorname{d}t={h_{\phi }}(iy){\int _{\mathbb{R}}}{e^{-ity}}f(t)\hspace{0.1667em}\operatorname{d}t,\hspace{2em}y\in \mathbb{R},\]

the assumption on ${h_{\phi }}$ can be interpreted as an invertibility assumption on ${\mathcal{A}_{\phi }}$, and its inverse ${\mathcal{A}_{\phi }^{-1}}$ takes the form $({\mathcal{A}_{\phi }^{-1}}f)(t)={\textstyle\int _{-\infty }^{t}}{g_{\phi }}(t-s)f(s)\hspace{0.1667em}\operatorname{d}s$. While it might sometimes be challenging to show that ${h_{\phi }}(z)\ne 0$ when $z\in {\mathbb{C}_{+}}$, as required by Theorem 1, a necessary and easy-to-check condition is that $\phi ([0,\infty ))<0$. Indeed, this follows from the facts that ${h_{\phi }}(0)=-\phi ([0,\infty ))$, ${h_{\phi }}$ is continuous, and ${h_{\phi }}(x)\to \infty $ as $x\to \infty $ (where $x\in [0,\infty )$).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.

##### Example 1 (CARMA processes).

Let $\lambda \in \mathbb{R}$ and, for a given $p\in \mathbb{N}$, suppose $\boldsymbol{b}\in {\mathbb{R}^{p}}$ and $\boldsymbol{A}\in {\mathbb{R}^{p\times p}}$ with a spectrum $\sigma (\boldsymbol{A})$ contained in $\{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)>0\}$. Consider the delay measure where ${\delta _{0}}$ is the Dirac measure at 0 and ${\boldsymbol{e}_{1}}$ is the first canonical basis vector of ${\mathbb{R}^{p}}$. Note that ${\boldsymbol{e}_{1}}$ is merely used as a normalization: the effect of replacing ${\boldsymbol{e}_{1}}$ by an arbitrary vector $\boldsymbol{c}\in {\mathbb{R}^{p}}$ can be incorporated in $\boldsymbol{b}$ and $\boldsymbol{A}$. By the assumption on $\sigma (\boldsymbol{A})$, all the entries of ${e^{-\boldsymbol{A}t}}$ are exponentially decaying as $t\to \infty $ and, thus, $|\phi |$ is a finite measure with moments of any order (in particular, (6) is satisfied). The function ${h_{\phi }}$ takes the form ${h_{\phi }}(z)=z+\lambda -{\boldsymbol{b}^{\top }}{(\boldsymbol{A}+z{\boldsymbol{I}_{p}})^{-1}}{\boldsymbol{e}_{1}}$, where ${\boldsymbol{I}_{p}}$ is the $p\times p$ identity matrix. By the fraction decomposition it follows that ${\boldsymbol{b}^{\top }}{(\boldsymbol{A}+z{\boldsymbol{I}_{p}})^{-1}}{\boldsymbol{e}_{1}}=R(z)/Q(z)$ for $z\in {\mathbb{C}_{+}}$, where $(Q,R)$ is the unique pair of real polynomials $Q,R:\mathbb{C}\to \mathbb{C}$ such that (i) In other words, ${({X_{t}})_{t\in \mathbb{R}}}$ is a causal and invertible CARMA process with autoregressive polynomial

##### (11)

\[ \phi (\operatorname{d}t)=-\lambda \hspace{0.1667em}{\delta _{0}}(\operatorname{d}t)+{\boldsymbol{b}^{\top }}{e^{-\boldsymbol{A}t}}{\boldsymbol{e}_{1}}\hspace{0.1667em}\operatorname{d}t,\]*Q*is monic and has no zeroes on ${\mathbb{C}_{+}}$, (ii) $\text{deg}(Q)>\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+\lambda )Q(z)-R(z)\ne 0\hspace{2em}\text{for all}\hspace{2.5pt}z\in {\mathbb{C}_{+}}\text{,}\]

there exists a unique stationary solution ${({X_{t}})_{t\in \mathbb{R}}}$ to (7), and it is given by (9) with ${g_{\phi }}$ satisfying ##### (12)

\[ {\int _{0}^{\infty }}{e^{-ity}}{g_{\phi }}(t)\hspace{0.1667em}\operatorname{d}t=\frac{Q(iy)}{P(iy)},\hspace{2em}y\in \mathbb{R}.\]*P*and moving average polynomial*Q*(see Section 3 or [8, Remark 4]). Conversely, given a moving average \[ {X_{t}}={\int _{-\infty }^{t}}g(t-s)\hspace{0.1667em}\operatorname{d}{L_{s}},\hspace{2em}t\in \mathbb{R},\]

with $g:[0,\infty )\to \mathbb{R}$ characterized by (12) for some polynomials *P*and*Q*having no zeroes in ${\mathbb{C}_{+}}$, and which satisfy $\text{deg}(P)=\text{deg}(Q)+1$, one can choose a unique constant $\lambda \in \mathbb{R}$ such that the polynomial $R(z):=(z+\lambda )Q(z)-P(z)$ meets $\text{deg}(R)<\text{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,\infty )\to \mathbb{R}$ with a rational Laplace transform $R/Q$ can always be represented as $f(t)={\boldsymbol{b}^{\top }}{e^{-\boldsymbol{A}t}}{\boldsymbol{e}_{1}}$ for suitable $\boldsymbol{b}\in {\mathbb{R}^{\deg (Q)}}$ and $\boldsymbol{A}\in {\mathbb{R}^{\deg (Q)\times \deg (Q)}}$ with $\sigma (\boldsymbol{A})\subseteq \{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)>0\}$. For more on the relation between solutions of SDDEs and CARMA processes, see [3, Section 4.3].##### Example 2 (Discrete delay).

Let $\lambda ,\xi \in \mathbb{R}$ and $\tau \in (0,\infty )$, and assume that $|\xi |\le {\tau ^{-1}}$. Consider the following SDDE written in differential form: To show existence of a unique stationary solution using Theorem 1, we must argue that $z+\lambda -\xi {e^{-z\tau }}\ne 0$ whenever $z\in {\mathbb{C}_{+}}$ or, equivalently, that the two equations cannot hold simultaneously if $x\ge 0$ and $y\in \mathbb{R}$. Since the only real solution of an equation of the form $u=\sin (\alpha u)$ is $u=0$ when $|\alpha |\le 1$, the second equation in (14) implies that $y=0$. If this is the case we have that since $u\mapsto u-\xi {e^{-\tau u}}$ is increasing on $[0,\infty )$. In view of the last part of Remark 2, we conclude that ${h_{\phi }}(z)\ne 0$ for all $z\in {\mathbb{C}_{+}}$ if and only if $\xi <\lambda $.

##### (13)

\[ \operatorname{d}{X_{t}}=(-\lambda {X_{t}}+\xi {X_{t-\tau }})\hspace{0.1667em}\operatorname{d}t+\operatorname{d}{L_{t}},\hspace{2em}t\in \mathbb{R}.\]##### (14)

\[ x+\lambda -\xi {e^{-\tau x}}\cos (\tau y)=0\hspace{2em}\text{and}\hspace{2em}y+\xi {e^{-\tau x}}\sin (\tau y)=0\]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

##### (15)

\[ {(-1)^{n}}\frac{{\operatorname{d}^{\hspace{0.1667em}n}}}{\operatorname{d}{t^{n}}}f(t)\ge 0\hspace{2em}\text{for all}\hspace{2.5pt}n\in {\mathbb{N}_{0}}\hspace{2.5pt}\text{and}\hspace{2.5pt}t>0\text{.}\]*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].##### Theorem 2.

##### Remark 3.

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

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.##### Theorem 3.

*Let*${h_{\phi }}$

*be defined as in*(8)

*and assume that*${h_{\phi }}(z)\ne 0$

*for all*$z\in {\mathbb{C}_{+}}$

*. Suppose further that ϕ is a nonnegative measure when restricted to*$(0,\infty )$

*, that is,*

##### (16)

\[ \phi (B\cap (0,\infty ))\ge 0\hspace{2em}\textit{for all Borel measurable sets}\hspace{2.5pt}B\textit{.}\]*Then*(i)

*–*(iv)

*of Theorem*

*2*

*are satisfied.*

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

##### (17)

\[ \operatorname{d}{X_{t}}=-\lambda {X_{t}}\hspace{0.1667em}\operatorname{d}t+{\int _{0}^{\infty }}{X_{t-s}}\hspace{0.1667em}\eta (\operatorname{d}s)\hspace{0.1667em}\operatorname{d}t+\operatorname{d}{L_{t}},\hspace{2em}t\in \mathbb{R},\]*η*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 )}}$.##### Remark 4.

As will, for instance, appear from Example 3, the condition (16) is not necessary for ${({X_{t}})_{t\in \mathbb{R}}}$ to be nonnegative. According to Theorem 2 it is necessary and sufficient that (15) is satisfied for $1/{h_{\phi }}$. The case $n=0$ is always true when ${h_{\phi }}(z)\ne 0$ for all $z\in {\mathbb{C}_{+}}$, and by relying on Faà di Bruno’s formula (see, e.g., [16, Theorem 2]) it can be checked that when $n\ge 1$, where ${\mathcal{L}_{\eta }^{(j)}}$ denotes the

##### (18)

\[ \begin{aligned}{}& {(-1)^{n}}\frac{{\operatorname{d}^{n}}}{\operatorname{d}{x^{n}}}\frac{1}{{h_{\phi }}(x)}=\\ {} & \sum \limits_{\alpha \in {\mathbb{N}_{0}^{n}}:{\textstyle\textstyle\sum _{j=1}^{n}}j{\alpha _{j}}=n}\frac{n!|\alpha |!{(1-{\mathcal{L}_{\eta }^{(1)}}(x))^{{\alpha _{1}}}}}{{\alpha _{1}}!{(1!)^{{\alpha _{1}}}}\cdots {\alpha _{n}}!{(n!)^{{\alpha _{n}}}}{h_{\phi }}{(x)^{|\alpha |+1}}}{\prod \limits_{j=2}^{n}}{({(-1)^{j}}{\mathcal{L}_{\eta }^{(j)}}(x))^{{\alpha _{j}}}}\end{aligned}\]*j*th derivative of the Laplace transform ${\mathcal{L}_{\eta }}(x):={\textstyle\int _{0}^{\infty }}{e^{-xt}}\hspace{0.1667em}\eta (\operatorname{d}t)$ of $\eta :=\phi (\hspace{0.2222em}\cdot \hspace{0.2222em}\cap (0,\infty ))$. Moreover, ${\alpha _{j}}$ refers to the*j*th entry of $\alpha \in {\mathbb{N}_{0}^{n}}$ and $|\alpha |={\textstyle\sum _{j=1}^{n}}{\alpha _{j}}$. From this expression it is easy to see that (15) is satisfied if*η*is nonnegative, since then ${\mathcal{L}_{\eta }}$ 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 ${({X_{t}})_{t\in \mathbb{R}}}$ is required to be nonnegative; for instance, the restriction for $n=1$ implies in particular that 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 $\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 \[ {h_{\phi }}{(x)^{3}}\frac{{\operatorname{d}^{2}}}{\operatorname{d}{x^{2}}}\frac{1}{{h_{\phi }}(x)}\bigg|{_{x=0}}={\xi ^{2}}+5\xi +2\]

when $\lambda =\tau =1$, and hence $1/{h_{\phi }}$ is not completely monotone on $[0,\infty )$ if $\xi =-0.8$.##### Fig. 1.

Simulations of ${X_{1}},\dots ,{X_{200}}$ from the model (17) with $\eta =\xi {\delta _{1}}$ when $\xi =0.2$ (left) and $\xi =-0.8$ (right). In both cases, $\lambda =1$

## 3 Nonnegative CARMA processes

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 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:

*P*has no zeroes on ${\mathbb{C}_{+}}$ (causality). Define the companion matrix \[ \boldsymbol{A}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}0& 1& 0& \cdots & 0\\ {} 0& 0& 1& \cdots & 0\\ {} \vdots & \vdots & \ddots & \ddots & \vdots \\ {} 0& 0& \cdots & 0& 1\\ {} -{a_{p}}& -{a_{p-1}}& \cdots & -{a_{2}}& -{a_{1}}\end{array}\right),\]

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##### (19)

\[ {Y_{t}}={\int _{-\infty }^{t}}{\boldsymbol{b}^{\top }}{e^{\boldsymbol{A}(t-s)}}{\boldsymbol{e}_{p}}\hspace{0.1667em}\operatorname{d}{L_{s}},\hspace{2em}t\in \mathbb{R}.\]##### (20)

\[ {\int _{0}^{\infty }}{e^{-ity}}g(t)\hspace{0.1667em}\operatorname{d}t=\frac{Q(iy)}{P(iy)},\hspace{2em}y\in \mathbb{R}.\] \[\begin{aligned}{}\operatorname{d}{\boldsymbol{X}_{t}}& =\boldsymbol{A}{\boldsymbol{X}_{t}}\hspace{0.1667em}\operatorname{d}t+{\boldsymbol{e}_{p}}\hspace{0.1667em}\operatorname{d}{L_{t}},\\ {} {Y_{t}}& ={\boldsymbol{b}^{\top }}{\boldsymbol{X}_{t}}.\end{aligned}\]

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 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( 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 ${Q^{\prime }}$ being the derivative of

*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##### (22)

\[ \begin{aligned}{}{\alpha _{1}},\dots ,{\alpha _{p}},{\beta _{1}},\dots ,{\beta _{q}}& \in \{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)<0,\hspace{0.1667em}\operatorname{Im}(z)=0\}\\ {} \text{and}\hspace{2em}{\sum \limits_{j=1}^{k}}{\alpha _{j}}& \ge {\sum \limits_{j=1}^{k}}{\beta _{j}}\hspace{2em}\text{for}\hspace{2.5pt}k=1,\dots ,q,\end{aligned}\]*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##### (23)

\[ {\int _{0}^{\infty }}{e^{-ity}}f(t)\hspace{0.1667em}\operatorname{d}t=\frac{R(iy)}{Q(iy)},\hspace{2em}y\in \mathbb{R}.\] \[ \lambda ={a_{1}}-{b_{p-2}}=-{\alpha _{p}}+{\sum \limits_{j=1}^{p-1}}({\beta _{j}}-{\alpha _{j}}).\]

Note that, in case the zeroes ${\beta _{1}},\dots ,{\beta _{q}}$ of *Q*are distinct, Cauchy’s residue theorem implies that##### (24)

\[ f(t)=-{\sum \limits_{j=1}^{q}}\frac{P({\beta _{j}})}{{Q^{\prime }}({\beta _{j}})}{e^{{\beta _{j}}t}},\hspace{2em}t\ge 0,\]*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.##### Theorem 4.

*Let*${({Y_{t}})_{t\in \mathbb{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\in \mathbb{R}}}$

*is nonnegative if*$R/Q$

*is completely monotone on*$[0,\infty )$

*or, equivalently, if the function f characterized by*(23)

*is nonnegative.*

##### Remark 5.

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 Then the CARMA($p,q$) process associated with the pair $({P_{1}}\times {P_{2}},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).

##### Example 3.

Suppose that ${({Y_{t}})_{t\in \mathbb{R}}}$ is a causal and invertible CARMA($2,1$) process, so that its polynomials On the other hand, from (24) we have $f(t)=-P(\gamma ){e^{\gamma t}}$, and so one would need to require that $P(\gamma )=(\gamma -\alpha )(\gamma -\beta )\le 0$ in order to apply Theorem 4. Imposing this condition is equivalent to assuming that While the first restriction in (26) is necessary, the second is not. In general, we see that (25) and (26) differ the most when $|\alpha -\beta |$ is small. In Figure 2 we have marked the feasible area for

*P*and*Q*are representable as $P(z)=(z-\alpha )(z-\beta )$ and $Q(z)=z-\gamma $ for suitable $\alpha ,\beta \in \{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)<0\}$ and $\gamma \in (-\infty ,0)$. It follows from (22) that for ${({Y_{t}})_{t\in \mathbb{R}}}$ to be nonnegative it is necessary and sufficient that##### (25)

\[ \operatorname{Im}(\alpha )=\operatorname{Im}(\beta )=0\hspace{2em}\text{and}\hspace{2em}\gamma \le \max \{\alpha ,\beta \}.\]##### (26)

\[ \operatorname{Im}(\alpha )=\operatorname{Im}(\beta )=0\hspace{2em}\text{and}\hspace{2em}\min \{\alpha ,\beta \}\le \gamma \le \max \{\alpha ,\beta \}.\]*γ*implied by (25) and (26), respectively, for a given polynomial*P*with zeroes $\alpha ,\beta \in (-\infty ,0)$.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.

##### Corollary 1.

*Let*${({Y_{t}})_{t\in \mathbb{R}}}$

*be a causal and invertible CARMA(*$3,2$

*) process with polynomials P and Q, and let*${\beta _{1}}$

*and*${\beta _{2}}$

*,*$\operatorname{Re}({\beta _{1}})\ge \operatorname{Re}({\beta _{2}})$

*, denote the zeroes of Q. Then the function f characterized by*(23)

*is nonnegative if and only if*$\operatorname{Im}({\beta _{1}})=\operatorname{Im}({\beta _{2}})=0$

*and one of the following conditions holds:*

*In particular, if either*(i)

*or*(ii)

*is satisfied, then*${({Y_{t}})_{t\in \mathbb{R}}}$

*is nonnegative.*

##### Example 4.

Let ${({Y_{t}})_{t\in \mathbb{R}}}$ be the causal and invertible CARMA($3,2$) associated to real polynomials (Here, as in (22), the zeroes are ordered such that ${\alpha _{1}}\ge {\alpha _{2}}\ge {\alpha _{3}}$ and ${\beta _{1}}\ge {\beta _{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 ${\alpha _{1}}$, ${\alpha _{2}}$, and ${\alpha _{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 ${\alpha _{1}}>{\alpha _{2}}={\alpha _{3}}$ and ${\beta _{1}}={\beta _{2}}$. If ${\beta _{1}}$ is larger than ${\alpha _{2}}$, but smaller than the local extremum ${x^{\ast }}$ of

*P*and*Q*representable as \[ P(z)=(z-{\alpha _{1}})(z-{\alpha _{2}})(z-{\alpha _{3}})\hspace{1em}\text{and}\hspace{1em}Q(z)=(z-{\beta _{1}})(z-{\beta _{2}}),\hspace{2em}z\in \mathbb{C},\]

for some ${\alpha _{1}},{\alpha _{2}},{\alpha _{3}}\in \{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)<0\}$ and ${\beta _{1}},{\beta _{2}}\in (-\infty ,0)$. If one seeks to apply the result of [1, 21] stated in (22), it must be required that ##### (27)

\[ \operatorname{Im}({\alpha _{j}})=0,\hspace{2em}{\beta _{1}}\le {\alpha _{1}},\hspace{2em}\text{and}\hspace{2em}{\beta _{1}}+{\beta _{2}}\le {\alpha _{1}}+{\alpha _{2}}.\]*P*located in the interval $({\alpha _{2}},{\alpha _{1}})$, it follows that both $P({\beta _{1}})\le 0$ and ${P^{\prime }}({\beta _{1}})\le 0$. It is easy to check that ${x^{\ast }}=(2{\alpha _{1}}+{\alpha _{2}})/3$, so we conclude that the conditions of Corollary 1 are satisfied while those of (27) are not whenever \[ \frac{{\alpha _{1}}+{\alpha _{2}}}{2}<{\beta _{1}}\le \frac{{\alpha _{1}}+{\alpha _{2}}}{2}+\frac{{\alpha _{1}}-{\alpha _{2}}}{6}.\]

On the other hand, a simple example where (27) applies, but Corollary 1 does not, is when ${\alpha _{3}}<{\beta _{1}}<{\alpha _{2}}$, since this implies $P({\beta _{1}})>0$.## 4 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\ge 2$ we let ${\boldsymbol{L}_{t}}={({L_{t}^{1}},\dots ,{L_{t}^{d}})^{\top }}$, $t\in \mathbb{R}$, be a 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 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

*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 \[ {\int _{0}^{\infty }}{t^{2}}\hspace{0.1667em}|{\phi _{jk}}|(\operatorname{d}t)<\infty ,\hspace{2em}j,k=1,\dots ,d.\]

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 ##### (28)

\[ {X_{t}^{j}}-{X_{s}^{j}}={\sum \limits_{k=1}^{d}}{\int _{s}^{t}}{\int _{[0,\infty )}}{X_{u-v}^{k}}\hspace{0.1667em}{\phi _{jk}}(\operatorname{d}v)\hspace{0.1667em}\operatorname{d}u+{L_{t}^{j}}-{L_{s}^{j}},\hspace{2em}s<t,\]*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##### (29)

\[ \operatorname{d}{\boldsymbol{X}_{t}}=-\boldsymbol{\lambda }{\boldsymbol{X}_{t}}\hspace{0.1667em}\operatorname{d}t+(\boldsymbol{\eta }\ast \boldsymbol{X})(t)\hspace{0.1667em}\operatorname{d}t+\operatorname{d}{\boldsymbol{L}_{t}},\hspace{2em}t\in \mathbb{R}.\]##### (30)

\[ {\boldsymbol{h}_{\phi }}(z)=z{\boldsymbol{I}_{d}}-{\int _{[0,\infty )}}{e^{-zt}}\hspace{0.1667em}\boldsymbol{\phi }(\operatorname{d}t),\hspace{2em}z\in {\mathbb{C}_{+}},\]##### (31)

\[ {\boldsymbol{X}_{t}}={\int _{-\infty }^{t}}{\boldsymbol{g}_{\phi }}(t-s)\hspace{0.1667em}\operatorname{d}{\boldsymbol{L}_{s}},\hspace{2em}t\in \mathbb{R},\]##### (32)

\[ {\int _{0}^{\infty }}{e^{-ity}}{\boldsymbol{g}_{\phi }}(t)\hspace{0.1667em}\operatorname{d}t={\boldsymbol{h}_{\phi }}{(iy)^{-1}},\hspace{2em}y\in \mathbb{R}.\]*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.##### Theorem 5.

*Let*${\boldsymbol{h}_{\phi }}$

*be defined as in*(30)

*and assume that*$\det ({\boldsymbol{h}_{\phi }}(z))\ne 0$

*for all*$z\in {\mathbb{C}_{+}}$

*. Suppose further that all the entries of*$\boldsymbol{\eta }$

*are nonnegative measures and that*$\boldsymbol{\lambda }$

*is an M-matrix. Then the entries of the unique solution*${({\boldsymbol{X}_{t}})_{t\in \mathbb{R}}}$

*to*(28)

*, which is given by*(9)

*, are nonnegative in the sense that*${X_{t}^{j}}\ge 0$

*for*$j=1,\dots ,d$

*almost surely for each fixed*$t\in \mathbb{R}$

*.*

##### Remark 6.

When comparing Theorem 5 to Theorem 3, it seems that we have to impose an additional assumption of $\boldsymbol{\lambda }$ being an

*M*-matrix, which was not needed in the univariate case. However, requiring that the one-dimensional quantity $\lambda =-\phi (\{0\})$ is an*M*-matrix would naturally be interpreted as a nonnegativity constraint, but this is automatically satisfied under the stated assumptions on ${h_{\phi }}$, 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.

##### Example 5.

Let $\boldsymbol{\lambda }\in {\mathbb{R}^{p\times p}}$ with a spectrum $\sigma (\boldsymbol{\lambda })$ contained in $\{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)>0\}$, and consider the Ornstein–Uhlenbeck equation:

\[ \operatorname{d}{\boldsymbol{X}_{t}}=-\boldsymbol{\lambda }{\boldsymbol{X}_{t}}\hspace{0.1667em}\operatorname{d}t+\operatorname{d}{\boldsymbol{L}_{t}},\hspace{2em}t\in \mathbb{R}.\]

The stationary solution ${({\boldsymbol{X}_{t}})_{t\in \mathbb{R}}}$ takes the form ${\boldsymbol{X}_{t}}={\textstyle\int _{-\infty }^{t}}{e^{-\boldsymbol{\lambda }(t-s)}}\hspace{0.1667em}\operatorname{d}{\boldsymbol{L}_{s}}$ and thus, to ensure that ${\boldsymbol{X}_{t}}$ is nonnegative (and without further restrictions on the driving Lévy process ${({\boldsymbol{L}_{t}})_{t\in \mathbb{R}}}$), it must be required that all entries of ${e^{-\boldsymbol{\lambda }t}}$ are nonnegative for every $t\ge 0$. According to Lemma 2 this is the case if and only if $\boldsymbol{\lambda }$ is an *M*-matrix, which is exactly the condition of Theorem 5 for nonnegativity.##### Example 6.

Let ${\boldsymbol{A}_{1}},{\boldsymbol{A}_{2}},{\boldsymbol{B}_{0}}\in {\mathbb{R}^{d\times d}}$ be given such that both

\[ z\mapsto \det ({\boldsymbol{I}_{d}}{z^{2}}+{\boldsymbol{A}_{1}}z+{\boldsymbol{A}_{2}})\hspace{1em}\text{and}\hspace{1em}z\mapsto \det ({\boldsymbol{I}_{d}}z+{\boldsymbol{B}_{0}})\]

have no roots on ${\mathbb{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 \[ {\boldsymbol{Y}_{t}}={\int _{-\infty }^{t}}\boldsymbol{g}(t-s)\hspace{0.1667em}\operatorname{d}{\boldsymbol{L}_{s}},\hspace{2em}t\in \mathbb{R},\]

where the Fourier transform of $\boldsymbol{g}:[0,\infty )\to {\mathbb{R}^{d\times d}}$ is \[ {\int _{0}^{\infty }}{e^{-ity}}\boldsymbol{g}(t)\hspace{0.1667em}\operatorname{d}t={(-{\boldsymbol{I}_{d}}{y^{2}}+{\boldsymbol{A}_{1}}iy+{\boldsymbol{A}_{2}})^{-1}}({\boldsymbol{I}_{d}}iy+{\boldsymbol{B}_{0}}),\hspace{2em}y\in \mathbb{R}.\]

(See also [12, Theorem 3.22].) By [3, Theorem 4.8], ${({\boldsymbol{Y}_{t}})_{t\in \mathbb{R}}}$ can also be identified as the unique solution to (29) with $\boldsymbol{\lambda }={\boldsymbol{A}_{1}}-{\boldsymbol{B}_{0}}$ and \[ \boldsymbol{\eta }(\operatorname{d}t)={e^{-{\boldsymbol{B}_{0}}t}}\hspace{0.1667em}\operatorname{d}t\hspace{0.1667em}[{\boldsymbol{B}_{0}}({\boldsymbol{A}_{1}}-{\boldsymbol{B}_{0}})-{\boldsymbol{A}_{2}}].\]

In view of Theorem 5 and Lemma 2 we conclude that the entries of ${({\boldsymbol{Y}_{t}})_{t\in \mathbb{R}}}$ are nonnegative if both ${\boldsymbol{B}_{0}}$ and ${\boldsymbol{A}_{1}}-{\boldsymbol{B}_{0}}$ are *M*-matrices, and all entries of the matrix ${\boldsymbol{B}_{0}}({\boldsymbol{A}_{1}}-{\boldsymbol{B}_{0}})-{\boldsymbol{A}_{2}}$ are nonnegative.## 5 Proofs

##### Proof of Theorem 2.

The equivalence between (i) and (ii) is an immediate consequence of stationarity of ${({X_{t}})_{t\in \mathbb{R}}}$. The fact that (iii) implies (ii) is rather obvious—this is also readily seen by representing the subordinator ${({L_{t}})_{t\in \mathbb{R}}}$ as a cumulative sum of its positive jumps plus a nonnegative drift. To see that (ii) implies (iii), note by [15] that ${X_{t}}$ is infinitely divisible with a Lévy measure ${\nu _{X}}$ given by Since ${({L_{t}})_{t\in \mathbb{R}}}$ is a subordinator and

\[ {\nu _{X}}(B):=\big(\text{Leb}\times \nu \big)\big(\{(s,x)\in [0,\infty )\times \mathbb{R}\hspace{0.1667em}:\hspace{0.1667em}x{g_{\phi }}(s)\in B\setminus \{0\}\}\big),\]

where $\text{Leb}(\hspace{0.2222em}\cdot \hspace{0.2222em})$ refers to the Lebesgue measure. Since ${X_{t}}\ge 0$ almost surely, ${\nu _{X}}$ is concentrated on $(0,\infty )$, and hence we deduce that ##### (33)

\[ 0={\nu _{X}}((-\infty ,0))\ge \text{Leb}(\{s\in [0,\infty )\hspace{0.1667em}:\hspace{0.1667em}{g_{\phi }}(s)<0\})\nu ((0,\infty )).\]*ν*is nonzero, $\nu ((0,\infty ))>0$. By combining this observation with (33) we conclude that ${g_{\phi }}$ is nonnegative almost everywhere. In view of Bernstein’s theorem on monotone functions [6], which states that ${g_{\phi }}$ is nonnegative if and only if its Laplace transform is completely monotone on $[0,\infty )$, we have that (iii) holds if and only if (iv) does, and this completes the proof. □##### Proof of Corollary 1.

Since

*f*takes the form (23) and $\deg (Q)=2$, it follows from [21, Remark 2] that, necessarily, $\operatorname{Im}({\beta _{1}})=\operatorname{Im}({\beta _{2}})=0$ (as we also discussed in relation to (22)). If ${\beta _{1}}>{\beta _{2}}$, we may use (24) and represent*f*explicitly as \[ f(t)=\frac{1}{{\beta _{1}}-{\beta _{2}}}\big(P({\beta _{2}}){e^{{\beta _{2}}t}}-P({\beta _{1}}){e^{{\beta _{1}}t}}\big),\hspace{2em}t\ge 0.\]

Here we have also used the fact that ${Q^{\prime }}({\beta _{1}})=-{Q^{\prime }}({\beta _{2}})={\beta _{1}}-{\beta _{2}}$. By considering $t=0$ and $t\to \infty $ we conclude that $f(t)\ge 0$ for all *t*if and only if both $P({\beta _{1}})\le P({\beta _{2}})$ and $P({\beta _{1}})\le 0$, and this thus establishes (i). If ${\beta _{1}}={\beta _{2}}=:\beta $, we observe initially that \[\begin{aligned}{}R(z)& =(z+{a_{1}}+2\beta ){(z-\beta )^{2}}-P(z)\\ {} & =2{\beta ^{3}}+{a_{1}}{\beta ^{2}}-{a_{3}}-(3{\beta ^{2}}+2\beta {a_{1}}-{a_{2}})z\\ {} & ={P^{\prime }}(\beta )\beta -P(\beta )-{P^{\prime }}(\beta )z.\end{aligned}\]

Since the function having Fourier transform $y\mapsto {(iy-\beta )^{-2}}$ is $t\mapsto t{e^{\beta t}}$, it follows from (23) and the above expression for *R*that \[ f(t)=({P^{\prime }}(\beta )\beta -P(\beta ))t{e^{\beta t}}-{P^{\prime }}(\beta )(\beta t{e^{\beta t}}+{e^{\beta t}})=-(P(\beta )t+{P^{\prime }}(\beta )){e^{\beta t}}.\]

By considering $t=0$ and $t\to \infty $ once again we deduce that *f*is nonnegative on $[0,\infty )$ if and only if $P(\beta )\le 0$ and ${P^{\prime }}(\beta )\le 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

\[ {(\boldsymbol{\mu }\ast \boldsymbol{f})_{jk}}=\sum \limits_{l}\int {f_{lk}}(\hspace{0.2222em}\cdot \hspace{0.2222em}-t)\hspace{0.1667em}{\mu _{jl}}(\operatorname{d}t)\]

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.##### Lemma 1.

*Suppose that*$\det ({\boldsymbol{h}_{\phi }}(z))\ne 0$

*for all*$z\in {\mathbb{C}_{+}}$

*and let*${\boldsymbol{g}_{\phi }}:[0,\infty )\to {\mathbb{R}^{d\times d}}$

*be characterized by*(32)

*. Then it holds that*

\[ {\boldsymbol{g}_{\phi }}(t)={\boldsymbol{g}_{\phi }}(t-s){\boldsymbol{g}_{\phi }}(s)+{\int _{s}^{t}}{\boldsymbol{g}_{\phi }}(t-u)(\boldsymbol{\phi }\ast ({\boldsymbol{g}_{\phi }}{\mathbb{1}_{[0,s]}}))(u)\hspace{0.1667em}\operatorname{d}u\]

*for arbitrary*$s,t\in [0,\infty )$

*with*$s<t$

*.*

##### Proof.

To lighten notation, and in line with Remark 4, we denote the Laplace transform by $\mathcal{L}$. Specifically, for a given signed measure where ${({\boldsymbol{g}_{\phi }}\ast \boldsymbol{\phi })_{jk}}:={\textstyle\sum _{l=1}^{d}}{\textstyle\int _{[0,\infty )}}{g_{\phi ,jl}}(\hspace{0.2222em}\cdot \hspace{0.2222em}-t)\hspace{0.1667em}{\phi _{lk}}(\operatorname{d}t)$. Since $\mathcal{L}[{\boldsymbol{g}_{\phi }}](z)={\boldsymbol{h}_{\phi }}{(z)^{-1}}$ by (32), it is easy to check that $\mathcal{L}[{\boldsymbol{g}_{\phi }}]$ and $\mathcal{L}[\boldsymbol{\phi }]$ commute, and hence For a fixed $s\in [0,\infty )$, it thus follows from (34) and (35) that By multiplying $\mathcal{L}[{\boldsymbol{g}_{\phi }}](z)$ from the left on both sides of (36) and using uniqueness of the Laplace transform, this shows that

*μ*and an integrable function $f:[0,\infty )\to \mathbb{R}$, set $\mathcal{L}[\mu ](z)=\mathcal{L}[\mu (\operatorname{d}t)](z):={\textstyle\int _{[0,\infty )}}{e^{-zt}}\hspace{0.1667em}\mu (\operatorname{d}t)$ for $z\in {\mathbb{C}_{+}}$ and $\mathcal{L}[f]:=\mathcal{L}[f(t)\hspace{0.1667em}\operatorname{d}t]$. From [3, Proposition 5.1] it follows that##### (34)

\[ {\boldsymbol{g}_{\phi }}(t)={\boldsymbol{g}_{\phi }}(s)+{\int _{s}^{t}}({\boldsymbol{g}_{\phi }}\ast \boldsymbol{\phi })(u)\hspace{0.1667em}\operatorname{d}u,\hspace{2em}t>s\ge 0,\]##### (35)

\[ {\boldsymbol{g}_{\phi }}\ast \boldsymbol{\phi }=\boldsymbol{\phi }\ast {\boldsymbol{g}_{\phi }}.\] \[\begin{aligned}{}& \mathcal{L}[{\boldsymbol{g}_{\phi }}{\mathbb{1}_{(s,\infty )}}](z)\\ {} & \hspace{1em}=\mathcal{L}[{\mathbb{1}_{(s,\infty )}}](z){\boldsymbol{g}_{\phi }}(s)+\mathcal{L}\Big[{\int _{0}^{\hspace{0.2222em}\cdot \hspace{0.2222em}}}(\boldsymbol{\phi }\ast {\boldsymbol{g}_{\phi }})(u){\mathbb{1}_{(s,\infty )}}(u)\hspace{0.1667em}\operatorname{d}u\Big](z)\\ {} & \hspace{1em}=\frac{1}{z}\big(\mathcal{L}[{\delta _{s}}](z){\boldsymbol{g}_{\phi }}(s)+\mathcal{L}[(\boldsymbol{\phi }\ast ({\boldsymbol{g}_{\phi }}{\mathbb{1}_{[0,s]}})){\mathbb{1}_{(s,\infty )}}](z)+\mathcal{L}[\boldsymbol{\phi }](z)\mathcal{L}[{\boldsymbol{g}_{\phi }}{\mathbb{1}_{(s,\infty )}}](z)\big)\end{aligned}\]

for $z\in \mathbb{C}$ with $\operatorname{Re}(z)>0$. After rearranging terms we obtain ##### (36)

\[ {\boldsymbol{h}_{\phi }}(z)\mathcal{L}[{\boldsymbol{g}_{\phi }}{\mathbb{1}_{(s,\infty )}}](z)=\mathcal{L}[{\delta _{s}}](z){\boldsymbol{g}_{\phi }}(s)+\mathcal{L}[(\boldsymbol{\phi }\ast ({\boldsymbol{g}_{\phi }}{\mathbb{1}_{[0,s]}})){\mathbb{1}_{(s,\infty )}}](z).\] \[\begin{aligned}{}{\boldsymbol{g}_{\phi }}(t)& =({\boldsymbol{g}_{\phi }}\ast {\delta _{s}})(t){\boldsymbol{g}_{\phi }}(s)+({\boldsymbol{g}_{\phi }}\ast (\boldsymbol{\phi }\ast ({\boldsymbol{g}_{\phi }}{\mathbb{1}_{[0,s]}})))(t)\\ {} & ={\boldsymbol{g}_{\phi }}(t-s){\boldsymbol{g}_{\phi }}(s)+{\int _{s}^{t}}{\boldsymbol{g}_{\phi }}(t-u)(\boldsymbol{\phi }\ast ({\boldsymbol{g}_{\phi }}{\mathbb{1}_{[0,s]}}))(u)\hspace{0.1667em}\operatorname{d}u\end{aligned}\]

for $t>s$, 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.##### Lemma 2.

*Let*$\boldsymbol{A}\in {\mathbb{R}^{d\times d}}$

*. If*$\boldsymbol{A}$

*is an M-matrix, then*${e^{-\boldsymbol{A}t}}$

*has nonnegative entries for all*$t\ge 0$

*. Conversely, if*${e^{-\boldsymbol{A}t}}$

*has nonnegative entries for all*$t\ge 0$

*and the spectrum*$\sigma (\boldsymbol{A})$

*of*$\boldsymbol{A}$

*belongs to*$\{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)>0\}$

*, then*$\boldsymbol{A}$

*is an invertible M-matrix.*

##### Proof.

If $\boldsymbol{A}$ is an

*M*-matrix, we have that ${e^{-\boldsymbol{A}t}}={e^{-\alpha t}}{e^{\boldsymbol{B}t}}$ for some $\alpha \in [0,\infty )$ and $\boldsymbol{B}\in {\mathbb{R}^{d\times d}}$ with nonnegative entries. Since ${e^{\boldsymbol{B}t}}$ has nonnegative entries, we conclude that the same holds for ${e^{-\boldsymbol{A}t}}$. If instead the entries of ${e^{-\boldsymbol{A}t}}$ are nonnegative and $\sigma (\boldsymbol{A})\subseteq \{z\in \mathbb{C}\hspace{0.1667em}:\hspace{0.1667em}\operatorname{Re}(z)>0\}$, it follows that the entries of \[ {\int _{0}^{\infty }}{e^{-xt}}{e^{-\boldsymbol{A}t}}\hspace{0.1667em}\operatorname{d}t={({\boldsymbol{I}_{d}}x+\boldsymbol{A})^{-1}}\]

are nonnegative for all $x\in [0,\infty )$ as well. By [14, Theorem 2] this implies that $\boldsymbol{A}$ is an invertible *M*-matrix. □##### Proof of Theorem 5.

For $\varepsilon \in (0,\infty )$ define the measure ${\boldsymbol{\phi }_{\varepsilon }}:=-\boldsymbol{\lambda }{\delta _{0}}+\boldsymbol{\eta }(\hspace{0.2222em}\cdot \hspace{0.2222em}\cap (\varepsilon ,\infty ))$. We start by observing that, when To see this suppose, for the sake of contradiction, there exist sequences ${({\varepsilon _{n}})_{n\ge 1}}\subseteq (0,\infty )$ and ${({z_{n}})_{n\ge 1}}\subseteq {\mathbb{C}_{+}}$ such that $\det ({\boldsymbol{h}_{{\phi _{{\varepsilon _{n}}}}}}({z_{n}}))=0$ for all $n\ge 1$ and ${\varepsilon _{n}}\to 0$ as $n\to \infty $. Note that the length of any entry of $\mathcal{L}[{\boldsymbol{\phi }_{\varepsilon }}](z)$ is bounded by the constant for all

*ε*is sufficiently small,##### (37)

\[ \det ({\boldsymbol{h}_{{\phi _{\varepsilon }}}}(z))\ne 0\hspace{2em}\text{for all}\hspace{2.5pt}z\in {\mathbb{C}_{+}}.\]*ε*and*z*, and hence ${\inf _{\varepsilon }}|\det ({\boldsymbol{h}_{{\phi _{\varepsilon }}}}(z))|\sim |z{|^{d}}$ as $|z|\to \infty $ by the Leibniz formula (the notation ∼ means that the ratio tends to one). This shows in particular that the sequence ${({z_{n}})_{n\ge 1}}$ must be bounded, and so it has a subsequence ${({z_{{n_{k}}}})_{k\ge 1}}$ which converges to a point ${z^{\ast }}\in {\mathbb{C}_{+}}$. However, this would imply that \[ \det ({\boldsymbol{h}_{\phi }}({z^{\ast }}))=\underset{k\to \infty }{\lim }\det ({\boldsymbol{h}_{{\phi _{{\varepsilon _{{n_{k}}}}}}}}({z_{{n_{k}}}}))=0,\]

which contradicts the original assumption that $\det ({\boldsymbol{h}_{\phi }}(z))\ne 0$ for all $z\in {\mathbb{C}_{+}}$. Consequently, the condition (37) is satisfied as long as *ε*is smaller than a certain threshold, say, ${\varepsilon ^{\ast }}$. Thus, for $\varepsilon \le {\varepsilon ^{\ast }}$ we can define the corresponding function ${\boldsymbol{g}_{{\phi _{\varepsilon }}}}:[0,\infty )\to {\mathbb{R}^{d\times d}}$ through (32) (with $\boldsymbol{\phi }$ replaced by ${\boldsymbol{\phi }_{\varepsilon }}$). By (34) and (35), \[\begin{aligned}{}{\boldsymbol{g}_{{\phi _{\varepsilon }}}}(t)& ={\boldsymbol{I}_{d}}+{\int _{0}^{t}}({\boldsymbol{\phi }_{\varepsilon }}\ast {\boldsymbol{g}_{{\phi _{\varepsilon }}}})(s)\hspace{0.1667em}\operatorname{d}s\\ {} & ={\boldsymbol{I}_{d}}-\boldsymbol{\lambda }{\int _{0}^{t}}{\boldsymbol{g}_{{\phi _{\varepsilon }}}}(s)\hspace{0.1667em}\operatorname{d}s\end{aligned}\]

for $t\in [0,\varepsilon ]$. By uniqueness of solutions of such differential equation it follows that for ${\boldsymbol{g}_{{\phi _{\varepsilon }}}}(t)={e^{-\boldsymbol{\lambda }t}}$ for $t\in [0,\varepsilon ]$, in which case we can rely on Lemma 2 to deduce that its entries are nonnegative. Now, given that ${\boldsymbol{g}_{{\phi _{\varepsilon }}}}(t)$ has nonnegative entries for $t\in [0,k\varepsilon ]$ for some positive integer $k\ge 1$, we can apply Lemma 1 with $s=k\varepsilon $ to deduce that this remains true when $t\in (k\varepsilon ,(k+1)\varepsilon ]$. Consequently, it follows by induction that ${\boldsymbol{g}_{{\phi _{\varepsilon }}}}(t)$ has nonnegative entries for all $t\in [0,\infty )$.Next, for $\varepsilon \in [0,{\varepsilon ^{\ast }}]$ define the Fourier transform Since ${\mathcal{F}_{{\phi _{\varepsilon }},jk}}$ converges pointwise to ${\mathcal{F}_{{\phi _{0}},jk}}$ (the Fourier transform of the $(j,k)$-th entry ${g_{\phi ,jk}}$ of ${\boldsymbol{g}_{\phi }}$) as $\varepsilon \to 0$, it follows by the Plancherel theorem and (38) that

\[ {\mathcal{F}_{{\phi _{\varepsilon }},jk}}(y):={({\boldsymbol{h}_{{\phi _{\varepsilon }}}}{(iy)^{-1}})_{jk}},\hspace{2em}y\in \mathbb{R},\]

of the $(j,k)$-th entry ${g_{{\phi _{\varepsilon }},jk}}$ of ${\boldsymbol{g}_{{\phi _{\varepsilon }}}}$. By Cramer’s rule it follows that ${\mathcal{F}_{{\phi _{\varepsilon }},jk}}(y)=\det (\widehat{{\boldsymbol{h}_{{\phi _{\varepsilon }}}}(iy)})/\det ({\boldsymbol{h}_{{\phi _{\varepsilon }}}}(iy))$, where $\widehat{{\boldsymbol{h}_{{\phi _{\varepsilon }}}}(iy)}$ is the matrix formed by replacing the *j*th column of ${\boldsymbol{h}_{{\phi _{\varepsilon }}}}(iy)$ by the*k*th canonical basis vector. By the same type of arguments as above (in particular, by relying on the Leibniz formula), \[\begin{aligned}{}& \underset{\varepsilon \le {\varepsilon ^{\ast }}}{\inf }|\det ({\boldsymbol{h}_{{\phi _{\varepsilon }}}}(iy))|\ge {c_{1}}(1\vee |y{|^{d}})\\ {} \text{and}\hspace{1em}& \underset{\varepsilon \le {\varepsilon ^{\ast }}}{\sup }|\det (\widehat{{\boldsymbol{h}_{{\phi _{\varepsilon }}}}(iy)})|\le {c_{2}}(1\vee |y{|^{d-1}})\end{aligned}\]

for suitably chosen constants ${c_{1}},{c_{2}}\in (0,\infty )$. Consequently, we conclude that ##### (38)

\[ \underset{0\le \varepsilon \le {\varepsilon ^{\ast }}}{\sup }|{\mathcal{F}_{{\phi _{\varepsilon }},jk}}(y)|\le {c_{1}^{-1}}{c_{2}}(1\wedge |y{|^{-1}}),\hspace{2em}y\in \mathbb{R}.\] \[ {\int _{0}^{\infty }}{({g_{{\phi _{\varepsilon }},jk}}(t)-{g_{\phi ,jk}}(t))^{2}}\hspace{0.1667em}\operatorname{d}t={\int _{\mathbb{R}}}|{\mathcal{F}_{{\phi _{\varepsilon }},jk}}(y)-{\mathcal{F}_{{\phi _{0}},jk}}(y){|^{2}}\hspace{0.1667em}\operatorname{d}y\to 0,\hspace{2em}\varepsilon \to 0.\]

From this we establish that ${g_{{\phi _{{\varepsilon _{n}}}},jk}}\to {g_{\phi ,jk}}$ almost everywhere for a suitable sequence ${({\varepsilon _{n}})_{n\ge 1}}\subseteq (0,{\varepsilon ^{\ast }}]$ with ${\varepsilon _{n}}\to 0$ as $n\to \infty $. This shows that ${g_{\phi ,jk}}(t)\ge 0$ for almost all $t\in [0,\infty )$ and, thus, completes the proof. □