Multifractal analysis of financial markets: a review

Consider a measure m embedded in a geometric support ${\mathcal{F}}$, whose density is $\dot{m}(t')$ at position $t'$. By definition, $\int_{t'\in{\mathcal{F}}} \dot{m}(t'){\rm d}t'=1$. The measure in the neighbourhood of $t'$ is $\dot{m}(t'){\rm d}t'$. Using the idea of the classic box-counting method [163], we cover the geometric support ${\mathcal{F}}$ using boxes of size s. The integrated measure $m(s,t)$ in the tth box ${\mathcal{B}}(s,t)$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m(s,t) = \int_{t'\in{\mathcal{B}}(s,t)} \dot{m}(t'){\rm d}t', \nonumber \end{align} \tag{ 14 }$$

where $\sum_t {\rm m}(s,t)=1$ for any s. The fractal dimension of ${\mathcal{F}}$ can be determined as follow

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_f := D_0= \lim_{s\to0}\frac{\ln\sum_t[m(s,t)]^0}{\ln{1/s}}, \nonumber \end{align} \tag{ 15 }$$

where the term $\sum_t[m(s,t)]^0$ gives the number of non-empty boxes needed to cover the support ${\mathcal{F}}$. Note that fractal dimension is also called similarity dimension or capacity dimension. It is obvious that, if there is no measure distributed in a box, the box should not be counted. In most cases in handling financial time series, we have D₀  =  1. In the empirical determination of D₀, we plot $\sum_t[m(s,t)]^0$ against s in log–log scales and the slope of the linear part in a proper scaling range is regarded as its estimate.

The information entropy or Shannon entropy of the measure can be defined as [182, 183]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I(s) = -\sum_t m(s,t)\ln[m(s,t)], \nonumber \end{align} \tag{ 16 }$$

and one can define the information dimension as follows

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma = \lim_{s\to0}\frac{I(s)}{\ln{(1/s)}}, \nonumber \end{align} \tag{ 17 }$$

which was introduced by Balatoni and Rényi in [184, 185] and used to study strange attractors [186, 187].

A third index used to characterize the measure is the correlation dimension [188, 189]. The correlation dimension was originally introduced to study the scaling behavior of the correlation integral of time series of length N:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C(s) = \lim_{N\to\infty}\frac{1}{N^2}\sum_{i,j=1}^N{\mathcal{H}}(s-|\vec {x}_i-\vec {x}_j|) \approx \sum_t[m(s,t)]^2, \nonumber \end{align} \tag{ 18 }$$

where ${\mathcal{H}}(x)$ is the Heaviside function. The correlation dimension is defined as follows

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu = \lim_{s\to0}\frac{\ln\sum_t[m(s,t)]^2}{\ln{1/s}}. \nonumber \end{align} \tag{ 19 }$$

These three dimensions D_f, $\sigma \ {\rm and} \ \nu$ can be unified into one framework of generalized dimensions [167–169]. We define the q-order Rényi entropy as [190, 191]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_q(s) = \lim_{p\to{q}}\frac{1}{1-p}\ln\sum_t[m(s,t)]^{\,p}, \label{Eq:RenyiEn} \nonumber \end{align} \tag{ 20 }$$

and the generalized dimension as [192]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_q &= \lim_{s\to0}\frac{I_q(s)}{\ln(1/s)} = \lim_{s\to0} \lim_{p\to{q}}\frac{1}{p-1}\frac{\ln{\chi(\,p,s)}}{\ln{s}}\nonumber \\ &=\left\{\begin{array}{@{}ll@{}} \lim_{s\to0} \frac{1}{q-1}\frac{\ln{\chi(q,s)}}{\ln{s}}, & q\neq1 \nonumber \\ \lim_{s\to0} \frac{\sum_t m(s,t)\ln[m(s,t)]}{\ln{s}}, & q=1 \end{array} \right. \label{Eq:MF-PF:Dq:Iq:s} \nonumber \end{align} \tag{ 21 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi(q,s)= \sum_t\left[m(s,t)\right]^q \label{Eq:MF-PF:chi:s} \nonumber \end{align} \tag{ 22 }$$

is the q-order partition function of the measure. It is easy to verify that $D_0=D_f$, $D_1=\sigma$ and $D_2=\nu$. In practice, we can compute D_q according to equation (21). For a given q, we plot I_q(s) against $\ln{s}$ for various box sizes s and perform linear regression in a proper scaling range to obtain the slope D_q.

We can rewrite equation (21) as follows

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (q-1)D_q = \lim_{s\to0} \frac{\ln{\chi(q,s)}}{\ln{s}} \label{Eq:MF-PF:Dq:chi} \nonumber \end{align} \tag{ 23 }$$

or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi(q,s) \sim s^{\tau(q)}, \label{Eq:MF-PF:tau:chi:s} \nonumber \end{align} \tag{ 24 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = (q-1)D_q. \label{Eq:MF-PF:tau:Dq} \nonumber \end{align} \tag{ 25 }$$

In practice, we can compute $\tau(q)$ according to equation (24). For a given q, we can compute $\chi(q,s)$ for various box sizes s and perform a linear regression of $\ln\chi(q,s)$ against $\ln{s}$ in a proper scaling range to obtain $\tau(q)$. Alternatively, $\tau(q)$ can be transformed from D_q using equation (25).

Recently, Xiong and Shang proposed a variance-weighted partition function approach and established the relationship between the corresponding multifractal functions [193]. Numerical simulations show that the variance-weighted partition function approach performs comparatively as the standard partition function approach.

Based on the box-counting idea, the geometric support is partitioned into non-overlapping boxes of size s. The local singularity strength $\alpha$ in the tth box ${\mathcal{B}}(s,t)$ is defined according to the following relationships [171]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m(s,t) \sim s^{\alpha}, \label{Eq:MF-PF:m:s:alpha} \nonumber \end{align} \tag{ 26 }$$

which may differ for different boxes. Let $N_s(\alpha)$ denote the number of boxes in which the singularity strengths are within $[\alpha,\alpha+d\alpha]$. Hence, the fractal dimension of the set is determined according to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle N_s(\alpha) = \rho(\alpha')s^{-f(\alpha')}{{\rm d}}\alpha', \label{Eq:MF-PF:Ns:s:f:alpha} \nonumber \end{align} \tag{ 27 }$$

where $\rho(\alpha)$ is the density of the singularity strength and $f(\alpha)$ is the fractal dimension of the boxes under consideration [171], which is usually called multifractal spectrum or singularity spectrum. One can directly obtain $\alpha$ from equation (26) and $f(\alpha)$ from equation (27) using local and pointwise singularity methods [194–196].

Inserting equation (26) into the partition function (22) and rewriting the sum into an integral over $\alpha$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi(q,s) &= \int{s}^{q\alpha'}\rho(\alpha')s^{-f(\alpha')}{{{\rm d}}}\alpha' \nonumber \\ &= \int \rho(\alpha')s^{q\alpha'-f(\alpha')}{{{\rm d}}}\alpha'. \label{Eq:MF-PF:chi:integral} \nonumber \end{align} \tag{ 28 }$$

Assume that $\rho(\alpha)$ is nonzero and non-singular. It therefore contributes a constant proportional factor to the leading behavior of the integral. Because s is small, according to the method of steepest descent, the integral will be dominated by the value of $\alpha(q)$ that makes the power $q\alpha'-f(\alpha')$ smallest. We thus replace $\alpha'$ by $\alpha(q)$, which is defined by the extremal condition:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle {{{{\rm d}}}\left[q\alpha'-f(\alpha')\right]}/{{{{\rm d}}}\alpha'}\big|_{\alpha'=\alpha(q)}=0 \label{Eq:MF-PF:MinCond1} \nonumber \end{align} \tag{ 29a }$$

and

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle {{{{\rm d}}}^2\left[q\alpha'-f(\alpha')\right]}/{{{{\rm d}}}\alpha'^2}\big|_{\alpha'=\alpha(q)}>0. \label{Eq:MF-PF:MinCond2} \nonumber \end{align} \tag{ 29b }$$

Then, equation (28) becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi(q,s) \sim s^{q\alpha-f(\alpha)}. \label{Eq:MF-PF:chi:s:qalpha:f} \nonumber \end{align} \tag{ 30 }$$

It follows from equation (29) that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f'(\alpha) = {{{{\rm d}}}f(\alpha)}/{{{{\rm d}}}\alpha} = q \label{Eq:MF-PF:df:dalpha:q} \nonumber \end{align} \tag{ 31a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f''(\alpha) < 0, \label{Eq:MF-PF:d2f:dalpha2} \nonumber \end{align} \tag{ 31b }$$

which means that the multifractal spectrum $f(\alpha)$ is a concave function and its slope at point $(\alpha(q),f(\alpha(q))$ is q.

Comparing equations (24) and (30), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = q\alpha - f(\alpha). \label{Eq:MF-PF:tau:alpha:f} \nonumber \end{align} \tag{ 32 }$$

Taking the derivative of equation (32) with respect to q, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{{{\rm d}}}\tau(q)}{{{{\rm d}}}q} = \alpha +q\frac{{{{\rm d}}}\alpha}{{{{\rm d}}}q}-\frac{{{{\rm d}}}f(\alpha)}{{{{\rm d}}}\alpha}\frac{{{{\rm d}}}\alpha}{{{{\rm d}}}q} = \alpha. \label{Eq:MF-X-PF:dtau:dq:alpha} \nonumber \end{align} \tag{ 33 }$$

Rewriting equations (32) and (33), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha = {{{{\rm d}}}\tau(q)}/{{{{\rm d}}}q} \label{Eq:MF-X-PF:alpha:dtau:dq} \nonumber \end{align} \tag{ 34a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) = q\alpha -\tau(q). \label{Eq:MF-PF:f:alpha} \nonumber \end{align} \tag{ 34b }$$

This shows that $f(\alpha)$ is the Legendre transform of $\tau(q)$ [161, 171]. Therefore, after obtaining $\tau(q)$, we can determine numerically $\alpha(q)$ using equation (34a) and $f(\alpha)$ using equation (34b). To qualify a multifractal measure, we can use either $(q,\tau)$ or $(\alpha,f)$, which are equivalent.

From the canonical perspective, we can obtain the $f(\alpha)$ function directly [197, 198]. We can define the canonical measures

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu(q,s,t) = \frac{[m(s,t)]^q}{\sum_t [m(s,t)]^q}. \label{Eq:MF-PF:mu:q} \nonumber \end{align} \tag{ 35 }$$

The singularity strength $\alpha(q)$ and its spectrum $f(q)$ can be computed by linear regressions in log–log scales using the following equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(q) = \lim_{s\to0} \frac{\sum_t \mu(q,s,t) \ln[m(s,t)]}{\ln{s}}, \label{Eq:MF-PF:alpha:mu:m} \nonumber \end{align} \tag{ 36a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle f(q) & \triangleq f(\alpha(q)) \nonumber \\ &= \lim_{s\to0}\frac{\sum_t \mu(q,s,t) \ln\left[\mu(q,s,t)\right]}{\ln{s}}. \label{Eq:MF-PF:fq:mu} \nonumber \end{align} \tag{ 36b }$$

To compute $\alpha(q)$, we can plot $\sum_t \mu(q,s,t) \ln{[m(s,t)]}$ against $\ln{s}$, determine the scaling range and obtain the slope by linear regression. The $f(q)$ function can be determined similarly based on equation (36b). Combining equations (24) and (36), we can verify equation (32) easily. The mass exponent function $\tau{(q)}$ can be obtained by using equation (32).

For financial volatilities, we can further define the inverse partition function [199]. We start with an illustration of the partition function approach using a high-frequency volatility time series.

Denote $\{I(t): t=1,\cdots, T\}$ the price time series of a financial asset. The logarithmic return $r(t)$ at the finest time resolution is calculated by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r(t) = \ln{I(t)} - \ln{I(t-1)}, \label{Eq:Def:return} \nonumber \end{align} \tag{ 37 }$$

where $t=2,\cdots,T$. The absolute return is utilized as a proxy for volatility according to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v(t) = |r(t)|. \label{Eq:Def:volatility} \nonumber \end{align} \tag{ 38 }$$

Figure 3(a) illustrates a segment of the 1 min volatility time series of the S&P 500 index, in which, according to the partition function method, the original volatility time series is divided into N boxes with identical size s  =  T/N. The box sizes are chosen such that the number of boxes of each size is an integer to cover the whole time series.

Zoom In Zoom Out Reset image size

We define the continuous volatility measure $u(t')$ for any $t'\in (t-1, t]$ by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u(t') = \frac{v(t)}{\sum_{t=1}^Tv(t)}, \label{Eq:continue:volatility} \nonumber \end{align} \tag{ 39 }$$

where $t=1,\cdots,T$. The cumulative function of the volatility measure is obtained as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U(t) = \int_0^t u(t'){\rm d}t'. \label{Eq:CumSum:measure} \nonumber \end{align} \tag{ 40 }$$

The function $U(t)$ of $v(t)$ in figure 3(a) is presented in figure 3(b). The measure in each box of size s is determined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_n(s) = U(ns)-U(ns-s), \label{Eq:direct:measure} \nonumber \end{align} \tag{ 41 }$$

which is also illustrated in figure 3(b). Usually, one computes the measure in a discrete way,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_n(s) = \sum_{t=(n-1)s+1}^{ns} u(t), \label{Eq:disctrete:measure} \nonumber \end{align} \tag{ 42 }$$

in which T/s must be an integer to avoid edge effects (see section 5.1.2 for more information). When the continuous measure $u(t)$ is used, more choices for s values are possible as N can be any integer.

For a given order q, the direct partition function $\chi_q(s)$ can be estimated by using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_q(s) = \sum_{n=1}^{N} [\mu_n(s)]^q. \label{Eq:dpf} \nonumber \end{align} \tag{ 43 }$$

In figure 3(d), we plot $[\chi_q(s)]^{1/q-1}$ as a function of s for different q's. The scaling exponent function $\tau(q)$ can be obtained through power-law regressions to the data points in the scaling range, which spans about three orders of magnitude. We find from figure 3(f) that $\tau(q)$ is a nonlinear function of q, confirming the presence of multifractality in the volatility measure.

We now turn to investigate the inverse partition function. For each threshold $\mu$, a sequence of exit times $\{s_j(\mu): j=1,\cdots,J\}$ can be determined successively by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{j=1}^J s_j = \inf \left\{t : U(t) \geqslant j \mu \right\}, \label{Eq:exittime} \nonumber \end{align} \tag{ 44 }$$

where $J = 1/\mu$ is an integer. Graphically, figure 3(c) shows how the exit times s_j can be determined. The inverse measure is defined as the normalized exit time

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mu^\dagger_j(\mu) = s_j/T, \label{Eq:imeasure} \nonumber \end{align} \tag{ 45 }$$

and the inverse partition function can be determined as follows:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \chi^\dagger(\mu,p) = \sum_{j=1}^J \left[\mu^\dagger_j(\mu)\right]^{\,p}. \label{Eq:MF-PF:iPF} \nonumber \end{align} \tag{ 46 }$$

Figure 3(e) shows the dependence of $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}[\chi^\dagger(\mu,p)]^{1/(\,p-1)}$ as a function of the thresholds $\mu$ for different p  values. Power-law scaling can be observed over about three orders of magnitude:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \chi^\dagger(\mu,p) \sim \mu^{\tau^\dagger(\,p)}. \label{Eq:MF-PF-chip2} \nonumber \end{align} \tag{ 47 }$$

The scaling exponent function $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\tau^\dagger(\,p)$ can be obtained through power-law regressions to the data points in the scaling range, as is shown in figure 3(f). The nonlinearity of the $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\tau^\dagger(\,p)$ function confirms the presence of multifractality in the exit time measure.

An important issue arises concerning the relationship between $\tau(q)$ and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\tau^\dagger(\,p)$. Mandelbrot and Riedi proved an elegant inversion formula mathematically for both discontinuous and continuous multifractal measures [200, 201]. Roux and Jensen independently proved the exact relation between the direct and inverse scaling exponents for the classical binomial measures [202]. Xu et al provided a simple 'proof' of the inversion formula for multinomial measures [203], which is presented below.

Let $\mu$ be a probability measure on $[0,1]$ with its integral function $M(t) = \mu([0,t])$. Then its inverse measure can be defined by

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mu^\dagger = M^\dagger(s) = \left\{\begin{array}{@{}lll@{}} \inf\{t:M(t)>s\}, && {{{\rm if~}}} s<1 \nonumber \\ 1, && {{{\rm if~}}} s=1 \end{array} \right., \label{Eq:MF-PF:invMu} \nonumber \end{align} \tag{ 48 }$$

where $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}M^\dagger(s)$ is the inverse function of $M(t)$. If $\mu$ is self-similar, then the relation $\mu = \sum_{i=1}^n p_i\mu(m_i^{-1}(\cdot))$ holds, where the m_i's are similarity maps with scale contraction ratios $r_i\in(0,1)$ and $\sum_{i=1}^np_i=1$ with p _i  >  0. The multifractal spectrum of measure $\mu$ is the Legendre transform $f(\alpha)$ of $\tau$, which is defined by the generating function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{i=1}^n p_i^qr_i^{-\tau}=1. \label{Eq:MF-PF:iPF:tau} \nonumber \end{align} \tag{ 49 }$$

It can be shown [200, 201] that the inverse measure $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\mu^\dagger$ is also self-similar with ratio $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}r_i^\dagger=p_i$ and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}p_i^\dagger=r_i$, whose multifractal spectrum $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}f^\dagger(\alpha^\dagger)$ is the Legendre transform of $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\tau^\dagger$, which is defined implicitly by

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \sum_{i=1}^n (\,p_i^\dagger)^{\,p} (r_i^\dagger)^{-\tau^\dagger}=1. \label{Eq:MF-PF:iPF:theta} \nonumber \end{align} \tag{ 50 }$$

It is easy to verify that the following inversion formula holds

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \left\{\begin{array}{@{}lll@{}} \tau(q) &=& -p \nonumber \\ \tau^\dagger(\,p) &=& -q \end{array} \right.. \label{Eq:MF-PF:Tau:Theta} \nonumber \end{align} \tag{ 51 }$$

Two equivalent testable formulae follow immediately

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \left\{\begin{array}{@{}lll@{}} \tau(q) = -\tau^{\dagger[-1]}(-q) \nonumber \\ \tau^\dagger(\,p) = -\tau^{[-1]}(-p) \end{array} \right., \label{Eq:MF-PF:InvForm} \nonumber \end{align} \tag{ 52 }$$

where the superscript ${[-1]}$ denotes the inverse function operator. Figure 3(f) verifies that the inversion formulae (52) hold for the high-frequency volatility time series of the S&P 500 index [199].

We can also relate the direct and inverse singularity strengths and the direct and inverse singularity spectra. From equation (51), we have

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \alpha(q)=\frac{{{{\rm d}}}\tau(q)}{{{{\rm d}}}q}=\frac{-{{{\rm d}}}p}{-{{{\rm d}}}\tau^\dagger(\,p)}=\frac{1}{\alpha^\dagger(\,p)}. \label{Eq:MF-PF:InvForm:alpha} \nonumber \end{align} \tag{ 53 }$$

Combining equations (51), (53) and the Legendre transform, we have

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle f(\alpha) &= q\alpha-\tau = -\tau^\dagger(\,p)\alpha+p = \alpha[-\tau^\dagger(\,p)+p\alpha^\dagger]\nonumber \\ &= f^\dagger(\alpha^\dagger)/\alpha^\dagger = \alpha f^\dagger(1/\alpha). \label{Eq:MF-PF:InvForm:falpha} \nonumber \end{align} \tag{ 54 }$$

Rigorous proofs of these results can also be found in [200, 201]. It is easy to obtain the relation between direct and inverse generalized dimensions [204]:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \left\{\begin{array}{@{}lll@{}} p &=& -(q-1)D(q) \nonumber \\ D^\dagger(\,p) &=& q/[1+(q-1)D(q)] \end{array} \right., \label{Eq:MF-PF:Dq:DP} \nonumber \end{align} \tag{ 55 }$$

where $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}D^\dagger(\,p)$ is the inverse generalized dimension function.

In the multifractal analysis of financial time series, almost all studies with very few exceptions focused on individual time series. One can further investigate an ensemble of financial time series of individual assets [205]. Ensemble averaging was developed to determine the multifractal dimensions of diffusion-limited aggregations [206–210], and unearthed a subtle discrete hierarchy quantified by complex fractal dimensions [211].

One defines the quenched and annealed mass exponents $\tau_\mathrm{quen} (q)$ and $\tau_\mathrm{ann} (q)$ as follows,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle \ln \chi_{q}(s) \right\rangle \sim -\tau_\mathrm{quen}(q) \ln s, \label{Eq:quenched} \nonumber \end{align} \tag{ 56 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln \left\langle \chi_{q}(s) \right\rangle \sim -\tau_\mathrm{ann}(q) \ln s, \label{Eq:Annealed} \nonumber \end{align} \tag{ 57 }$$

where the angular brackets $\langle \cdot \rangle$ signify the ensemble average over all time series. The annealed exponents are more sensitive to rare samples of the ensemble with unusual values of $\chi_q(s)$, while the quenched exponents are more characteristic of typical members of the ensemble [212].

The reasoning under the ensemble averaging lies in the consideration that one performs the measurement many times of the dynamics of an 'ensemble' of financial instruments. The ensemble averaging method provides an alternative way to measure market risks from individual equities other than the market index. It is also different from directly averaging the risk measures of many individual equities. In addition, although the method was initially developed for the partition function approach, its extension to other multifractal analysis methods is straightforward.

Another important statistical quantity of turbulent fields is the structure function of velocity increments [172, 173]. Classical experiments have been carried out to measure the structure function and its nonlinear scaling behaviors [174]. Such nonlinear scaling behavior is also termed as multifractality [161, 175]. The structure function approach has been also employed to investigate financial time series [18, 177, 213]. The qth-order structure function is also called the qth-order height-height correlation function in the literature [214, 215].

Consider a time series $\{X(i):i=1,\cdots, N\}$. The increments or innovations over time scale s are defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta{X}(i,s) = X(i)-X(i-s). \label{Eq:MF-SF:dX} \nonumber \end{align} \tag{ 58 }$$

For financial assets, $X(i)$ is the logarithmic price and $\Delta{X}$ is the return over the time interval of duration s. The qth-order structure function is defined as the qth-order moment of the increment distribution:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K(q,s) = \frac{\frac{1}{N-s}\sum_{i=s+1}^N\left|X(i)-X(i-s)\right|^q}{\frac{1}{N}\sum_i[X(i)]^q}. \label{Eq:MF-SF:Kqs} \nonumber \end{align} \tag{ 59 }$$

We stress that $q\geqslant0$ for time series in turbulence [174], finance [177, 216], and other fields, because $\Delta{X}(i,s)$ can be zero such that the moments of negative orders are not defined. When q  =  0, we have $\newcommand{\e}{{\rm e}} K(0,s)\equiv 1$, which is independent of the scale s. When q  =  2, $K(2,s)$ is proportional to the autocorrelation function $\langle{X(i)X(i-s)}\rangle$ [177, 216].

For self-similar time series, we expect to have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle K(q,s) \sim s^{\zeta(q)} \triangleq s^{qH(q)}, \label{Eq:MF-SF:Kq:s} \nonumber \end{align} \tag{ 60 }$$

where $H(q)$ is the generalized Hurst exponent. The scaling function is obtained as [217]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = \zeta(q)-1 = qH(q)-1. \label{Eq:MF-SF:tau:Hq} \nonumber \end{align} \tag{ 61 }$$

When q  =  0, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(0) = -1. \label{Eq:MF-SF:tau:q=0} \nonumber \end{align} \tag{ 62 }$$

In practice, for each fixed q, we calculate the $K(q,s)$ values with varying s values and perform a linear least-squares regression of $\ln{K(q,s)}$ against $\ln{s}$ in a properly chosen scaling range $[s_{\min},s_{\max}]$ to obtain $qH(q)$ and $\tau(q)$. It follows that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = \frac{\sum_s\ln{s}\sum_s\ln{K(q,s)}-N\sum_s\ln{s}\ln{K(q,s)}} {\left(\sum_s\ln{s}\right)^2-N\sum_s(\ln{s)}^2}-1. \label{Eq:MF-SF:tau:Hq:LSE} \nonumber \end{align} \tag{ 63 }$$

We can further determine the singularity function $\alpha(q)$ and the singularity spectrum $f(\alpha)$ through the Legendre transform. When $X(t)$ is monofractal, $H(q)=H$ is a constant. When $X(t)$ has a multifractal nature, $H(q)$ decreases with q.

According to the definition of $K(q,s)$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{{{\rm d}}}\ln{K(q,s)}}{{{{\rm d}}}q} = \frac{1}{K(q,s)}\frac{{{{\rm d}}}K(q,s)}{{{{\rm d}}}q} = \frac{\sum_tv^q\ln{v}}{\sum_t{v^q}}, \nonumber \end{align} \tag{ 64 }$$

where $\sum_t$ is the sum for $s+1\leqslant{t}\leqslant{T}$ with $v>0$. Denote k the proportionality coefficient of the scaling relation (60). It is obvious that k  >  0. Combining equations (60) and (61), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau'(q)k\ln{s} =\frac{{{{\rm d}}}\ln{K(q,s)}}{{{{\rm d}}}q} = \frac{\sum_tv^q\ln{v}}{\sum_tv^q}. \nonumber \end{align} \tag{ 65 }$$

Let $g(x)=x\ln{x}$. We have $g'(x)>0$ when $x\in[1/e,+\infty)$ and $g'(x)<0$ when $x\in(0,1/e]$.

Figure 4 presents the results of the multifractal analysis of the daily time series of the logarithmic price $X(t)$ of the Dow Jones Industrial Average index (DJIA, 1896–2015) based on the structure function approach.

Zoom In Zoom Out Reset image size

In a seminal study [213], Müller et al presented a statistical analysis of financial data documenting scaling laws in the form of equation (60) for q  =  1 and q  =  2, though they did not mention 'multifractal', 'multiscaling' or 'structure function'. Their analysis included intra-day prices over 3 years and daily prices over 15 years of four foreign exchange spot rates (DEM/USD, JPY/USD, CHF/USD, USD/GBP) and intra-day prices over 3 years of gold (XAU/USD). Similar analyses for the first- and second-order moments have been adopted by other economists [86].

The method of fluctuation analysis is a classical approach to extract the Hurst exponent of time series [218], which has been applied extensively in econophysics [219–223]. The fluctuation function can be calculated as follows,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z(2,s) = \sigma_2^2(s) = \left\langle\left[\Delta{X(i,s)} - \langle\Delta{X(i,s)}\rangle\right]^2 \right\rangle \sim s^ {2H} \label{Eq:MF-FA:Sigma:s} \nonumber \end{align} \tag{ 66 }$$

where H is the fluctuation analysis (FA) exponent, i.e. approximatively the Hurst index. A value of the Hurst exponent H larger than 0.5 means that the time series is correlated, For H  <  0.5, the time series is anti-correlated, and for H  =  0.5, it is uncorrelated.

One can extend the fluctuation analysis in equation (66) to higher orders as follows [224–226],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z(q,s) = \left\langle|\Delta{X(i,s)} - \langle\Delta{X(i,s)}\rangle|^q \right\rangle \sim s^{\zeta(q)} \label{Eq:MF-FA:Zq:s} \nonumber \end{align} \tag{ 67 }$$

which enables us to understand the multifractal nature in the dynamics of markets. The relationship between the scaling exponent $\zeta(q)$ and the generalized Hurst exponent $H(q)$ can be described by equation (61). When q  =  2, $H = H(2)$ is the Hurst exponent presented in equation (67). Using the Legendre transform, we can obtain the singularity strength $\alpha$ and its spectrum $f(\alpha)$.

Analogous to the inverse measure in the partition function approach [201, 227], Jensen [228] investigated the moments of exit time or first passage time. For a given increment $\Delta{X}$, using definition (58), the exit time scale s at time i is defined by

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle s(i, \Delta{X}) = \inf\{s^{\dagger}: \Delta{X}(i,s^{\dagger})\geqslant \Delta{X}\}, \nonumber \end{align} \tag{ 68 }$$

which is the minimal time needed for the fluctuation to exceed the threshold $\Delta{X}\geqslant0$. The first passage time has important usage in finance. For instance, Cho and Frees adopted it to construct an asymptotically unbiased volatility estimator [229]. They argued that the natural volatility estimator focuses on how much the price changes, whereas the 'temporal' estimator based on the first passage time focuses on how quickly the price changes [229]. The inverse statistics have also some connections to the persistence probability [230–233]. The variable $s(i, \Delta{X})$ and derived quantities have natural interpretations in finance, such the 'time to recovery' after a drawdown and are part of metrics used in a standard way in evaluating performance of financial strategies and hedge-funds. This explains that inverse statistics have been extensively applied in finance, as well as in other fields like in the study of heart beat rate and cardiac diseases [234].

If $X(i)$ is a Brownian motion, the distribution of s is solved analytically as a Gamma distribution [235, 236]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Pr(s) = \frac{1}{\sqrt{\pi}}\frac{a}{s^{3/2}}\exp\left(-\frac{a^2}{s}\right) \label{Eq:FPT:PDF:BM} \nonumber \end{align} \tag{ 69 }$$

where a is proportional to $\Delta{X}$. Since asset prices are not Brownian motions, Simonsen et al [237] proposed a generalized shifted Gamma distribution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Pr(s) = \frac{\nu}{\Gamma\left(\frac{\delta}{\nu}\right)}\frac{\beta^{2\delta}}{(s+s_0)^{\delta+1}}\exp\left[-\left(\frac{\beta^2}{s+s_0}\right)^\nu\right], \label{Eq:FPT:PDF:GS:Gamma} \nonumber \end{align} \tag{ 70 }$$

which has a power-law tail $\Pr(s)\sim s^{-(\delta+1)}$. The Gamma distribution in equation (69) is recovered if $\nu=1$, $\delta=1/2$, s₀  =  0 and $\beta=a$. They analyzed the wavelet filtered logarithmic daily closing prices of the DJIA (26 May 1896 to 5 June 2001) and found that the distribution (70) with $\delta\approx0.5$ fits the data well. This observation seems quite ubiquitous for uniformly spaced time series of indices and stocks in different markets [238–240], although there are studies reporting that $\delta$ can deviate from 0.5 [241, 242]. For tick-by-tick data at the transaction level, the tail exponent can be larger than 0.5, as found for five highly liquid stocks (AstraZeneca, GlaxoSmithKline, Lloyds TSB Group, Shell, and Vodafone) in the electronic market (SETS) of the London Stock Exchange (LSE) during the year 2002 [239] and the foreign exchange rate of DEM against USD for the full year of 1998 ($\delta=1.4$) [243].

Simonsen et al [237] argued that the most probable first passage time s^* is the optimal investment horizon and found that s^* scales with $\Delta{X}$ as a power law,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s^* \sim (\Delta{X})^\gamma, \label{Eq:FPT:OIH:PL} \nonumber \end{align} \tag{ 71 }$$

where $\gamma\approx1.8$, deviating from $\gamma=2$ for Brownian motions. It is found that $\gamma <2$ holds for regularly spaced samples in most markets and the $\gamma$ value is larger in developed stock markets than emerging markets [238, 244]. For tick-by-tick data that are unevenly sampled, it was reported that $\gamma>2$ [245].

For asset prices, it is natural to also consider the minimal time needed for the price to depreciate by a certain return $\Delta{X}\leqslant 0$. Then, the exit time s at time i is defined by [246]

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle s(i, \Delta{X}) = \inf\{s^{\dagger}: \Delta{X}(i,s^{\dagger})\leqslant \Delta{X}\}. \nonumber \end{align} \tag{ 72 }$$

A lot of research has shown that the optimal investment horizon is longer for gains than for losses, which is termed the gain-loss asymmetry [244, 246–254]. Several models have been proposed to reproduce the gain-loss asymmetry [255–257]. Nevertheless, there are also assets that do not show significant gain-loss asymmetry [241, 258].

An analysis in the frequency space based on the discrete wavelet transform showed that the gain-loss asymmetry is introduced mainly by the low-frequency content of the price series and the asymmetry disappears if enough of the low-frequency content is removed [257]. Alternatively, the gain-loss asymmetry is found to be caused by the non-Pearson-type autocorrelations in the time series [259]. Furthermore, the gain-loss asymmetry vanishes if the temporal dependence structure is destroyed by shuffling the time series [260].

The gain-loss asymmetry can be related to the leverage effect as follows. Recall first that the leverage effect is the fact that, after a large negative return, the volatility increases significantly and then relaxes back, usually roughly exponential in time [261–263]. In contrast, after a positive return, there is no significant change of volatility. And the causality is from (negative) returns to future volatility and not the reverse, in the sense that a variation of volatility does not produce any measurable change of the expected return in the future. The leverage effect can be interpreted as due to the impact of a loss on the risk perception of the firm, which increases as the equity over debt ratio has decreased as a result of the loss [261]. It also reflects a behavioral response of investors, who after a loss frantically reassess their risk exposure and readjust their portfolios, leading to large price moves. The leverage effect allows one to account for both (i) the asymmetry in the average time for a price drop of a fixed percentage versus a price gain of the same value (for instance 10 d for a loss of 5% and 20 d for a gain of 5% on the DJIA) and (ii) the fact that this asymmetry in waiting times for positive versus negative objectives is stronger for indices and weaker or absent for individual stocks [257, 260].

Consider a fixed drop target of $-5\%$. Such a drop is achieved by a succession of daily returns, more negative than positive. When some negative loss occurs on one day, the leverage effect leads to an increase of the amplitude of the next daily return. Because we are calculating a waiting time conditional to a cumulative drop of $-5\%$, the following returns will tend to be negative and with larger amplitude due to the leverage effect. In contrast, for the waiting time to reach a positive gain of $+5\%$, a majority of the daily moves will be positive, and the leverage effect will be weaker or absent if there are no or weak negative returns along the price path. As a consequence, as the amplitude of the daily returns tend to be small for a sequence of daily returns involved in a positive level ($+5\%$), the average waiting time to reach a positive return ($+5\%$) is larger than for a negative return goal. The asymmetry of a single stock is small because the leverage effect has a rather short memory and the difference between the number of up and down daily moves for reaching the target threshold is not large. In contrast, for an index, or equivalently a portfolio of N stocks, one must account for the ubiquitous existence of correlations between the returns of the stocks constituting the portfolio. These correlations, which are in general positive (most stocks move approximately together on average), produce an amplification of the gain-loss asymmetry, because the leverage effect is stronger by the averaging over the idiosyncratic residuals of the constituting stocks. This also leads to predict that the gain-loss asymmetry for a portfolio is the stronger, the stronger is the average correlation coefficients of its constituting stocks.

The moments of the exit time, called the distance structure functions [228] or inverse structure functions [264, 265], are defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_p(\Delta{X}) \equiv \langle{s^{\,p}}(\Delta{X})\rangle. \label{Eq:MF-SFA:ISF} \nonumber \end{align} \tag{ 73 }$$

Due to the duality between the structure function and the inverse structure function, one can intuitively expected that there is a power-law scaling stating that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_p({\Delta{X}})\sim {\Delta{X}}^{\phi(\,p)}, \nonumber \end{align} \tag{ 74 }$$

where $\phi(\,p)$ is a nonlinear concave function [228].

Synthetic data for the GOY shell model of turbulence showed perfect power-law dependence of the inverse structure functions on the velocity threshold [228]. For two-dimensional turbulence, the inverse structure functions exhibit well-defined multifractal properties [265]. In contrast, for three-dimensional turbulence, the inverse structure functions of an experimental time series at high Reynolds number do not exhibit clear power law scaling [264]. Nevertheless, different experiments show that the inverse structure functions of three-dimensional turbulence exhibit a more general scaling behavior called extended self-similarity (see next subsection) [266–268]. To our knowledge, the behavior of inverse structure functions of financial asset prices has not been studied.

The inversion formula relating the direct scaling exponents $\zeta(q)$ of direct structure functions and the inverse scaling exponents $\phi(\,p)$ of inverse structure functions can also be derived [269, 270] and give

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}lll@{}} \zeta(q) &=& -p \nonumber \\ \phi(\,p) &=& -q \end{array} \right.. \label{Eq:ZetaTheta} \nonumber \end{align} \tag{ 75 }$$

These expressions are verified by the simulated velocity fluctuations from the shell model [202]. However, this prediction (75) cannot be confirmed by wind-tunnel turbulence experiments (Reynolds numbers ${{{\rm Re}}}=400{{\rm \mbox{--}}}1100$) [267], which is not surprising because the inverse structure function does not scale as a power law of the velocity increment thresholds [264, 267].

To investigate the scaling properties of the inverse structure functions, one can define a set of relative exponents using the framework of extended self-similarity (ESS) [271]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K(q,s) \sim [K_{q_0}(s)]^{\xi(q,q_0)}, \label{Eq:MF-SF:ESS:Kq:Kq0} \nonumber \end{align} \tag{ 76 }$$

where q₀ is the order taken as a reference value. If equation (60) holds, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi(q,q_0) = \frac{qH(q)}{q_0H(q_0)}. \label{Eq:MF-SF:ESS:xi:Hq:Hq0} \nonumber \end{align} \tag{ 77 }$$

When the time series is monofractal, $H(q)$ is independent of q, that is, H(q)  =  H(q₀). It follows for monofractal time series that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi(q,q_0) = q/q_0. \label{Eq:MF-SF:ESS:xi:q:q0} \nonumber \end{align} \tag{ 78 }$$

In the case of velocity structure functions in fluid mechanics, q₀  =  3 is a natural choice based on the exact Kolmogorov's four-fifth law [268, 272]. For financial time series, we suggest to choose q₀  =  2 because $H(2)$ is the Hurst index. Generally, the ESS approach provides a wider scaling range for the extraction of scaling exponents.

Figure 5(a) presents the log–log plots of $K(q,s)$ against K₂(s) for $q= 1, 2, \cdots, 6$ with $s \in [1,2^{11}]$. The straight lines hold over more than two orders of magnitude and over at least 3.5 orders of magnitudes for $q\leqslant3$, showing the existence of extended self-similarity in the structure functions. The scaling range for small q's seems to be broader than for large q's. The ESS scaling exponents $\xi(q,2)$ are shown in figure 5(b). There is a weak indication that $\xi(q,2)$ has a nonlinear dependence as a function of q, with a downward curvature making the curve depart from the linear dependence $\xi(q,2) = q/2$ observes for small q's.

Zoom In Zoom Out Reset image size

Whether a time series possesses monofractal or multifractal behaviors is related to the innovation distributions. For simplicity of notation, we denote $v=|\Delta{X}(t,s)|$. By definition, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K(q,s) = \int_{0}^{\infty} v^q \Pr(s;v){{{\rm d}}}v = \left[\sigma(s)\right]^q \int_{0}^{\infty} {{{\rm d}}}x~\Phi(s;x) x^q, \label{Eq:MF-SF:ESS:PDF1} \nonumber \end{align} \tag{ 79 }$$

where $\Pr(s;v)$ is the distribution of $v$, $\sigma(s)$ is the standard deviation of $v$, $x=v/\sigma$ is the normalized volatility, and $\Phi(s;x)$ is the distribution of x. Eliminating $\sigma(s)$ using q  =  2, it follows that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K(q,s) = [K_2(s)]^{q/2}~ \frac{\int_{0}^{\infty} {{{\rm d}}}x~\Phi(s;x) ~x^q}{\left[\int_{0}^{\infty} {{{\rm d}}}x~\Phi(s;x) ~x^2 \right]^{q/2}}. \label{Eq:MF-SF:ESS:PDF2} \nonumber \end{align} \tag{ 80 }$$

If $\Phi(s;x)$ is universal and independent of s, then the last term in equation (80) is a number independent of s (and thus of K₂(s)) and the monofractal behavior $\xi(q,2) = q/2$ follows. Figure 5(c) shows the distributions of the normalized volatility $|\Delta{X}(t,s)|/\sigma(s)$ at three different scales. Although the three curves well overlap in their tails, they show clear differences in the bulk around x  =  1. This suggests that the ESS cannot be nonfractal, which is consistent with the presence of nonlinear curvature in figure 5(b).

We note that the ESS analysis is not constrained to the structure functions. Rather, it can be used for the moments of other multifractal analysis approaches. In validating and qualifying the multiscaling behaviour of the recurrence intervals of financial volatility, the relative dependence of recurrence interval moments of different orders of empirical time series has been investigated [273], which is in essence an ESS analysis [274].

The wavelet transform is widely used as a mathematical microscope to analyze time series [275, 276]. In particular, it can be used to analyze the singular structure of fractal and multifractal time series [277–279].

The wavelet transform of function $X(t)$ is [275, 276]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:WT:Def} C_{\psi}(b,s) = \frac{1}{s}\int_{-\infty}^{+\infty}X(t)\psi\left(\frac{t-b}{s}\right){{\rm d}}t, \nonumber \end{align} \tag{ 81 }$$

where $b\in{\mathbb{R}}$ is the position parameter, $s \in{\mathbb{R}}^+$ is the dilation parameter, $\psi(t)$ is the analyzing 'mother' wavelet. Gaussian wavelets are widely adopted [280] and are defined as derivatives of the Gaussian function:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_n(t) = g^{(n)}(t) = c_n\frac{{{\rm d}}^n}{{{\rm d}}t^n}{{\rm e}}^{-t^2/2}, \label{Eq:WT:GaussianWavelets} \nonumber \end{align} \tag{ 82 }$$

where c_n is the normalization coefficient. The second order (n  =  2) Gaussian wavelet is known as the Mexican hat.

To perform wavelet transform on a time series $\{X(i)\}_{i=1}^N$, one usually discretizes equation (81) in which b can take any value along the sampling spacing and s takes n_s scales forming a geometric sequence:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s_j=\lambda^{\,j-1}s_{\min}, ~~1\leqslant{j}\leqslant{n_s}, \label{Eq:WT:aj} \nonumber \end{align} \tag{ 83 }$$

where $\lambda>1$. The discretised wavelet transform of a given time series $Y(i)$ is then

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W(s,i) = \frac{1}{s} \sum_{j = 1}^N X(\,j) \psi \left(\frac{j-i}{s}\right),~i=1,\cdots, N. \label{Eq:MFXWT:WTTransform} \nonumber \end{align} \tag{ 84 }$$

The wavelet transform is used to identify different types of singularities in the signals in the time-scale plane, depending on the order n of derivative of the mother wavelet (82). More generally, choosing a mother wavelet $\psi$ satisfying $\int x^{m+1} \psi(s) {{\rm d}} x = 0$ implies that the wavelet transform is blind to polynomial dependence up to order m and only higher order powers are detected. For instance, for m  =  0, the wavelet transform is independent of the level of the signal and informs essentially on the local slope. For m  =  1, the wavelet transform is independent of the level of the signal and of its local slope, and it informs on the local curvature. We should also mention the existence of orthogonal wavelet transforms [280–282] and biorthogonal wavelet transform [283], which can be used directly to determine the Hurst index of a long-range correlated time series [284–286].

To investigate the multifractal nature of a time series, one can calculate the qth order moments of the wavelet coefficients [287, 288]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_q(s)=\frac{1}{N}\sum_{i=1}^N|W(s,i)|^q \label{Eq:WT:BO:Mq:a} \nonumber \end{align} \tag{ 85 }$$

and we expect that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_q(s) \sim s^{qH(q)}, ~s\gg0. \label{Eq:WT:BO:Mq} \nonumber \end{align} \tag{ 86 }$$

We note that M_q(s) might diverge for $q\leqslant0$ because some values of $|W(i,k)|$ might be close or equal to 0.

One can obtain the scaling exponents through ordinary least-squares regression. Or, we can adopt the weighted least-squares regression by putting more weight on the points with smaller variance [286]. Denoting

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:WT:BO:LqN} {\mathcal{L}}_{j,q}=\ln{M_q}/q, \nonumber \end{align} \tag{ 87 }$$

the corresponding weight w_j is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:WT:BO:wj} w_j \sim {{\rm Var}}({\mathcal{L}}_{j,q})^{-1} \nonumber \end{align} \tag{ 88 }$$

and the objective function is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:WT:BO:Obj} \sum_{j}w_j{\mathrm{E}}\left[({\mathcal{L}}_{j,q}-H(q)\ln{s_j}-A)^2\right]. \nonumber \end{align} \tag{ 89 }$$

The estimate of $H(q)$ is determined as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:WT:BO:hatH} \hat{H}(q)=\frac{\sum_{j}w_j\sum_{j}w_j\ln{s_j}{\mathcal{L}}_{j,q}-\sum_jw_j\ln{s_j}\sum_{j}w_j{\mathcal{L}}_{j,q}} {\sum_{j}w_j\sum_{j}w_j\ln^2{s_j}-\left(\sum_{j}w_j\ln{s_j}\right)^2}. \nonumber \end{align} \tag{ 90 }$$

A more elegant approach, which removes the redundancy contained in continuous wavelet transforms, is based on the wavelet transform modulus maxima (WTMM) [284, 288–290]. At each scale s_j, there exists several local maxima of |W(s_j, i)|. These maxima form several maxima lines in the $(i,s)$ plane. The WTMM is capable of detecting all the singularities of a large class of time series [291]. The traditional method based on partition functions may face problems in dealing with the situation of $q\leqslant0$. However, we will find that the multifractal analysis based on WTMM can overcome this difficulty because the wavelet coefficients on the WTMM line are local maxima that are non 0 [290].

Assume that there are N_j local maxima at scale s_j and $W(i_{j,k},s_j)$ is the kth local maximum located at i_j ,k. The qth order moments defined on the WTMM are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M(q,s_j)=\frac{1}{N_j}\sum_{k=1}^{N_j}|W(i_{j,k},s_j)|^q \label{Eq:MF-WTMM:Mq:a} \nonumber \end{align} \tag{ 91 }$$

or alternatively the partition functions read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF-WTMM:Zq:a} Z(q,s_j)=\sum_{k=1}^{N_j}|W(i_{j,k},s_j)|^q. \nonumber \end{align} \tag{ 92 }$$

For multifractal time series, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z(q,s)\sim s^{\tau(q)}. \label{Eq:MF-WTMM:tau:q} \nonumber \end{align} \tag{ 93 }$$

The generalized dimensions and singularity spectrum can then be calculated straightforwardly.

Analogous to the thermodynamic approximation in the direct determination of singularities and singularity spectrum in section 2.1.3, we can define the canonical measure [289],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_k(q,s_j)=\frac{|W(i_{j,k},s_j)|^q}{Z_q(s_j)}, \label{Eq:MF-WTMM:DC:mu} \nonumber \end{align} \tag{ 94 }$$

such that the singularity strength $\alpha$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(q)= \lim_{s_j\to0}\frac{1}{\ln{s_j}}\sum_{i=1}^{N_j}\left[\mu_k(q, s_j)\ln{|W(i_{j,k},s_j)|}\right] \label{Eq:MF-WTMM:DC:alpha} \nonumber \end{align} \tag{ 95 }$$

and the singularity spectrum $f(\alpha)$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) =\lim_{s_j\to0}\frac{1}{\ln{s_j}}\sum_{i=1}^{N_j}\left[\mu_k(q,s_j)\ln{\mu_k(q,s_j)}\right] \label{Eq:MF-WTMM:DC:f}. \nonumber \end{align} \tag{ 96 }$$

