CARMA Approximations and Estimation

CARMA(p, q) processes are compactly defined through a stochastic differential equation (SDE) involving q + 1 derivatives of the Lévy process driving the noise, despite this latter having in general no differentiable paths. We replace the Lévy noise with a continuously differentiable process obtained via stochastic convolution. The solution of the new SDE is then a stochastic process that converges to the original CARMA in an L2 sense. CARMA processes are largely applied in finance, and this article mathematically justifies their representation via SDE. Consequently, the use of numerical schemes in simulation and estimation by approximating derivatives with discrete increments is also justified for every time step. This provides a link between discrete-time ARMA models and CARMA models in continuous-time. We then analyze an Euler method and find a rate of convergence proportional to the square of the time step for discretization and to the inverse of the convolution parameter. These must then be adjusted to get accurate approximations.


INTRODUCTION
Continuous-time autoregressive moving-average processes, CARMA in short, represent the continuous-time version of the well-known ARMA models. Because of their linear structure, they are convenient for modeling empirical data, hence they find application in several fields, particularly in finance. Indeed, with a Lévy process as a driver for the noise, one obtains a rich class of possibly heavy-tailed continuous-time stationary processes, and such heavy tails are frequently observed in finance [1]. In Barndorff-Nielsen and Shepard [2], for example, a Lévydriven Ornstein-Uhlenbeck process (corresponding to a CAR(1)) is employed in a stochastic volatility model, and [3] extend the idea to a wider class of Lévy-driven CARMA processes with non-negative kernel, allowing for non-negative, heavy-tailed processes with a larger range of autocovariance functions.
In view of financial applications, Brockwell et al. [4] derive efficient estimates for the parameters of a non-negative Lévy-driven CAR(1) process, and Brockwell et al. [5] generalize the technique to higher-order CARMA processes with a non-negative kernel. In Benth et al. [6] and Benth and Benth [7], the authors exploit the link between discrete-time ARMA and continuous-time ARMA to model the observed processes in a continuous manner, starting from discrete-time observations-based estimates, in the context of energy markets and weather derivatives. It is also worth mentioning that CARMA processes allow us to construct a wide range of long-memory processes, as shown, e.g., in Marquardt [8] and Marquardt and James [9].
Among the several financial applications, in Paschke and Prokopczuk [10,11] a CARMA process is employed for the long-term equilibrium component in the two-factors spot price model introduced by Schwartz and Smith [12] (which uses an Ornstein-Uhlenbeck process), while Benth and Koekebakker [13] use a CARMA process for the shortterm factor instead. In Benth et al. [14] and García et al. [15], a CARMA process models the spot price in the electricity market, while in Benth et al. [6] and Benth and Benth [7] it models weather variables, such as temperature and wind speed. Finally, Andresen et al. [16] apply CARMA processes to interest rates modeling and [17] to the dynamics of regional ocean freight rates.
Formally, CARMA(p, q) processes are compactly defined through a stochastic differential equation (SDE) involving q + 1 derivatives of the Lévy process driving the noise, despite the latter having in general no differentiable paths. Numerical schemes, such as the Euler scheme, are also applied both for simulations and estimation of CARMA processes. These schemes approximate derivatives with discrete increments with respect to a time step , and their convergence for approaching zero is indeed a well-studied topic. However, this opens up for a discrepancy since it is not mathematically clear what happens when the time step approaches zero, as, in general, the derivatives of a Lévy process are not well-defined.
We give a mathematical justification for the representation via SDE of a CARMA process, and, consequently, for the discrete approximations used for simulating and estimating the continuous-time dynamics of CARMA processes. This is done by replacing the Lèvy process L driving the noise with a process smooth enough so that the SDE makes mathematical sense. For ε > 0 we construct a continuously differentiable function L ε , which converges in the L 2 sense to the original Lévy process in the space [0, T] × , for T < +∞. More precisely, the process L ε is constructed via stochastic convolution of L with a smooth sequence of functions {φ ε } ε>0 depending on the convolution parameter ε > 0 and with first derivative converging to the delta-Dirac function when ε approaches zero. We show a rate of convergence of order 1.
For the L ε continuous differentiable, the framework developed solves the differentiability issue for CAR(p) processes in which no more than the first derivative of the driving noise process is required. We can, however, make use of the result in Brockwell et al. [5,Proposition 2]. This states that, under some conditions on the coefficients of the SDE, a CARMA(p, q) process can be decomposed into a sum of p dependent CAR(1) processes. The approach can then be extended to CARMA(p, q) processes of every order q ≥ 0 by considering p different approximations based on the continuous differentiable L ε and taking their sum. This justifies an L ε -only continuous differentiable, as the general case q ≥ 0 is still covered.
We thus replace the Lèvy process L with L ε in the SDE defining the CARMA process X. This gives a different SDE, which depends on ε, whose solution is the new process X ε . We show that X ε converges to X in L 2 ([0, T] × ) and that the convergence is of order 2.
The ε-approximation of the Lévy process justifies the use of an Euler scheme for every time step . We present a possible approach for simulating L ε starting from a simulated path of the Lèvy process L and then X ε starting from L ε . This, for a fixed time interval , gives a new noise process L ε and a third process X ε . We show a rate of convergence for both L ε and X ε proportional to 2 and inversely so to ε. This highlights that the two parameters and ε must be chosen in a way to get accurate simulations and prevent them from exploding.
We finally consider a possible Euler scheme approach for estimating the CARMA dynamics starting from the simulated time series. It is worth mentioning that the process L ε is not a Lèvy process. Due to the convolution operator, it indeed fails on the crucial property of independent increments. We show in particular that the increments depend on the ratio ε , and, as this ratio becomes smaller, the increments become more dependent on each other. This has a negative drawback when it comes to estimation. We find indeed that the Euler scheme recovers the true dynamics of the CARMA process only if ε is smaller than .
All the results are deeply investigated and are also confirmed by numerical experiments. In particular, we consider two different Lèvy processes, a Brownian motion and a normal inverse Gaussian (NIG) process, and CARMA processes of two different orders, a Gaussian CAR(1) and a Gaussian CARMA (3,2). Crucial to all the numerical simulations is the convolution property of both the Gaussian and the NIG distribution function, which allows us to compare path-wise different time grids by simulating only one path of the Lèvy process.
Since CARMA processes have significant application within finance, from the modeling of the electricity spot price to the modeling of wind speed and temperature, with the present study, we justify its mathematical formulation and application in a rigorous way. Moreover, this may open for further applications in the context of sensitivity analysis with respect to paths of the driving noise process. Indeed, since the new process L ε admits derivatives, one could potentially analyze the sensitivity of the CARMA process, and hence of the variable of interest (e.g., the spot price), with respect to the Lévy noise in a path-wise sense. It might also be that when observing a signal, the driving-noise process is more regular than, e.g., a Brownian motion. This might then be modeled by the regular process here introduced and studied.
The rest of the paper is organized as follows. In section 2 we construct the approximating process L ε by stochastic convolution, and we study the rate of convergence. In section 3 we introduce CARMA processes and their spectral representation. We then show the convergence of the new process X ε to the original CARMA process, and we briefly recover the general case q > 0. In section 4 we study an Euler scheme to simulate the process X ε starting from the simulations of L ε , and then to calibrate the parameters of the CARMA dynamics starting from the discrete time series. Appendix A in Supplementary Material contains the proofs of the main results.

