We consider the superpositions of Ornstein–Uhlenbeck (supOU) processes which are an important class of ambit stochastic processes [11] and play a fundamental role in finance and insurance, degradation modelling, diffusion, turbulence, astrophysics and transport, see [20, 52, 30, 51] and references therein. SupOU processes are proving to be invaluable tools for modeling complex phenomena with long-range dependence, intermittency, and non-Gaussian behaviors.

In finance, the supOU processes can be incorporated in different ways. The popular continuous time stochastic volatility model for financial assets d S ( t ) = μ S ( t ) d t + σ ( t ) S ( t ) d B ( t ) , t ≥ 0 , generalizes the classical Black–Scholes model by replacing the volatility parameter σ by the volatility process σ ( t ) in the form of a superposition of positive OU processes, where B ( t ) is the standard Brownian motion that is independent of σ ( t ), see the fundamental work [8]. If σ ( t ) is a supOU process with inverse gamma marginals, then log-returns have approximately the Student distribution, see [8, p. 170] for a general correspondence.

The fractal activity time geometric Brownian motion model introduced in [35] has the form S ( t ) = S 0 e μ t + θ T ( t ) + σ B ( T ( t ) ) , t ≥ 0 , where S ( t ) is the asset price and the fractal activity time T ( t ) is a positive, nondecreasing stochastic process in the form of the linear spline of a discrete integrated supOU process with positive marginals that is independent of B ( t ). If T ( t ) is an integrated supOU process with inverse gamma marginals, then log-returns have exactly the Student distribution, see [27, 28] for other marginals.

For modeling high-frequency financial data, supOU processes were explored in [37] by demonstrating their efficacy in capturing the stochastic volatility of stock prices with improved accuracy over traditional models like GARCH. The authors leveraged the superposition property to model multiscale temporal dependencies, which proved particularly effective for intraday trading strategies. Furthermore, supOU processes were applied to stock indexes [36], zero-coupon bonds [61], and option pricing [48, 59].

SupOU processes are well suited for modeling migrating fish counts [62] and intricate hydrological phenomena, such as river discharge time series [65, 63] and water characteristics [64]. In the context of wind energy systems, the utility of graph supOU models in capturing the spatial and temporal variability of wind capacity factors across a European electricity network was demonstrated in [50]. In astrophysics, the supOU processes were applied to modelling the X-ray light curves generated by black holes [42].

For the main properties of supOU processes, we refer to several works by Barndorff-Nielsen and his collaborators [5, 9, 14, 6]. In particular, for every self-decomposable (SD) distribution, there is a stationary supOU process with this SD distribution as marginal distribution and quite flexible class of covariance structures including long-range dependence (LRD), that follows from the following theorem.

We note that Theorem 1 is a reformulation of Theorem 3.9.3 in [41]. Let us now give the definition of supOU processes.

The quadruple ( b , σ 2 , μ Z , π ) is referred to as the characteristic quadruple and it completely determines the supOU process (1) with any given marginal SD distribution with Lévy triplet ( b Y , σ Y 2 , μ Y ), that is, (3) κ Y ( z ) = i z b Y − 1 2 z 2 σ Y 2 + ∫ R ∖ { 0 } ( e i z x − 1 − i z x 1 [ − 1 , 1 ] ( x ) ) μ Y ( d x ) . Recall that the ID random variable V = d Y ( 0 ) with characteristic function ϕ ( z ) = exp ( κ Y ( z ) ) is self-decomposable if for every constant c ∈ ( 0 , 1 ) there is a characteristic function ϕ c ( z ) such that ϕ ( z ) = ϕ ( c z ) ϕ c ( z ), z ∈ R. The ID random variable V is self-decomposable if and only if its Lévy measure has the form μ Y ( d x ) = g ( x ) | x | d x, where g ( x ) is increasing in ( − ∞ , 0 ) and decreasing on ( 0 , ∞ ), see [58, 54]. We also note that κ Z ( z ) = z d d z κ Y ( z ).

