Universality Classes of Fluctuation Dynamics in Hierarchical Complex Systems

TL;DR

This paper introduces a unified framework for modeling the short-term fluctuation statistics of multiscale complex systems, revealing two universality classes characterized by distinct tail behaviors in their distributions.

Contribution

It develops a hierarchical stochastic model that derives two fundamental classes of fluctuation dynamics, supported by empirical validation in turbulence and financial data.

Findings

Two universality classes identified: power law and stretched exponential tails.

The model's predictions align well with empirical data from turbulence and markets.

Distributions are expressed using Meijer G-functions for analytical clarity.

Abstract

A unified approach is proposed to describe the statistics of the short time dynamics of multiscale complex systems. The probability density function of the relevant time series (signal) is represented as a statistical superposition of a large time-scale distribution weighted by the distribution of certain internal variables that characterize the slowly changing background. The dynamics of the background is formulated as a hierarchical stochastic model whose form is derived from simple physical constraints, which in turn restrict the dynamics to only two possible classes. The probability distributions of both the signal and the background have simple representations in terms of Meijer G-functions. The two universality classes for the background dynamics manifest themselves in the signal distribution as two types of tails: power law and stretched exponential, respectively. A detailed analysis of empirical data from classical turbulence and financial markets shows excellent agreement with the theory.

Equations

$$P(x|\epsilon )=\frac{1}{\sqrt{2\pi \epsilon }}\exp \left(-\frac{x^2}{2\epsilon }\right).$$

$$P(x)={\int }_0^{\infty }P(x|\epsilon )f(\epsilon )d\epsilon ,$$

$$d{\epsilon }_i=F({\epsilon }_i,{\epsilon }_{i-1})dt+G({\epsilon }_i,{\epsilon }_{i-1})d{W}_i,$$

$$d{\epsilon }_i=-{\gamma }_i\left({\epsilon }_i-{\epsilon }_{i-1}\right)dt+{\kappa }_i{\epsilon }_i^s{\epsilon }_{i-1}^{1-s}d{W}_i,i=1,...,N,$$

$$f({\epsilon }_i|{\epsilon }_{i-1})=\frac{{({\beta }_i{\epsilon }_{i-1})}^{{\beta }_i+1}}{\Gamma ({\beta }_i+1)}{\epsilon }_i^{-{\beta }_i-2}{e}^{-{\beta }_i{\epsilon }_{i-1}/{\epsilon }_i},$$

$$f({\epsilon }_i|{\epsilon }_{i-1})=\frac{{({\beta }_i/{\epsilon }_{i-1})}^{{\beta }_i}}{\Gamma ({\beta }_i)}{\epsilon }_i^{{\beta }_i-1}{e}^{-{\beta }_i{\epsilon }_i/{\epsilon }_{i-1}},$$

$$f_N({\epsilon }_N)=\int d{\epsilon }_{N-1}\dots \int d{\epsilon }_{1}f({\epsilon }_N|{\epsilon }_{N-1})\dots f({\epsilon }_1|{\epsilon }_0),$$

$$f_{N}(\varepsilon_{N})=\frac{1}{\varepsilon_{0}\omega\Gamma(\boldsymbol{\beta}+{\bf 1})}G_{N,0}^{0,N}\left(\begin{array}[]{c}{-\boldsymbol{\beta}-{\bf 1}}\\ -\end{array}\bigg{|}\frac{\varepsilon_{N}}{\varepsilon_{0}\omega}\right),$$

$$f_{N}(\varepsilon_{N})=\frac{\omega}{\varepsilon_{0}\Gamma(\boldsymbol{\beta})}G_{0,N}^{N,0}\left(\begin{array}[]{c}{-}\\ {\boldsymbol{\beta}-{\bf 1}}\end{array}\bigg{|}\frac{\omega\varepsilon_{N}}{\varepsilon_{0}}\right).$$

$$P_N(x)=\frac{1}{\sqrt{2\pi }}{\int }_0^{\infty }\exp \left(-\frac{x^2}{2{\epsilon }_N}\right){\epsilon }_N^{-1/2}{f}_N({\epsilon }_N)d{\epsilon }_N,$$

$$P_{N}(x)=\frac{1}{\sqrt{2\pi\omega\varepsilon_{0}}\Gamma(\boldsymbol{\beta}+{\bf 1})}G_{N,1}^{1,N}\left(\begin{array}[]{c}{-\boldsymbol{\beta}-{\bf 1/2}}\\ 0\end{array}\bigg{|}\frac{x^{2}}{2\omega\varepsilon_{0}}\right).$$