In applications, we use different scales s_j and perform (weighted) linear regression on the log–log scales.

Using respectively the partition function approach and wavelet analysis to detect the multifractal nature of a single time series, it is found that $\tau^{{{\rm PF}}}(q)=\tau^{{{\rm WT}}}(q)+q$ and $\alpha^{{{\rm PF}}}(q)=\alpha^{{{\rm WT}}}(q)+1$ [292–294].

A more recent multifractal formalism is based on wavelet leaders [295–300]. Wavelet leaders, introduced by Jaffard [295, 296], are defined based on discrete wavelet coefficients, which decompose the time series $X(i)$ on the bases $\{\psi_{j,k}\}_{j \in Z, k \in Z}$ composed of discrete wavelets $\psi_{j, k}$. The integers $j \in Z$ and $k \in Z$ define the scale s  =  2^j and location b  =  k2^j. The discrete wavelets $\{\psi_{j,k}\}_{j \in Z, k \in Z}$ are space-shifted and scale-dilated templates of an analyzing wavelet $\psi(t)$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi_{j,k}(t) = \frac{1}{2^{\,j}} \psi \left(\frac{i-k2^{\,j}}{2^{\,j}} \right), \label{Eq:MF-WL:OrthoBasis} \nonumber \end{align} \tag{ 97 }$$

