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 } .