Nothing
# private function
# Estimates an appropiate number of boxes for using in the box assisted
# algorithm using the time series being analyzed and the radius of the grid.
# @details
# The estimation takes the difference between the maximum and the minimum of the
# time series and divides it by the radius of the grid. If the number of boxes
# is too large (over a kMaximumNumberBoxes), the number of boxes is fixed to
# kMaximumNumberBoxes. If the number of boxes is too small (below
# kMinimumNumberBoxes), the number of boxes is fixed to kMinimumNumberBoxes.
# @param time.series the original time.series from which the surrogate data is
# generated
# @param radius width of each of the rectangles that will be used to build the
# grid in the box assisted algorithm
# @return the number of boxes to used in the box assisted algorithm
# @author Constantino A. Garcia
estimateNumberBoxes = function(time.series, radius){
kMinimumNumberBoxes = 10
kMaximumNumberBoxes = 500
number.boxes = as.integer((max(time.series) - min(time.series)) / radius)
if (number.boxes > kMaximumNumberBoxes) number.boxes = kMaximumNumberBoxes
if (number.boxes < kMinimumNumberBoxes) number.boxes = kMinimumNumberBoxes
number.boxes
}
# private function
# checks if takens' vectors v1 and v2 are neighbours
# the neighbourhood is determined using the max norm and an radius radious.
# embedding.dim is not strictly necessary but it is used to accelerate the
# computations
isNeighbour = function(v1,v2,embedding.dim,radius){
for (i in 1:embedding.dim) {
if (abs(v1[[i]] - v2[[i]]) >= radius) {
return(FALSE)
}
}
TRUE
}
#' Build the Takens' vectors
#' @description
#' This function builds the Takens' vectors from a given time series. The set
#' of Takens' vectors is the result of embedding the time series in a
#' m-dimensional space. That is, the \eqn{n^{th}} Takens' vector is defined as
#' \deqn{T[n]=\{time.series[n], time.series[n+ timeLag],...time.series[n+m*timeLag]\}.}
#' Taken's theorem states that we can then reconstruct an equivalent dynamical
#' system to the original one (the
#' dynamical system that generated the observed time series) by using the
#' Takens' vectors.
#' @param time.series The original time series.
#' @param embedding.dim Integer denoting the dimension in which we shall embed
#' the time.series.
#' @param time.lag Integer denoting the number of time steps that will be use
#' to construct the Takens' vectors.
#' @return A matrix containing the Takens' vectors (one per row). The resulting
#' matrix also contains information about the time lag and the embedding
#' dimension used (as attributes).
#' @references H. Kantz and T. Schreiber: Nonlinear Time series Analysis
#' (Cambridge university press)
#' @author Constantino A. Garcia and Gunther Sawitzki.
#' @examples
#' \dontrun{
#'# Build the Takens vector for the Henon map using the x-coordinate time series
#' h = henon(n.sample= 3000,n.transient= 100, a = 1.4, b = 0.3,
#' start = c(0.73954883, 0.04772637), do.plot = FALSE)
#' takens = buildTakens(h$x,embedding.dim=2,time.lag=1)
#' # using the x-coordinate time series we are able to reconstruct
#' # the state space of the Henon map
#' plot(takens)}
#' @export buildTakens
buildTakens = function(time.series, embedding.dim, time.lag) {
N = length(time.series)
#vector of increments
maxjump = (embedding.dim - 1) * time.lag
jumpsvect = seq(0, maxjump, time.lag)
#matrix that will store the takens' vectos. One vector per row
numelem = N - maxjump
takens = matrix(nrow = numelem, ncol = embedding.dim)
#build takens' vectors
for (i in 1:numelem) {
takens[i, 1:embedding.dim] = time.series[jumpsvect + i]
}
#
#class(takens) = "takens"
id = deparse(substitute(time.series))
attr(takens,"embedding.dim") = embedding.dim
attr(takens,"time.lag") = time.lag
attr(takens,"id") = id
takens
}
#' Estimate an appropiate time lag for the Takens' vectors
#' @description
#' Given a time series (time.series), an embedding dimension (m) and a
#' time lag (timeLag), the \eqn{n^{th}}
#' Takens' vector is defined as
#' \deqn{T[n]=\{time.series[n], time.series[n+ timeLag],...time.series[n+m*timeLag]\}.}
#' This function estimates an appropiate time lag by using the autocorrelation
#' function or the average mutual information .
#' @details
#' A basic criteria for estimating a proper time lag is based on the following
#' reasoning: if the time lag used to build the Takens' vectors is too small,
#' the coordinates will be too highly temporally correlated and the embedding
#' will tend to cluster around the diagonal in the phase space. If the time lag
#' is chosen too large, the resulting coordinates may be almost uncorrelated and
#' the resulting embedding will be very complicated. Thus, the autocorrelation
#' function can be used for estimating an appropiate time lag of a time series.
#' However, it must be noted that the autocorrelation is a linear statistic,
#' and thus it does not take into account nonlinear dynamical correlations. To
#' take into account nonlinear correlations the average mutual information (AMI)
#' can be used. Independently of the technique used to compute the correlation,
#' the time lag can be selected in a variety of ways:
#' \itemize{
#' \item Select the time lag where the autocorrelation/AMI function decays to 0
#' (\emph{first.zero} selection method). This
#' method is not appropriate for the AMI function, since it only takes
#' positive values.
#' \item Select the time lag where the autocorrelation/AMI function decays to
#' 1/e of its value at zero (\emph{first.e.decay} selection method).
#' \item Select the time lag where the autocorrelation/AMI function reaches
#' its first minimum (\emph{first.minimum} selection method).
#' \item Select the time lag where the autocorrelation/AMI function decays to
#' the value specified by the user (\emph{first.value} selection method and
#' \emph{value} parameter).
#' }
#' @param time.series The original time series.
#' @param technique The technique that we shall use to estimate the time lag
#' (see the Details section). Allowed values are \emph{"acf"} and \emph{"ami"}.
#' @param selection.method Method used for selecting a concrete time lag.
#' Available methods are \emph{"first.zero"}, \emph{"first.e.decay"} (default),
#' \emph{"first.minimum"} and \emph{"first.value"}.
#' @param value Numeric value indicating the value that the autocorrelation/AMI
#' function must cross in order to select the time lag. It is used only with
#' the "first.value" selection method.
#' @param lag.max Maximum lag at which to calculate the acf/AMI.
#' @param do.plot Logical value. If TRUE (default value), a plot of the
#' autocorrelation/AMI function is shown.
#' @param main A title for the plot.
#' @param ... Additional parameters for the \emph{acf} or the
#' \emph{mutualInformation}.
#' @return The estimated time lag.
#' @examples
#' \dontrun{
#' sx = sinaiMap(a=0.3,n.sample=5000,start=c(0.23489,0.8923),do.plot=FALSE)$x
#' timeLag(sx, technique="ami",
#' n.partitions = 20, units = "Bits")
#' timeLag(sx, technique="acf") }
#' @note If the autocorrelation/AMI function does not cross the specifiged value,
#' an error is thrown. This may be solved by increasing the \emph{lag.max} or
#' selecting a higher value to which the autocorrelation/AMI function may decay.
#' @seealso \code{\link{mutualInformation}}
#' @references H. Kantz and T. Schreiber: Nonlinear Time series Analysis
#' (Cambridge university press)
#' @author Constantino A. Garcia
#' @export timeLag
timeLag = function(time.series, technique=c("acf","ami"),
selection.method=c("first.e.decay", "first.zero",
"first.minimum", "first.value"),
value = 1/exp(1), lag.max=NULL, do.plot=TRUE,
main = NULL, ...) {
technique = match.arg(technique)
if (is.null(main)) {
main = switch(technique,
acf = "Autocorrelation function",
ami = "Average Mutual Information (AMI)")
}
if (is.null(lag.max) && (technique == "acf")) {
lag.max = max(20, length(time.series)/2)
# if technique=="ami", mutualInformation function will handle the NULL value
}
fx = switch(technique,
ami = mutualInformation(
time.series, lag.max = lag.max,
do.plot = do.plot, main = main, ...
)$mutual.information,
acf = stats::acf(
time.series, lag.max = lag.max, plot = do.plot,
type = "correlation", main = main, ...
)$acf
)
selection.method = match.arg(selection.method)
lag = switch(match.arg(selection.method),
first.e.decay = get.position.decays(fx,fx[1]/exp(1),technique),
first.zero = get.position.decays(fx,0,technique),
first.minimum = get.first.minimum(fx,technique),
first.value = get.position.decays(fx,value,technique) )
if (do.plot) {
if (selection.method == "first.e.decay") {
abline(h = fx[1] / exp(1), lty = 3, col = 3, lwd = 2)
}else if (selection.method == "first.value") {
abline(h = value, lty = 3, col = 3, lwd = 2)
}
}
lag
}
# internal function. It will calculate the first time lag where
# the vector x decays to a given value
get.position.decays = function(x, value, method ) {
cross.position = which(x <= value)
# if x does not cross the value specified
# by the user, a warning is given
if (length(cross.position) == 0) {
stop(paste("The ",method," function does not cross ", value,
". Choose another \"cross\" value!\n",sep = ""))
} else{
cross.position = cross.position[[1]]
}
# If possible: convert from positions to time lags by substracting 1
if ((cross.position) > 1) {
cross.position = cross.position - 1
}
cross.position
}
# internal function. It will calculate the first time lag where
# the x vector has a minimum
get.first.minimum = function(x,method){
derivative = diff(x)
derivative = ifelse( derivative < 0, -1, derivative)
derivative = ifelse( derivative == 0, 0, derivative)
derivative = ifelse( derivative > 0, 1, derivative)
minimums = which(diff(derivative) >= 1)
# if the autocorrelation function is monotonically decreasing, an error is
# given
if (length(minimums) == 0) {
stop(paste("The ", method, " function does not have a minimum\n",sep = ""))
}else{
#take into account that we have loss some samples after diff
first.minimum = minimums[[1]] + 1
}
# convert from positions to time lags by substracting 1
(first.minimum - 1)
}
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.