$$P_N(x)\sim \sum _{i=1}^N\frac{c_i}{x^{2{\beta }_i+3}},\text{mbox}for|x|\to \infty ,$$

$$P_{N}(x)=\frac{\omega^{1/2}}{\sqrt{2\pi\varepsilon_{0}}\Gamma(\boldsymbol{\beta})}G_{0,N+1}^{N+1,0}\left(\begin{array}[]{c}-\\ {\boldsymbol{\beta}-{\bf 1/2}},0\end{array}\bigg{|}\frac{\omega x^{2}}{2\varepsilon_{0}}\right).$$

$$P_N(x)\sim x^{2\theta }\exp \left[-(N+1)(\omega x^2/2{\epsilon }_0{)}^{1/(N+1)}\right],$$

$$d{\epsilon }_1(t)={\gamma }_1{\epsilon }_0>0,$$

$$y=\frac{{\epsilon }_N}{{\epsilon }_0}=\prod _{j=1}^N{\xi }_j,$$

$$g(y)={\int }_0^{\infty }\prod _{j=1}^N{g}_j({\xi }_j)d{\xi }_j\delta (y-{\xi }_1\cdots {\xi }_N),$$

$$g_j({\xi }_j)=\frac{{\beta }_j^{{\beta }_j+1}}{\Gamma ({\beta }_j+1)}{\xi }_j^{-{\beta }_j-2}{e}^{-{\beta }_j/{\xi }_j},$$

$$g_j({\xi }_j)=\frac{{\beta }_j^{{\beta }_j}}{\Gamma ({\beta }_j)}{\xi }_j^{{\beta }_j-1}{e}^{-{\beta }_j{\xi }_j}.$$

$$\mathcal{M}[g;s]\equiv {\int }_0^{\infty }dy{y}^{s-1}g(y).$$

$$\mathcal{M}[g;s]=\prod _{j=1}^N\mathcal{M}[g_j;s],$$

$$\mathcal{M}[g_j;s]={\beta }_j^{s-1}\frac{\Gamma ({\beta }_j-s+2)}{\Gamma ({\beta }_j+1)}$$

$$\mathcal{M}[g_j;s]=\frac{\Gamma ({\beta }_j+s-1)}{{\beta }_j^{s-1}\Gamma ({\beta }_j)}$$

$$\mathcal{M}[g;s]=\frac{{\alpha }^{-s}{\prod }_{j=1}^m\Gamma (s+b_j){\prod }_{j=1}^n\Gamma (1-s+a_j)}{{\prod }_{j=m+1}^q\Gamma (1-s+b_j){\prod }_{j=n+1}^p\Gamma (s+a_j)}$$

$$g(y)=G_{p,q}^{m,n}\left(\begin{array}[]{l}a_{1},\ldots,a_{p}\\ b_{1},\ldots,b_{q}\end{array}\bigg{|}\alpha y\right).$$

Full text

I Introduction

Complex phenomena are known to exhibit stationary time series which often show large deviations from Gaussian statistics. A common procedure to produce non-Gaussian tails is to allow for the violation of one of the assumptions of the central limit theorem, usually the condition that the increments have finite variance which leads to stable Lévy distributions levy . The physically undesirable consequence that such distributions have infinite variance is usually dealt with, a posteriori, by imposing some sort of truncation, generating the so-called truncated Lévy distributions truncatedlevy . In other cases, heavy-tailed distributions can be accounted for by a superposition of two statistics—a procedure known in mathematics as compounding compounding and in physics as superstatistics beck . Notwithstanding their empirical successes in several areas appliedlevy ; andrews_1989 ; castaing_PhysD90 ; chabaud_PRL94 ; beck , such ad hoc treatments of non-Gaussian effects are not completely satisfactory in that they offer no clear physical explanation for the underlying dynamics driving the fluctuation phenomena. Here we propose a general stochastic dynamical framework whereby heavy-tailed distributions naturally emerge. It is shown in particular that only two classes of distributions are allowed as characterized by the nature of the tails: power law and stretched exponential, respectively.

