# Joint Multifractal Analysis and Source Testing of River Level Records Based on Multifractal Detrended Cross-Correlation Analysis

Authors: Liang Wu; Manling Wang; Tongzhou Zhao
Journal: Complexity
 (2020)

---

## Abstract

The joint multifractal analysis is usually conducted in two different variables for their cross-correlations but rarely used for two records of one variable collected at two different places. It is important for the detection of change in multifractality in space. Besides, the cross-correlations in two analyzed series make the analysis of sources of joint multifractality difficult. There are few studies on the source of joint multifractality. We focus on the two issues for two level records at pairs of adjacent sites along one river and carry out an extension of our previous work which is about the single multifractality of one record with the same data set. The data set is collected from 10 observation stations of a northern China river and contains about two million high-frequency river level records. Results of joint multifractal analysis via multifractal detrended cross-correlation analysis show that the change in joint multifractality at pairs of adjacent sites caused by weak cross-correlations can be detected by comparing the single generalized Hurst exponent with the joint scaling exponent function and reveal the effects of human activities on joint multifractality. This analysis provides an approach for detecting the change in multifractality. Following the idea of our previous work, two robust hypothesis tests via a set of pairs of surrogate series are proposed for the source testing of joint multifractality. The analysis of the effects of cross-correlations is carried out via a proposed simultaneously half-shifting technique which can both minimize the cross-correlations between original series and make full use of records. Results of source analysis show not only the effects of autocorrelations in series and probability distribution of river levels but also the effects of cross-correlations between series.


---

## Body

## 1. Introduction

The multifractality has been studied in many areas of applied sciences since its concept was first introduced by Mandelbrot [1]. The structure function method [2] is a classical multifractal analysis method which was generalized from box-counting algorithms [3–5] and once widely used for multifractal analysis of time series. Until now, many multifractal analysis methods were developed, such as multifractal detrended fluctuation analysis (MF-DFA) [6] and wavelet transform modulus maxima (WTMM) [7] and its discrete version, namely, wavelet leaders multifractal analysis (WLMF) [8–11]. These methods have been widely applied to many areas of applied sciences instead of the structure function method [12–22]. The MF-DFA is a multifractal generalization of the detrended fluctuation analysis (DFA) that is used to estimate Hurst exponent of monofractal process [17, 23, 24]. Similarly, the monofractal versions of wavelet-based methods can be also used for Hurst exponent [25–31]. Besides, many complex systems usually have the joint multifractal nature which is exhibited by simultaneous records of two variables. There are also many methods for this issue such as multifractal detrended cross-correlation analysis (MF-DCCA) [32], multifractal cross wavelet analysis (MF-X-WT) [33], and joint multifractal analysis based on wavelet leaders (MF-X-WL) [34].The MF-DCCA which was first proposed by Zhou [32] is a generalization of the MF-DFA [6] for two nonstationary time series. As MF-DFA, the MF-DCCA is easy to implement and robust when time series is short. Besides, it can also eliminate some polynomial trends that may exist in nonstationary time series. The MF-DCCA has been a commonly used joint multifractal analysis method in many fields including hydrology [35–38]. Therefore, we conduct the MF-DCCA for the joint multifractal analysis of river level records. The river levels we analyze are collected from the whole 10 observation stations of a tributary of Haihe River in North China and contain about two million records. The sampling frequency is 6 minutes. The multifractal analysis of these data and their source testing has been conducted in [39, 40]. Such many and high-frequency records are rarely seen in previous literatures. Venugopal et al. [41–43] studied the multifractality of high-resolution temporal rainfall data using wavelet methods.This paper follows up our previous work [39, 40] with the same data set and focuses on the joint multifractality and its source testing of river levels at adjacent stations. The joint multifractality is exhibited by two simultaneous records and distinguished from the single multifractality of one record studied in [40]. The joint multifractal analysis is usually conducted in two different variables for their cross-correlations but rarely used for two records of one variable collected at two different places like the data we analyze in this paper. It is very important for the detection of change in multifractality in space. In our study, the joint multifractal analysis of two river level records at adjacent stations is conducted using MF-DCCA. The obtained joint multifractal results are used to analyze the change in multifractality along the river and can detect the change in multifractality affected by human activities. It provides an approach for detecting the change in multifractality. Besides, there are also few studies on the source of joint multifractality. Jiang et al. [33] studied the source of joint multifractality in financial series using MF-X-WT. It is known that the single multifractality of one record has two main sources: the fatness of the probability distribution of original time series and the different correlations in small and large fluctuations. The two sources can be usually distinguished by analyzing the corresponding surrogate series [44–48]. In addition to the two sources that cause single multifractality, the cross-correlations that exist in two analyzed records also affect the joint multifractality. The cross-correlations in two analyzed records make some surrogate series methods ineffective and the analysis of sources of joint multifractality difficult. Following the idea in [40], we propose two hypothesis tests for the source of joint multifractality, which are based on the empirical distributions of scaling exponent function estimated from some pairs of surrogate series. The pairs of surrogate series are generated via shuffling method and rank-ordered remapping technique which are still effective for joint multifractality. For further testing the effect of cross-correlations, we also propose a simultaneously half-shifting technique that can minimize cross-correlations. Testing results show that the joint multifractality in water levels is mainly caused by the different correlations within and between original series and is also related to the probability distribution of river levels. The simultaneously half-shifting technique further shows the effects of cross-correlations. The MF-DCCA and related comparison results not only show the effects of cross-correlations but also can detect the change of multifractality.The paper is organized as follows: in Section2, we introduce the MF-DCCA method and propose the source testing method and the simultaneously half-shifting technique. The data and related preliminaries are described in Section 3. In Section 4, we report the results and give detailed analysis. Our work is concluded in Section 5.


## 2. Methods

### 2.1. Multifractal Detrended Cross-Correlation Analysis