which usually uses the Daubechies wavelet of order 1 for the mother wavelet. The discrete wavelet coefficients are defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W(\,j, k) = \int_i X(i) \frac{1}{2^{\,j}} \psi \left(\frac{i-k2^{\,j}}{2^{\,j}} \right). \label{Eq:MF-WL:DCT} \nonumber \end{align} \tag{ 98 }$$

One defines a dyadic interval $\lambda(\,j, k)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda(\,j,k) = \left[k2^{\,j}, (k+1)2^{\,j}\right), \label{Eq:MF-WL:Lambda:jk} \nonumber \end{align} \tag{ 99 }$$

and denotes the union of interval $\lambda$ and its two adjacent neighbors as $3\lambda$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 3\lambda(\,j,k) = \lambda(\,j, k-1) \cup \lambda(\,j, k) \cup \lambda(\,j, k+1). \label{Eq:MF-WL:3Lambda:jk} \nonumber \end{align} \tag{ 100 }$$

The wavelet leader $\newcommand{\e}{{\rm e}} \ell(\,j, k)$ is defined by [295, 296]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ell(\,j, k) = \sup_{\lambda' \subset 3\lambda(\,j, k) } |W(\lambda')|, \label{Eq:MF-WL:Leader} \nonumber \end{align} \tag{ 101 }$$

which corresponds to the largest value of the absolute wavelet coefficients $|W(\,j',k')|$ calculated on intervals $(k-1)2^{\,j} \leqslant 2^{\,j'} k' < (k+2)2^{\,j}$ with $0 < j' \leqslant j$. Note that, to compute the wavelet leaders, all the fine scales $2^{\,j'} \leqslant 2^{\,j}$ must be considered.

The multifractal analysis based on wavelet leaders can be described as follows [295–300]. As usual, we can define the moments of order q by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M(q,j) = \frac{1}{n_j}\sum_{k=1}^{n_j} \ell(\,j,k)^q, \label{Eq:MF-WL:SF} \nonumber \end{align} \tag{ 102 }$$

where n_j is the number of wavelet leaders at scale s  =  2^j. One can also expect the following scaling behavior if the underlying processes are jointly multifractal,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M(q,j) \sim 2^{\,j \zeta(q)}=s^{\zeta(q)} \label{Eq:MF-WL:SF:Scaling} \nonumber \end{align} \tag{ 103 }$$

where $\zeta(q)$ is related to $\tau(q)$ by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \zeta(q) = \tau(q)+1 = qH(q). \nonumber \end{align} \tag{ 104 }$$

We can also provide a canonical approach to determine the singularity strength function $\alpha(q)$ and the multifractal spectrum $f(\alpha(q))$ directly, as for the WTMM case in section 2.3.2.

An interesting application of wavelet leaders is the pointwise wavelet leaders entropy [301], which is defined over different scales (j ) for a given point (k) by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E(k) = -\sum_{j=1}^{m} \rho_{j,k}\log_2\rho_{j,k}, \label{Eq:MF-WL:pL:Entropy} \nonumber \end{align} \tag{ 105 }$$

where $\rho_{j,k}\log_2\rho_{j,k}=0$ if $\rho_{j,k}=0$ and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_{j,k} = \left\{\begin{array}{@{}lll@{}} \frac{W(\,j,k)^2}{\sum_{j=1}^{m}W(\,j,k)^2} & {\mathrm{if}}~W(\,j,k)\neq0 \nonumber \\ 0 & {\mathrm{if}}~W(\,j,k)=0. \end{array} \right. \label{Eq:MF-WL:pL:Entropy:Rho} \nonumber \end{align} \tag{ 106 }$$

When applied to the DJIA index (1928–2011), the resulting wavelet leader entropy is able to identify historical large market fluctuations [301].

Recently, wavelet p -leaders are introduced to characterize the p -exponents of negative pointwise regularity [302]. Based on the wavelet coefficients calculated in equation (98), the wavelet p -leaders are defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ell^{(\,p)}(\,j, k) = \left(\sum_{j'\leqslant{j},\lambda' \subset 3\lambda(\,j, k)} |W(\,j,k)|^{\,p}2^{\,j'-j}\right)^{1/p}. \label{Eq:MF-WL:p-Leader} \nonumber \end{align} \tag{ 107 }$$

The p -leader multifractal formalism relies on the sample moments of the p -leaders [303]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M^{(\,p)}(q,j) = \frac{1}{n_j}\sum_{k=1}^{n_j} \left[\ell^{(\,p)}(\,j,k)\right]^q. \label{Eq:MF-WL:pL:SF} \nonumber \end{align} \tag{ 108 }$$

It is found that MF-DFA characterizes local regularity via the two-exponent, not the Hölder exponent, and the wavelet p -leader multifractal formalism can be viewed as an wavelet-based extension of MF-DFA [303]. The p -leaders suffer from finite-resolution effects, which can be resolved through an explicit and universal closed-form correction [304].

Consider a time series $\{X(i)|i =1,2,\cdots,N\}$, whose data points are usually the cumulative sums of a time series of innovations. Let $\widetilde{X}(t)$ be the local trend function, which is usually determined locally around the data points using some method. The detrended residuals are defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon(i) = X(i)-\widetilde{X}(i). \label{Eq:MF:Detrend:epsilon} \nonumber \end{align} \tag{ 109 }$$

We divide $\newcommand{\e}{{\rm e}} \{\epsilon(i)\}$ into $N_s=\mathrm{int}[N/s]$ non-overlapping segments with boxes of size s, where $\mathrm{int}[y]$ is the largest integer that is no larger than y . We denote the $v$th segment as $\newcommand{\e}{{\rm e}} S_v=\{\epsilon((v-1)s+j)| j=1,2, \cdots, s\}$. The local detrended fluctuation function $F_v(s)$ in the $v$th box is defined as the r.m.s. of the detrended residuals:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[F_v(s)\right]^2 = \frac{1}{s}\sum_{i=1}^{s} \left[\epsilon((v-1)s+j)\right]^2. \label{Eq:MF:Detrend:Fv:s} \nonumber \end{align} \tag{ 110 }$$

The qth-order overall detrended fluctuation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q(s) = \left\{\frac{1}{N_s}\sum_{v=1}^{N_s} {F_v^q(s)}\right\}^{\frac{1}{q}}, \label{Eq:MF:Detrend:Fqs} \nonumber \end{align} \tag{ 111 }$$

where q can take any real value except q  =  0. When q  =  0, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln[F_0(s)] = \frac{1}{N_s}\sum_{v=1}^{N_s}{\ln[F_v(s)]}, \label{Eq:MF:Detrend:Fq0} \nonumber \end{align} \tag{ 112 }$$

according to L'Hôspital's rule.

When the whole series $\newcommand{\e}{{\rm e}} \{\epsilon(i)\}$ cannot be completely covered by N_s boxes, one can use 2N_s boxes to cover the series from both ends of the series, where $N_s=\mathrm{int}[N/s]$ and a short part at each end of the residuals may remain uncovered [178]. We denote the $v$th segment partitioned from left to right as $\newcommand{\e}{{\rm e}} S_v=\{\epsilon((v-1)s+j)| j=1,2, \cdots, s\}$ for $v=1,2, \cdots, N_s$ and the $v$th segment partitioned from left to right as $\newcommand{\e}{{\rm e}} S_v=\{\epsilon(N-(v-N_s)s+j)| j=1,2, \cdots, s\}$ for $v=N_s+1,N_s+2, \cdots, 2N_s$. The local detrended fluctuation function $F_v(s)$ in the $v$th box is defined as the r.m.s. of the detrended residuals,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[F_v(s)\right]^2 = \frac{1}{s}\sum_{i=1}^{s} \left[\epsilon((v-1)s+j)\right]^2 \label{Eq:MF:Detrend:Fv:s:L2R} \nonumber \end{align} \tag{ 113a }$$

for $v=1,2, \cdots, N_s$ and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[F_v(s)\right]^2 = \frac{1}{s}\sum_{i=1}^{s} \left[\epsilon(N-(v-N_s)s+j)\right]^2 \label{Eq:MF:Detrend:Fv:s:R2L} \nonumber \end{align} \tag{ 113b }$$

for $v=N_s+1,N_s+2, \cdots, 2N_s$. The qth-order overall detrended fluctuation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q(s) = \left\{\frac{1}{2N_s}\sum_{v=1}^{2N_s} \left[F_v(s)\right]^q\right\}^{\frac{1}{q}}, \label{Eq:MF:Detrend:Fqs:LR} \nonumber \end{align} \tag{ 114 }$$

where q can take any real value except for q  =  0. When q  =  0, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln[F_0(s)] = \frac{1}{2N_s}\sum_{v=1}^{2N_s}{\ln[F_v(s)]}, \label{Eq:MF:Detrend:Fq0:LR} \nonumber \end{align} \tag{ 115 }$$

according to L'Hôspital's rule.

Varying the values of the segment size s, we can determine the power-law relation between the function F_q(s) and the size scale s,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q(s) \sim {s}^{H(q)}. \label{Eq:MF:Detrend:hq} \nonumber \end{align} \tag{ 116 }$$

According to the standard multifractal formalism, the multifractal scaling exponent $\tau(q)$ can be used to characterize the multifractal properties, which reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q)=qH(q)-D_f, \label{Eq:tau:hq} \nonumber \end{align} \tag{ 117 }$$

where D_f is the fractal dimension of the geometric support of the multifractal measure [178]. For time series analysis, we have D_f  =  1. If the scaling exponent function $\tau(q)$ is a nonlinear function of q, the time series is regarded to have a multifractal nature.

Following the approach determining directly $f(\alpha)$ from the partition function [197, 198, 305], the singularity strength and its spectrum in the detrended fluctuation approaches can also be determined directly [306]. Defining the canonical measure $\mu(q,s,v)$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu(q,s,v) = \frac{F_v^q(s)}{\sum_{v=1}^{N_s} F_v^q(s)}, \label{Eq:MF-DMA-CN:mu} \nonumber \end{align} \tag{ 118 }$$

we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(q) = \lim_{s\to0} \frac{\sum_{v=1}^{N_s} \mu(q,s,v) \ln F_v(s)}{\ln{s}}, \label{Eq:MF-DMA-CN:alpha} \nonumber \end{align} \tag{ 119a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha(q)) = \lim_{s\to0}\frac{\sum_{v=1}^{N_s} \mu(q,s,v) \ln\left[\mu(q,s,v)\right]}{\ln{s}}. \label{Eq:MF-DMA-CN:fq} \nonumber \end{align} \tag{ 119b }$$

Numerical experiments based on the p -model and fractional Brownian motions showed that MF-DMA and its direct determination variant have comparative performance to unveil the fractal and multifractal properties [306].

The detrended fluctuation analysis (DFA) was invented originally to investigate the long-range dependence in coding and noncoding DNA nucleotide sequences [307]. The multifractal extension of DFA, MF-DFA, was developed shortly thereafter [308, 309]. Since the work of Kantelhardt et al [178], MF-DFA has become one of the most important methods for multifractal analysis in diverse fields.

In the MF-DFA method, the trend function is a polynomial [178], and the linear trend is most generally adopted. For each point i, its trend value $\widetilde{X}(i)$ is determined together with other points in the same box $S_v$, using a polynomial fitting of order l to the s data points in $i\in S_v$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \widetilde{X}(i) = \sum_{k=0}^\ell a_k i^k, ~~l=1,2, \cdots. \label{Eq:MF:Trend:Polynomial} \nonumber \end{align} \tag{ 120 }$$

The detrending process for $\newcommand{\e}{{\rm e}} \ell=1$ is shown in figure 6. The figure shows a case with N_s  <  N/s in which the data points at each end of the time series should not be included to form a box. Otherwise, one obtains incorrect results with $f_{\max}=f(q=0)>1$, as will be elaborated in section 5.1.2. In the standard procedure implementing MF-DFA, researchers remove the constant shift $\langle{\Delta{X}}\rangle$ from the increments $\Delta{X}$ and work on the cumulative profile

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x(i) = \sum_{j=1}^{i} \left[\Delta{X}(j)-\langle{\Delta{X}}\rangle\right],~i = 1, 2, \cdots, N, \label{Eq:cumsum} \nonumber \end{align} \tag{ 121 }$$

where $\langle{\Delta{X}}\rangle$ is the sample mean of the $\Delta{X}(i)$ series. However, when the trend function is polynomial, this manipulation does not have any impact on the results.

Zoom In Zoom Out Reset image size

A subtle issue of the MF-DFA method is the choice of polynomial order. Using classic mathematical models (fractional Brownian motion, Lévy process and binomial measure) and real-world time series (foreign exchange rate and literary texts), Oświęcimka et al found that the calculated singularity spectra could be very sensitive to the order of the detrending polynomial and this sensitivity depends on the analyzed time series [310]. Without a deep understanding of the underlying mechanism driving the dynamics, it is hard to determine which polynomial order should be used. Anyway, for financial time series such as returns, we can check if the value of the Hurst index $H(2)$ is close to 0.5, based on the fundamental stylized fact of the absence of long memory in financial returns [94], which however might not hold at the transaction level [130].

A possible solution is to use 'feasible trend functions' in the multifractal feasibly detrended fluctuation analysis (MFFDFA) [311]. One presets a priori a set of any trend functions, say, $\{ax^2+bx+c, a\sin(x^2)+bx+c, ax^3+bx+c\}$, and fit each segment with all the preset trend functions. For each segment, only one detrended function is chosen which has the maximum coefficient of determination R². Hence, for each time scale s, the trend functions for different segments are not necessary the same. Numerical experiments on fractional Brownian motions and binomial measures show that the MFFDFA method outperforms the MF-DFA method [311]. For binomial measures, the improvement is partially due to the fact that the fluctuation functions F_q(s) for MFFDFA have much smaller log-periodic oscillation amplitudes, which reduces the impact of the choice of scaling ranges [312].

Based on the moving average (MA) technique for estimating the Hurst exponent of self-affinity signals [313], the detrending moving average (DMA) analysis considers the difference between the original signal and its moving average function [314, 315]. The DMA method can also be easily implemented to estimate the correlation properties of non-stationary series without any a priori assumptions, which is widely applied to the analysis of real-world time series [316–323] and synthetic signals as well [324]. Extensive numerical experiments unveil that the performance of DMA is comparable to DFA with slightly different priorities under different situations and both DFA and DMA remain 'The Methods of Choice' in determining the Hurst index of time series [325–327].

The DMA method was generalized to MF-DMA for the investigation of the multifractal nature of time series [328]. The main difference between the MF-DFA and MF-DMA algorithms is about the detrending procedure. In the MF-DMA approach, one calculates the moving average function $\widetilde{X}(i)$ in a moving window of size s [321],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \widetilde{X}(i)=\frac{1}{s}\sum_{k=-\lfloor(s-1)\theta\rfloor}^{\lceil(s-1)(1-\theta)\rceil}X(i-k), \label{Eq:MF-DMA:trend} \nonumber \end{align} \tag{ 122 }$$

where $\lfloor{y}\rfloor$ is the largest integer not greater than y , $\lceil{y}\rceil$ is the smallest integer not smaller than y , and $\theta\in[0,1]$ is the position parameter. In other words, one considers $\lceil(s-1)(1-\theta)\rceil$ data points in the past and $\lfloor(s-1)\theta\rfloor$ points in the future to determine $\widetilde{X}(i)$. When i is close to each of the endpoints, the size of the box might be less than s. One can either abandon those points or calculate the moving averages directly. Usually, we consider three special cases, corresponding to three representative versions of MF-DMA:

• MF-BDMA. The first case $\theta=0$ refers to the backward moving average [325], in which the moving average function $\widetilde{X}(i)$ is calculated over all the past s  −  1 data points of the signal.
• MF-CDMA. The second case $\theta=0.5$ corresponds to the centered moving average [325], where $\widetilde{X}(i)$ contains half past and half future information in each window. This version is also termed MF-CMA [329].
• MF-FDMA. The third case $\theta=1$ is called the forward moving average, where $\widetilde{X}(i)$ considers the trend information of s  −  1 data points in the future.

We note that the size of the window used to calculate the moving averages must be identical to the segment size in calculating the detrended functions [328, 329]. The detrending process is illustrated in figure 7.

Zoom In Zoom Out Reset image size

Although determining the intrinsic trend embedded in non-stationary time series is difficult, there are many different options [330]. Because the intrinsic trend of a time series is usually unknown, methods have different performances for different types of time series. Hence, quite a few variants of MF-DFA and MF-DMA have been proposed. Local trends can be further classified into discontinuous ones and continuous ones. MF-DFA and MF-DMA belong respectively to these two classes (see also figures 6 and 7).

Zhou and Leung proposed the multifractal temporally weighted detrended fluctuation analysis (MF-TWDFA) in moving windows [331]. For each i, they fit a linear function to the data points $(X(i-s), \cdots, X(i), \cdots, X(i+s))$ in the moving window of size 2s  +  1, using weighted least-squares regression. The objective function is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal{O}}(i; s) = \sum_{j=-s}^s\left[X(i+j)-\sum_{k=0}^la_k(i+j)^k\right]^2w_{ij}, \nonumber \end{align} \tag{ 123 }$$

where $|i-j|\leqslant{s}$, w_ij's are weights that decrease with increasing $|i-j|$. In other words, the deviation $\left[X-a_0-a_1\times(i+j)\right]^2$ has a larger weight when j  is far from i. The weights w_ij in MF-TWDFA are defined by [331]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w_{ij} = \left[1-\left(\frac{i-j}{s}\right)^2 \right]^2. \nonumber \end{align} \tag{ 124 }$$

For the equally weighted case,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w_{ij} = 1, \nonumber \end{align} \tag{ 125 }$$

we obtain the multifractal moving window detrended fluctuation analysis (MF-MWDFA) [331]. These two methods belong actually to the continuous class, since the trend changes mildly without jumps. Indeed, when l  =  0, we recover the centered MF-DMA method [331].

The empirical mode decomposition (EMD) is an innovative method of analysis for nonlinear and non-stationary time series, which decomposes the time series $\{X((v-1)s+i)\}_{i=1}^s$ into a finite number of intrinsic mode functions (IMFs) [332]. Let $\newcommand{\tr}{{\rm tr}} \{x(i)\}_{i=1}^s\triangleq\{X((v-1)s+i)\}_{i=1}^s$. The decomposition is a sifting process, which has six steps [330]: (1) Identify all extrema of $x(i)$; (2) Interpolate the local maxima to form an upper envelope $U(i)$; (3) Interpolate the local minima to form a lower envelope $L(i)$; (4) Calculate the mean envelope:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu(i)=[U(i)+L(i)]/2; \label{Eq:EMD:mu} \nonumber \end{align} \tag{ 126 }$$

(5) Extract the mean from the signal

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_k(i) = X(i)-\mu(i); \label{Eq:EMD:Ck} \nonumber \end{align} \tag{ 127 }$$

(6) Check whether C_k(i) satisfies the IMF conditions. If the IMF conditions are fulfilled, C_k(i) is an IMF, and the sifting stops. Otherwise, let x(i)  =  C_k(i) and keep sifting. Finally, we obtain the EMD-based local trend of the time series segment $\{X((v-1)s+i)\}_{i=1}^s$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle r_n\triangleq\widetilde{X}((v-1)s+i)= X((v-1)s+i)-\sum_{k=1}^K C_k(i). \label{Eq:EMD:rn} \nonumber \end{align} \tag{ 128 }$$

Qian et al used this local trend in the MF-DFA method and developed the EMD-based MF-DFA method [333]. Numerical experiments on multifractal binomial measures indicate that the EMD-based MF-DFA method has comparable performance with respect to the MF-DFA method except that the EMD-based method performs slightly better for small singularities [333]. The method has been applied to analyze stock market indices [333, 334], foreign exchange rates [335], and certified emission reduction (CER) and European Union Allowances (EUA) futures prices [336], all of which confirm the presence of multifractality in financial returns.

In order to handle time series with highly non-stationary sinusoidal and polynomial overlying trends, Horvatic et al proposed to detrend with polynomials of varying orders [337]. In this method, the polynomial order $\newcommand{\e}{{\rm e}} \ell$ is not fixed but increases with the scale s and the detrended fluctuations at different scales should be normalized to secure an identical intercept. Although an objective criterion for the determination of the varying order is lacking, a subjective choice ($\newcommand{\e}{{\rm e}} \ell=1$ for s  =  4, $\newcommand{\e}{{\rm e}} \ell=2$ for s  =  6, $\newcommand{\e}{{\rm e}} \ell=4$ for s  =  10, $\newcommand{\e}{{\rm e}} \ell=8$ for s  =  18, and $\newcommand{\e}{{\rm e}} \ell=16$ for s  >  18) seemingly works well [337]. The underlying logic is that we need high-order polynomials to detrend properly the oscillations.

Du et al constructed the dual-tree complex wavelet transform based multifractal detrended fluctuation analysis (DTCWT-MFDFA) for non-stationary time series [338]. Based on the anti-aliasing and nearly shift-invariance of the dual-tree complex wavelet transform [339, 340], the raw time series is decomposed through the pyramid algorithm to extract the scale-dependent trends and the fluctuations from the wavelet coefficients. The length of the non-overlapping segment on each time scale is computed through the Hilbert transform and each of the extracted fluctuations is partitioned into non-overlapping segments with identical size. The detrended fluctuation function for each segment can be calculated. The DTCWT-MFDFA method is found to outperform the MF-DFA method for multifractal binomial measures [338].

In order to perform the MF-DFA, Manimaran et al proposed to detrend using discrete wavelet transform [341]. The discrete wavelets belonging to the Daubechies family are adopted to analyze the raw time series up to level s. The resulting low-pass coefficients capture the polynomial trends in the data. One then reconstructs a time series using these low-pass coefficients to represent the local trends at different scales. The method has been tested for synthetical time series and applied to NASDAQ and BSE indices [342, 343]. Ghosh et al used Db4 and Db6 basis sets to isolate respectively the local linear and quadratic trends at different scales from the closing prices of fifteen large-capitalization companies listed on the New York Stock Exchange (NYSE) [344]. Multifractality is observed in the investigated stock indices and individual stocks.

We note that other detrending methods used to preprocess the raw time series in section 5.3.6 can also be used in the local detrending at different scales to construct variants of MF-DFA.

The concept of the multiplier was introduced by Novikov to describe the intermittency and scale self-similarity in turbulent flows [164]. In econophysics, one can apply this concept to study volatility multipliers. For high-frequency financial data, we first construct an additive measure in the time interval $[t_1,t_2]$, which is the sum of absolute returns:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu([t_1,t_2]) = \sum_{t = t_1}^{t_2}|r(t)|. \label{Eq:MF:VolatilityMeasure} \nonumber \end{align} \tag{ 129 }$$

The absolute return is actually a measure of the volatility on $[t_1,t_2]$ [85]. The volatility time series $|r(t)|$ is divided into boxes of identical size s. Each mother box is further partitioned into a daughter boxes, in which a is called a base. The multiplier m_a,s is determined by the ratio of the measure on a daughter box to that on her mother box [347, 348].

If the multiplier is scale invariant, for any two bases a and b, the multiplier density functions $p_a(m_a)$ and $p_b(m_b)$ are related through Mellin transform $\mathbb{M}$ as follows [348]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle [\mathbb{M} \{\,p_a(m_a) \}]^{1 / \ln a} = [\mathbb{M} \{\,p_b(m_b) \}]^{1 / \ln{b}}, \label{Eq:Mellin} \nonumber \end{align} \tag{ 130 }$$

which leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\ln \int_0^1 m_a^q p_a(m_a) d m_a}{\ln a} = \frac{\ln \int_0^1 m_b^q p_b(m_b) {{{\rm d}}} m_b}{\ln{b}}. \label{Eq:MF:multiplier:base} \nonumber \end{align} \tag{ 131 }$$

