  We introduce a pipeline including multifractal detrended cross-correlation
analysis (MF-DXA) modified by either singular value decomposition or the
adaptive method to examine the statistical properties of the pulsar timing
residual ($PTR$) induced by a gravitational wave (GW) signal. We propose a new
algorithm, the so-called irregular-MF-DXA, to deal with irregular data
sampling. Inspired by the quadrupolar nature of the spatial cross-correlation
function of a gravitational wave background, a new cross-correlation function,
$\bar{\sigma}_{\times}$, derived from irregular-MF-DXA will be introduced. We
show that, this measure reveals the quadrupolar signature in the $PTRs$ induced
by stochastic GWB. We propose four strategies based on the $y$-intercept of
fluctuation functions, the generalized Hurst exponent, and the width of the
singularity spectrum to determine the dimensionless amplitude and power-law
exponent of the characteristic strain spectrum as
$\mathcal{H}_c(f)\sim\mathcal{A}_{yr}(f/f_{yr})^{\zeta}$ for stochastic GWB.
Using the value of Hurst exponent, one can clarify the type of GWs. We apply
our pipeline to explore 20 millisecond pulsars observed by Parkes Pulsar Timing
Array. The computed scaling exponents confirm that all data are classified into
a nonstationary class implying the universality feature. The value of the Hurst
exponent is in the range $H\in [0.56,0.87]$. The $q$-dependency of the
generalized Hurst exponent demonstrates that the observed $PTRs$ have
multifractal behavior, and the source of this multifractality is mainly
attributed to the correlation of data which is another universality of the
observed datasets. Multifractal analysis of available $PTRs$ datasets reveals
an upper bound on the dimensionless amplitude of the GWB, $\mathcal{A}_{yr}<
2.0\times 10^{-15}$.