In this section, we introduce the MF-DCCA [32, 49] which is used for joint multifractal analysis of river level fluctuations in this work. The MF-DCCA is a generalization of the MF-DFA [6] for two nonstationary time series. The method is described as follows.LetXt and Yt be the two processes to be analyzed. First, the profile is computed as follows:(1)X′i≔∑k=1iXk−X¯,Y′i≔∑k=1iYk−Y¯,i=1,…,L,where L is the minimum length of Xt and Yt and X¯ and Y¯ are the means of Xt and Yt, respectively.Dividing the profilesX′t and Y′t into Ns:=intL/s nonoverlapping segments with equal length s and then computing the detrended covariance for each segment, we get(2)FXYv,s=1s∑i=1sX′v−1s+i−xviY′v−1s+i−yvi,v=1,…,Ns,where xvi and yvi are, respectively, the best polynomial fit of X′v−1s+i and Y′v−1s+i in each of the Ns segments v. In this step, the trends are removed via the polynomial fit of profile. The m-order polynomial used in the fitting procedure can eliminate the polynomial trend of order m−1 in the original series. This ensures the fluctuation analysis available for data affected by trends or other nonstationarities. Note that the detrended covariance computed by equation (2) uses absolute values of fluctuations for fluctuation strength, which is different from that of original paper [32] using absolute values of FXYv,s for fluctuation itself. Besides, the operation of absolute value here can also avoid no obvious power-law scaling when q<0 [32, 37], which exists in this study.Average the fluctuation functionFXYv,s over all the segments, given by(3)FXYq,s≔1Ns∑v=1NsFXYv,sq/21/q,where q can be any real values except 0. For q=0,(4)FXY0,s≔exp12Ns∑v=1NslnFXYv,s.The final step is calculatinghXYq via the slope of the log-log plot of FXYq,s versus s, which is based on the power-law [32]:(5)FXYq,s∼shXYq,for larges.The functionhXYq is called the scaling exponent function [32]. When q=2, hXYq is the bivariate Hurst exponent [50].WhenX=Y, MF-DCCA degenerates to MF-DFA and hXYq is called the generalized Hurst exponent [6]. Let hXq and hYq be the generalized Hurst exponents of Xt and Yt estimated using MF-DFA, respectively. For binomial multifractal measures and some multifractal random walks, it is validated numerically in the case of dependent pairs [32, 51]:(6)hXYq=hXq+hYq2.Furthermore, more general relationship is derived [37]:(7)hXYq≤hXq+hYq+logsk2,where k is a positive constant. When scale s⟶+∞,(8)hXYq≤hXq+hYq2.Since the computation ofFXYq,s in equation (3) is determined by FXYv,s from the segments with large fluctuation for q>0, hXYq of q>0 shows the scaling behavior of the segments with large fluctuations. For q<0, hXYq shows that with small fluctuations [6, 17].The joint multifractal spectrumfXYα is obtained via the Legendre transform of Renyi scaling exponent τXYq:(9)α=τXY′q,fXYα=qα−τXYq,where α is the singularity strength and τXYq can be obtained by calculating τXYq=qhXYq−1 [35, 38]. The strength of joint multifractality of records can be characterized by ΔαXY=αmax−αmin.The dynamics of multifractal system are described by a continuous spectrum of exponents (multifractal spectrum), rather than a single exponent (fractal dimension). This indicates that the exponentα (singularity strength or Hölder exponent) values continuously in one range. The width of this range is the width of multifractal spectrum ΔαXY. It can reflect the complexity of exponent value and characterize the strength of multifractality. Inherently, the multifractal formalism is originally used to characterize the complexity of invariant measures of nonlinear dynamical systems.Ifhq satisfies the generalized binomial multifractal model [18](10)hXYq=1q−lnaq+bqqln2,then the τq and Δα can be calculated via the following formulas:(11)τXYq=−lnaq+bqln2,ΔαXY=lnb−lnaln2.


### 2.2. Testing Method for Source of Joint Multifractality

This section focuses on the testing for the source of joint multifractality. As shown in previous articles [33, 44, 47, 48, 51], the source of multifractality can be clarified by analyzing the corresponding surrogate series. Unlike the source testing of single series, we need to consider cross-correlations between series which also play an important role in generating joint mulitfractality [33, 34]. In view of the above, we propose two hypothesis tests using pairs of surrogate series for the following two different sources:(1)
The different correlations within and between original series.(2)
The fatness of probability distribution of original series.The following two surrogate methods [51] can generate pairs of surrogate series which keep one of the two sources and destroy another. Let Xt and Yt be the two series to be analyzed (t=1,…,L).(1)
For testing Type (1) source, the pair of surrogate seriesXsg1t,Ysg1t|t=1,2,…,L is generated via shuffling the two series Xt and Yt independently.(2)
For testing Type (2) source, the pair of surrogate seriesXsg2t,Ysg2t|t=1,2,…,L is generated using rank-ordered remapping technique: let Zt|t=1,2,…,L be a sequence of random numbers generated from the standard normal distribution. We rearrange Zt such that the rearranged series ZXt|t=1,2,…,L has the same rank ordering as the original series Xt. The series ZYt|t=1,2,…,L is generated using the same technique for the original series Yt. Let Xsg2t,Ysg2t=σXZXt,σYZYt, where σX and σY are the sample standard deviations of Xt and Yt, respectively.The pair of surrogate seriesXsg1t,Ysg1t preserves the distribution of original series and removes the correlations within and between original series. If the joint multifractality is caused in part by different correlations within and between original series, their joint multifractality will be significantly different from that of two original series. The pair of surrogate series Xsg2t,Ysg2t preserves the correlations within and between original series, and both have the Gaussian distribution. If the joint multifractality is caused in part by the fat-tailed distribution of original series, they will show a weaker joint multifractality than original series. Then, the two types of joint multifractality can be distinguished by comparing hXYq of two original series with that of corresponding pair of surrogate series.Following the idea in [40], we propose two hypothesis tests for the source of joint multifractality, which are based on the empirical distribution of scaling exponent function estimated from 1000 pairs of surrogate series. It can not only achieve more robust results than the previous clarification method using one pair of surrogate series but also show more details about the source of joint multifractality such as the effects of large and small fluctuations.The null hypotheses for the two hypothesis tests are described as follows:(i)
H0sg1: the joint multifractality is not due to the correlations within and between series.(ii)
H0sg2: the joint multifractality is not due to the fatness of probability distribution.Monte Carlo simulation is adopted for the acceptance region of null hypotheses [52]. The proposed procedure for the two hypothesis tests is described as follows:(i)
Assume that the significance level of test isα. Generate 1000 pairs of surrogate series and calculate the scaling exponent function hsg1q (or hsg2q) using MF-DCCA to obtain the empirical probability distribution of the scaling exponent function, where hsg1q and hsg2q denote the scaling exponent function estimated from Xsg1t,Ysg1t and Xsg2t,Ysg2t, respectively.(ii)
Obtain the acceptance regionhlqi,huqi for null hypothesis at the significance level αi=α/N, based on the empirical probability distribution. In detail, hlqi is estimated by the (αi/2)-quantile of the empirical probability distribution and huqi is estimated by the (1−αi/2)-quantile, i=1,…,N, where N is the number of the analyzed qs.(iii)
Connecthlqis and huqis, respectively, using straight line. Since hlqis or huqis varies continuously, it is suitable to connect using straight line. The area surrounded by these connection lines from q1 to qN is the acceptance region.(iv)
Calculate the scaling exponent functionhXYq of original series using MF-DCCA.(v)
Accomplish the hypothesis testing via comparing thehXYq with the acceptance region.Take testing for the null hypothesisH0sg1 as an example. If hq falls into the acceptance region for H0sg1, we accept that the joint multifractality is not mainly due to different correlations within and between series. If hq falls out of the acceptance region for H0sg1, we believe that the joint multifractality is caused in part by different correlations within and between series.Note that the significance level for eachq is αi=α/N because this procedure is a multiple-testing procedure so that Bonferroni correction is conducted. In this paper, N=31 and we use 5% as its significance level. So the significance level for each q is αi≈0.16%.