In this work we also address another important problem that one often faces when dealing with fluctuation phenomena in complex systems, namely the fact the empirical data can almost equally well be fitted with different probability distributions, making it difficult to select between competing models beck-cohen-swinney ; sornette_review . We tackle this model selection problem by introducing a joint fitting procedure whereby we simultaneously fit the empirical distributions of both the measured quantity (signal) and an internal variable (background), which represents the changing local environment. In our formalism, the background is described by a hierarchical stochastic model whose form is severely constrained by basic physical requirements, yielding only two ‘universality classes’ of possible dynamics. Analytical expressions are obtained for the distributions of both the signal and the background in terms of certain special transcendental functions—the Meijer $G$-functions. The availability of closed form solutions within the same family of special functions not only offers a unified theoretical treatment to the problem, but more importantly it allows one to implement the joint fitting procedure described above, which in turn helps to discriminate between the allowed distributions and also provides a direct check on the assumptions of the model. Excellent agreement with the theory is found in a detailed analysis of turbulence data as well as of financial asset price fluctuations. Owing to its generality, our theory could have a wide range of applications across different areas and disciplines.

II The Multiscale Approach

Consider a multiscale complex system in a stationary state where the probability distribution of the relevant observable, denoted by $x$, depends on the spatial or temporal scale with which the system is observed. The large-scale distribution, i.e., measured at a time scale $\tau_{0}$ above which fluctuations in $x$ are essentially uncorrelated, is assumed to be known and given by $P(x|\varepsilon_{0})$, where $\varepsilon_{0}$ is a parameter that characterizes the large-scale stationary state (global equilibrium) of the system. To be specific, we shall assume that at large scales our system follows a Gaussian distribution with variance $\epsilon_{0}$. (A Gibbsian distribution with mean $\epsilon_{0}$ could also be chosen with similar results.) We suppose furthermore that the observed quantity $x$ has a much faster dynamics than that of its local environment. This means that over short time periods (during which the environment does not change appreciably) the system follows the same distribution $P(x|\varepsilon)$, but where $\varepsilon$ now characterizes the ‘local equilibrium’. Under these assumptions, we can write

[TABLE]

Here we assume that $\varepsilon$ is a fluctuating quantity but one that varies more slowly (in time and space) than $x$. Physically, the quantity $\varepsilon$ can be identified, for example, with the local energy flux in a turbulent flow SV1 or as a ‘local temperature’ in a hierarchical system in thermal equilibrium SV2 . If we sample the system at short intervals but over a long time span (comparable to $\tau_{0}$), then the statistics of $x$ will be described by the marginal distribution

[TABLE]

where $f(\varepsilon)$ is the stationary probability density function of the background variable $\varepsilon$.

One distinctive aspect of our approach is that we wish to determine the possible distributions for $f(\varepsilon)$ from general physical arguments, rather than prescribe it a priori as is normally done in the compound and superstatistics formalisms. Another important ingredient is that we seek a theory that incorporates multiple time scales—a common feature of complex systems—and that clearly exhibits the connection between the local equilibrium variable $\varepsilon$ and its large scale counterpart $\varepsilon_{0}$. To this end, the Salazar-Vasconcelos model SV1 ; SV2 recently introduced to describe intermittency in multiscale fluctuation phenonema is a natural starting point. Our approach is also akin in spirit to the Palmer-Stein-Abrahams-Anderson (PSAA) hierarchical model for relaxation in spin glass anderson .

III The Dynamical Model

Let us assume that the system has $N$ well-separated time scales $\tau_{i}$, $i=1,...,N$, in addition to the large scale $\tau_{0}$, and let us order them from smallest to largest, i.e., $\tau_{i}\ll\tau_{i-1}$. Our arguments will be presented in a general setting, but to fix the ideas one may think of the system as a particle diffusing in a slowly changing environment. The relevant observable $x$ (say, the particle velocity) is measured at the shortest time scale $\tau_{N}$, while the environment seen by the particle is described by a slower degree of freedom (related to the local diffusion coefficient) represented by a stochastic variable $\varepsilon_{N}$ with distribution $f_{N}(\varepsilon_{N})$. In order to find $f_{N}(\varepsilon_{N})$, we use a hierarchical dynamical model where at each level $i$ $(i=1,...,N$) of the hierarchy the corresponding variable $\varepsilon_{i}$ is described by the following stochastic differential equation (SDE):

[TABLE]

