A causal continuous-time stochastic model for the turbulent energy cascade in a helium jet flow
We discuss continuous cascade models and their potential for modelling the energy dissipation in a turbulent flow. Continuous cascade processes, expressed in terms of stochastic integrals with respect to Lévy bases, are examples of ambit processes. These models are known to reproduce experimentally observed properties of turbulence: The scaling and self-scaling of the correlators of the energy dissipation and of the moments of the coarse-grained energy dissipation. We compare three models: a normal model, a normal inverse Gaussian model and a stable model. We show that the normal inverse Gaussian model is superior to both, the normal and the stable model, in terms of reproducing the distribution of the energy dissipation; and that the normal inverse Gaussian model is superior to the normal model and competitive with the stable model in terms of reproducing the self-scaling exponents. Furthermore, we show that the presented analysis is parsimonious in the sense that the self-scaling exponents are predicted from the one-point distribution of the energy dissipation, and that the shape of these distributions is independent of the Reynolds number.
Since the pioneering work of Kolmogorov  and Oboukhov , where the turbulent energy dissipation is assumed to be log-normally distributed, the small-scale intermittency of the energy dissipation in turbulence has received much attention [24, 44]. The small scale intermittency is primarily expressed in terms of multifractal and universal scaling of inertial range statistics, including extended self-similarity , scaling and self-scaling of correlators , and the statistics of breakdown coefficients .
Early attemps to model the rapid variation of the turbulent velocity field include [45, 9, 23, 28, 12, 1] (among many others). Such phenomenological approaches are sometimes called “synthetic turbulence” and can be divided into two classes. The first direction starts from modelling the velocity field and derives the model for the energy dissipation by taking squared small scale increments. The second line of investigation focuses on modelling the energy dissipation field and derives the velocity field by various, partly ad hoc, manipulations. The approach presented here focuses on modelling the energy dissipation as the fundamental field, not derived from an a priori velocity field.
In , an iterative, geometric multi-affine model for the one-dimensional velocity process is contructed and some of the basic, global statistical quantities of the energy dissipation field are derived. However, this discrete, dyadic approach does not allow to give explicit expressions for more specific statistical quantities. In the approach discussed here, all main statistical quantities, like -point correlations and probability densities can be derived analytically since our approach defines the energy dissipation field as an explicit closed expression that is mathematically tractable and does not involve an iterative procedure.
Another dyadic, iterative approach for the construction of the velocity field is discussed in . Their model is based on a wavelet decomposition of the velocity field combined with a multiplicative cascading structure for the wavelet coefficients. The energy dissipation field is then derived from small scale increments, again leaving only limited room for an exact analytical treatment of higher order statistics. However, as discussed in , such wavelet approaches are superior over discrete geometric approaches as they allow to model stationarity in a mathematical more rigorous way. The approach discussed here does not suffer from problems related to mathematical rigour and no iterative limit arguments are needed for the construction. A related and interesting wavelet-based approach is discussed in , which allows for a sequential construction of the field. A further wavelet-based approach  builds on random functions and their orthogonal wavelet fransform. The authors show that to each such random function there is an associated cascade on a dyadic tree of wavelet coefficients. The performance of the model is illustrated by numerical examples, with little analytical insight.
The models [45, 9] fail to incorporate skewness for the velocity increments , a basic property of turbulent fields. As an alternative approach,  proposes a combination of a multiplicative cascade for the energy dissipation, the use of Kolmogorov’s refined similarity hypothesis  and an appropriate summation rule for the increments to construct the velocity field. Here, again, only discrete iterative procedures are employed which make analytical statistical statements very difficult.
It is important to mention that the continuous cascade fields discussed in this paper are potentially useful in the above cited works [9, 23, 28, 12, 1] as a candidate for a closed continuous and mathematically tractable and flexible version of the discrete multiplicative procedures employed there.
Discrete and continuous random cascade processes have proved useful in describing phenomenologically the small-scale behaviour of the turbulent energy dissipation [31, 25, 11, 38, 32, 27, 18, 26]. In  the surrogate energy dissipation is modelled as a discrete random multiplicative cascade process. Choosing the law of the cascade generators to be log-normal yields the Kolmogorov-Oboukhov model. A continuous analogue to the discrete multiplicative cascade processes is formulated in terms of integrals with respect to Lévy bases and has been shown [4, 41, 39] to be computationally tractable and to accurately describe the two- and three-point statistics of the energy dissipation.
In the cited works, focus is on the modelling of -point statistics of the energy dissipation, not the distribution of the energy dissipation itself. Indeed,  concludes with a remark that the law of the Lévy basis driving the cascade model should be inferred and its dependency on the Reynolds number should be investigated.
Both, discrete and continuous multiplicative cascade processes, suggest that the law of the logarithm of the energy dissipation should be infinitely divisible. Infinite divisibility is necessary for the cascade models discussed here to be defined in a mathematical rigorous way. Furthermore, it greatly simplifies analytic calculations of some of the statistical properties of the models. Among the infinitely divisible distributions are the normal, stable, and normal inverse Gaussian distributions. These three classes of distributions each have their own tail behaviour.
The use of stable Lévy bases for modelling of the energy dissipation has been investigated in  and it is concluded by analysing the breakdown coefficients that “except for the log-normal limit, this leaves no room for the log-stable modelling of the turbulent energy cascade.” The present paper investigates the alternative of using a normal inverse Gaussian Lévy basis to model the energy dissipation and addresses the one-point distributions, multifractality of the coarse-grained energy dissipation, and two-point statistics.
The use of normal inverse Gaussian distributions in turbulence modelling is not new. In , the parameters of the normal inverse Gaussian distribution cascade generator are estimated from scaling exponents and cumulants, which are moment estimates, notorious for their sensitivity to outliers. Indeed, estimation of the normal inverse Gaussian parameters (and those of other distributions) may not be feasible from sample moments. In this paper we will apply likelihood methods which do not suffer from the same problems as the moment based methods.
Our motivation for the use of the normal inverse Gaussian distribution is not based on physical arguments. We mainly exploit that these distributions are flexible, yet simple, and capable of attaining a wide range of shapes. Some shortcomings of the normal inverse Gaussian distribution will be discussed. Being a normal mean-variance mixture with an inverse Gaussian as mixing distribution, the tail behavior of the mixture is related to that of the mixing distribution . Hence one obtains a recipe for the construction of distributions with a prescribed tail behaviour. The infinite divisibility is necessary for the calculus of the stochastic integrals used in this paper.
The She-Leveque-Dubrulle model [42, 21, 43], which has been shown to accurately predict the scaling exponents of the coarse-grained energy dissipation and the structure functions of the velocity increments, prescibes a Poisson distribution for the log-energy dissipation. We find that the Poisson distribution provides a poor fit and that the normal inverse Gaussian distribution is a clear improvement. It should be noted that, being infinitely divisible, the Lévy-Itô decomposition expresses the normal inverse Gaussian distribution as a linear combination of Poisson distributions.
The paper is organised as follows. Section 2 provides some background on the data analysed in this paper. Section 3 recalls the construction of continuous cascade processes in terms of integrals with respect to Lévy bases. Section 4 applies the theory to the data and shows how the distribution of the surrogate energy dissipation determines the scaling and self-scaling exponents of the two-point correlators of the energy dissipation and the coarse-grained energy dissipation. Section 5 concludes. The two appendices provide necessary background on the normal inverse Gaussian distribution and integration with respect to Lévy bases.
2 Background on the data
We analyse thirteen data sets, each consisting of approximately sixteen million one-point time records of the velocity component in the mean stream direction in helium gas jet flow . The time series can be assumed to be stationary. In , eleven of the thirteen data sets are used in an analysis of the intermittency exponent of the turbulent energy cascade. In particular, [19, table I] summarises useful information about the data sets. Data sets no. 6 and 9 are not considered in .
Let denote the velocity component in the mean stream direction, and let denote the mean stream velocity. Since the flow can be assumed to be homogeneous and isotropic, we use the surrogate energy dissipation as a proxy for the energy dissipation and henceforth omit the “surrogate” predicate. Here denotes the position along the mean stream direction, and denotes the viscosity. We apply the Taylor frozen flow hypothesis to express the energy dissipation in terms of the measured time series. We do not invoke an “instantaneous Taylor correction”  since it introduces spurious effects in correlators of the energy dissipation. Since any change of the energy dissipation by a multiplicative constant is inconsequential for the conclusions, we scale the energy dissipation to have unit mean. Finally, the derivative is calculated from the discrete samples using interpolation with third-order splines.
The resolution scale (mean velocity times sampling time ) is – times the Kolmogorov length , and the Taylor-microscale based Reynolds number varies from to . Table 1 lists in addition to parameters that will be explained later.
The spectral density of the velocity component in the mean stream direction is estimated using Welsh’s overlapping segment averages with a Hanning taper, a block length of , and an overlap of , see e.g.  for details about the method. Figure 1 shows the spectral density for data sets no. 1 and 13. The inertial ranges are found at frequencies from approximately to and to , respectively, after which the transition to the dissipation range occurs. The slope in the inertial range is approximately , obeying Kolmogorov’s -law. For most of the data sets, the spectral density attains an almost constant value at the highest frequencies. We interpret this as instrument noise. As shown by fig. 2 (top) this does not appear to distort the data significantly. For some of the data sets, the resolution scale is almost five times the Kolmogorov length. The effect of this becomes apparent in fig. 2 (bottom) where the time series appears much less smooth than in fig. 2 (top) where the resolution scale is close to the Kolmogorov length. If data set no. 1 is subsampled by retaining every fourth record, a picture similar to fig. 2 (bottom) is obtained. Having no means to improve the resolution scale of the measured data sets, we consider downsampling the data to improve the smoothness so that calculation of the derivative becomes reliable. As will be remarked later, the influence of the downsampling on the shape of the one-point distribution of the energy dissipation is negligible and it significanly decreases the scatter of all quantities derived from moments estimates. Henceforth all graphs and tables refer to data that has been downsampled by a factor of two. We found little or no change when downsampling by a factor of three or four compared to a factor of two.
To assess the existence of moments, we consider, for a given time series , the relative absolute sample moment of order and sample size , where . Figure 3 shows the relative absolute sample moments of of orders 6 and 8 for data set no. 7. We see that the sample moments are corrupted by outliers and that removal of outliers increases the reliability of the moment estimates. We will assume that has moments of order at least 8. The removal of outliers does not distort the parameters of the one-point distribution of the energy dissipation.
3 A Lévy based continuous multiplicative cascade model
In this section we present the model for the energy dissipation. The model provides a link between the distribution and the two-point correlators of the energy dissipation.
where is called the ambit set. We will assume that the Lévy basis is homogeneous (see appendix B), and that the ambit sets are defined by for some bounded set . Thus we ensure that is stationary in space and time. The model (2) is an example of a random multiplicative cascade process in continuous space and time. Since only one-point time series are available for our analysis, we will ignore the spatial dependency and consider the energy dissipation as a function of time alone and at a fixed position in space,
where . However, it must be emphasised that the model is equally capable of modelling specific spatial statistics in one or more spatial dimensions.
3.1 One-point distributions
where denotes the logarithm of the Laplace transform of the random variable , denotes the Lévy seed corresponding to the Lévy basis , and denotes the volume (or area) of the set . By (3), the distribution of the energy dissipation is closely related to the distribution of the Lévy seed . For brevity we define and .
3.2 Coarse-grained energy dissipation
Let denote the coarse-grained energy dissipation,
Kolmogorov’s refined similarity hypothesis states that
where is a random variable independent of the energy dissipation. This hypothesis provides a relation between the scaling exponents of the moments of longitudinal velocity differences and the scaling exponents of the moments of the coarse-grained energy dissipation, namely
The exponent is the intermittency exponent which quantifies the deviation from Kolmogorov’s 1941 theory where .
for integral values of where is a constant depending only on the size and shape of the ambit set . It follows that the self-scaling exponent is independent of the ambit set and hence only depends on the distribution of the Lévy seed .
3.3 Two-point correlators
The two-point correlator of order of the energy dissipation is defined as
We observe that the exponent in (9) is expressed as a product where the first factor depends only on the Lévy basis and the order of the correlator, and the second factor depends only on the overlap of the ambit sets. This provides a way of modelling a wide range of correlators, since the shape of the ambit set, under suitable assumptions, can be determined from the correlator, see  for details or subsec. 3.4 for a particular example.
in a range of comparable to the inertial range of the velocity structure functions. It is straightforward to show that the scaling exponent is equal to the intermittency exponent . This follows from the fact that the second derivative of with respect to behaves as .
In  it is shown that the two-point correlators also enjoy the property of self-scaling,
for an even wider range of , in analogy to extended self-similarity  of structure functions. Here is the self-scaling exponent.
is the self-scaling exponent. As noted in , the self-scaling property is independent of the shape of the ambit set and thus scaling of the correlators is not necessary for self-scaling of the correlators. The Lévy seed determines through (3) the distribution of the energy dissipation and therefore the self-scaling exponents. Both, the distribution and the self-scaling exponents, can be estimated from data and hence compared with the model.
The correlators are moment estimates. Therefore they may not exist beyond a certain order. When they exist, they are sensitive to noise and outliers, particularly at high orders. By Hölder’s inequality we have
It follows that the correlator exists provided the velocity derivative has finite moment of order . Under the model (2) it follows by (3) that exists if and only if the Lévy seed has exponential moments of order .
3.4 Ambit sets and scaling of correlators
While the self-scaling of the correlators and the self-scaling of the coarse-grained energy dissipation is independent of the shape of the ambit set , the scaling property (10) requires to be specified appropriately according to (9). In this subsection we consider one such specification. We define for , , and the function by
In the limit we have that
The ambit set is given as
It follows that
Hence is the scaling exponent of . In the limit we have perfect scaling,
The parameter is merely a tuning parameters to account for imperfect scaling of the correlators. In view of (9), where the overlap of the ambit sets determines the correlator, may be interpreted as the decorrelation time, and is a time scale near and below which the scaling behaviour has terminated.
3.5 Three Lévy bases
We consider three distributions of the Lévy seed :
The normal Lévy seed yields a log-normal model for the energy dissipation, consistent with the Kolmogorov-Oboukhov theory . The stable Lévy seed is included to enable comparison with . The use of a normal inverse Gaussian Lévy seed can be motivated in several ways, though none of them is based on physical arguments. The normal inverse Gaussian distribution is computationally tractable and allows a wide range of distributions (see fig. 11). Furthermore, the representation (42) shows that the normal inverse Gaussian distributions form a natural generalisation of the normal distributions by allowing the variance of the normal distribution to be random. Since the variance is positive, and if there is a “typical variance”, it is natural to model the variance using an unimodal distribution on the positive real line. A flexible distribution of this type is the inverse Gaussian law as shown in fig. 12.
While other infinitely distributions could be considered as well, the three classes of distributions above exhaust a wide range of unimodal distributions with a specific tail behaviour: The normal distributions have “light tails”, the stable distributions have “heavy tails”, decaying algebraically, and the normal inverse Gaussian distributions have “semi-heavy tails”, decaying exponentially, see (36).
The common symbols for the parameters have been chosen due to the similarity of their interpretations. In all three cases, is a location parameter, is a scale parameter, and are shape parameters, and determines the asymmetry of the distribution. The domain of the parameters is and in all three cases; and in the stable case; and in the normal inverse Gaussian case. The parametrisation of the stable distribution is chosen to follow [37, eq. (1.1.6)], so that the log-characteristic function is given by
We see that the shape parameters of the distribution of are, in all three cases, identical to the shape parameters of the distribution of the Lévy seed , and that the location and scale parameters of are multiplied with a factor determined by the size of the ambit set.
The normal distribution has exponential moments of all orders. Provided , the stable distribution has exponential moments of all non-negative orders. We will therefore assume that in the stable case. The normal inverse Gaussian distribution has exponential moments of order if
In these cases the log-Laplace transform of the Lévy seed is given by
It follows that
Note that under the normal inverse Gaussian model, we must require that to ensure the existence of the correlator of total order . Note also that depends only on the order and the shape parameters. The self-scaling exponents of the correlators are now given by
In a similar manner it follows that the self-scaling exponents of the coarse-grained energy dissipation depend only on the order and the shape parameters of the Lévy seed. We therefore conclude that both kinds of self-scaling exponents are uniquely determined by the orders under consideration and the shape of the one-point distribution of the energy dissipation.
4 Data analysis
In this section we apply the model from section 3 to the thirteen data sets. The one-point distribution of the energy dissipation is used to predict the scaling and self-scaling exponents of the two-point correlators and the coarse-grained energy dissipation.
4.1 One-point distributions
Figure 5 shows that the distribution of has a distinct, non-Gaussian, asymmetric shape. The shape is independent of the downsamling applied to smooth the data. Four parametric distributions are fitted to the data using likelihood methods: normal, stable, normal inverse Gaussian, and normal inverse Gaussian constrained to have finite exponential moment of order according to (25). Clearly, the normal inverse Gaussian distribution provides an excellent fit for all amplitudes. The left tail is overestimated by the stable distribution. The unconstrained normal inverse Gaussian distribution is also able to accurately fit the distribution, though it slightly underestimates the left tail. In what follows, any reference to the normal inverse Gaussian distribution will imply the constrained version.
It is necessary to constrain the shape parameters of the normal inverse Gaussian distribution to ensure the existence of the correlators up to a certain total order. This constrainment is a delicate issue. Firstly, it introduces the problem of determining an appropriate upper order. And secondly, it implies a slightly poorer fit. The choice of as the upper order is somewhat arbitrary but provides a compromise between fitting the distribution of the energy dissipation accurately and predicting the self-scaling exponents of the correlators and the coarse-grained energy dissipation. We conjecture the existence of a more appropropriate infinitely divisible distribution, perhaps within the class of normal mean-variance mixtures. For the present paper, the normal inverse Gaussian will suffice.
Table 1 summarises the estimated parameters for each of the thirteen data sets. Note that for both, the normal inverse Gaussian distribution and the stable distribution, the estimated shape parameters do not depend on the Taylor micro-scale Reynolds number, while the location and scale parameters vary only slightly. This is a clear indication of universality of the distribution of the energy dissipation, at least within the thirteen data sets considered here. Downsampling the data by a factor of two changes the normal inverse Gaussian shape parameters and by only around and the stable shape parameter by around .
4.2 Coarse-grained energy dissipation
The behaviour of the moments as a function of the lag does not reveal any pronounced scaling range, see fig. 6 (top) as an example. Compare with fig. 8 which demonstrates a clear scaling range for the two-point correlators. However, as shown by fig. 6 (bottom) the moments display a clear self-scaling behaviour, and the self-scaling exponents can be extracted.
The estimated self-scaling exponents can now be compared with the predictions of the model (2) and with the She-Leveque-Dubrulle model. The former predicts
while the latter predicts
Figure 7 (top) shows that the normal inverse Gaussian and the stable model agree well with both data and the She-Leveque-Dubrulle model for the self-scaling exponent . For fig. 7 (bottom) displays a somewhat greater discrepancy between the models. The She-Leveque-Dubrulle model agrees best with data for the data sets with the lowest while the stable model agrees best for the other data sets. The normal inverse Gaussian model slightly overpredicts the self-scaling exponents, but in general it provides a prediction of comparable accuracy.
The normal inverse Gaussian model can be forced to align better with the stable model and the She-Leveque-Dubrulle model by requiring exponential moments of order to exist. This, however, compromises the prediction of the self-scaling exponents of the two-point correlators (subsec. 4.4).
4.3 Two-point correlators: scaling
Since the sample moments of the velocity derivative become increasingly unreliable with increasing order, we consider only the orders for which and are positive half-integers with . This leads to 36 non-trivial combinations in the analysis of self-scaling of the correlators. The self-scaling exponents satisfy
and it is therefore sufficient to consider eight combinations.
Figure 8 shows that the correlator for data set no. 7 exhibits scaling for . The figure also shows that the correlator (18) determined by (15) provides a very good fit to the data for lags . The parameter determines the behaviour of the correlator for lags close to the decorrelation time . At lags , the correlators are corrupted by surrogacy effects . Lags have therefore been excluded from the fit. Similar results hold for the other data sets and other orders of the correlators.
The parameter , which determines the extend of the scaling range of the correlators, cannot be identified since surrogacy effects corrupt the correlator at small lags. We have chosen as it allows for clear scaling at all lags for all thirteen data sets.
where is a factor which depends only on and . Since cannot be reliably identified by the data available, (33) does not allow us to predict the correlator scaling exponents from the distribution of the energy dissipation alone.
The scaling exponents may still be estimated directly by fitting (15) and (18) to the correlators. Table 1 shows the estimated intermittency exponent where it is also compared to the estimate found in . The two estimates are in reasonable agreement. The differences are likely due to the differences in the estimation procedures. We note that for all the data sets, the intermittency exponent is approximately half of the prediction under the She-Leveque-Dubrulle model where .
As a final consistency check, we note that, according to (19), the scaling exponents satisfy a relation of the form
By (26), we may choose . The other unknowns are found by solving the linear least squares system corresponding to (34). The root mean square of the relative error is below except for data sets no. 1, 2, and 5 where it is , , and , respectively. We conclude that (34) is fulfilled to a satisfying degree.
4.4 Two-point correlators: self-scaling
The self-scaling of the correlators is confirmed by fig. 9. Under model (2), the self-scaling exponent is given by the ratio . It depends only on the shape parameters of the distribution of and on the order of the involved correlators. In particular, the self-scaling exponents are independent of the ambit set and the location and scale parameters of . Therefore, the self-scaling exponents are predicted directly from the estimated shape parameters listed in tab. 1, and they inherit the -independence from .
Figure 10 shows the estimated self-scaling exponents for all the data sets. In most cases the predictions of the normal inverse Gaussian model and the stable model are within the variation of the data point, though the normal inverse Gaussian model tends to predict larger values than the stable model. Both models are clearly superior to the normal model. The stable model appears to be slightly more accurate than the normal inverse Gaussian model. The case is exceptional as all models fail to predict the observed self-scaling exponents by approximately .
At low Reynolds numbers, the self-scaling exponents tend to deviate from the self-scaling exponents at higher Reynolds numbers. At low Reynolds numbers the inertial range is shorter. Whether this affects the estimation of the self-scaling exponents remains to be investigated.
Further observations may be drawn from fig. 10. Firstly, the predicted self-scaling exponents are in all three cases almost independent of the Taylor micro-scale Reynolds number. This follows from the fact that the estimated shape parameters of the distribution of are almost independent of (tab. 1). As noted previously, this independence hints at a universal character of the self-scaling exponents and the distribution of . Secondly, the shape parameters are derived from one-point statistics, while the self-scaling exponents are derived from two-point statistics. The predictability of the latter from the former hints at a parsimonious description of the energy dissipation as the exponential of an integral with respect to a Lévy basis: Higher order statistics are reasonably well predicted from lower order statistics.
5 Conclusion and outlook
Thirteen time series of one-point measurements of the velocity component in the mean stream direction in a helium jet are analysed from the point of view of the (surrogate) energy dissipation. The distribution of the logarithm of the energy dissipation is shown to be well approximated by normal inverse Gaussian distributions, a property possessed by neither the normal distribution nor the stable distribution. Furthermore, the shape of the distribution is apparently universal, i.e., it does not depend on the Reynolds number, at least within the considered data sets.
The two-point correlators of the energy dissipation show scaling and self-scaling. By modelling the energy dissipation as the exponential of an integral with respect to a Lévy basis, a connection is established between the shape of the distribution of the logarithm of the energy dissipation and the scaling and self-scaling exponents of the two-point correlators. Thus the model is parsimonious: the two-point correlators are all, to good accuracy, determined by the one-point distribution. The normal inverse Gaussian distributions are compared to the stable distributions. While the stable distributions are also capable of modelling the scaling and self-scaling exponents, they are far from capable of accurately modelling the distribution of the energy dissipation (fig. 5). In particular, the stable model implies infinite variance of , yet the data suggest that has finite moments at least up to order 6. The use of normal inverse Gaussian distributions allows accurate modelling of both the correlators and the distribution of the energy dissipation.
The model also allows prediction of the self-scaling exponents of the moments of the coarse-grained energy dissipation. Comparison to the predictions of the She-Leveque-Dubrulle model is made.
While the normal inverse Gaussian distributions allow fairly accurate modelling of the one-point distribution, two-point correlators of the energy dissipation, and the self-scaling exponents of the moments of the coarse-grained energy dissipation, a slight disadvantage related to exponential moments is identified. We conjecture the existence of a more appropriate infinitely divisible distribution without this disadvantage.
The energy dissipation exhibits stylised features beyond the scaling and self-scaling of the two-point correlators. In , breakdown coefficients and Kramers-Moyal coefficients are employed, in particular, to evaluate the use of stable distributions in the modelling. A similar analysis remains to be performed in the case of the normal inverse Gaussian distributions.
Appendix A The normal inverse Gaussian distribution
The normal inverse Gaussian distributions form a four-parameter family of probability distributions on the real line. They are a special case of the generalised hyperbolic distributions introduced in  to describe the law of the logarithm of the size of sand particles (see also [3, 2, 6, 13, 22]). The generalised hyperbolic distributions are applied in many areas of science, see e.g. [35, 7] and the references therein.
The probability density function of a normal inverse Gaussian distribution is given by
where , , and denotes the modified Bessel function of the second kind with index . The real parameter determines the location, and the positive parameter determines the scale. The parameters and are shape parameters and lie within the shape cone: .
From the asymptotic property as it follows that the NIG distribution has semi-heavy tails, specifically
as . This illuminates the role of and in determining the tails of the distribution.
The cumulant function of a random variable with distribution is given by
and the radius of convergence for the cumulant function is .
By differentiating (37) we obtain the following expressions for the first four cumulants,
where . Hence, the standardised third and fourth cumulants are
Equations (38) and (39) further illuminate the roles of the four parameters: and determine location and scale, respectively; is related to the skewness, specifically the tail asymmetry (if , the distribution is symmetric); and is related to the kurtosis.
It follows immediately from (37) that if are independent normal inverse Gaussian variables with common parameters and , but with individual location and scale parameters and (), then the distribution of the sum is where and . Therefore, the normal inverse Gaussian distributions are infinitely divisible, see also .
It is often desirable to describe the NIG distributions in terms of location-scale invariant parameters. By letting and we have that and are invariant under change of location and scale. For the purpose of interpreting the shape parameters it is sometimes advantageous to express the shape in terms of the location-scale invariant steepness and asymmetry , defined by
where , and is the alternative asymmetry. These parameters are within the shape triangle, defined by
Figure 11 shows the shape of the normal inverse Gaussian distributions for various values of the asymmetry and steepness . A wide range of shapes is possible. (See  for details on the shape of the family of generalised hyperbolic distributions).
A useful property of the normal inverse Gaussian distribution is the representation in terms of a mean-variance mixture of a normal distribution with the mixing distribution being an inverse Gaussian distribution, hence the name. Specifically, if is normal with mean and variance , and if is endowed with an independent inverse Gaussian distribution with parameters and , then follows a normal inverse Gaussian distribution with parameters . In short we may write
where is a normal distribution with zero mean and unit variance. For reference, the probability density function of an inverse Gaussian distribution with parameters and is
for . Figure 12 shows the logarithm of the probability density functions of the inverse Gaussian distribution corresponding to the three symmetric normal inverse Gaussian distributions in fig. 11. As the steepness increases, the probability that the random variance in (42) will attain large values increases.
To estimate the normal inverse Gaussian parameters from data one may apply maximum (or pseudo) likelihood methods. The maximum likelihood estimation is non-trivial since the likelihood function is very flat near the optimum. The computer program “hyp”  implements numerical maximisation of the likelihood function, but, in general, non-linear optimisation algorithms may also be applied. The approach of expectation-maximisation developed in  may be used in conjunction with other optimisation algorithms.
Appendix B Integration with respect to Lévy bases
The stochastic processes to be considered in the present paper are expressed in terms of integrals of deterministic functions with respect to Lévy bases on . A Lévy basis on is an infinitely divisible, independently scattered random measure, i.e., to each bounded Borel subset of , an infinitely divisible random variable is associated, the random variables associated to disjoint subsets are independent, and the random variable associated to a disjoint union is almost surely equal to the sum of the random variables associated to each subset, provided the union is a bounded Borel subset. For more details and mathematical rigour we refer to .
The stochastic integral of a deterministic function with respect to a Lévy basis is defined in three steps. Firstly, for an indicator function we define . Secondly, by requiring that the integral is linear in the integrand, the integral is extended to simple functions, i.e., linear combinations of indicator functions. Finally, since a measurable function may be approximated by a sequence of simple functions, the integral is defined to be the limit in probability of the sequence of integrals of the simple functions, provided this limit exists.
An important class of Lévy bases has the property that, informally, the distribution of the random variable associated to a subset does not depend on the location of the subset. In this case, we have the following fundamental representations of the Laplace transform and the characteristic function of the integral of a deterministic function with respect to a Lévy basis . Let and denote the logarithm of the Laplace transform and the characteristic function of the random variable , respectively. Then
where is a random variable (called the Lévy seed) whose cumulant function is related to the Lévy basis by
see  for more details. It follows that the distribution of the stochastic integral is determined by the function and the log-characteristic function of the Lévy seed .
The authors wish to thank B. Chabaud for granting permission to use the data. The authors furthermore acknowledge useful comments from Ole E. Barndorff-Nielsen and M. Greiner.
-  A. Arneodo, E. Bacry, and J. F. Muzy. Random cascades on wavelet dyadic trees. J. Math. Phys., 39:4142–4164, 1998.
-  O E Barndorff-Nielsen. Hyperbolic distributions and distributions on hyperbolae. Scand. J. Statist., 5(3):151–157, 1978.
-  O E Barndorff-Nielsen and Christian Halgreen. Infinite divisibility of the hyperbolic and generalized inverse Gaussian distributions. Z. Wahrscheinlichkeitstheorie verw. Gebiete, 38(4):309–311, 1977.
-  O E Barndorff-Nielsen and J Schmiegel. Lévy-based spatial-temporal modelling, with applications to turbulence. Russian Math. Surveys, 59(1):65, 2004.
-  Ole E Barndorff-Nielsen. Exponentially decreasing distributions for the logarithm of particle size. Proc. R. Soc. Lond. A, 353(1674):401–419, 1977.
-  Ole E Barndorff-Nielsen and P Blæsild. Hyperbolic distributions and ramifications: Contributions to theory and application. In C Taillie, G Patil, and B Baldessari, editors, Statistical Distributions in Scientific Work, volume 4, pages 19–44. Reidel, Dordrecht, 1981.
-  Ole E Barndorff-Nielsen, P Blæsild, and Jürgen Schmiegel. A parsimonious and universal description of turbulent velocity increments. Eur. Phys. J. B, 2004.
-  Ole E Barndorff-Nielsen, John Kent, and M Sørensen. Normal variance-mean mixtures and distributions. Int. Stat. Rev., 50(2):145–159, 1982.
-  R Benzi, L. Biferale, A. Crisanti, G. Paladin, M. Vergassola, and A. Vulpiani. A random process for the construction of multiaffine fields. Physica D, 65:352–358, 1993.
-  R Benzi, S Ciliberto, C Baudet, G.R. Chavarria, and R Tripiccione. Extended self-similarity in the dissipation range of fully developed turbulence. Europhys. Lett., 24(4):275–279, 1993.
-  R Benzi, G Paladin, G Parisi, and A Vulpiani. On the multifractal nature of fully developed turbulence and chaotic systems. J. Phys. A, 17(18):3521, 1984.
-  L. Biferale, G. Boffetta, A. Celani, A. Crisanti, and A. Vulpiani. Mimicking a turbulent signal: Sequential multiaffine processes. Phys. Rev. E, 57(6):R6261–R6264, 1998.
-  P Blæsild. The shape cone of the d-dimensional hyperbolic distribution. Technical report, Dep. of Theoretical Statistics, Inst. of Mathematics, Univ. of Aarhus, 1990.
-  P Blæsild and M K Sørensen. Hyp: A Computer Program for Analyzing Data by Means of the Hyperbolic Distribution. Technical report, Dep. of Theoretical Statistics, Inst. of Mathematics, Univ. of Aarhus, 1992.
-  O. Chanal, B. Chabaud, B Castaing, and B. Hébral. Intermittency in a turbulent low temperature gaseous helium jet. Eur. Phys. J. B, 17(2):309–317, 2000.
-  J. Cleve, M. Greiner, and K. R. Sreenivasan. On the effects of surrogacy of energy dissipation in determining the intermittency exponent in fully developed turbulence. Europhys. Lett., 61(6):756, 2003.
-  Jochen Cleve, Thomas Dziekan, Jürgen Schmiegel, Ole E Barndorff-Nielsen, Bruce R Pearson, Katepalli R Sreenivasan, and Martin Greiner. Finite-size scaling of two-point statistics and the turbulent energy cascade generators. Phys. Rev. E, 71(2):1–12, 2005.
-  Jochen Cleve and Martin Greiner. The markovian metamorphosis of a simple turbulent cascade model. Phys. Lett. A, 273:104 – 108, 2000.
-  Jochen Cleve, Martin Greiner, Bruce R Pearson, and Katepalli R Sreenivasan. Intermittency exponent of the turbulent energy cascade. Phys. Rev. E, 69(6):1–6, 2004.
-  Jochen Cleve, Jürgen Schmiegel, and Martin Greiner. Apparent scale correlations in a random multifractal process. Eur. Phys. J. B, 63(1):109–116, 2008.
-  Bérengère Dubrulle. Intermittency in fully developed turbulence: Log-Poisson statistics and generalized scale covariance. Phys. Rev. Lett., 73(7):959–962, 1994.
-  Ernst Eberlein and Ernst August v. Hammerstein. Generalized hyperbolic and inverse gaussian distributions: limiting cases and approximation of processes. In R. C. Dalang, M. Dozzi, and F. Russo, editors, Seminar on Stochastic Analysis, Random Fields and Applications IV, volume 58 of Progress in Probability, pages 221–264. Birkhäuser Verlag, 2004.
-  F.W. Elliot, Jr., A.J. Majda, D.J. Horntrop, and R.M McLaughlin. Hierarchical Monte Carlo Methods for Fractal Random Fields. J. Stat. Phys., 81(3/4):717–736, 1995.
-  U Frisch. Turbulence. Cambridge Univ. Press, Cambridge, 1995.
-  Uriel Frisch, Pierre-Louis Sulem, and Mark Nelkin. A simple dynamical model of intermittent fully developed turbulence. J. Fluid Mech., 87(04):719–736, 1978.
-  B. Jouault, M. Greiner, and P. Lipa. Fix-point multiplier distributions in discrete turbulent cascade models. Physica D, 136:125 – 144, 2000.
-  Bruno Jouault, Peter Lipa, and Martin Greiner. Multiplier phenomenology in random multiplicative cascade processes. Phys. Rev. E, 59:2451–2454, 1999.
-  A. Juneja, D.P. Lathrop, K. R. Sreenivasan, and G. Stolovitzky. Synthetic turbulence. Phys. Rev. E, 49(6):5179–5194, 1994.
-  Dimitris Karlis. An EM type algorithm for maximum likelihood estimation of the normal-inverse Gaussian distribution. Statist. Probab. Lett., 57(1):43–52, 2002.
-  A N Kolmogorov. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. J. Fluid Mech, 13(1):82–85, 1962.
-  Benoit B. Mandelbrot. Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier. J. Fluid Mech., 62(02):331–358, 1974.
-  Charles Meneveau and K. R. Sreenivasan. The multifractal nature of turbulent energy dissipation. J. Fluid Mech., 224:429–484, 1991.
-  A M Oboukhov. Some specific features of atmospheric tubulence. J. Fluid Mech., 13(01):77–81, 1962.
-  Donald B Percival and Andrew T Walden. Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques. Cambridge University Press, 1993.
-  K Prause. The Generalized Hyperbolic Model: Estimation, Financial Derivatives and Risk Measures. PhD thesis, Albert-Ludwigs University, 1999.
-  Balram S. Rajput and Jan Rosiński. Spectral representations of infinitely divisible processes. Probab. Th. Rel. Fields, 82(3):451–487, 1989.
-  Gennady Samorodnitsky and Murad S. Taqqu. Stable non-Gaussian random processes. Chapman & Hall, 1994.
-  D Schertzer and S Lovejoy. Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes. J. Geophys. Res, 92(D8):9693–9714, 1987.
-  Jürgen Schmiegel. Self-scaling of turbulent energy dissipation correlators. Phys. Lett. A, 337(4-6):342–353, 2005.
-  Jürgen Schmiegel, Ole E Barndorff-Nielsen, and H Eggers. A class of spatio-temporal and causal stochastic processes, with application to multiscaling and multifractality. S. Afr. J. Sci., 101:513–519, 2005.
-  Jürgen Schmiegel, Jochen Cleve, Hans C. Eggers, Bruce R. Pearson, and Martin Greiner. Stochastic energy-cascade model for -dimensional fully developed turbulence. Phys. Lett. A, 320(4):247–253, 2004.
-  Zhen-Su She and Emmanuel Leveque. Universal scaling laws in fully developed turbulence. Phys. Rev. Lett., 72(3):336–339, 1994.
-  Zhen-Su She and Edward Waymire. Quantized Energy Cascade and Log-Poisson Statistics in Fully Developed Turbulence. Phys. Rev. Lett., 74(2):262–265, 1995.
-  K. R. Sreenivasan and R. Antonia. The phenomenology of small-scale turbulence. Annu. Rev. Fluid Mech., 29:435–472, 1997.
-  T. Vicsek and A.-L. Barabási. Multi-affine model for the velocity distribution in fully developed turbulent flows. J. Phys. A, 24:L845–L851, 1991.