CONSTRUCTION AND CONVERGENCE RESULTS
Let ( , F , P) be a probability space. We consider a onedimensional Lèvy process L, which we assume to be a squareintegrable martingale with characteristic triplet (κ, , ν), ν being the Lévy measure. From Cont and Tankov [18,Proposition 3.13], for every t > 0 the variance of L(t) is of the form Var L(t) = tσ 2 with We now consider a sequence of measurable functions {φ ε } ε>0 :[0, T] −→ R such that, for every ε > 0, φ ε is twice continuously differentiable on the interval [0, T] and φ ε (0) = 0. We then construct the process {L ε (t), t ≥ 0} by stochastic convolution of the Lévy process L with φ ε : L ε is well-defined in [0, T] for every T ∈ ]0, +∞[, and, in particular, the condition φ ε (0) = 0 is used to prove the following result, where the symbol * denotes the convolution operator. For ψ ε (x) : = φ ′ ε (x), since φ ε (0) = 0, we write that φ ε (t) = t 0 ψ ε (u)du so that φ ε is twice continuously differentiable for any continuously differentiable ψ ε . We introduce the function ψ(x) : = 2 √ 2π exp − 1 2 x 2 , and for every ε > 0 we define Let us notice that R ψ(x)dx = R ψ ε (x)dx = 2, while the total area of the mollifier function usually equals 1. However, we construct the process {L ε (t), t ≥ 0} by stochastic convolution restricted to the interval [0, t], and, denoting by the standard Gaussian cumulative distribution function, the area of ψ ε restricted to [0, t] is so that the unitary total integral property for ψ ε is preserved in the limit. By differentiating Equation (2.4) with respect to t, according to the Leibniz rule, we get This shows that L ε ∈ C 1 ([0, T]), though L ε / ∈ C n ([0, T]) for any n > 1, because of the term L(t) which is not differentiable, unless φ ′ ε (0) = ψ ε (0) = 0, which is not our case. Hence, we start with the Lévy process L, which is not differentiable, and we end up with L ε , which is continuously differentiable. For a CARMA(p, q) process, we need ε being the m-th derivative of φ ε . However, for our particular choice of the convolution function ψ ε , these conditions of null derivatives in zero are not satisfied. Another possible choice would have been, for example, to take a sequence of gamma densities with scale parameter θ ε converging to 0 as ε converges to 0 and shape With L ε constructed above, we thus address the differentiability issue in the stochastic differential equation defining the CARMA process only for q = 0, namely when the process is a CAR(p). However, Brockwell et al. [5,Proposition 2] shows that, under some conditions, we can express a CARMA(p, q) process as the sum of p dependent and possibly complex-valued CAR(1) processes. The idea is then to construct the ε-approximation with L ε in Equation (2.3) for each of the p processes of CAR(1)-type and obtain the ε-CARMA(p, q) process by summing them up. We can thus deal with q > 0, but it must be done separately. With the gamma densities mentioned above, for example, the special treatment for q > 0 would not be necessary. However, we started this analysis from CAR(p) processes, and in this setting the scaled Gaussian densities ψ ε make the job. We shall analyze the case q > 0 in section 3.2.
A first important fact about L ε in Equation (2.3) is that the convolution operation does not preserve the properties of the Lévy process L.
is not a Lévy process.
Proof: See Appendix A.2 in Supplementary Material. By Proposition 2.2, for L = B, a Brownian motion, the approximated process B ε is no longer a Brownian motion. We can, however, show that the distribution of B is preserved.
has Gaussian distribution with zero mean and variance Proof: B ε is defined in Equation (2.3) as the stochastic integral of φ ε with respect to the Brownian motion B. Then B ε has Gaussian distribution with zero mean and variance t 0 φ 2 ε (s)ds, which can be explicitly calculated.
Let us point out that Lemma 2.3 is valid for every choice of the functions {φ ε } ε>0 , with the obvious modification for the variance formula t 0 φ 2 ε (s)ds. Moreover, from Lemma 2.3, we get that lim ε↓0 Var(B ε (t)) = t, as one would expect. Figure 1 shows indeed the typical bell-shape for a Gaussian distribution. We sampled the process B ε at time t = 1 for 5, 000 times and different values of ε. The empirical (black line) and the theoretical (red line) distributions (according to Lemma 2.3) are also shown. For ε = 1, there is a certain discrepancy between empirical and theoretical distributions, which may be due to the fact that for big values of ε the increments of B ε are more dependent from each other than for smaller values of ε. The same phenomena will be encountered in section 4.3.

Convergence of the Lévy Approximation
We prove now the convergence of the process L ε to the Lévy process L.  shows graphically the convergence stated in Theorem 2.4 for a Brownian motion and for a normal inverse Gaussian process with parameters (α, β, δ, µ) = (1, 0, true , 0), true being the time step used in the Euler discretization (see section 4). As ε decreases, the error between L ε and L decreases. Even with a large ε (ε = 1), the process L ε is able to capture the shape of L. It seems however to be slightly shifted to the right.
We compute explicitly the squared error between L ε and L at the terminal point t = T. Lemma 2.5. For every ε > 0 and T < +∞, the squared error between L ε and L is given by for σ 2 in Equation (2.1) and ϕ and respectively the density and the cumulative distribution functions of a standard Gaussian random variable.
Proof: The proof applies the integration by parts formula several times.
The following lemma concerning the cumulative Gaussian distribution function is useful. Lemma 2.6. The function can be expanded by integration by parts into the series !! denoting the double factorial. Moreover, for large x, the asymptotic expansion is also valid We now use the expansion (2.8) to calculate the rate of convergence of L ε to the Lévy process L, starting from Lemma 2.5. The expansion (2.9) will be used further in the paper.
Proposition 2.7. For every ε > 0 and T < +∞, there exists a constant K 1 > 0 such that Similarly to Lemma 2.3, the rate of convergence found in Proposition 2.7 can be calculated explicitly for every choice of the functions {φ ε } ε>0 , possibly obtaining a different order. The L 2 norm of the approximation error studied in Proposition 2.7 as a function of ε is in Figure 3A for the Brownian motion and in Figure 3B for the NIG Lévy process. In the next section we shall introduce CARMA processes in a more rigorous way and study their ε-approximation.

CARMA REPRESENTATION AND CONVERGENCE RESULTS
For p > q and t ≥ 0, a second-order Lévy-driven continuoustime ARMA process of order (p, q) is formally defined by the stochastic differential equation where P and Q are the characteristic polynomials of the CARMA process of the form with a j ≥ 0, j = 1, . . . , p, and a p > 0, b q = 1 and b j = 0 for q < j < p. Equation (3.1) is, however, usually interpreted by means of its state-space representation

3)
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org where ⊤ denotes the transpose operator, and By Itô's rule, the solution of the first equation in (3.3) at time t ≥ s is given by then the process {Y(t), t ≥ 0} is strictly stationary and the CARMA(p, q) process is given by The function g is referred to as the kernel function associated to the CARMA process. We also notice that the process X is a causal function of L, in the sense that X(t) is independent of {L(s) − L(t), s ≥ t} for all t ≥ 0. We shall assume that condition (3.4) is satisfied, referring to it as the causality condition or stationarity condition. In a similar way, we introduce the state-space representation for X ε . This is the process obtained by replacing the Lévy process L with the smoother L ε , namely . (3.7) The solution to Equation (3.7) can be expressed in terms of the kernel function g by We point out that the kernel function g in Equation (3.6) is continuously differentiable. Its value at x = 0 is given in the following lemma, where δ denotes the Kronecker delta.
Lemma 3.1. For p > q, the value of the kernel function at x = 0 is equal to g(0) = δ p,q+1 .
Proof: The proof follows by direct verification.
From the basic algebra analysis, if A has distinct eigenvalues, then the matrix exponential e Ax can be decomposed into the matrix product P ′ being the derivative of the polynomial P. We then find the following spectral representations for the kernel function g and its first derivative h(x) : = g ′ (x).  These representations will be used for proving the convergence results in the next section.

