text:
# On comparing the writhe of a smooth curve to the writhe of an inscribed polygon

## 1. Introduction

The writhing number measures the wrapping and coiling of space curves. Writhe has proved useful in molecular biology, where it is used to study the geometry of tangled strands of DNA ; often with the famous Călugăreanu-White formula for a curve $`C`$ in space with a normal field $`V`$ :

$$\mathrm{Lk}(C,C+ϵV)=\mathrm{Tw}(C,V)+\mathrm{Wr}(C).$$

In these applications, and in numerical simulations performed by biologists and mathematicians, it is often required to compute writhing numbers using numerical methods.

Several authors have presented algorithms for computing the exact writhing number of an $`n`$-edge polygonal curve in a finite number of steps . The fastest of these algorithms runs in time $`O(n^{1.6})`$, while earlier methods use time $`O(n^2)`$.

Careful implementations of such algorithms provide acceptable accuracy in computing writhe for polygonal curves. But reliably computing the writhe of smooth curves requires another step: we must be able to bound the error introduced in approximating a smooth curve by an inscribed polygonal curve. The purpose of this paper is to prove:

###### Theorem 1.

Suppose $`C(t)`$ is a simple, closed curve of class $`𝒞^4`$. We assume $`C(t)`$ is parametrized so that $`|C^{}(t)|1`$, and that we have upper bounds $`B_1,\mathrm{},B_4`$ on $`|C^{}(t)|,\mathrm{}|C^{(4)}(t)|`$. Let $`C_n(t)`$ be any $`n`$-edge polygonal curve inscribed in $`C`$ with maximum edge length $`x`$ and $`1/x>5B_2`$.

If the ribbon formed by joining $`C_n(t)`$ to $`C(t)`$ for every $`t`$ is embedded,

(1)
$$|\mathrm{Wr}(C)\mathrm{Wr}(C_n)|<\alpha nx^3+nO(x^4).$$

where $`\alpha `$ is a numerical constant less than $`B_2(5B_2^2+B_3)`$.

That is, if the lengths of the edges of $`C_n`$ are approximately constant, the error is bounded by a multiple of $`1/n^2`$.

The proof is based on Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula, which gives the difference in writhing number between two curves as the spherical area of the ribbon bounded by the curves on $`S^2`$ swept out by their unit tangent vectors . (Following Bruce Solomon , we will refer to such curves as tantrices, though they are classically referred to as tangent indicatrices.)

We begin by defining the writhing number in Section 2. Sections 3 and 4 then introduce the original form of Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula. In Sections 5 and 6 we extend Fuller’s formula to the case where one curve is polygonal and the other is of class $`𝒞^2`$ using a natural geometric idea: the tantrix of a polygonal curve should be defined to be the chain of geodesic segments on $`S^2`$ joining the (isolated) tangent vectors of the curve (this was pointed out by Chern in ). In the process, we discover a surprising fact: the writhe of a polygonal curve is equal to the writhe of any smooth curve obtained by carefully rounding off its corners!

Section 7 contains the remainder of our work: estimating the terms in our improved version of the $`\mathrm{\Delta }\mathrm{Wr}`$ formula to obtain Theorem 9. We test our error bounds in Section 8 by computing the writhe of a collection of polygonal curves inscribed in a smooth curve of known writhe.

The last section contains a discussion of some open problems inspired by the present work. We state the most important of them now: Like most of the theory of writhing numbers, the proof of our main theorem depends essentially on the fact that $`C`$ is closed. Can these methods be extended to open curves?

## 2. Definitions

The writhing number of a space curve is defined by:

###### Definition 2.

The writhe of a piecewise differentiable curve $`C(s)`$ is given by:

(2)
$$\mathrm{Wr}(C)=\frac{1}{4\pi }_{C\times C}\frac{C^{}(s)\times C^{}(t)(C(s)C(t))}{|C(s)C(t)|^3}dsdt,$$

Definition 2 is inspired by the Gauss formula for the linking number of two space curves, $`A(s)`$ and $`B(s)`$ (see Epple  for a fascinating discussion of the history of this formula):

(3)
$$\mathrm{Lk}(A,B)=\frac{1}{4\pi }_{A\times B}\frac{A^{}(s)\times B^{}(t)(A(s)B(t))}{|A(s)B(t)|^3}dsdt.$$

When the two curves $`A`$ and $`B`$ become a single curve, their linking number becomes the writhing number. This introduces a potential singularity on the diagonal of $`C\times C`$, but a careful calculation shows that the integral still converges. In fact, the integrand of Equation 2 approaches $`0`$ on the diagonal of $`C\times C`$, even when the curve $`C`$ has a corner.

From now on, we’ll assume that $`C`$ is simple. With this assumption, another way to look at the integral of Definition 2 is to observe that the integrand is the pullback of the area form on $`S^2`$ under the Gauss map $`C\times CS^2`$ defined by

(4)
$$(C(s),C(t))\frac{C(s)C(t)}{|C(s)C(t)|}.$$

From this perspective, we can see that the (signed) multiplicity of the Gauss map at any point $`p`$ on $`S^2`$ is just the number of self-crossings of the projection of $`C`$ in direction $`p`$.

## 3. Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula

Suppose we have a differentiable curve $`C(t)`$, with unit tangent vector $`T(t)`$. As we mentioned in Section 1, the curve $`T(t)`$ on the unit sphere is known as the tantrix of $`C`$. This curve divides the unit sphere into a number of cells. Within each cell, the signed crossing number of the projection of $`C`$ is constant: changing projection directions within the cell amounts to altering the projection of the knot by a regular isotopy consisting of Reidemeister moves of type II and III (pictured below). Neither of these moves changes the signed crossing number of the knot.

Figure 1. Changing the projection direction within a cell can only alter the diagram by one of these two moves. Neither changes the signed crossing number of the diagram, as we can see by counting the $`+`$ and $``$ markers at the crossings of $`C`$.

This observation motivates the idea that the writhe of a closed space curve is related to the fraction of the sphere’s area enclosed by its tantrix. In 1978, Brock Fuller stated the following:

###### Theorem 3.

(Fuller’s Spherical Area Formula) For any closed space curve $`C(s)`$ of class $`𝒞^3`$, let $`A`$ be the spherical area enclosed by the tantrix of $`C`$. Then

(5)
$$1+\mathrm{Wr}(C)=\frac{A}{2\pi }mod2.$$

Fuller used this formula to conclude that the difference in writhe between two curves $`X_0`$ and $`X_1`$ whose tantrices $`T_0`$ and $`T_1`$ are sufficiently close is given by a certain formula, which represents the spherical area of the ribbon between $`T_0`$ and $`T_1`$.

To be more specific, suppose that $`X_0`$ and $`X_1`$ are simple closed space curves of class $`𝒞^2`$, with regular parametrization (that is, parametrized so that $`X_0^{}`$ and $`X_1^{}`$ never vanish), and unit tangent vectors $`T_0`$ and $`T_1`$. Let $`F:S^1\times [0,1]𝐑^3`$ be a continuous deformation of $`X_0`$ into $`X_1`$, where $`F(t,\lambda )=X_\lambda (t)`$ and the $`X_\lambda `$ are simple curves of class $`𝒞^1`$, with unit tangent vectors $`T_\lambda (t)`$ continuous in $`(t,\lambda )`$.

###### Theorem 4.

(Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ Formula) If $`T_1(t)`$ and $`T_\lambda (t)`$ are not antipodal for all $`(t,\lambda )`$, then

(6)
$$\mathrm{Wr}(X_1)\mathrm{Wr}(X_0)=\frac{1}{2\pi }_C\frac{T_0(t)\times T_1(t)}{1+T_0(t)T_1(t)}\left[T_0^{}(t)+T_1^{}(t)\right]dt.$$

We observe that this formula does not require an arc-length parametrization of $`X_0`$ and $`X_1`$.

## 4. Justifying Fuller’s interpretion of the $`\mathrm{\Delta }\mathrm{Wr}`$ formula

While Fuller stated both these theorems in 1978, he did not provide complete proofs for either. The first rigorous proofs of Theorems 3 and 4 were given by Aldinger, Tabor, and Klapper  in 1995. While these authors proved both theorems as stated, they did not show that the formula in Theorem 4 represents the spherical area of the ribbon between $`T_0`$ and $`T_1`$ (in , the right-hand side of Equation 6 describes the difference between the twist of two frames on $`X_0`$ and $`X_1`$.)

In the spirit of their paper, we now justify Fuller’s original intuition about Equation 6.

###### Proposition 5.

Given two curves $`T_0(t),T_1(t):[0,1]S^2`$ where $`T_0(t)`$ and $`T_1(t)`$ are never antipodal, the area of the spherical region $`R`$ bounded by $`T_0`$, $`T_1`$ and the great circle arcs joining their endpoints is given by

(7)
$$\mathrm{Area}(R)=\frac{T_0(t)\times T_1(t)}{1+T_0(t)T_1(t)}\left(T_0^{}+T_1^{}\right)dt.$$

###### Proof.

We let

$$u(\theta ,t)=\mathrm{cos}\theta T_0(t)+\mathrm{sin}\theta T_1(t),$$

and parametrize the region $`R`$ by

$$v(\theta ,t)=\frac{u(\theta ,t)}{|u(\theta ,t)|}$$

where $`\theta `$ ranges from $`0`$ to $`\pi /2`$. Plugging this parametrization into the area form on $`S^2`$, and using the properties of the triple product, we find

$$\mathrm{d}\mathrm{Area}=\frac{1}{|u|^3}\left(\frac{u}{\theta }\times \frac{u}{s}u\right)\mathrm{d}\theta \mathrm{d}t.$$

Using the definition of $`u(\theta ,t)`$, this simplifies to

$$\mathrm{d}\mathrm{Area}=T_0\times T_1\left(\frac{\mathrm{cos}\theta }{(1+2\mathrm{cos}\theta \mathrm{sin}\theta T_0T_1)^{\frac{3}{2}}}T_0^{}+\frac{\mathrm{sin}\theta }{(1+2\mathrm{cos}\theta \mathrm{sin}\theta T_0T_1)^{\frac{3}{2}}}T_1^{}\right)\mathrm{d}\theta \mathrm{d}t.$$

Using the formula $`\mathrm{sin}2\theta =2\mathrm{cos}\theta \mathrm{sin}\theta `$, and the fact that the definite integrals of each of the trigonometric expressions above from $`0`$ to $`\pi /2`$ are equal, we have

$$\mathrm{Area}(R)=_0^1T_0\times T_1\left[_0^{\pi /2}\frac{\mathrm{cos}\theta }{(1+\mathrm{sin}2\theta T_0T_1)^{3/2}}d\theta \right](T_0^{}+T_1^{})dt.$$

This can be solved by the general integration formula

(8)
$$\frac{\mathrm{cos}\theta }{(1+a\mathrm{sin}2\theta )^{3/2}}d\theta =\frac{a\mathrm{cos}\theta \mathrm{sin}\theta }{(a^21)\sqrt{1+a\mathrm{sin}2\theta }},$$

which yields the formula in the statement of the Proposition. ∎

## 5. Extending Fuller’s formula to polygonal curves: I

To measure the difference in writhe between a smooth curve and a polygonal curve inscribed in the smooth curve, we must extend Theorem 4 to polygonal curves. To do so, we intend to approximate each polygonal curve with a family of smooth curves so that the writhe of the smooth curves converges to the writhe of the polygonal curve.

Examining Definition 2, it might seem that this result follows from general principles. For instance, one might conjecture that $`\mathrm{Wr}`$ was continuous in the $`C^1`$ norm on curves, and hope to obtain an approximating family using standard techniques. Unfortunately, the situation is not so simple; as the example in Figure 2 shows, writhe is not continuous in any $`C^k`$ norm on curves. Thus, our proof depends explicitly on the hypothesis that the limit curve is polygonal; it cannot be easily extended to the case where the limit curve is merely piecewise $`C^2`$.