where $F(\varepsilon_{i},\varepsilon_{i-1})$ represents the deterministic driving term, $G(\varepsilon_{i},\varepsilon_{i-1})$ is the noise amplitude, and the $dW_{i}$’s are independent Wiener processes. That the functions $F$ and $G$ depend only on $\varepsilon_{i}$ and $\varepsilon_{i-1}$ encodes the hierarchical nature (local interactions) of the system. It expresses the fact that the physical constraints imposed on the system at the large scale are not directly felt at the small scales but rather are transferred down the hierarchy through the intervening scales.

The possible functional forms of $F(\varepsilon_{i},\varepsilon_{i-1})$ and $G(\varepsilon_{i},\varepsilon_{i-1})$ are severely constrained by three general physical requirements: (i) equilibrium condition, which states that $\left\langle\varepsilon_{i}(t)\right\rangle=\varepsilon_{0}$ for $t\to\infty$; (ii) invariance under change of scale, which requires that $F_{i}(\lambda\varepsilon_{i-1},\lambda\varepsilon_{i})=\lambda F_{i}(\varepsilon_{i-1},\varepsilon_{i})$ and $G_{i}(\lambda\varepsilon_{i-1},\lambda\varepsilon_{i})=\lambda G_{i}(\varepsilon_{i-1},\varepsilon_{i})$; and (iii) positivity of $\varepsilon_{i}$, meaning that ${\rm Prob}(\varepsilon_{i}(t)<0)=0,\;\forall t,\;{\rm if}\;\varepsilon_{j}(t=0)\geq 0,\forall j$, which in turn entails that $G_{i}(0,0)=0$.

The most general model that satisfies the three requirements above is of the form

[TABLE]

where $\gamma_{i}$ and $\kappa_{i}$ are positive constants, and $s$ is, in principle, an arbitrary real number $0\leq s\leq 1$. The positivity of $\varepsilon_{i}(t)$ is proved in Appendix A. The stationary distribution $f(\varepsilon_{i}|\varepsilon_{i-1})$ can be computed from the Fokker-Planck equation associated with the SDE (4), for $\varepsilon_{i-1}$ kept fixed. If we now impose the (physically reasonable) condition that the corresponding Fokker-Planck equation should have analytic coefficients, then this further restricts the model to only two possible cases: $s=1$ or $s=1/2$.

In the case $s=1$ (which was first treated in SV1 ), the stationary conditional distribution $f(\varepsilon_{i}|\varepsilon_{i-1})$ is given by an inverse-gamma distribution

[TABLE]

whereas for $s=1/2$ the result is a gamma distribution:

[TABLE]

where $\beta_{i}={2\gamma_{i}}/{k^{2}_{i}}$ in both cases.

III.1 Stationary Solutions

Explicit stationary solutions for $f(\varepsilon_{N})$ can be obtained for both models above in the regime of large separation of time scales, i.e., when $\gamma_{N}\gg\gamma_{N-1}\gg\ldots\gg\gamma_{1}$. In this case we can write

[TABLE]

where $f(\varepsilon_{i}|\varepsilon_{i-1})$ is given by either (5) or (6). A remarkable fact, not noticed before SV1 , is that this multiple integral can be carried out exactly for both cases, with the result being expressed in terms of the Meijer $G$-function meijer (a detailed derivation is shown in Appendix B). For the inverse-gamma case ($s=1$) one finds

[TABLE]

whereas for the gamma case ($s=1/2$) the integral yields

[TABLE]

Here $\omega=\prod_{j=1}^{N}\beta_{j}$ and we have introduced the vector notation ${\boldsymbol{\beta}}\equiv(\beta_{1},\dots,\beta_{N})$ and $\Gamma({\bf a})\equiv\prod_{j=1}^{N}\Gamma(a_{j}).$

It is perhaps worth noting at this point that the lognormal distribution naturally arises from our multiscale model in the limit of infinitely many scales. This follows from a direct application of the central limit theorem (CLT) to the variable $\ln\varepsilon_{N}$, where $\varepsilon_{N}=\xi_{N}\xi_{N-1}\cdots\xi_{1}$, with $\xi_{i}={\varepsilon_{i}}/{\varepsilon_{i-1}}$. To apply CLT, we need however to ensure that the second moment of $\varepsilon_{N}$ remains finite for $N\to\infty$, which is the case for both $s=1$ and $s=1/2$ if we take the limit $\beta_{i}=\beta\to\infty$ in such a way that $\sigma^{2}=N/\beta$ stays finite SV1 . We thus see that the lognormal model often used in the compound approach castaing_PhysD90 ; chabaud_PRL94 ; beck-cohen-swinney is a limiting case of our multiscale formalism when infinitely many scales are considered.