Convergence of the CARMA Approximation
We prove the convergence of X ε to X. We first start with the following preliminary result.  Here h(x) = g ′ (x) and h ε (x) : = g(0)ψ ε (x) +   (3.11) and prove the first convergence result. In Proposition 3.4, we do not require the limit to hold in x = 0. However, this does not affect the convergence in L 2 since X and X ε are both integral functions, and {x = 0} is a set of measure 0. Indeed, we can prove the required convergence for X ε .  Figure 4 shows the convergence proven in Theorem 3.5 for a Gaussian CAR(1) process with parameter a = 0.6. We notice that X ε captures the shape of X even with ε large, and we observe a right-shifting of X ε with respect to X, as noticed already for the Lévy process approximation.

From Lemma 3.2, we also derive the spectral representation
We now study explicitly the rate of convergence of X ε to the CARMA process X. Proposition 3.6. For every ε > 0 and T < +∞, there exists a constant K 2 > 0 such that  The rate of convergence obtained in Proposition 3.6 is shown in Figure 3C for a CAR(1). We point out that the results in this section have been stated and proved for a general CARMA(p, q) process. The convergence results are indeed valid for every q ≥ 0, even when L ε is not smooth enough to make the SDE having mathematical sense. We then present them for the general case. Since the same ideas apply to other families of convolution functions, such as the Gamma densities previously mentioned, this is useful for possible further applications.
We end this section with a distribution property holding for X ε when the original noise is driven by a Brownian motion. If the CARMA process X is driven by a Brownian motion B, then X has Gaussian distribution. We show that the same holds for X ε . This is indeed not a straightforward result since the process B ε is no longer a Brownian motion by means of Proposition 2.2. A rigorous proof is thus necessary. Proposition 3.7. For every ε > 0, the process {X ε (t), t ≥ 0} has Gaussian distribution.
Proof: See Appendix A.9 in Supplementary Material. Frontiers in Applied Mathematics and Statistics | www.frontiersin.org , 0 ≤ t ≤ 5} (red line) for different 's generated from a single path with ε = 0.01, for a Brownian motion and an NIG process with parameters (α, β, δ, µ) = (1, 0, true , 0). The black bar lines in the four rightmost images represent the difference between L ε and L ε .

