Fractal Analysis of Stock Index and Electrocardiograph

Fractal Analysis of Stock Index and Electrocardiograph

Sy-Sang Liaw,1 _{Feng-Yuan Chiu,}1 _{Cheng-Yen Wang,}1 _{and You-Hsien Shiau}2

1

Department of Physics, National Chung-Hsing University, 250 Guo-Kuang Road, Taichung 402, Taiwan

2

Graduate Institute of Applied Physics, National Chengchi University, Taipei 11605, Taiwan

(Received July 13, 2010)

We use the modified inverse random midpoint displacement (mIRMD) method to calculate the fractal dimensions of time sequences of arbitrary length. We show that monofractals, crossover-fractals, and the superposition of a fractal, with either a differentiable or periodic function, can be easily identified from the results calculated by the mIRMD. We compare our method with the method of detrended fluctuation analysis (DFA) and point out their differences. We apply our method to real financial and physiological data by considering the movements of the S&P 500 stock index and the electrocardiograph of a patient with ventricular fibrillation.

PACS numbers: 05.45.Tp, 05.45.Df, 89.65.Gh

I. INTRODUCTION

The hallmark of fractals is their structure under unlimited magnification. Thus, to be precise, the fractal dimension of a curve has to be defined within the limits of an infinitely small scale. In practice, however, a real time sequence consists of only a finite number of points. When one calculates the fractal dimension of a finite time sequence, one assumes implicitly that the sequence will vary similarly at scales less than the minimal interval of the sequence. In other words, one assumes that the time sequence can be interpolated by a monofractal. This assumption is, in general, false. Real non-stationary data normally varies differently at different scales, so that it generally has scale-dependent fractal dimensions. Therefore, one can characterize a time sequence better by a fractal spectrum, showing how the sequence varies at different scales, rather than by one single “average” value.

Since the invention of fractal geometry [1], many methods have been proposed for the calculation of fractal dimensions. It is surprising that none of the existing practical methods, such as the box-counting dimension, correlation dimension, and information di-mension, is related directly to the geometric properties of the fractals. If we take the simple case of curves that can be represented by fractal functions as an example, without any cal-culation we know intuitively that the fractal dimension of the curve is small if the curve looks smoother, but it is large if the curve makes turns more dramatically. None of the above mentioned methods meets this intuitive understanding. Recently, Liaw and Chiu [2] proposed a method for calculating the fractal dimensions of a time series from a geometric point of view by evaluating the turning angles at every point. This method, called IRMD, can also evaluate midpoint displacements between any two points on the curve, so that the calculation procedure is equivalent to the inverse of the procedure in the well-known

ran-http://PSROC.phys.ntu.edu.tw/cjp 814 2010 THE PHYSICAL SOCIETYc OF THE REPUBLIC OF CHINA

0

s

s

s

G

G

G

FIG. 1: Piecewise interpolation Lk of a function f (x) (red) at level 0 (dotted), 1 (dashed), and 2 (solid). The midpoint displacements Gjk and scales sk are as indicated.

dom midpoint displacement (RMD) [3] method used for generating fractal curves. In this report, a modified version of IRMD, called mIRMD, which is more suitable for calculating the scale-dependent fractal dimensions of time sequences of arbitrary length, has been used. This paper has been organized as follows: the IRMD method and its modified ver-sion, mIRMD, are described in Sec. II; our method is tested on some known monofractal functions, including random walks and fractional Brownian motions (fBms), in Sec. III; in Sec. IV we compare the mIRMD method with the popular detrended fluctuation analysis (DFA) method for calculating Hurst exponents; in Sec. V we discuss the relation between the Hurst exponent and the fractal dimension, and show explicitly that curves with Hurst exponents between 1 and 2 are actually differentiable but not twice-differentiable functions; we then describe how to identify an additive noise or a trend from a fractal sequence based on the mIRMD results in Sec. VI; and in Sec. VII we apply the mIRMD to a movement of the S&P 500 stock index and the electrocardiograph of a patient with ventricular fibrillation. Sec. VIII is the conclusion.