Figure 2. The family of almost-planar curves on the left converge in any $`C^k`$ norm to the planar figure eight curve on the right. However, the writhe of the curves on the left approaches one, while the writhe of the planar figure eight is zero. This shows that writhe is not continuous in any $`C^k`$ norm on curves.

To prepare for the proof, we establish some notation for polygonal curves. Let $`C(t)`$ be a polygonal curve with corners at cyclically ordered parameter values $`t_0<t_1<\mathrm{}<t_n=t_0`$. We let $`T(t)`$ denote the unit tangent to $`C`$, and set up the convention that $`T(t_i)`$ will be the tangent vector leaving $`C(t_i)`$.

We now construct a family of smooth curves approximating our polygonal curve.

###### Proposition 6.

Given an embedded polygonal curve $`C`$ with corners at $`t_0,\mathrm{},t_{n1},t_n=t_0`$, there exists a family of smooth curves $`C_i`$ converging pointwise to $`C`$ with

1. $`C_i=C`$ outside a neighborhood of each corner point $`C(t_j)`$ of radius $`1/i`$.
2. Near each corner, the tangent vectors of $`C_i`$ interpolate between $`T(t_{j1})`$ and $`T(t_j)`$.
3. $`\mathrm{Wr}(C_i)\mathrm{Wr}(C)`$.

###### Proof.

It is easy to construct a family of $`C_iC`$ obeying conditions (1) and (2) by rounding off each corner of $`C`$. We claim that this can be done in such a way that the writhe integrand has a uniform upper bound on all the $`C_i`$. Since condition (1) implies that the $`C_iC`$ pointwise in the $`𝒞^1`$ norm, the bounded convergence theorem \[15, p.81\] will then yield condition (3).

Since any pair of adjacent edges is planar, we can choose the $`C_i`$ so that the region of each $`C_i`$ approximating a pair of adjacent edges is also planar. This means that for some universal $`ϵ`$, the writhe integrand of each $`C_i`$ vanishes in an $`ϵ`$-neighborhood of the diagonal of $`C_i\times C_i`$.

Since $`C`$ has no self-intersections and the angle at each corner of $`C_i`$ is positive, the distance between any pair of non-adjacent edges of $`C`$ is bounded below by some constant. Since the $`C_i`$ converge to $`C`$ pointwise, we may assume the same for the portions of the $`C_i`$ approximating any pair of disjoint edges. Throwing away finitely many of the $`C_i`$ if necessary, this means that for any $`\delta >0`$, there exists a universal lower bound (depending on $`\delta `$) on the distance between any pair of points in $`C_i\times C_i`$ outside an $`\delta `$-neighborhood of the diagonal.

But for any pair of points on $`C_i`$, the writhe integrand is bounded above by the inverse square of the distance between them. Thus, our lower bound on self-distances yields a universal upper bound on the writhe integrand for $`C`$ and all the $`C_i`$ outside a $`\delta `$-neighborhood of the diagonal. Choosing $`\delta <ϵ`$, this completes the proof of the proposition. ∎

## 6. Extending Fuller’s formula to polygonal curves: II

We now state our extension of Fuller’s theorem. Our formula will apply to the following situation (c.f. Section 3): Suppose that $`X_0`$ and $`X_1`$ are simple closed space curves, with $`X_0`$ of class $`𝒞^2`$ and $`X_1`$ polygonal, with regular parametrization (that is, parametrized so that $`X_0^{}`$ and $`X_1^{}`$ never vanish where they are defined), and unit tangent vectors $`T_0`$ and $`T_1`$.

Let $`F:S^1\times [0,1]𝐑^3`$ be a $`C^0`$ deformation of $`X_0`$ into $`X_1`$, where $`F(t,\lambda )=X_\lambda (t)`$, so that the $`X_\lambda `$ are simple curves of class $`C^1`$ for $`\lambda [0,1)`$, with unit tangent vectors $`T_\lambda (t)`$ continuous in $`(t,\lambda )`$. As above, we take the corners of $`X_1`$ to be at parameter values $`t_0,t_1,\mathrm{},t_n=t_0`$. We let $`T_1`$ denote the unit tangent vector to $`X_1`$, and let $`T_1(t_i)`$ be the tangent vector leaving $`X_1(t_i)`$.

###### Theorem 7.

If each corner angle of $`X_1`$ is strictly greater than $`\pi /2`$, and each $`T_1(t)`$ and $`T_\lambda (t)`$ are at an angle less than $`\pi /2`$, then

$$\mathrm{Wr}(X_1)\mathrm{Wr}(X_0)=\frac{1}{2\pi }\underset{i=1}{\overset{n}{}}\mathrm{Area}R(T_0(t_i),T_0(t_{i+1}),T_1(t_i))+\mathrm{Area}\mathrm{}T_0(t_i)T_1(t_{i1})T_1(t_i),$$

where $`R(T_0(t_i),T_0(t_{i+1}),T_1(t_i))`$ is the spherical region bounded by geodesics from $`T_1(t_i)`$ to $`T_0(t_i)`$ and $`T_0(t_{i+1})`$ and the portion of $`T_0`$ between $`t_i`$ and $`t_{i+1}`$, $`\mathrm{}T_0(t_i)T_1(t_{i1})T_1(t_i)`$ is the spherical triangle with these three vertices, and $`\mathrm{Area}`$ represents oriented area on $`S^2`$.

###### Proof.

Construct a sequence of smooth curves $`C_jX_1`$ using Proposition 6. For large enough $`j`$, each of these curves can be homotoped to $`X_1`$ through a family of simple $`𝒞^1`$ curves with a continuous family of tangent vectors, as in the setup for the statement of this theorem above.

Joining these homotopies to the homotopy from $`X_1`$ to $`X_0`$ assumed by our hypotheses generates a family of (non-smooth) homotopies from the $`X_0`$ to each of the $`C_j`$. We wish to smooth each of these to obtain homotopies from $`X_0`$ to $`C_j`$ which obey the conditions of Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula (Theorem 4).

We first prove that the tangent vectors of each of the intermediate curves in each homotopy from $`X_0`$ to $`C_j`$ are never antipodal to the corresponding tangent vectors $`T_j`$ of $`C_j`$. By hypothesis, for each $`t`$ and $`\lambda `$, $`\mathrm{}T_\lambda (t),T_1(t)<\pi /2`$. On the other hand, since the difference between the tangent vectors to $`X_1`$ at any corner is less than $`\pi /2`$, for large enough $`j`$, $`\mathrm{}T_1(t),T_j(t)<\pi /2`$. Putting these equations together, we see that $`\mathrm{}T_\lambda (t),T_j(t)<\pi ,`$ and so these vectors are never antipodal.

It is easy to smooth the combined homotopy from $`X_0`$ to $`C_j`$ so that each of the intermediate curves is of class $`𝒞^1`$ while preserving this condition. Since the smoothed homotopy satisfies the hypotheses of Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula (Theorem 4), Proposition 5 tells us that the difference between $`\mathrm{Wr}(X_0)`$ and $`\mathrm{Wr}(C_j)`$ is the spherical area of the ribbon joining $`T_0`$ and $`T_j`$.

For each $`i`$, the contribution to the spherical area from the straight part of $`C_j`$ between $`t_i`$ and $`t_{i+1}`$ comes from the ribbon between $`T_1(t_i)`$ and the portion of $`T_0`$ with $`t(t_i+1/j,t_{i+1}1/j)`$. As $`j\mathrm{}`$, this area converges to the area of the ribbon between the portion of $`T_0`$ with $`t(t_i,t_{i+1})`$ and $`T_1(t_i)`$. This is the first term in our sum above.

At each vertex $`t_i`$ of $`X_1`$, the contribution to our spherical area from the curved part of $`C_j`$ comes from the ribbon between the great circle arc connecting $`T_1(t_{i1})`$ and $`T_1(t_i)`$ and a portion of $`T_0`$ of parameter length $`2/j`$. As $`j\mathrm{}`$, the area of this ribbon converges to the area of the spherical triangle with vertices $`T_0(t_i)`$, $`T_1(t_i)`$, $`T_1(t_{i1})`$. This is the second term in our sum above. Figure 3 shows both these terms on the unit sphere.

Figure 3. This figure shows the two types of regions in the sum in the statement of Theorem 7. The top (dotted) curve shows the great circle arcs joining the tangent vectors $`T(t_i)`$ of the polygonal curve $`X_1`$. The bottom curve shows the continuous curve of unit tangents $`T_0`$ to the smooth curve $`X_0`$. The light gray regions show the first terms in the sum, while the dark gray spherical triangles show the second terms.

We have shown that the right-hand side of the statement of the Theorem is equal to the limit $`lim_j\mathrm{}(\mathrm{Wr}(C_j)\mathrm{Wr}(X_0))`$. However, by Proposition 6, $`lim_j\mathrm{}\mathrm{Wr}(C_j)=\mathrm{Wr}(X_1)`$. Thus

(9)
$$\underset{j\mathrm{}}{lim}\mathrm{Wr}(C_j)\mathrm{Wr}(X_0)=\mathrm{Wr}(X_1)\mathrm{Wr}(X_0),$$

which is the left-hand side in the statement of the Theorem. This completes the proof. ∎

We now make a surprising observation: Since the tantrices of the $`C_j`$ differ as curves on $`S^2`$ only in parametrization, the area between each of these curves and the tantrix of $`X_0`$ is constant. Thus, by Fuller’s formula, each $`C_j`$ has the same writhe! And since (by Proposition 6) these writhing numbers converge to the writhe of $`X_1`$, each $`\mathrm{Wr}(C_j)`$ is equal to $`\mathrm{Wr}(X_1)`$ as well! So we have the following corollary:

###### Corollary 8.

If $`C_n`$ is a polygonal curve, and $`C`$ is a smooth curve obtained by rounding off the corners of $`C_n`$ under the conditions of Proposition 6, then

(10)
$$\mathrm{Wr}(C_n)=\mathrm{Wr}(C).$$

## 7. Bounding the $`\mathrm{\Delta }\mathrm{Wr}`$ Formula

We now prove our main theorem by finding asymptotic bounds for Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula. Our theorem deals with the following situation: Assume that $`C(t)`$ is a simple closed curve of class $`𝒞^4`$, parametrized so that $`|C^{}(t)|1`$. (Given any initial parametrization, this can be accomplished by rescaling.) Further, assume we have upper bounds $`B_1`$, …, $`B_4`$ on the norms of the first four derivatives of $`C`$.

Let $`C_n(t)`$ be any $`n`$-edge polygonal curve inscribed in $`C`$. We assume that the maximum edge length of $`C`$ is bounded by $`x`$.

###### Theorem 9.

If the ribbon formed by joining $`C_n(t)`$ to $`C(t)`$ for every $`t`$ is embedded, and $`1/x>5B_2`$,

(11)
$$|\mathrm{Wr}(C)\mathrm{Wr}(C_n)|<\alpha nx^3+nO(x^4).$$

where $`\alpha `$ is a numerical constant less than $`B_2(5B_2^2+B_3)`$.

We make a few comments on this theorem before diving into the proof. First, we observe that if the lengths of the edges of $`C_n`$ are all of the same order of magnitude, the difference between the writhe of $`C`$ and the writhe of $`C_n`$ is of order $`1/n^2`$.

Next, we observe that the form of our theorem was chosen to be of maximal use in applications. In particular, we did not require that $`C`$ be parametrized by arclength and state our bounds in terms of curvature and torsion because in practice it is very difficult to obtain an arc-length parametrization of a given curve, while it is comparatively easy to obtain values for the derivative bounds given above.

Last, we discuss the role of the additional hypotheses in the statement above; that the ribbon between $`C`$ and $`C_n`$ be embedded and that $`1/x`$ be greater than $`5B_2`$. Both are intended to exert enough control over the approximation to guarantee the existence of a homotopy from $`C`$ to $`C_n`$ obeying the requirements of Theorem 7.

