Response and Uncertainty of the Parabolic Variance PVAR to Non-Integer Exponents of the Power Law
11 Responses and Degrees of Freedom of PVARfor a Continuous Power-Law PSD
Franc¸ois Vernotte, Enrico Rubiola, Siyuan Chen
Abstract —This paper is devoted to the use of the ParabolicVariance (PVAR) to characterize continuous power-law noises.First, from a theoretical calculation, we give the responses ofPVAR extended to continuous power-law noises. Second, froman empirical study, we provide an approximated but relevantexpression of the uncertainties of PVAR estimates. This simpleexpression easily applies to continuous power-law noises.
Keywords —frequency stability; fractional noise; uncertaintyassessment; degrees of freedom
I. I
NTRODUCTION
We generally assume that the FM noises affecting an oscil-lator are five and can be modeled by a PSD following 5 integerpower-laws, from f − to f +2 . It is sometimes necessary toconsider fractional or even real exponents, whether due to theprocess itself or to the needs of the estimation. In our case, themillisecond pulsar timing analysis, we are faced to both cases:on the one hand, we are searching for a red noise which isexpected to follow a f − / FM slope [1], [2]; and on the otherhand, we aim to estimate the exponent of this red noise by aBayesian statistics analysis of the PVAR estimates. We havethen to consider this exponent as a real number around whichwe have to place a confidence interval. For these reasons,we need an expression of the PVAR response for continuouspower-law noises as well as its uncertainty, i.e. the number ofdegrees of freedom of a PVAR estimate versus the power-lawexponent. II. C
ONTINUOUS
PVAR
RESPONSE
The response of a variance XVAR (X may be A for AVAR,M for MVAR, H for HVAR, P for PVAR, etc.) for a f α noiseis given by the following relationshipXVAR ( τ ) = (cid:90) ∞ | H X ( f ) | h α f α (1)where H X ( f ) is the transfer function associated to XVAR and h α the f α noise level.Following the pioneer work of Walter [3], nothing forces α to be an integer and we can calculate the response of avariance for continuous power-law noise. Thus, using Eq. (1)leads to the continuous response of AVAR:AVAR ( τ, α ) = (cid:0) − α +1 − (cid:1) Γ( α −
1) sin( πα/ πτ ) α +1 h α (2) F. Vernotte, and E. Rubiola are with FEMTO-ST, Time and FrequencyDepartment, Observatory THETA, UBFC, Besanc¸on, FranceSiyuan Chen is with Station Radioastro. Nanc¸ay, PSL, Nanc¸ay, France,FEMTO-ST, T-F Department, UBFC, Besanc¸on, France and LPC2E, Univer-sit´e d’Orl´eans, Orl´eans, France π /3 26 π /35 2 ln(2)2[7-ln(16)]/5 1/2 3/5 V a r i an c e r e s pon s e Exponent α AVARPVAR
Fig. 1. Continuous response of AVAR and PVAR compared to the knownresponses for α ∈ {− , − , } . which is equal to [3, Eq. (14)] since we don’t consider a highcut-off frequency in (1).Similarly, by using the transfer function of PVAR [4, Eq.(17)], we foundPVAR ( τ, α ) = 9 · − α (cid:2) α − α − − α ( α − (cid:3) × Γ( α −
5) sin( πα/ πτ ) α +1 h α . (3)Since PVAR converges from f − to f +2 FM noise, we assumethat this expression is valid for α ∈ ] − , +3[ .Fig. 1 shows that this relationship gives results which arein perfect accordance with the PVAR responses for integerpower-law noises [4, Table I].III. D EGREES OF F REEDOM OF
PVAR
ESTIMATES
First, we have to find a simplified expression of the numberof degrees of freedom (dof) of PVAR estimates for integerpower-law noises. Since this equation was solved for a whitePM noise (see Eq. (24) in [4]), we assume that the followingexpression should have the same form: ν ≈ A ( α ) m/M − B ( α )( m/M ) (4)with m = τ /τ , i.e. the integration time τ divided by thesampling time τ , M is the number of summations in thecomputation of the variance estimates, i.e. M = N − m + 2 ( N is the total number of samples in the data-run) for PVARand A ( α ) , B ( α ) are coefficients to be determined for each α exponent. Thanks to [4, Eq. (24)], we already know that A (+2) = 23 and B (+2) = − . a r X i v : . [ phy s i c s . d a t a - a n ] M a y We determined these coefficients by using massive Monte-Carlo simulations and verified the results by comparing themto the dof computed for continuous power-law noises.
A. Determination of the coefficients from Monte-Carlo simu-lations
The Monte-Carlo simulation was performed by computing10 000 sequences of frequency deviations for each α ∈{− , − , , +1 , +2 } and for each data-run length N ∈{ , , } , i.e. 150 000 simulated sequences. Forgiven α , N and τ , we deduced the dof from the averagesand the variances of the PVAR for the corresponding set ofsequences by using the following well-known property of χ ν distributions [4]: ν = 2 E [ PVAR ( τ )] V [ PVAR ( τ )] , (5)where E [ · ] and V [ · ] stands respectively for the mathematicalexpectation and the variance of the quantity between thebrackets. By using a least square fit, we obtained the followingresults: • A ( − ≈ , A ( − ≈ , A (0) ≈ , A (+1) ≈ , A (+2) = 23 • B ( α ) ≈ for all α .We have then modeled A ( α ) by the following rd orderpolynomial and assumed that B ( α ) = B is constant: (cid:26) A ( α ) = 27 + α + α − α B = 12 . (6)
1 10 100 1000 10000 D eg r ee s o f F r eedo m ν Normalized integration time m f -2 FMf -1 FMWhite FMf FMf FM-70-60-50-40-30-20-10 0 10 20 1 10 100 1000 10000 E rr o r / % Normalized integration time m f -2 FMf -1 FMWhite FMf FMf FM Fig. 2. Above: comparison of the empirical dof (crosses) and the approxima-tions (lines) given by Eq. (4) and (6) for all types of noise. Below: relativedifference (in %) between the empirical dof and the approximations.
Thanks to Eq. (4) and (6), we are now able to assess the dofof all PVAR estimates whatever the integration time or thenumber of samples.The upper plot of Fig. 2 compares the dof obtained by theMonte-Carlo simulations and by Eq. (4) and (6) for all integertypes of noise. The pretty good agreement is confirmed bythe lower plot which shows that the discrepancies are within ± except for the very first m -values ( m = 1 , ).Therefore, the approximated dof assessment by using Eq.(4) and (6) can then be applied to all integer power-law noises.In order to check if it remains valid for non-integer power-lawnoises when α ∈ ] − , let us compute the continuous dof ofPVAR. B. Verification for continuous power-law noises
The dof may be computed from Eq. (5). The mathematicalexpectation is the response of PVAR given above in (3) andthe variance can be computed from (21) and (22) of [4]: V [ PVAR ( τ )] = 2 M M − (cid:88) i =0 M − (cid:88) j =0 (cid:20) m τ m − (cid:88) k =0 m − (cid:88) l =0 (cid:18) m − − k (cid:19) (cid:18) m − − l (cid:19)(cid:110) R x [( i + k − j − l ) τ ] − R x [( i + k − j − m − l ) τ ] − R x [( i + m + k − j − l ) τ ] (cid:111)(cid:21) (7)where R x ( τ ) is the autocorrelation function of the phase-timesamples x ( t ) = (cid:82) t y ( θ ) d θ , i.e. R x ( τ ) = E { x ( t ) x ( t + τ ) } .We used the following continuous expression of R x ( τ ) versusthe power-law exponent α (see [5], [3]): R x ( mτ ) = h α π ) α τ α − Γ( m − α/ α − m + α/ α/ − α/ . However, since this expression involves Γ -functions with argu-ments of the order of m , we have limited these computationsto N = 128 samples ( Γ(128) = 3 · !). Thanks tothis equation, we have computed the theoretical variance ofPVAR ( τ ) versus continuous α and deduced the dof from (5).Let us call P ν ( α, m, M ) = Mmν . From Eq. (4), we see that P ν ( α, m, M ) ≈ A ( α ) − Bm/M . The top plot of Fig. 3 shows P ν ( α, m, M ) computed from (7) (crosses) and approximatedfrom (6) (solid lines) versus the noise power-law α for m ∈ { , , } (we prefer to plot P ν ( α, m, M ) than ν forreasons of readability of the graph). The agreement is quitegood for m = 11 and m = 16 , but there is a notable differencefor m = 4 and α < − . The lower plot of Fig. 3 showsthat this discrepancy is around 20% at worst but generallyremains within ± % in most cases (all α for m > andall m for α > − ). This agreement is satisfactory to getan acceptable assessment of the PVAR uncertainties sincethe relative uncertainties are proportional to / √ ν : they aretherefore always below 10% and mostly within ± . . P ν ( α , m , M ) Noise power law exponent α m= 4m=11m=32-25-20-15-10-5 0 5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 E rr o r / % Noise power law exponent α m=4m=5m=6m=7m=8m=10m=11m=13m=16m=19m=23m=27m=32 Fig. 3. Above: comparison of P ν ( α, m, M ) computed from (7) ( × , + , ∗ )and approximated by (4) (solid lines) for N = 128 data. The blue squaresand the green circles are respectively the values obtained for m = 4 and m = 32 from the Monte-Carlo simulations. Below: Error (in %) between theapproximated values of P ν ( α, m, M ) and the computed values. C. Case of the largest integration times
The approximation given by Eq. (4) and (6) remains enoughclose to the empirical dof ν if m ≤ N/ . Moreover, we knowthat ν = 1 for m = N/ . On the other hand, we can note thatthe approximation diverges beyond N/ (see dashed lines inthe upper plot of Fig. 4). However, it is of importance to assessthe uncertainties within this interval, particularly if N is nota power of 2.In order to fill this gap, we decided to interpolate the dofwithin round (cid:0) / N/ (cid:1) ≤ m ≤ round (cid:0) − / N/ (cid:1) , i.e.between m ≈ round (1 . N/ and m ≈ round (0 . N/ ,by using the following semi-logarithmic fit ν ( m ) = a ln( m ) + b with a = ν ( m ) − m ) − ln( m ) b = ln( m ) − ν ( m ) ln( m )ln( m ) − ln( m ) . For m ≥ m , the dof are set to 1.Fig. 4 (top) is an enlargement of the highest 2 decades of m , i.e. m ∈ [4096 , fo N = 32768 data, to focus on theresult of the semi-logarithmic fit. The bottom plot shows theerror between the fit and the dof computed from the Monte-Carlo simulations. We can see than most of these errors arewithin ± except for the case of the white FM which rangesfrom +5% to − and even reaches -24% for m = 14 766 .However, this fit is sufficient to ensure an estimation of the D eg r ee s o f F r eedo m ν Normalized integration time mf -2 FMWhite FMf FM-30-25-20-15-10-5 0 5 10 100004000 20000 E rr o r / % Normalized integration time mf -2 FMf -1 FMWhite FMf FMf FM Fig. 4. Above: comparison of the empirical dof (crosses) and the semi-logarithmic fits (solid lines) for N/ ≤ m ≤ N/ and for random walk FM,white FM and white PM. The dashed lines represent the approximations givenby Eq. (4) and (6). In this example N = 32768 samples and the logarithmicincrement of the m -values is / within [ N/ , N/ . Below: Error (in %)between the empirical dof and the semi-logarithmic fits for all types of noise. PVAR uncertainty for the highest τ value within ∼ atworst. IV. C ONCLUSION
From a theoretical calculation, we have determined theresponse of PVAR for continuous power-law noises. FromMonte-Carlo simulations, we have got a simplified expressionproviding the dof of the the PVAR estimates within 10 %.We have proved that this expression remains valid for non-integer power-law noises. Finally, we have shown that a simpleinterpolation is efficient to fit the dof for the highest octaveof integration times. Thanks to these results, we will be ableto use PVAR to analyze millisecond pulsar timings and toestimate the non-integer exponent of a red noise if it isdetected. On the other hand, this paper make it possible togeneralize the use of PVAR to process any signal with non-integer power-law noise.A
CKNOWLEDGEMENT
This work was partially funded by the ANR Programmesd’Investissement d’Avenir (PIA) Oscillator IMP (Project 11-EQPX-0033) and FIRST-TF (Project 10-LABX-0048).R
EFERENCES[1] E. S. Phinney, “A practical theorem on gravitational wave backgrounds,”August 2001, arXiv:astro-ph/0108028. [2] S. Chen, A. Sesana, and W. Del Pozzo, “Efficient computation of thegravitational wave spectrum emitted by eccentric massive black hole bina-ries in stellar environments,”
Monthly Notices of the Royal AstronomicalSociety , vol. 470, no. 2, pp. 1738–1749, 05 2017.[3] T. Walter, “Characterizing frequency stability: A continuous power-lawmodel with discrete sampling,”
IEEE Transactions on Instrumentationand Measurement , vol. 43, no. 1, pp. 69–79, February 1994.[4] F. Vernotte, M. Lenczner, P.-Y. Bourgeois, and E. Rubiola, “The parabolicvariance (PVAR), a wavelet variance based on least-square fit,”