Nothing
################################################
# FRACTAL Hurst coefficient estimators and
# related statistical functions
#
# Functions:
#
# dispersion
# DFA
# HDEst
# hurstACVF
# hurstBlock :: aggAbs, aggVar, diffVar, Higu
# hurstSpec :: pgramReg, pgramRegMod, RobInt
# RoverS
#
################################################
##
# dispersion
##
"dispersion" <- function(x, front=FALSE)
{
checkScalarType(front,"logical")
n.sample <- length(x)
if (n.sample < 32)
stop("Time series must contain at least 32 points")
scale <- as.integer(logScale(scale.min=1, scale.max=n.sample, scale.ratio=2))
sd <- unlist(lapply(scale, function(m, x, n.sample, front){
nscale <- n.sample %/% m
dyad <- nscale * m
n.offset <- ifelse1(front, n.sample - dyad, 0)
stdev(aggregateData(x[seq(1 + n.offset, length=dyad)],
by=m, FUN=mean))
},
x=x,
n.sample=n.sample,
front=front))
list(sd=sd, scale=scale)
}
##
# DFA
##
"DFA" <- function(x, detrend="poly1", sum.order=0, overlap=0,
scale.max=trunc(length(x)/2), scale.min=NULL, scale.ratio=2,
verbose=FALSE)
{
# define local functions
"polyfit.model" <- function(polyfit.order){
if (polyfit.order < 0)
stop("Polynomial fit order must be positive")
if (polyfit.order == 0)
return("x ~ 1")
if (polyfit.order == 1)
return("x ~ 1 + t")
else{
return(paste(c("x ~ 1", "t", paste("I(t^", seq(2, polyfit.order), ")", collapse = " + ", sep = "")),
collapse=" + ", sep = ""))
}
}
"regression.poly" <- function(x, model=polyfit.model(1)){
order <- length(attr(terms(formula(x~1)),"order"))
t <- x@positions
x <- x@data
if (order == 0)
return(sum((x - mean(x))^2))
if (order == 1){
x <- x - mean(x)
t <- t - mean(t)
slope <- sum(t*x)/sum(t^2)
return (sum((x - slope*t)^2))
}
fit <- lm(model, data=data.frame(list(t=t, x=x)))
return(sum(fit$residuals^2))
}
"regression.bridge" <- function(x, ...){
x <- x@data
N <- length(x)
bridge <- seq(x[1], x[N], length=N)
x <- x - bridge
return(sum(x^2))
}
"regression.none" <- function(x, ...)
return(sum(x@data^2))
# check input argument types and lengths
checkScalarType(detrend,"character")
checkScalarType(scale.ratio,"numeric")
checkScalarType(verbose,"logical")
checkScalarType(overlap,"numeric")
# obtain series name
data.name <- deparse(substitute(x))
# convert sequence to signalSeries class
x <- wmtsa::create.signalSeries(x)
# check overlap argument
if ((overlap < 0) | (overlap >= 1))
stop("Overlap factor must be in the range [0,1)")
# map the detrending method
detrend <- lowerCase(detrend)
if (substring(detrend,1,4) == "poly"){
polyfit.order <- as.integer(substring(detrend,5))
if (is.na(polyfit.order)){
stop("Improperly formed detrending string")
}
else if (polyfit.order < 0){
stop("Polynomial fit order must be positive")
}
# form detrending model
model <- formula(polyfit.model(polyfit.order))
if (verbose) {
cat("Detrending model: ");
print(model)
}
# set regression function
regressor <- regression.poly
modstr <- as.character(model)[2:3]
regress.str <- paste(modstr,collapse=" ~ ")
}
else if (charmatch(detrend, "bridge", nomatch=FALSE)){
regressor <- regression.bridge
regress.str <- "Bridge detrended"
}
else if (charmatch(detrend, "none", nomatch=FALSE)){
regressor <- regression.none
regress.str <- "None"
}
else stop("Detrending method is not supported")
# coerce orders to integers
sum.order <- trunc(sum.order)
# perform cumulative summation(s) or difference(s)
if (sum.order > 0){
for (i in seq(sum.order))
x <- cumsum(x)
}
else if (sum.order < 0){
for (i in seq(sum.order))
x <- diff(x)
}
# obtain sequence length (after possibly
# taking differences)
N <- length(x)
# create scale for data analysis
if (is.null(scale.min))
scale.min <- ifelse1(substring(detrend,1,4) == "poly", 2*(polyfit.order + 1), min(N/4, 4))
checkScalarType(scale.min, "numeric")
checkScalarType(scale.max, "numeric")
# DFA requires integer scales
scale.min <- trunc(scale.min)
scale.max <- trunc(scale.max)
if (scale.min > scale.max)
stop("Scale minimum cannot exceed scale maximum")
if (any(c(scale.min, scale.max) < 0))
stop("Scale minimum and maximum must be positive")
if (any(c(scale.min, scale.max) > N))
stop(paste("Scale minimum and maximum must not exceed length",
"of (differenced/cummulatively summed) time series"))
# initialize output
scale <- logScale(scale.min, scale.max, scale.ratio=scale.ratio, coerce=trunc)
scale <- scale[scale > 1]
rmse <- rep(0.0, length(scale))
# calculate DFA
for (i in seq(along=scale)){
sc <- scale[i]
noverlap <- trunc(overlap * sc)
if (verbose)
cat("Processing scale", sc,"\n")
cumsum <- 0.0
count <- 0
index <- seq(sc)
# sum the mse over all blocks
while(index[sc] <= N){
rmse[i] <- rmse[i] + regressor(x[index], model=model)
count <- count + 1
index <- index + sc - noverlap
}
# normalize by the number points
# in all blocks. this normalization
# factor may differ from the number
# of points in the time series (N)
# since N may not be evenly divisible
# by the current block size. finally,
# take the sqrt() to form the root
# mean square error
rmse[i] <- sqrt(rmse[i] / (count * sc))
}
# fit the log-log data
# weight the coefficients starting from 1 down to 1/N
xx <- log(scale)
yy <- log(rmse)
ww <- 1/seq(along=xx)
w <- NULL # quell R CMD check scoping issue
logfit <- lm(y ~ 1 + x, data=data.frame(x=xx, y=yy, w=ww), weights=w)
exponent <- logfit$coefficients["x"]
# return the result as a fractal block exponent
fractalBlock(
domain ="time",
estimator ="Detrended Fluctuation Analysis",
exponent =exponent,
exponent.name="H",
scale =scale,
stat =rmse,
stat.name ="RMSE",
detrend =regress.str,
overlap =overlap,
data.name =data.name,
sum.order =sum.order,
series =asVector(x),
logfit =logfit)
}
##
# HDEst
##
"HDEst" <- function(NFT, sdf, A=0.3, delta=6/7)
{
checkScalarType(NFT,"integer")
checkScalarType(A,"numeric")
checkScalarType(delta,"numeric")
# define frequencies
omega <- ((2 * pi)/NFT) * (1:NFT)
# notation as in reference: construct X matrix
L <- round(A * NFT^delta)
c1 <- rep(1, L)
c2 <- log(abs(2 * sin(omega[1:L]/2)))
c3 <- omega[1:L]^2/2
X <- cbind(c1, c2, c3)
# construct least squares product matrix and find 3rd row:
TX <- t(X)
ProdX <- solve(TX %*% X) %*% TX
b <- ProdX[3,]
# evaluate K.hat, C.hat, and optimum m, as in Reference;
# use first L non-zero frequencies:
K.hat <- as.numeric(b %*% log(sdf[1:L]))
if (K.hat == 0)
stop("Method fails: K.hat=0")
C.hat <- ((27/(128 * pi^2))^0.2) * (K.hat^2)^-0.2
round(C.hat * NFT^0.8)
}
##
# hurstACVF
##
"hurstACVF" <- function(x, Ascale=1000, lag.min=1, lag.max=64)
{
# check input argument types and lengths
checkScalarType(Ascale,"numeric")
checkScalarType(lag.min,"integer")
checkScalarType(lag.max,"integer")
# estimates long memory parameters by linear regression of scaled asinh
# of ACVF vs log(lag) over intermediate lag values.
x <- wmtsa::create.signalSeries(x)@data
n.sample <- length(x)
# evaluate autocovariance function (ACVF)
acfx <- acf(x, n.sample, "covariance")$acf
# include lag=0 point
yy <- c(1, acfx)
# linear regression over lag.min to lag.max
fitcoef <- lsfit(log(lag.min:lag.max), asinh(Ascale * yy[(lag.min + 1):
(lag.max + 1)])/asinh(Ascale))$coef
# relate slope to long memory parameters (beta=slope of ACVF,
# alpha=parameter in corresponing PPL model, HG=corresponding
# Hurst parameter for stationary FGN model.
beta <- fitcoef[2] * asinh(Ascale)
alpha <- -1 - beta
HG <- 1 + beta/2
drop(list(beta=beta, alpha=alpha, HG=HG))
}
##
# hurstBlock
##
"hurstBlock" <- function(x, method="aggAbs", scale.min=8, scale.max=NULL,
scale.ratio=2, weight=function(x) rep(1,length(x)), fit=lm)
{
# obtain series name
data.name <- deparse(substitute(x))
# convert input to signalSeries class
x <- wmtsa::create.signalSeries(x)
# initialize variables
n.sample <- length(x)
if (is.null(scale.max))
scale.max <- n.sample
# check inputs
checkScalarType(method, "character")
method <- match.arg(lowerCase(method), c("aggabs","aggvar","diffvar","higuchi"))
checkScalarType(scale.ratio, "integer")
if (scale.ratio < 2)
stop("scale.ratio must be greater than one")
checkScalarType(scale.min, "integer")
if (scale.min < 2)
stop("scale.min must be greater than one")
checkScalarType(scale.max, "integer")
if (scale.max > n.sample)
stop("scale.max cannot exceed series length")
if (!is.function(weight))
stop("weight must be a function")
if (!is.function(fit))
stop("fit must be a linear regression function")
fitstr <- deparse(substitute(fit))
supported <- c("lm","lmsreg","ltsreg")
if (!is.element(fitstr,supported ))
stop("Supported linear regression functions are:", supported)
lmfit <- is.element(fitstr,"lm")
# form the scale vector
scale <- logScale(scale.min, scale.max, scale.ratio=scale.ratio, coerce=trunc)
# define block statistic
statfunc <- ifelse1(method == "aggabs",
function(x) mean(abs(x)),
function(x) variance(x, unbiased=FALSE))
statistic <- unlist(lapply(scale,
function(blockWidth, x, statfunc, n.sample, higuchi){
if (higuchi){
return(mean(abs(aggregateData(x, by=1, FUN=mean, moving=blockWidth))) *
((n.sample - 1) / blockWidth))
}
nBlocks <- n.sample %/% blockWidth
nUpper <- nBlocks * blockWidth
xcut <- as.numeric(x[1:nUpper])
# reshape recentered time series as a matrix
tsMat <- matrix(xcut - mean(xcut), blockWidth, nBlocks)
# evaluate the statistic of each column's/block's mean value
statfunc(colMeans(tsMat))
},
x=x@data, statfunc=statfunc, n.sample=n.sample,
higuchi=is.element(method,"higuchi")))
# set up for linear regression of variability vs group size:
if (is.element(method, c("aggabs","aggvar","higuchi"))){
xvar <- scale
yvar <- statistic
}
else if (is.element(method,"diffvar")){
# take first difference, and exclude any points
# with non-positive differences
dStat <- -diff(statistic)
xvar <- scale[dStat > 0][-length(scale)]
yvar <- dStat[dStat > 0]
}
w <- NULL # quell R CMD check scoping issue
# fit the log-log data weight the coefficients
logfit <- ifelse1(lmfit,
fit(y ~ 1 + x, data=data.frame(x=log(xvar), y=log(yvar), w=weight(yvar)), weights=w),
fit(y ~ 1 + x, data=data.frame(x=log(xvar), y=log(yvar))))
exponent <- logfit$coefficients["x"]
# estimate H from the slope of the regression fit:
H <- switch(method,
aggabs =1 + exponent,
aggvar =1 + exponent/2,
diffvar=1 + exponent/2,
higuchi=exponent + 2)
# return the result as a fractal block exponent
fractalBlock(
domain ="time",
estimator ="Hurst coefficient via Aggregated Series",
exponent =H,
exponent.name="H",
scale =xvar,
stat =yvar,
stat.name =method,
detrend =NULL,
overlap =0,
data.name =data.name,
sum.order =0,
series =asVector(x),
logfit =logfit)
}
##
# hurstSpec
##
"hurstSpec" <- function(x, method="standard", freq.max=0.25, dc=FALSE, n.block=NULL,
weight=function(x) rep(1,length(x)), fit=lm, sdf.method="direct", ...)
{
# check input arguments
checkScalarType(method,"character")
method <- match.arg(lowerCase(method),c("standard","smoothed","robinson"))
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)")
if (!is.function(fit))
stop("fit must be a linear regression function")
fitstr <- deparse(substitute(fit))
supported <- c("lm","lmsreg","ltsreg")
if (!is.element(fitstr,supported ))
stop("Supported linear regression functions are:", supported)
lmfit <- is.element(fitstr,"lm")
# obtain series name
data.name <- deparse(substitute(x))
# if SDF smoothing is requested, force removal of DC component
if (is.element(method,"smoothed"))
dc <- FALSE
# calculate single-sided SDF estimate with a normalized
# spectral resolution of 1/N
sdf <- sapa::SDF(x, method=sdf.method, single.sided=TRUE, npad=length(x), ...)
attr(sdf,"series.name") <- data.name
# associate freq.max normalized frequency
# with corresponding index into SDF
# frequency vector
df <- attr(sdf, "deltaf")
freq.indx.max <- ceiling(freq.max/df)
freq.indx.min <- ifelse1(dc, 1, 2)
if (is.element(method,c("standard","smoothed"))){
# set up variables for linear regression:
xvar <- log10(freq.indx.min:freq.indx.max)
yvar <- log10(sdf[freq.indx.min:freq.indx.max])
if (is.element(method,"smoothed")){
# check n.block, the number of logarithmic paritions
# of the SDF. if null, use the floor(log2(n.sample))
# to maintain similarity to DWT decompositions
if (is.null(n.block))
n.block <- as.integer(floor(logb(length(x),base=2)))
checkScalarType(n.block,"integer")
# partition SDF into nonoverlapping dyadic length
# blocks and obtain block averages
lxvar <- length(xvar)
delx <- (xvar[lxvar] - xvar[1])/n.block
xbox <- c(xvar[1] + ((1:n.block) - 0.5) * delx)
xbdy <- as.numeric(c(xvar[1], (xbox + 0.5 * delx)))
ybox <- vector(mode="double", length=n.block)
for(k in 1:n.block){
ybox[k] <- mean(yvar[(xvar >= xbdy[k]) & (xvar <= xbdy[k + 1])])
}
xvar <- xbox
yvar <- ybox
}
w <- NULL # quell R CMD check scoping issue
logfit <- ifelse1(lmfit,
fit(y ~ 1 + x, data=data.frame(x=xvar, y=yvar, w=weight(yvar)), weights=w),
fit(y ~ 1 + x, data=data.frame(x=xvar, y=yvar)))
# perform trimmed least squares fit to get 1 - 2 H
exponent <- logfit$coefficients["x"]
H <- (1 - exponent) / 2
}
else if (is.element(method,"robinson")){
n.freq <- length(sdf)
# form cumulative sum of periodogram
#spec.low <- sdf[freq.indx.min:freq.indx.max]
spec.low <- sdf[1:freq.indx.max]
csumspec <- cumsum(spec.low)
# Set initial value of q, and a large enough delH to ensure
# at least one iteration
parmq <- 0.6
delH <- 1
# Evaluate Robinson estimate of H
while(abs(delH) > 0.001){
qm <- floor(parmq * freq.indx.max)
NewH1 <- 1 - log(csumspec[qm]/csumspec[freq.indx.max])/(2 * log(parmq))
if (NewH1 >= 0.75){
# Since H >= 0.75, use suboptimal expression;
# take avg of derivative over three nearest frequencies
NewH1 <- 1 - (((pi * qm)/n.freq) * mean(spec.low[(
qm - 2):qm]))/csumspec[qm]
if (NewH1 >= 0.75){
# terminate with current value;
# it is the best we can do for now
AvgH <- NewH1
# set flag
parmq <- 1
}
else{
# do another iteration using a smaller q
AvgH <- NewH1 - 1
parmq <- 0.9 * parmq
}
}
else{
# Since H < 0.75, we can iterate
# to optimal q(H). First calculate H(Old q)
NewH2 <- 0.75 - exp((0.01829 - parmq)/0.1292)
# average previous H and H(Old q)
AvgH <- mean(c(NewH1, NewH2))
# obtain new q=q(AvgH)
parmq <- 0.01829 - 0.1292 * log(0.75 - AvgH)
}
delH <- NewH1 - AvgH
xvar <- yvar <- logfit <- NULL
}
H <- NewH1
# estimate stdev of the inferred H
#StdofH <- sqrt(Hvar(H, parmq)/freq.indx.max)
# And drop H, optimum parameter q,
# and the estimated std of H
#drop(list(Hparm=NewH1, Qopt=parmq, StdH=StdofH))
} # robinson
# return Hurst coefficient estimate
names(H) <- "H"
# return the result as a fractal block exponent
fractalBlock(
domain = "frequency",
estimator = "Hurst coefficient via regression of nonparametric SDF estimate",
exponent = H,
exponent.name = "H",
scale = xvar,
stat = yvar,
stat.name = switch(method,standard="SDF",smoothed="mean of SDF blocks",robinson="Robinson Integration"),
detrend = NULL,
overlap = NULL,
data.name = data.name,
sum.order = 0,
series = asVector(x),
logfit = logfit,
sdf = sdf)
}
##
# RoverS
##
"RoverS" <- function(x, n.block.min=2, scale.ratio=2, scale.min=8)
{
checkScalarType(n.block.min,"integer")
checkScalarType(scale.ratio,"numeric")
checkScalarType(scale.min,"numeric")
# define local functions
"Rstat" <- function(cSumSeries, timeValue, shiftValue){
# Provide storage for the differenced cumulative sum:
difSum <- vector(mode="double", length=shiftValue + 1)
# Evaluate the difSum vector and the scalar difference dSumk:
difSum <- cSumSeries[timeValue:(timeValue + shiftValue)] - cSumSeries[timeValue]
dSumk <- cSumSeries[timeValue + shiftValue] - cSumSeries[timeValue]
# Evaluate max and min range:
RMax <- max(difSum - (c(timeValue:(timeValue + shiftValue)) * dSumk)/shiftValue)
Rmin <- min(difSum - (c(timeValue:(timeValue + shiftValue)) * dSumk)/shiftValue)
# return the R-statistic
return(RMax - Rmin)
}
# convert data, if necessary, to signalSeries class:
x <- wmtsa::create.signalSeries(x)
# find length of series:
n.sample <- length(x@data)
# and find number of data blocks that can be used:
nBlocksMax <- floor(n.sample/n.block.min)
# find max number of times number of blocks can be increased,
# and provide storage for result vectors:
kmax <- floor(1 + log(n.sample/(n.block.min * scale.min))/log(scale.ratio))
RoS <- vector(mode="double", length=kmax * nBlocksMax)
scale <- vector(mode="double", length=kmax * nBlocksMax)
# do cumulative sum of time series once and for all:
cSumser <- cumsum(x@data)
# set up initial values for while loop:
numBlocks <- n.block.min
k <- 1
indx <- 0
# and begin loop to generate R/S statistic vs k:
while(k <= kmax) {
nBlocks <- floor(n.sample/numBlocks)
nUpper <- nBlocks * numBlocks
imax <- floor(nBlocks - 1 - 1/numBlocks)
for (i in (1:imax)){
indx <- indx + 1
t <- numBlocks * i + 1
R <- Rstat(cSumser, t, numBlocks)
S <- variance(x[(t + 1):(t + numBlocks)], unbiased=FALSE)
S <- sqrt(S)
scale[indx] <- numBlocks
RoS[indx] <- R/S
}
# increase numBlocks and k; then (maybe) go around again:
numBlocks <- numBlocks * scale.ratio
k <- k + 1
}
# do linear regression over all points that have positive R/S:
xvar <- log10(scale[RoS > 0])
yvar <- log10(RoS[RoS > 0])
# remove NA values
xbad <- which(is.na(xvar))
ybad <- which(is.na(yvar))
bad <- union(xbad, ybad)
if (length(bad)){
xvar <- xvar[-bad]
yvar <- yvar[-bad]
}
fitcoef <- lsfit(xvar, yvar)$coef[1:2]
# and evaluate H:
RoverSH <- fitcoef[2]
RoverSH
}
###
## pgramReg
###
#
#"pgramReg" <- function(sdf, freq.indx.min=2, freq.indx.max=32, plot=FALSE)
#{
# # Convert data, if necessary, to signalSeries class:
# x <- create.signalSeries(sdf)
# spec <- x@data
#
# # set up variables for linear regression:
# xvar <- log10(freq.indx.min:freq.indx.max)
# yvar <- log10(spec[freq.indx.min:freq.indx.max])
#
# fitcoef <- lsfit(xvar, yvar)$coef[1:2]
#
# # return Hurst coefficient estimate
# (1 - fitcoef[2]) / 2
#}
#"pgramRegMod" <- function(sdf, freq.indx.min=2, freq.indx.max=32, Nbox=8, plot=FALSE)
#{
# # convert data, if necessary, to signalSeries class:
# x <- create.signalSeries(sdf)
# spec <- x@data
#
# # set up variables for linear regression:
# xvar <- log10(freq.indx.min:freq.indx.max)
# yvar <- log10(sdf[freq.indx.min:freq.indx.max])
# ybox <- vector(mode="double", length=Nbox)
#
# # assign to boxes and get box averages:
# lxvar <- length(xvar)
# delx <- (xvar[lxvar] - xvar[1])/Nbox
# xbox <- c(xvar[1] + ((1:Nbox) - 0.5) * delx)
# xbdy <- c(xvar[1], (xbox + 0.5 * delx))
# for(k in 1:Nbox){
# ybox[k] <- mean(yvar[(xvar >= as.numeric(xbdy[k])) & (xvar <=
# as.numeric(xbdy[k + 1]))])
# }
#
# # perform trimmed least squares fit to get 1 - 2 H; solve for H
# fitcoef <- ltsreg(xbox, ybox)$coef[1:2]
#
# if(plot == T) {
# plot(xbox, ybox, xlab="log10(Frequency Index)", ylab =
# "log10(SDF)")
# panel.abline(fitcoef)
# }
#
# # return Hurst coefficient estimate
# (1 - fitcoef[2]) / 2
#}
##
# RobInt
##
#"RobInt" <- function(sdf, freq.max=0.2)
#{
# #TODO: what is Hvar(q=1) supposed to be?
# # define local functions
# "Hvar" <- function(H, q)
# {
# # Evaluates variance estimate of distribution of parameter H
# # Ref: I. Lobato and P. M. Robinson,
# temp1 <- (1 + 1/q - 2 * (q^(1 - 2 * H)))/((log(q))^2)
# (temp1 * ((1 - H)^2))/(3 - 4 * H)
# }
#
# if (!is(sdf, "SDF"))
# stop("sdf must be an object of class \"SDF\"")
# if (!attr(sdf,"single.sided"))
# stop("sdf must be single-sided")
#
# 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)")
#
# # associate freq.max normalized frequency
# # with corresponding index into SDF
# # frequency vector
# df <- attr(sdf, "deltaf")
# freq.indx.max <- ceiling(freq.max/df)
# spec <- as.vector(sdf)
# n.freq <- length(spec)
#
# #TODO: does RobInt need to be a periodogram?
# # form cumulative sum of periodogram
# spec.low <- spec[seq(freq.indx.max)]
# csumspec <- cumsum(spec.low)
#
# # Set initial value of q, and a large enough delH to ensure
# # at least one iteration
# parmq <- 0.6
# delH <- 1
#
# # Evaluate Robinson estimate of H
# while(abs(delH) > 0.001){
#
# qm <- floor(parmq * freq.indx.max)
# NewH1 <- 1 - log(csumspec[qm]/csumspec[freq.indx.max])/(2 * log(parmq))
#
# if (NewH1 >= 0.75){
#
# # Since H >= 0.75, use suboptimal expression;
# # take avg of derivative over three nearest frequencies
# NewH1 <- 1 - (((pi * qm)/n.freq) * mean(spec.low[(
# qm - 2):qm]))/csumspec[qm]
#
# if (NewH1 >= 0.75){
# # terminate with current value;
# # it is the best we can do for now
# AvgH <- NewH1
#
# # set flag
# parmq <- 1
# }
# else{
#
# # do another iteration using a smaller q
# AvgH <- NewH1 - 1
# parmq <- 0.9 * parmq
# }
# }
# else{
#
# # Since H < 0.75, we can iterate
# # to optimal q(H). First calculate H(Old q)
# NewH2 <- 0.75 - exp((0.01829 - parmq)/0.1292)
#
# # average previous H and H(Old q)
# AvgH <- mean(c(NewH1, NewH2))
#
# # obtain new q=q(AvgH)
# parmq <- 0.01829 - 0.1292 * log(0.75 - AvgH)
# }
#
# delH <- NewH1 - AvgH
# }
#
#
# # estimate stdev of the inferred H
# StdofH <- sqrt(Hvar(NewH1, parmq)/freq.indx.max)
#
# # And drop H, optimum parameter q,
# # and the estimated std of H
# drop(list(Hparm=NewH1, Qopt=parmq, StdH=StdofH))
#}
###
## disWhittle
###
#
#"disWhittle" <- function(sdf, NFT, freq.indx.min=2, freq.indx.max=floor((NFT + 1)/2))
#{
# # convert datato signalSeries class
# x <- create.signalSeries(sdf)
# Sx <- x@data
#
# # define function to evaluate the weighted sum of the spectrum
# Isumd <- function(delta, Sx, n.freq, Mindx, Maxdx)
# {
# xatt <- attributes(Sx)
# f <- xatt$frequency
# N <- xatt$n.sample
# Np <- (N - 1) %/% 2
#
# # f <- c(1:(Maxdx - Mindx + 1)) / n.freq
# # Nstar <- length(f)
#
# Isum1 <- mean(as.numeric(Sx) * (4 * (sin(pi * f))^2)^delta)
# log(Isum1) + 2 * delta * (((N/(2 * Np)) - 1) *
# log(2) - log(2 * N)/(2 * Np))
# }
#
# # call optimize() function to find the best value of delta
# dfinal <- optimize(Isumd, lower=0, upper=0.5, tol=0.0001,
# maximum=F, Sx=Sx[freq.indx.min:freq.indx.max], n.freq=NFT,
# Mindx=freq.indx.min, Maxdx=freq.indx.max)
#
# # return Hurst coefficient estimate
# dfinal$minimum + 0.5
#}
#
#######################################################################################
#
#"Whittle" <- function(sdf, NFT, freq.indx.min=2, freq.indx.max=floor((NFT + 1)/2))
#{
# # Convert data, if necessary, to signalSeries class:
# x <- create.signalSeries(sdf)
# spec <- x@data
#
# # Here is the function to evaluate the weighted sum of the spectrum:
# Isumd <- function(dparm, spec, n.freq, Mindx, Maxdx)
# {
# f <- c(1:(Maxdx - Mindx + 1))/n.freq
# Isumd <- sum(spec * (4 * (sin(pi * f))^2)^dparm)/n.freq
# Isumd
# }
#
# # Just call the S-Plus "optimize" function to find the best
# # value of delta:
# dfinal <- optimize(Isumd, lower=0, upper=0.5, tol=0.001,
# maximum=F, spec=spec[freq.indx.min:freq.indx.max], n.freq=NFT,
# Mindx=freq.indx.min, Maxdx=freq.indx.max)
#
# # return Hurst coefficient estimate
# dfinal$minimum + 0.5
#}
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.