We can guarantee that $`C_n`$ satisfies the first hypothesis by proving that $`C_n`$ lies in an embedded tubular neighborhood of $`C`$. Since $`C`$ is of class $`𝒞^4`$, and has no self-intersections, such a neighborhood is guaranteed to exist: for a discussion of how to compute the radius of this tube (which is known as the thickness of $`C`$), see the literature on ropelength of knots (*e.g.* ).

###### Proof.

We begin by reparametrizing our curve by arclength. This forces us to recompute our bounds for the derivatives of $`C(t)`$ (a standard computation), arriving at

(12)
$$|C^{}(s)|=1,|C^{\prime \prime }(s)|<K:=2B_2,|C^{\prime \prime \prime }(s)|<T:=2B_3+10B_2^2,$$

while $`C^{(4)}(s)`$ is again bounded above. To remind ourselves of the connection between these bounds and the curvature and torsion of our curve, we will refer to the bound for the second derivative as $`K`$, and the bound for the third derivative as $`T`$. Further, we note that the curvature $`\kappa (s)`$ of our curve is bounded above by $`K`$, and that our hypotheses imply that $`1/x>(5/2)K`$.

We also establish the convention that the corners of $`C_n`$ are at parameter values cyclically ordered as $`s_0,\mathrm{},s_{n1},s_n=s_0`$.

By smoothing the linear interpolation between $`C`$ and $`C_n`$, we can construct a homotopy between $`C`$ and $`C_n`$ according to the conditions of Theorem 7 as long as:

1. the ribbon joining $`C`$ to $`C_n`$ is embedded,
2. the angle at each corner of $`C_n`$ is at least $`\pi /2`$,
3. the angle between $`T(s)`$ and $`T_n(s)`$ is at most $`\pi /2`$ for any $`s`$.

Borrowing from Lemma 12 (below), we see that our assumption that $`1/x>(5/2)K`$ is enough to bound the angle in (3) by $`0.20402<\pi /4`$. At any corner $`s_i`$, the same Lemma implies that the corner angle is the supplement of at most twice $`0.20402`$, so this is enough to ensure that condition (2) holds as well.

Theorem 7 now tells us that

(13)
$$\begin{array}{c}|\mathrm{Wr}(C)\mathrm{Wr}(C_n)|\frac{1}{2\pi }\underset{i=1}{\overset{n}{}}\left|\mathrm{Area}R(T\left(s_i\right),T\left(s_{i+1}\right),T_n\left(s_i\right))\right|\hfill \\ \hfill +\left|\mathrm{Area}\mathrm{}T\left(s_i\right)T_n\left(s_{i1}\right)T_n\left(s_i\right)\right|,\end{array}$$

where the first term is the area of the spherical region bounded by the geodesics from $`T_n(s_i)`$ to $`T(s_i)`$ and $`T(s_{i+1})`$ and the portion of $`T`$ between $`s_i`$ and $`s_{i+1}`$, and the second term is the area of the spherical triangle. Our job now is to estimate the areas of these regions. To do so, we first invoke Taylor’s Theorem, in the form commonly used in numerical analysis (c.f. , Thm.1.4):

###### Theorem 10.

(Taylor’s Theorem) Suppose $`C(s)`$ is a curve of class $`𝒞^4`$, with fourth derivative bounded by $`B_4^{}`$. Then (choosing coordinates so that $`C(0)`$ is at the origin),

(14)
$$C(s)=sC^{}(0)+\frac{s^2}{2}C^{\prime \prime }(0)+\frac{s^3}{6}C^{\prime \prime \prime }(0)+R_4(s),$$

where $`|R_4(s)|<s^4B_4^{}`$.

We will use this expression for $`C(s)`$ frequently in our work below.

###### Lemma 11.

For any $`s`$, we have

(15)
$$|s|C(s)||<\frac{K^2}{24}|s^3|+\frac{1}{120}|s^5|,\text{ and }|C(s)||s|.$$

Further, for any edge of $`C_n`$, the difference $`|s_{i+1}s_i|`$ is at most $`1.01x`$.

###### Proof.

We assume without loss of generality that $`s`$ is positive. By Schur’s lemma (), since the curvature of $`C`$ is bounded above by $`K`$, $`|C(x)|`$ is at least the length of a chord across an arc of length $`s`$ on a circle of radius $`1/K`$, or $`(2/K)\mathrm{sin}(K/2)s`$. This means that we have

$$\frac{2}{K}\mathrm{sin}\frac{K}{2}s=s\frac{K^2}{24}s^3+R_5(s),$$

where $`R_5(s)`$ is the term of order $`s^5`$ which comes from the usual Taylor expansion of $`\mathrm{sin}s`$. In particular,

$`|s|C(s)||`$ $`<`$ $`s{\displaystyle \frac{2}{K}}\mathrm{sin}{\displaystyle \frac{K}{2}}s`$
$`<`$ $`{\displaystyle \frac{K^2}{24}}s^3R_5(s),`$
where $`R_5(s)<\frac{1}{120}s^5`$. The upper bound on $`|C(s)|`$ comes from the fact that $`C`$ is unit-speed.

The second statement is another Schur’s lemma calculation; this time invoking our hypothesis that $`x>(5/2)K`$ and observing that $`1.01\mathrm{sin}y>y`$ for $`y`$ between $`0`$ and $`1/5`$. ∎

We will also need an upper bound on the angle between $`T(s)`$ and $`C_n(s)`$.

###### Lemma 12.

The angle between the tangent vector $`T(s)`$ and the corresponding tangent vector $`T_n(s)`$ to $`C_n`$ is bounded above by

(16)
$$\mathrm{}T(s)T_n(s)<0.51005Kx.$$

###### Proof.

Assume that $`s`$ is between $`s_i`$ and $`s_{i+1}`$. Then

(17)
$$\mathrm{sin}\mathrm{}T(s)T_n(s)=\frac{|[C(s_{i+1})C(s_i)]\times T(s)|}{|C(s_{i+1})C(s_i)|}.$$

But we have

$$C(s_{i+1})C(s_i)=_{s_i}^{s_{i+1}}T(t)𝑑t,$$

and for any $`t`$, we have

$$T(t)=T(s)+_s^tT^{}(u)𝑑u.$$

This means that

(18) $`[C(s_{i+1})C(s_i)]\times T(s)`$ $`=`$ $`{\displaystyle _{s_i}^{s_{i+1}}}T(t)\times T(s)𝑑t`$
(19) $`=`$ $`{\displaystyle _{s_i}^{s_{i+1}}}{\displaystyle _s^t}T^{}(u)\times T(s)𝑑u𝑑t.`$
Since $`|T^{}(u)\times T(s)||T^{}(u)||T(s)|\kappa (u)<K`$, and $`s`$ is between $`s_i`$ and $`s_{i+1}`$, a small computation reveals that this integral is bounded by $`\frac{K}{2}(s_{i+1}s_i)^2`$.

Since the length $`|C(s_{i+1})C(s_i)|`$ is bounded below by $`(1/1.01)(s_{i+1}s_i)`$ by Lemma 11, we get

(20)
$$\mathrm{sin}\mathrm{}T(t)T_n(t)<\frac{1.01}{2}Kx.$$

Since $`1/x>(5/2)K`$, this is always bounded above by $`1.01/5`$, and so

(21)
$$\mathrm{}T(t)T_n(t)<\frac{1.01^2}{2}Kx.$$

We are now ready to embark on the main work of the proof: estimating the areas in Equation 13. We begin with the first term: the area bounded by the portion of $`T(s)`$ between $`s_i`$ and $`s_{i+1}`$, together with the great circle arcs joining $`T(s_i)`$ and $`T(s_{i+1})`$ to $`T_n(s_i)`$. Without loss of generality, we may assume that $`i=0`$, that $`s_0=0`$, and that $`C(0)=\mathrm{𝟎}`$, and apply the Taylor expansion of Equation 14 to $`C`$ at $`0`$. Our strategy is to prove that this region is contained in a neighborhood of the great circle arc joining $`T(0)`$ and $`T(s_1)`$. Suppose $`s`$ is between $`0`$ and $`s_1`$. We want to bound the height of $`T(s)`$ above the $`T(0)`$, $`T(s_1)`$ plane, or

(22)
$$h(s):=\frac{C^{}(s)C^{}(0)\times C^{}(s_1)}{|C^{}(0)\times C^{}(s_1)|}.$$

First, we have

$`C^{}(s_1)`$ $`=`$ $`C^{}(0)+s_1C^{\prime \prime }(0)+{\displaystyle \frac{s_1^2}{2}}C^{\prime \prime \prime }(0)+R_3(s_1).`$
$`C^{}(s)`$ $`=`$ $`C^{}(0)+sC^{\prime \prime }(0)+{\displaystyle \frac{s^2}{2}}C^{\prime \prime \prime }(0)+R_3(s).`$
Using the triple product identities, we can rewrite $`h(s)`$ in terms of the inner product of $`C^{}(0)`$ and the cross product of these vectors. For the triple product, we get

(23)
$$\left[\frac{s_1^2s}{2}\frac{s^2s_1}{2}\right]C^{}(0)C^{\prime \prime }(0)\times C^{\prime \prime \prime }(0)+C^{}(0)\left[R_3(s_1)\times C^{}(s)+C^{}(s_1)\times R_3(s)\right].$$

Expanding the last term, we see that is the sum of a term of order $`s_1s^3`$ and a term of order $`ss_1^3`$. Thus, to leading order, the norm of the entire triple product is bounded above by

(24)
$$|h(s)|<H:=\frac{KT}{2|C^{}(0)\times C^{}(s_1)|}s_1^3+O(s_1^4),$$

since $`s[0,s_1]`$. We now consider the height of $`T_n(0)`$ above the $`T(0)`$, $`T(s_1)`$ plane. Since $`T_n(0)`$ is the normalization of $`C(s_1)C(0)=C(s_1)`$, this height is given by

(25)
$$\frac{C(s_1)}{|C(s_1)|}\frac{C^{}(0)\times C^{}(s_1)}{|C^{}(0)\times C^{}(s_1)|}.$$

As before, we get

(26)
$$C^{}(0)\times C^{}(s_1)=s_1C^{}(0)\times C^{\prime \prime }(0)+\frac{s_1^2}{2}C^{}(0)\times C^{\prime \prime \prime }(0)+O(s_1^3).$$

Taking the dot product with the Taylor expansion of $`C(s_1)`$, we get only terms of order $`O(s_1^4)`$ and higher. Thus, to leading order, this region is contained in a rectangle based on the great circle arc joining $`C^{}(0)`$ and $`C^{}(s_1)`$ of height $`H`$. We now estimate the area of this rectangle.

First, we note that the length of the great circle joining $`C^{}(0)`$ and $`C^{}(s_1)`$ is given by the angle $`\theta `$ between $`C^{}(0)`$ and $`C^{}(s_1)`$. Since $`s_1<1.01x`$ by Lemma 11, this length is bounded above by $`1.01Kx`$, which is less than $`0.404`$ by our hypotheses on $`x`$. Since $`H`$ is small compared to $`s`$, we may assume that the entire rectangle is contained within a spherical disk of radius $`0.5`$.

We project the rectangle to the plane by central projection: this map is increasing on lengths and areas, and increases length by at most a factor of $`1.01`$. The area of the rectangle in the plane is overestimated by the product $`1.01\theta H`$. On the other hand, we have $`|C^{}(0)\times C^{}(s_1)|=\mathrm{sin}\theta `$. And for $`\theta <0.404`$, $`1.02\mathrm{sin}\theta >\theta `$. Keeping track of the various constants involved, and using the fact that $`s_1<1.01x`$ again, the area of this spherical region is overestimated by

(27)
$$\mathrm{Area}R(T(s_i),T(s_{i+1}),T_n(s_i))<KTx^3+O(x^4),$$

We now turn to the second term in the Equation 13: the area of the spherical triangle bounded by $`T(s_i)`$, $`T_n(s_{i1})`$ and $`T_n(s_i)`$. Without loss of generality we assume that $`i=1`$, that $`s_1=0`$, and that $`C(0)=\mathrm{𝟎}`$, and we expand $`C`$ around $`0`$ using Equation 14. We wish to compute

