State Spaces of Multifactor Approximations of Nonnegative Volterra Processes

1. Introduction

Multifactor Markovian approximations for fractional Brownian motion were initially introduced by Carmona and Coutin (1998) and revisited more recently by Abi Jaber and El Euch (2019b) in the context of nonnegative stochastic Volterra equations, motivated by rough and Volterra Heston models of El Euch and Rosenbaum (2019) and Abi Jaber et al. (2019b). Since then, substantial literature has emerged on such multifactor processes for numerical approximation methods (Harms 2019; Chevalier et al. 2022; Bayer and Breneis 2023a, b, 2024; Alfonsi and Kebaier 2024), deep learning approaches (Papapantoleon and Rou 2025), modeling (Abi Jaber 2019), and optimal control (Abi Jaber et al. 2021). It should be noted that such approximations are also heavily used in physics, chemistry, and other fields; see, for example, Baczewski and Bond (2013) and Bochud and Challet (2007).

The starting point is a nonnegative solution to the stochastic Volterra equation

where the kernel K is (approximated by) a weighted sum of exponentials of the formwith positive nodes $\mathit{x}={(x_i)}_{i=1,\dots ,N}$ and weights $\mathit{w}={(w_i)}_{i=1,\dots ,N}$. Such series in terms of exponential functions are sometimes known as Prony series. The coefficients $b,\sigma :\mathbb{R}\to \mathbb{R}$ are continuous and satisfy the boundary conditionsto ensure that the process Y remains nonnegative for any $Y_0\ge 0$.

Then, the nonnegative Volterra process Y can be written in the form $Y_t={\sum }_{i=1}^N{w}_i{Y}_t^{(i)}$, where $\mathit{Y}={(Y^{(i)})}_{i=1,\dots ,N}$ is the solution to the N-dimensional Stochastic Differential Equation (SDE)

with initial values $Y_0^{(i)}$, such that ${\sum }_{i=1}^N{w}_i{Y}_0^{(i)}=Y_0$. More generally, we will consider the N-dimensional SDEwith $y_0^{(i)}\in \mathbb{R}$, which are often, but not necessarily, chosen to coincide with the initial values $Y_0^{(i)}$, for $i=1,\dots ,N$.

The aim of the paper is to determine a state space of the multifactor Markovian process $\mathit{Y}$. That is, we want to determine a set $\mathcal{D}\subseteq {\mathbb{R}}^N$ such that for every starting value ${\mathit{Y}}_0\in \mathcal{D}$, there exists a $\mathcal{D}$-valued solution $\mathit{Y}$ to (1)—that is ${\mathit{Y}}_t\in \mathcal{D}$ for all $t\ge 0$ almost surely.

Beyond the mathematical importance of defining the state space of the Markovian process $\mathit{Y}$, the knowledge of the state space is crucial for several practical applications, some of which have been considered so far:

• Modeling and Calibration: The Multifactor Model (1) can serve as a model in its own right (and not solely as an approximation of Volterra models), usually for stochastic volatility factors, as seen in the lifted Heston model of Abi Jaber (2019). Here, the knowledge of the state space is crucial for calibrating the initial values of ${\mathit{Y}}_0$ directly to market data.

• Simulation Accuracy: The identification of a valid state space allows for more precise simulation schemes for $\mathit{Y}$. For instance, Bayer and Breneis (2024) generalized a simulation scheme for the square-root process of Lileika and Mackevičius (2021) to simulate paths from the dynamics in (1) with $\sigma (z)=\sqrt{z}$. However, the authors were not able to show that their simulation scheme is well-defined because they did not know the state space. Indeed, when simulating from (1), great care has to be taken to ensure that the aggregated process Y does not become negative, or if it does, one has to determine how to proceed with the square root term $\sqrt{Y}$.

• Efficient Domain Meshing for Partial Differential Equation (PDE) Solutions: In Papapantoleon and Rou (2025), for example, the authors use a rejection algorithm that discards simulations failing to satisfy the condition ${\sum }_i{w}_i{Y}_0^{(i)}\ge 0$; this approach is not precise and becomes inefficient in high dimensions.

For all these reasons, the geometry of the state space for multifactor processes $\mathit{Y}$ quickly became a central focus in the associated literature for nonnegative Volterra processes. For the lifted Heston model of Abi Jaber (2019), simulations demonstrated that although individual processes $Y^{(i)}$ might take negative values, their aggregated sum Y remains nonnegative. In this context, two key papers provided abstract characterizations of possible state spaces: one based on the resolvent of the first kind of the kernel by Abi Jaber and El Euch (2019a) and the other using the resolvent of the second kind of the kernel by Cuchiero and Teichmann (2020) (we refer to Section 5 for further details on these state spaces and their connections). Although these spaces are valid for a wide range of kernels, they remain somewhat abstract and challenging to make explicit, making it almost impossible to determine the good conditions on the initial values $\mathit{Y}\ge 0$ of the process $\mathit{Y}$. Finally, using the simulation algorithm for the rough Heston model due to Bayer and Breneis (2024), we visually represent the process’s domain by a sample plot; see Figure 1—replicating a similar plot already presented in that paper. One can clearly recognize that the sample paths do not only seem to lie in the half-plane $Y=w_1{Y}^{(1)}+w_2{Y}^{(2)}\ge 0$ marked with the downward-oriented black line, but also in an even smaller cone seemingly below the upward-oriented line $\{\mathit{y}\in {\mathbb{R}}^2:y_1\ge y_2\}$.

1.1. Main Contributions

The main question we are interested in can be summarized as follows:

What constitutes a suitable state space for the multifactor Markovian process $\mathit{Y}$ in (1)?

Our main results in Theorem 2 and Corollary 1 establish that this state space can be represented as a linear transformation of ${\mathbb{R}}_{+}^N$, and we provide an explicit form for this transformation.

As a first application of this result, we prove that the weak simulation scheme for the rough Heston process proposed by Bayer and Breneis (2024) is well-defined in the sense that the variance process always stays nonnegative; see Section 3. In Section 4, we derive the corresponding pricing PDE on the transformed domain ${\mathbb{R}}_{+}^N$ and solve it numerically by the finite element method after truncation of the domain. Naturally, knowing the PDE’s precise domain is crucial for accurate numerical approximations. Finally, in Section 5, we show how the explicit formula for the domain compares with general, abstract characterizations given in the literature by Abi Jaber and El Euch (2019a) and Cuchiero and Teichmann (2020).

1.2. Notation and Conventions

We denote by ${\mathit{e}}_i$ the i-th unit vector, which has a 1 in the i-th component and 0 in every other component; $\mathbf{1}≔{(1,1,\dots ,1)}^{\top }$ is the vector with 1 in every component; $\text{Id}$ is the identity matrix; $\text{diag}(\mathit{a})$ for a vector $\mathit{a}\in {\mathbb{R}}^N$ is the diagonal matrix with entries $\mathit{a}$ in the diagonal; and $\overline{w}≔{\mathbf{1}}^{\top }\mathit{w}={\sum }_{i=1}^N{w}_i.$ Throughout, italic letters a denote real numbers and bold letters $\mathit{a}$ denote vectors, where we write $\mathit{a}={(a_i)}_{i=1}^N$ for the components of $\mathit{a}$. An exception are stochastic processes, where components are denoted by ${\mathit{Y}}_t={(Y_t^{(i)})}_{i=1}^N$ (due to the time variable in the subscript).