### 2.3. Simultaneously Half-Shifting Technique

As shown in [33], the shift of two series does not change the correlations within each series but can weaken cross-correlations between them. For further analyzing the effects of cross-correlations on joint multifractality, we propose the following technique in MF-DCCA for both minimizing the cross-correlations between original series and making full use of records:(i)
LetXt and Yt be the two series to be analyzed (t=1,…,L). tm:=intL/2.(ii)
LetY0t=Y1,…,Ytm, Y1t=Ytm+1,…,YL and X0t=X1,…,XL−tm, X1t=XL−tm+1,…,XL. Two constructed pairs of shifting series are X0t,Y1t and X1t,Y0t. They keep the autocorrelations in each series and minimizing the cross-correlations between series.(iii)
For making full use of records, the joint multifractality without cross-correlations is obtained by simultaneously conducting two pairs of shifting series using MF-DCCA. For details, in the step of computation ofFXYq,s, for each s, the FXYq,ss is computed via averaging the fluctuations of both X0t,Y1t and X1t,Y0t. The other steps are the same.This technique simultaneously uses two pairs of series that are shifted half their length relative to each other. We thus call it simultaneously half-shifting technique (SHST). The scaling exponent functionhXYSq obtained via simultaneously half-shifting technique is compared with the scaling exponent function estimated from two original series and the acceptance region for H0sg1.For MF-DFA of single seriesXt, there are no cross-correlations. In this case, the simultaneously half-shifting technique degenerates the operation of half-dividing; that is, the Xt is divided into two series X0t and X1t which are defined above. Further for multifractal detrended fluctuation analysis of the two divided series, FXq,ss is computed via averaging the fluctuations of both X0t and X1t. The operation of half-dividing should not change the multifractality since it can keep the autocorrelations in series and the distribution of series. We can analyze the change in multifractality of each series (Xt and Yt) after the operation of half-dividing for validating the simultaneously half-shifting technique in keeping autocorrelations.


## 2.1. Multifractal Detrended Cross-Correlation Analysis

In this section, we introduce the MF-DCCA [32, 49] which is used for joint multifractal analysis of river level fluctuations in this work. The MF-DCCA is a generalization of the MF-DFA [6] for two nonstationary time series. The method is described as follows.LetXt and Yt be the two processes to be analyzed. First, the profile is computed as follows:(1)X′i≔∑k=1iXk−X¯,Y′i≔∑k=1iYk−Y¯,i=1,…,L,where L is the minimum length of Xt and Yt and X¯ and Y¯ are the means of Xt and Yt, respectively.Dividing the profilesX′t and Y′t into Ns:=intL/s nonoverlapping segments with equal length s and then computing the detrended covariance for each segment, we get(2)FXYv,s=1s∑i=1sX′v−1s+i−xviY′v−1s+i−yvi,v=1,…,Ns,where xvi and yvi are, respectively, the best polynomial fit of X′v−1s+i and Y′v−1s+i in each of the Ns segments v. In this step, the trends are removed via the polynomial fit of profile. The m-order polynomial used in the fitting procedure can eliminate the polynomial trend of order m−1 in the original series. This ensures the fluctuation analysis available for data affected by trends or other nonstationarities. Note that the detrended covariance computed by equation (2) uses absolute values of fluctuations for fluctuation strength, which is different from that of original paper [32] using absolute values of FXYv,s for fluctuation itself. Besides, the operation of absolute value here can also avoid no obvious power-law scaling when q<0 [32, 37], which exists in this study.Average the fluctuation functionFXYv,s over all the segments, given by(3)FXYq,s≔1Ns∑v=1NsFXYv,sq/21/q,where q can be any real values except 0. For q=0,(4)FXY0,s≔exp12Ns∑v=1NslnFXYv,s.The final step is calculatinghXYq via the slope of the log-log plot of FXYq,s versus s, which is based on the power-law [32]:(5)FXYq,s∼shXYq,for larges.The functionhXYq is called the scaling exponent function [32]. When q=2, hXYq is the bivariate Hurst exponent [50].WhenX=Y, MF-DCCA degenerates to MF-DFA and hXYq is called the generalized Hurst exponent [6]. Let hXq and hYq be the generalized Hurst exponents of Xt and Yt estimated using MF-DFA, respectively. For binomial multifractal measures and some multifractal random walks, it is validated numerically in the case of dependent pairs [32, 51]:(6)hXYq=hXq+hYq2.Furthermore, more general relationship is derived [37]:(7)hXYq≤hXq+hYq+logsk2,where k is a positive constant. When scale s⟶+∞,(8)hXYq≤hXq+hYq2.Since the computation ofFXYq,s in equation (3) is determined by FXYv,s from the segments with large fluctuation for q>0, hXYq of q>0 shows the scaling behavior of the segments with large fluctuations. For q<0, hXYq shows that with small fluctuations [6, 17].The joint multifractal spectrumfXYα is obtained via the Legendre transform of Renyi scaling exponent τXYq:(9)α=τXY′q,fXYα=qα−τXYq,where α is the singularity strength and τXYq can be obtained by calculating τXYq=qhXYq−1 [35, 38]. The strength of joint multifractality of records can be characterized by ΔαXY=αmax−αmin.The dynamics of multifractal system are described by a continuous spectrum of exponents (multifractal spectrum), rather than a single exponent (fractal dimension). This indicates that the exponentα (singularity strength or Hölder exponent) values continuously in one range. The width of this range is the width of multifractal spectrum ΔαXY. It can reflect the complexity of exponent value and characterize the strength of multifractality. Inherently, the multifractal formalism is originally used to characterize the complexity of invariant measures of nonlinear dynamical systems.Ifhq satisfies the generalized binomial multifractal model [18](10)hXYq=1q−lnaq+bqqln2,then the τq and Δα can be calculated via the following formulas:(11)τXYq=−lnaq+bqln2,ΔαXY=lnb−lnaln2.


