Fractals and Stochastic Calculus

Just another WordPress.com weblog

R/s Analysis to estimate the Hurst exponent

R/s analysis (or Rescaled Range analysis) was originally devised by Harold Edwin Hurst in its studies of the Nile discharge in 1951. It is a rather simple method, easily implemented in a program and that provides a direct estimation of the Hurst Exponent which is a precious indicator of the state of randomness of a time-series. It is especially interesting in revealing the existence of long-term dependence, which prevents, when it exists, the time-series to be reasonably modelised by a random walk.

I wish to expose here the minimal information for understanding the method, I will also provide a few references for those who wish to deepen their understanding of the matter.

Given a time-series with n elements X_1, X_2,...,X_n , the R/s statistic is defined as:

R/s(n)= \frac{1}{s}[Max_k{\sum_{i=1}^{k}(x_i-\bar{x})}-Min_k{\sum_{i=1}^{k}(x_i-\bar{x})}]

Where 1\leq k \leq n,
\bar{x} is the arithmetic mean
And s=\sqrt{\frac{1}{n}\sum(x_i-\bar{x})^2} is the standard deviation from the mean.

With this R/s value, Hurst found a generalization of a result found by Einstein in 1905 (Investigations on the Theory of the Brownian Movement ) as equation (11) (in the cited paper) in the following formula:

E[R/s(n)]=Cn^H \text{ as }n\rightarrow\infty \,\,\,\,(1)

Where H is the Hurst exponent.
From there, it is clear that we can obtain an estimation of the Hurst exponent pretty easily from an R/s analysis.
Several sites and articles propose a detailed methodology to implement R/s Analysis, I primarily use the approach exposed in a paper by O. Rose from February 1996:
Estimation of the Hurst Parameter of Long-Range Dependent Time Series
With a slight difference, however, I shall only plot one value of R/s for each value of d, in a manner similar to the following site: Estimating the Hurst Exponent
Anyway, I feel both articles are not very clear in their notations, and I therefore will detail the analysis I wish to implement.

I- RESCALED RANGE ANALYSIS

Considering the time series above X_j (j=1,..,N)
We divide the time series into K^u (*) non-overlapping blocks of length d^u=\frac{N}{K^u}
And we fix: t_i=d^u(i-1)+1

Next we get a new time series W(i,k):

W(i,k)=\sum_{j=1}^{k}{[X_{t_i+j-1}-\frac{1}{d^u}\sum_{v=1}^{d^u}{X_{t_i+v-1}}]},\,k=1,..,d^u

From there, we get the following rescaled range:
R/s(i,u)=\frac{R(i,d^u)}{s(i,d^u)}

With:

R(i,d^u)=Max\{0,W(i,1),...,W(i,d^u)\}-Min\{0,W(i,1),...,W(i,d^u)\}
And
s(i,d^u)=\sqrt{\frac{1}{d^u}\sum_{j=1}^{d^u}[X_{t_i+j-1}-\frac{1}{d^u}\sum_{v=1}^{d^u}X_{t_i+v-1}]^2}

Taking the mean over i , we then get R/s(d^u):

R/s(d^u)=\frac{1}{K^u}\sum_{i=1}^{K^u}R/s(i,d^u)

Considering equation (1):            log(R/s(d^u))=log(C) +Hlog(d^u)

We can plot log(R/s(d^u)) vs log(d^u) for u varying, H is then the slope of the regression line which we simply get from the linear least squares method.
Fixing:
x_u=log(d^u) and y_u=log(R/s(d^u))

We get:

\boxed{H=\frac{U\sum_{u}x_uy_u -(\sum_{u}x_u)(\sum_{u}y_u)}{U(\sum_{u}x_{u}^2)-(\sum_{u}x_u)^2}} \text{ with }u \text{ varying from }1 \text{ to } U

(*): N and K^u are chosen adequately so that d^u is always an integer.

II- IMPLEMENTATION

I implemented this method within an indicator for Metatrader 4 in order to compute a Fractalised Simple Moving Average for FOREX fluctuations. The implementation is not very interesting, partly because the computation time on this platform are not very good, I was therefore constraint to use the R/s method on a very limited number of data.
Anyway, the implementation can be seen on my other blog: Rescaled Range Analysis

October 14, 2009 Posted by | Uncategorized | Leave a comment

Standard Deviation of Fractional Brownian Motion

Let B_{H}(t) be the zero-mean Fractional Brownian Motion process of Hurst parameter H (0 \le H \le 1), such a process is defined as having the following covariance structure:

E[B_{H}(t)B_{H}(s)]=\frac{1}{2}(|t|^{2H}+|s|^{2H}-|t-s|^{2H}) (1)

Then(by taking s=t in (1))

Var(B_{H}(t))=E[(B_{H}(t)-E(B_{H}(t)))^2]=E[B_{H}^{2}(t)]=|t|^{2H}

Which gives the standard deviation as |t|^H (and for the Wiener Brownian Motion, we indeed get a standard deviation of \sqrt{t})

May 5, 2009 Posted by | Fractional Brownian Motion | Leave a comment

Maximum of Wiener Brownian Motion

Let W_0 be the standard Wiener Brownian Motion process, and M(t) defined as M(t)=\max_{0 \leq s \leq t} W_{0}(t), we then have the following result:

P(M(t) \geq x)= 2(1-\phi(\frac{x}{\sqrt{t}}))

where \phi denotes the standard normal cumulative distribution function:

\phi(z)=\frac{1}{\sqrt{2\pi}}e^{-z^{2}/2} \forall{z}\in\mathbb{R}

This result is simply a direct application of Theorem 2 from this paper with \alpha=0

May 4, 2009 Posted by | Wiener Brownian Motion | Leave a comment