The mass scaling exponent $\tau(q)$ of the moment of m_a can be obtained as [347, 348]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = - D_0 - \frac{\ln\langle m_a^q \rangle}{\ln a}, \label{Eq:MF:multiplier:tau} \nonumber \end{align} \tag{ 132 }$$

where D₀ is the fractal dimension of the support of the measure. For time series, we have D₀  =  1.

The local singularity exponent $\alpha$ and its spectrum $f(\alpha)$ are related to $\tau(q)$ through the Legendre transform. It follows that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(q) = -\frac{\langle m_a^q \ln m_a \rangle}{\langle m_a^q \rangle \ln a} \label{Eq:MF:multiplier:alpha} \nonumber \end{align} \tag{ 133 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) = \frac{\langle m_a^q \rangle \ln \langle m_a^q \rangle - \langle m_a^q \ln m_a^q \rangle}{\langle m_a^q \rangle \ln a}. \label{Eq:MF:multiplier:falpha} \nonumber \end{align} \tag{ 134 }$$

Equations (131)–(134) imply that, for each q, each of these characteristic multifractal quantities converges to a constant in the scaling range and is independent of the base a.

Note that, for financial volatility, q  >  −1 [349], which is quite analogous to the situation in turbulence [267]. Because the returns can be zero, the multiplier m can be zero and the probability density $p_a(0)\neq0$ at m_a  =  0. When m is smaller than a small number $\delta$, $p(m)=p(0)$ is a constant. For a given small $\delta$, $\int_{\delta}^1m^qp(m)dm = C$ is finite. However, the moment

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_0^1m^qp(m) {{{\rm d}}}m & =\int_0^{\delta}m^qp(m) {{{\rm d}}}m+ C \nonumber \\ &=\left.p(0)\frac{1}{q+1}m^{q+1}\right|_0^{\delta}+C \nonumber \end{align} \tag{ 135 }$$

diverges when $q\leqslant-1$. This property is verified empirically [349], based on the idea that the integrand $m_a^qp_a(m_a)$ diverges for a given q [268, 350].

The scale-invariant multiplier distribution is perceived to be more fundamental than the standard multifractal spectrum $f(\alpha)$ [347, 348]. Moreover, it needs exponentially less computational time and is more accurate than the partition function approach [347, 348].

Using 1 min high-frequency data of the S&P 500 index, Jiang and Zhou reported that the volatility multiplier is scale invariant and the volatility has the property of multifractality [349]. Niu and Wang confirmed the presence of scale invariance in the multiplier distributions in the Shanghai SSEC index and in the simulated data from the model of 'voter interacting dynamic system' [351].

The Hilbert–Huang transform combines the EMD and the Hilbert spectral analysis [332, 352], in which the Hilbert transform is applied to each component C_k(i) in equation (127). The overall trend can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \widetilde{X}(i) &= X(i)-\sum_{k=1}^nC_k(i) \nonumber \\ &= X(i)-\sum_{k=1}^n\mathcal{A}_k(i)\exp\left({\mathbf{j}}\int\omega_k(\ell){{{\rm d}}}\ell\right), \label{Eq:HHT:expression} \nonumber \end{align} \tag{ 136 }$$

where ${\mathbf{j}}^2=-1$. For each mode, the Hilbert–Huang spectrum is defined as the squared amplitude

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H(\omega,i)=\mathcal{A}^2(\omega,i). \label{Eq:HHT:Spectrum} \nonumber \end{align} \tag{ 137 }$$

The Hilbert–Huang marginal spectrum of the original time series is then written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h(\omega)=\int H(\omega,i) \, {\rm d}i, \label{Eq:HHT:MarginalSpectrum} \nonumber \end{align} \tag{ 138 }$$

which corresponds to an energy density at frequency $\omega$. This method is known as EMD-based Hilbert spectral analysis (EMD-HSA), which we term the Hilbert–Huang spectral analysis (HHSA) for simplicity.

Huang, Schmitt, Lu and Liu proposed the multifractal Hilbert–Huang spectral analysis (MF-HHSA) by generalizing the Hilbert–Huang transform from the second order to arbitrary order [353, 354]. The Hilbert–Huang spectrum $H(\omega,i)$ represents the original signal at the local level and can be used to define the joint probability density function $p(\omega,\mathcal{A})$ of the frequency $\omega_k$ and amplitude $\mathcal{A}_k$. The Hilbert–Huang marginal spectrum in equation (138) can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h(\omega)=\int_0^{\infty} p(\omega,\mathcal{A}) \mathcal{A}^2 {{{\rm d}}} \mathcal{A}, \label{Eq:HHT:MarginalSpectrum:2} \nonumber \end{align} \tag{ 139 }$$

which corresponds to the second-order moment. Huang et al generalize equation (139) into arbitrary-order moments

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{L}_q(\omega)=\int_0^{\infty} p(\omega,\mathcal{A}) \mathcal{A}^q {{{\rm d}}} \mathcal{A}, \label{Eq:HHT:Moments} \nonumber \end{align} \tag{ 140 }$$

where $q \geqslant 0$ and $h(\omega)=\mathcal{L}_2(\omega)$. In the inertial range, we assume the following scaling relation,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{L}_q(\omega) \sim \omega^{-\tau(q)}, \label{Eq:HHT:Moments:tau:q} \nonumber \end{align} \tag{ 141 }$$

where $\tau(q)$ is the corresponding scaling exponent function in the amplitude-frequency space.

Li and Huang applied this MF-HHSA method to the CSI 300 index of the Chinese stock market and confirmed the presence of multifractality [355].

Welter and Esquef proposed an EMD-based dominant amplitude multifractal formalism (EMD-DAMF), which is implemented as follows [356]:

• Step 1:  
  Perform EMD on $X(i)$ to obtain K IMFs C_k(i).
• Step 2:  
  Determine the amplitude a_k(i) and characteristic timescale s for C_k(i). The characteristic timescale s could be estimated in each mode by the Hilbert spectral analysis (HSA) as the reciprocal of the mean frequency, or simply by twice the mean time interval among all adjacent zero crossings.
• Step 3:  
  Determine the dominant amplitude coefficients

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_{k,j}:=\sup_{k'\leqslant{k}}\{\max[|a_{k'}(i\in{I}_{k,j})|]\}, \label{Eq:MF:EMD:DA:ukj} \nonumber \end{align} \tag{ 142 }$$

  where $k=1,2,\cdots,K$, $j=1,2,\cdots, n_k$, n_k is the number of local maxima of a_k(i), and I_k,j is a time support around the j th maxima of a_k(i).
• Step 4:  
  Calculate the qth order moment for different q values:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_q(s_k)=\frac{1}{n_k}\sum_{i=1}^{n_k}(u_{k,j})^q,~k=1,2,\cdots, K, \label{Eq:MF:EMD:DA:Sk} \nonumber \end{align} \tag{ 143 }$$

  where the characteristic time scale s_k can be obtained by HSA in each mode as the reciprocal of the mean (linear) frequency or simply by the mean time interval among all adjacent zero crossings.
• Step 5:  
  Estimate the scaling function

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_q(s_k) \sim s_k^{\tau(q)+1-q}. \label{Eq:MF:EMD:DA:PL} \nonumber \end{align} \tag{ 144 }$$

In real applications, Welter and Esquef suggested that the number of siftings is fixed to ten (that is, K  =  10) and signal envelopes are computed through cubic spline interpolation [356]. They validated the method with fractional Brownian motions [357], Lévy processes [358], Multifractal random wavelet cascades [359, 360], and multifractal multiplicative cascades from the p  model [361]. Note that the EMD-DAMF method has a small number n_k of amplitude maxima at larger scale modes and a limited number of scales [356]. The latter drawback might impact the determination of the scaling range.

Alternatively, Zhou et al proposed the so-called EMD-DFA method [362]. The integrated time series $X(i)$ is decomposed into a set of intrinsic mode functions and a residual or trend component. All intrinsic time scales s are identified by computing the number of data points between each neighboring local minima throughout each intrinsic mode function (IMF). For a given s, all IMFs are resampled through replacing by 0's all data points corresponding to scales not equal to s and the fluctuation time series X_s(i) is obtained by summing all resampled IMFs. The overall fluctuation at scale s is the root mean square of X_s(i):

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle F_2(s) = \left[\frac{1}{N}\sum_{i=1}^N X_s^2(i) \right]^{1/2}, \nonumber \end{align*}$$

which produces better scaling than the DFA method [362]. The EMD-DFA method can be extended for multifractal time series by working on the q-order moments:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q(s) = \left[\frac{1}{N}\sum_{i=1}^N |X_s(i)|^q \right]^{1/q}. \nonumber \end{align} \tag{ 145 }$$

Note that the EMD-DFA method does not have an explicit detrending form like equation (109).

Gieraltowski et al proposed the multiscale multifractal analysis (MMA) based on the MF-DFA method [363]. After calculating all values of the overall fluctuations F_q(s) by the MF-DFA method, rather than estimating $H(q)$ from equation (116) in the scaling range, they proceeded to estimate the local generalized Hurst exponents $H(q,s)$ in moving windows around different scales. In this method, two free parameters are introduced, the moving window size and the step. Although the original MMA method is designed for MF-DFA, it can obviously be applied to other methods of multifractal analysis.

The MMA method has been applied to investigate the multifractal nature of daily returns of the Dow Jones Industrial Average (DJIA), New York Stock Exchange Index (NYSE), S&P 500 Index, Heng Seng Index (HSI), Shanghai Stock Exchange Composite Index and Shenzhen Stock Exchange Component Index (SZCI) from 12 May 1992 to 8 May 2012 [364, 365] and eight stock market indices (CAC 40, DAX, FTSE 100, HSI, NIKKEI 225, SSEC, NASDAQ, S&P 500) from 3 January 2000 to 1 October 2014 [366].