## 2.2. Testing Method for Source of Joint Multifractality

This section focuses on the testing for the source of joint multifractality. As shown in previous articles [33, 44, 47, 48, 51], the source of multifractality can be clarified by analyzing the corresponding surrogate series. Unlike the source testing of single series, we need to consider cross-correlations between series which also play an important role in generating joint mulitfractality [33, 34]. In view of the above, we propose two hypothesis tests using pairs of surrogate series for the following two different sources:(1)
The different correlations within and between original series.(2)
The fatness of probability distribution of original series.The following two surrogate methods [51] can generate pairs of surrogate series which keep one of the two sources and destroy another. Let Xt and Yt be the two series to be analyzed (t=1,…,L).(1)
For testing Type (1) source, the pair of surrogate seriesXsg1t,Ysg1t|t=1,2,…,L is generated via shuffling the two series Xt and Yt independently.(2)
For testing Type (2) source, the pair of surrogate seriesXsg2t,Ysg2t|t=1,2,…,L is generated using rank-ordered remapping technique: let Zt|t=1,2,…,L be a sequence of random numbers generated from the standard normal distribution. We rearrange Zt such that the rearranged series ZXt|t=1,2,…,L has the same rank ordering as the original series Xt. The series ZYt|t=1,2,…,L is generated using the same technique for the original series Yt. Let Xsg2t,Ysg2t=σXZXt,σYZYt, where σX and σY are the sample standard deviations of Xt and Yt, respectively.The pair of surrogate seriesXsg1t,Ysg1t preserves the distribution of original series and removes the correlations within and between original series. If the joint multifractality is caused in part by different correlations within and between original series, their joint multifractality will be significantly different from that of two original series. The pair of surrogate series Xsg2t,Ysg2t preserves the correlations within and between original series, and both have the Gaussian distribution. If the joint multifractality is caused in part by the fat-tailed distribution of original series, they will show a weaker joint multifractality than original series. Then, the two types of joint multifractality can be distinguished by comparing hXYq of two original series with that of corresponding pair of surrogate series.Following the idea in [40], we propose two hypothesis tests for the source of joint multifractality, which are based on the empirical distribution of scaling exponent function estimated from 1000 pairs of surrogate series. It can not only achieve more robust results than the previous clarification method using one pair of surrogate series but also show more details about the source of joint multifractality such as the effects of large and small fluctuations.The null hypotheses for the two hypothesis tests are described as follows:(i)
H0sg1: the joint multifractality is not due to the correlations within and between series.(ii)
H0sg2: the joint multifractality is not due to the fatness of probability distribution.Monte Carlo simulation is adopted for the acceptance region of null hypotheses [52]. The proposed procedure for the two hypothesis tests is described as follows:(i)
Assume that the significance level of test isα. Generate 1000 pairs of surrogate series and calculate the scaling exponent function hsg1q (or hsg2q) using MF-DCCA to obtain the empirical probability distribution of the scaling exponent function, where hsg1q and hsg2q denote the scaling exponent function estimated from Xsg1t,Ysg1t and Xsg2t,Ysg2t, respectively.(ii)
Obtain the acceptance regionhlqi,huqi for null hypothesis at the significance level αi=α/N, based on the empirical probability distribution. In detail, hlqi is estimated by the (αi/2)-quantile of the empirical probability distribution and huqi is estimated by the (1−αi/2)-quantile, i=1,…,N, where N is the number of the analyzed qs.(iii)
Connecthlqis and huqis, respectively, using straight line. Since hlqis or huqis varies continuously, it is suitable to connect using straight line. The area surrounded by these connection lines from q1 to qN is the acceptance region.(iv)
Calculate the scaling exponent functionhXYq of original series using MF-DCCA.(v)
Accomplish the hypothesis testing via comparing thehXYq with the acceptance region.Take testing for the null hypothesisH0sg1 as an example. If hq falls into the acceptance region for H0sg1, we accept that the joint multifractality is not mainly due to different correlations within and between series. If hq falls out of the acceptance region for H0sg1, we believe that the joint multifractality is caused in part by different correlations within and between series.Note that the significance level for eachq is αi=α/N because this procedure is a multiple-testing procedure so that Bonferroni correction is conducted. In this paper, N=31 and we use 5% as its significance level. So the significance level for each q is αi≈0.16%.


## 2.3. Simultaneously Half-Shifting Technique

As shown in [33], the shift of two series does not change the correlations within each series but can weaken cross-correlations between them. For further analyzing the effects of cross-correlations on joint multifractality, we propose the following technique in MF-DCCA for both minimizing the cross-correlations between original series and making full use of records:(i)
LetXt and Yt be the two series to be analyzed (t=1,…,L). tm:=intL/2.(ii)
LetY0t=Y1,…,Ytm, Y1t=Ytm+1,…,YL and X0t=X1,…,XL−tm, X1t=XL−tm+1,…,XL. Two constructed pairs of shifting series are X0t,Y1t and X1t,Y0t. They keep the autocorrelations in each series and minimizing the cross-correlations between series.(iii)
For making full use of records, the joint multifractality without cross-correlations is obtained by simultaneously conducting two pairs of shifting series using MF-DCCA. For details, in the step of computation ofFXYq,s, for each s, the FXYq,ss is computed via averaging the fluctuations of both X0t,Y1t and X1t,Y0t. The other steps are the same.This technique simultaneously uses two pairs of series that are shifted half their length relative to each other. We thus call it simultaneously half-shifting technique (SHST). The scaling exponent functionhXYSq obtained via simultaneously half-shifting technique is compared with the scaling exponent function estimated from two original series and the acceptance region for H0sg1.For MF-DFA of single seriesXt, there are no cross-correlations. In this case, the simultaneously half-shifting technique degenerates the operation of half-dividing; that is, the Xt is divided into two series X0t and X1t which are defined above. Further for multifractal detrended fluctuation analysis of the two divided series, FXq,ss is computed via averaging the fluctuations of both X0t and X1t. The operation of half-dividing should not change the multifractality since it can keep the autocorrelations in series and the distribution of series. We can analyze the change in multifractality of each series (Xt and Yt) after the operation of half-dividing for validating the simultaneously half-shifting technique in keeping autocorrelations.


