Joint multifractal analysis based on the partition function approach: Analytical analysis, numerical simulation and empirical application

Many complex systems generate multifractal time series which are long-range cross-correlated. Numerous methods have been proposed to characterize the multifractal nature of these long-range cross correlations. However, several important issues about these methods are not well understood and most methods consider only one moment order. We study the joint multifractal analysis based on partition function with two moment orders, which was initially invented to investigate fluid fields, and derive analytically several important properties. We apply the method numerically to binomial measures with multifractal cross correlations and bivariate fractional Brownian motions without multifractal cross correlations. For binomial multifractal measures, the explicit expressions of mass function, singularity strength and multifractal spectrum of the cross correlations are derived, which agree excellently with the numerical results. We also apply the method to stock market indexes and unveil intriguing multifractality in the cross correlations of index volatilities.


Introduction
Measurements of a complex evolving system from different view angles provide us many time series that are usually long-range cross-correlated and exhibit multifractal nature. In turbulent flows, there are velocity field, temperature field and concentration field embedded in the same spatial domain. One can measure these quantities at fixed locations to obtain time series, which are mutually correlated [1,2]. In financial markets, there are also many pairs that are cross-correlated, such as market index volatilities, price returns of different markets, price returns of different equities, different quantities of a same equity [3][4][5][6][7][8][9][10]. Moreover, examples come from very diverse fields, including agronomy [11,12], seismic data [13], meteorology [14][15][16], medical science [17,18], geophysics [19], transportation [20][21][22], to list a few.
The joint multifractal analysis is a classic method and has been applied to study the joint multifractal nature between different pairs of time series recorded in natural and social sciences [3,11,12,14,15,17,26]. Due to its elegant geometric nature, many important properties can be derived, which is however very difficult in the frameworks of other methods mentioned above. For instance, although there is numerical evidence and analytical results for the relationship between the cross-multifractal spectrum f xy (α x , α y ) and the multifractal spectra f x (α x ) and f y (α y ) of individual time series [26,28,35,41,63], the problem is not solved. Moreover, the original MF-X-PF method is important because it handles moments with two different orders, while recent methods for multifractal cross-correlation analysis focus only on one order.
In this work, we recover the uni-order MF-X-PF method [26] and propose a direct determination approach for the multifractal spectrum using the idea from the bi-order MF-X-PF framework [2]. Based on this framework, we are able to derive important geometric properties of the uni-order MF-X-PF method. We perform numerical simulations using different mathematical models and explain the results of multifractal binomial measures analytically. Finally, we apply the bi-order MF-X-PF method to stock market indices.

Joint multifractal analysis based on partition function approach
In this section, we first present the joint multifractal analysis based on partition function approach with two moment orders [2], abbreviated MF-X-PF(p, q), and then derive the uni-order method MF-X-PF(q) that was independently proposed recently [26]. Although the joint partition function χ xy (q, s) of the uni-order method can be directly recovered from the joint partition function χ xy (p, q, s) of the bi-order method by posing p = q, we will show that the nexus between the multifractal properties of the two methods is not obvious, which is caused by the application of the steepest descent approach.