(28)
$$\mathrm{Area}\mathrm{}(\frac{C(s_0)}{|C(s_0)|},\frac{C(s_2)}{|C(s_2)|},C^{}(0))=\left|\left(\frac{C(s_0)}{|C(s_0)|}C^{}(0)\right)\times \left(\frac{C(s_2)}{|C(s_2)|}C^{}(0)\right)\right|.$$

If we factor out $`1/|C(s_0)||C(s_2)|`$, we are left with the norm of the cross product of two terms:

$`C(s_0)|C(s_0)|C^{}(0)`$ $`=`$ $`(s_0|C(s_0)|)C^{}(0)+{\displaystyle \frac{s_0^2}{2}}C^{\prime \prime }(0)+{\displaystyle \frac{s_0^3}{6}}C^{\prime \prime \prime }(0)+R_4(s_0)`$
$`C(s_2)|C(s_2)|C^{}(0)`$ $`=`$ $`(s_2|C(s_2)|)C^{}(0)+{\displaystyle \frac{s_2^2}{2}}C^{\prime \prime }(0)+{\displaystyle \frac{s_2^3}{6}}C^{\prime \prime \prime }(0)+R_4(s_2).`$
Using Lemma 11, we see that $`|s|C(s)||<(K^2/24)s^3+O(s^5)`$, and we see that the leading term of this expression contains fifth powers of of $`s_0`$ and $`s_2`$, and is bounded by:

(29)
$$s_0^2s_2^2\left(\frac{K^3}{48}+\frac{KT}{12}\right)(s_0+s_2)$$

However, we must still divide by $`|C(s_0)||C(s_2)|`$. By Lemma 11, we see that the ratios $`s_0/|C(s_0)|`$ and $`s_2/|C(s_2)|`$ are bounded above by $`1.01`$. Thus, using the same Lemma to conclude that $`s_2`$ and $`s_0`$ are less than $`1.01x`$, and making a central projection argument as before, we are left with

(30)
$$\mathrm{Area}\mathrm{}(T_n(s_i),T_n(s_{i1}),T(s_i))<\frac{K^3+KT}{3}x^3+O(x^4).$$

Summing over $`i`$, and dividing by $`2\pi `$, then writing $`K`$ and $`T`$ in terms of $`B_2`$ and $`B_3`$, we obtain the statement of the theorem. Note that we have overestimated the numerical constants to simplify the resulting formula. ∎

If a curve has a small region of high curvature, and larger regions of low curvature, it may be desirable to approximate the curve more carefully in the regions of high curvature in order to save time in the computation of writhe. Since our error bound is additive along the curve, these methods are well suited to this case. We have

###### Corollary 13.

Suppose $`C`$ is a $`𝒞^4`$ curve and $`C_n`$ is a curve inscribed in $`C`$ so that $`C`$ and $`C_n`$ obey the hypotheses of Theorem 9.

If $`C`$ and $`C_n`$ are divided into regions $`R_i`$, each containing $`n_i`$ edges which are bounded above in length by $`x_i`$, and so that the derivatives of $`C`$ are bounded by $`B_{1i}`$, …, $`B_{4i}`$ and $`1/x_i>5B_{2i}`$, then

$$|\mathrm{Wr}(C)\mathrm{Wr}(C_p)|<\underset{i}{}\alpha _in_ix_i^3+n_iO(x_i^4).$$

where each $`\alpha _i`$ is a numerical constant less than $`B_{2i}(5B_{2i}^2+B_{3i})`$.

We make one more observation:

###### Proposition 14.

Let $`C`$ be a simple, closed space curve of class $`𝒞^2`$, and $`C_p`$ be a polygonal approximating curve as in Theorem 9 or Corollary 13.

If the arc joining the endpoints of a sequence of $`n`$ edges of $`C_p`$ is planar, then the $`n2`$ edges interior to this region contribute nothing to the error bound in the Theorem.

In particular, this means that the derivative bounds in both statements can be taken to be bounds on the derivatives of the non-planar regions of the curve $`C`$.

###### Proof.

On these edges, the tantrix of the smooth curve and the polygonal curve parametrize the same great circle arc on $`S^2`$. Thus, the ribbon between these curves has zero area. ∎

## 8. Example Computations

We are now prepared to test Theorem 1 by computing the writhing numbers of various polygonal approximations of a smooth curve, and comparing the results to the exact writhe of the smooth curve. To control the numerical error introduced in these calculations, all of these computations were performed using an arbitrary-precision implementation of Banchoff’s formula for the writhing number of a polygonal curve. The initial runs were performed with $`45`$ decimal digits of precision. They were checked against runs performed with $`54`$ digits of precision. Since the results agreed, we feel confident that roundoff error does not affect the computations reported on below.

The curve whose writhe we computed is an example of Fuller:

|  |  |
| --- | --- |

Figure 4. This example of Fuller’s “closed helix” is composed of $`3`$ turns of a helix of radius $`1`$ with pitch angle $`0.33`$, with ends joined by a planar curve.

Using Theorem 3, and the Călugăreanu-White formula, it is easy to see that the writhe of this curve is $`3(1\mathrm{sin}0.33)2.0278709`$. After all, the area enclosed by the tantrix of this curve $`C`$ is that of a hemisphere, plus $`3`$ enclosures of a spherical cap of radius $`\pi /20.33`$. Thus the writhe of the curve is equal to $`1\mathrm{sin}0.33mod2`$. To complete the computation, one sets up a frame on the curve, and computes its twist and linking number. (Details for this computation can be found in .)

We now take a series of polygonal approximations to $`C`$, and compare the difference between their writhing numbers and the writhe of $`C`$ to the bounds of Theorem 9.

We begin by finding bounds on the derivatives of $`C`$ and the edge length of our approximations. By Proposition 14, it suffices to find derivative bounds for the helical region of $`C`$. Since the helix has unit radius, both $`B_2`$ and $`B_3`$ can be taken to be one. The curve is parametrized so that $`|C^{}(s)|1`$.

Here are the results of computing writhe with various numbers of edges:

| $`n`$ | $`\mathrm{Wr}(C_n)`$ | $`|\mathrm{Wr}(C_n)\mathrm{Wr}(C)|`$ | $`x`$ | $`\alpha nx^3`$ |
| --- | --- | --- | --- | --- |
| 100 | $`2.00541`$ | $`0.02246`$ | $`0.506`$ | 77.73 |
| 250 | $`2.02434`$ | $`0.00353`$ | $`0.203`$ | 12.55 |
| 500 | $`2.02697`$ | $`0.0009`$ | $`0.101`$ | 3.09 |
| 1000 | $`2.02763`$ | $`0.00024`$ | $`0.051`$ | 0.786 |

It is worth examing a graph of these results.

Figure 5. This graph shows a log-log plot of the actual error in computing the writhing number for one of Fuller’s “closed helices” with various numbers of edges (lower solid line), together with our error bounds (upper dotted line). The fact that the lines are parallel shows that the convergence is of order $`n^2`$, as predicted by Theorem 9.

## 9. Further Directions

In this paper, we have given a set of asymptotic error bounds which allow us to compute the writhe of a closed space curve with defined accuracy by computing the writhe of a polygonal approximation to this curve. The example we computed in Section 8 shows that our bounds are of the right order of magnitude: roughly speaking, the writhe converges quadratically in the number of edges of the approximation. Our work leaves open several directions for further inquiry.

First, it is puzzling that our approximation theorem for curves with corners (Proposition 6) should depend on the hypothesis that the limit curve is polygonal. To sketch an extension of this result, we recall a definition from Chern ():

###### Definition 15.

The tantrix of a piecewise $`C^1`$ curve $`C(s)`$ with positive corner angles is the image of $`T(s)`$ on the unit sphere, together with the great circle arcs joining the pairs of tangent vectors at each corner of the curve.

We note that our Theorem 7 shows that Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula holds for polygonal curves with the definition of tantrix extended as above. We further suspect that:

###### Conjecture 16.

Fuller’s Spherical Area formula (Theorem 3) and Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula (Theorem 4) hold for piecewise $`C^2`$ curves with the extended definition of tantrix given by Definition 15.

The proofs of both of these theorems depend on the Călugăreanu-White formula, which only applies to closed curves. Thus all of our results are restricted to closed curves. This leaves open a much more important problem:

###### Problem 17.

Extend all these theorems (the Călugăreanu-White formula, Fuller’s spherical area formula, and Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula) to open curves.

In particular, extending the results of this paper to open curves would be useful for applications in biology, where the curves of interest are not neccesarily closed. We note that while Fuller’s $`\mathrm{\Delta }\mathrm{Wr}`$ formula makes sense for open curves, computational examples show that it does not give the correct answer: boundary terms must be added to account for the ends of the curves.

## 10. Acknowledgements

I am grateful to Herbert Edelsbrunner, Herman Gluck, Issac Klapper, John Maddocks, Kathleen Rogers, and David Swigon, among many others, for fruitful conversations. This work was supported by both a National Defense Science and Engineering Graduate Fellowship and an NSF Postdoctoral Research Fellowship (DMS-99-02397).


format: text

text:
# Multifractal Detrended Fluctuation Analysis of Nonstationary Time Series

## I Introduction

In recent years the detrended fluctuation analysis (DFA) method  has become a widely-used technique for the determination of (mono-) fractal scaling properties and the detection of long-range correlations in noisy, nonstationary time series . It has successfully been applied to diverse fields such as DNA sequences , heart rate dynamics , neuron spiking , human gait , long-time weather records , cloud structure , geology , ethnology , economics time series , and solid state physics . One reason to employ the DFA method is to avoid spurious detection of correlations that are artifacts of nonstationarities in the time series.

Many records do not exhibit a simple monofractal scaling behavior, which can be accounted for by a single scaling exponent. In some cases, there exist crossover (time-) scales $`s_\times `$ separating regimes with different scaling exponents , e. g. long-range correlations on small scales $`ss_\times `$ and another type of correlations or uncorrelated behavior on larger scales $`ss_\times `$. In other cases, the scaling behavior is more complicated, and different scaling exponents are required for different parts of the series . This occurs, e. g., when the scaling behavior in the first half of the series differs from the scaling behavior in the second half. In even more complicated cases, such different scaling behavior can be observed for many interwoven fractal subsets of the time series. In this case a multitude of scaling exponents is required for a full description of the scaling behavior, and a multifractal analysis must be applied.

In general, two different types of multifractality in time series can be distinguished: (i) Multifractality due to a broad probability density function for the values of the time series. In this case the multifractality cannot be removed by shuffling the series. (ii) Multifractality due to different long-range (time-) correlations of the small and large fluctuations. In this case the probability density function of the values can be a regular distribution with finite moments, e. g. a Gaussian distribution. The corresponding shuffled series will exhibit non-multifractal scaling, since all long-range correlations are destroyed by the shuffling procedure. If both kinds of multifractality are present, the shuffled series will show weaker multifractality than the original series.

The simplest type of multifractal analysis is based upon the standard partition function multifractal formalism, which has been developed for the multifractal characterization of normalized, stationary measures . Unfortunately, this standard formalism does not give correct results for nonstationary time series that are affected by trends or that cannot be normalized. Thus, in the early 1990s an improved multifractal formalism has been developed, the wavelet transform modulus maxima (WTMM) method , which is based on wavelet analysis and involves tracing the maxima lines in the continuous wavelet transform over all scales. Here, we propose an alternative approach based on a generalization of the DFA method. This multifractal DFA (MF-DFA) does not require the modulus maxima procedure, and hence does not involve more effort in programming than the conventional DFA.

The paper is organized as follows: In Section II we describe the MF-DFA method in detail and show that the scaling exponents determined via the MF-DFA method are identical to those obtained by the standard multifractal formalism based on partition functions. In Section III we introduce several multifractal models, where the scaling exponents can be calculated exactly, and compare these analytical results with the numerical results obtained by MF-DFA. In Section IV, we show how the comparison of the MF-DFA results for original series to the MF-DFA results for shuffled series can be used to determine the type of multifractality in the series. In Section V, we compare the results of the MF-DFA with those obtained by the WTMM method for nonstationary series and discuss the performance of both methods for multifractal time series analysis.