## 3. Data and Preliminaries

The river level records we consider have been analyzed for single multifractality in previous papers [39, 40]. The records are collected in real time from 10 observations along a tributary of Haihe River in North China. The locations of 10 observations are shown in Figure 1. We mark these locations from upstream to downstream with Site 1-Site 10. The records are of high-frequency, collected every 6 minutes at 10 water level observations of this river from April 2011 to September 2013. The location of river has a significant continental monsoon climate, which is cold and dry in the winter and warm and humid in the summer. More information about records can be seen in Section 3 of [40] and Table 1. The sudden change in records in Table 1 at Site 2 and Site 7 is led by water storage and discharge at the two sites. Besides, there are some tributaries located downstream of Site 7 that flows into this river. The main tributary of this river is located at Site 9.Figure 1
Site 1-Site 10 locations.Table 1
Information about records.
Statistics for daily average obs.Obs.DaysMeanMedianStd.Max.Min.Site 1182855808481.132481.0390.314482.826480.495Site 21523076681.0761.2480.3821.6220.073Site 31899678091.2751.0700.5772.6150.060Site 41909598091.4491.2480.6222.9020.230Site 51880748061.0661.1750.3551.5220.010Site 61855577851.1091.0560.5052.1650.100Site 71908408090.2160.1190.2171.0900.020Site 81904458090.5140.2360.5492.4710.084Site 91858228011.8851.8810.6112.9690.292Site 101116005151.4951.4140.6032.5040.530The preprocessing of records considers three issues: (1) the noise effects led by high-frequency data; (2) the missing data led by abnormal work of data-collecting equipment; and (3) seasonal trend in the records led by monsoon climate. Following the literature [40], we conduct the preaverage of one day’s record and the linear interpolation, respectively, for the first two issues. The preaverage of one day’s record allows us to use the daily average river levels for analysis. More detail about the preprocessing can be seen in Section 3 of [40].For the third issue, it is known that MF-DCCA cannot remove the seasonal periodic trends in the records, which cause cross-overs in fluctuation function and further give incorrect results [53–58]. There are some robust methods for the time series with periodic trends such as Fourier detrended fluctuation analysis (F-DFA) [54, 55], singular value decomposition (SVD) [56, 57, 59], and the polynomial of varying order l [60]. The method of SVD can achieve the same goal as F-DFA [36]. A useful review on this issue can refer to the literature [58]. In this work, we use F-DFA for removing the seasonal periodic trends. After truncating the first 38 coefficients of low Fourier frequencies, we obtain the preprocessed daily average river levels for analysis, which show no cross-overs in fluctuation function (see Figure 2).Figure 2
Log-scaling plots ofFXYq,s versus s for two preprocessed daily average river levels at adjacent sites. The scaling range of s is selected from 25.1 to 27.2. The range of q is chosen every 0.4 from −6 to 6, which corresponds to lines from bottom to top in each subfigure.Note that since we calculate joint multifractality of two records, the two records used for analysis should keep the same in time. This results in less available data than the sampling days.


## 4. Results and Discussion

### 4.1. Multifractal Results