When q > 0
The process L ε is continuously differentiable but is not twice continuously differentiable. However, at the right-hand side of the stochastic differential equation (3.1), we require q + 1 derivatives of the driving noise process. For the general case q > 0 we need the following result. Here, the X (r) are dependent and possibly complex-valued CAR(1) processes with α r : = Q(λ r ) P ′ (λ r ) .
stochastic differential equation (3.12) that will be used for numerical purposes in section 4. In Equation (3.12), we write ε in parentheses since it applies to both X (r) and X (r) ε . In Figure 5 we consider a Gaussian CARMA(3, 2) process with parameters (a 1 , a 2 , a 3 ) = (2.043, 1.339, 0.177) taken from [6] where a CAR(3) is used to model the daily average temperature dynamics in Stockholm (Sweden) and (b 0 , b 1 , b 2 ) = (2.0, 0.8, 1.0). In Figure 3D, we observe the rate of convergence as a function of ε.

NUMERICAL SCHEMES AND RESULTS
We consider a possible approach to simulate X and X ε for some CARMA parameters, and we then estimate these parameters starting from the simulated time series. In particular, we need first to simulate L, then L ε , and only afterwards can we simulate X ε . For n > 0, we shall consider the partition (n) : = {0 = s (n) 0 < s (n) 1 < · · · < s (n) n = T} such that s (n) j : = s (n) j − s (n) j−1 is constant, that is s (n) j = T/n = (n) = : , and s (n) j = j for each j = 0, . . . , n.