II. THE METHOD OF IRMD AND ITS MODIFICATION

Consider a monofractal f sampled by a sequence of N = 2M _{+ 1 points [i, f (i)], i =}

0, 1, 2, · · · , N − 1. We construct a set of the piecewise linear functions Lk, k =

0, 1, 2, · · · , M that connect the points [jsk, f (jsk)] and [(j + 1)sk, f ((j + 1)sk)], where

sk = 2M −k = 2M/2k, j = 0, 1, 2, · · · , 2k − 1, and define ∆k as the sum of the

ab-solute values of the differences Gj_k between Lk and Lk+1 at all points [jsk+1, f (jsk+1)],

j = 0, 1, 2, · · · , 2k+1 (Fig. 1). We found that [2]

∆k∝ (sk)1−D, (1)

where D > 1 is the fractal dimension of the monofractal f . The IRMD method calculates ∆k

and uses Eq. (1) for determining D. We have shown [2] that the IRMD is fast and accurate in calculating D for monofractals such as Weierstrass functions, random walk trajectories,

1

3

5

7

6

4

2

FIG. 2: For a sequence of N = 9 points (big black dots), mIRMD averages over N − s = 7 values of midpoint displacement G(s) (solid black) at the scale s = 2, while IRMD sums over (N − 1)/s = 4 non-overlapped intervals (#1, 3, 5, and 7) only.

and for fBms of arbitrary Hurst exponents smaller than 1. There is, however, one drawback for the IRMD. In practice we found that one normally needs at least N = 29_{+ 1 = 513}

points in IRMD in order to obtain satisfactory results. This is because for N = 2M + 1, we have only M points in the log-log plot of ∆k versus sk, and the first 6 values of ∆k,

k = 0, 1, 2, · · · , 5 may not be reliable because they are a sum of not enough numbers (at most 25) of Gj_k values. Besides, real data do not always have exactly 2M + 1 points, so we are forced to cut the data length short. For example, for a sequence of 500 points, we can only calculate the fractal dimension of a partial sequence with N = 28_{+ 1 = 257, which}

might yield unsatisfactory results.

Notice that we add 2k values of Gj_k to obtain ∆k. Thus we have

∆k= 2khGki, (2)

where hGki is the average of the values Gj_k. To find D, we can equivalently calculate hGki

which is proportional to (sk)2−D:

hGki ∝ (sk)2−D. (3)

Now, hGki stands for the average displacement of the points [x, f (x)] from the midpoints

of the line segments connecting [x − sk+1, f (x − sk+1)] and [x + sk+1, f (x + sk+1)] by

restricting x to be any of 2k possible values: x = jsk, j = 1, 3, 5, · · · , 2k+1− 1, so that no

intervals [x − sk+1, x + sk+1] of length sk are overlapped. However, from the point of view

of calculating the average midpoint displacements of a time sequence, it would be more reasonable to consider the displacements of all intervals of sk, no matter whether or not

they are overlapped. There are in total N − skintervals of length sk(Fig. 2). Note that the

value of N − sk is larger than N/2 for all values of k except 0. The reliability of the average

value hGki of these N − sk(k 6= 0) values of displacement is more enhanced than the value

0.1 0.01 0.001 0.0001 0.00001 <G(s)> 2 4 6 8 20 200 2000 s

FIG. 3: Comparison between the results of mIRMD (dot) and those of IRMD for white noise (circle), random walk (diamond), and the double differentiable function f (x) = sin(100x) (square). The slopes of the fitting lines for these three sequences were close to 0, 0.5, and 2, indicating that their fractal dimensions were 2, 1.5, and 1,, respectively.

restrict the interval lengths to be an integral power of 2. An interval has to be a multiple of 2 though, so that it can have a midpoint x for calculating the displacement. Accordingly, N does not have to be in the form N = 2M+ 1 either. Summing up these changes, we see that the modified method of the IRMD, called mIRMD for short, is applicable to any sequence of arbitrary length N . In mIRMD we calculate the average midpoint displacement: hG(s)i, where s = 2, 4, 6, · · · , N/2 stands for the length of the interval, by

hG(s)i = 1 N − s N −1−s/2 X x=s/2     f (x) −f (x − s/2) + f (x + s/2) 2     ∝ sα. (4)

There will be a large number (N /4, as compared to log₂(N ) in IRMD) of calculated values of hG(s)i each of which is an average of more than N /2 midpoint displacements. The slope α of log(hG(s)i) versus log(s) plot for a monofractal of dimension D > 1 is equal to 2 − D (Eq. (3)). For the special case when f is a twice differentiable function [2] (D = 1), α = 2.

III. TESTS ON RANDOM WALKS AND FBMS

We test the mIRMD on the trajectory of a random walk, a twice differentiable func-tion, f (x) = sin(100x), and a white noise, and compare the results with that of IRMD. Fig. 3 shows the results for time sequences of length N = 2049. We see that in a small-scale range the IRMD, although having a smaller number of calculated values, is as good as mIRMD in determining the fractal dimensions of the sequences. In a large-scale range, re-sults from IRMD deviate slightly from that of mIRMD. Since the fractal dimensions of the twice differentiable functions, random walks, and white noises are 1, 1.5, and 2 respectively, we obtain in mIRMD, as expected, the slopes α of their log(hG(s)i) versus log(s) plots to be close to 2, 0.5, and 0, respectively.

We next check how well the mIRMD works on fBms with arbitrary dimensions be-tween 1 and 2. There exist several different algorithms for generating fBms. Here we use the algorithm proposed originally by Mandelbrot [4]. First, we create a time sequence of increments ξi by the formula

ξi(H, M, m) = m−_H Γ(H + 0.5)    m X j=1 jH−0.5θ_{1+m(M +i)−j} + m(M −1) X j=1 (m + j)H−0.5_{− j}H−0.5_{ θ} 1+m(M −1+i)−j    , (5)

where Γ is the Gamma function defined by Γ(p) =R∞ 0 xp−1e

−_x

dx. The Hurst exponent H is a constant between 0 and 1, and {θk} is a set of random numbers of which the probability

distribution is Gaussian with a unit variance and a zero mean. M is an integer that should be moderately large, and m indicates the number of the fractional steps in each unit time. The larger the values of M and m are, the more likely the distribution of {ξi} is close to a

Gaussian distribution. An fBm xl l = 1, 2, · · · , N , is given by the accumulation of these

increments: xl= l X i=1 (ξi− hξi), (6) where hξi = 1 N N X i=1 ξi.

We used M = 1000 and m = 50 and generated 30 fBms of length N = 500 for each of H = 0.1, 0.2, . . . , 0.9. We used mIRMD to find the average fractal dimension D of these 30 fBms and used the relation D = 2 − H to obtain the calculated value: HmIRMD. The

results are plotted against H in Fig. 4. We see that HmIRMD matches H very well, except

for the small values of H. The value of HmIRMD can be improved for small H by increasing

the value of m in generating ξi (Eq. (5)). For example, for the case of H = 0.1, when m =

1000 was used, we obtained HmIRMD= 0.13 instead of 0.18 as shown in Fig. 4 where m =

50 was used. However, it took a longer cpu time to create fBms when a larger m was used.

IV. COMPARISON WITH DFA

The detrended fluctuation analysis (DFA) was first introduced [5] in 1994 and has become a popular method in analyzing the long-term correlation of time series. In DFA one divides the time sequence into (N − 1)/b non-overlapped boxes of size b. Note that the box size can be any integer b = 1, 2, 3, · · · so that it can be related to the interval s

H

H

FIG. 4: Average Hurst exponents calculated by mIRMD and DFA, respectively, for 30 fBms of lengthN = 500. Standard deviations of the mIRMD results are shown in error bars. Error bars for DFA, which are about the same size as those for mIRMD, are not shown for clarity.

we used in mIRMD by s = 2b. One then defines the local trend in each box by fitting the data in the box with a polynomial function of the order n. The detrending effects of using different values of n have been investigated [6, 7]. Here we use n = 1 for simplicity. Let tl(s) be the local linear trend in the box where xl is located. One defines the detrended

sequence {Xl(s)} by subtracting tl(s) from {xl}:

Xl(s) = xl− tl(s). (7)

The mean square displacements F2_{(s) are then calculated by}

F2(s) = 1 N N X l=1 [Xl(s)]2. (8)

F (s) has a power relation to the scale s:

F (s) ∝ sH. (9)

The scaling factor is similar to the Hurst exponent, and we simply use H for it in this paper.

The DFA results for the fBms we had generated above are also plotted in Fig. 4. We see the results of DFA and mIRMD have about the same accuracy in calculating the H of the fBms. Coronado and Carpena [8] have shown that the DFA can provide accurate results for data with length as short as N = 28 _{for most cases. Fig. 4 shows that the}

mIRMD is as good as the DFA for N = 500. At this point, we would like to point out two differences between the mIRMD and DFA. The mIRMD calculates the average midpoint

D=2-slope = 1.800 H=slope = 0.203 0.6 0.5 0.4 0.3 0.2 0.1 <G(s)> or F(s) 2 4 68 20 200 2000 s

FIG. 5: Results of the mIRMD (circle) and DFA (square) in determining, respectively, the D and H of the Weierstrass function: f (x) = P∞

n=02

−n(2−D)

cos(2n_{x), D = 1.8. The mIRMD results} fluctuated at the large scales, while the DFA results deviated from a monofractal property in the small scales.

displacements hG(s)i and uses the relation Eq. (3) to obtain the fractal dimension D, while the DFA calculates the mean variances F2_{(s) and uses the relation Eq. (9) to obtain the}

Hurst exponent H. The first difference between these two methods is that one gives D and the other yields H. Direct comparison of their results may not be possible if the relation between D and H for time sequences is not known. It is generally believed that the relation D = 2 − H is good for a series whose return distributions are Gaussian, such as for the cases of random walks and fBms. We had assumed the relation in the comparison shown in Fig. 4. The second important difference between the mIRMD and DFA is that the calculation of hG(s)i in mIRMD is linear in f (x), while the calculation of F (s) in DFA involves quadratures of f (x). The linear formulation in mIRMD not only speeds up the calculation but more importantly makes the analysis of non-monofractals, such as crossover-fractals and monofractals with a trend, much easier.

Let us compare the results of mIRMD and DFA when applying the Weierstrass function [9, 10] W (x) = P∞

n=02

−n(2−D)

cos(2nx) with dimension D = 1.8. We generate N = 216_{+ 1 values of W (x}

i) by adding the first 51 terms of the series at xi = i/(N − 1),

i = 0, 1, 2, · · · , N − 1. The calculated values of hG(s)i and F (s) are plotted in Fig. 5. The value of hG(s)i fluctuates along the linear fitted line with a slope 0.200, showing that the calculated fractal dimension is DmIRMD = 1.800. The fitted line for the values of F (s),

when neglecting the first 5 points, yields H = 0.203. If the relation D = 2 − H is used, we have DDFA = 1.793. Both mIRMD and DFA work well for the Weierstrass function in

determining the overall fractal dimension. Examined closely, we see that the DFA results deviate from a monofractal property in the small scales—an intrinsic nuisance of DFA that

has been reported [11]. The mIRMD results, on the other hand, fluctuate in the large scales. The fluctuation is actually due to the quasi-periodic property of the Weierstrass function [12].

Note that the formulation of the mIRMD has automatically included the effect of the first order detrending. To see this, let us assume the first order trend in an interval of size s is given by the function ax +b. The midpoint displacement calculated using the detrended data is given by

[f (x) − (ax + b)] −[f (x − s/2) − (a(x − s/2) + b)] + [f (x + s/2) − (a(x + s/2) + b)] 2

=f (x) − f (x − s/2) + f (x + s/2)

2 ,

(10) which is equal to the midpoint displacement without detrending. We therefore see that the mIRMD results are normally similar to the results of DFA1 (DFA with 1storder detrending). However, when time sequences are periodic (e.g., sinusoidal functions) or quasi-periodic (e.g., the Weierstrass functions), mIRMD can reveal their periods (see Sec. VI) while DFA1 tends to average out the periodic variation and obtains a whitenoise-like result at large scales.

In practice, because only a few operations are needed in calculating the midpoint displacement for an interval, mIRMD runs very fast and therefore is capable of carrying out calculations at all possible scales s = 2, 4, 6, . . . N /2 and of averaging over all overlapping windows of the same size within a short time period. The efficiency advantage makes mIRMD a better candidate than DFA for analyzing two- or three-dimensional data.

V. COMPARISON WITH THE CALCULATION OF THE HURST EXPONENTH1 For more than 20 years, the Hurst exponent has been calculated for many non-stationary time sequences to reveal the correlation properties of the sequences. In addition to the value calculated by the DFA as described above, the Hurst exponent can also be calculated by other formulas. The simplest form of the Hurst exponent is denoted by H1

and calculated by the formula

Y = 1 N − s N −s X x=0 |f (x + s) − f (x)| ∝ sH1_{. (11)}

Notice that the Eq. (11) formula is involved with the absolute values of the linear form of the function f (x), similar to our formulation in Eq. (4). Using the argument presented in Ref. [2], one can show that for a fractal f (x) of dimension 1 < D ≤ 2, H1 is given by

H1 = 2 − D, and is equal to α. However, the difference between Eq. (11) and Eq. (4) can

be easily seen. At the limit s → 0, Eq. (11) is summing 1st derivatives of f (x) at all x’s, while Eq. (4) is summing 2nd derivatives. The values of H1 and α differ when f (x) is a

differentiable function of which the 1st _{derivative, f}0

(x), is a fractal. Time sequences that can be approximated by differentiable function but have 2nd derivative nowhere have been observed in heartbeat time series [13] and wind-speed data [14]. For this kind of function H1 is always 1, while α ranges from 1 to 2, and a fractal dimension of f0(x) is given by

D = 3 − α (Fig. 6). Note that the Hurst exponent H for the differentiable but not twice-differentiable functions calculated by DFA also have values between 1 and 2. No connection between 1 < H < 2 and the fractal property of the function or its derivative has been made previously. 0 0.2 0.4

f(x)

_x

f’(x)

_x

FIG. 6: Sequence with H1= 1.0, α = 1.4 can be interpolated by a differentiable function (above). However, the derivative of the sequence was a fractal with dimension D = 3 − α = 1.6 (below).

We generated differentiable fBms with no 2nd derivative anywhere, using the method described in Sec. II, by choosing the value of H in the range of [1, 2]. These fBms are of dimension 1. The value of H1 calculated by Eq. (11) is always 1 for all values of H in

the range of [1, 2]. The value of α calculated by Eq. (4) is given by α = 2 − H. Another example of a differentiable function with no 2nd derivatives anywhere is given by

V (x) =X∞

n=02

−n(3−D)

sin(2nx), 1 ≤ D ≤ 2. (12)

The 1st derivative of V (x) is nothing but the well-known Weierstrass function W (x) of dimension D , as given in Sec. IV. The values H1 and α for V (x) calculated by Eqs. (11)

and (4), respectively, are H1 = 1 (independent of the value D given in Eq. (12)) and

α = 2 − D.

VI. NOISE AND TRENDS

Often, real non-stationary sequences are either contaminated by white noise or are superposed with a trend which is differentiable. Because Eq. (4) is linearly dependent on the function f (x), we easily identified the dominant part of the superposed components from the log-log plot of hG(s)i versus s. For example, an additive white noise can be

easily identified at small scales because it has a constant contribution to hG(s)i, while contributions from other components having fractal dimensions less than 2 decrease to zero as s goes to 0. On the other hand, since the α value of a differentiable trend (1 ≤ α ≤ 2) is larger than that of a fractal (0 ≤ α < 1), we could easily tell whether or not a superposed trend existed from the log-log plot of hG(s)i versus s at large scales. Fig. 7 is the log-log plot of hG(s)i versus s for a monofractal superposed, respectively, with white noise and a twice-differentiable trend. A crossover was clearly seen in each case.

α = 0.76 α = 0.03 0.01 0.02 0.05 0.1 0.15 <G(s)> 2 4 6 8 20 200 2000 s α = 1.97 α = 0.64 1 0.1 0.01 0.001 <G(s)> 2 4 6 8 20 200 2000 s

FIG. 7: Log-log plot of hG(s)i versus s for a monofractal of dimension 1.3 superposed, respectively, with a white noise (α = 0) (left) and a twice-differentiable trend (α = 2) (right).

A periodic trend, often observed in environmental and physiological data, requires special attention, since its average midpoint displacement behaves differently from a power law. For a periodic function f (x) of period T , the midpoint displacement at the scale s = 2T is zero:

f (x) −f (x − s/2) + f (x + s/2)

2 = f (x) −

f (x − T ) + f (x + T )

2 = f (x) − f (x) = 0. (13) It also vanishes at the scale of any integral multiple of 2T . Thus, one half of the interval of s between two consecutive minima of hG(s)i is equal to the period T . Besides, these minimal values of hG(s)i are equal. Consider now that f (x) is the superposition of a monofractal of dimension D = 1.4 and a periodic trend sin(1000x). The log-log plot of hG(s)i versus s for f (x) is plotted in Fig. 8. The results indicated that for scales less than T , the series behaved as a fractal of dimension D , as expected. For scales larger than T , the average midpoint displacement began to oscillate, and it reached a minimum at every scale that was an integral multiple of 2T , just as for a pure periodic function. However, the minimal values were not the same. These minimal points in the plot actually lined up with a slope equal to that for the monofractal (Fig. 8), since the contributions from a periodic trend is zero at these scales. That is, the fractal property could still be identified at large scales when a periodic trend was superposed onto a fractal by looking at the minimal values of hG(s)i.

= 0.60 α = 0.61

2T

FIG. 8: Log-log plot of hG(s)i versus s for the superposition of a monofractal of dimension D = 1.4 and a periodic double-differentiable trend of period T .

VII. APPLICATIONS TO REAL NON-STATIONARY SEQUENTIAL DATA

In this section, two sets of real time sequences have been calculated as a demonstration of the usage of mIRMD. The stock index has been studied by many physicists recently. Here, we have considered the S&P 500 stock index at one minute intervals during the period from 1986 to 1992 [15]. The second series of data we considered was the electrocardiographic recordings (ECG) of a patient with ventricular fibrillation. A better understanding of ECG is no doubt very important in the diagnosis in heart disease. Here we have described some characteristics of this data based on mIRMD. A more complete analysis and discussion will be published elsewhere.

VII-1. Financial data: S&P 500

We calculated the fractal dimensions of the S&P 500 stock index at one minute intervals. The time series of the stock index have been shown [16] previously to have a crossover-fractal property [17–20]. In Fig. 9, we show the results of mIRMD for the stock index sequence in September 1987. There were, in total, N = 7725 points in the sequence. According to the calculations of mIRMD, the S&P had a fractal dimension of 1.05 for a time scale shorter than 25 minutes, and 1.36 for longer scales. The Hurst exponent calculated by DFA had a crossover at 33 minutes, below which H = 0.93, and above which H = 0.65. Qualitatively, mIRMD and DFA had consistent results in showing the crossover-fractal property of the stock index at one minute intervals. Interestingly, the crossover behavior was present and qualitatively the same for almost any period of the S&P 500 stock index during 1986 and 1992, whether the data length was for a year, a month, a week, or even a day. On the other hand, the crossover property was not universal for all stock indexes

around the world. For example, we found that the German stock index XDAX [15], during the period from September 2008 to June 2009, was well approximated by a monofractal instead. We are currently investigating the possible implication of the crossover behavior found in stock index sequences.

slope = 0.95 slope = 0.64 slope = 0.93 slope = 0.65 1 0.1 5 0.5 0.05 <G(s)> or F(s) 2 4 68 20 200 2000 s

FIG. 9: Time sequence of the S&P500 stock index for September, 1987 at one min intervals was a crossover-fractal. It had a fractal dimension 1.05 for a time scale shorter than 25 min and 1.36 for larger scales according to mIRMD (circle). The Hurst exponents calculated by DFA (square) were 0.93 and 0.65, respectively, before and after the crossover at 33 min.

VII-2. Physiological data: ECG

Fig. 10 describes an ECG recording before/after the onset of ventricular fibrillation (VF), which was obtained from the public internet site, www.physionet.org. The cu ven-tricular tachyarrhythmia database (cudb) contains ECG recordings of VF patients and the sampling rate in this database is 250 Hz. Fig. 10 shows the monomorphic, rather than the polymorphic VF pattern.

We calculated the fractal spectra for a typical ECG sequence when the patient was in a normal condition and another when VF occurred. The log-log plots of hG(s)i versus s are shown in Figs. 11(a) and 11(b), respectively. The ECG of a person in a normal condition (Fig. 11(a)) had two fractal dimensions of 1.02 and 1.80, respectively, in the small time scales separated at the scale s = 0.4 sec. The fact that the ECG was periodic was readily seen from the quasi-periodic behavior of the mIRMD results at large scales. From the intervals between two consecutive minima in Fig. 11(a), we found that its dominant period was about 0.85 sec. Furthermore, the fractal dimension of 1.64 at a large scale could be read from the increasing slope of the minimal values log(hG(s)i). Thus, according to the basic properties of mIRMD as described in Sec. VI, the normal ECG was a fractal of dimension 1.64 superimposed with a periodic trend of period 0.85 sec in the large scale. When VF

–1 0 1 2 ECG Signal (mV) 0 1 2 3 4 Time (sec)

FIG. 10: The ECG signal of a patient’s heartbeat changes pattern when ventricular fibrillation occurs. = 0.36 ∆ α = 0.98 α

s =424

s =90

FIG. 11: Log-log plot of hG(s)i versus sfor a typical patient’s ECG signal. Fractal dimensions at different scales were determined from the slopes of the interpolated lines shown. The unit of the scale s was 0.004 sec. (a) Results for normal condition. The dominant period at large scale was T = 424/2 × 0.004 = 0.85 sec. (b) Results when ventricular fibrillation occurred. The period of the signal was T = 90/2 × 0.004 = 0.18 sec.

occurred, the behavior of the ECG changed abruptly. As shown in Fig. 11(b), the ECG changed into a fractal with dimension 1.41 superposed with a periodic differentiable function of period 0.18 sec.

VIII. CONCLUSION

We have introduced a modified version of IRMD, called mIRMD, to determine the fractal dimensions of time sequences by calculating the midpoint displacements of the data. We have shown that both IRMD and mIRMD have satisfactory accuracy in determining the fractal dimensions of monofractals. IRMD calculated the midpoint displacements at scales limited to integral powers of 2. mIRMD was free from this constraint, and so revealed more details of the fractal spectra than the IRMD, which is important for scale-dependent real-time sequences. The mIRMD results were also more reliable, particularly for time sequences with short lengths of data, since all midpoint displacements at the same scale, rather than the non-overlapped ones as in the case of IRMD, were included in calculating the average. Of course, when computation time is a concern for cases with very long time sequences, IRMD is favored.

We have compared the results of mIRMD with those of the popular method of DFA for some monofractal functions and a typical real-time sequence with crossover-fractal struc-ture. Based on our comparison, we found that, qualitatively, mIRMD and DFA worked equally well in determining the scaling exponents and revealing the fractal structures of time sequences. However, the newly developed mIRMD method had one important fea-ture: its formulation was linear in the values of time sequences. This feature not only made calculations easier and faster than DFA, but also helped in the analyzing of scale-dependent fractal structures of non-stationary time sequences. In particular, we showed that an additive white noise or a differentiable trend was easily identified by the results of mIRMD.

We have applied mIRMD to two sets of real non-stationary time sequences. The results for the S&P500 stock index showed that during the period from 1986 to 1992 the S&P500 stock index had a fractal dimension close to 1 at scales less than about 25 minutes, and that the other fractal dimension approached 1.5, similar to a random walk at larger scales. The fractal spectrum of the ECG of a patient with ventricular fibrillation was more complex. The normal ECG had two distinct fractal dimensions in the small time scale region. At large scales there was a fractal superposed with a periodic trend. When ventricular fibrillation occurred, the fractal dimension of the ECG decreased. At small scales, the sequence changed from a fractal into a differentiable (but not twice-differentiable) curve. The period of the superposed trend at large scales also decreased at the VF stage.

Acknowledgments

This work was supported by grants from the National Science Council under grant number NSC97-2112-M005-001, and the National Center for Theoretical Sciences of Taiwan.

References

[1] B. B. Mandelbrot, Fractals, Form, Chance and Dimension (Freeman, San Francisco, 1977). [2] S.-S. Liaw and F.-Y. Chiu, Physica A388, 3100 (2009).

[3] L. Carpenter, Computer Graphics 14 (3), 109 (1980).

[4] N. Vandewalle and M. Ausloos, Phys. Rev. E 58, 6832 (1998). [5] C.-K. Peng et al., Phys. Rev. E49, 1685 (1994).

[6] K. Hu, P.Ch. Ivanov, Z. Chen, P. Carpena, and H. E. Stanley, Phys. Rev. E64, 011114 (2001). [7] J. W. Kantelhardt, E. Koscielby-Bunde, H. H. A. Rego, S. Havlin, and A. Bunde, Physica

A295, 441 (2001).

[8] A. V. Coronado and P. Carpena, J. Biol. Phys. 31, 121 (2005). [9] M. V. Berry and Z. V. Lewis, Proc. R. Soc. Lond. A370, 459 (1980).

[10] B. B. Mandelbrot, The Fractal Geometry of Nature, W. H. Freeman and Co., New York (1982). [11] A. Bashan, R. Bartsch, J.W. Kantelhardt, and S. Havlin, Physica A387, 5080 (2008).

[12] If we write the Weierstrass function as W (x) =P∞

n=02−n(2−D)cos(2

n_{πx +}π

2), we see that at x = 1/2m_{, all terms in the series with n > m vanish. Because of this, hG(s)i has many minima} at scales s = N (₂m11 −

1 2m2).

[13] C.-K. Peng, S. Havlin, H. E. Stanley, and A. L. Goldberger, CHAOS 5, 82 (1995).

[14] R. G. Kavasseri and R. Nagarajan, IEEE Trans. Circuits Syst., Part I: Fundamental Theory and Applications 51, 2255 (2004).

[15] Olsen Financial Technologies: http://www.olsendata.com/

[16] Y. Liu, P. Cizeau, M. Meyer, C.-K. Peng, and H. E. Stanley, Physica A245, 437 (1997). [17] T. Penzel, J. W. Kantelhardt, H. F. Becker, J. H. Peter, and A. Bunde, Comput. Cardiol. 30,

307 (2003).

[18] S. Havlin et al., Physica A274, 99 (1999).

[19] N. Scafetta, A. Ray, and B.J. West, Physica A359, 1 (2006). [20] J. W. Kantelhardt et al., Geophys. Res. 111, D01106 (2008).