2. State Spaces of the Multifactor Markovian Process

Fix $N\ge 1$. We consider the N-dimensional SDE

where $b,\sigma :\mathbb{R}\to \mathbb{R}$ are continuous, satisfy the linear growth conditionand the Boundary Conditions (1). The speeds of mean-reversion $\mathit{x}={(x_i)}_{i=1,\dots ,N}$ are positive and ordered—that is, $0<x_1\le x_2\le \cdots \le x_N$—the weights $\mathit{w}={(w_i)}_{i=1,\dots ,N}$ are positive, W is a one-dimensional Brownian motion, and where ${\mathit{Y}}_0,{\mathit{y}}_0\in {\mathbb{R}}^N$ may be different. This corresponds to (1) written in vector form.

The aim of this section is to determine a state space of the multifactor process $\mathit{Y}$. That is, we want to determine a set $\mathcal{D}\subseteq {\mathbb{R}}^N$ such that for every starting value ${\mathit{Y}}_0\in \mathcal{D}$, there exists a $\mathcal{D}$-valued weak solution $\mathit{Y}$ to (6)—that is, ${\mathit{Y}}_t\in \mathcal{D}$ for all $t\ge 0$ almost surely. In particular, the domain $\mathcal{D}$ should be a subset of the half-plane $\{\mathit{y}\in {\mathbb{R}}^N:{\mathit{w}}^{\top }\mathit{y}\ge 0\}$ to ensure nonnegativity of the aggregated weighted process $Y≔{\mathit{w}}^{\top }\mathit{Y}$, which for the specific case ${\mathit{y}}_0={\mathit{Y}}_0$ would correspond to the Volterra Process (1) for the weighted sum of exponential kernel K given in (1). As illustrated in Figure 1, and following the abstract characterizations of domains of (possibly infinite-dimensional) lifts of nonnegative Volterra processes in Abi Jaber and El Euch (2019a) and Cuchiero and Teichmann (2020), we suspect the domain $\mathcal{D}$ to be a cone.

2.1. Main Result

We prove that the domain $\mathcal{D}$ is a cone characterized by the set $\mathcal{Q}$ of admissible matrices, which is defined as follows.

Definition 1.

A matrix $Q\in {\mathbb{R}}^{N\times N}$ is called admissible if it satisfies the following assumptions:

1. Q is invertible,

2. ${\mathit{e}}_N^{\top }Q={\mathit{w}}^{\top }$,

3. $Q\mathbf{1}=\overline{w}{\mathit{e}}_N$, where $\overline{w}≔{\mathit{w}}^{\top }\mathbf{1}$,

4. ${(Q\text{diag}(\mathit{x})Q^{-1})}_{i,j}\le 0$ for $i,j\in \{1,\dots ,N\}$ with $i\ne j$.

We denote by $\mathcal{Q}$ the set of all admissible matrices.

Before stating our main theorem, we first show that the set of admissible matrices $\mathcal{Q}$ is nonempty by providing an explicit example of an admissible matrix. However, $\mathcal{Q}$ is not reduced to a singleton, as shown in Example 3 below.

Theorem 1.

The matrix $Q={(q_{i,j})}_{i,j=1,\dots N}\in {\mathbb{R}}^{N\times N}$ given by

$$q_{i,j}=w_j,\hspace{1em}j\le i,\hspace{1em}\hspace{1em}{q}_{i,i+1}=-{\displaystyle \sum _{j=1}^i{w}_j},\hspace{1em}i=1,\dots ,N-1,$$

and zeros elsewhere is admissible. In particular, $\mathcal{Q}$ is nonempty.

Proof.

The proof is given in Section 2.4. □

We are now ready to state our main theorem that gives state spaces of the Multifactor Markovian Process (2).

Theorem 2.

Let $b,\sigma :\mathbb{R}\to \mathbb{R}$ be continuous functions satisfying the Linear Growth Condition (2) and the Boundary Conditions (1). Let $Q\in \mathcal{Q}$ be an admissible matrix in the sense of Definition 1 and suppose that

$${\mathit{y}}_0=\mu \text{diag}{(\mathit{x})}^{-1}\mathbf{1}\hspace{1em}\mathit{\text{for}} \mathit{\text{some}} \mu \ge 0.$$

Set $\mathcal{D}=Q^{-1}{\mathbb{R}}_{+}^N$. Then, for each ${\mathit{Y}}_0\in \mathcal{D}$, there exists a $\mathcal{D}$-valued weak solution $\mathit{Y}$ to (2).

Proof.

The proof is given in Section 2.5. □

Example 1.

For the admissible matrix Q given in Theorem 1, the set $\mathcal{D}=Q^{-1}{\mathbb{R}}_{+}^N$ corresponds to the set of $\mathit{y}\in {\mathbb{R}}^N$ such that ${\mathit{w}}^{\top }\mathit{y}\ge 0$ and

$${\displaystyle \sum _{j=1}^i{w}_j}{y}_j\ge {\displaystyle \sum _{j=1}^i{w}_j}{y}_{i+1}\hspace{1em}\hspace{1em}\text{for}\hspace{1em}\hspace{1em}i=1,\dots ,N-1.$$

Remark 1.

The domain $Q^{-1}{\mathbb{R}}_{+}^N$ is not unique; see the appendix.

In practice—for instance, for the multifactor approximations of Volterra processes—we often have ${\mathit{Y}}_0={\mathit{y}}_0$. Hence, it may be interesting to know whether ${\mathit{y}}_0$ as given in (2) is in $Q^{-1}{\mathbb{R}}_{+}^N$. We treat this question in a slightly more general context in the following lemma.

Lemma 1.

Let Q be an admissible matrix and let ${\mathit{y}}_0\in {\mathbb{R}}^N$ satisfy (2). Then, ${\mathit{y}}_0\in Q^{-1}{\mathbb{R}}_{+}^N.$

Proof.

We prove the equivalent statement $Q{\mathit{y}}_0\in {\mathbb{R}}_{+}^N$. First, note that

$$Q{\mathit{y}}_0=\mu Q\text{diag}{(\mathit{x})}^{-1}\mathbf{1}=\mu Q\text{diag}{(\mathit{x})}^{-1}{Q}^{-1}Q\mathbf{1}=\mu \overline{w}Q\text{diag}{(\mathit{x})}^{-1}{Q}^{-1}{\mathit{e}}_N.$$

