Nothing
######################################################
## FRACTAL fractionally differenced process functions
##
## FDSimulate
## FDWhittle
##
######################################################
###
# FDSimulate
##
"FDSimulate" <- function(delta, innovations.var=1, method="ce", seed=0)
{
# define local functions
"FDCirculantEmbedding" <- function(delta, innovations.var=1, n.sample=16, seed=0)
{
# check input parameters
checkScalarType(delta,"numeric")
checkScalarType(innovations.var,"numeric")
checkScalarType(n.sample,"integer")
checkScalarType(seed,"integer")
seed <- as.integer(abs(seed))
if (n.sample < 2)
stop("n.sample must be an integer greater than unity")
# initialize variables
ncumsum <- max((2 * delta - 1) %/% 2 + 1, 0)
delta <- delta - ncumsum
# form fdp avcs
fdp.acvs <- lmACF(lmModel("fdp",delta=delta, innov=innovations.var),
lag.max=n.sample, type="covariance")@data
# concatenate reversed fdp acvs
fdp.acvs <- c(fdp.acvs, rev(fdp.acvs[2:n.sample]))
# take DFT
S <- Re(fft(fdp.acvs, inverse=FALSE))
# use CE method create bootstrap surrogate
z <- as.vector(.Call("RS_fractal_bootstrap_circulant_embedding",
S, seed,
COPY=rep(FALSE,2), CLASSES=c("matrix","integer"),
PACKAGE="ifultools"))[seq(n.sample)]
if (!ncumsum)
return(z)
for (i in seq(ncumsum))
z <- cumsum(z)
z
}
# check input argument types and lengths
checkVectorType(delta,"numeric")
checkScalarType(method,"character")
if (!is.numeric(innovations.var))
stop("innovations.var must be a numeric object")
checkScalarType(seed,"integer")
seed <- as.integer(abs(seed))
# initialize variables
n.innov <- length(innovations.var)
n.delta <- length(delta)
# perform length check
if (n.delta > n.innov)
innovations.var <- c(innovations.var, rep(innovations.var[n.innov],
n.delta - n.innov))
# check method
method <- match.arg(method, c("ce","cholesky"))
#TODO: axe this probably: we may only support ce method
if (method == "cholesky"){
z <- as.vector(.Call("RS_fractal_fdp_simulate", delta, innovations.var,
COPY = rep(FALSE,2), CLASSES= c("numeric", "numeric"),
PACKAGE="ifultools"))
}
else{
# find unique delta levels and match against
# original delta vector
delta.levels <- sort(unique(delta))
idelta <- match(delta, delta.levels)
# create matrix of FD realizations, each column
# corresponds to a unique delta
z <- apply(matrix(delta.levels), MARGIN=1, function(x, n.delta, seed, FDCirculantEmbedding)
FDCirculantEmbedding(delta=x, innovations.var=1, n.sample=n.delta),
n.delta=n.delta, seed=seed, FDCirculantEmbedding=FDCirculantEmbedding)
z <- z[cbind(seq(n.delta), idelta)] * innovations.var
}
oldClass(z) <- "FDSimulate"
attr(z, "delta") <- delta
attr(z, "innov") <- innovations.var
attr(z, "method") <- "circulant embedding"
z
}
##
# FDWhittle
##
"FDWhittle" <- function(x, method="continuous", dc=FALSE,
freq.max=0.5, delta.min=-1,delta.max=2.5, sdf.method="direct", ...)
{
# define local functions
Isumd <- function(delta, Sx, discrete)
{
# function to evaluate the weighted sum of an FD process spectrum
f <- Sx$frequency
N <- Sx$n.sample
Np <- (N - 1) %/% 2
Q <- 2 * mean(Sx$data * (4 * (sin(pi * f))^2)^delta)
ifelse1(!discrete, Q, log(Q) + 2 * delta * (((N/(2 * Np)) - 1) *
log(2) - log(2 * N)/(2 * Np)))
}
# check input arguments
checkScalarType(method,"character")
method <- match.arg(lowerCase(method),c("discrete","continuous"))
checkScalarType(freq.max,"numeric")
if (freq.max <= 0.0 || freq.max > 0.5)
stop("freq.max must be on the normalized frequency range (0,0.5)")
# calculate single-sided SDF estimate with a normalized
# spectral resolution of 1/N
sdf <- SDF(x, method=sdf.method, single.sided=TRUE, npad=length(x), ...)
# associate freq.max normalized frequency
# with corresponding index into SDF
# frequency vector
xatt <- attributes(sdf)
freq.indx.max <- ceiling(freq.max/xatt$deltaf)
freq.indx.min <- ifelse1(dc, 1, 2)
ifreq <- seq(freq.indx.min, freq.indx.max)
# pack relevant spectral information into a list.
# limit the range of sdf and corresponding frequency
# according to the input specifications
Sx <- list(data=sdf[ifreq], frequency=xatt$frequency[ifreq], n.sample=xatt$n.sample)
# use optimize() function to find the best value of delta
# return Hurst coefficient estimate
#dfinal$minimum + 0.5
z <- optimize(Isumd, lower=delta.min, upper=delta.max, tol=0.001,
maximum=FALSE, Sx=Sx, discrete=is.element(method,"discrete"))$minimum
names(z) <- "delta"
z
}
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.