In this section, we study the joint multifractality between the preprocessed daily average river levels using MF-DCCA and mainly focus on the change in joint multifractality along the river. So the joint multifractal properties of two records at Sitek and Site k+1 are obtained, where k=1,…,9. Parameters for MF-DCCA are selected following [39, 40]. They have been proved to be robust and useful for the same records. The third-order polynomial is adopted to calculate the best polynomial fit of the profile in each segment, which can eliminate the second-order polynomial trend in the original series. Limited by the series length and following [39, 40], the scaling range of s is selected from 25.1 to 27.2 every 0.1 power.The length of river levels is short and only about 800. For avoiding inaccurate results at largeq caused by the finite-size effects [61–64], the q-range cannot be large. The range of q is chosen −6≤q≤6 carefully. The range is large enough to contain accurate hXYqs for our analysis.Figure2 shows the log-scaling plots of FXYq,s versus s for two preprocessed daily average river levels at Site k and Site k+1. It can be seen that although there is some volatility in log2FXYq,s which is led by short series, log2FXYq,s increases approximately linearly with log2s for all qs at all adjacent sites. Therefore, there exists the cross power-law behavior of two preprocessed daily average river levels at adjacent sites.The joint multifractal results between two preprocessed daily average river levels at Sitek and Site k+1 are shown in Figure 3. The errors bars of hXYq are based on the linear least square fit. It is shown in Section 2.1 that there exist some relationships between hXYq and hXq+hYq/2. We also plot the values of hXq+hYq/2 in Figure 3 for comparison. From Figure 3, we can find the values of hXq+hYq/2 deviate obviously from the values of hXYq and fall out of the errors bars for q<0 at all adjacent sites except Site 3-Site 4 and Site 9-Site 10. The deviations of all pairs of adjacent sites become smaller as they get away from Site 2 and Site 7. As mentioned in Section 2.1, it has been validated numerically in the case of dependent pairs [32, 51]:(12)hXYq=hXq+hYq2.Figure 3
Comparison of values ofhXYq and hXq+hYq/2 for two preprocessed daily average river levels at adjacent sites. Different colors and symbols indicate the corresponding different results: blue triangles for values of hXYq and red circles for values of hXq+hYq/2. The curves are obtained by fitting the two-parameter binomial model (equation (10)) to points of hXYq. The resulting model parameters are reported in Table 2. The errors bars of hXYq are based on the linear least square fit. All the errors bars reveal 95% confidence interval.For dependent pairs, the values ofhXq+hYq/2 should not differ much from that of hXYq. So the deviations of hXYq from hXq+hYq/2 can reflect the strength of cross-correlations. The small deviations correspond to the strong cross-correlations. The results shown in Figure 3 are consistent with this conclusion and indicate weak cross-correlations between two daily average river levels of these adjacent sites except Site 3-Site 4 and Site 9-Site 10. In fact, we note that there are water storage and discharge activities at Site 2 and Site 7, which can lead to the weak cross-correlations between them and their upstream or downstream (Site 1, Site 3, Site 6, and Site 8). Site 4 and Site 5 are close to Site 7 (see Figure 1) and also affected by this activity. The hXYq is usually used to characterize the joint multifractality. Thus, we conclude that the human activities (water storage and discharge) at Site 2 and Site 7 can change the corresponding joint multifractality of river levels (the part caused by cross-correlations) by weakening cross-correlations. The significant deviations at Site 8-Site 9 may be due to the lateral inflow from the main tributary at Site 9 which can also weaken cross-correlations.It can be also seen in Figure3 that the shapes of hXYq for all pairs of adjacent sites suggest the two-parameter binomial model. The fitting results are shown in Figure 3 and reported in Table 2. The results show that the hXYq given by MF-DCCA is fitted well with two-parameter binomial model (see equation (10)), and the values of R2 are all larger than 0.95. For two-parameter binomial model, the strength of joint multifractality can be characterized by Δα=lnb−lna/ln2. From Table 2, the values of Δα are larger than 0.6, which indicates that there is strong joint multifractality at all adjacent sites.Table 2
Fitting results of binomial multifractal model.
hXY2hXY2 errors barsabΔαR2Site 1-site 20.0540.0210.6881.1560.7480.993Site 2-site 30.0500.0160.7131.1520.6920.996Site 3-site 40.0700.0150.5351.2071.1730.979Site 4-site 50.0740.0190.6361.1740.8840.974Site 5-site 60.0670.0160.7281.1350.6400.998Site 6-site 70.0620.0170.6751.1670.7910.994Site 7-site 80.0530.0170.6011.2031.0020.993Site 8-site 90.0550.0200.6271.1900.9240.995Site 9-site 100.1260.0270.6471.1310.8060.999


### 4.2. Source of Joint Multifractality

In this section, we test the source of joint multifractality obtained from river level records at pairs of adjacent sites in the section above. The testing method is described in Section2.2. All the testing results are based on 1000 pairs of surrogate series.Figure4 shows the testing results for the null hypothesis H0sg1: the joint multifractality is not due to the correlations within and between series. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series Xsg1t,Ysg1t. For testing of the null hypothesis H0sg1, if the joint multifractality is only due to the correlations within and between series, the pair of surrogate series Xsg1t,Ysg1t should show monofractality with hXYq=0.5, and the acceptance region should be around 0.5 for all qs. In this case, our testing method will be invalid for hXYq near 0.5 since hXYq near 0.5 will always fall into the acceptance region and accept H0sg1. This suggests that we should focus on the hXYq away from 0.5 when conducting this testing. From Figure 4, triangles indicate hXYqs of two preprocessed daily average river levels at adjacent sites. It can be seen that the acceptance region is always around 0.5 and most values of hXYq fall out of the shadow area at all pairs of adjacent sites except those near 0.5; that is, we should reject H0sg. So we can conclude that the joint multifractality is mainly due to the correlations within and between series. Note that most values of hXYq are smaller than the acceptance region, which indicates strong negative correlations within and between records.Figure 4
Testing via pair of surrogate seriesXsg1t,Ysg1t. X-axis denotes the range of q values. Triangles indicate hXYSqs of two preprocessed daily average river levels at adjacent sites. Circles indicate hXYSqs obtained via simultaneously half-shifting technique. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series for the null hypothesis H0sg1: the joint multifractality is not due to the correlations within and between series.For further analyzing the effects of cross-correlations on joint multifractality, we obtain thehXYSq (see circles in Figure 4) via simultaneously half-shifting technique which is described in Section 2.3 for minimizing the cross-correlations between original series. It can be seen that the difference between hXYSq obtained via simultaneously half-shifting technique and hXYq without this technique is large at Site 3-Site 4, Site 4-Site 5, Site 8-Site 9, and Site 9-Site 10. This indicates that the cross-correlations have significant effects on the joint multifractality at those pairs of adjacent sites. The insignificant effects of the cross-correlations on joint multifractality for Site 1-Site 2, Site 2-Site 3, Site 5-Site 6, Site 6-Site 7, and Site 7-Site 8 are consistent with the previous analysis of deviations between hXYq and hXq+hYq/2. It is due to the weak cross-correlations caused by human activities (water storage and discharge) at Site 2 and Site 7. We can also see that the values of hXYSq obtained via simultaneously half-shifting technique fall out of the shadow area at all pairs of adjacent sites except those near 0.5. It indicates that the autocorrelations within each series also have significant effects on the joint multifractality at all pairs of adjacent sites. Note that most values of hXYSq is smaller than the acceptance region, which indicates strong negative autocorrelations in records. In summary, for Site 1-Site 2, Site 2-Site 3, Site 5-Site 6, Site 6-Site 7, and Site 7-Site 8, the joint multifractality is mainly caused by the autocorrelations within each series; For other pairs of adjacent sites, both autocorrelations and cross-correlations have significant effects on the joint multifractality.Besides, we also compare the values of the generalized Hurst exponenthXq, respectively, obtained from the preprocessed daily average river levels and its half-dividing series at all sites for validating the simultaneously half-shifting technique in keeping autocorrelations. The results are shown in Figure 5. It shows that the values of hXq change little after the operation of half-dividing. This validates the effectiveness of simultaneously half-shifting technique in keeping autocorrelations.Figure 5
Comparison of values of the generalized Hurst exponenthXq, respectively, obtained from the preprocessed daily average river levels and its half-dividing series. Different symbols indicate the corresponding different results: triangles for results of the preprocessed daily average river levels and circles for its half-dividing series. The errors bars for the preprocessed daily average river levels are based on the linear least square fit. All the errors bars reveal 95% confidence interval.Figure6 shows the testing results for the null hypothesis H0sg2: the joint multifractality is not due to the fatness of probability distribution. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series Xsg2t,Ysg2t. It can be seen that almost all values of hXYq fall out of the shadow area; that is, we should reject H0sg2. But the difference between hXYq and acceptance region is small. This confirms that the probability distribution of river levels also has effects on the joint multifractality, but the effects are small.Figure 6
Testing via pair of surrogate seriesXsg2t,Ysg2t. X-axis denotes the range of q values. Triangles indicate hXYSqs of two preprocessed daily average river levels at adjacent sites. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series for the null hypothesis H0sg2: the joint multifractality is not due to the fatness of probability distribution.