Recall that in linear algebra, an M-matrix is a square matrix with nonpositive off-diagonal entries and with eigenvalues whose real parts are nonnegative. Clearly, the matrix $Q\text{diag}(\mathit{x})Q^{-1}$ is an M-matrix: the nonpositivity of the off-diagonal entries holds by assumption, and its eigenvalues are (real and) positive, because it is just the matrix $\text{diag}(\mathit{x})$ written in a different basis. Now, it is well-known that the inverse of an M-matrix has nonnegative entries (in fact, this property characterizes M-matrices). Therefore,

$${(Q\text{diag}(\mathit{x})Q^{-1})}^{-1}=Q\text{diag}{(\mathit{x})}^{-1}{Q}^{-1},$$

has only nonnegative entries. In particular,

$$Q{\mathit{y}}_0=\mu \overline{w}Q\text{diag}{(\mathit{x})}^{-1}{Q}^{-1}{\mathit{e}}_N\in {\mathbb{R}}_{+}^N,$$

proving the lemma. □

We remark that Condition (2) on ${\mathit{y}}_0$ can easily be dropped by using affine transformations of ${\mathbb{R}}_{+}^N$ instead of linear ones. We give the specific details in the next corollary.

Corollary 1.

Let $b,\sigma :\mathbb{R}\to \mathbb{R}$ be continuous functions satisfying the Linear Growth Condition (2) and the Boundary Conditions (1). Let $Q\in \mathcal{Q}$ be an admissible matrix in the sense of Definition 1, and assume that ${\mathit{w}}^{\top }{\mathit{y}}_0\ge 0$. Let ${\tilde{\mathit{y}}}_0$ be chosen according to (2) with ${\mathit{w}}^{\top }{\tilde{\mathit{y}}}_0={\mathit{w}}^{\top }{\mathit{y}}_0.$ Define the set

$$\mathcal{D}=Q^{-1}{\mathbb{R}}_{+}^N+({\mathit{y}}_0-{\tilde{\mathit{y}}}_0).$$(2)

Then, for any ${\mathit{Y}}_0\in \mathcal{D}$, there exists a $\mathcal{D}$-valued weak solution $\mathit{Y}$ to (2).

Proof.

Denote by $\tilde{\mathit{Y}}$ a weak solution to

$$\text{d}{\tilde{\mathit{Y}}}_t=-\text{diag}(\mathit{x})\left({\tilde{\mathit{Y}}}_t-{\tilde{\mathit{y}}}_0\right)\hspace{0.17em}\text{d}t+b({\mathit{w}}^{\top }{\tilde{\mathit{Y}}}_t)\mathbf{1}\hspace{0.17em}\text{d}t+\sigma ({\mathit{w}}^{\top }{\tilde{\mathit{Y}}}_t)\mathbf{1}\hspace{0.17em}\text{d}{W}_t,$$

with initial condition ${\tilde{\mathit{Y}}}_0≔{\mathit{Y}}_0+{\tilde{\mathit{y}}}_0-{\mathit{y}}_0.$ Note that ${\tilde{\mathit{Y}}}_0\in Q^{-1}{\mathbb{R}}_{+}^N$, so Theorem 2 implies the existence of such a solution $\tilde{\mathit{Y}}$ that stays in $Q^{-1}{\mathbb{R}}_{+}^N$. Define the process $\mathit{R}≔\tilde{\mathit{Y}}+{\mathit{y}}_0-{\tilde{\mathit{y}}}_0$, and note that ${\mathit{w}}^{\top }\mathit{R}={\mathit{w}}^{\top }\tilde{\mathit{Y}}$. Thus, $\mathit{R}$ satisfies ${\mathit{R}}_0={\mathit{Y}}_0$ and

$$\text{d}{\mathit{R}}_t=-\text{diag}(\mathit{x})\left({\mathit{R}}_t-{\mathit{y}}_0\right)\hspace{0.17em}\text{d}t+b({\mathit{w}}^{\top }{\mathit{R}}_t)\mathbf{1}\hspace{0.17em}\text{d}t+\sigma ({\mathit{w}}^{\top }{\mathit{R}}_t)\mathbf{1}\hspace{0.17em}\text{d}{W}_t.$$

But this means that $\mathit{R}$ is a solution to (2) that stays in $\mathcal{D}$. The corollary follows immediately. □

We now give the specific result for the multifactor square-root process.

Example 2.

Consider the multifactor square-root model

$$\text{d}{\mathit{V}}_t^N=-\text{diag}(\mathit{x})\left({\mathit{V}}_t^N-{\mathit{v}}_0\right)\hspace{0.17em}\text{d}t+\left(\theta -\lambda {\mathit{w}}^{\top }{\mathit{V}}_t^N\right)\mathbf{1}\hspace{0.17em}\text{d}t+\nu \sqrt{{\mathit{w}}^{\top }{\mathit{V}}_t^N}\mathbf{1}\hspace{0.17em}\text{d}{W}_t.$$

Let Q be an admissible matrix, and assume that ${\mathit{w}}^{\top }{\mathit{v}}_0\ge 0$. Then,

$$\mathcal{D}=Q^{-1}{\mathbb{R}}_{+}^N+\left({\mathit{v}}_0-\frac{{\mathit{w}}^{\top }{\mathit{v}}_0}{{\mathit{w}}^{\top }\text{diag}{(\mathit{x})}^{-1}\mathbf{1}}\text{diag}{(\mathit{x})}^{-1}\mathbf{1}\right).$$

2.2. Link with Nonnegative Volterra Processes

As an application of our result, one can obtain the existence of nonnegative solutions to Volterra equations with kernels of the Form (1). We note that such existence can be obtained by working directly on the level of the Volterra equation as done in Abi Jaber et al. (2019b, theorem 3.6 and example 3.7). Here, our result provides another alternative, as illustrated in the following corollary.

Corollary 2.

Let $b,\sigma :\mathbb{R}\to \mathbb{R}$ be continuous functions satisfying the Linear Growth Condition (2) and the Boundary Conditions (1). Let the kernel K be given by a weighted sum of exponentials as in (1). Then, for each $Y_0\ge 0$, the Stochastic Volterra Equation (1) admits a nonnegative weak solution Y.

Proof.

Fix $Y_0\ge 0$, and let ${\mathit{y}}_0\in {\mathbb{R}}^N$ be such that ${\mathit{w}}^{\top }{\mathit{y}}_0=Y_0$. Let ${\tilde{\mathit{y}}}_0$ be chosen according to (2) with ${\mathit{w}}^{\top }{\tilde{\mathit{y}}}_0={\mathit{w}}^{\top }{\mathit{y}}_0$. Let $Q\in \mathcal{Q}$ be an admissible matrix—for instance, given by Theorem 1. Then, it follows from Lemma 1 that ${\tilde{\mathit{y}}}_0\in Q^{-1}\mathbb{R}$. Hence, ${\mathit{y}}_0={\tilde{\mathit{y}}}_0+({\mathit{y}}_0-{\tilde{\mathit{y}}}_0)\in \mathcal{D}$, with $\mathcal{D}$ given by (2). An application of Corollary 1, with the starting value $Y_0={\mathit{y}}_0\in \mathcal{D}$, yields the existence of a $\mathcal{D}$-valued weak solution $\mathit{Y}$ to the Equation (2). Thanks to the variation of constants formula, we can rewrite the equation in the form