More generally, this method is equivalent to calculate the local logarithmic slopes of the F_q(s) functions:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h(q,s) = \lim_{s'\to{s}}\frac{{\rm d}\ln{F_q(s')}}{{\rm d}\ln{s'}}, \label{Eq:MF:Detrend:hq:Local} \nonumber \end{align} \tag{ 146 }$$

which is widely utilized in the study of structure functions in turbulence [265, 268], as well as in DFA [367]. The idea to estimate local logarithmic slopes is especially informative when the scaling range is narrow and it can serve as a method for the determination of scaling ranges. Different numerical methods can be adopted to estimate numerically the local logarithmic slopes $H(q,s)$ and the linear fitting in a moving window suggested in the MMA method is one of them. Empirical analyses have been performed on financial time series, such as the WTI crude oil market [368] and stock market indices [369], without referring to the name of multiscale multifractal analysis.

To quantify the evolution of a multifractal spectrum, one can perform a multifractal analysis of time series in moving windows or rolling windows. We would like to stress that short window sizes will reduce the estimation accuracy of multifractal quantities and increase the uncertainty of the results. Time-varying multifractal quantities have practical implications, as discussed in section 8.2.4.

Alternatively, Xiong and Zhang proposed the time-varying multifractal spectrum distribution (TM-MFSD) [370]. The key idea is to perform multifractal analysis on the instantaneous autocorrelation function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r(i, \delta) = {\mathbb{E}}[X^*(i-\Delta/2)X(i+\Delta/2)], \nonumber \end{align} \tag{ 147 }$$

rather than the raw time series $X(i)$, where the time-delayed conjugation X^*(i) of the analyzed series $X(i)$ is selected as the window function. To obtain the time-varying multifractal spectrum distribution of a time series, one can adopt different multifractal analysis methods described in previous sections, such as WTMM [370], wavelet leaders [371, 372], MF-DMA [373], and MF-DFA [374]. Essentially, the multifractal analysis is performed with these methods on moving windows.

A classic method in nonlinear dynamics is to embed properly the time series into a phase space and estimate the generalized dimensions D_q of the data points in the phase space [167, 188, 189, 375], which is however less applied in econophysics. Based on the delay embedding technique [376], the time series $X(i)$ can be converted into a sequence of vectors:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overrightarrow{X}(\,j) = [X(\,j),X(\,j+\tau), \cdots, X(\,j+(d-1)\tau)], \label{Eq:Enbedding:Vector} \nonumber \end{align} \tag{ 148 }$$

where $j=1,\cdots,N_v$, $N_v$ is the number of points (or vectors) in the phase space, $\tau$ is a suitably chosen time delay, and d is the embedding dimension of the phase space. The embedding dimension d and the delay $\tau$ can be determined in different ways [189, 377–381].

The relative number of data points within a distance R to point i is calculated by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_{i}(s) = \frac{1}{N_v} \sum_{j=1,j\neq i}^{N_v} {\mathcal{H}}\left(s-\left|\overrightarrow{X}(i)-\overrightarrow{X}(\,j)\right|\right), \label{Eq:Enbedding:ci:R} \nonumber \end{align} \tag{ 149 }$$

where ${\mathcal{H}}$ is the Heaviside function. The generalized correlation integral C_q(s) is [382]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{e.5} C_{q} (s) = \left\{\frac{1}{N_{v}} \sum_{i}^{N_{v}} \left[c_i (s)\right]^{q-1}\right\}^{1/(q-1)}, \nonumber \end{align} \tag{ 150 }$$

which is the generalization of the correlation integral [167, 188, 189]. The generalized dimensions are given by [375]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{e.6} D_{q} \equiv \lim_{s \rightarrow 0} \frac{\log C_q (s)}{\log s}. \nonumber \end{align} \tag{ 151 }$$

Lee investigated the multifractal characteristics of the 1 min volatility of the Korea composite stock price index KOSPI from 20 March 1992 to 28 February 2007 [383]. The volatility time series is first detrended with a wavelet filter and then a multifractal analysis in the phase space with $\tau=1$ and d  =  1 is performed to confirm the presence of multifractality.

Inspired by the stylized fact that the correlation between two assets is usually asymmetric [384–387], Alvarez-Ramirez et al developed the asymmetric multifractal detrended fluctuation analysis (A-MF-DFA) to detect asymmetric multifractal scaling in individual time series [388]. In this method, the trend function is linear such that a_i  =  0 for all i  >  1 in equation (120). The main idea is to calculate the upwards (downwards) fluctuation functions in which the removed linear trends have non-negative (negative) slopes. Mathematically, the qth-order upwards overall detrended fluctuation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q^+(s) &= \left\{\frac{1}{\sum_{v=1}^{2N_s} \frac{1+k_v}{2}}\sum_{v=1}^{2N_s} \left[\frac{1+k_v}{2}F_v(s)\right]^q\right\}^{\frac{1}{q}} \nonumber \\ & = \left\langle[F_v(s)]^q\right\rangle|_{a_1\geqslant0}, \label{Eq:MF:A-MF-DFA:Fqs:LR:p} \nonumber \end{align} \tag{ 152 }$$

and the qth-order downwards overall detrended fluctuation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q^-(s) = \left\langle[F_v(s)]^q\right\rangle|_{a_1<0}, \label{Eq:MF:A-MF-DFA:Fqs:LR:n} \nonumber \end{align} \tag{ 153 }$$

where $k_v=1$ if the slope of the linear trend (120) is nonnegative, that is $a_1\geqslant0$, and $k_v=-1$ if the slope a₁  <  0. When q  =  0, we apply L'Hôspital's rule as in equation (115). The respective generalized Hurst exponents can be determined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q^{\pm}(s) \sim {s}^{H^{\pm}(q)}. \label{Eq:MF:A-MF-DFA:hq} \nonumber \end{align} \tag{ 154 }$$

Note that the original A-MF-DFA method uses the mean absolute value of the residuals in each segment $v$ [388].

Empirical studies using the A-MF-DFA method investigated and confirmed the asymmetric multifractal behavior in the daily returns of the DJIA index from 23 May 1980 to 25 August 2008 [389], the WTI oil price from 2 January 1986 to 14 December 2010 [389], the SSEC index from 19 December 1990 to 27 April 2012 and the SZCI index from 2 April 1991 to 27 April 2012 [390], and the DJIA, NASDAQ, NYSE, and S&P500 indices from 1 January 1991 to 31 December 2015 [391].

The A-MF-DFA method can be easily extended to design the asymmetric multifractal detrending moving average analysis (A-MF-DMA), in which moving averages serve as the trend function and the direction (or 'sign') is determined exactly as for the A-MF-DFA method [392]. The asymmetric multiscale detrended fluctuation analysis is obtained by calculating the local exponents at different time scales, which has been applied to investigate the hourly returns of the California electricity spot prices during years 1999 and 2000 [393].

Cao et al developed the asymmetric multifractal detrended cross-correlation analysis (MF-ADCCA) method [394], as an asymmetry extension of the MF-X-DFA method (see section 3.4). The method has been applied to identify asymmetric cross-correlations between the daily returns of the SSEC index and six foreign exchange rates from 22 July 2005 to 13 January 2012 [394], the 5 min returns of the CSI 300 index spot and futures returns from 30 June 2011 to 7 June 2013 [395], the daily returns of the gold spot prices in the Shanghai Gold Exchange (Au 99.95) and the New York Gold Exchange from 2 January 2003 to 27 April 2012 [396], the EU ETS for carbon emission trading products versus energy futures—Brent crude oil (Brent), Richards bay coal (Coal), UK base electricity (Electricity) and UK natural gas (Gas)—from 1 January 2008 to 31 December 2012 [397], the daily returns of the WTI price and eight foreign exchange rates (CAD, MXN, NOK, GBP, JPY, AUD, EUR and KRW) against the USD from 4 January 2000 to 31 December 2014 [398], the Shanghai SSEC and four other indices (US S&P 500, Germany DAX, India BSESN, and Brazil BVSP) from 1 January 2002 to 26 September 2014 [399], the WTI crude oil prices and eight stock merket indices (Brazil IBOVESPA, Chile IPSA, Argentina MERVAL, Mexico IPC, Colombia COLCAP, Peru SPBVL, S&P500 and DJIA) from the 28 January 2005 to 29 December 2016 [400], and the daily data of Bitcoin till 25 August 2017 [401].

Chen and Zheng proposed the asymmetric joint multifractal analysis based on partition functions [402], as an extension of the MF-X-PF method (see section 3.1). This extension becomes feasible for financial volatility because one can use the associated returns to define the signs in the segments.

The diffusion entropy analysis (DEA) for time series $\{X(i)\}_{i=1}^N$ is based on the Shannon entropy [403–406], which can be determined by taking asset prices as an example. For a given diffusion time or time scale s, we calculate the return series

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Delta{X}(i,s) = X(i)-X(i-s), ~i=1,2,\cdots, N-s, \nonumber \end{align*}$$

which denotes the total diffusion distance of s jump steps. For each s, the values of $\{\Delta{X}(i,s)\}_{i=1}^N$ are covered by N_s non-overlapping intervals of length $\newcommand{\e}{{\rm e}} \epsilon(s)$. The interval length $\newcommand{\e}{{\rm e}} \epsilon(s)$ is chosen as a fraction of the square root of the variance of the fluctuations $\{\Delta{X}(i,1)\}$, which is independent of the time scale s. We count the number N_j(s) of the returns $\Delta{X}(i,s)$ that fall in the j th interval and the probability (or frequency) by means of

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle p_{j}(s) \equiv \frac{N_j(s)}{N-s}. \nonumber \end{align*}$$

The Shannon entropy of the diffusion process is determined as

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle I(s) = - \sum_{j=1}^{N_s} p_j(s) \ln p_j(s). \nonumber \end{align*}$$

The scaling behavior is characterized by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle I(s) \sim \delta \ln{s}, \nonumber \end{align*}$$

where $\delta$ is the scaling exponent.

For fractional Brownian motions, the DEA scaling exponent is equal to the Hurst index, that is, $\delta=H$ [405], which does not necessary hold for other processes such as Lévy processes [405, 406]. For financial assets, it seems that one should work on returns to obtain approximately $\delta=H$; otherwise, $\delta$ will be greater than 0.9 [407–409]. There is also an evident finite-size effect that the entropy $I(s)$ deviates from the 'theoretical' value for small scales and converges to the theoretical value for large s [406]. It is possible to overcome this drawback using the idea of integrated moments [410] (see section 5.1.5 for details).

Using DEA, Palatella et al have provided very interesting discussions of multifractality, non-linear time correlation and fat tails for financial market activities [69, 70]. An alternative generalization of DEA for multifractal analysis is to replace the Shannon entropy by the Rényi entropy expressed in equation (20), leading to the multifractal diffusion entropy analysis (MF-DEA), whose main output is a spectrum of scaling exponents $\delta(q)$ [411]. Analogous to the multiscale multifractal analysis based on the MF-DFA method [363] (see section 2.5.4), the multiscale multifractal diffusion entropy analysis can be designed [412]. We notice that $\delta(q)$ is not related directly to the well-known generalized dimensions D_q. Hence, the MF-DEA method unveils multifractality in certain mappings of the original time series, but not in the original time series per se.

We should however mention that the MF-DEA method has some drawbacks [413], as all approaches have weaknesses accompanying their positive properties. First, the MF-DEA method for negative values of q needs to be re-defined. Second, in the DEA or MF-DEA method, the estimation of the probabilities p _j(s) has a crucial impact on the estimation of the Rényi entropy and consequently the $\delta(q)$ spectrum [414]. As a result, an optimal choice of $\newcommand{\e}{{\rm e}} \epsilon(s)$ is required [414]. Third, in general, the MF-DEA method does not allow one to investigate correlations between extremely small and extremely large fluctuations. Fourth, the Rényi entropy has a particular property of conserving the scaling term.

Xu and Beck applied symbolic dynamics techniques to financial returns to determine the Rényi dimensions or the generalized dimensions D_q [415]. There are different symbolization techniques that can convert a financial time series $\{X(i)\}_{i=1}^N$ into a symbol sequence. Li and Wang proposed an approach by considering the speed of price changes, which is quantified by an angle [416, 417],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \theta_i = \arctan[(X(i+\ell)-X(i))/\ell], \label{Eq:MF:Symb:Angle} \nonumber \end{align} \tag{ 155 }$$

which falls within the interval $[-90^{\circ},90^{\circ}]$. To symbolize the raw time series, we partition $[-90^{\circ},90^{\circ}]$ into 2n subintervals, which are usually almost symmetric with respect to $\theta=0^{\circ}$. Li and Wang considered four subintervals $(-90^{\circ},-45^{\circ})$, $[-45^{\circ},0^{\circ})$, $[0^{\circ},45^{\circ})$, and $[45^{\circ},90^{\circ})$ and four corresponding symbols D, d, r and R, representing respectively fast falls, slow falls, slow rises and fast rises of prices.

The symbolization method of Xu and Beck is similar [415]. They suggested to partition the return domain $\Delta{X_i}\in(-1,1)$ into 2n subintervals, mapped into 2n symbols. When n  =  1, the two subintervals are $(-1,0)$ and $[0,\infty)$ and two corresponding symbols d and u, representing respectively falling down and rising of the prices. When n  =  2, the four subintervals are (−1, −  c₁), [−c₁, 0], [0, c₂) and $[c_2,\infty]$, where $c_1\in(0,1)$ and $c_2\in[0,\infty]$ are parameters. Note that Xu and Beck used open intervals, which may discard some data points. In addition, a return cannot be less than  −1. For stock markets with price limit rules, the subintervals should be adjusted accordingly.

For any manner of symbolization, if $\theta_i$ or $\Delta{X}_i$ is an element of a certain subinterval, the ith symbol S_i is determined by the corresponding symbol. In this way, the raw time series is mapped to a symbol sequence {S_i}. The symbol sequence is divided into N_k segments of equal length k. For any given k, we obtain $\omega(k)=(2n){}^k$ allowed symbol subsequences or configurations. The occurrence probability of each configuration C_j of length k can be determined, denoted as $p_j = p(C=C_j)$. Xu and Beck recalled the well-known Rényi information [191]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_{q} =\frac{1}{q-1} \ln \sum^{\omega(k)}_{j=1}p_j^q \label{renin} \nonumber \end{align} \tag{ 156 }$$

and the Rényi dimension

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_q=\lim_{\varepsilon \to 0} \frac{I_q}{\ln \varepsilon}=\lim_{\varepsilon \to 0} \frac{1}{\ln \varepsilon} \frac{1}{q-1} \ln \sum^{\omega(k)}_{j=1}p_{j}^{q}, \nonumber \end{align} \tag{ 157 }$$

where $\varepsilon=(2n){}^{-k}$.

The method has been applied to the daily and 1 min price time series of several US stocks (Alcoa, Bank of America, General Electric, Intel, Johnson & Johnson, Coca Cola and WalMart) from January 1998 to May 2013 and the generalized dimensions were calculated and compared for different values of k [415].

For the integrated measures m_x(s,i) and m_x(s,i) in the ith box of size s, the local singularity strengths $\alpha_x$ and $\alpha_y$ are defined as [418]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m_x(s,i) \sim s^{\alpha_x} ~~~{{{\rm and}}}~~~ m_y(s,i) \sim s^{\alpha_y}. \label{Eq:MF-X-PF:mxy:s:alpha} \nonumber \end{align} \tag{ 158 }$$

Denoting $N_s(\alpha_x,\alpha_y)$ as the number of boxes needed to cover the neighborhoods of points with singularities $\alpha_x$ and $\alpha_y$, the fractal dimension of the set is determined according to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle N_s(\alpha_x,\alpha_y) \sim s^{-f(\alpha_x,\alpha_y)}, \label{Eq:MF-X-PF:Ns:s:f:alpha:xy} \nonumber \end{align} \tag{ 159 }$$

in which $f(\alpha_x,\alpha_y)$ is the joint distribution of the two singularities [305], or the joint multifractal spectrum.

Through the joint partition function [418]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{xy}(\,p,q,s)= \sum_i\left[m_x(s,i)\right]^{\,p/2}[m_y(s,i)]^{q/2}, \label{Eq:MF-X-PF:pq:chi:s} \nonumber \end{align} \tag{ 160 }$$

we can obtain the joint mass exponent function $\tau_{xy}(\,p,q)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{xy}(\,p,q,s) \sim s^{\tau_{xy}(\,p,q)}. \label{Eq:MF-X-PF:pq:tauxy} \nonumber \end{align} \tag{ 161 }$$

Inserting equation (158) into the joint partition function, rewriting the sum into a double integral over $\alpha_x$ and $\alpha_y$, and then applying the steepest descent approach to estimate the integral at small s values, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_{xy}(\,p,q) = \frac{1}{2}(\,p\alpha_x+q\alpha_y) - f(\alpha_x,\alpha_y), \label{Eq:MF-X-PF:tau:alpha:f} \nonumber \end{align} \tag{ 162 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial f(\alpha_x,\alpha_y)}{\partial \alpha_x} = \frac{p}{2} ~~~{{{\rm and}}}~~~ \frac{\partial f(\alpha_x,\alpha_y)}{\partial \alpha_y} = \frac{q}{2}. \label{Eq:MF-X-PF:df:alphaxy:q} \nonumber \end{align} \tag{ 163 }$$

Taking the partial derivative of equation (162) over p  or q, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial\tau_{xy}(\,p,q)}{\partial p} = \frac{\alpha_x}{2} ~~~{{{\rm and}}}~~~ \frac{\partial\tau_{xy}(\,p,q)}{\partial q} = \frac{\alpha_y}{2}. \nonumber \end{align} \tag{ 164 }$$

We can then obtain the double Legendre transforms:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_x(\,p,q) = \frac{2\partial \tau(\,p,q)}{\partial p} ~{{{\rm and}}}~ \alpha_y(\,p,q) = \frac{2\partial \tau(\,p,q)}{\partial q}, \label{Eq:MF-X-PF:alphaxy:tau} \nonumber \end{align} \tag{ 165a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha_x,\alpha_y) = \frac{p\alpha_x(\,p,q)+q\alpha_y(\,p,q)}{2} -\tau_{xy}(\,p,q). \label{Eq:MF-X-PF:f:alpha} \nonumber \end{align} \tag{ 165b }$$

From the canonical perspective, we can obtain the $f(\alpha_x,\alpha_y)$ function directly by defining two canonical measures

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_{xy}(\,p,q,s,i) = \frac{[m_x(s,i)]^{\,p/2}[m_y(s,i)]^{q/2}}{\sum_t [m_x(s,i)]^{\,p/2}[m_y(s,i)]^{q/2}}. \label{Eq:MF-X-PF:pq:muxy} \nonumber \end{align} \tag{ 166 }$$

The two singularity strengths $\alpha_x(\,p)$ and $\alpha_x(\,p)$ and the joint multifractal spectrum f _xy(p ,q) can be computed by linear regressions in log–log scales using the following equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\alpha_x(\,p,q)} = \lim_{s\to0} \frac{\sum_i \mu_{xy}(\,p,q,s,i) \ln{m_x(s,i)}}{\ln{s}}, \label{Eq:MF-X-PF:pq:alphax:mu:mx} \nonumber \end{align} \tag{ 167a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\alpha_y(\,p,q)} = \lim_{s\to0}\frac{\sum_t \mu_{xy}(\,p,q,s,i) \ln{m_y(s,i)}}{\ln{s}}, \label{Eq:MF-X-PF:pq:alphay:mu:my} \nonumber \end{align} \tag{ 167b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {f_{xy}(\,p,q)} = \lim_{s\to0}\frac{\sum_t \mu_{xy}(\,p,q,s,i) \ln\left[\mu_{xy}(\,p,q,s,i)\right]}{\ln{s}}. \label{Eq:MF-X-PF:pq:fpq:mu} \nonumber \end{align} \tag{ 167c }$$

The joint mass exponent function can be obtained by using equation (162).

Xie et al derived analytical results for binomial measures with multipliers p _x and p _y [418]. They found that the two integrated measures m_x(s,i) and m_y(s,i) in boxes of size s  =  2^l are related explicitly by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m_x(s,i) &=C(s) \left[m_y(s,i)\right]^{\beta} \nonumber \\ &= {\rm e}^{-\gamma L}s^{{\gamma}/{\ln 2}} \left[m_y(s,i)\right]^{\beta}, \label{Eq:MF-X-PF:analytic:mxmy} \nonumber \end{align} \tag{ 168 }$$

where 2^L is the length of the time series,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \beta=\frac{\ln p_x-\ln (1-p_x)}{\ln p_y-\ln (1-p_y)}, \label{Eq:MF-X-PF:analytic:beta} \nonumber \end{align} \tag{ 169 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \gamma=\beta\ln (1-p_y) -\ln (1-p_x). \label{Eq:MF-X-PF:analytic:gamma} \nonumber \end{align} \tag{ 170 }$$

When $p_x+p_y=1$, we have $\beta=-1$. When $p_x=p_y$, we have $\beta=1$ and $C(s)=1$. When both p _x and p _y are greater than 0.5 or less than 0.5, that is, $(\,p_x-0.5)(\,p_y-0.5)>0$, we have $\beta>0$; Otherwise, when $(\,p_x-0.5)(\,p_y-0.5)<0$, we have $\beta<0$.

The joint mass exponent function is expressed as follows [418]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_{xy}(\,p,q) &=\frac{p\gamma}{2\ln 2}+\tau_y(Q) \nonumber \\ &=\frac{p\gamma}{2\ln 2}-\frac{\ln[\,p_y^{Q}+(1-p_y)^{Q}]}{\ln{2}} \label{Eq:MF-X-PF:analytic:tauxy:p:q} \nonumber \end{align} \tag{ 171 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Q = \frac{\beta p}{2}+\frac{q}{2} \label{Eq:MF-X-PF:analytic:chi:p:q:c} \nonumber \end{align} \tag{ 172 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_y(Q) =- \ln[\,p_y^{Q}+(1-p_y)^{Q}]/\ln{2}. \label{Eq:MF-X-PF:analytic:tauy} \nonumber \end{align} \tag{ 173 }$$

It follows that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_x = \frac{\gamma}{\ln 2}-\frac{\beta}{\ln 2}\frac{p_y^{Q}\ln p_y+(1-p_y)^{Q}\ln (1-p_y) }{p_y^{Q}+(1-p_y)^{Q}} \label{Eq:MF-X-PF:analytic:alphax} \nonumber \end{align} \tag{ 174 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_y = -\frac{1}{\ln 2}\frac{p_y^{Q}\ln p_y+(1-p_y)^{Q}\ln (1-p_y) }{p_y^{Q}+(1-p_y)^{Q}}. \label{Eq:MF-X-PF:analytic:alphay} \nonumber \end{align} \tag{ 175 }$$

We obtain immediately the relationship between $\alpha_x$ and $\alpha_y$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_x =\frac{\gamma}{\ln 2}+\beta\alpha_y, \label{Eq:MF-X-PF:analytic:alphaxalphay} \nonumber \end{align} \tag{ 176 }$$

which shows that $f_{xy}(\alpha_x,\alpha_y)$ is a curve along this line rather than a surface and the line segment (176) is the projection of $f_{xy}(\alpha_x,\alpha_y)$ onto the $(\alpha_x,\alpha_y)$ plane.

Xie et al obtained that [418]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} \alpha_{y,\min} & = \min\left\{-\frac{\ln p_y}{\ln 2},-\frac{\ln(1-p_y)}{\ln 2}\right\} \nonumber \\ \alpha_{y,\max} & = \max\left\{-\frac{\ln p_y}{\ln 2},-\frac{\ln(1-p_y)}{\ln 2}\right\} \end{array} \right. \label{Eq:MF-X-PF:properties:analytic:max:min:alphay} \nonumber \end{align} \tag{ 177 }$$

such that the width of the singularity spectrum of $\alpha_y$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta \alpha_y =\frac{\left|\ln (1-p_y) - \ln p_y\right|}{\ln 2}. \label{Eq:MF-X-PF:properties:analytic:deltaalphay} \nonumber \end{align} \tag{ 178 }$$

When p _y  =  0.5, $\Delta \alpha_y=0$. In this case, the measure is neither multifractal nor monofractal since it is uniformly distributed on the support. Similarly, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} \alpha_{x,\min} & = \min\left\{-\frac{\ln p_x}{\ln 2},-\frac{\ln(1-p_x)}{\ln 2}\right\} \nonumber \\ \alpha_{x,\max} & = \max\left\{-\frac{\ln p_x}{\ln 2},-\frac{\ln(1-p_x)}{\ln 2}\right\}. \end{array} \right. \label{Eq:MF-X-PF:properties:analytic:max:min:alphax} \nonumber \end{align} \tag{ 179 }$$

Further, the multifractal spectrum $f_{xy}(\alpha_x,\alpha_y)$ can be expressed as [418]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{xy}(\alpha_x,&\alpha_y) =\left\{Q\left(\frac{1-p_y}{p_y}\right)^{Q} \ln\left(\frac{p_y}{1-p_y}\right)\right. \nonumber \\ &+\left.\left[1+\left(\frac{1-p_y}{p_y}\right)^{Q}\right]\ln\left[1+\left(\frac{1-p_y}{p_y}\right)^{Q}\right]\right\} \nonumber \\ & ~~\left\{{\ln 2}{\left[1+\left(\frac{1-p_y}{p_y}\right)^{Q}\right]}\right\}^{-1}. \label{Eq:MF-X-PF:analytic:fpq} \nonumber \end{align} \tag{ 180 }$$

It follows immediately that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{xy}(Q=0) = 1 ~{{\rm and}}~f_{xy}(Q)=f_{xy}(-Q), \nonumber \end{align} \tag{ 181 }$$

where $\newcommand{\tr}{{\rm tr}} f_{xy}(Q)\triangleq f_{xy}(\alpha_x,\alpha_y; Q)$. It implies that $f_{xy}(\alpha_x,\alpha_y)$ is symmetric with respect to the line Q  =  0. In addition, it is found that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{Q\to\pm\infty} f_{xy}(\,p,q) = 0. \label{Eq:MF-X-PF:properties:analytic:fpq:infty} \nonumber \end{align} \tag{ 182 }$$

When Q  <  0, ${{\rm d}}f_{xy}(Q)/{{\rm d}}Q>0$ so that f _xy(Q) is a monotonically increasing function of Q. When Q  >  0, ${{\rm d}}f_{xy}(Q)/{{\rm d}}Q>0$ so that f _xy(Q) is a monotonically decreasing function of Q. Therefore, the maximum of f _xy(Q) is 1 and its minimum is 0.

The multifractal cross-correlation analysis based on statistical moments (MFSMXA) proposed in [419] is actually a special case of MF-X-PF$(\,p,q)$ when p   =  q. We call it MF-X-PF$(q)$ for consistence [418]. In this case, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_{xy}(q) = \frac{q}{2}(\alpha_x+\alpha_y) - f_{xy}(\alpha_x,\alpha_y), \label{Eq:MF-X-PF:q:tau:ax:ay:f} \nonumber \end{align} \tag{ 183 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial f_{xy}(\alpha_x,\alpha_y)}{\partial \alpha_x} = \frac{\partial f_{xy}(\alpha_x,\alpha_y)}{\partial \alpha_y} = \frac{q}{2}. \label{Eq:MF-X-PF:q:df:dalphax:q} \nonumber \end{align} \tag{ 184 }$$

Taking the derivative of equation (183) over q and applying equation (184), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{{\rm d}}\tau_{xy}(q)}{{{\rm d}}q} = \frac{\alpha_x+\alpha_y}{2}. \label{Eq:MF-X-PF:q:dtau:dq:alphax:alphay} \nonumber \end{align} \tag{ 185 }$$

Defining that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \alpha_{xy} \triangleq \frac{1}{2}[\alpha_{x}(q) + \alpha_{y}(q)], \label{Eq:MF-X-PF:q:alpha:alphax:alphay} \nonumber \end{align} \tag{ 186 }$$

we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_{xy} = \frac{{{\rm d}}\tau_{xy}(q)}{{{\rm d}}q}, \label{Eq:MF-X-PF:q:dtau:dq:alpha} \nonumber \end{align} \tag{ 187a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle f_{xy}(\alpha_{xy}(q)) \triangleq f_{xy}(\alpha_x,\alpha_y) = q \alpha_{xy}(q) -\tau_{xy}(q), \label{Eq:MF-X-PF:q:f:alpha} \nonumber \end{align} \tag{ 187b }$$

which is the Legendre transform.

It is derived that [418]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_{xy}(q)= \frac{1}{2}[\tau_x(q) + \tau_y(q)], \label{Eq:MF-X-PF:q:tauxy:taux:tauy} \nonumber \end{align} \tag{ 188 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{xy}(q)= \frac{1}{2}[\,f_x(q) + f_y(q)]. \label{Eq:MF-X-PF:q:fxy:fx:fy} \nonumber \end{align} \tag{ 189 }$$

We define the generalized joint Hurst exponent H_xy(q) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_{xy}(q)= qH_{xy}(q)-1. \label{Eq:MF-X-PF:tauxy:Hxy} \nonumber \end{align} \tag{ 190 }$$

It follows that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{xy}(q)= \frac{1}{2}[H_x(q)+H_y(q)], \label{Eq:MF-X-PF:Hxy:Hx:Hy} \nonumber \end{align} \tag{ 191 }$$

where H_x(q) and H_y(q) are the generalized Hurst exponent of m_x and m_y. These relations were observed numerically using the MF-X-DFA method [180], the MF-X-DMA method [144] and the MF-X-PF method with p   =  q [419].

If we can define the canonical measures

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_{xy}(q,s,i) = \frac{[m_x(s,i)m_y(s,i)]^{q/2}}{\sum_i [m_x(s,i)m_y(s,i)]^{q/2}} \label{Eq:MF-X-PF:q:mu:p:q} \nonumber \end{align} \tag{ 192 }$$

the two singularity strengths $\alpha_x(\,p)$ and $\alpha_x(\,p)$ and the joint multifractal spectrum f _xy(p ,q) can be computed directly:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_{xy}(q) &= \lim_{s\to0} \frac{\sum_i \mu_{xy}(q,s,i) \ln{[m_x(s,i)m_y(s,i)]^{1/2}}}{\ln{s}}\nonumber \\ &= \frac{\alpha_{x}(q)+\alpha_{y}(q)}{2}, \label{Eq:MF-X-PF:q:alpha:mu:mx} \nonumber \end{align} \tag{ 193a }$$

where equations (167a) and (167b) are used in the second equality, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{xy}(\alpha_{xy}(q)) = \lim_{s\to0}\frac{\sum_i \mu_{xy}(q,s,i) \ln\left[\mu_{xy}(q,s,i)\right]}{\ln{s}}. \label{Eq:MF-X-PF:q:fpq:mu} \nonumber \end{align} \tag{ 193b }$$

The joint mass exponent function can be obtained by using equation (187b).

Xiong and Shang proposed the weighted MF-X-PF$(q)$ method (W-MFSMXA) [429]. Their numerical experiments on two-exponent ARFIMA processes, binomial measures and NBVP time series illustrate that W-MFSMXA slightly outperforms MF-X-PF$(q)$ for relatively short time series.

The joint multifractal analysis based on the structure function approach, termed MF-X-SF, also has its root in turbulence, in which the joint structure function of the velocity and temperature fluctuations was introduced to study the cross-correlation between the temperature and velocity dissipation fields in a heated turbulent jet [430, 431]. For two self-similar time series $X(i)$ and $Y(i)$ with $i=1,2,\cdots,N$, the $(\,p,q)$-order joint structure function is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{xy}(\,p,& q,s) = \left\langle \Delta{X(i,s)}^{\,p/2} \Delta{Y(i,s)}^{q/2} \right\rangle \nonumber \\ &=\frac{1}{N-s+1}\sum_{i=s}^{N} \Delta{X(i,s)}^{\,p/2} \Delta{Y(i,s)}^{q/2}, \label{Eq:MF-X-SF:Kxy:pq} \nonumber \end{align} \tag{ 194 }$$

where $\Delta{X(i,s)}=X(i)-X(i-s)$ and $\Delta{Y(i,s)} \,=$ $Y(i)-Y(i-s)$. The special case of p   =  1 is called by some researchers as the analogous multifractal cross-correlation analysis (AMF-XA) [432], in which the authors were seemingly unaware of previous works [430, 433]. Because the innovations $\Delta{X}$ and $\Delta{Y}$ could be negative, the structure functions K_xy(p ,q,s) are usually defined on positive integers of p  and q. Therefore, it is more convenient to use absolute innovations to define the joint structure function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{xy}&(\,p,q,s) = \left\langle \left|\Delta{X(i,s)}\right|^{\,p/2} |\Delta{Y(i,s)}|^{q/2} \right\rangle \nonumber \\ &=\frac{1}{N-s+1}\sum_{i=s}^{N} |\Delta{X(i,s)}|^{\,p/2} |\Delta{Y(i,s)}|^{q/2}, \label{Eq:MF-X-SF:Kxy:pq:absolute} \nonumber \end{align} \tag{ 195 }$$

which is widely adopted in the structure function approach, as well as in the multifractal height cross-correlation analysis (MFHXA) [433]. If there is a scaling relation between the joint structure function K_xy(p ,q,s) and the scale s, we can define the scaling exponent function $\tau_{xy}(\,p,q)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{xy}(\,p,q,s) \sim s^{\zeta_{xy}(\,p,q)} = s^{\tau_{xy}(\,p,q)+1}, \label{Eq:MF-X-SF:Kxy:pq:s} \nonumber \end{align} \tag{ 196 }$$

which is an analog of the MF-X-PF$(\,p,q)$ approach.

We further define the singularity functions $\alpha_x(\,p,q)$ and $\alpha_y(\,p,q)$ as the partial derivatives of $\tau_{xy}(\,p,q)$ with respect to p  and q, respectively, as in equation (165a). The singularity spectrum function $\newcommand{\tr}{{\rm tr}} f_{xy}(\,p,q)\triangleq f_{xy}(\alpha_x,\alpha_y)$ can be obtained as in equation (165b). We note that it is hard to define a generalized Hurst index H_xy(p , q) in a consistent way. The joint structure function approach has been considered for the multifractal random walks [434].

A special case of the MF-X-SF$(\,p,q)$ approach when p   =  q was proposed by Kristoufek, which was termed the multifractal height cross-correlation analysis [433]. For consistency, we call it the MF-X-SF$(q)$ approach. In this case, the qth order joint structure function is defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{xy}(q,s) &= \left\langle \left|\Delta{X(i,s)}\Delta{Y(i,s)}\right|^{q/2} \right\rangle \nonumber \\ &=\frac{1}{N-s+1}\sum_{i=s}^{N} |\Delta{X(i,s)}\Delta{Y(i,s)}|^{q/2}. \label{Eq:MF-X-SF:Kxy:q} \nonumber \end{align} \tag{ 197 }$$

If there is a scaling relation between the joint structure function K_xy(p , q, s) and the scale s, we can define the scaling exponent function $\tau_{xy}(q)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{xy}(q,s) \sim s^{\tau_{xy}(q)+1} = s^{qH_{xy}(q)}, \label{Eq:MF-X-SF:Kxy:q:s} \nonumber \end{align} \tag{ 198 }$$

where H_xy(q) are the generalized bivariate Hurst exponents or the generalized joint Hurst exponents. Under certain conditions, we have [433]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{xy}(q) = \frac{H_x(q)+H_y(q)}{2}. \label{Eq:MF-X:Hxy:Hx:Hy} \nonumber \end{align} \tag{ 199 }$$

We define the joint singularity function and its spectrum through the Legendre transform:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_{xy} = \frac{{{{\rm d}}}\tau_{xy}(q)}{{{{\rm d}}}q}, \label{Eq:MF-X-SF:q:dtau:dq:alpha} \nonumber \end{align} \tag{ 200a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{xy}(\alpha_{xy}) = q \alpha_{xy} -\tau_{xy}. \label{Eq:MF-X-SF:q:f:alpha} \nonumber \end{align} \tag{ 200b }$$

If equation (199) holds, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_{xy}(q) = \frac{\tau_x(q)+\tau_y(q)}{2} \label{Eq:MF-X:tayxy:taux:tauy} \nonumber \end{align} \tag{ 201 }$$

according to the definitions that $\tau_x(q)=qH_x(q)-1$, $\tau_y(q)=qH_y(q)-1$ and $\tau_{xy}(q)=qH_{xy}(q)-1$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_{xy}(q) = \frac{\alpha_x(q)+\alpha_y(q)}{2} \label{Eq:MF-X:alphaxy:alphax:alphay} \nonumber \end{align} \tag{ 202 }$$

according to the definitions that $\alpha_x(q)={{{\rm d}}}\tau_x(q)/{{{\rm d}}}q$, $\alpha_y(q)={{{\rm d}}}\tau_y(q)/{{{\rm d}}}q$ and $\alpha_{xy}(q)={{{\rm d}}}\tau_{xy}(q)/{{{\rm d}}}q$. It follows from equation (200b) that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{xy}(\alpha_{xy}) = \frac{f_{x}(\alpha_x)+f_y(\alpha_y)}{2}. \label{Eq:MF-X:fxy:fx:fy} \nonumber \end{align} \tag{ 203 }$$

The cross wavelet transform (XWT) of two time series $\{X_i\}_{i=1}^N$ and $\{Y_i\}_{i=1}^N$ was introduced by Hudgins et al to study the highly intermittent pattern in atmospheric turbulence [435], which is defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal{W}}_{xy}(s,i) = {\mathcal{W}}_x(s,i){\mathcal{W}}_y(s,i). \label{Eq:MF-X-WT:Wxy} \nonumber \end{align} \tag{ 204 }$$

The cross-wavelet power is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{xy}(s) = \|{\mathcal{W}}_{xy}(s,i)\|. \label{Eq:MF-X-WT:PWxy} \nonumber \end{align} \tag{ 205 }$$

The wavelet coherence $\rho_{xy}^2(s)$ is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_{xy}^2(s)& = \frac{\|{\mathcal{W}}_{xy}(s,i)\|}{\|{\mathcal{W}}_x(s,i)\|\|{\mathcal{W}}_y(s,i)\|} \nonumber \\ & = \frac{\sum_{i=1}^N|W_x(s,i)W_y(s,i)|}{\left(\sum_{i=1}^NW_x^2(s,i)\sum_{i=1}^NW_y^2(s,i)\right)^{1/2}}. \label{Eq:MF-X-WT:WC} \nonumber \end{align} \tag{ 206 }$$

The wavelet coherence is able to uncover the co-movement between two time series in the time-frequency domain, which is widely applied in finance and economics [436–438].

Similar to the MF-X-PF method, Jiang et al designed the multifractal cross wavelet analysis [439]. It generalizes the concept of wavelet coherence and constructs the following joint partition function:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{xy}(\,p,q,s)= \sum_{i=1}^N|W_x(s,i)|^{\,p/2}|W_y(s,i)|^{q/2}. \label{Eq:MF-X-WT:Chi:p:q:s} \nonumber \end{align} \tag{ 207 }$$

The definition of the partition function allows us to uncover the more intricate relationship between the coherency and the scale under different scopes, which corresponds to a cross-multifractal behavior. If the underlying processes are jointly multifractal, the result is a scaling behavior:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{xy}(\,p,q,s) \sim s^{{\mathcal{T}}_{xy}(\,p, q)}, \label{Eq:MF-X-WT:Chi:Scaling} \nonumber \end{align} \tag{ 208 }$$

where ${\mathcal{T}}_{xy}(\,p, q)$ is the joint mass exponent function, which can be estimated by regressing $\ln \chi_{xy} (\,p,q,s)$ against $\ln s$ in the scaling range for a given pair $(\,p,q)$.

Analogous to the double Legendre transform in the joint multifractal analysis based on the partition function approach MF-X-PF$(\,p,q)$ [418], we define the joint singularity strength function h_x and h_y

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_x(\,p,q) = \frac{2\partial {\mathcal{T}}_{xy}(\,p, q)}{\partial p}, \label{Eq:MF-X-WT:hx} \nonumber \end{align} \tag{ 209 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_y(\,p,q) = \frac{2\partial {\mathcal{T}}_{xy}(\,p, q)}{\partial q}, \label{Eq:MF-X-WT:hy} \nonumber \end{align} \tag{ 210 }$$

and the multifractal spectrum $D_{xy}(h_x, h_y)$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_{xy}(h_x, h_y) = \frac{1}{2}(\,p h_x + q h_y) - {\mathcal{T}}_{xy}. \label{Eq:MF-X-WT:Dh} \nonumber \end{align} \tag{ 211 }$$

When p   =  0 and q  =  0, the joint partition function in equation (207) is equal to the number of wavelet coefficients, or equivalently the total number of data points of the original series. In other words, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{xy}(0,0,s) = N, \label{Eq:MF-X-WT:Chi:0:0} \nonumber \end{align} \tag{ 212 }$$

which results in

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal{T}}_{xy}(0, 0) = 0. \label{Eq:MF-X-WT:Txy:0:0} \nonumber \end{align} \tag{ 213 }$$

In contrast, in the MF-X-PF$(\,p, q)$ method, $\tau_{xy} (0,0)=-1$. Therefore,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal{T}}_{xy}(\,p, q) \neq \tau_{xy}(\,p,q). \label{Eq:MF-X-WT:Txy:neq:tau} \nonumber \end{align} \tag{ 214 }$$

It follows immediately that h_x(p , q), h_y(p , q) and $D_{xy}(h_x, h_y)$ obtained from the MF-X-WT$(\,p,q)$ method are different from the joint singularity strengths $\alpha_x(\,p,q)$, $\alpha_y(\,p,q)$ and the joint multifractal spectrum $f_{xy}(\alpha_x, \alpha_y)$ obtained from the MF-X-PF$(\,p,q)$ method.

Analogous to the partition function approach [197, 198], we can directly estimate the joint singularity strength h_x and h_y and the joint multifractal spectrum D_xy(p , q) using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_x(\,p,q) = \lim_{s \rightarrow 0} \sum_i \frac{\mu_{xy}(\,p,q,s,i) \ln |w_x(s,i)|}{\ln s}, \label{Eq:MF-X-WT:DE:hx} \nonumber \end{align} \tag{ 215 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_y(\,p,q) = \lim_{s \rightarrow 0} \sum_i \frac{\mu_{xy}(\,p,q,s,i) \ln |w_y(s,i)|}{\ln s}, \label{Eq:MF-X-WT:DE:hy} \nonumber \end{align} \tag{ 216 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_{xy}(\,p,q) = \lim_{s \rightarrow 0} \sum_i\frac{\mu_{xy}(\,p,q,s,i) \ln \mu_{xy}(\,p,q,s,i)}{\ln s} \label{Eq:MF-X-WT:DE:Dxy} \nonumber \end{align} \tag{ 217 }$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \mu_{xy} (\,p,q,s,i) = \frac{|w_x(s,i)|^{\,p/2} |w_y(s,i)|^{q/2}}{\chi_{xy}(\,p,q,s)}. \nonumber \end{align*}$$

Thus we can directly determine the joint singularity strength functions h_x(p , q) and h_y(p , q) and the joint multifractal function D_xy(p , q) from equations (215)–(217) through linear regressions.

Extensive numerical experiments on two multifractal binomial measures unveil connections between the theoretical joint multifractal formulas of MF-X-PF$(\,p,q)$ derived in [418] and the empirical joint multifractal features of MF-X-WT$(\,p,q)$ [439]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathcal{T}}_{xy}(\,p, q) + \frac{p}{2} + \frac{q}{2}-1 = \tau_{xy}(\,p, q), \label{Eq:pModel:WTTau:PFTau} \nonumber \end{align} \tag{ 218 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_x(\,p, q) + 1 = \alpha_x(\,p, q), \label{Eq:pModel:WThx:PFalphax} \nonumber \end{align} \tag{ 219 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_y(\,p, q) + 1 = \alpha_y(\,p, q), \label{Eq:pModel:WThy:PFalphay} \nonumber \end{align} \tag{ 220 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_{xy} (h_x, h_y) + 1 = f_{xy}(\alpha_x, \alpha_y), \label{Eq:pModel:WTDh:PFfalpha} \nonumber \end{align} \tag{ 221 }$$

where $\tau_{xy}(\,p,q)$, $\alpha_x(\,p,q)$, $\alpha_y(\,p,q)$ and $f_{xy}(\alpha_x, \alpha_y)$ are given in equations (171), (174), (175) and (180).

Although the MF-X-WT$(\,p,q)$ method is able to successfully identify the monofractal or multifractal cross correlations in time series, the results show non-negligible deviations from the theoretical results [439]. Nevertheless, empirical applications to financial returns and financial volatilities uncover evident cross-multifractality [439].

The implementation of a joint multifractal analysis based on WTMM is not straightforward. We can conveniently determine the two sets of modulus maxima lines ${\mathcal{L}}_x$ and ${\mathcal{L}}_y$ for $X(i)$ and $Y(i)$. However, there is no one-to-one correspondence between ${\mathcal{L}}_x$ and ${\mathcal{L}}_y$ because there are different numbers of modulus maxima for the two time series at a given scale. Lin and Sharif found a nice solution to overcome this difficulty by properly pairing up the modulus maxima lines in ${\mathcal{L}}_x$ and ${\mathcal{L}}_y$ [440, 441]. Two modulus maxima lines in ${\mathcal{L}}_x$ and ${\mathcal{L}}_y$ are paired up if they are close to each other in the time coordinate. The joint partition function can thus be defined as [440]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z_{p,q}(s_j)= \sum_{k=1}^{N_j}|W_x(i_{j,k},s_j)|^{\,p}|W_y(i_{j,k},s_j)|^q , \label{Eq:MF-X-WTMM:Zpq} \nonumber \end{align} \tag{ 222 }$$

which scales as a power law with respect to the scale s:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z_{p,q}(s) \sim s^{\tau_{xy}(\,p,q)}. \label{Eq:MF-X-WTMM:Zpq:s:Tau} \nonumber \end{align} \tag{ 223 }$$

Note that we use p  and q instead of p /2 and q/2 following [440]. In this case, the Legendre transform is expressed as [440]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\alpha_x}=\frac{\partial\tau(\,p,q)}{\partial p},~~ {\alpha_y}=\frac{\partial\tau(\,p,q)}{\partial q}, \label{Eq:MF-X-WTMM:alpha:dtau:dq} \nonumber \end{align} \tag{ 224a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f({\alpha_x},{\alpha_y})=p{\alpha_x}(\,p,q) +q{\alpha_y}(\,p,q)-\tau(\,p,q), \label{Eq:MF-X-WTMM:f:alpha} \nonumber \end{align} \tag{ 224b }$$

where $p=\partial{f}/\partial{\alpha_x}$ and $q=\partial{f}/\partial{\alpha_y}$.

The validity of the MF-X-WTMM method is confirmed by numerical experiments on coupled random binomial cascades [440]. Coupled random binomial cascades are characterized by a coupling strength parameter $g\in[0,1]$, which g  =  0 and g  =  1 stand respectively for uncoupled and fully coupled time series [305]. Their numerical result that the $f(\alpha_x, \alpha_y)$ spectrum collapses into a one-dimensional curve in the fully coupled case is in excellent agreement with the analytical result of equation (176) [418].

Due to the inherent drawbacks associated with the multifractal wavelet transform analysis, the MF-X-WT$(\,p,q)$ is not recommended for multifractal cross-correlation analysis. In addition to the MF-X-WTMM method, a suitable alternative candidate is the wavelet leaders [295–300]. Jiang et al generalized the MF-WL method to the joint multifractal analysis based on wavelet leaders with two moment orders p  and q, termed MF-X-WL$(\,p, q)$ [442].

One first obtains the wavelet leaders $\newcommand{\e}{{\rm e}} \ell_x(\,j,k)$ and $\newcommand{\e}{{\rm e}} \ell_y(\,j,k)$ of $\{X_i\}_1^N$ and $\{Y_i\}_1^N$ at different scales s  =  2^j, and then defines the joint moment with orders p  and q as [442]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{xy}(\,p, q, j) = \frac{1}{n_j}\sum_{k=1}^{n_j} \ell_x(\,j,k)^{\,p/2}\ell_y(\,j,k)^{q/2}, \label{Eq:MF-X-WL:SF} \nonumber \end{align} \tag{ 225 }$$

where n_j is the number of wavelet leaders at scale s  =  2^j. When X  =  Y and p   =  q, we recover the traditional multifractal formalism based on wavelet leaders presented in section 2.3.3. Further scaling analyses follow the conventional way as other similar methods.

Through numerical experiments, it is found that the MF-X-WL$(\,p,q)$ method can successfully detect joint multifractality in binomial measures and monofractality in bivariate fractional Brownian motions [442]. The method has also been used to investigate the joint multifractality between the daily returns of the DJIA and NASDAQ indices, and between their daily volatilities as well. The MF-X-WL$(\,p,q)$ method works for both positive and negative $(\,p,q)$ values. However, the method seemingly cannot provide very accurate results.

Multiplicative cascade models generate multinomial measures that are constructed in a recursive manner [162]. Multiplicative cascades exist in financial markets since there is evidence of a causal information cascade from large scales to fine scales [478, 479], which is quite similar to turbulence. It is from this insight that the MRW and many other models emerged. The multiplicative cascade process starts with an interval $[0,1]$ on with the measure is uniformly distributed. Without loss of generality, assume that the total measure on the interval is 1. At the first iteration, the interval is divided into b segments whose lengths are $\{s_i:i=1,\cdots,b\}$ and the measure is redistributed on each segment such that the measure on the ith segment is m_i. At the second iteration, each segment is divided into b smaller segments and the measure on it is redistributed in the same way. This procedure goes to infinity, resulting in a multiscale multinomial measure. The mass exponent function can be solved from the following equation [171]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{i=1}^b \frac{m_i^q}{s_i^{\tau(q)}} = 1. \label{Eq:MF-PF:Gamma} \nonumber \end{align} \tag{ 249 }$$

In this iterative process, the daughter segments are not necessary to cover the mother segment. Therefore, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{i=1}^b s_i \leqslant 1 ~~~{{\rm and}}~~~~ \sum_{i=1}^b m_i = 1. \nonumber \end{align} \tag{ 250 }$$

It can be proven that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}l@{}} {{{\rm d}}}\tau(q)/{{{\rm d}}}q>0 \nonumber \\ {{{\rm d}}}^2\tau(q)/{{{\rm d}}}q^2\leqslant0 \end{array}\right., \nonumber \end{align} \tag{ 251 }$$

where the equality '=' holds if and only if all $\ln{m_i}/\ln{s_i}$ terms are identical. In this situation, the measure is uniformly distributed on the support and is monofractal such that $\tau(q)$ is a linear function of q. According to equation (25), we have $\newcommand{\e}{{\rm e}} D_q \equiv D_0 = D_f$, indicating that D_q is independent of q. Since $\alpha(q)=\tau'(q)$, we have $\alpha'(q)= \tau''(q)= 0$, suggesting that the singularity width $\alpha_{\max}-\alpha_{\min}=0$.

It is easy to verify that, when q  =  0 and q  =  1,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(0) = -1 ~{\mathrm{and}}~ \tau(1)=0 \label{Eq:MF-PF:tau:0} \nonumber \end{align} \tag{ 252 }$$

are respectively solutions to equation (249). Because $\tau(q)$ is a monotonically increasing function of q for multifractals, these solutions are the unique solutions. For the D_q function, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{{\rm d}}}D_q/{{{\rm d}}}q\leqslant0, \label{Eq:MF-PF:dDq:dq} \nonumber \end{align} \tag{ 253 }$$

where the equality holds if and only if all $\ln{m_i}/\ln{s_i}$ terms are identical. In addition, the limits $D_{+\infty}$ and $D_{- \infty}$ exist when $q\to\pm\infty$, which read

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle D_{+\infty}\triangleq\lim_{q\to+\infty}D_q=\min_i\left\{\ln{m_i}/\ln{s_i}\right\} \label{Eq:MF-PF:Dq:infty:p} \nonumber \end{align} \tag{ 254a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle D_{-\infty}\triangleq\lim_{q\to-\infty}D_q=\max_i\left\{\ln{m_i}/\ln{s_i}\right\}. \label{Eq:MF-PF:Dq:infty:n} \nonumber \end{align} \tag{ 254b }$$

For the $\alpha(q)$ function, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rm d}}\alpha(q)/{{\rm d}}q\leqslant0, \label{Eq:MF-PF:dalpha:dq} \nonumber \end{align} \tag{ 255 }$$

where the equality holds if and only if all $\ln{m_i}/\ln{s_i}$ terms are identical. In addition, the limits $\alpha(+\infty)$ and $\alpha(- \infty)$ exist when $q\to\pm\infty$, which read [162]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \alpha_{\min}=\alpha_{+\infty}\triangleq\lim_{q\to+\infty}\alpha(q)=\min_i\left\{\ln{m_i}/\ln{s_i}\right\} \label{Eq:MF-PF:alpha:infty:p} \nonumber \end{align} \tag{ 256a }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \alpha_{\max}=\alpha_{-\infty}\triangleq\lim_{q\to-\infty}\alpha(q)=\max_i\left\{\ln{m_i}/\ln{s_i}\right\}. \label{Eq:MF-PF:alpha:infty:n} \nonumber \end{align} \tag{ 256b }$$

Note that the $\alpha(q)$ function has a very similar shape to the D_q function. For the $f(\alpha(q))$ function, the relationship ${{\rm d}}f(\alpha)/{{\rm d}}\alpha=q$ still holds and $f_{\max} = f(q=0) =1$. Moreover, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}l@{}} f(\alpha(q))\geqslant0 \nonumber \\ {{\rm d}}^2f(\alpha)/{{\rm d}}\alpha^2\leqslant0 \end{array} \right., \nonumber \end{align} \tag{ 257 }$$

where the equality in the second formula holds if and only if all $\ln{m_i}/\ln{s_i}$ terms are identical. The proofs of these properties can be found in [82] and references therein. The multifractal measures discussed here are conservative. Non-conservative multifractal measures exhibit different properties [480].

When s_i  =  1/b for all i's, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = -\frac{\ln\sum_{i=1}^b m_i^q}{\ln{b}}. \label{Eq:MF:Model:Deterministic:Multinomial:tau} \nonumber \end{align} \tag{ 258 }$$

The generalized dimension function D_q and the generalized Hurst exponent function $H(q)$ can be explicitly expressed due to equations (12) and (13). According to equation (11a), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \alpha(q) \triangleq \frac{{{\rm d}}\tau(q)}{{{\rm d}}q}= -\frac{\sum_{i=1}^b m_i^q\ln{m_i}}{\sum_{i=1}^b m_i^q\ln{b}}. \label{Eq:MF:Model:Deterministic:Multinomial:alpha} \nonumber \end{align} \tag{ 259 }$$

The multifractal spectrum $f(\alpha)$ can also be explicitly expressed from the parameter q using equation (11b).

The most widely used multifractal model in econophysics is the p -model that generates multifractal binomial measures in a recursive way [361], It is also frequently used for testing the performance of different multifractal methods. The p -model is a special case of the deterministic multiplicative cascade model with b  =  2 and $s_1=s_2=1/2$. To construct the binomial measure, one starts with the zeroth iteration k  =  0, where the data set $z(i)$ consists of one value, z⁽⁰⁾(1)  =  1. In the kth iteration, the data set $\{z^{(k)}(i): i = 1, 2, \cdots, 2^k\}$ is obtained from

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}rcl@{}} z^{(k)}(2i-1)&=& p_z z^{(k-1)}(i) \nonumber \\ z^{(k)}(2i) &=& (1-p_z)z^{(k-1)}(i) \end{array} \label{Eq:pModel} \nonumber \end{align} \tag{ 260 }$$