The Lévy Simulations
By definition, for a fixed time interval , the increments of the Lèvy process L(s (n) j ) : = L(s (n) j ) − L(s (n) j−1 ), j = 1, . . . , n, are independent and identically distributed random variables following a certain distribution, which depends only on the time increment . For example, for L = B, the increments B(s (n) j ) have all Gaussian distribution with zero mean and variance . A possible approach to simulate L ε in an exact way is to find the distribution of L ε (t), starting from the distribution of L(t). For example, for L = B, by means of Lemma 2.3, we can simulate B ε sampling directly from a Gaussian distribution with zero mean and variance in Lemma 2.3. However, with this approach, L and L ε are from two different ω's in , and comparing the two paths we would not get the error measure required. We then need an alternative approach that simulates L ε and L from the same ω ∈ .
From Equation (2.3), L ε is the stochastic integral of φ ε with respect to the Lévy process L. From the definition of the stochastic integral as a limit of approximating Riemann-Stieltjes sums [20], for t ∈ (n) such that t = s (n) m = m , for m = m(t) and 1 ≤ m ≤ n, we approximate L ε (t) by approximating the integral with a sum, namely The approximation scheme (4.1) adds an error to the simulation. It, however, constructs L ε starting from the simulations of L, allowing to compare the two paths for different 's. In the next proposition, we study the approximation error. Proposition 4.1. For every > 0, ε > 0 and T < +∞, there exists a constant K 3 > 0 such that the simulation error is bounded by Proof: See Appendix A.10 in Supplementary Material.
With Proposition 4.1, we estimated an upper bound for the error added by the Euler scheme for simulating the process L ε by means of the process L ε . In particular, it shows that, for a fixed ε > 0, the simulation error decreases at the decreasing of the time step , as we would expect. In Kloeden et al. [21, section 9.3], the authors study the rate of convergence for the Euler scheme when applied to an Itô process, with respect to the absolute error criterion, finding a strong rate of convergence equal to 1/2. Here, in line with the rest of the paper, we consider the squared error in Proposition 4.1 instead. This gives a rate of convergence equal to 2, corresponding to a strong order 1. From Proposition 4.1, we also notice that the bound depends on the inverse of ε, meaning that, for ε going to zero, the simulation error might explode. Figure 5 shows the convergence of the Euler scheme (4.1) for a Brownian motion and for the normal inverse Gaussian process introduced in section 2.1. Since an exact simulation is not possible, to generate B ε , we consider the scheme (4.1) with a small time step true = 2 −10 . Moreover, since we want to compare B ε and B ε path-wise for different values of , B ε and B ε must be related to the same ω ∈ , that is, to the same realization of the Brownian motion B. To get this, we exploit the following well-known property: given X i ∼ N (0, σ 2 i ), for i = 1, ..., k, independent Gaussian random variables, their sum follows the distribution k i=1 X i ∼ N (0, k i=1 σ 2 i ). Thus, we start by simulating the true process B ε for a certain time step true << 1. For this, we sample the increments B(s (n) j ) from a Gaussian distribution with zero mean and variance true and obtain B ε as in Equation (4.1). Then, we choose a time step multiple of true , namely, = k true for a certain k ∈ N. The Brownian motion path to be used in Equation (4.1) for simulating B ε is then obtained by summing up the increments of B (above simulated) at groups of k elements. These increments have indeed Gaussian distribution with zero mean and variance k i=1 true = k true = , as required. In Figure 6, we considered true = 2 −10 and ∈ {2 −7 , 2 −5 , 2 −3 , 2 −1 }. We see that, for a fixed ε, as decreases, the error between B ε and B ε decreases. Even with a large ( = 2 −1 ), the process B ε captures the shape of B ε . In Figure 3E we show the rate of convergence of the Euler scheme proven in Proposition 4.1 as function of time step and for a fixed ε (ε = 0.01).
Similarly, for the NIG Lévy process, to simulate and compare L ε with L ε path-wise, we use the same idea as in the Brownian motion case. We can indeed exploit the convolution property of the NIG distribution (see [22]): given X i ∼ NIG(α, β, δ i , µ i ), for i = 1, ..., k, independent normal inverse Gaussian random variables with common parameters α and β, but individual scalelocation parameters δ i and µ i , their sum follows In section 2.1, we introduced NIG distribution with parameters (1, 0, true , 0), for true being the time step used in the Euler discretization. As above, we thus simulate the true process L ε from Equation (4.1) by sampling the increments of L from a NIG distribution with parameters (1, 0, true , 0) for a certain true << 1. Then, we obtain L ε by choosing so that = k true for a certain k ∈ N, and summing up the increments of L (simulated above) at groups of k elements. These increments have indeed NIG distribution with parameters (1, 0, k i=1 true , k i=1 0) = (1, 0, , 0) as required. We consider the same 's as for the Brownian motion. We notice that, for a fixed ε, L ε captures the shape of L ε even with a large . In Figure 3F we show the rate of convergence for the Euler scheme proved in Proposition 4.1 as function of the time step and for a fixed ε (ε = 0.01).
From Proposition 2.7 and 4.1, we calculate the squared error for L ε with respect to L. Proposition 4.2 shows that the squared error between L ε and L decreases when decreases. However, because of the term proportional to the inverse of ε, it might also explode for small values of ε. We then need to tune the parameters and ε to find a combination that minimizes the error preventing it from taking large values.

