*=======================================================================
* dfa.f - zero-order detrended fluctuation analysis
* search for long term correlations in a data sequence
*
* by Andy Allinger, 2012, released to the public domain
* This program may be used by any person for any purpose.
*
* refer to:
* Jan W. Kantelhardt
* Fractal and Multifractal Time Series
* 4 April 2008
* published at: arxiv.org
*
* COMMENTS: Detrended fluctuation analysis (DFA) computes a term called
* fractal dimension. "Fractal dimension" has been defined in several ways.
* The value DFA yields is not the probably more interesting Hausdorff
* dimension. It merely demonstrates that the variance of short segments
* of a time series, compared to longer segments, is different from what
* would be expected by chance if there were no long-range correlations.
* The present program calculates the 0-order DFA, which can use the
* short formula for variance. For my purposes, the program returns a
* meaningless number, but it sure does it fast. Maybe you know what it's
* good for.
*
*=======================================================================
* I/O LIST
* Name Type In/Out Description
* X[n] REAL in time series of data points
* n INTEGER in length of the series
* a REAL out fractal dimension
* fault INTEGER out error code 0 = no errors
*
*=======================================================================
SUBROUTINE dfa(X, n, a, fault)
IMPLICIT NONE
INTEGER n, fault
REAL X, a
DIMENSION X(n)
* constant
INTEGER rlim ; PARAMETER (rlim = 9) ! how many segment lengths
* local variables
INTEGER
& i, j, ! counters
& jlim, ! index of slim in R[]
& p, ! # segments compared each length
& R(rlim), ! segment lengths
& s, slim ! length of segment to compare
DOUBLE PRECISION
& ave, sum, sum2, ! average, sum, sum of squares
& F(rlim), v, ! fluctuation
& var, ! variance in a segment
& x1, sy, sxy, ! least squares terms
& Y(n), y1, ! cumulative X
& CY(0:n), CY2(0:n) ! cumulative Y
DATA R / 10, 20, 40, 80, 160, 320, 640, 1280, 2560 /
*=======================================================================
* begin
*=======================================================================
fault = 0
IF (n < 30) THEN ! never mind the short stuff
a = 0.
fault = -1
RETURN
END IF
slim = n / 4 ! limit of statistical significance?
* Let X have zero mean
sum = 0.
DO j = 1, n ; sum = sum + X(j) ; END DO
ave = sum / n
DO j = 1, n ; X(j) = X(j) - ave ; END DO
*=======================================================================
* accumulate
*=======================================================================
sum = 0.
DO j = 1, n
sum = sum + X(j)
Y(j) = sum
END DO
sum = 0. ; sum2 = 0.
CY(0) = 0. ; CY2(0) = 0.
DO j = 1, n
y1 = Y(j)
sum = sum + y1
sum2 = sum2 + y1 * y1
CY(j) = sum
CY2(j) = sum2
END DO
*=======================================================================
* detrended fluctuation analysis
*=======================================================================
DO j = 1, rlim ! initialize F and P
IF (R(j) < slim) jlim = j ! index of maximum segment length
END DO
DO j = 1, jlim ! for each segment length
s = R(j)
v = 0.
p = 0
DO i = 0, n-s, s ! for each segment
sum = CY(i+s) - CY(i)
sum2 = CY2(i+s) - CY2(i)
ave = sum / s ! average
var = sum2 / s - ave * ave ! variance
IF (var < 0) PRINT *, 'affine: negative variance!'
v = v + var ! fluctuation
p = p + 1 ! count
END DO ! next i
DO i = n, s, -s ! the same thing, backwards
sum = CY(i) - CY(i-s)
sum2 = CY2(i) - CY2(i-s)
ave = sum / s
var = sum2 / s - ave * ave
v = v + var
p = p + 1
END DO
F(j) = SQRT(v / p) ! fluctuation
END DO ! next j
* slope of a least squares line to find a, where
* F ~ s^a
sum = 0. ; sum2 = 0. ; sy = 0. ; sxy = 0.
DO j = 1, jlim
s = R(j)
x1 = FLOAT(s)
x1 = LOG(x1)
y1 = LOG(F(j))
sum = sum + x1
sum2 = sum2 + x1 * x1
sy = sy + y1
sxy = sxy + x1 * y1
END DO
x1 = jlim * sum2 - sum * sum ! solve the normal equations
IF (x1 < 1E-7) THEN
fault = 1
RETURN
END IF
y1 = jlim * sxy - sum * sy
a = SNGL(y1 / x1)
END ! of affine