MF-X-PF(p, q)
Based on the box-counting idea, the geometric support is partitioned into boxes of size s. We consider two integrated measures m x (s, t) and m y (s, t) in the t-th box. The local singularity strengths α x and α y are defined according to the following relationships: and m y (s, t) ∼ s αy .
Let N s (α x , α y ) denote the number of boxes of size s needed to cover the set of points in which the singularity strengths are around α x and α y with bands dα x and dα y . Hence, the fractal dimension of the set is determined according to [64] N s (α x , α y ) ∼ s −fxy(αx,αy) , in which f xy (α x , α y ) is the joint distribution of the two singularity strengths [2] or the joint multifractal spectrum. We consider the joint partition function This definition is slightly different from that in Ref. [2], in which the orders are p and q rather than p/2 and q/2. In this setting, we recover the traditional partition function when m x = m y and p = q [27]. The joint mass exponent function τ xy (p, q) can be obtained from the following relation χ xy (p, q, s) ∼ s τxy(p,q) .
In practice, for a given pair (p, q), we compute χ xy (p, q, s) for a various of box sizes s and perform linear regression of ln χ xy (p, q, s) against ln s in a proper scaling range to obtain τ xy (p, q).
We insert the two relations in Eq. (1) into the joint partition function, rewrite the sum into a double integral over α x and α y , and then apply the steepest descent approach to estimate the integral at small s values, which leads to where ∂f xy (α x , α y )/∂α x = p/2, and ∂f xy (α x , α y )/∂α y = q/2.
Similar derivation can be done over q and one can obtain the double Legendre transforms Therefore, after obtaining τ xy (p, q), we can numerically determine α x using Eq. (8a), α y using Eq. (8b), and f xy (α x , α y ) using Eq. (8c). From the canonical perspective, one can obtain the f xy (α x , α y ) function directly [2,65,66]. Defining the canonical measures as follows the two singularity strengths α x (p, q) and α y (p, q) and the joint multifractal spectrum f xy (p, q) can be computed by linear regressions in log-log scales using the following equations: The joint mass exponent function can be obtained by using Eq. (5).
Taking derivative of Eq. (12) over q and using Eq. (13), we have Defining that Eq. (14) and Eq. (12) can be rewritten as follows where f xy (α xy ) f xy (α x , α y ). We notice that Eq. (16) has the same form of the Legendre transform [27].
Combining Eqs. (17) and (18) and using τ x (0) = τ y (0) = −1, we have C = 0 and thus Inserting Eq. (15) and Eq. (19) into Eq. (16b), we obtain that We note that Eq. (19) and Eq. (20) still hold when D 0 = 1. In this case, we use τ xy (0) = τ x (0) = τ y (0) = −D 0 to conduct the derivation. These relations were observed numerically using the MF-X-DFA method [28], the MF-X-DMA method [41] and the MF-X-PF(q) method [26]. As shown in Eq. (11), the problem is to handle a measure [m x (s, t)m y (s, t)] 1/2 . From the canonical perspective, we can obtain the f xy (α xy ) function directly [2,65,66]. We can define the canonical measures The two singularity strengths α x (p) and α 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: where Eq. (10a) and Eq. (10b) are used in the second equality, and The joint mass exponent function can be obtained by using Eq. (16b).

Numerical analysis applying MF-X-PF(p, q)
We perform joint multifractal analysis numerically of two binomial measures [67]. We use p x = 0.3 and p y = 0.4 and generate two binomial measures of length 2 20 . Figure 1(a) shows on log-log scales the dependence of χ xy (p, q, s) against box size s for different q with fixed p = 2. It is obvious that the curves for different q exhibit excellent power law relationships. The power-law exponents obtained by linear regressions of ln χ xy (p, q, s) against ln s are estimates of the mass exponents τ xy (p, q), whose contour plot is shown in Fig. 1(e). We find that τ xy (p, q) increases with p and q. Adopting the double Legendre transform in Eq. (8), we obtain numerically the singularity functions α x (p, q) and α y (p, q) and the multifractal spectrum f xy (p, q), whose contour plots are illustrated in Fig. 1(fh) respectively. We find that α x (p, q) and α y (p, q) are decreasing functions of p and q, while f xy (p, q) has a saddle shape. An intriguing feature is that the contour lines are parallel to each other for α x (p, q), α y (p, q) and f xy (p, q). Figure 1(i) plots the singularity spectrum f xy (α x , α y ), which is not a surface but a curve. We also calculate the multifractal functions using the direct determination approach presented in Eq. (10) for comparison. In Fig. 1 and t µ xy (2, q, s, t) ln[µ xy (2, q, s, t)] against ln s for different q with fixed p = 2. The singularity strength functions α x (p, q) and α y (p, q) and the multifractal spectrum f xy (p, q) are computed from the slopes of the lines in these three plots. The corresponding contour plots are presented in Fig. 1(j) to Fig. 1(l), which are the same as those in Fig. 1(f) to Fig. 1(h). The numerical results presented in to Fig. 1 can be derived analytically.