If E Y 2 ( t ) < ∞, the correlation function of the supOU process (1) is given by (4) r ( τ ) = Corr ( Y ( t ) , Y ( t + τ ) ) = η ∫ 0 ∞ ξ − 1 e − τ ξ π ( d ξ ) , τ ≥ 0 , which can exhibit either short-range dependence (SRD) or LRD depending on the measure π.

For explanation of the concept of the superposition, let us introduce the Lévy-driven Ornstein–Uhlenbeck process X ( t ) = ∫ − ∞ t e − ξ ( t − s ) d Z ( ξ s ) , t ∈ R , where Z ( t ) is a two-sided Lévy process with Lévy–Khintchine triplet ( b , σ 2 , μ Z ).

The correlation function of the OU process X ( t ) is of the form r ( τ ) = Corr ( X ( t ) , X ( t + τ ) ) = e − ξ τ , τ > 0 , provided E X 2 ( t ) < ∞.

Thus, if π is a probability measure, then the measure π provides a randomization of the rate parameter ξ of the Lévy-driven OU process X ( t ) and hence the correlation structure (4) of the supOU process is more flexible and can exhibit LRD. The popular example of the measure π is the gamma distribution, some other examples can be found in [7]. Although the OU process X ( t ) is a linear process, the supOU process is nonlinear in general, that is, Y ( t ) cannot be represented in the form ∫ − ∞ t a ( t − s ) d Z ( s ) for some function a ( · ), see [23, 5, 47].

In [38–40] self-decomposability was shown for the gamma distribution, the inverse gamma distribution, the inverse Gaussian distribution, the Student distribution, the hyperbolic cosine distribution, the Gumbel distribution and the Bessel distribution. Also, self-decomposability was shown for the generalized inverse Gaussian distribution in [34], for the generalized Gaussian distribution in [24], for the positive α-stable distribution in [9], and for the generalized Linnik distribution in [3].

Self-decomposability of the log-normal distribution follows from [18], where it is proven that the log-normal distribution belongs to the class of generalized gamma convolutions, which is a subclass of SD distributions. Finally, self-decomposability of the Rosenblatt distribution has been proven in [49].

Let us present the fundamental result which explains the structure of supOU processes.

The representation of ν in the proof of Theorem 6 is most convenient for us among several series representations of random measures which are reviewed in [56].

The most famous example of the application of Theorem 2 is the gamma process, however, the relation (6) does not work for many other marginal distributions, see Section 2. We note that the equation (7) is called the Bondesson–Rosinski representation due to [16, 56], which is not unique, if we change the intensity θ, then the shot-height distribution F ( x ) should be scaled. If R k = λ for all k, then the formula (7) becomes the representation of the usual OU process. Since S k are jump times of the Poisson process with intensity θ, we have S k − S k − 1 ∼ Exp ( θ ) , that is, S k − S k − 1 has the exponential distribution with rate θ (implying mean 1 / θ ). Thus, the Bondesson–Rosinski representation (7) has in average θ summands per any unit interval. We note that S k / η are jumps of the Poisson process with intensity θ / η. In Figure 1 we show a typical sample path of the supOU process Y ( t ), which consists of jumps with different heights and mixed-exponential decays with different rates.

Overall, we can see that the supOU processes form a wide class of stochastic processes and the representation (7) enables to independently specify the marginal distribution by the choice of the shot-height distribution and the correlation structure by the choice of the measure π.

We note that the Lévy measure μ Y for the specified marginal distribution of the process Y ( t ) can be derived from the one-to-one correspondence between the Lévy density and the characteristic function as stated in the following result.