## II Multifractal DFA

### A Description of the method

The generalized multifractal DFA (MF-DFA) procedure consists of five steps. The first three steps are essentially identical to the conventional DFA procedure (see e. g. ). Let us suppose that $`x_k`$ is a series of length $`N`$, and that this series is of compact support, i.e. $`x_k=0`$ for an insignificant fraction of the values only.

$``$ Step 1: Determine the “profile”

$$Y(i)\underset{k=1}{\overset{i}{}}\left[x_kx\right],i=1,\mathrm{},N.$$
(1)
Subtraction of the mean $`x`$ is not compulsory, since it would be eliminated by the later detrending in the third step.

$``$ Step 2: Divide the profile $`Y(i)`$ into $`N_s\mathrm{int}(N/s)`$ non-overlapping segments of equal length $`s`$. Since the length $`N`$ of the series is often not a multiple of the considered time scale $`s`$, a short part at the end of the profile may remain. In order not to disregard this part of the series, the same procedure is repeated starting from the opposite end. Thereby, $`2N_s`$ segments are obtained altogether.

$``$ Step 3: Calculate the local trend for each of the $`2N_s`$ segments by a least-square fit of the series. Then determine the variance

$$F^2(s,\nu )\frac{1}{s}\underset{i=1}{\overset{s}{}}\left\{Y[(\nu 1)s+i]y_\nu (i)\right\}^2$$
(2)
for each segment $`\nu `$, $`\nu =1,\mathrm{},N_s`$ and

$$F^2(s,\nu )\frac{1}{s}\underset{i=1}{\overset{s}{}}\left\{Y[N(\nu N_s)s+i]y_\nu (i)\right\}^2$$
(3)
for $`\nu =N_s+1,\mathrm{},2N_s`$. Here, $`y_\nu (i)`$ is the fitting polynomial in segment $`\nu `$. Linear, quadratic, cubic, or higher order polynomials can be used in the fitting procedure (conventionally called DFA1, DFA2, DFA3, $`\mathrm{}`$) . Since the detrending of the time series is done by the subtraction of the polynomial fits from the profile, different order DFA differ in their capability of eliminating trends in the series. In (MF-)DFA$`m`$ \[$`m`$th order (MF-)DFA\] trends of order $`m`$ in the profile (or, equivalently, of order $`m1`$ in the original series) are eliminated. Thus a comparison of the results for different orders of DFA allows one to estimate the type of the polynomial trend in the time series .

$``$ Step 4: Average over all segments to obtain the $`q`$th order fluctuation function

$$F_q(s)\left\{\frac{1}{2N_s}\underset{\nu =1}{\overset{2N_s}{}}\left[F^2(s,\nu )\right]^{q/2}\right\}^{1/q},$$
(4)
where, in general, the index variable $`q`$ can take any real value except zero . For $`q=2`$, the standard DFA procedure is retrieved. We are interested in how the generalized $`q`$ dependent fluctuation functions $`F_q(s)`$ depend on the time scale $`s`$ for different values of $`q`$. Hence, we must repeat steps 2 to 4 for several time scales $`s`$. It is apparent that $`F_q(s)`$ will increase with increasing $`s`$. Of course, $`F_q(s)`$ depends on the DFA order $`m`$. By construction, $`F_q(s)`$ is only defined for $`sm+2`$.

$``$ Step 5: Determine the scaling behavior of the fluctuation functions by analyzing log-log plots $`F_q(s)`$ versus $`s`$ for each value of $`q`$. Several examples of this procedure will be shown in Section III. If the series $`x_i`$ are long-range power-law correlated, $`F_q(s)`$ increases, for large values of $`s`$, as a power-law,

$$F_q(s)s^{h(q)}.$$
(5)
In general, the exponent $`h(q)`$ may depend on $`q`$. For stationary time series, $`h(2)`$ is identical to the well-known Hurst exponent $`H`$ (see e. g. ). Thus, we will call the function $`h(q)`$ generalized Hurst exponent.

For monofractal time series with compact support, $`h(q)`$ is independent of $`q`$, since the scaling behavior of the variances $`F^2(s,\nu )`$ is identical for all segments $`\nu `$, and the averaging procedure in Eq. (4) will give just this identical scaling behavior for all values of $`q`$. Only if small and large fluctuations scale differently, there will be a significant dependence of $`h(q)`$ on $`q`$: If we consider positive values of $`q`$, the segments $`\nu `$ with large variance $`F_s^2(\nu )`$ (i. e. large deviations from the corresponding fit) will dominate the average $`F_q(s)`$. Thus, for positive values of $`q`$, $`h(q)`$ describes the scaling behavior of the segments with large fluctuations. Usually the large fluctuations are characterized by a smaller scaling exponent $`h(q)`$ for multifractal series . On the contrary, for negative values of $`q`$, the segments $`\nu `$ with small variance $`F_s^2(\nu )`$ will dominate the average $`F_q(s)`$. Hence, for negative values of $`q`$, $`h(q)`$ describes the scaling behavior of the segments with small fluctuations, which are usually characterized by a larger scaling exponent.

However, the MF-DFA method can only determine positive generalized Hurst exponents $`h(q)`$, and it already becomes inaccurate for strongly anti-correlated signals when $`h(q)`$ is close to zero. In such cases, a modified (MF-)DFA technique has to be used. The most simple way to analyze such data is to integrate the time series before the MF-DFA procedure. Hence, we replace the single summation in Eq. (1), which is describing the profile from the original data $`x_k`$, by a double summation,

$$\stackrel{~}{Y}(i)\underset{k=1}{\overset{i}{}}\left[Y(k)Y\right].$$
(6)
Following the MF-DFA procedure as described above, we obtain generalized fluctuation functions $`\stackrel{~}{F}_q(s)`$ described by a scaling law as in Eq. (5), but with larger exponents $`\stackrel{~}{h}(q)=h(q)+1`$,

$$\stackrel{~}{F}_q(s)s^{\stackrel{~}{h}(q)}=s^{h(q)+1}.$$
(7)
Thus, the scaling behavior can be accurately determined even for $`h(q)`$ which are smaller than zero (but larger than $`1`$) for some values of $`q`$. We note that $`\stackrel{~}{F}_q(s)/s`$ corresponds to $`F_q(s)`$ in Eq. (5). If we do not subtract the average values in each step of the summation in Eq. (6), this summation leads to quadratic trends in the profile $`\stackrel{~}{Y}(i)`$. In this case we must employ at least the second order MF-DFA to eliminate these artificial trends.

### B Relation to standard multifractal analysis

For stationary, normalized records with compact support the multifractal scaling exponents $`h(q)`$ defined in Eq. (5) are directly related, as shown below, to the scaling exponents $`\tau (q)`$ defined by the standard partition function-based multifractal formalism.

Suppose that the series $`x_k`$ of length $`N`$ is a stationary, normalized sequence. Then the detrending procedure in step 3 of the MF-DFA method is not required, since no trend has to be eliminated. Thus, the DFA can be replaced by the standard Fluctuation Analysis (FA), which is identical to the DFA except for a simplified definition of the variance for each segment $`\nu `$, $`\nu =1,\mathrm{},N_s`$, in step 3 \[see Eq. (2)\]:

$$F_{\mathrm{FA}}^2(s,\nu )[Y(\nu s)Y((\nu 1)s)]^2.$$
(8)
Inserting this simplified definition into Eq. (4) and using Eq. (5), we obtain

$$\left\{\frac{1}{2N_s}\underset{\nu =1}{\overset{2N_s}{}}|Y(\nu s)Y((\nu 1)s)|^q\right\}^{1/q}s^{h(q)}.$$
(9)
For simplicity we can assume that the length $`N`$ of the series is an integer multiple of the scale $`s`$, obtaining $`N_s=N/s`$ and therefore

$$\underset{\nu =1}{\overset{N/s}{}}|Y(\nu s)Y((\nu 1)s)|^qs^{qh(q)1}.$$
(10)
This already corresponds to the multifractal formalism used e. g. in . In fact, a hierarchy of exponents $`H_q`$ similar to our $`h(q)`$ has been introduced based on Eq. (10) in .

In order to relate also to the standard textbook box counting formalism , we employ the definition of the profile in Eq. (1). It is evident that the term $`Y(\nu s)Y((\nu 1)s)`$ in Eq. (10) is identical to the sum of the numbers $`x_k`$ within each segment $`\nu `$ of size $`s`$. This sum is known as the box probability $`p_s(\nu )`$ in the standard multifractal formalism for normalized series $`x_k`$,

$$p_s(\nu )\underset{k=(\nu 1)s+1}{\overset{\nu s}{}}x_k=Y(\nu s)Y((\nu 1)s).$$
(11)
The scaling exponent $`\tau (q)`$ is usually defined via the partition function $`Z_q(s)`$,

$$Z_q(s)\underset{\nu =1}{\overset{N/s}{}}|p_s(\nu )|^qs^{\tau (q)},$$
(12)
where $`q`$ is a real parameter as in the MF-DFA above. Sometimes $`\tau (q)`$ is defined with opposite sign (see e. g. ).

Using Eq. (11) we see that Eq. (12) is identical to Eq. (10), and obtain analytically the relation between the two sets of multifractal scaling exponents,

$$\tau (q)=qh(q)1.$$
(13)
Thus, we have shown that $`h(q)`$ defined in Eq. (5) for the MF-DFA is directly related to the classical multifractal scaling exponents $`\tau (q)`$. Note that $`h(q)`$ is different from the generalized multifractal dimensions

$$D(q)\frac{\tau (q)}{q1}=\frac{qh(q)1}{q1},$$
(14)
that are used instead of $`\tau (q)`$ in some papers. While $`h(q)`$ is independent of $`q`$ for a monofractal time series with compact support, $`D(q)`$ depends on $`q`$ in that case. Our assumption of compact support of the series $`x_k`$ can be directly observed in Eq. (14), since the fractal dimension of the support is $`D(0)\tau (0)=1`$.

Another way to characterize a multifractal series is the singularity spectrum $`f(\alpha )`$, that is related to $`\tau (q)`$ via a Legendre transform ,

$$\alpha =\tau ^{}(q)\mathrm{and}f(\alpha )=q\alpha \tau (q).$$
(15)
Here, $`\alpha `$ is the singularity strength or Hölder exponent, while $`f(\alpha )`$ denotes the dimension of the subset of the series that is characterized by $`\alpha `$. Using Eq. (13), we can directly relate $`\alpha `$ and $`f(\alpha )`$ to $`h(q)`$,

$$\alpha =h(q)+qh^{}(q)\mathrm{and}f(\alpha )=q[\alpha h(q)]+1.$$
(16)
## III Four Illustrative examples

### A Example 1: monofractal uncorrelated and long-range correlated series

As a first example we apply the MF-DFA method to monofractal series with compact support, for which the generalized Hurst exponent $`h(q)`$ is expected to be independent of $`q`$,

$$h(q)=H\mathrm{and}\tau (q)=qH1.$$
(17)
Such series have been discussed in the context of conventional DFA in several studies before, see e.g. . Long-range correlated random numbers are usually generated by the Fourier transform method, see e. g. . Using this method we can generate long-range anti-correlated ($`0<H<0.5`$), uncorrelated ($`H=0.5`$), or (positively) long-range correlated ($`0.5<H<1`$) series. The latter are characterized by a power-law decay of the autocorrelation function $`C(s)x_kx_{k+s}s^\gamma `$ for large scales $`s`$ with $`\gamma =22H`$ if the series is stationary. Alternatively, all stationary long-range correlated series can be characterized by the power-law decay of their power spectra, $`S(f)f^\beta `$ with frequency $`f`$ and $`\beta =2H1`$. Note that $`H`$ corresponds to the Hurst exponent of the integrated time series here.