Analytical results for MF-X-PF(p, q)
Let us start with two multifractal binomial measures of length 2 L . Consider two integrated measures m x (s, t) and m y (s, t) in boxes of size s = 2 l . There are n types For the t-the box, we have where k ∈ {1, ..., n}. It follows that Inserting Eq. (25) into Eq. (24a), we get and Note that β and γ depend only on p x and p y . When p x + p y = 1, we have β = −1.
When p x = p y , we have β = 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 β > 0; Otherwise, when (p x − 0.5)(p y − 0.5) < 0, we have β < 0. Combining Eq. (3) and Eq. (26), we obtain where Because m y is a multifractal measure, we have where τ y (Q) has an analytical expression [27]: The joint partition function can be rewritten as follows Comparing Eq.(4) and Eq. (33), we obtain the joint mass exponent function: It follows that and We obtain immediately the relationship between α x and α y α x = γ ln 2 + βα y .
This relationship explains the observation in Fig. 1(i) that f xy (α x , α y ) is a curve along this line rather than a surface and the line segment (37) is the projection of f xy (α x , α y ) onto the (α x , α y ) plane.
We now derive the main geometric properties of α x (Q) and α y (Q). We find that α y (Q) is a monotonically decreasing function of Q, because We can prove that the limits of α y exist when Q → ±∞. We rewrite Eq. (36) as follows We can obtain that Therefore, the solution of Eq. (36) exists and is unique if and only if α y ∈ [α y,min , α y,min ]. The explicit form of the solution is Further, the width of the singularity spectrum of α y is These results explain the parallel observation of the contour lines in Fig. 1(g). When p y = 0.5, ∆α y = 0. In this case, the measure is neither multifractal nor monofractal since it is uniformly distributed on the support. According to Eq. (37), we have which suggests that α x is a strictly monotonic function of Q. Moreover, it is easy to show that Therefore, the solution of Eq. (35) exists, which is unique if and only if α x ∈ [α x,min , α x,max ]. Due to the symmetry between the two measures m x and m y , the results for α x are obvious, provided that we know the geometric properties of α y .
We now turn to investigate the geometric properties of the multifractal spectrum f xy (α x , α y ), which has the following form: It is easy to find that where f xy (Q) f xy (α x , α y ; Q). It indicates that f xy (α x , α y ) is symmetric with respect to the line Q = 0, as numerically shown in Fig. 1(h). Furthermore, we obtain lim Q→±∞ f xy (p, q) = 0.
Taking derivative of f xy (α x , α y ) with respect to Q, we have When Q < 0, df xy (Q)/dQ > 0 so that f xy (Q) is a monotonically increasing function of Q. When Q > 0, df xy (Q)/dQ < 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. These properties explain the parallel feature of the contour lines in Fig. 1(h). We note that the numerical results are in excellent agreement with the analytical results for τ xy (p, q) in Eq. (34), α x (p, q) in Eq. (35), α y (p, q) in Eq. (36), and f xy (p, q) in Eq. (45). Combining Eq. (41) and Eq. (45), we find that f xy (α x , α y ) is a univariate function of α y , or of α x by using Eq. (37).