Let us now reformulate Theorem 3 for the supOU processes Y ( t ) with positive marginals using the Laplace transform ψ ( s ) = E e − s Y ( t ) = ∫ 0 ∞ e − s x p ( x ) d x, s > 0, where p ( x ) is the density of the marginal distribution. Following [17, p. 15], the density p ( x ) is infinitely divisible if and only if ζ ( s ) = log ψ ( s ) = − a s + ∫ 0 ∞ ( e − s x − 1 ) m ( x ) d x , where a ≥ 0 is the left-extremity of the distribution. Taking the first and second derivatives of the above identity with a = 0, we obtain (10) ζ ′ ( s ) = − ∫ 0 ∞ e − s x x m ( x ) d x , (11) ζ ″ ( s ) = ∫ 0 ∞ e − s x x 2 m ( x ) d x , which can be interpreted as the Laplace transforms of x m ( x ) and x 2 m ( x ), respectively. Therefore, the Lévy density m ( x ) can be found by the inversion of either (10) or (11) using one of several methods: the Bromwich integral, the Post–Widder formula, or Fredholm equations, see [21] for details. In particular, the Lévy density for the log-normal distribution was numerically computed by the method of Fredholm equations in [18]. In our numerical examples, where the Lévy density is not known analytically, we use the method based on the Bromwich integral.

The following theorem establishes another useful relation for the OU-type processes, one that does do not involve the Lévy measure and is motivated by the theory of shot-noise processes [58].

From Theorem 4 for the exponential impulse response h ( t ) = e − R t , we have log ϕ ( z ) = θ R ∫ − ∞ ∞ K ( z x ) d F ( x ) , K ( u ) = ∫ 0 u e i s − 1 s d s .

We will concentrate on simulation of supOU processes with fat-tail marginals such as the inverse gamma distribution and the inverse Gaussian distribution. We note that supOU processes with fat-tail marginals usually have the Lévy density m ( x ) such that x m ( x ) is singular at zero that disables the direct use of (6).

We propose to adapt the Bondesson–Rosinski representation (7) for simulation of the supOU processes which may have infinitely many jumps. In numerical studies, we usually simulate a stochastic process at equidistant points t j = j Δ with time step Δ and, therefore, the contribution of frequent small jumps between points t 1 , t 2 , … can be disregarded that can be viewed as the truncation of the Lévy measure.

Indeed, consider the Lévy measure μ Y on the interval [ 0 , ∞ ). This measure on the interval [ 0 , ϵ ] corresponds to small jumps of the supOU process Y ( t ). We approximate Y ( t ) by a stochastic process Y ϵ ( t ) with the Lévy measure μ Y truncated to the interval [ ϵ , ∞ ). Specifically, we construct Y ϵ ( t ) as follows.

For the given Lévy density m ( x ) of the process Y ( t ), from the relation (6), we obtain F ( x ) = 1 − x θ m ( x ) , which gives negative values near zero if x m ( x ) is singular at zero, and F ( x ) is an increasing function with F ( x ) → 1 as x → ∞. Solving F ( ϵ ) = 0 with respect to θ for given ϵ ∈ ( 0 , ∞ ), we get the solution θ ϵ = ϵ m ( ϵ ). By truncating F ( x ) to the interval [ ϵ , ∞ ), we define the function (14) F ϵ ( x ) = 1 − x θ ϵ m ( x ) , x ≥ ϵ , 0 , x < ϵ , which is a cumulative distribution function. Now we will state a main theorem which will be used for simulation.

Since the probability density function (pdf) is the inverse Fourier transform of the exponent of the cumulant function, the pdf of the marginal distribution of Y ϵ ( t ) converges to that of Y ( t ). Moreover, the finite-dimensional distributions of the process Y ϵ ( t ) converge to those of Y ( t ) and, consequently, the statistical properties of sample paths of Y ϵ ( t ) approximate those of Y ( t ) for small ϵ. Since the supOU process Y ϵ ( t ) can be viewed as a surrogate model of Y ( t ), the sampling error for modelling the stochastic process Y ( t ) on a finite interval is larger than the error due to the replacement of the model Y ( t ) by Y ϵ ( t ) with small ϵ, see Supplementary Materials in [45] for details. From a practical point of view, we can take ϵ = 0.001. For other Lévy processes, the convergence under truncation of the Lévy measure is investigated in [29].

In summary, the simulation algorithm for a supOU process with the given measure π and the Lévy density m ( x ) of the given marginal distribution is as follows.