for $i = 1, 2, \cdots, 2^{k-1}$. An early construction of binomial measures can be attributed to de Wijs in 1951 [163, 481, 482].

When $k\to\infty$, z^(k)(i) approaches a binomial measure. According to equation (249), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = -\ln[m_1^q+(1-m_1)^q]/\ln{2}. \label{Eq:pModel:tau} \nonumber \end{align} \tag{ 261 }$$

The analytical expressions of D_q and $H(q)$ can be obtained using equations (12) and (13) respectively. According to the Legendre transform expressed in equations (11a) and (11b), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(q) = -\frac{m_1^q\ln{m_1}+(1-m_1)^q\ln(1-m_1)}{[m_1^q+(1-m_1)^q]\ln{2}} \label{Eq:pModel:alpha} \nonumber \end{align} \tag{ 262 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) &= -q\frac{m_1^q\ln{m_1}+(1-m_1)^q\ln(1-m_1)}{[m_1^q+(1-m_1)^q]\ln{2}} \nonumber \\ &~~+ \frac{\ln[m_1^q+(1-m_1)^q]}{\ln{2}}\nonumber \\ &= -\frac{\alpha_{\max}-\alpha}{\alpha_{\max}-\alpha_{\min}}\log_2\left(\frac{\alpha_{\max}-\alpha}{\alpha_{\max}-\alpha_{\min}}\right)\nonumber \\ &~~-\frac{\alpha-\alpha_{\min}}{\alpha_{\max}-\alpha_{\min}}\log_2\left(\frac{\alpha-\alpha_{\min}}{\alpha_{\max}-\alpha_{\min}}\right). \label{Eq:pModel:f:alpha:2} \nonumber \end{align} \tag{ 263 }$$

Equation (263) was derived by Calvet et al [483].

We now perform a multifractal analysis of a binomial measure with m₁  =  0.3 and length N  =  2²⁰. Figure 8(a) presents the first 1000 data points of the measure. The partition function approach and the direct determination approach are utilized and compared. We first work on the partition function approach. Figure 8(b) illustrates the power-law dependence of $\chi(q,s)$ on scale s for different q values from  −10 to 10. All the curves for different q values have nice power-law scaling in a broad scaling range. The slopes obtained by linear regressions of $\ln\chi(q,s)$ against $\ln{s}$ are estimates of $\tau(q)$, which are shown in figure 8(e). We obtain D_q using equation (25), $\alpha(q)$, $f(q)$ and $f(\alpha)$ using the Legendre transform, which are presented in figures 8(f)–(i). Based on the direct determination approach, figure 8(c) plots the dependence of $\sum_i \mu(q,s,i) \ln{[m(s,i)]}$ against s and figure 8(d) the dependence of $\sum_i\mu(q,s,i)\ln[\mu(q,s,i)]$ against s in linear-log coordinates, where $m(s,i)$ is the total measure in the ith box of size s and $\mu(q,s,i)$ is the corresponding canonical measure. We also observe nice linearity. The slopes of the linear fits in figures 8(c) and (d) are estimates of $\alpha(q)$ and $f(q)$ respectively, which are also shown in figures 8(g) and (h). We obtain $\tau(q)$ by the Legendre transform and D_q according to equation (25), which are depicted in figures 8(e) and (f), respectively. In figure 8(e), we also mark $\tau(0)=-1$ and $\tau(1)=0$, while in figure 8(i), we show that $f'(\alpha)|_{q=0}=0$ and $f'(\alpha)|_{q=1}=1$. In figures 8(e)–(i), we also show the corresponding analytical functions. It is evident that both the classical partition function approach and the direct determination approach unveil perfectly the multifractal nature of the binomial measure with very high accuracy. Note that the choice of s values crucially impacts the results, because of the log-periodic oscillations in deterministic binomial measures [312], which will be discussed in section 5.

Zoom In Zoom Out Reset image size

We can introduce randomness into multiplicative processes to construct stochastic multinomial models [162, 165, 166, 484]. We focus on the uni-scale case that a segment is divided into b sub-segments of identical length. At the first step, a unity segment of measure $\mu_0=1$ is divided into b non-overlapping subsegments. With probability $p'_i$, measures $m'_{i,1},\cdots,m'_{i,b}$ are assigned to the subsegments. There are k probabilities such that the construction matrix is expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bf{CRM}= \left( \begin{array}{@{}ccc|c@{}} m'_{1,1} &\cdots & m'_{1,b} & p'_1 \nonumber \\ \cdots &\cdots &\cdots &\cdots \nonumber \\ m'_{k,1} &\cdots & m'_{k,b} & p'_k \nonumber \\ \end{array} \right). \label{Eq:SDMF:CRM1} \nonumber \end{align} \tag{ 264 }$$

At the second step, each of the b sub-segments is divided into b segments of length 1/b². For the ith sub-segment, with probability $p_j'$, the measures are $m'_{i,1}m'_{j,1},\cdots,m'_{i,b}m'_{j,b}$. This process goes to infinity, resulting in stochastic multinomial measures.

Let $m_{(i-1)b+j}=m'_{ij}$, $p_{(i-1)b+j}=p'_i$ and n  =  kb, we rewrite equation (264) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:SDMF:CRM2} \bf{CRM}= \left( \begin{array}{@{}ccccccc@{}} m_1 & \cdots & m_{(i-1)b+j} & \cdots & m_n \nonumber \\ p_1 & \cdots & p_{(i-1)b+j} & \cdots & p_n \nonumber \\ \end{array} \right). \nonumber \end{align} \tag{ 265 }$$

Note that $\sum_{i=1}^{n}p_i=b$. The mass exponent function $\tau(q)$ is the solution of the following equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:SDMF:Gamma} \Gamma(q,\tau)=\sum_{i=1}^n p_i {m_{ij}^q}{n^{\tau}} = 1. \nonumber \end{align} \tag{ 266 }$$