$${\mathit{Y}}_t={\mathit{y}}_0+{\displaystyle {\int }_0^t\text{exp}}(-\text{diag}(\mathit{x})(t-s))\mathbf{1}(b({\mathit{w}}^{\top }{\mathit{Y}}_s)\hspace{0.17em}\text{d}s+\sigma ({\mathit{w}}^{\top }{\mathit{Y}}_s)\hspace{0.17em}\text{d}{W}_s),$$

so that the process Y defined by $Y={\mathit{w}}^{\top }\mathit{Y}$ solves the equation

$$Y_t=Y_0+{\displaystyle {\int }_0^t{\displaystyle \sum _{i=1}^N{w}_i}}{e}^{-x_i(t-s)}(b(Y_s)\hspace{0.17em}\text{d}s+\sigma (Y_s)\hspace{0.17em}\text{d}{W}_s),$$

which is precisely the Volterra Equation (1) with the kernel K given by (1). It remains to argue that, for all $t\ge 0$, $Y_t$ remains nonnegative by using the fact that ${\mathit{Y}}_t\in \mathcal{D}$. Indeed, using Definition 1, condition (2) and the fact that ${\mathit{w}}^{\top }{\tilde{\mathit{y}}}_0={\mathit{w}}^{\top }{\mathit{y}}_0$, we obtain that

$${\mathit{w}}^{\top }\mathcal{D}={\mathit{e}}_N^{\top }Q{Q}^{-1}{\mathbb{R}}_{+}^N+({\mathit{w}}^{\top }{\mathit{y}}_0-{\mathit{w}}^{\top }{\tilde{\mathit{y}}}_0)={\mathit{e}}_N^{\top }{\mathbb{R}}_{+}^N={\mathbb{R}}_{+}.$$

Hence, for all $t\ge 0$, $Y_t={\mathit{w}}^{\top }{\mathit{Y}}_t\in {\mathit{w}}^{\top }\mathcal{D}={\mathbb{R}}_{+}$, which ends the proof. □

2.3. On Admissible Matrices for $N\in \{2,3\}$

In this section, we give examples of admissible matrices.

Example 3.

In the case $N=2$, in Definition 1, conditions (2) and (3) imply that we are looking for a matrix of the form

$$\begin{array}{l}Q=\left(\begin{array}{cc}q& -q\\ w_1& w_2\end{array}\right),\hfill \end{array}$$

for some $q\ne 0$ to ensure invertibility. Then,

$$\begin{array}{l}Q\text{diag}(\mathit{x})Q^{-1}=\frac{1}{\overline{w}}\left(\begin{array}{cc}{w}_1{x}_2+w_2{x}_1& (x_1-x_2)q\\ w_1{w}_2(x_1-x_2)q^{-1}& w_1{x}_1+w_2{x}_2\end{array}\right).\hfill \end{array}$$

Because $x_1\le x_2$, the last condition in Definition 1 is satisfied for any $q>0$, and indeed, the domain $Q^{-1}{\mathbb{R}}_{+}^2$ is independent of the precise choice of q given by

For the case of the Multifactor Square-Root Process (2), the resulting sample paths of ${\mathit{V}}^2$ and $\mathit{U}≔Q{\mathit{V}}^2$ are illustrated in Figure 2. Note that we chose the large maturity $T=100$ to give the process more time to explore its domain. Thereby, it is more clearly visible that the domain of $\mathit{U}$ is indeed ${\mathbb{R}}_{+}^2$ than if we had set $T=1$.

For $N=3$, a similar computation—relegated to the appendix due to its length—gives multiple choices of domains. See the appendix for details.

2.4. Proof of Theorem 1

In preparation for the proof, we introduce the matrix $R={(r_{i,j})}_{i,j=1,\dots ,N}$ defined by

that will turn out to be the inverse of Q. For example, for $N=4$, we have

Proof of Theorem 1.

In Definition 1, conditions (2) and (3) are readily satisfied by construction. To argue condition (1), we will prove that R is actually the inverse of Q—that is, $QR=\text{Id}$. Indeed, consider first the diagonal elements. Here, we have