III.2 Marginal Distributions

The marginal distribution $P_{N}(x)$ of the signal $x$ measured at the time scale $\tau_{N}$ can now be computed from the superposition integral (2), with $P(x|\varepsilon)$ as shown in (1). More specifically, we have

[TABLE]

where $f_{N}(\varepsilon_{N})$ is as given in (8) or (9). It is again remarkable that this integral can be performed exactly for both cases, yielding the following two classes of distributions:

i) Power-law class. This corresponds to the case $s=1$. Upon substituting (8) into (10) and using known properties of the $G$-functions meijer , the resulting integral can be expressed also in terms of a $G$-function:

[TABLE]

This formula is an alternative representation of the generalized hypergeometric function ${{}_{N}}F_{0}$ reported in SV1 . We note, in particular, that the distributions in this family all have power-law tails SV1 :

[TABLE]

where the $c_{i}$’s are constants.

ii) Stretched exponential class. This class corresponds to $s=1/2$, in which case the superposition integral (10) computed with (9) yields

[TABLE]

This new class of distributions is a generalization of the $K$-distribution ($N=1$) which is a known compound distribution with applications to, e.g., scattering in random media radar and turbulence andrews_1989 ; guhr2015 . From the asymptotic expansion of the function $G_{m,n}^{p,q}$, with $m=q=0$ and $n=p=N+1$, one finds that the tail of the distribution in this case is given by a modified stretched exponential

[TABLE]

where $\theta=(\sum_{i=1}^{N}\beta_{i}-N)/(N+1)$. It is interesting to note that a stretched exponential also appears (as the long time asymptotics) in the PSAA model for relaxation in spin glass anderson .

III.3 The Joint Fitting Procedure.

As a first application of our theory we have analyzed two sets of data: i) velocity measurements in a turbulent fluid and ii) financial asset prices. Both systems display hierarchical structures and hence are natural candidates for applications of our formalism. In turbulence frisch one has an energy cascade from large to small length scales, so that our background variables $\epsilon_{i}$ correspond to the fluctuating energy fluxes between adjacent scales. In financial markets there is an information cascade from long to short temporal scales nature1996 , with the background variables representing the stochastic volatilities at the different time scales. Furthermore, the statistics of both velocity and price fluctuations have been modeled by different distributions (often with comparable degrees of agreement, see, e.g., sornette_review ; guhr2015 ; beck_2016 ), and so they are a good test for our joint fitting procedure (see below) whereby the distributions of both the signal and the background are fitted according to Eqs. (11) and (8), respectively. (Because the empirical distributions in both cases show a power-law trend it suffices to consider the case $s=1$.)

In the turbulence data we consider, our variable $x$ represents velocity increments in a low temperature gaseous helium jet french_group , i.e., $x(t)=u(t+\tau)-u(t)$, where $u$ is the velocity measured (by a hot wire probe) on the axis of the jet and $\tau$ is the inverse of the sampling rate. Here we analyze a data set of $10^{7}$ points obtained from an experimental run performed at Reynolds number $Re=295000$ and $\tau^{-1}=271.5$ kHz french_group . In Fig. 1(a) we show the symmetric part of the empirical distribution for the velocity increments together with the best fit by the theoretical formula $P_{N}(x)$, given in (11), for various values of $N$. In performing the fits we set $\beta_{i}=\beta$ and thus are left with only one parameter to adjust for any given $N$. One sees from Fig. 1(a) that there appears to be valid solutions for different values of $N$, and so an unambiguous choice of model distribution cannot be made on the basis of this comparison alone.