Numerical analysis applying MF-X-PF(q)
We also apply the MF-X-PF(q) method to the same mathematical example. The results are shown in Fig. 2. We find that the three theoretical relationships in Eq. (19), Eq. (15), and Eq. (20) are nicely verified. In addition, we observe again that the results from the classic partition function approach and the direct determination approach agree with each other. We note that this is also the case for other mathematical and empirical examples investigated in this work. Thus we will not show the results obtained from the direct determination approach in the rest of this paper.
These properties are indicators of monofractality. The mathematical model used here is bivariate fractional Brownian motions (BFBMs). The two components x(t) and y(t) of the BFBM are two univariate fractional Brownian motions with Hurst indices H xx and H yy , respectively. The basic properties of multivariate fractional Brownian motions have been comprehensively studied [69][70][71]. Extensive numerical experiments of other MF-DCCA algorithms have been conducted using bivariate fractional Brownian motions [41,57]. The two Hurst indexes H xx and H yy of the two univariate FBMs and their cross-correlation coefficient ρ are input arguments of the simulation algorithm. By using the simulation procedure described in Refs. [70,71], we have generated as an example a realization of BFBM with H xx = 0.1, H yy = 0.5 and ρ = 0.5. The length of the BFBM is 2 16 . The joint multifractal analysis of the BFBM using the MF-X-PF(p, q) algorithm is presented in Fig. 3.
The corresponding power-law dependence of the joint partition function χ xy (p, q, s) with respect to the box size s for different q's and fixed p = 2 is shown in Fig. 3(a). The scaling ranges span over two orders of magnitude. The slopes of the lines give the estimates of τ xy (p, q), where p and q vary from −10 to 10 with a spacing of 0.1. The resulting mass exponents τ xy (p, q) are shown in the contour plot of Fig. 3(b). We observe that τ xy (p, q) increases with p and q, the contour curves are parallel lines, and the parallel lines are evenly spaced. These features suggest that τ xy (p, q) is a linear function of p and q, which is an indicator of monofractality.
In order to further show the performance of the MFXPF algorithm, we calculate the errors between the estimated exponents τ xy (p, q) and the theoretical exponents as ∆τ xy (p, q) = τ xy (p, q) − [p/2 + q/2 − 1]. Fig. 3(c) shows the dependence of ∆τ xy (p, q) with respect to p and q. All the ∆τ xy (p, q) values are less than 0.15, implying that the algorithm gives good estimates.
By adopting the double Legendre transform in Eq. (8) numerically, we get the singularity strength functions α x (p, q) and α y (p, q) and the multifractal spectrum f xy (α x , α y ), whose contour plots are shown in Fig. 3(d,e,f). The singularity strength functions α x (p, q) and α y (p, q) are close to 1, indicating that the functions α x (p, q) and α y (p, q) are independent of the order p and q. Although there is a trend in each function α x (p, q) and α y (p, q), the theoretical functions α x (p, q) = 1 and α y (p, q) = 1 are basically satisfied. Hence, the MF-X-PF algorithm is able to correctly capture the monofractal nature of the BFBMs. Fig. 3(g) plots the singularity spectrum f xy (α x , α y ), which is a surface and the contour lines are closed curves. It is easy to find that the vast majority of the surface is nearly equal to the theoretical function f xy (α x , α y ) = 1. We observe that the errors ∆τ xy (p, q) is equal to the difference between f xy (p, q) and 1, as shown by the Legendre transform.
We point out that the results using the direct determination approach are exactly the same as shown in Fig. 3. We thus summarize that the theoretical analysis is well verified by the numerical results.

