Zhao et al.[4] have introduced MF-DXA method based on the MF-DFA method[1],[2] and analysed the cross-correlation between two non-stationary series quantitatively.
The steps for the MF-DXA method are as follows.

1. Let $x(i)$ and $y(i)$ are two data series for $i = 1,2,\ldots,N$, of length $N$. The mean of
these series is calculated as $\bar{x} = \frac{1}{N}\sum_{i=1}^{N} x(i)$ and $\bar{y} = \frac{1}{N}\sum_{i=1}^{N} y(i)$ respectively.
Then accumulated deviation series for $x(i)$ and $y(i)$, are calculated as per the below equations respectively.
\begin{eqnarray}
X(i) = \sum_{k=1}^{i} [x(k)-\bar{x}], i = 1,2,\ldots,N \nonumber \\
Y(i) = \sum_{k=1}^{i} [y(k)-\bar{y}], i = 1,2,\ldots,N \nonumber
\end{eqnarray}
This subtraction of the mean($\bar{x}$ and $\bar{y}$) from the data series, is a standard way of removing noisy data from the input data series. The effect of this subtraction would be eliminated by the detrending in the fourth step [2].

2. Both $X(i)$ and $Y(i)$ are divided into $N_s$ non-overlapping segments, where $N_s$ = $int(N/s)$, $s$ is the length of the segment. In our experiment $s$ varies from $16$ as minimum to $512$ as maximum value in log-scale.

3. For each $s$, we denote a particular segment by $v$($v = 1,2,\ldots,N_s$).
For each segment least-square fit is performed to obtain the local trend of the particular segment~[4].
Here $x_v(i)$ and $y_v(i)$ denote the least square fitted polynomials for the segment $v$ in $X(i)$ and $Y(i)$ respectively. $x_v(i)$ and $y_v(i)$ are calculated as per the equations $x_v(i) = \sum_{k=0}^{m} {C_{x_k}}{(i)^{m-k}}$ and $y_v(i) = \sum_{k=0}^{m} {C_{y_k}}{(i)^{m-k}}$, where $C_{x_k}$ and $C_{y_k}$ are the $k$th coefficients of the fit polynomials with degree $m$. For fitting linear, quadratic, cubic or higher $m$-order polynomials may be used~[1],[2],[4]. For this experiment $m$ is taken as $1$ [4].

4. To detrend the data series, we have to subtract the polynomial fit from the data series. There is presence of slow varying trends in natural data series. To quantify the scale invariant structure of the variation around the trends, detrending is required.
Here for each $s$ and segment $v$, $v = 1,2,\ldots,N_s$, detrending is done by subtracting the least-square fits $x_v(i)$ and $y_v(i)$ from the part of the data series $X(i)$ and $Y(i)$ respectively, for the segment $v$. The covariance of the these residuals, denoted by $f^2_{xy}(s,v)$ for a particular $s$ and $v$, is then calculated as per the following equation [4].
\begin{eqnarray}
f^2_{xy}(s,v) = \frac{1}{s}\sum_{i=1}^{s} \{X[(v-1)s+i]-x_v(i)\} \times \nonumber \\
\{Y[(v-1)s+i]-y_v(i)\}, \nonumber
\end{eqnarray}
for each segment $v$, $v = 1,2,\ldots,N_s$.

5. Then the $q$th-order detrended covariance, denoted by $F_{xy}(q,s)$, is calculated by averaging $f^2_{xy}(s,v)$ over all the segments($v$) generated for a particular $s$ and $q$, as per the equation below [1],[2],[4].
\begin{eqnarray}
F_{xy}(q,s) = \left\{\frac{1}{N_s}\sum_{v=1}^{N_s} [f^2_{xy}(s,v)]^{\frac{q}{2}}\right\}^{\frac{1}{q}}, \nonumber
\end{eqnarray}
for $q \neq 0$ because in that case $\frac{1}{q}$ would blow up.

6. The above process is repeated for different values of $s \in 16,32,\ldots,512$ and it can be seen that for a specific $q$, $F_{xy}(q,s)$ increases with increasing $s$. If the series are long range power correlated,
the $F_{xy}(q,s)$ versus $s$ for a particular $q$, will show power-law behaviour as below [4].
\begin{eqnarray}
F_{xy}(q,s) \propto s^{h_{xy}(q)} \nonumber
\end{eqnarray}
If this kind of scaling exists, $\log_{2} [F_{xy}(q,s)]$ would depend linearly on $\log_{2} s$, where $h_{xy}(q)$ is the slope and represents the degree of the cross-correlation between the the data series $x(i)$ and $y(i)$.

In general, $h_{xy}(q)$ depends on $q$.
$q$ ranges from negative to positive values. $h_{xy}(q)$-values calculated with negative $q$-s, are influenced by scaling behaviour of the segments($v$) with smaller fluctuations, whereas, $h_{xy}(q)$-values calculated with positive $q$-s, are influenced by scaling behaviour of the segments($v$) with larger fluctuations. For $q = 2$, method boils down to standard method of DXA [4].

As confirmed from several experiments done by Zhao et al. [4], if $h_{xy}(q)=0.5$ then there is no cross-correlation. Further, if $h_{xy}(q) > 0.5$ then there is persistent long-range cross-correlations, where the large value of one variable is likely to be followed by a large value of another variable in the series, whereas, $h_{xy}(q) < 0.5$ then there is anti-persistent long-range cross-correlations, where a large value of one variable is most likely to be followed by a small value and vice versa, in the series. $h_{xy}(q)$ for $q=2$, i.e. $h_{xy}(2)$ is the DXA exponent. As per Podobnik et al. the cross-correlation exponent between two non-stationary series, denoted by $\gamma_x$, is calculated as per the equation $\gamma_x = 2-2\{h_{xy}(2)\}$ [3]. For uncorrelated data series, $\gamma_x=1$, and the lower the value of $\gamma_x$, the more correlated the data series are.

References
1. J. , S. , E.-B. , S. , A. and H. , "Multifractal detrended fluctuation analysis of nonstationary time series", "Physica A: Statistical Mechanics and its Applications" "316"(0000), no. "1", 87---114.
2. J. , S. , E.-B. , S. , A. and H. , "Multifractal detrended fluctuation analysis of nonstationary time series", "Physica A: Statistical Mechanics and its Applications" "316"(0000), no. "1", 87---114.
3. B. Podobnik and H. Stanley, Detrended cross-correlation analysis: A new method for analyzing two nonstationary time series, Physical Review Letters 100(2008), no. 8,
4. X. WANG and C. YANG, FRACTAL PROPERTIES OF PARTICLES IN PHASE SPACE FROM URQMD MODEL, International Journal of Modern Physics E 22(2013), no. 04, 1350021.