The CARMA Simulations
It is possible to construct X and X ε starting from Equations (3.5) and (3.8) and considering the same kind of numerical approximation as done for L ε , that is, by replacing the integral with a suitable sum. However, to not add a second simulation error to the one already due to L ε , which is approximated by L ε , we start from the state-space representations instead, Equations (3.3) and (3.7), respectively. By using the Euler scheme, according to [21,Chapter 9], we get (here we write ε in parenthesis since the equation holds both for X and X ε ): For the reasons already discussed, we apply the scheme (4.2) only for q = 0, that is, when b = e 1 . For q > 0, by means of Proposition 3.8, we instead apply the Euler scheme to the stochastic differential equation (3.12) for each of the p components X (r) (ε) of X (ε) , and we get X (ε) as the sum of them, namely, When simulating X ε , we use the approximation L ε in Equation (4.1), and this adds an error. We study the rate of convergence for X ε to the true process X ε in the following proposition. |λ j | 2 , there exist two constants K 4 , K 5 > 0 such that the simulation error for the Euler scheme (4.2) is bounded by Proof: See Appendix A.11 in Supplementary Material.
Let us notice that the bound on in Proposition 4.3 has been introduced to obtain the estimate for the squared error. The scheme (4.2) works however for every > 0.
In Figure 7 we show the convergence for the simulations of a Gaussian CAR(1) and a Gaussian CARMA (3,2). In both cases we notice that for decreasing, the Euler approximation X ε converges to X ε , and, even with a large value of ( = 2 −1 ), X ε captures the shape of X ε . The underlined digits are the common digits inâ andâε . The symbol * marks when condition ε < is satisfied.
In Figures 3G,H, we show the rate of convergence for the two processes as function of 2 and for a fixed value of ε (ε = 0.01).
From Proposition 3.6 and 4.3, we estimate the squared error for X ε with respect to X. Proposition 4.4. For every > 0 and ε > 0 small enough, the following bound applies: Proof: This is a direct consequence of Proposition 3.6 and 4.3 by the inequality (a + b) 2 ≤ 2(a 2 + b 2 ) for a = X ε (T) − X ε (T) and b = X ε (T) − X(T).