## 4.1. Multifractal Results

In this section, we study the joint multifractality between the preprocessed daily average river levels using MF-DCCA and mainly focus on the change in joint multifractality along the river. So the joint multifractal properties of two records at Sitek and Site k+1 are obtained, where k=1,…,9. Parameters for MF-DCCA are selected following [39, 40]. They have been proved to be robust and useful for the same records. The third-order polynomial is adopted to calculate the best polynomial fit of the profile in each segment, which can eliminate the second-order polynomial trend in the original series. Limited by the series length and following [39, 40], the scaling range of s is selected from 25.1 to 27.2 every 0.1 power.The length of river levels is short and only about 800. For avoiding inaccurate results at largeq caused by the finite-size effects [61–64], the q-range cannot be large. The range of q is chosen −6≤q≤6 carefully. The range is large enough to contain accurate hXYqs for our analysis.Figure2 shows the log-scaling plots of FXYq,s versus s for two preprocessed daily average river levels at Site k and Site k+1. It can be seen that although there is some volatility in log2FXYq,s which is led by short series, log2FXYq,s increases approximately linearly with log2s for all qs at all adjacent sites. Therefore, there exists the cross power-law behavior of two preprocessed daily average river levels at adjacent sites.The joint multifractal results between two preprocessed daily average river levels at Sitek and Site k+1 are shown in Figure 3. The errors bars of hXYq are based on the linear least square fit. It is shown in Section 2.1 that there exist some relationships between hXYq and hXq+hYq/2. We also plot the values of hXq+hYq/2 in Figure 3 for comparison. From Figure 3, we can find the values of hXq+hYq/2 deviate obviously from the values of hXYq and fall out of the errors bars for q<0 at all adjacent sites except Site 3-Site 4 and Site 9-Site 10. The deviations of all pairs of adjacent sites become smaller as they get away from Site 2 and Site 7. As mentioned in Section 2.1, it has been validated numerically in the case of dependent pairs [32, 51]:(12)hXYq=hXq+hYq2.Figure 3
Comparison of values ofhXYq and hXq+hYq/2 for two preprocessed daily average river levels at adjacent sites. Different colors and symbols indicate the corresponding different results: blue triangles for values of hXYq and red circles for values of hXq+hYq/2. The curves are obtained by fitting the two-parameter binomial model (equation (10)) to points of hXYq. The resulting model parameters are reported in Table 2. The errors bars of hXYq are based on the linear least square fit. All the errors bars reveal 95% confidence interval.For dependent pairs, the values ofhXq+hYq/2 should not differ much from that of hXYq. So the deviations of hXYq from hXq+hYq/2 can reflect the strength of cross-correlations. The small deviations correspond to the strong cross-correlations. The results shown in Figure 3 are consistent with this conclusion and indicate weak cross-correlations between two daily average river levels of these adjacent sites except Site 3-Site 4 and Site 9-Site 10. In fact, we note that there are water storage and discharge activities at Site 2 and Site 7, which can lead to the weak cross-correlations between them and their upstream or downstream (Site 1, Site 3, Site 6, and Site 8). Site 4 and Site 5 are close to Site 7 (see Figure 1) and also affected by this activity. The hXYq is usually used to characterize the joint multifractality. Thus, we conclude that the human activities (water storage and discharge) at Site 2 and Site 7 can change the corresponding joint multifractality of river levels (the part caused by cross-correlations) by weakening cross-correlations. The significant deviations at Site 8-Site 9 may be due to the lateral inflow from the main tributary at Site 9 which can also weaken cross-correlations.It can be also seen in Figure3 that the shapes of hXYq for all pairs of adjacent sites suggest the two-parameter binomial model. The fitting results are shown in Figure 3 and reported in Table 2. The results show that the hXYq given by MF-DCCA is fitted well with two-parameter binomial model (see equation (10)), and the values of R2 are all larger than 0.95. For two-parameter binomial model, the strength of joint multifractality can be characterized by Δα=lnb−lna/ln2. From Table 2, the values of Δα are larger than 0.6, which indicates that there is strong joint multifractality at all adjacent sites.Table 2
Fitting results of binomial multifractal model.
hXY2hXY2 errors barsabΔαR2Site 1-site 20.0540.0210.6881.1560.7480.993Site 2-site 30.0500.0160.7131.1520.6920.996Site 3-site 40.0700.0150.5351.2071.1730.979Site 4-site 50.0740.0190.6361.1740.8840.974Site 5-site 60.0670.0160.7281.1350.6400.998Site 6-site 70.0620.0170.6751.1670.7910.994Site 7-site 80.0530.0170.6011.2031.0020.993Site 8-site 90.0550.0200.6271.1900.9240.995Site 9-site 100.1260.0270.6471.1310.8060.999


## 4.2. Source of Joint Multifractality