Figure 1 shows the generalized fluctuation functions $`F_q(s)`$ for all three types of monofractal series ($`H=0.75,0.5,0.25`$) and several $`q`$ values. On large scales $`s`$, we observe the expected power-law scaling behavior according to Eq. (5), which corresponds to straight lines in the log-log plot. In Fig. 1(d), the scaling exponents $`h(q)`$ determined from the slopes of these straight lines are shown versus $`q`$. Although a slight $`q`$ dependence is observable, the values of $`h(q)`$ are always very close to the $`H`$ of the generated series that has been analyzed. The degree of the $`q`$ dependence observed for this monofractal series allows to estimate the usual fluctuation of $`h(q)`$ to be expected for monofractal series in general.

Next, we analyze multifractal series for which $`\tau (q)`$ can be calculated exactly, and compare the numerical results with the expected scaling behavior.

### B Example 2: binomial multifractal series

In the binomial multifractal model , a series of $`N=2^{n_{\mathrm{max}}}`$ numbers $`k`$ with $`k=1,\mathrm{},N`$ is defined by

$$x_k=a^{n(k1)}(1a)^{n_{\mathrm{max}}n(k1)},$$
(18)
where $`0.5<a<1`$ is a parameter and $`n(k)`$ is the number of digits equal to 1 in the binary representation of the index $`k`$, e. g. $`n(13)=3`$, since 13 corresponds to binary 1101.

The scaling exponents $`\tau (q)`$ can be calculated straightforwardly. According to Eqs. (11) and (18) the box probability $`p_{2s}(\nu )`$ in the $`\nu `$th segment of size $`2s`$ is given by

$`p_{2s}(\nu )=p_s(2\nu 1)+p_s(2\nu )`$
$`=[(1a)/a+1]p_s(2\nu )=p_s(2\nu )/a.`$
Thus, according to Eqs. (12) and (18),

$`Z_q(s)={\displaystyle \underset{\nu =1}{\overset{N/s}{}}}[p_s(\nu )]^q={\displaystyle \underset{\nu =1}{\overset{N/2s}{}}}[p_s(2\nu 1)]^q+[p_s(2\nu )]^q`$
$`=\left[{\displaystyle \frac{(1a)^q}{a^q}}+1\right]{\displaystyle \underset{\nu =1}{\overset{N/2s}{}}}[p_s(2\nu )]^q`$
$`=[(1a)^q+a^q]{\displaystyle \underset{\nu =1}{\overset{N/2s}{}}}[p_{2s}(\nu )]^q=[(1a)^q+a^q]Z_q(2s)`$
and according to Eqs. (12) and (13),

$`\tau (q)={\displaystyle \frac{\mathrm{ln}[a^q+(1a)^q]}{\mathrm{ln}(2)}},`$ (19)
$`h(q)={\displaystyle \frac{1}{q}}{\displaystyle \frac{\mathrm{ln}[a^q+(1a)^q]}{q\mathrm{ln}(2)}}.`$ (20)
Note that $`\tau (0)=1`$ as required. There is a strong non-linear dependence of $`\tau (q)`$ upon $`q`$, indicating multifractality. The same information is comprised in the $`q`$ dependence of $`h(q)`$. The asymptotic values are $`h(q)\mathrm{ln}(a)/\mathrm{ln}(2)`$ for $`q+\mathrm{}`$ and $`h(q)\mathrm{ln}(1a)/\mathrm{ln}(2)`$ for $`q\mathrm{}`$. They correspond to the scaling behavior of the largest and weakest fluctuations, respectively. Note that $`h(q)`$ becomes independent of $`q`$ in the asymptotic limit, while $`\tau (q)`$ approaches linear $`q`$ dependences.

Figure 2 shows the MF-DFA fluctuation functions $`F_q(s)`$ for the binomial multifractal model with $`a=0.75`$. The results for MF-DFA1 and MF-DFA4 are compared in parts (a) and (b). Fig. 2(c) shows the corresponding slopes $`h(q)`$ for three values of $`a`$ together with the exact results obtained from Eq. (20). The numerical results are in good agreement with Eq. (20), showing that the MF-DFA correctly detects the multifractal scaling exponents. Figures 2(d) and (e) show the corresponding exponents $`\tau (q)=qh(q)1`$ \[see Eq. (13)\] and the corresponding $`f(\alpha )`$ spectrum calculated from $`h(q)`$ using the modified Legendre transform (16). Both are also in good agreement with Eq. (19). We have also checked that the results for the binomial multifractal model remain unchanged if the double summation technique \[see Eq. (6)\] is applied. We obtain slopes $`\stackrel{~}{h}(q)=h(q)+1`$ as expected in Eq. (7). Note that there is no need to use this modification, except if $`h(q)`$ is close to zero or has negative values.

### C Example 3: dyadic random cascade model with log-Poisson distribution

For another independent test of the MF-DFA, we employ an algorithm based on random cascades on wavelet dyadic trees proposed in  (see also ). This algorithm builds a random multifractal series by specifying its discrete wavelet coefficients $`c_{n,m}`$, defined recursively,

$`c_{1,1}=1,c_{n,2m1}=Wc_{n1,m},c_{n,2m}=Wc_{n1,m},`$
where $`n=2,\mathrm{},n_{\mathrm{max}}`$ (with $`N=2^{n_{\mathrm{max}}}`$) and $`m=1,\mathrm{},2^{n2}`$. The values of $`W`$ are taken from a log-Poisson distribution, $`|W|=\mathrm{exp}(P\mathrm{ln}\delta +\gamma )`$, where $`P`$ is Poisson distributed with $`P=\lambda `$. There are three independent parameters, $`\lambda `$, $`\delta `$, and $`\gamma `$. Inverse wavelet transform is applied to create the multifractal random series $`x_k`$ once the wavelet coefficients $`c_{n,m}`$ are known,

$$x_k=\underset{n=1}{\overset{n_{\mathrm{max}}}{}}\underset{m=1}{\overset{2^{n1}}{}}c_{n,m}\psi _{n,m}(k),$$
(21)
where $`\psi _{n,m}(k)`$ is a set of wavelets forming an orthonormal wavelet basis. Here, we employ the Haar wavelets, $`\psi _{n,m}(k)2^{(nn_{\mathrm{max}}1)/2}\psi [2^{nn_{\mathrm{max}}1}km]`$ with $`\psi (x)1`$ for $`0<x0.5`$, $`\psi (x)1`$ for $`0.5<x1`$ and $`\psi (x)0`$ otherwise. For this model the multifractal scaling exponents are given by

$`\tau (q)={\displaystyle \frac{\lambda (1\delta ^q)\gamma q}{\mathrm{ln}2}}1,`$ (22)
$`h(q)=[\lambda (1\delta ^q)\gamma q]/(q\mathrm{ln}2).`$ (23)
Figure 3 shows the MF-DFA fluctuation functions $`F_q(s)`$ for the dyadic random cascade model. The numerically determined slopes $`h(q)`$ for three sets of parameters are compared with the exact results obtained from Eq. (23) and the good agreement shows that the MF-DFA correctly detects the multifractal scaling exponents. Large deviations occur only for very small moments ($`q<10`$), indicating that the range of $`q`$ values should not exceed $`10`$.

### D Example 4: uncorrelated multifractal series with power-law distribution function

The examples discussed in the previous three subsections were based on series involving long-range correlations. In the present example we want to apply the MF-DFA method to an uncorrelated series, that nevertheless exhibits multifractal scaling behavior due to the broad distribution of its values. We denote by $`P(x)`$ the probability density function of the values $`x_k`$ in the series. The distribution $`P(x)`$ does not affect the multifractality of a series on large scales $`s`$, if all moments

$$|x|^q_{\mathrm{}}^{\mathrm{}}|x|^qP(x)𝑑x$$
(24)
are finite. Here we choose a (normalized) power-law probability distribution function,

$$P(x)=\alpha x^{(\alpha +1)}\mathrm{for}\mathrm{\hspace{0.33em}1}x<\mathrm{}\mathrm{with}\alpha >0$$
(25)
and $`P(x)=0`$ for $`x<1`$, where already the second moment diverges if $`\alpha 2`$. In this case, the series exhibits multifractal scaling behavior on all scales. Note, that Eq. (25) becomes identical to a Levy distribution of class $`\alpha `$ for large values of $`x`$. The parameter $`\alpha `$ is not related to the Hölder exponent $`\alpha `$ in Eq. (15). The scaling properties of random walks with Levy distributed steps (Levy flights and Levy walks) have been analyzed in . The multifractal nature of Levy processes has been investigated in .

In order to derive the multifractal spectrum, let us consider $`s`$ uncorrelated random numbers $`r_k`$, $`k=1,\mathrm{},s`$, distributed homogeneously in the interval $`[0,1]`$. Obviously, the typical value of the minimum of the numbers, $`r_{\mathrm{min}}(s)\mathrm{min}_{k=1}^sr_k`$, will be $`r_{\mathrm{min}}(s)=1/s`$. It can be easily shown that the numbers $`r_k`$ are transformed into numbers $`x_k`$ distributed according to the power-law probability distribution function (25) by $`r_kx_k=r_k^{1/\alpha }`$. Thus, the typical value of the maximum of the $`x_k`$ will be $`x_{\mathrm{max}}(s)\mathrm{max}_{k=1}^sx_k=[r_{\mathrm{min}}(s)]^{1/\alpha }=s^{1/\alpha }`$.

If $`\alpha 2`$, the fluctuations of the profile $`Y(i)`$ \[Eq. (1)\] and the corresponding DFA variance $`F^2(s,\nu )`$ \[Eq. (2)\] will be dominated by the square of the largest value $`x_{\mathrm{max}}^2(s)=s^{2/\alpha }`$ in the segment of $`s`$ numbers, since the second moment of the distribution (25) diverges. Now the whole series consists of $`N_s\mathrm{int}(N/s)`$ segments of length $`s`$ and not just of one segment. For some segments $`\nu `$, $`[F^2(s,\nu )]^{1/2}`$ is larger than its typical value $`x_{\mathrm{max}}(s)=s^{1/\alpha }`$, since the maximum within the whole series of length $`N`$ is $`x_{\mathrm{max}}(N)=N^{1/\alpha }`$. In order to calculate $`F_q(s)`$ \[Eq. (4)\], we need to take into account the whole distribution $`P_s(y)`$ of the values $`y[F^2(s,\nu )]^{1/2}`$. Since each of the maxima in the $`N_s`$ segments corresponds to an actual number $`x_k`$ and these $`x_k`$ are random numbers from the power-law distribution (25), it becomes obvious, that the distribution of the maxima will have the same form, i. e. $`P_s(y)P(x=y)`$ for large $`y`$. Small values of $`y`$ are excluded because of the maximum procedure, but the large $`x_k`$ values are very likely to be identical to the maxima of the corresponding segments. Since the smallest maxima for segments of length $`s`$ are of the order of $`x_{\mathrm{max}}(s)=s^{1/\alpha }`$, the lower cutoff for $`P_s(y)`$ must be proportional to $`s^{1/\alpha }`$. From the normalization condition $`_{As^{1/\alpha }}^{\mathrm{}}P_s(y)𝑑y=1`$ (with an unimportant prefactor $`A<1`$) we get

$$P_s(y)=A^\alpha \alpha sy^{(\alpha +1)}.$$
(26)
Now $`F_q(s)`$ \[Eq. (4)\] can be calculated by integration from the minimum value $`As^{1/\alpha }`$ of $`y[F^2(s,\nu )]^{1/2}`$ to the maximum value $`N^{1/\alpha }`$. For $`sN`$ we obtain

$`F_q(s)\left[{\displaystyle _{As^{1/\alpha }}^{N^{1/\alpha }}}y^qP_s(y)𝑑y\right]^{1/q}`$
$`\left|A^\alpha sN^{q/\alpha 1}A^qs^{q/\alpha }\right|^{1/q}\{\begin{array}{cc}s^{1/q}(q>\alpha )\hfill & \\ s^{1/\alpha }(q<\alpha )\hfill & \end{array}.`$
Comparing with Eq. (5), we finally get

