Nothing
################################################################################
#' Detrended Fluctuation Analysis
#' @description
#' Functions for performing Detrended Fluctuation Analysis (DFA), a widely used
#' technique for detecting long range correlations in time series. These functions
#' are able to estimate several scaling exponents from the time series being analyzed.
#' These scaling exponents characterize short or long-term fluctuations, depending of
#' the range used for regression (see details).
#' @details The Detrended Fluctuation Analysis (DFA) has become a widely used
#' technique for detecting long range correlations in time series. The DFA procedure
#' may be summarized as follows:
#' \enumerate{
#' \item Integrate the time series to be analyzed. The time series resulting from the
#' integration will be referred to as the profile.
#' \item Divide the profile into N non-overlapping segments.
#' \item Calculate the local trend for each of the segments using least-square regression.
#' Compute the total error for each of the segments.
#' \item Compute the average of the total error over all segments and take its root square. By repeating
#' the previous steps for several segment sizes (let's denote it by t), we obtain the
#' so-called Fluction function \eqn{F(t)}.
#' \item If the data presents long-range power law correlations: \eqn{F(t) \sim t^\alpha}{F(t) proportional t^alpha} and
#' we may estimate using regression.
#' \item Usually, when plotting \eqn{\log(F(t))\;Vs\;log(t)}{log(F(t)) Vs log(t)} we may distinguish two linear regions.
#' By regression them separately, we obtain two scaling exponents, \emph{\eqn{\alpha_1}{alpha1}}
#' (characterizing short-term fluctuations) and \emph{\eqn{\alpha_2}{alpha2}} (characterizing long-term fluctuations).
#' }
#' Steps 1-4 are performed using the \emph{dfa} function. In order to obtain a estimate
#' of some scaling exponent, the user must use the \emph{estimate} function specifying
#' the regression range (window sizes used to detrend the series).
#' @param time.series The original time series to be analyzed.
#' @param npoints The number of different window sizes that will be used to estimate
#' the Fluctuation function in each zone.
#' @param window.size.range Range of values for the windows size that will be used
#' to estimate the fluctuation function. Default: c(10,300).
#' @param do.plot logical value. If TRUE (default value), a plot of the Fluctuation function is shown.
#' @return A \emph{dfa} object.
#' @author Constantino A. Garcia
#' @examples
#' \dontrun{
#' noise = rnorm(5000)
#' dfa.analysis = dfa(time.series = noise, npoints = 10,
#' window.size.range=c(10,1000), do.plot=FALSE)
#' cat("Theorical: 0.5---Estimated: ", estimate(dfa.analysis),"\n")
#' }
#' @rdname dfa
#' @export dfa
#' @exportClass dfa
dfa=function(time.series, window.size.range=c(10,300), npoints=20,do.plot=TRUE){
f="y~x"
window.sizes=c()
fluctuation.function=c()
#compute windows' sizes
log.window.sizes=seq(log10(window.size.range[[1]]),log10(window.size.range[[2]]),len=npoints)
window.sizes=10^log.window.sizes
# window.sizes must be an integer
window.sizes=unique(floor(window.sizes))
#integrate the time series time.series: this is called the "profile" in Penzel et al.
Y=calculateProfile(time.series)
#Calculate the local trend for
#each of the segments by a least-square fit
#of the time.series and compute the error
for (i in 1:npoints){
fluctuation.function[[i]]=fluctuationFunction(Y,f,window.sizes[[i]])
}
#create dfa object
dfa.structure = list(window.sizes=window.sizes,fluctuation.function=fluctuation.function)
class(dfa.structure) = "dfa"
#plot
if (do.plot){
plot(dfa.structure)
}
return(dfa.structure)
}
#' @return The \emph{getWindowSizes} function returns the windows sizes used
#' to detrend the time series
#' @rdname dfa
#' @export getWindowSizes
#'
getWindowSizes= function(x){
return (x$window.sizes)
}
#' @return The \emph{getFluctuationFunction} function returns the Fluctuation function
#' of the DFA
#' @rdname dfa
#' @export getFluctuationFunction
getFluctuationFunction= function(x){
return (x$fluctuation.function)
}
#' @rdname dfa
#' @method plot dfa
#' @param ... Additional parameters.
#' @S3method plot dfa
#'
plot.dfa = function(x, ...){
plot(x$window.sizes,x$fluctuation.function,log="xy",main="Detrended fluctuation analysis\n",xlab="Window size: t",ylab="Fluctuation function: F(t)")
}
#' @rdname dfa
#' @method estimate dfa
#' @param x A \emph{dfa} object.
#' @param regression.range Vector with 2 components denoting the range where the function will perform linear regression.
#' @S3method estimate dfa
#'
estimate.dfa = function(x, regression.range = NULL,do.plot=FALSE,...){
log.fluctuation.function = log10(x$fluctuation.function)
log.window.sizes = log10(x$window.sizes)
if (is.null(regression.range)){ regression.range = range(x$window.sizes)}
#prepare interpolations to calculate the exponents
index.alpha = which(x$window.sizes>=regression.range[[1]]&x$window.sizes<=regression.range[[2]])
#get alpha
interp1=lm(log.fluctuation.function[index.alpha]~log.window.sizes[index.alpha])
alpha=interp1$coefficients[[2]]
if (do.plot){
plot(x)
lines(x$window.sizes[index.alpha],10^interp1$fitted.values,col="blue")
legend("bottomright",col="blue",lty=c(1),lwd=c(2.5), legend=paste("Scaling exponent estimate =", sprintf("%.2f",alpha)))
}
return(alpha)
}
# private
# Calculate the local trend for
# each of the segments by a least-square fit
# of the time.series and compute the error
# Y: the profile
# fint: interpolation formula
# l: size of the window
fluctuationFunction=function(Y,fint,l){
#preliminary definitions
ldata=length(Y)
ngroups=floor(ldata/l)
x=1:l
# Divide the profile in ngroups non-overlaping groups. Interpolate and calculate the residual error
F2=c()
for (i in 1:ngroups){
beg=(i-1)*l+1; en=i*l;
y=Y[beg:en]
#interpolate
fit=lm(as.formula(fint))
#get error
F2[i]=sum(fit$residuals^2)/l
}
return (sqrt(sum(F2)/ngroups))
}
# private
# Determine the “profile” of the timeseries time.series
calculateProfile=function(time.series){
lendata=length(time.series)
time.series=time.series-mean(time.series)
#profile
Y=c(time.series[1])
for (i in (2:lendata)){
Y[i]=time.series[i]+Y[i-1]
}
return(Y)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.