In this section, we test the source of joint multifractality obtained from river level records at pairs of adjacent sites in the section above. The testing method is described in Section2.2. All the testing results are based on 1000 pairs of surrogate series.Figure4 shows the testing results for the null hypothesis H0sg1: the joint multifractality is not due to the correlations within and between series. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series Xsg1t,Ysg1t. For testing of the null hypothesis H0sg1, if the joint multifractality is only due to the correlations within and between series, the pair of surrogate series Xsg1t,Ysg1t should show monofractality with hXYq=0.5, and the acceptance region should be around 0.5 for all qs. In this case, our testing method will be invalid for hXYq near 0.5 since hXYq near 0.5 will always fall into the acceptance region and accept H0sg1. This suggests that we should focus on the hXYq away from 0.5 when conducting this testing. From Figure 4, triangles indicate hXYqs of two preprocessed daily average river levels at adjacent sites. It can be seen that the acceptance region is always around 0.5 and most values of hXYq fall out of the shadow area at all pairs of adjacent sites except those near 0.5; that is, we should reject H0sg. So we can conclude that the joint multifractality is mainly due to the correlations within and between series. Note that most values of hXYq are smaller than the acceptance region, which indicates strong negative correlations within and between records.Figure 4
Testing via pair of surrogate seriesXsg1t,Ysg1t. X-axis denotes the range of q values. Triangles indicate hXYSqs of two preprocessed daily average river levels at adjacent sites. Circles indicate hXYSqs obtained via simultaneously half-shifting technique. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series for the null hypothesis H0sg1: the joint multifractality is not due to the correlations within and between series.For further analyzing the effects of cross-correlations on joint multifractality, we obtain thehXYSq (see circles in Figure 4) via simultaneously half-shifting technique which is described in Section 2.3 for minimizing the cross-correlations between original series. It can be seen that the difference between hXYSq obtained via simultaneously half-shifting technique and hXYq without this technique is large at Site 3-Site 4, Site 4-Site 5, Site 8-Site 9, and Site 9-Site 10. This indicates that the cross-correlations have significant effects on the joint multifractality at those pairs of adjacent sites. The insignificant effects of the cross-correlations on joint multifractality for Site 1-Site 2, Site 2-Site 3, Site 5-Site 6, Site 6-Site 7, and Site 7-Site 8 are consistent with the previous analysis of deviations between hXYq and hXq+hYq/2. It is due to the weak cross-correlations caused by human activities (water storage and discharge) at Site 2 and Site 7. We can also see that the values of hXYSq obtained via simultaneously half-shifting technique fall out of the shadow area at all pairs of adjacent sites except those near 0.5. It indicates that the autocorrelations within each series also have significant effects on the joint multifractality at all pairs of adjacent sites. Note that most values of hXYSq is smaller than the acceptance region, which indicates strong negative autocorrelations in records. In summary, for Site 1-Site 2, Site 2-Site 3, Site 5-Site 6, Site 6-Site 7, and Site 7-Site 8, the joint multifractality is mainly caused by the autocorrelations within each series; For other pairs of adjacent sites, both autocorrelations and cross-correlations have significant effects on the joint multifractality.Besides, we also compare the values of the generalized Hurst exponenthXq, respectively, obtained from the preprocessed daily average river levels and its half-dividing series at all sites for validating the simultaneously half-shifting technique in keeping autocorrelations. The results are shown in Figure 5. It shows that the values of hXq change little after the operation of half-dividing. This validates the effectiveness of simultaneously half-shifting technique in keeping autocorrelations.Figure 5
Comparison of values of the generalized Hurst exponenthXq, respectively, obtained from the preprocessed daily average river levels and its half-dividing series. Different symbols indicate the corresponding different results: triangles for results of the preprocessed daily average river levels and circles for its half-dividing series. The errors bars for the preprocessed daily average river levels are based on the linear least square fit. All the errors bars reveal 95% confidence interval.Figure6 shows the testing results for the null hypothesis H0sg2: the joint multifractality is not due to the fatness of probability distribution. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series Xsg2t,Ysg2t. It can be seen that almost all values of hXYq fall out of the shadow area; that is, we should reject H0sg2. But the difference between hXYq and acceptance region is small. This confirms that the probability distribution of river levels also has effects on the joint multifractality, but the effects are small.Figure 6
Testing via pair of surrogate seriesXsg2t,Ysg2t. X-axis denotes the range of q values. Triangles indicate hXYSqs of two preprocessed daily average river levels at adjacent sites. The shadow area indicates acceptance region of hXYq at 5% significance level estimated from 1000 pairs of surrogate series for the null hypothesis H0sg2: the joint multifractality is not due to the fatness of probability distribution.


## 5. Conclusion

The joint multifractal analysis is usually conducted in two different variables for their cross-correlations but rarely used for two records of one variable collected at two different places. Besides, there are also few studies on the source of joint multifractality. In this paper, we focused on the two issues for river level records of a tributary of Haihe River in North China and proposed two hypothesis tests and simultaneously half-shifting technique for the source analysis of joint multifractality. This study is an extension of our previous work on the single multifractality of one record with the same data set.The joint multifractal analysis was conducted using MF-DCCA for two river level records at pairs of adjacent sites along the river. The obtainedhXYqs was compared with hXq+hYq/2 at each pair of adjacent sites. Results of comparison showed that human activities (water storage and discharge) can change the joint multifractality to deviate from the average of two single multifractality. The change is due to the weak cross-correlations caused by human activities. It provides an approach for detecting the change in joint multifractality caused by cross-correlation.For the source of joint multifractality, we proposed two hypothesis tests, which are based on empirical distributions of the scaling exponent function estimated from 1000 pairs of surrogate series. It can achieve more robust results than the previous clarification method using one pair of surrogate series. Results of source testings showed that the joint multifractality in river level records is mainly caused by the correlations within and between series and is also related to the probability distribution of river levels. The simultaneously half-shifting technique was proposed for further analyzing the effects of the cross-correlations. This technique can both minimize the cross-correlations between original series and make full use of records. The further analysis showed more details about the source of joint multifractality. For Site 1-Site 2, Site 2-Site 3, Site 5-Site 6, Site 6-Site 7, and Site 7-Site 8, the joint multifractality is mainly caused by the autocorrelations within each series. For other pairs of adjacent sites, both autocorrelations and cross-correlations have significant effects on the joint multifractality. These results are partly consistent with results of joint multifractal analysis.The multifractality of hydrologic dynamics and its source is an important topic in the fields of hydrology and meteorology. The hydrologic complex system has strong nonlinear correlations which cause the multifractal feature of dynamics. The source of multifractality is the direct influence factor of hydrologic dynamics. Thus, the join multifractal analysis of river levels in this study is useful for the theory and simulation of hydrological phenomena. Especially, the analysis of effects of human activities on joint multifractality can be a reference for the detection and control of human activities in hydrological and meteorological environment. From the research method, this study provides an approach for detecting the change in joint multifractality and some techniques for source analysis of joint multifractality. These methods can be extended to other applications of joint multifractal analysis.



---

*Source: 1532805-2020-12-07.xml*