$$h(q)\{\begin{array}{cc}1/q(q>\alpha )\hfill & \\ 1/\alpha (q\alpha )\hfill & \end{array}.$$
(27)
Note that $`\tau (q)`$ follows a linear $`q`$ dependence, $`\tau (q)=q/\alpha 1`$ for $`q<\alpha `$, while it is equal to zero for $`q>\alpha `$ according to Eq. (13). Hence, the series of uncorrelated power-law distributed values has rather bi-fractal  instead of multifractal properties. Since $`h(2)=1/2`$ holds exactly for all values of $`\alpha `$, it is not possible to recognize the multifractality due to the broad power-law distribution of the values if only the conventional DFA is applied. The second moment shows just the uncorrelated behavior of the values. In a very recent preprint  this behavior has been interpreted as a failure of the DFA and corresponding non-detrending methods for series with a broad distribution, and another method to determine the exponent $`1/\alpha `$ has been proposed. We believe that a multifractal description with more than one exponent is required to characterize this kind of series, and thus any method calculating just one exponent will be insufficient for a full characterization.

Figure 4(a) shows the MF-DFA3 fluctuation functions for series of independent random numbers $`x_k[1,\mathrm{})`$ distributed according to Eq. (25) with $`\alpha =1`$. Since the scaling exponents $`h(q)`$ become very close to zero asymptotically for large positive values of $`q`$ according to Eq. (27), we must use the modified MF-DFA technique involving the double sum as described in the last paragraph of Subsection II.A. Hence, for this technical reason, $`\stackrel{~}{F}_q(s)/s`$ is calculated instead of $`F_q(s)`$. The corresponding slopes $`\stackrel{~}{h}(q)1`$ are identical to $`h(q)`$, see Eq. (7). In Fig. 4(b) the slopes $`h(q)`$ for series with $`\alpha =0.5`$, 1.0, and 2.0 are compared with the theoretical result Eq. (27), and nice agreement is observed.

## IV Comparison of the multifractality for original and shuffled series

### A Distinguishing the two types of multifractality

As already mentioned in the introduction, two different types of multifractality in time series can be distinguished. Both of them require a multitude of scaling exponents for small and large fluctuations. (i) Multifractality of a time series can be due to a broad probability density function for the values of the time series, and (ii) multifractality can also be due to different long-range correlations for small and large fluctuations. The example discussed in Subsection III.D, the uncorrelated multifractal series with a power-law probability density function, is of type (i), while the examples discussed in Subsections III.A – III.C are of type (ii), where the probability density function of the values is a regular distribution with finite moments .

Now, we would like to distinguish between these two types of multifractality. The most easy way to do so is by analyzing also the corresponding randomly shuffled series. In the shuffling procedure the values are put into random order, and thus all correlations are destroyed. Hence the shuffled series from multifractals of type (ii) will exhibit simple random behavior, $`h_{\mathrm{shuf}}(q)=0.5`$, i. e. non-multifractal scaling like in Fig. 1(b). For multifractals of type (i), on the contrary, the original $`h(q)`$ dependence is not changed, $`h(q)=h_{\mathrm{shuf}}(q)`$, since the multifractality is due to the probability density, which is not affected by the shuffling procedure. If both kinds of multifractality are present in a given series, the shuffled series will show weaker multifractality than the original one.

The effect of the shuffling procedure is illustrated in Fig. 5(a), where the MF-DFA2 fluctuation functions $`F_{10}^{\mathrm{shuf}}(s)`$ and $`F_{10}^{\mathrm{shuf}}(s)`$ are shown for shuffled series for three of the multifractal examples taken from the previous section. Random behavior, $`h_{\mathrm{shuf}}(q)=0.5`$, is observed for the series that were long-range correlated or generated from the dyadic random cascade model before the shuffling procedure \[upper four curves in Fig. 5(a)\]. In contrast, we observe the original multifractal scaling for the shuffled multifractal series with power-law probability density function $`P(x)x^2`$ \[lower two curves in Fig. 5(a)\]. The $`h_{\mathrm{shuf}}(q)`$ dependences are shown in Fig. 6, which can be compared with the corresponding slopes shown in Figs. 1(d), 3(b), and 4(b). Thus, the fluctuation analysis of the shuffled series, $`F_q^{\mathrm{shuf}}(s)`$, directly indicates the presence of type (i) multifractality, which is due to a broad probability distribution, by deviations from $`h_{\mathrm{shuf}}(q)=0.5`$.

Now we want to determine directly the magnitude of the (ii) multifractality, which is due to correlations. For that purpose we compare the fluctuation function for the original series, $`F_q(s)`$, with the result for the corresponding shuffled series, $`F_q^{\mathrm{shuf}}(s)`$. Differences between these two fluctuation functions directly indicate the presence of correlations in the original series. These differences can be observed best in a plot of the ratio $`F_q(s)/F_q^{\mathrm{shuf}}(s)`$ versus $`s`$ . Since the anomalous scaling due to a broad probability density affects $`F_q(s)`$ and $`F_q^{\mathrm{shuf}}(s)`$ in the same way, only multifractality due to correlations will be observed in $`F_q(s)/F_q^{\mathrm{shuf}}(s)`$. This is illustrated in Fig. 5(b) for the same three multifractal examples as in Fig. 5(a). In order not to have increased statistical errors in the results when considering the ratio $`F_q(s)/F_q^{\mathrm{shuf}}(s)`$ instead of $`F_q(s)`$ itself, $`F_q^{\mathrm{shuf}}(s)`$ can be calculated by averaging over a large number of randomly shuffled series generated from the same original series.

The scaling behavior of the ratio is

$$F_q(s)/F_q^{\mathrm{shuf}}(s)s^{h(q)h_{\mathrm{shuf}}(q)}=s^{h_{\mathrm{cor}}(q)}.$$
(28)
Note that $`h(q)=h_{\mathrm{shuf}}(q)+h_{\mathrm{cor}}(q)`$. If only distribution multifractality \[type (i)\] is present, $`h(q)=h_{\mathrm{shuf}}(q)`$ depends on $`q`$ and $`h_{\mathrm{cor}}(q)=0`$. On the other hand, deviations of $`h_{\mathrm{cor}}(q)`$ from zero indicate the presence of correlations, and a $`q`$ dependence of $`h_{\mathrm{cor}}(q)`$ indicates correlation multifractality \[type (ii)\]. If only correlation multifractality is present, $`h_{\mathrm{shuf}}(q)=0.5`$ and $`h(q)=0.5+h_{\mathrm{cor}}(q)`$. If both, distribution multifractality and correlation multifractality are present, both, $`h_{\mathrm{shuf}}(q)`$ and $`h_{\mathrm{cor}}(q)`$ depend on $`q`$.

### B Significance of the results

In Figs. 1-6 we have shown the results of the MF-DFA for single configurations of long time series. Now we address the significance and accuracy of the MF-DFA results for short series. How much do the numerically determined exponents $`h(q)`$ vary from one configuration (sample series) to the next, and how close are the average values to the theoretical values? In other words, how large are the statistical and systematical deviations of exponents practically determined by the MF-DFA for finite series? These questions are particularly important for short series, where the statistics is poor. If the values of $`h(q)`$ are determined inaccurately, the multifractal properties will be reported inaccurately or even false conclusions on multifractal behavior might be drawn for monofractal series.

To address the significance and accuracy of the MF-DFA results we generate, for each of the three examples considered already in Fig. 5, 100 series of length $`N=2^{13}=8192`$ and calculate $`h(10)`$, $`h(+10)`$, $`h_{\mathrm{shuf}}(10)`$, and $`h_{\mathrm{shuf}}(+10)`$ for each of these series. The corresponding histograms are shown in Fig. 7. For the long-range power-law correlated series with $`H=0.75`$ we find the following mean values and standard deviations of the generalized Hurst exponents:

$`h(10)=0.80\pm 0.03,h_{\mathrm{shuf}}(10)=0.56\pm 0.02,`$
$`h(+10)=0.72\pm 0.04,h_{\mathrm{shuf}}(+10)=0.48\pm 0.02.`$
The mean values for the original series are rather close to, but not identical to the theoretical value $`H=0.75`$. The mean value for $`q=10`$ is about two standard deviations larger than 0.75, while the value for $`q=+10`$ is slightly smaller. These deviations, though, certainly cannot indicate multifractality, since we analyzed monofractal series. Instead, they are due to the finite, random series, where parts of the series have slightly larger and slightly smaller scaling exponent just by statistical fluctuations. By considering negative values of $`q`$ we focus on the parts with small fluctuations, which are usually described by a larger scaling exponent . For positive values of $`q`$ we focus on the parts with large fluctuations usually described by a smaller value of $`h`$. Thus for short records we always expect a slight difference between $`h(10)`$ and $`h(+10)`$ even if the series are monofractal. If this difference is weak, one has to be very careful with conclusions about multifractality. Practically it is always wise to compare with generated monofractal series with otherwise similar properties before drawing conclusions regarding the multifractality of a time series. In addition to the statistical fluctuations of the $`h`$ values, the average $`h(10)`$ is usually determined slightly too large, while $`h(+10)`$ is slightly too small. The same behavior is obtained if the WTMM method is used instead of the MF-DFA, as we will show in Subsection V.C.

The same kind of difference is also observed for the average $`h_{\mathrm{shuf}}(10)`$ and $`h_{\mathrm{shuf}}(+10)`$ values. After all correlations have been destroyed by the shuffling, $`h_{\mathrm{shuf}}=0.5`$ is expected since the probability density is Gaussian with all finite moments. The deviations from $`h_{\mathrm{shuf}}=0.5`$ we observe for the finite random series are characteristic for monofractal series of this length ($`N=8192`$). Only for the second moment we obtain $`h_{\mathrm{shuf}}(2)=0.5`$ exactly if a sufficient number of series is considered.

For multifractal series generated from the dyadic random cascade model, Fig. 7(c,d) shows the histograms of the scaling exponents $`h(10)`$, $`h(+10)`$, $`h_{\mathrm{shuf}}(10)`$, and $`h_{\mathrm{shuf}}(+10)`$. Their averages and standard deviations,

$`h(10)=0.69\pm 0.04,h_{\mathrm{shuf}}(10)=0.57\pm 0.02,`$
$`h(+10)=0.54\pm 0.02,h_{\mathrm{shuf}}(+10)=0.48\pm 0.02.`$
have to be compared with the theoretical values from Eq. (23), $`h(10)=0.743`$ and $`h(+10)=0.567`$. Surprisingly, the mean $`h(10)`$ is smaller than the theoretical value in this example, but for the mean $`h(+10)`$ the deviation is similar to the deviation observed for the monofractal data in the previous example. Again, similar results are obtained with the WTMM method. For the shuffled series, the mean generalized Hurst exponents are practically identical to those for the shuffled monofractal series \[the average $`h_{\mathrm{shuf}}(10)`$ is larger by half the standard deviation\], and both are evidently consistent with monofractal uncorrelated behavior, $`h(q)=0.5`$, as discussed above. Hence, the series from the dyadic random cascade model show no signs of distribution multifractality and are characterized by correlation multifractality only.

The histograms of the scaling exponents for our last example, the power-law distributed random numbers with $`P(x)x^2`$, are shown in Fig. 7(e,f). The corresponding mean values and standard deviations,

$`h(10)=1.24\pm 0.09,h_{\mathrm{shuf}}(10)=1.26\pm 0.09,`$
$`h(+10)=0.11\pm 0.03,h_{\mathrm{shuf}}(+10)=0.11\pm 0.04.`$
show obviously no differences between original and shuffled series as expected for uncorrelated series. This indicates that the multifractality is due to the broad probability density function only. The values have to be compared with $`h(10)=1`$ and $`h(+10)=0.1`$ from Eq. (27). As usual, the average value of $`h(10)`$ is too large because we analyzed short series.

## V Comparison to the wavelet transform modulus maxima method

### A Brief description of the wavelet transform modulus maxima method

The wavelet transform modulus maxima (WTMM) method  is a well-known method to investigate the multifractal scaling properties of fractal and self-affine objects in the presence of nonstationarities. It is an application of the wavelet transform with continuous basis functions. One defines the wavelet-transform of a series $`x_k`$ of length $`N`$ by

$$W(n,s)=\frac{1}{s}\underset{k=1}{\overset{N}{}}x_k\psi [(kn)/s].$$
(29)
Note that in this case the series $`x_k`$ are analyzed directly instead of the profile $`Y(i)`$ defined in Eq. (1). Here, the function $`\psi (x)`$ is the analyzing wavelet and $`s`$ is, as above, the scale parameter. The wavelet is chosen orthogonal to the possible trend. If the trend can be represented by a polynomial, a good choice for $`\psi (x)`$ is the $`m`$-th derivative of a Gaussian, $`\psi ^{(m)}(x)=d^m(e^{x^2/2})/dx^m`$. This way, the transform eliminates trends up to $`(m1)`$th order.

Now, instead of averaging over all values of $`W(n,s)`$, one averages, within the modulo-maxima method, only the local maxima of $`|W(n,s)|`$. First, one determines for a given scale $`s`$, the positions $`n_i`$ of the local maxima of $`|W(n,s)|`$ as function of $`n`$, so that $`|W(n_i1,s)|<|W(n_i,s)||W(n_i+1,s)|`$ for $`i=1,\mathrm{},i_{\mathrm{max}}`$. Then one sums up the $`q`$th power of these maxima,

$$Z(q,s)=\underset{i=1}{\overset{i_{\mathrm{max}}}{}}|W(n_i,s)|^q.$$
(30)
The reason for this maxima procedure is that the absolute wavelet coefficients $`|W(n,s)|`$ can become arbitrarily small. The analyzing wavelet $`\psi (x)`$ must always have positive values for some $`x`$ and negative values for other $`x`$, since it has to be orthogonal to possible constant trends. Hence there are always positive and negative terms in the sum (29), and these terms might cancel. If that happens, $`|W(n,s)|`$ can become close to zero. Since such small terms would spoil the calculation of negative moments in Eq. (30), they have to be eliminated by the maxima procedure. In the MF-DFA, the calculation of the variances $`F^2(s,\nu )`$ in Eq. (2), i. e. the deviations from the fits, involves only positive terms under the summation. The variances cannot become arbitrarily small, and hence no maximum procedure is required for series with compact support.

In addition, the MF-DFA variances will always increase if the segment length $`s`$ is increased, because the fit will always be worse for a longer segment. In the WTMM method, in contrast, the absolute wavelet coefficients $`|W(n,s)|`$ need not increase with increasing scale $`s`$, even if only the local maxima are considered. The values $`|W(n,s)|`$ might become smaller for increasing $`s`$ since just more (positive and negative) terms are included in the summation (29), and these might cancel even better. Thus, an additional supremum procedure has been introduced in the WTMM method in order to keep the dependence of $`Z(q,s)`$ on $`s`$ monotonous: If, for a given scale $`s`$, a maximum at a certain position $`n_i`$ happens to be smaller than a maximum at $`n_i^{}n_i`$ for a lower scale $`s^{}<s`$, then $`W(n_i,s)`$ is replaced by $`W(n_i^{},s^{})`$ in Eq. (30). There is no need for such a supremum procedure in the MF-DFA.

Often, scaling behavior is observed for $`Z(q,s)`$, and scaling exponents $`\widehat{\tau }(q)`$ can be defined that describe how $`Z(q,s)`$ scales with $`s`$,

$$Z(q,s)s^{\widehat{\tau }(q)}.$$
(31)
The exponents $`\widehat{\tau }(q)`$ characterize the multifractal properties of the series under investigation, and theoretically they are identical to the $`\tau (q)`$ defined in Eq. (12)  and related to $`h(q)`$ in Eq. (13).

### B Examples for series with nonstationarities

Since the WTMM method has been developed to analyze multifractal series with nonstationarities, such as trends or spikes, we will compare its performance with the performance of the MF-DFA for such nonstationary series. In Fig. 8 the MF-DFA fluctuation function $`F_q(s)`$ and its scaling behavior are compared with the rescaled WTMM partition sum $`Z(q,s)`$ for the binomial multifractal described in Subsection III.B. To test the detrending capability of both methods, we have added linear as well as quadratic trends to the generated multifractal series. The trends are removed by both methods, if a sufficiently high order of detrending is employed. The deviations from the theoretical values of the scaling exponents $`h(q)`$ \[given by Eq. (20)\] are of similar size for the MF-DFA and the WTMM method. Thus, the detrending capability and the accuracy of both methods is equivalent.

We also obtain similar results for a monofractal long-range correlated series with additional spikes (outliers) that consist of large random numbers and replace a small fraction of the original series in randomly chosen positions. The spikes lead to multifractality on small scales $`s`$, while the series remains monofractal on large scales. Thus, the effects of the spikes are eliminated neither by the WTMM method nor by the MF-DFA, but both methods again give rather equivalent results.

### C Significance of the results

The last problem we address is a comparison of the significance of the results obtained by the MF-DFA and the WTMM method. The significance of the MF-DFA results has already been discussed in detail in Subsection IV.B. Here we will compare the significance of both methods for short and long series.

We begin with the significance of the results for random series involving neither correlations nor a broad distribution \[as in Fig. 1(b)\]. Fig. 9 shows the distribution of the multifractal Hurst exponents $`h(10)`$ and $`h(+10)`$ calculated by the MF-DFA as well as by the WTMM using the relation $`h(q)=[\widehat{\tau }(q)+1]/q`$ based on Eq. (13). Similar to the results presented in Fig. 7, we have analyzed 100 generated series of uncorrelated random numbers. In addition, we compare the results for the (relatively short) series length $`N=2^{13}=8192`$ and for $`N=2^{16}=65532`$. Ideally, both, $`h(10)`$ and $`h(+10)`$, should be equal to the Hurst exponent of the uncorrelated monofractal series, $`H=0.5`$. The histograms show that similar deviations as well as remarkable fluctuations of the exponents occur for both methods, as discussed in Subsection IV.B for the MF-DFA. We find the following mean values and standard deviations,

$`h(10)=\{\begin{array}{cc}0.55\pm 0.03\hfill & \text{for MF-DFA (}N=8k\text{)}\hfill \\ 0.52\pm 0.02\hfill & \text{for MF-DFA (}N=64k\text{)}\hfill \\ 0.58\pm 0.05\hfill & \text{for WTMM (}N=8k\text{)}\hfill \\ 0.56\pm 0.03\hfill & \text{for WTMM (}N=64k\text{)}\hfill \end{array}`$
$`\mathrm{and}h(+10)=\{\begin{array}{cc}0.49\pm 0.03\hfill & \text{for MF-DFA (}N=8k\text{)}\hfill \\ 0.49\pm 0.02\hfill & \text{for MF-DFA (}N=64k\text{)}\hfill \\ 0.46\pm 0.04\hfill & \text{for WTMM (}N=8k\text{)}\hfill \\ 0.48\pm 0.02\hfill & \text{for WTMM (}N=64k\text{)}\hfill \end{array}.`$
As already discussed in Subsection IV.B, the deviations of the average $`h(q)`$ values from $`H=0.5`$ do not indicate multifractality. For the WTMM method and short series, one has to be very careful in order not to draw false conclusions from results like $`h(10)=0.58`$ and $`h(+10)=0.46`$. The corresponding results of the MF-DFA are closer to the theoretical value.

Figure 10 shows the distribution of the multifractal scaling exponents $`h(10)`$ and $`h(+10)`$ calculated for generated multifractal series from the binomial model with $`a=0.75`$ described in Subsection III.B. Like for Fig. 9, 100 generated series have been analyzed for each of the histograms. Now, the differences between the distributions of $`h(10)`$ and $`h(+10)`$ are much larger, indicating multifractality. We find

$`h(10)=\{\begin{array}{cc}1.88\pm 0.06\hfill & \text{for MF-DFA (}N=8k\text{)}\hfill \\ 1.89\pm 0.03\hfill & \text{for MF-DFA (}N=64k\text{)}\hfill \\ 1.86\pm 0.05\hfill & \text{for WTMM (}N=8k\text{)}\hfill \\ 1.89\pm 0.02\hfill & \text{for WTMM (}N=64k\text{)}\hfill \end{array}`$
$`\mathrm{and}h(+10)=\{\begin{array}{cc}0.50\pm 0.02\hfill & \text{for MF-DFA (}N=8k\text{)}\hfill \\ 0.51\pm 0.01\hfill & \text{for MF-DFA (}N=64k\text{)}\hfill \\ 0.46\pm 0.01\hfill & \text{for WTMM (}N=8k\text{)}\hfill \\ 0.47\pm 0.01\hfill & \text{for WTMM (}N=64k\text{)}\hfill \end{array}.`$
These values must be compared with the theoretical values $`h(10)=1.90`$ and $`h(+10)=0.515`$ from Eq. (20). Again, the MF-DFA results turn out to be slightly more significant than the WTMM results. The MF-DFA seems to have slight advantages for negative $`q`$ values and short series, but in the other cases the results of the two methods are rather equivalent. Besides that, the main advantage of the MF-DFA method compared with the WTMM method lies in the simplicity of the MF-DFA method.

## VI Conclusion

We have generalized the DFA, widely recognized as a method to analyze the (mono-) fractal scaling properties of nonstationary time series. The MF-DFA method allows a reliable multifractal characterization of multifractal nonstationary time series. The implementation of the new method is not more difficult than that of the conventional DFA, since just one additional step, a $`q`$ dependent averaging procedure, is required. We have shown for stationary signals that the generalized (multifractal) scaling exponent $`h(q)`$ for series with compact support is directly related to the exponent $`\tau (q)`$ of the standard partition function-based multifractal formalism. Further, we have shown in several examples that the MF-DFA method can reliably determine the multifractal scaling behavior of the time series, similar to the WTMM method which is a more complicated procedure for this purpose. For short series and negative moments, the significance of the results for the MF-DFA seems to be slightly better than for the WTMM method.

Contrary to the WTMM method, the MF-DFA method as described in Subsection II.A requires series of compact support, because the averaging procedure in Eq. (4) will only work if $`F^2(s,\nu )>0`$ for all segments $`\nu `$. Although most time series will fulfill this prerequisite, it can be overcome by a modification of the MF-DFA technique in order to analyze data with fractal support: We restrict the sum in Eq. (4) to the local maxima, i. e. to those terms $`F^2(s,\nu )`$ that are larger than the terms $`F^2(s,\nu 1)`$ and $`F^2(s,\nu +1)`$ for the neighboring segments. By this restriction all terms $`F^2(s,\nu )`$ that are zero or very close to zero will be disregarded, and series with fractal support can be analyzed. The procedure reminds slightly of the modulus maxima procedure in the WTMM method (see Subsection V.A). There is no need, though, to employ a continuously sliding window or to calculate the supremum over all lower scales for the MF-DFA, since the variances $`F^2(s,\nu )`$, which are determined by the deviations from the fit, will always increase when the segment size $`s`$ is increased. In the maxima MF-DFA procedure the generalized Hurst exponent $`h(q)`$ defined in Eq. (5) will depend on $`q`$ and even diverge for $`q0`$ for monofractal series with non-compact support. Thus, it is more appropriate to consider the scaling exponent $`\tau (q)`$, calculating

$$\underset{F^2(s,\nu 1)<F^2(s,\nu )F^2(s,\nu +1)}{}[F^2(s,\nu )]^{q/2}s^{\tau (q)}.$$
(32)
This extended MF-DFA procedure will also be applicable for data with fractal support.

In a later work we will apply the MF-DFA method to a range of physiological and meteorological data.

### Acknowledgements

We would like to thank Yosef Ashkenazy for useful discussions and the Deutsche Akademischer Austauschdienst, the Deutsche Forschungsgemeinschaft, the German Israeli Foundation, and the Minerva Foundation for financial support.