CARMA Estimation
We introduce an approach to estimate the parameters (a 1 , . . . , a p , b 0 , . . . , b q ) starting from the time series for {X(t), 0 ≤ t ≤ T} and {X ε (t), 0 ≤ t ≤ T} by applying the Euler scheme. Following Benth et al. [7, section 4.3], we focus on CAR models only, that is, on the process X(t) = Y 1 (t), for Y 1 (t) being the first coordinate of the vector Y(t) in Equation an Ornstein-Uhlenbeck process and Equation (4.3) becomes which coincides with the Euler scheme (4.2) when p = 1. For p = 3, from Equation (4.3) we get From the definition for z i j , rearranging the terms, the last equation can be rewritten as with θ 1 : = 3 − a 1 , θ 2 : = 2a 1 − a 2 2 − 3, Then, one first estimates the AR coefficients θ 1 , θ 2 , θ 3 and afterwards recover the CAR coefficients a 1 , a 2 , a 3 by This scheme holds for every Lévy process. To recover the case q > 0 different approaches have to be considered, see, e.g., [23] and [24].
In Tables 1, 2 we report the parameters estimated for a Gaussian CAR(1) and a Gaussian CAR(3) for different ε's and different 's. In particular,â is the vector of parameters estimated starting from the time series for X, andâ ε is the vector of parameters estimated starting from the time series for X ε . The underlined digits are the common digits betweenâ andâ ε , with respect to the same ε and the same . The estimates get more accurate for a small ε, but they are not accurate if ε is bigger than . We marked the cases corresponding to ε < with the symbol * , and we notice that as soon as ε gets smaller than , the accuracy improves and gets better for smaller values of ε. We summarize this empirical result with a proposition. A more rigorous study should be performed regarding the condition in Proposition 4.5, analysing in details the estimation technique for ARMA processes, which is the starting point to estimate the CARMA parameters. This is, however, beyond the purposes of the current article.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.