Application to stock market indexes
We now apply the MF-X-PF(p, q) algorithm to investigate the long-range power-law cross correlations of the daily volatility time series of the Dow Jones Industrial Average (DJIA) and the National Association of Securities Dealers Automated Quotations (NASDAQ) index. The daily volatility is defined as the absolute value of the logarithmic difference of daily closing prices: where P (t) is the closing price on day t and has been retrieved for the DJIA and NASDAQ indices. The time period of the samples is from 5 February 1971 to 25 January 2011, containing 10084 data points. The daily return time series of the two indexes are shown in Figure S1 (New J. Phys. online).  Fig. 4(a) shows on log-log scales the dependence of the joint partition function χ xy (p, q, s) with respect to the box size s for different q's and fixed p = 2. We observe nice power-law scaling over about 1.5 orders of magnitude. The contour plot of the exponents τ xy (p, q) is shown in Fig. 4(b), where p and q vary from −10 to 10 with a spacing of 0.1. The contour curves are not straight lines and the spacings between neighboring curves are not equidistant. Fig. 4(c) and Fig. 4(d) illustrate respectively the contour plots of the singularity strength functions α x (p, q) and α y (p, q), which are obtained numerically from τ xy (p, q). We observe that the values of the singularity strength range from 0.6 to 1.2, which are well dispersed. In addition, the singularity strength functions are not monotonic with respect to p or q. Fig. 4(e) illustrates the multifractal function f xy (p, q) obtained from the Legendre transform, whose values range from 0 to 1. The maximum f xy (p, q) = 1 is reached at point (p, q) = (0, 0). Within the investigated intervals of p and q, the small f xy (p, q) values concentrated in the region with large values of p and q. In Fig. 4(f), we present the singularity spectrum f xy (α x , α y ). These empirical findings suggest that the cross correlations between daily volatilities of DJIA and NASDAQ possess multifractal nature, which is consistent with previous results using the MF-X-DFA, MF-X-DMA and MF-X-PF(q) methods [26,28,41,72]. To reveal whether the joint multifractality between the daily volatilities of the two indices remains or changes along time, we perform the MF-X-PF(p, q) analysis in moving windows on a decade basis with a step of one year. The results are presented in Figure S2 (New J. Phys. online). We show six plots in Fig. 5. We find that the joint multifractal singularity spectrum f xy (α x , α y ) changes over time. Moreover, the inclusion or exclusion of financial turmoils (high volatile periods) has a significant impact on the shape of f xy (α x , α y ). In the sample period under investigation, there were two infamous market crises, the Black Monday in 1987 and the latest crisis in 2008. During relatively calm periods, the f xy (α x , α y ) contour looks roughly like an American football. However, when one of the crisis is included, the contours are significantly stretched to the southwest. In other words, the singularity strengthes α x and α y have much smaller values during turmoil periods. This is actually not surprising because this feature is well-documented for ordinary multifractals [27].
We repeat the same analysis for two stocks Du Pont (NYSE:DD) and Exxon Mobil (NYSE:XOM) over time period from 05-Jan-1970 to 01-Sep-2015, containing 11522 data points. The daily return time series of the two stocks are shown in Figure S3 (New J. Phys. online). The results are illustrated in Figure S4 (New J. Phys. online). As expected, very similar results are observed. Two more pairs of financial time series are investigated and the results are presented in Figure S5 to Figure S8 of the Supplementary data (New J. Phys. online). One pair is about crude oil commodities, Arab Light to USA and WTI Cushing. The sample period from is 03-Jan-1991 to 18-Dec-2012, containing 5510 data points. Another pair is about Special Drawing Rights (SDRs) per currency unit for the U.K. pound sterling (GBP) and the U.S. dollar (USD) over time period from 05-Jan-1994 to 01-Sep-2015, containing 5452 data points.
Compared with the results of binomial measures and fractional Brownian motions, the multifractal function and the multifractal singularity spectrum exhibit different shapes for different data sets studied. For example, in Fig. 4(f) for the financial market data there is a pronounced asymmetry, and the spectrum exhibits a stretched shape, in sharp contrast to Fig. 3(f) for the artificial BFBM data. These features reflect the irregular nonlinear traits of financial indexes. Roughly, the spectrum contour parallels to the diagonal α x = α y (cf. Eq. (37)), which is due to the fact that the DJIA and NASDAQ indexes comove along time so that the volatilities fulfill Eq. (26) to certain extend. A direct conjecture is that the correlation coefficient ρ(α x , α y ) is greater if the correlation coefficient ρ(R x (t), R y (t)) is greater. This is validated by Fig. S9 in the Supplementary Data (New J. Phys. online).

Conclusions
We have studied the properties of joint multifractal analysis based on partition function with two moment orders, termed MF-X-PF(p, q). The uni-order method MF-X-PF(q) has then been derived. The main properties of these methods have been obtained analytically. For instance, for the MF-X-PF(q) method, we have obtained the relationship between the joint mass exponent function and the individual mass exponent functions, τ xy (q) = [τ x (q) + τ y (q)]/2, which was numerically and empirically observed in the literature.
We applied the MF-X-PF(p, q) method to multifractal binomial measures. The expressions of mass function, singularity strength and multifractal spectrum of the cross correlations have been derived, which agree excellently with the numerical results. We further validated the performance of the method by using bivariate fractional Brownian motions without multifractal cross correlations. When applied to the daily volatility time series of two stock market indexes, intriguing multifractality in the cross correlations is confirmed. The multifractal properties of these examples are found to be the same when we use the conventional determination approach and the direct determination approach.
Multifractal cross-correlation analysis has been applied in many fields, especially in Econophysics. Although there are numerous methods, most of them consider only one moment order. It is natural that bi-order methods such as MF-X-PF(p, q) can be developed for other uni-order methods. We expect that such bi-order methods will unveil new stylized facts in the analysis of financial time series, which can serve to calibrate agent-based models [73]. In addition, the joint multifractal nature extracted from two long-range cross-correlated time series has potential applications. One possibility is to construct a multi-scale cross-correlation measure, analogous to other DCCA coefficients [38,61,62,74,75]. Another possibility is to construct a measure quantifying market efficiency [76][77][78][79]. A related possibility is to quantitatively characterize the degree of market unrest other than the volatility measure [80].