It follows immediately that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = -\frac{\ln\sum_{i=1}^n p_im_i^q}{\ln{n}}. \label{Eq:MF:Model:Stochastic:Multinomial:tau} \nonumber \end{align} \tag{ 267 }$$

The generalized dimension function D_q and the generalized Hurst exponent function $H(q)$ can be explicitly expressed from equations (12) and (13). According to equation (11a), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \alpha(q) \triangleq \frac{{\rm d}\tau(q)}{{\rm d}q}= -\frac{\sum_{i=1}^np_im_i^q\ln{m_i}}{\sum_{i=1}^n p_im_i^q\ln{n}}. \label{Eq:MF:Model:Stochastic:Multinomial:alpha} \nonumber \end{align} \tag{ 268 }$$

Hence, all the multifractal quantities have analytical expressions.

In figure 9, we illustrate an example of a stochastic binomial measure, together with its multifractal properties. An interesting feature is the presence of negative dimensions $f(\alpha)<0$. We will come back to discuss this issue later in section 5.2.3.

Zoom In Zoom Out Reset image size

When the multipliers have continuous distributions, we obtain a large class of stochastic multiplicative cascade models [162]. The well-known lognormal model originated in fluid mechanics [172, 485]. For lognormal binomial measures [143, 483, 486], the multipliers follow the lognormal distribution in which the logarithmic multipliers have mean $\lambda$ and variance $\sigma^2$. The mass exponent function is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) =-\frac{\sigma^2\ln2}{2}q^2+\lambda q -1 \nonumber \end{align} \tag{ 269 }$$

and the singularity spectrum is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) = 1-\frac{(\alpha-\lambda)^2}{2\ln2\sigma^2}. \nonumber \end{align} \tag{ 270 }$$

For log-Poisson binomial measures [143, 483, 486], the multipliers follows the Poisson distribution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(m) = \frac{{{\rm e}}^{-\gamma}\gamma^m}{m!} \nonumber \end{align} \tag{ 271 }$$

and the singularity spectrum is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) = 1-\frac{\gamma}{\ln2}+\alpha\frac{\ln(\gamma{{\rm e}}/\alpha)}{\ln2}. \nonumber \end{align} \tag{ 272 }$$

For log-Gamma binomial measures [143, 483, 486], the multipliers follows the Gamma distribution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p(m) = \beta^{\gamma}m^{\gamma-1}{{\rm e}}^{-\beta{m}}/\Gamma(\gamma), \nonumber \end{align} \tag{ 273 }$$

where $\beta>0$ and $\gamma>0$, and the singularity spectrum is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) = 1+\gamma\frac{\ln(\alpha\beta/\gamma)}{\ln2}+\frac{\gamma-\alpha\beta}{\ln2}. \nonumber \end{align} \tag{ 274 }$$

Similar random cascading processes such as the $\mathcal{W}$-cascades on the dyadic trees of their orthogonal wavelet coefficients can also generate stochastic multifractal time series with explicit expressions of $\tau(q)$ and $f(\alpha)$ [360]. For lognormal $\mathcal{W}$-cascades [360], we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = -\frac{\sigma^2}{2\ln2}q^2-\frac{\mu}{\ln2}q-1 \label{Eq:MF:MathModel:W:cascades:LogNormal:tau:q} \nonumber \end{align} \tag{ 275 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) = 1-\frac{(\alpha+\mu/\ln2)^2\ln2}{2\sigma^2}, \label{Eq:MF:MathModel:W:cascades:LogNormal:f:alpha} \nonumber \end{align} \tag{ 276 }$$

where $\mu$ and $\sigma^2$ are respectively the mean and the variance of $\ln|W|$. The minimum and maximum singularities are obtained when $f(\alpha)=0$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha_{\min,\max} = \pm \frac{\sqrt{2}\sigma}{\sqrt{\ln2}}-\frac{\mu}{\ln2}. \label{Eq:MF:MathModel:W:cascades:LogNormal:alpha:min:max} \nonumber \end{align} \tag{ 277 }$$

For log-Poisson $\mathcal{W}$-cascades [360], we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q) = \frac{\lambda(1-\delta^q)-\gamma{q}}{\ln2}-1 \label{Eq:MF:MathModel:W:cascades:LogPoisson:tau:q} \nonumber \end{align} \tag{ 278 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) &= \left(\frac{\alpha}{\ln\delta}+\frac{\gamma}{\ln2\ln\delta}\right)\left[\ln\left(\frac{\alpha\ln2+\gamma}{-\lambda\ln\delta}\right)-1\right] \nonumber \\ &~~-\frac{\lambda}{\ln2}+1, \label{Eq:MF:MathModel:W:cascades:LogPoisson:f:alpha} \nonumber \end{align} \tag{ 279 }$$

where the law of $\ln|W|$ is the same as $P\ln\delta+\gamma$, and $\lambda$ is the mean and the variance of the Poisson variable P. The differences of the multifractal expressions between the $\mathcal{W}$-cascades and the conventional cascade models stem from the different expressions of the multiplier distributions.

The discussions about stochastic continuous multifractals with some other multiplier distributions can be found in [82, 487, 488], in which anomalous multifractal behavior emerges.

Mandelbrot et al proposed a seminal multifractal model of asset returns (MMAR) [217], which is also termed Brownian motions in multifractal time (BMMT) [489]. In the MMAR model, the logarithmic price is assumed to be a subordinate process of a Brownian motion and a multifractal process of trading time. The MMAR model contains three assumptions:

Assumption 1. $X (t)$ is a subordinate process

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle X (t) \equiv B_H[\theta(t)], \label{Eq:MFModels:MMAR:Compound} \nonumber \end{align} \tag{ 280 }$$

where t is the clock time, $B(t)$ is a fractional Brownian motion with self-affine index H and $\theta(t)$ is a stochastic non-decreasing trading time.

Assumption 2. The trading time (or 'subordinator') $\theta(t)$ is the cumulative distribution function of a multifractal measure $\mu$ defined on $[0, T]$ such that it is a multifractal process with continuous, non-decreasing paths, and stationary increments.

Assumption 3. The processes $B(t)$ and $\theta(t)$ are independent.

The trading time $\theta(t)$ plays a crucial role in the MMAR model, which presumably causes the price $X(t)$ to be multifractal. It is shown that, under these three assumptions, the process $X(t)$ is multifractal with stationary increments and scaling function [217].

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_X(q) = \tau_\theta(Hq). \label{Eq:MFModels:MMAR:tauX:tauTheta} \nonumber \end{align} \tag{ 281 }$$

When $\theta(t)$ is a canonical multifractal, the increments of $X(t)$ have fat tails [217, 490]. When H  =  1/2, the price process displays both uncorrelated increments and long memory in volatility [217]. Therefore, the MMAR model with H  =  1/2 is able to reproduce the main stylized facts of financial returns [94]. A natural choice of the canonical multifractal measure $\mu$ in assumption 2 is from the random multiplicative cascade models [143, 162, 217, 483, 486].

Recent empirical results show that the intertrade durations do have multifractal nature [491–495], which provides strong supports for assumption 2. Calvet et al developed in detail the theory of multiplicative multifractals that is closely relevant to the $\theta(t)$ process [483]. Fisher et al unveiled the presence of multifractal scaling laws in the DEM/USD (bid/ask quotes from 31 October 1992 to 1 September 1993 and daily data from June 1973 to December 1996) and JPY/USD exchange rates data [496]. They found that the simulated multifractal data from the MMAR generating process is consistent with the DEM/USD data, but do not appear completely consistent with the JPY/USD data [496]. Calvet and Fisher confirmed the presence of multifractality and the consistency with the MMAR model in the DEM/USD data, and in the daily data from 1962 to 1998 of the value-weighted NYSE-AMEX-NASDAQ index ('CRSP Index'), Archer Daniels Midland (ADM), General Motors (GM), Lockheed-Martin, Motorola, and United Airlines (UAL) as well [143].

Günay investigated the comparative performance of the lognormal MMAR model with respect to four benchmark econometric models (GARCH, EGARCH, FIGARCH and MRS-GARCH) using daily stock index returns of four emerging markets (Croatia, Greece, Poland and Turkey) from 4 January 2000 to 3 July 2014 [497]. He estimated the parameters of the five models and simulated each model with 1000 runs. He found that the deviation of the average mass scaling exponents of the lognormal MMAR model from those of the empirical results is the smallest.

Using the daily returns of 12 foreign exchange rates against the USD from January 1993 to February 2012 as an example, Goddard and Onali proposed a joint test for long memory and multifractality in the framework of the lognormal MMAR model [498].

Batten et al modified the MMAR model by representing the trading time in assumption 2 with the series of volume ticks [499], that is, the trading time $\theta(t)$ is given by the normalized volume ticks $\overline{V}_t$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{V}_t = V_t/\sum_{t=1}^TV_t, \label{Eq:MFModels:MMAR:Volume:as:Theta} \nonumber \end{align} \tag{ 282 }$$

where $V_t\geqslant 0$ denotes the number of volume ticks at time t. The rational of this model is based on the fact that trading volumes possess multifractal nature that is universal in different markets [222, 500–504]. They studied the 5 min EUR/USD spot quotes and trading ticks from 5 January 2006 to 31 December 2007. They found that the model can produce admissible VaR forecasts at the 12 h forecast horizon and outperform the historical simulation approach and the benchmark GARCH(1,1) location-scale VaR model at the daily horizon.

A caveat about the construction of some versions (including the initial one [217]) of the MMAR model should be stressed. Assumption 2 for the subordinator $\theta(t)$ defines it in general without paying attention to the condition of causality. In general, a multifractal measure is a geometrical construction that exists over all times, without a specific definition or dependence on the arrow of time. In other words, the multifractal measure pre-exists time. The fact that the subordinator $\theta(t)$ can mimic a time flow results from the construction consisting in taking the cumulative distribution function of this multifractal measure. However, the measure at a given instant already 'knows' its future values in the constructions that are for instance purely geometric without attention to causality. While these classes of MMAR models may appear good candidates to describe financial price time series, they should be only considered as toys, with no serious relevance to real markets. Indeed, the arguably most important of all properties that models of financial markets should obey is that of causality. Otherwise, spurious conclusions will derive in portfolio allocation and optimisation, in derivative pricing, investment strategies and so on.

In assumption 2, when the trading time $\theta(t)$ is the c.d.f. of a Poisson multifractal measure $\mu$ defined on $[0, T]$, we obtain the Poisson MMAR model or PMM for short [505]. The PMM improves the MMAR in the sense that the resulting measure is grid-free since its partitioned sub-intervals do not have identical length b^−k and it can be viewed as a stochastic volatility model with a Markov latent vector [505]. The PPM leads to the Markov-switching multifractal models for volatility.

The Markov-switching multifractal (MSM) model designed by Calvet and Fisher is a discrete-time Markov process with multi-frequency stochastic volatility [505, 506]. The return process is specified as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_i = \sigma_i \epsilon_i, \label{Eq:MFModels:MSM:ri:sigmai:epsiloni} \nonumber \end{align} \tag{ 283 }$$

where the random variables $\newcommand{\e}{{\rm e}} \epsilon_i$ are i.i.d. standard Gaussians ${\mathcal{N}}(0,1)$ and $\sigma_i$ is a stochastic volatility with $\bar{k}$ volatility components $M_{1,i}, M_{2,i}, \cdots, M_{\bar{k},i}$ decaying at heterogeneous frequencies $\gamma_1, \gamma_2, \cdots, \gamma_{\bar{k}}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_i^2 = \sigma^2 M_{1,i} M_{2,i} \cdots M_{\bar{k},i} = \sigma^2 \prod_{k=1}^{\bar{k}}M_{k,i} \label{Eq:MFModels:MSM:ri:sigmai2:Mki} \nonumber \end{align} \tag{ 284 }$$

where $\sigma$ is the unconditional standard deviation of the return r_i under the assumption that the multipliers M_1,i, M_2,i, ⋯, $M_{\bar{k},i}$ are statistically independent. The random volatility components M_k,i are random multipliers that are persistent, non-negative and canonical such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\mathbb{E}}(M_{k,i})=1 \label{Eq:MFModels:MSM:E:Mki} \nonumber \end{align} \tag{ 285 }$$

for all i's. It is further assumed for simplicity that the multipliers $M_{1,i}, M_{2,i}, \cdots, M_{\bar{k},i}$ at a given time i are statistically independent.

In the MSM model, one stacks the volatility components in period i into a $1\times\bar{k}$ vector

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\overrightarrow{M}}_i= [M_{1,i}, M_{2,i}, \cdots, M_{\bar{k},i}]. \label{Eq:MFModels:MSM:Mi:vector} \nonumber \end{align} \tag{ 286 }$$

The vector ${\overrightarrow{M}}_i$ is assumed to be a first-order Markov process and thus defines the volatility states. When the state vector ${\overrightarrow{M}}_{i-1}$ is determined in period i  −  1, the i period multiplier M_k,i for each $k\in\{1,2,\cdots,\bar{k}\}$ is either drawn from a fixes distribution $P(M)$ with probability $\gamma_k$ or remains unchanged at its current state:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{k,i} = \left\{\begin{array}{@{}lll@{}} m_{k,i} & {\mathrm{with~probability~}} \gamma_k, \nonumber \\ m_{k,i-1} & {\mathrm{with~probability~}} 1-\gamma_k, \nonumber \\ \end{array} \right. \label{Eq:MFModels:MSM:Mki:Markov} \nonumber \end{align} \tag{ 287 }$$

where m_k,i is drawn from $P(M)$. The MSM model allows for a wide variety of specifications for the distribution $P(M)$ and any discrete or continuous distribution with ${\mathbb{E}}(M)=1$ defined on $M\in(0,\infty)$ could be utilized.

Introducing two parameters $\gamma_{\bar{k}}\in(0,1)$ and $b\in(1,\infty)$, the transition probabilities are specified as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \gamma_k=1-(1-\gamma_{\bar{k}})^{b^{k-\bar{k}}}, \label{Eq:MFModels:MSM:gammak} \nonumber \end{align} \tag{ 288 }$$

or equivalently

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \gamma_k = 1-(1-\gamma_1)^{b^{k-1}}. \label{Eq:MFModels:MSM:gammak:bis} \nonumber \end{align} \tag{ 289 }$$

Obviously, $\gamma_k$ is a monotonically increasing function of k so that M_k,i with small k values correspond to low-frequency components and M_k,i with large k values correspond to high-frequency components. For small values of $\gamma_1$ and k, a Taylor polynomial approximation of equation (289) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \gamma_k \approx \gamma_1 b^{k-1}. \label{Eq:MFModels:MSM:gammak:TaylorApprox} \nonumber \end{align} \tag{ 290 }$$

Hence, the transition probabilities of low-frequency components grow approximately at geometric rate b and the growth rate of the transition probabilities of high-frequency components slows down.

The MSM model can produce the main stylized facts of financial asset returns. It contains multifrequency volatility and the price process is a multifractal jump-diffusion process [507]. The changes of low-frequency multipliers induce discontinuous volatility and long memory in volatility, high-frequency multipliers produce substantial outliers and hence fat tails in returns, and the hierarchical cascade of the multiplicative process ensures the emergence of multifractality in returns. The multifractal spectrum of the price process is fully determined by the probability distribution $P(M)$ of the multipliers [143].

Calvet and Fisher proposed a parsimonious setup in which the probability distribution $P(M)$ is a binomial distribution with the random variables taking values $m_0\in[1,2]$ or $2-m_0\in[0,1]$ with equal probability, which is a binomial MSM model [506]. For the binomial MSM model and other MSM models with discrete multiplier distributions, the model parameters ($\sigma$, b, $\gamma_{\bar{k}}$ and parameters in $P(M)$) can be estimated via the maximum likelihood estimation (MLE) method, because the close-form likelihood can be expressed and the number of frequencies $\bar{k}$ can be determined using the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) by comparing the likelihoods of different MSM($\bar{k}$) models that have different numbers of free parameters [506].

Using the daily exchange rates of the Deutsche Mark (DEM), Japanese Yen (JPY), British Pound (GBP), and Canadian dollar (CAD) against the US Dollar for about 30 years, Calvet and Fisher estimated binomial MSM models with four parameters and more than a thousand states. They found that the MSM model outperforms benchmark econometric models (GARCH, Markov-switching GARCH and FIGARCH) both in-sample and out-of-sample and considerable gains in forecasting accuracy are obtained at horizons of 10 to 50 d [506]. Using daily data of two stock indices (the Dow Jones Composite 65 Average Index and NIKKEI 225 Average Index) from January 1969 to October 2004, two foreign exchange rates (GBP and AUD against USD) from March 1973 to February 2004, and two US treasury bond rates (TB-1y and TB-2y) from June 1976 to October 2004, it is found that the estimated model captures the temporal dependence of the data, through estimating and comparing the generalized Hurst exponents $H(q=1)$ and $H(q=2)$ for both empirical data and simulated data of the binomial MSM model [508].

Lux studied the lognormal MSM model [509], in which multipliers are determined by random draws from a lognormal distribution with parameters $\lambda$ and S:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(M) = \frac {1}{MS {\sqrt {2\pi}}} \exp\left(-{\frac {\left(\ln M+\lambda \right)^{2}}{2S^2}}\right). \label{Eq:MFModels:MSM:LN:P(M)} \nonumber \end{align} \tag{ 291 }$$

According to the canonical condition (285), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S^2 = 2\lambda. \label{Eq:MFModels:MSM:LN:CanonicalCond} \nonumber \end{align} \tag{ 292 }$$

Hence, the lognormal SMS model has four parameters. Lux proposed a generalized method of moments (GMM) estimator together with linear forecasts for model calibration, which can be applied to any continuous distribution with any number of volatility components. His numerical experiments show that the GMM estimator performs well for the binomial and lognormal models.

The mass exponent function $\tau(q)$ defined through structure functions of the lognormal MSM model can be explicitly expressed as [263]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(q)&=\frac{q}{2}-1-\left(\frac{q^2}{4}-\frac{q}{2}\right)(\lambda-1) \nonumber \\ &=-\frac{\lambda-1}{4}q^2+\frac{\lambda}{2}q-1. \label{Eq:MFModels:MSM:LN:tau:q} \nonumber \end{align} \tag{ 293 }$$

The parameter $\lambda$ can be determined by fitting equation (293) to the empirically obtained $\tau(q)$ curve. Based on the DAX data, Eisler and Kertész showed that the MF-DFA method provides a nice estimate of $\lambda$, which is only slightly worse than the GMM estimate [263]. It follows from the Legendre transform that (see [510])

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha(q) = -\frac{\lambda-1}{2}q+\frac{\lambda}{2} \label{Eq:MFModels:MSM:LN:alpha:q} \nonumber \end{align} \tag{ 294 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(\alpha) = 1- \frac{(\lambda-2\alpha)^2}{4(\lambda-1)}. \label{Eq:MFModels:MSM:LN:f:alpha} \nonumber \end{align} \tag{ 295 }$$

We see that $\alpha(q)$ is a linear function of q without limits. When $q\geqslant\lambda/(\lambda-1)$, we have $\alpha\leqslant 0$.

Empirical analysis of five exchange rates (DEM, GBP, CND, JPY and CHF against USD) from 1 January 1979 to 31 December 2004 shows that the lognormal MSM model performs comparably to the binomial MSM model [509]. Using 1 min high frequency data of the Polish WIG20 index sampled over the period 1999–2004, Oświęcimka et al found that the lognormal MSM model is able to reproduce some main stylized facts including the fat-tailed distribution, the absence of correlations and multifractality in the return time series [510]. Based on the same data sets as in [508], it is found that the lognormal MSM model performs almost identically to the binomial MSM model in reproducing the generalized Hurst exponents $H(1)$ and $H(2)$ of logarithmic prices, indicating that the parsimonious discrete specification of the binomial MSM model is flexible enough for empirical analysis [511]. Similar results are obtained for absolute returns and squared returns for small q values, but not for large q values [512]. The determination of the generalized Hurst exponents $H(q)$ in these works is based on the structure function approach. The situation might be improved if one uses other multifractal analyses such as MF-DFA [263].

The time series generated from the lognormal MSM model has a vanishing leverage autocorrelation [263], which is inconsistent with the well-known leverage effect in financial time series [94]. Using the idea in [513], Eisler and Kertész extended the MSM model by introducing leverage autocorrelations [263]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_i = \exp\left(-\sum_{i'<i}{\mathrm{sign}}(\epsilon_{i'})K(i-i')\right)\sigma_i \epsilon_i, \label{Eq:MFModels:MSM:ri:sigmai:epsiloni:Leverage} \nonumber \end{align} \tag{ 296 }$$

where $K(i)$ is a kernel function, which is qualitatively the leverage autocorrelation function in an appropriate limit [263, 513]. Obviously, this extension can be applied to other stochastic volatility models. Alternatively, incorporating the MSM volatility into the mean-reversion term in the square-root volatility process of the Heston model [514], the leverage effect emerges naturally due to the correlation between the Gaussian innovations of the volatility and the Gaussian innovations of the stock price in the Heston model [515]. Note that a large class of stochastic volatility processes can produce the leverage effect [262, 516].

The univariate MSM model can be extended to bivariate and multivariate MSM models to investigate volatility comovement between two or many financial assets [517]. For multivariate MSM models with discrete multiplier distributions, it is possible to derive the closed-form likelihoods, and parameter estimation can be conducted by maximum likelihood for state spaces of moderate size and by simulated likelihood through a particle filter in high-dimensional cases [517]. The GMM estimation approach proposed in [509] can also be designed for multivariate MSM models with discrete or continuous multiplier distributions [518]. This model is homogeneous, which can be improved by allowing correlations between volatility components to form a non-homogeneous bivariate MSM model to gain better risk prediction power [519].

Other extensions include the Copula-MSM model for characterizing spot and weekly futures price dynamics and more generally other pairs of financial assets [520], the MSM-t model in which the innovations $\newcommand{\e}{{\rm e}} \epsilon_i$ follow the Student distribution [521], the MSM-skewed t model in which the innovations follow a skewed Student distribution [522], the level-MSM model for interest rates which incorporates the well-known level effect observed in interest rates [523], the modified Hull–White model of interest rates in which the volatility of the short term rate is driven by an MSM model [524], the cyclic MSM model for insurance claims arising from climatic events in which the term $\newcommand{\e}{{\rm e}} \epsilon_i$ is set to the smoothed average of claims observed during the corresponding month of the year [525], the realized volatility lognormal MSM (RV-LMSM) model in which the non-RV $\sigma_i$ in equation (284) is replaced by the realized volatility [526], and the MSM duration (MSMD) model for the analysis of intertrade durations in financial markets [527–529].

Overall, different specifications of MSM-type models have comparative performance and on average outperform benchmark econometric models in capturing the multifractal behavior of financial variables, in forecasting financial volatility, and in measuring market risks.

Bacry et al introduced the multifractal random walk (MRW) to model stochastic volatility in which the volatility can be reduced as an exponential of a long memory process [445]. The MRW model was inspired by the seminal idea of causal cascades in financial markets, in which a multivariate multifractal description of a portfolio of assets was also introduced and studied [478, 479]. The MRW model possesses solvable multifractal properties [434, 445, 530], which belongs to the class of lognormal continuous cascade models of asset returns [531–533]. The MRW can be extended to the multivariate MRW [434]. Bacry and Muzy further proved the equivalence between the two definitions of MRW and MMAR [531].

Consider a multifractal process $X(t)$, whose discrete process are sampled at discrete times $k\Delta{t}$ is $X_{\Delta{t}}(t)=X_{\Delta{t}}(i\Delta{t})$, where $i=1,\cdots,N$. The increment $\delta_{\Delta{t}}X_{\Delta{t}}(i\Delta{t})$ of the MRW is constructed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:dX} \delta_{\Delta{t}}X_{\Delta{t}}(i\Delta{t}) = \varepsilon_{\Delta{t}}(i){{\rm e}}^{\omega_{\Delta{t}}(i)} \nonumber \end{align} \tag{ 297 }$$

and the corresponding MRW is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:X} X_{\Delta{t}}(t)= \sum_{i=1}^N\delta_{\Delta{t}}X_{\Delta{t}}(i\Delta{t}) = \sum_{i=1}^N\varepsilon_{\Delta{t}}(i){{\rm e}}^{\omega_{\Delta{t}}(i)} \nonumber \end{align} \tag{ 298 }$$