To solve this quandary, we need to extract the background series directly from the experimental data and examine its distribution. For that, we first divide our original time series in intervals of size $M$ and for each such interval compute a variance estimator, $\epsilon(t)=\frac{1}{M}\sum_{j=0}^{M-1}[x(t-j\delta t)-\bar{x}(t)]^{2}$, where $\bar{x}(t)=\frac{1}{M}\sum_{j=0}^{M-1}x(t-j\delta t)$, thus generating a new time series. Next, we numerically compound the background series $\epsilon(t)$ with a Gaussian for various $M$, and select the value of $M$ for which the corresponding superposition integral best fits the distribution of the original series. Excellent agreement is found for $M=19$; see inset of Fig. 1(b). In the main plot of Fig. 1(b) we compare the distribution of $\epsilon(t)$ for this optimal window size with the distribution $f_{N}(\varepsilon_{N})$ given in (8) for the same parameters as in Fig. 1(a). Taking into account both figures 1(a) and 1(b), we conclude that the solution with $N=7$ and $\beta=15.5$ gives the best overall fit to the turbulence data.

The financial data we analyzed correspond to the logarithmic intraday returns of the Ibovespa index of the São Paulo Stock Exchange, that is, $x(t)=\ln S(t+\tau)-\ln S(t)$, where $S(t)$ is the Ibovespa value at time $t$ and $\tau=30$s. We analyzed quotes from the period of November 2002 to March 2004, corresponding to a total of about 53000 points. The empirical distribution of the returns and respective fits using (11) are shown in Fig. 2(a). In this case, the optimal window size for the variance series is $M=5$, as indicated in the inset of Fig. 2(b), with the fits to the variance distribution being shown in the main plot of Fig. 2(b). Using the joint fitting procedure outlined above, we conclude that the best combined fit occurs for $N=3$ and $\beta=1.2$. It is interesting to note that for both sets of data considered here we have found $N>1$, thus indicating that fluctuation processes across multiple scales do indeed take place in these systems. (Recall that $N$ is an estimate of the number of relevant scales in the problem.)

IV Conclusions

Describing fluctuation phenomena in multiscale complex systems is an admittedly difficult task because (among other reasons) one does not usually have direct access to the interscale dynamics, and hence indirect inferences have to be made about its effect on the measured quantities. For instance, non-Gaussian statistics is usually seen as an evidence of complex interactions between scales, but a general dynamical framework to explain such deviations from Gaussianity has not yet been established. Here we have shown, from a rather minimal set of assumptions on the interscale stochastic dynamics, that there exist two general classes of heavy-tailed distributions for the statistics of multiscale fluctuations. The distributions in both classes are given in terms of the same family of special functions (Meijer $G$-function) but differ regarding the nature of the tail: power law and modified stretched exponential, respectively. Good agreement was found with experimental data on classical fluid turbulence as well as financial data—both sets of data analyzed here were shown to belong to the power-law class. Further development of the theory presented here and additional applications, including of the stretched-exponential class, will be discussed in forthcoming publications.

Acknowledgements.

We are grateful to B. Chabaud and P. E. Roche for sharing their turbulence data with us and to the São Paulo Stock Exchange for providing the financial data. This work was supported in part by the Brazilian agencies CNPq and FACEPE.

Appendix A Positivity of $\varepsilon_{i}(t)$

By hypothesis, the large scale variable $\varepsilon_{0}$ is positive and $\varepsilon_{1}(t)$ is a continuous stochastic process. Therefore, if a given realization of (4) were to give a negative value of $\varepsilon_{1}$ at some time $t^{\prime}$, then it would have reached the value zero at an earlier time $t<t^{\prime}$. Setting $\varepsilon_{1}(t)=0$ in (4), we get

[TABLE]

which implies that $\varepsilon_{1}>0$ for the whole process. Iterating this argument, we conclude that if $\varepsilon_{j-1}(t)>0$ then $\varepsilon_{j}(t)>0$ for all $t$. We close by remarking that the case $\varepsilon_{0}=0$ yields the trivial fixed point $\varepsilon_{j}(t)=0$, $\forall j,t$.

Appendix B Background Distributions

We start by introducing the variable

[TABLE]

where $\xi_{j}=\varepsilon_{j}/\varepsilon_{j-1}$, then $f_{N}(\varepsilon_{N})=g(y)/\varepsilon_{0}$ and

[TABLE]

where from (5) we get

[TABLE]

whilst from (6) we obtain

[TABLE]

Now we apply the Mellin transform, defined as

[TABLE]

to both sides of (17). We find

[TABLE]

where

[TABLE]

is the Mellin transform of (18) and

[TABLE]

is the Mellin transform of (19). Next, we use the following property of the Meijer $G$ function meijer . If the Mellin transform of $g(y)$ is

[TABLE]

then

[TABLE]

Using (21), (22), (23) and (24) we obtain (8) and (9) respectively.