$$\begin{array}{l}{(QR)}_{NN}={\displaystyle \sum _{k=1}^N{q}_{Nk}}{r}_{kN}={\displaystyle \sum _{k=1}^N{w}_k}\frac{1}{\overline{w}}=1,\\ {(QR)}_{ii}={\displaystyle \sum _{k=1}^N{q}_{ik}}{r}_{ki}={\displaystyle \sum _{k=1}^i{w}_k}\frac{w_{i+1}}{\sum _{\mathcal{\ell }=1}^i{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{i+1}{w}_{\mathcal{\ell }}}+\left(-{\displaystyle \sum _{\mathcal{\ell }=1}^i{w}_{\mathcal{\ell }}}\right)\frac{-1}{\sum _{\mathcal{\ell }=1}^{i+1}{w}_{\mathcal{\ell }}}=1,\end{array}$$

for $i=1,\dots ,N-1$. Next, consider off-diagonal elements. We have

$${(QR)}_{Nj}={\displaystyle \sum _{k=1}^N{q}_{Nk}}{r}_{kj}={\displaystyle \sum _{k=1}^j{w}_k}\frac{w_{j+1}}{\sum _{\mathcal{\ell }=1}^j{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}+w_{j+1}\frac{-1}{\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}=0,$$

for $j\le N-1$,

$${(QR)}_{iN}={\displaystyle \sum _{k=1}^N{q}_{ik}}{r}_{kN}={\displaystyle \sum _{k=1}^i{w}_k}\frac{1}{\overline{w}}+\left(-{\displaystyle \sum _{\mathcal{\ell }=1}^i{w}_{\mathcal{\ell }}}\right)\frac{1}{\overline{w}}=0,$$

for $i\le N-1$,

$${(QR)}_{ij}={\displaystyle \sum _{k=1}^N{q}_{ik}}{r}_{kj}={\displaystyle \sum _{k=1}^i{w}_k}\frac{w_{j+1}}{\sum _{\mathcal{\ell }=1}^j{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}+\left(-{\displaystyle \sum _{\mathcal{\ell }=1}^i{w}_{\mathcal{\ell }}}\right)\frac{w_{j+1}}{\sum _{\mathcal{\ell }=1}^j{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}=0,$$

for $i<j\le N-1$, and

$${(QR)}_{ij}={\displaystyle \sum _{k=1}^N{q}_{ik}}{r}_{kj}={\displaystyle \sum _{k=1}^j{w}_k}\frac{w_{j+1}}{\sum _{\mathcal{\ell }=1}^j{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}+w_{j+1}\frac{-1}{\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}=0,$$

for $j<i\le N-1$. In particular, this proves that $R=Q^{-1}$.

Finally, we verify Definition 1, condition (4) by direct computations. First, it is easily verified that we have $Q\text{diag}(\mathit{x})=S≔{(s_{i,j})}_{i,j=1,\dots ,N}$ with

$$s_{ij}=q_{ij}{x}_j=w_j{x}_j,\hspace{1em}\hspace{1em}j\le i,\hspace{1em}\hspace{1em}{s}_{i,i+1}=-x_{i+1}{\displaystyle \sum _{\mathcal{\ell }=1}^i{w}_{\mathcal{\ell }}},\hspace{1em}\hspace{1em}i=1,\dots ,N-1.$$

Now, let us compute $Q\text{diag}(\mathit{x})Q^{-1}=T≔{(t_{i,j})}_{i,j=1,\dots ,N}.$

We have

$$t_{i,N}={\displaystyle \sum _{k=1}^N{s}_{i,k}}{r}_{k,N}=\frac{1}{\overline{w}}\left({\displaystyle \sum _{k=1}^i{w}_k}{x}_k-x_{i+1}{\displaystyle \sum _{\mathcal{\ell }=1}^i{w}_{\mathcal{\ell }}}\right)\le 0,$$

for $i=1,\dots ,N-1$, because the $x_i$ are ordered increasingly. Similarly,

$$t_{N,j}={\displaystyle \sum _{k=1}^N{s}_{N,k}}{r}_{k,j}={\displaystyle \sum _{k=1}^j{x}_k}{w}_k\frac{w_{j+1}}{\sum _{\mathcal{\ell }=1}^j{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}+x_{j+1}{w}_{j+1}\frac{-1}{\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}\le 0,$$

for $j=1,\dots ,N-1$. Next,

$$t_{i,j}={\displaystyle \sum _{k=1}^N{s}_{i,k}}{r}_{k,j}=\left({\displaystyle \sum _{k=1}^i{w}_k}{x}_k-x_{i+1}{\displaystyle \sum _{\mathcal{\ell }=1}^i{w}_{\mathcal{\ell }}}\right)\frac{w_{j+1}}{\sum _{\mathcal{\ell }=1}^j{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}\le 0,$$

for $i<j\le N-1$. Finally,

$$t_{i,j}={\displaystyle \sum _{k=1}^N{s}_{i,k}}{r}_{k,j}={\displaystyle \sum _{k=1}^j{w}_k}{x}_k\frac{w_{j+1}}{\sum _{\mathcal{\ell }=1}^j{w}_{\mathcal{\ell }}\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}+w_{j+1}{x}_{j+1}\frac{-1}{\sum _{\mathcal{\ell }=1}^{j+1}{w}_{\mathcal{\ell }}}\le 0,$$

for $j<i\le N-1$. This verifies Definition 1, condition (4) and proves the theorem. □

2.5. Proof of Theorem 2

Fix an admissible matrix $Q\in \mathcal{Q}$. The main idea of the proof is to reduce the study to the process $\mathit{Z}=Q\mathit{Y}$ and prove that its associated SDE admits an ${\mathbb{R}}_{+}^N$-valued solution.

We start by writing the SDE for $\mathit{Z}$. For this we first observe that due to (2), $\mathit{Y}$ satisfies

where $b_{\mu }(z)=b(z)+\mu$. Using Q as a transformation of basis (in the sense $\mathit{Z}=Q\mathit{Y}$), we get, thanks to the invertibility of Q, the following SDE

Using the Definition 1 admissibility conditions (2) and (3), we have that ${\mathit{w}}^{\top }{Q}^{-1}={\mathit{e}}_N^{\top }Q{Q}^{-1}={\mathit{e}}_N^{\top }$ and $Q\mathbf{1}=\overline{w}{\mathit{e}}_N$, which simplifies the equation to

Recall that $Z^{(N)}$ is the N-th component of $\mathit{Z}$.

In order to prove Theorem 2, it suffices to prove that for each ${\mathit{Z}}_0\in {\mathbb{R}}_{+}^N$, there exists an ${\mathbb{R}}_{+}^N$-valued $\mathit{Z}$ weak solution to (3). In particular, this would hold for any initial value of the form ${\mathit{Z}}_0=Q{\mathit{Y}}_0$ with $Y_0\in \mathcal{D}=Q^{-1}{\mathbb{R}}_N^{+}$, and setting $\mathit{Y}=QZ$, one obtains a $\mathcal{D}$-valued weak solution $\mathit{Y}$ to (6) started at ${\mathit{Y}}_0$.

Hence, this boils down to establish that the set ${\mathbb{R}}_{+}^N$ is stochastically viable with respect to Equation (3). Viability and invariance theory for stochastic differential equations have been extensively studied in the literature in various contexts and with different assumptions on the domain and the coefficients; we refer to Abi Jaber et al. (2019a), Da Prato and Frankowska (2004, 2007), and the references therein.

For the nonnegative orthant ${\mathbb{R}}_{+}^N$, the characterization in terms of the coefficients is very simple and means that, at boundary points, the diffusive coefficient has to be tangential to the boundary and the drift inward pointing. This is summarized in the following lemma.

Lemma 2.

Let $\tilde{b},\tilde{\sigma }:{\mathbb{R}}^N\to {\mathbb{R}}^N$ be continuous satisfying the growth conditions

$$\Vert \tilde{b}(\mathit{z})\Vert +\Vert \tilde{\sigma }(\mathit{z})\Vert \le L(1+\Vert \mathit{y}\Vert ),\hspace{1em}\hspace{1em}\mathit{z}\in {\mathbb{R}}^N,$$

and the boundary conditions, for all $\mathit{z}\in {\mathbb{R}}_{+}^N$,

$$z_i=0\hspace{0.17em}\Rightarrow \hspace{0.17em}{\mathit{e}}_i^{\top }\tilde{b}(\mathit{z})\ge 0 \text{and} {\mathit{e}}_i^{\top }\tilde{\sigma }(\mathit{z})=0,\hspace{1em}i=1,\dots ,N,$$(4)

then, for each ${\tilde{\mathit{Z}}}_0\in {\mathbb{R}}_{+}^N$, there exists a weak ${\mathbb{R}}_{+}^N$-valued solution to the following SDE

$$\hspace{0.17em}\text{d}{\tilde{\mathit{Z}}}_t=\tilde{b}\left({\tilde{\mathit{Z}}}_t\right)\hspace{0.17em}\text{d}t+\tilde{\sigma }\left({\tilde{\mathit{Z}}}_t\right)\hspace{0.17em}\text{d}{W}_t.$$

Proof.

See for instance Da Prato and Frankowska (2007, example 2.7). □

We now proceed to the proof of Theorem 2.

Proof of Theorem 2.

It remains to apply Lemma 2 on Equation (3). For this, we define

$$\tilde{b}(\mathit{z})=-Q\text{diag}(\mathit{x})Q^{-1}\mathit{z}+\overline{w}{b}_{\mu }(z_N){\mathit{e}}_N\hspace{1em}\text{and}\hspace{1em}\tilde{\sigma }(\mathit{z})=\overline{w}\sigma (z_N){\mathit{e}}_N,\hspace{1em}z\in {\mathbb{R}}^N.$$

Then, it readily follows from the continuity and growth conditions of b and $\sigma$ that $\tilde{b},\tilde{\sigma }$ are also continuous with at most linear growth conditions. As for the Boundary Conditions (4), we fix $\mathit{z}\in {\mathbb{R}}_{+}^N$ such that $z_i=0$ for some $i=1,\dots ,N$.

• For the diffusion term, we have

  $${\mathit{e}}_i^{\top }\tilde{\sigma }(\mathit{z})=\overline{w}\sigma (z_N){\mathit{e}}_i^{\top }{\mathit{e}}_N=0,$$
  because ${\mathit{e}}_i^{\top }{\mathit{e}}_N=0$ if $i<N$ and $\sigma (z_N)=\sigma (0)=0$ if $i=N$, where we used the boundary condition on $\sigma$ in (3).

• For the drift term, we first observe that for the same reason $b_{\mu }(z_N){\mathit{e}}_i^{\top }{\mathit{e}}_N=(b(z_N)+\mu ){\mathit{e}}_i^{\top }{\mathit{e}}_N\ge 0$, because $b(0)+\mu \ge 0$ thanks to the boundary condition on b in (1) and the fact that $\mu \ge 0$, so that we can write

  $$\begin{array}{ll}{\mathit{e}}_i^{\top }\tilde{b}(\mathit{z})\hfill & =-{\mathit{e}}_i^{\top }Q\text{diag}(\mathit{x})Q^{-1}\mathit{z}+\overline{w}{b}_{\mu }(z_N){\mathit{e}}_i^{\top }{\mathit{e}}_N\hfill \\ \hfill & \ge -{\displaystyle \sum _{j\ne i}(Q\text{diag}(}\mathit{x})Q^{-1}{)}_{ij}{z}_j\hfill \\ \hfill & \ge 0,\hfill \end{array}$$

where the first inequality follows from $z_i=0$ and the second inequality follows from the Definition 1, admissibility condition (4) in for the matrix Q and the fact that $z_j\ge 0$.

This shows that the Boundary Conditions (4) are satisfied by $\tilde{b},\tilde{\sigma }$, so that an application of Lemma 2 yields the existence of an ${\mathbb{R}}_{+}^N$-valued solution $\mathit{Z}$ to (3) for any initial condition ${\mathit{Z}}_0\in {\mathbb{R}}_{+}^N$. In particular, it holds for the initial value ${\mathit{Z}}_0=Q{\mathit{Y}}_0$ with $Y_0\in \mathcal{D}=Q^{-1}{\mathbb{R}}_{+}^N$. Setting $\mathit{Y}=QZ$, one obtains a $\mathcal{D}$-valued weak solution $\mathit{Y}$ to (2) started at ${\mathit{Y}}_0$ and completes the proof of theorem. □

3. The Weak Scheme Is Cone-Preserving

In Section 2, we determined the state space $\mathcal{D}\subseteq {\mathbb{R}}^N$ of the multifactor square-root process ${\mathit{V}}^N$ given by (2). Assume now that we approximate the process ${\mathit{V}}^N$ using the weak simulation scheme proposed in Bayer and Breneis (2024). The goal of this section is to prove that the resulting approximation has the same viable domain $\mathcal{D}$ as ${\mathit{V}}^N$.

Let us first start by recalling the weak simulation scheme of Bayer and Breneis (2024). First, the SDE in (2) is split into two parts, one containing the drift and the other the diffusion. Denote by $D(\mathit{z},h)≔{\mathit{Z}}_h≔{(Z_h^i)}_{i=1}^N$ the solution at time h of the ordinary differential equation (ODE)

and by $S(\mathit{y},h)≔{\mathit{Y}}_h≔{(Y_h^i)}_{i=1}^N$ the solution at time h of the SDE

Then, the ODE (3) is linear and can hence be solved exactly. Therefore, the simulation scheme $\widehat{D}$ for the ODE is simply given by

where

We now recall the simulation scheme for the SDE (3). Note that the right-hand side of (3) is the same for all i. Thus, after multiplying (3) with $\mathit{w}$, we get

where $\overline{w}≔{\mathbf{1}}^{\top }\mathit{w}.$ This is now a one-dimensional SDE, which was already studied in Lileika and Mackevičius (2021), where a second-order simulation scheme was given. This scheme is based on matching the first six centralized moments (up to errors of order $O(h^3)$), while preserving the nonnegativity of Y. Define the quantities

Then, we define ${\widehat{Y}}_h$ to be the random variable which is $x_i$ with probability $p_i$, $i=1,2,3$, noting, in particular, that $p_1+p_2+p_3=1$.

We can now reconstruct an approximation $\widehat{\mathit{Y}}$ from $\widehat{Y}$. Indeed, because the right-hand side of (3) is the same for all $i=1,\dots ,N$, the solution of (3) must be of the form

for some scalar random variable R. Taking the inner product of (3) with $\mathit{w}$, we get

Hence, we set

Finally, we use Strang splitting to get the scheme

for approximating ${\mathit{V}}_h$ given $\mathit{v}$. Therefore, we get a simulation algorithmwhere $0=t_0<t_1<\cdots <t_M=T.$ The only problem that could occur is that the square root in (6) or (7) is not well-defined. However, note that if we can prove that ${\mathit{V}}^{N,M}$ does not leave $\mathcal{D}$, where $\mathcal{D}$ is the same cone as in Theorem 2, then in particular, $x={\mathit{w}}^{\top }\mathit{y}$ in (5) will always be nonnegative, and, hence, the square roots in (6) and (7) are always well-defined. Proving that ${\mathit{V}}^{N,M}$ stays in $\mathcal{D}$ is the aim of the following theorem.

Theorem 3.

Let Q be an admissible matrix, and let ${\mathit{v}}_0$ be chosen according to (2). Then, for all $\mathit{v}\in Q^{-1}{\mathbb{R}}_{+}^N$ and $h\ge 0$, the weak simulation algorithm $A^{\text{CIR}}$ described above satisfies $A^{\text{CIR}}(\mathit{v},h)\in Q^{-1}{\mathbb{R}}_{+}^N$. In particular, $A^{\text{CIR}}$ is well-defined.

Proof.

Given $\mathit{z},\mathit{y}\in Q^{-1}{\mathbb{R}}_{+}^N$ and $h\ge 0$, we want to show that $D(\mathit{z},h)\in Q^{-1}{\mathbb{R}}_{+}^N$ and $\widehat{S}(\mathit{y},h)\in Q^{-1}{\mathbb{R}}_{+}^N$. This will prove the theorem.

Consider first the algorithm $\widehat{S}$. Recall that $\widehat{S}(\mathit{y},h)=\mathit{y}+R\mathbf{1}$ for some scalar random variable R. We have to verify that $Q\widehat{S}(\mathit{y},h)=Q\mathit{y}+RQ\mathbf{1}\in {\mathbb{R}}_{+}^N$. The last component of this vector is given by

$${(Q\widehat{S}(\mathit{y},h))}_N={\mathit{w}}^{\top }\mathit{y}+R\overline{w},$$

and we recall that this was given by the random variable ${\widehat{Y}}_h$ in Section 3, which by definition is nonnegative, as verified in Lileika and Mackevičius (2021). Conversely, for $i=1,\dots ,N-1$, we have

$${(Q\widehat{S}(\mathit{y},h))}_i={(Q\mathit{y})}_i+0\ge 0,$$

by the assumption that $Q\mathit{y}\in {\mathbb{R}}_N^{+}$. Hence, $\widehat{S}$ leaves the domain $Q^{-1}{\mathbb{R}}_{+}^N$ invariant.

Next, consider the algorithm D. Recall that $D(\mathit{z},h)$ was given as the exact solution at time h of the ODE

$$\text{d}{\mathit{Z}}_t=-\text{diag}(\mathit{x})({\mathit{Z}}_t-{\mathit{v}}_0)\hspace{0.17em}\text{d}t+(\theta -\lambda {\mathit{w}}^{\top }{\mathit{Z}}_t)\mathbf{1}\hspace{0.17em}\text{d}t,\hspace{1em}\hspace{1em}{\mathit{Z}}_0=\mathit{z}.$$

Note that due to (2), $\text{diag}(\mathit{x}){\mathit{v}}_0=\mu \mathbf{1}$ for some $\mu \ge 0$. Defining $\tilde{\mathit{Z}}≔Q\mathit{Z}$, we have

$$\text{d}{\tilde{\mathit{Z}}}_t=(-Q\text{diag}(\mathit{x})Q^{-1}\widetilde{{\mathit{Z}}_t}+(\theta +\mu -\lambda {\tilde{Z}}_t^N)Q\mathbf{1})\text{d}t,\hspace{1em}\hspace{1em}{\tilde{\mathit{Z}}}_0=Q\mathit{z}\in {\mathbb{R}}_{+}^N,$$

and we have to show that ${\tilde{\mathit{Z}}}_t\in {\mathbb{R}}_{+}^N$. We prove this by invoking Lemma 2, where we note that

$$b(\mathit{z})=-Q\text{diag}(\mathit{x})Q^{-1}\mathit{z}+(\theta +\mu -\lambda z^N)Q\mathbf{1},\hspace{1em}\hspace{1em}\sigma \equiv 0.$$

In particular, we have to show that $b_i(\mathit{z})\ge 0$ for $\mathit{z}\in {\mathbb{R}}_{+}^N$ with $z^i=0$.

We start with $i=N$. Here, we have

$$b_N(\mathit{z})=-{(Q\text{diag}(\mathit{x})Q^{-1}\mathit{z})}_N+(\theta +\mu )\overline{w}.$$

Of course, $(\theta +\mu )\overline{w}\ge 0$. Moreover, because of Definition 1, $-{(Q\text{diag}(\mathit{x})Q^{-1}\mathit{z})}_N$ is a linear combination of $z^i$ where all the coefficients are nonnegative, with the exception of the coefficient of $z^N$. However, because $z^N=0$, this implies that $b_N(\mathit{z})\ge 0$.

Next, consider $i=1,\dots ,N-1$. Then,

$$b_i(\mathit{z})=-{(Q\text{diag}(\mathit{x})Q^{-1}\mathit{z})}_i.$$

As before, $-{(Q\text{diag}(\mathit{x})Q^{-1}\mathit{z})}_i$ is again a linear combination of the $z^j$, where all coefficients are nonnegative, with the exception of the coefficient of $z^i$. But because $z^i=0$, this implies that $b_i(\mathit{z})\ge 0$, proving the theorem. □

4. Solving PDEs

As an application of the domain, we want to solve PDEs. Recall that the multifactor square-root process ${\mathit{V}}^N$ is given by

see (2). After a transformation of variables using $\mathit{Z}≔Q{\mathit{V}}^N$ and ${\mathit{z}}_0≔Q{\mathit{v}}_0$, where Q is the matrix in Theorem 1, we have

Let $f:{\mathbb{R}}_{+}^N\to \mathbb{R}$ be a “nice” payoff function. Then, we define the value function $u:{\mathbb{R}}_{+}^N\times [0,T]\to \mathbb{R}$,

Then, u satisfies the PDE

with the boundary condition $u(\mathit{z},T)=f(\mathit{z})$, $\mathit{z}\in {\mathbb{R}}_{+}^N$.

For numerical approximation, we then need to truncate the domain in space and impose appropriate boundary conditions. For simplicity, we will instead fabricate an appropriate source term such that the PDE has an explicit, given solution, which we then also impose as Dirichlet boundary condition on the boundary of the truncated domain.

Specifically, suppose that we want the exact solution to have the form

Plugging this formula into (4), we obtain a source term

with $g_{ij}={(Q\text{diag}(\mathit{x})Q^{-1})}_{ij}$; that is, u satisfiesnow with the terminal condition $u(\mathbf{z},T)=\tilde{u}(\mathbf{z},T)$. The precise parameters chosen are summarized in Table 1. In dimension $N=2$, we choose the admissible matrix Q given by Example 3. We furthermore choose $\mathbf{\alpha }=(3,\hspace{0.17em}4)$, $\beta =1.6$, and the terminal time $T=2$.

Table 1. Parameters of the Stochastic Volatility of the Lifted Rough Heston Model Used for the Numerical Example

After truncation of the domain, we solve the PDE by the finite element method, using the package FEniCSx; see Baratta et al. (2023), and compare against the exact solution $\tilde{u}$. In Table 2, we present the $L^2$-errors on the truncated domain for three choices of truncated domains, each of side-length four:

Table 2. $L^2$ Errors over the Truncated Domain for the Approximate Finite Element Method (FEM) Solution to the PDE for $n_t=n_x=n$ in Dimension $N=2$

1. $D={[0,4]}^2$, corresponding to a truncation in v-space which respects the cone-shaped actual domain of the process;

2. $D={[-0.5,3.5]}^2$ corresponding to a truncation in v-space, which respects neither the cone-shaped actual domain nor the nonnegativity condition;

3. $D=[-0.5,3.5]\times [0,4]$ corresponding to a domain truncation, which does not respect the cone-shaped actual domain in v-space but does respect the nonnegativity.

We use first-order Lagrange-type finite elements, with $n_t$ time-steps as well as mesh-size $n_x=n_t$ in each space dimension. (We refer to https://github.com/bayerc2/domain_multifactor_volterra for more details.)

We would like to emphasize that, although it might seem trivial to choose ${[0,4]}^2$ as the domain instead of, say, $[-0.5,3.5]\times [0,4]$, this choice is based on correctly identifying the matrix Q and the domain $\mathcal{D}=Q^{-1}{\mathbb{R}}_{+}^2$, which is appropriately truncated here to $Q^{-1}{[0,4]}^2$. Without knowledge of Q, one would need to guess $\mathcal{D}$ to truncate the domain, a task that becomes increasingly nontrivial in higher dimensions.

When nonnegativity of the variance process is preserved (cases 1 and 3), the numerical method empirically exhibits second-order convergence, with slightly smaller error when the computational domain is a subset of the invariant domain of the process (case 1). On the other hand, when nonnegativity of the variance process is not preserved on the computational domain (case 2), the error explodes due to the instability of the heat equation backward in time.

We also provide an example in dimension $N=3$. In this case, we follow the construction outlined in the appendix. Note that the construction is not unique, so we tested different solutions to the system of inequalities leading to different admissible Q. Specifically, we choose

• $a=1$, $b=2$, the default choice suggested in the appendix.

• $a=2.04$, $b=0.88$, an admissible choice obtained by minimizing $\Vert Q\text{diag}(\mathbf{x})Q^{-1}\Vert$—related to the Lipschitz constant of the transformed drift—over all admissible choices of parameters a, b.

• $a=0.51$, $b=-0.03$, obtained by maximizing $\Vert Q\text{diag}(\mathbf{x})Q^{-1}\Vert$.

It turns out that there was no significant difference in the numerical results based on the different choices of transformation. Hence, we only report the results for the first choice. We report our numerical results in Table 3. Note that we had to limit the grid sizes due to the severely increased computational time. Hence, the accuracies reported may seem disappointing. However, note that the $L^2$-error is unnormalized here. A normalized error would be obtained by dividing by the total volume of the (computational) domain, which is $4^3=64$ in this case. We observe an expected error decay when the domain is respected and highly erratic, diverging behavior when positivity is not preserved.

Table 3. $L^2$ Errors over the Truncated Domain for the Approximate FEM Solution to the PDE for $n_t=n_x=n$ in Dimension $N=3$

5. Remark on the Link Between the Sets $\mathcal{E}$ and $\mathcal{G}$

In this section, we argue that the two abstract “invariance” sets that appeared in the literature in Abi Jaber and El Euch (2019a) and Cuchiero and Teichmann (2020) are equal. This part is valid for more general locally square-integrable kernels K beyond the weighted sum of exponential case.

We introduce the following notations. For suitable functions f, g and measure L, we denote their convolution by *:

The shift operator ${\mathrm{\Delta }}_h$ with $h\ge 0$ maps any function f on ${\mathbb{R}}_{+}$ to the function ${\mathrm{\Delta }}_{h}f$ given by

If the function f on ${\mathbb{R}}_{+}$ is right-continuous and of locally bounded variation, the measure induced by its distributional derivative is denoted $\text{d}f$, so that $f(t)=f(0)+{\int }_{[0,t]}\text{d}f(s)$ for all $t\ge 0$. By convention, df does not charge $\{0\}$.

Two sets appeared so far in the literature to characterize the nonnegativity of solutions to stochastic Volterra equations:

1. Set of Cuchiero and Teichmann (2020, equation (4.7), definition 4.12, and theorem 4.17(i)):¹

   $$\mathcal{E}={\displaystyle \underset{\eta >0}{\cap }{\mathcal{E}}^{\eta }} \text{with} {\mathcal{E}}^{\eta }≔\left\{g_0:[0,T]\to \mathbb{R} \text{such} \text{that} g_0-R^{\eta }*g_0\ge 0\right\}.$$
   Here, $R^{\eta }(t)$ is the resolvent of the second kind of the kernel $(\eta K)$ defined by
   $$R^{\eta }=\eta K-\eta K*R^{\eta }=\eta K-R^{\eta }*\eta K.$$

2. Set of Abi Jaber and El Euch (2019a, equations (2.4)–(2.5) and theorem 2.1):

   $$\begin{array}{ll}\mathcal{G}=\{g_0:\hfill & [0,T]\to \mathbb{R} \text{such} \text{that} \hfill \\ \hfill & {\mathrm{\Delta }}_h{g}_0-({\mathrm{\Delta }}_{h}K*L)(0)g_0-d({\mathrm{\Delta }}_{h}K*L)*g_0\ge 0 \text{and} g_0(0)\ge 0.\}\hfill \end{array},$$

where $L(\text{d}t)$ is the resolvent of the first kind of the kernel²

One can argue that the two sets are equal:

because both conditions that appear in the set are necessary and sufficient conditions for the nonnegativity of the linear Volterra equationfor $\eta >0$. Indeed, on the one hand, the solution of (8) can be expressed in terms of the resolvent of the second kind in the formwhich is exactly the form that appears in $\mathcal{E}$. On the other hand, by relying on the properties of the resolvent of the first kind—see, for instance, Abi Jaber and El Euch (2019a, the proof of theorem A.2)—one can write that

We note that the first line is exactly the condition that appears in the set $\mathcal{G}$. We sketch why $g_0\in \mathcal{G}$ is equivalent to the nonnegativity of $f^{\eta }$ for all $\eta >0$, under suitable assumptions on the nonnegative kernel K (for instance, when K is a weighted sum of exponentials; see Abi Jaber and El Euch 2019a, $(H_0)$ and $(H_1)$ for precise conditions). The implication $\Rightarrow$ follows from Abi Jaber and El Euch (2019a, theorem A.2) with $b(x)=-\eta x$ and $\sigma (x)=0$ therein. For the converse direction, fix $\eta >0$ and assume that $f^{\eta }\ge 0$. Clearly, $g_0(0)=f^{\eta }(0)\ge 0$, and from (5), we have

To conclude that $g_0\in \mathcal{G}$, it suffices to note that the right-hand side tends to zero as $\eta \to \infty$, because $f^{\eta }\to 0$ in this limit (this follows from the Laplace transform $\widehat{f^{\eta }}$ of (8), given by $\widehat{f^{\eta }}=\frac{\widehat{g_0}}{1+\eta \widehat{K}}$, which goes to 0 as $\eta \to \infty$).

In conclusion, the nonnegativity of the Linear Volterra Equation (8), for any $\eta >0$, is equivalent to the condition that appears in $\mathcal{E}$, as well as the one that appears in $\mathcal{G}$, which shows that the two sets are equal.

In principle, to establish a link with our cone $\mathcal{D}$, one should restrict to kernels that are weighted sums of exponentials of the form

and input curves of the form

Then, the resolvents of the second kind and first kind for such kernels must be computed and plugged into the conditions defining the sets $\mathcal{E}$ and $\mathcal{G}$. Even in dimension $N=2$, this leads to highly cumbersome and nontrivial computations, and it is not clear how to explicitly determine a suitable domain, as, for instance, our cone $\mathcal{D}$, for $Y_0^i$ from $\mathcal{E}$ and $\mathcal{G}$. This makes the approach in the current paper particularly crucial.