where $X(0)=0$ and $N=T/\Delta{t}$. Here, $\varepsilon_{\Delta{t}}$ and $\omega_{\Delta{t}}$ are independent, $\varepsilon_{\Delta{t}}$ is a Gaussian white noise with variance $\sigma^2{\Delta{t}}$, while $\omega_{\Delta{t}}$ is Gaussian and its covariance is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:w:cov} C_{\omega}(\Delta{t},s)={{\rm cov}}({\omega_{\Delta{t}}(i+s)},{\omega_{\Delta{t}}(i)})=\lambda^2\ln\rho_{\Delta{t}}[s] \nonumber \end{align} \tag{ 299 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:rho} \rho_{\Delta{t}}[s]= \left\{\begin{array}{@{}ll@{}} \frac{T}{(s+1)\Delta{t}}, & s+1\leqslant{T}/\Delta{t} \nonumber \\ 1, & {{\rm otherwise}}. \end{array} \right. \nonumber \end{align} \tag{ 300 }$$

Thus the variance is ${{\rm Var}}({\omega_{\Delta{t}}(i)})=\lambda^2\ln(T/\Delta{t})$. To ensure that the variance of $X_{\Delta{t}}(t)$ converges when $\Delta{t}\to0$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:w:mean} {\mathbb{E}}({\omega_{\Delta{t}}(i)})=-{{\rm Var}}({\omega_{\Delta{t}}(i)})=-\lambda^2\ln(T/\Delta{t}). \nonumber \end{align} \tag{ 301 }$$

It follows that the variance of $X_{\Delta{t}}(t)$ is $\sigma^2t$, which is independent of $\Delta{t}$. Because $\omega_{\Delta{t}}$ is Gaussian, this model is the lognormal MRW model and can be extended to other distributions [531, 534].

It has been proved that the scaling exponents of the structure functions of $X_{\Delta{t}}(t)$ is [445]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:zeta:q:1} \zeta(q)=-\frac{\lambda^2}{2}q^2+\left(\frac{1}{2}+\lambda^2\right)q \nonumber \end{align} \tag{ 302 }$$

when the scale s is less than N. When s  >  N, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:zeta:q:2} \zeta(q)=q/2. \nonumber \end{align} \tag{ 303 }$$

Speaking differently, in the latter case, $\omega_{\Delta{t}}$ is a Gaussian white noise so that $\delta_{\Delta{t}}X_{\Delta{t}}(i\Delta{t})$ is also a Gaussian white noise. Hence, $X_{\Delta{t}}(t)$ is a Brownian motion. Ludeña provides approximate confidence intervals of $\zeta(q)$ for empirical data [535]. One can estimate the MRW model using the generalized method of moments [532, 533] or an approximated maximum likelihood method [536].

Bacry et al [445] showed that this construction above has a well-defined and unique time-continuous limit, making the MRW the single quadratic multifractal generalisation of the fractional Brownian process with a bona-fide continuous limit (similar to the continuous limit of a discrete random walk that defines the Wiener process). Sornette et al [537] presented a convenient alternative representation of the MRW in terms of an exponential of a long memory process with a specific asymptotic $1/\sqrt{t}$ decaying kernel, which has been generalised to arbitrary power law kernel [538, 539] and to endogenous self-excited processes [540].

The MRW model is able to capture the multifractal nature of financial returns, such as the daily DJIA returns from 1900 to 2003 [541] and intraday S&P 500 future index over the period of 1988–1999 [532]. It also predicts the power-law tailed distribution of financial returns [532, 542].

If $\varepsilon_{\Delta{t}}$ is a fractional Gaussian noise, the MRW becomes the multifractal fractional random walk (MFRW) and the scaling exponents are [445]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:MRW:zeta:qH} \zeta_{q,H}= \left\{\begin{array}{@{}lll@{}} [qH-q(q-2)\lambda^2]/2, &{{\rm if}}& s\leqslant{T} \nonumber \\ qH, &{{\rm if}}& s>T. \end{array} \right. \nonumber \end{align} \tag{ 304 }$$

Mathematical treatments for the MFRW can be found in [543–546]. The application of the MFRW model is limited in econophysics since financial returns are usually uncorrelated such that $H\approx 0.5$. Rypdal and Løvsletten proposed two mean-reverting MRW models in which the anti-correlations are modeled either through a drift term as in an Ornstein–Uhlenbeck process or by using fractional Gaussian noise with H  <  0.5 for the innovations [547]. They applied the models to estimate the multifractal behavior in the 1 h electricity spot prices in Norwegian Kroner market from 4 May 1992 to 27 August 2011.

The MRW is a useful model in finance which is validated with respect to certain properties on the memory of volatility but is not validated for a fully faithful description of the stock market returns [548].

To capture the leverage effect that past price returns and future volatilities are negatively correlated, Pochart and Bouchaud proposed the skewed multifractal random walk (SMRW) [513]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:SMRW:dX} \delta_{\Delta{t}}X_{\Delta{t}}(i\Delta{t}) = \varepsilon_{\Delta{t}}(i)\exp\left[\omega_{\Delta{t}}(i)-\sum_{j<i}K(\,j,i)\varepsilon_{\Delta{t}}(\,j)\right] \nonumber \end{align} \tag{ 305 }$$

where the kernel $K(i,j)$ decays as a power law

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:SMRW:Kij} K(\,j,i) = K_0(i-j)^{-\beta}\Delta{t}^{-\beta'}~~~~~(\,j<i). \nonumber \end{align} \tag{ 306 }$$

They derived that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Eq:MF:SMRW:zeta} \left\{\begin{array}{@{}lll@{}} \zeta(2q) = q-2q(q-1)\lambda^2 \nonumber \\ \zeta(2q+1) = q+1-2q^2\lambda^2-\beta. \end{array} \right. \nonumber \end{align} \tag{ 307 }$$

The SMRW model has been applied to understanding the volatility smiles in option pricing [513]. Further mathematical treatments and empirical validation using the daily data of five stock market indices (CAC 40, DAX, FTSE 100, S&P 500 and DJIA) from 3 December 1990 to 15 February 2010 can be found in [549].

Sornette et al provided an alternative representation of the MRW [537],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta_{\Delta{t}}X(t) =\int_{t-\Delta{t}}^t{{{\rm d}}}W_1(t'){{\rm e}}^{\omega_{\Delta{t}}(t')} \label{Eq:MF:MRW:Sornette:dX} \nonumber \end{align} \tag{ 308 }$$

where $\delta_{\Delta{t}}X(t)=X(t)-X(t-\Delta{t})$ is the increment, W₁(t) is a Wiener process with unit diffusion coefficient, and $\omega_{\Delta{t}}(t)$ is the log-volatility which obeys an autoregressive equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_{\Delta{t}}(t) = \omega^0_{\Delta{t}}+\int_{-\infty}^t \eta(t)~h_{\Delta{t}}(t-t') {{{\rm d}}}t', \label{Eq:MF:MRW:Sornette} \nonumber \end{align} \tag{ 309 }$$

where $\omega^0_{\Delta{t}}$ is a constant, $\newcommand{\e}{{\rm e}} \eta(t)$ denotes the process of information flow in the form of a standardized Gaussian white noise and the memory kernel $h_{\Delta{t}}(\cdot)$ is a causal function ensuring that the system is not anticipative. Thus $\omega(t)$ represents the response of the market to incoming information $\newcommand{\e}{{\rm e}} \eta(t)$ up to the time t. At time t, the distribution of $\omega_{\Delta{t}}(t)$ is Gaussian with mean $\omega^0_{\Delta{t}}$ and variance $V_{\Delta{t}} = \int_0 ^\infty{h}^2_{\Delta{t}}(t'){{{\rm d}}}t' = \lambda^2 \ln \left( \frac{Te^{3/2}}{\Delta t} \right)$. The covariance of $\omega_{\Delta{t}}(t)$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_{\Delta{t}}(s\Delta{t}) = \int_0 ^\infty h_{\Delta{t}}(t) h_{\Delta{t}}(t+s|\Delta{t}|) {{{\rm d}}}t, \label{Eq:MF:MRW:Sornette:CoVar} \nonumber \end{align} \tag{ 310 }$$

specifies entirely the random process. When $\Delta{t}$ is sufficiently small, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_{\Delta{t}}(t') \sim h_0 \sqrt{\frac{\lambda^2 T}{t'}} ~~~ {{\rm for}}~ {\Delta{t}} \ll t' \ll T. \label{Eq:MF:MRW:Sornette:h} \nonumber \end{align} \tag{ 311 }$$

The long-range dependence and multifractality of the stochastic volatility are caused by the slow power-law decay of the memory kernel in equation (311).

Different from the techniques and evidence of multifractality in the domain of $f(\alpha)$ reported in this review, it can be derived from the above representation of the MRW that there is also a nice temporal signature of multifractality, in the sense that there is a full spectrum of exponent $p(A)$ describing the 1/t^p (A) relaxation of the volatility from a peak of amplitude A [537]. Such multifractality has also been found in the multifractal stress activation model of earthquakes [550, 551], which has been confirmed on many catalogs [552, 553], and exogenous shocks and on endogenous multifractal responses in financial markets [58, 59].

Saichev and Sornette generalized the MRW model to a large class of continuous stochastic processes [538]. Similar to equation (308), the increment process is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta_{\Delta{t}}X(t) =\int_{t-\Delta{t}}^t \kappa{{\rm e}}^{\omega(t')} {{{\rm d}}}t', \label{Eq:MF:QuasiMF:dX} \nonumber \end{align} \tag{ 312 }$$

where $\omega(t)$ is expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega(t) = \int_{-\infty}^t h(t-t'){{{\rm d}}}W(t') \label{Eq:MF:QuasiMF:w} \nonumber \end{align} \tag{ 313 }$$

through a Wiener process $W(t)$ and a power-law kernel

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h(t) \sim h_0 \left(1+\frac{t}{\ell}\right)^{-0.5-\varphi}{\mathcal{H}}(t), \label{Eq:MF:QuasiMF:h} \nonumber \end{align} \tag{ 314 }$$

where $\varphi\geqslant0$ and the microscopic characteristic scale $\newcommand{\e}{{\rm e}} \ell$ regularizes the singularity of the power law in the propagator $h(t)$ at t  =  0 and ${\mathcal{H}}(t)$ is the Heaviside function ensuring the condition of causality. The special case $\varphi=0$ recovers the MRW model with a truncation at an integral scale. The spectrum of the quasi-multifractal process can be computed and it is found that multifractality holds over a large range of dimensionless scales for $\varphi>0$ [538, 554].

Saichev and Filimonov considered the discrete version of the quasi-multifractal model and proposed a new method for the determination of the quasi-multifractal spectrum [539].

Filimonov and Sornette introduced the self-excited multifractal (SEMF) model, in which the absolute increments of the process are expressed as exponentials of a long memory of past increments [540]. The SEMF model has a self-excitation mechanism combined with exponential nonlinearity such that future values of the process explicitly depends on past ones. The self-excitation captures the microscopic origin of the emergent endogenous self-organization properties of financial markets, such as the cascade in information flows, the triggering of aftershocks by previous large fluctuations and the 'reflexive' interactions.

Without loss of generality, we denote the time increment $\Delta t \to 1$. In discrete time $t=0, 1, 2, \cdots, n$, the SEMF process reads $X_n=\sum_{i=0}^n \Delta{X}_i$, where its increments are given by the following relation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta{X}_n = \sigma \eta_n\exp\left\{-\frac{\omega_n}{\sigma}\right\}, \label{Eq:MF:SEMF:dX} \nonumber \end{align} \tag{ 315 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_n=\sum_{i=0}^{n-1}\Delta{X}_ih_{n-i-1} \label{Eq:MF:SEMF:w} \nonumber \end{align} \tag{ 316 }$$

in which the kernel has a power-law form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_n=h_0 n^{-0.5-\varphi}. \label{Eq:MF:SEMF:h} \nonumber \end{align} \tag{ 317 }$$

The random variables $\newcommand{\e}{{\rm e}} \eta_n$ is an i.i.d. Gaussian with zero-mean and unit variance representing a flow of external news, the parameter $\sigma$ stands for the impact amplitude of the external news and sets the dimension and scale of $\Delta{X}_n$ and X_n. The sum in equation (316) ensures that the dynamics is self-excited.

The discrete SEMF process can be viewed as the simplest multifractal generalization of the GARCH process [133], and in fact it is strongly related to EGARCH and other variants [131, 137] defined in section 1.3. The external flow of news $\newcommand{\e}{{\rm e}} \eta_n$ can be either positive or negative, controlling the signs of the returns. The absolute returns (or volatilities) are determined by all past returns with a decaying weight as a function of their distance to the present. The SEMF process can reproduce all the standard stylized facts found in financial time series, such as heavy tails of the return distribution, absence of autocorrelation of the returns, long-range correlation of the squared increments, multifractality and the leverage effect, which are robust to the specification of the parameters and the shape of the memory kernel [540].

The continuous-time SEMF process can be formally defined by the stochastic integro-differential equation in the Ito sense [540]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{{\rm d}}}X(t)=\sigma\exp\left\{-\frac1\sigma\int_{-\infty}^t h(t-t'){{{\rm d}}}X(t') \right\}{\rm d}W(t), \label{Eq:MF:SEMF:dX:continous} \nonumber \end{align} \tag{ 318 }$$

where $dW(t)$ is the increment of the regular Wiener process and $h(t)$ is a memory kernel function. The expression (318) is justified as the formal notation of the discrete formulation (315) in the limit when $\Delta{t}\to0$. When there is no memory in the kernel such that $h(t)=0$, we retrieve the standard random walk.

Following Scher and Montroll's valley model [556], the unconditional interevent duration density $\psi(t)$ is expressed based on the conditional density of the priority $\newcommand{\e}{{\rm e}} \epsilon$ (or volume for transactions):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi(t) = \int_{-\infty}^{\infty} \psi(t|\epsilon)\rho(\epsilon){{{\rm d}}}\epsilon, \nonumber \end{align} \tag{ 319 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi(t|\epsilon)=s_0^{-1} {\rm e}^{-\beta\epsilon}\exp\left(-s_0^{-1} {\rm e}^{-\beta\epsilon} t \right) \nonumber \end{align} \tag{ 320 }$$

is the conditional density with $s_0, \beta>0$ and $\newcommand{\e}{{\rm e}} \rho(\epsilon)$ is the density of $\newcommand{\e}{{\rm e}} \epsilon$ [561, 562]. This framework is in essence a mixture of distributions hypothesis from subordinated processes [564].

When the priority $\newcommand{\e}{{\rm e}} \epsilon$ has different density $\newcommand{\e}{{\rm e}} \rho(\epsilon)$, $\psi(t)$ will have different forms [561]. When $\newcommand{\e}{{\rm e}} \epsilon$ is a constant, $\psi(t)$ is an exponential. When $\newcommand{\e}{{\rm e}} \epsilon$ is distributed uniformly, $\psi(t)\sim t^{-1}$. When $\newcommand{\e}{{\rm e}} \epsilon$ follows the Laplace distribution, $\psi(t)$ has a power-law tail with arbitrary exponent for $t\gg s_0$.

Within the MF-CTRW formalism [561], the q-moment is defined as follows

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle{t^q}\rangle = \int_0^{\infty} t^q\psi(t){{{\rm d}}}t = \int_0^{\infty} \langle{t^q|\epsilon}\rangle\rho(\epsilon){{{\rm d}}}\epsilon, \nonumber \end{align} \tag{ 321 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle{t^q|\epsilon}\rangle = \int_0^{\infty} t^q\psi(t|\epsilon){{{\rm d}}}t=\Gamma(1+q)s_0^q{\rm e}^{q\beta\epsilon}. \nonumber \end{align} \tag{ 322 }$$

It is found that none of the three $\newcommand{\e}{{\rm e}} \rho(\epsilon)$ densities discussed above could produce a multifractal structure [561]. Rather, a proper selection of weight $\newcommand{\e}{{\rm e}} \rho(\epsilon)$ to produce multifractality is in the form of the exponential or Weibull family:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho(\epsilon) \sim \exp\left[-\left(\frac{|\epsilon-\mu|}{\sigma}\right)^\alpha\right], \nonumber \end{align} \tag{ 323 }$$

where $\alpha>0$ and $\sigma>0$.

The moment $\langle{t^q}\rangle$ exists if $\alpha>1$ and q  >  −1. When $\alpha=2$, one has

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle{t^q}\rangle \sim L^{q^2}, \nonumber \end{align} \tag{ 324 }$$

where $\newcommand{\e}{{\rm e}} L=\exp(\alpha^2\beta^2/4)$. When $\alpha\neq2$, one has

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle{t^q}\rangle \sim L^{|q|^{\alpha/(\alpha-1)}}, \nonumber \end{align} \tag{ 325 }$$

where $\newcommand{\e}{{\rm e}} L=\exp[(\alpha-1)(\beta\sigma/\alpha){}^{\alpha/(\alpha-1)}]$. Kasprzak et al further derived the generalized partition function in a closed form within the MF-CTRW model and unveiled higher-order phase transitions in financial markets [562].

Percolation models for financial markets are based on stochastic formation of opinion clusters with many variants of the rules of cluster merging and splitting [565, 566]. The multifractal analysis of returns using the structure function approach illustrates that the scaling exponent function $\xi(q)$ varies with model parameters and that properly determined parameters are able to produce results comparable to those of the DJIA index [567]. Several percolation models have been studied aiming at exploring the emergence of multifractal behavior in artificial financial time series, such as the stochastic cellular automata model [568], the particle-cluster aggregation model on scale-free networks [569, 570], the percolation model on the Sierpinski carpet lattice [571], and other variants of the percolation model [572–577].

Spin models, especially Ising models, have also been adopted to model the price formation process of financial assets [578]. In a sophisticated extension of the Ising model, Sornette and Zhou incorporated global news for all agents, imitation among neighboring agents and idiosyncratic preference of individual agents with an irrational interpretation of the predictive power of global news [58, 59]. The model can generate multifractal financial time series in a self-organized way. There are also Ising models for stock markets with market makers [579–581], on Sierpinski carpet lattice [582], and on small-world networks [583, 584]. These models also succeeded in generating multifractal returns.

Mike and Farmer introduced an empirical behavioral model for order-driven markets [585], known as the Mike-Farmer model [586], in which the basic rules (the long memory in order directions and a Student distribution in relative order prices) are based on the statistical regularities extracted from the order flow data of real stocks. Gu and Zhou improved the Mike-Farmer model by considering further the long memory in relative order prices [587]. The Gu-Zhou model is able to produce multifractality in the generated return series [587] and in the recurrence intervals of returns [588].

Thompson and Wilson found that the zero-intelligence model of Farmer et al [589] is not able to produce multifractality in the generated returns, but the positive-intelligence model can [590]. Crepaldi et al showed that, although the usual minority games do not produce multifractal returns, the non-synchronous minority games in the non-ergodic phase can generate multifractal returns [591]. Multifractal behavior of financial time series is also observed in other types of models, such as an artificial dealer market model [592], a double auction model [593], the stochastic voter model [594, 595], the stochastic interacting contact model [596, 597], the epidemic model [598], and strategic models with fundamentalist agents [599, 600].

In most multifractal analyses, we need to calculate high-order moments. The problem of computational precision may arise. Consider the partition function approach as an example. When $m(s,t)\ll1$ and $q\gg1$, the corresponding values of the partition function $\chi$ will be too small such that the computer is out of memory. To overcome this problem, we can compute the logarithm of the partition function ln$\chi_{q}(s)$ rather than the partition function itself. With a simple manipulation, we obtain [205]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln\chi_{q}(s) &= \ln\sum_{t}\left[\frac{m(s,t)}{m_{\max}}\times{m_{\max}}\right]^q \nonumber \\ &= \ln\sum_{t}\left[\frac{m(s,t)}{m_{\max}}\right]^q +q\ln{m_{\max}}, \nonumber \end{align} \tag{ 326 }$$

where $m_{\max}=\max \limits_{t}\{m(s,t)\}$ is the maximum of $m(s,t)$ in the boxes.

Multifractal analysis deals with the scaling behavior of time series at different scales. The most common way is to partition the time series of length N into nonoverlapping windows of scale s. Hence, it is not uncommon that N/s is not an integer. If one uses a window with size less than s, one may obtain incorrect results such as $f(\alpha)>1$, which happened sometimes in empirical studies. A usual strategy is to cover the time series from both sides and keep only the $\mathrm{int}[N/s]$ boxes of size s for each covering [178].

Alternatively, we can adopt the idea of the sand box method [601, 602], in which, instead of covering the times series, one uses M_s boxes of size s with the centers of the boxes randomly distributed on the time series [603]. Taking the family of MF-DFA as an example, the qth overall fluctuation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_q(s) = \left\{\frac{1}{M_s}\sum_{v=1}^{M_s} \left[F_v(s)\right]^q\right\}^{\frac{1}{q}}. \nonumber \end{align} \tag{ 327 }$$

This approach has an additional advantage to reduce fluctuations around the power-law scaling, especially when there are intrinsic log-periodic oscillations such as in binomial measures [312].

A preliminary condition for multifractal analysis is the estimation accuracy of the moments. When $|q|$ is large, the estimated moment is usually unreliable since the empirical integrand $x^q\Pr(x)$ for the moment diverges, where x is the measure and $\Pr(x)$ is its probability distribution in the partition function approach [350]. The divergence phenomenon has been empirically shown in turbulence data [268, 350] and financial data [349]. Similarly, when $|q|$ is larger, the estimated derivatives ${\rm d}F_q(s)/{\rm d}s$ deviate significantly from the theoretical values at small scales and fluctuate a lot at large scales in MF-DMA [604]. We believe that the divergence phenomenon is ubiquitous in multifractal analysis.

The estimation accuracy of moments is worse when $|q|$ is larger and the time series length N is shorter. Unfortunately, there are no theoretical results of the dependence on the length N of the maximum q for which the convergence of the integrand is ensured. With turbulence data, Zhou et al estimated that $N(q=8) \approx 6 \cdot 10^{6}$ and $N(q=9) \approx 5 \cdot 10^{7}$ [268]. For low-frequency financial time series, such as daily returns, we suggest that the range of q should be $|q|\leqslant6$, or even smaller.

If a time series is fractal or multifractal, there are power-law scaling relationships between moment-like quantities and scales. In practice, taking $F_q(s) \sim {s}^{h(q)}$ in equation (116) as an example, the scaling exponents are usually estimated using linear least-squares regressions of $\ln{F_q(s)}$ against $\ln{s}$ in a scaling range $[s_{\min},s_{\max}]$. In other words, the scaling laws hold only on the scaling range. The scaling range should be wide enough to ensure the presence of scale-invariant characteristics [605–608]. In many cases, the change of scaling range causes dramatic variations of the estimates and sometimes very strange singularity spectrum shapes [609, 610], some of which are called 'inverted multifractality' [329]. Hence, as the scaling range plays a crucial role in the estimation of scaling exponents, many efforts have been made for its determination. Nevertheless, it remains a difficult task that has not been fully solved [327]. The mostly adopted method is actually eye balling.

Wang et al argued that scaling laws hold in the scaling range but not in the two ranges outside the scaling range such that one can fit the field with three power laws by minimizing the RMS of residuals, called the three-segment method by some researchers [611]. Because the assumption of two extra power laws is merely empirical, this method is not widely adopted. However, when F_q(s) exhibits a crossover phenomenon as in many financial time series [135, 205, 496, 612–615], the idea can be used to simultaneously determine the crossover scale $s_{\times}$ and the two scaling exponents H₁ and H₂ [610, 616]. We rewrite the bi-power law for q  =  2 as follows