#############################################################
## WMTSA fractionally differenced process functionality
##
## Functions:
##
## wavFDPBand
## wavFDPBlock
## wavFDPSDF
## wavFDPTime
##
## Constructor Functions and methods:
##
## wavFDP
##
## eda.plot.wavFDP
## plot.wavFDP
## print.wavFDP
## print.summary.wavFDP
## summary.wavFDP
##
#############################################################
###
# wavFDPBand
##
"wavFDPBand" <- function(delta=1/4, method="bandpass", scaling=TRUE,
levels=1:5, n.sample=n.sample <- 2^(max(levels)+1))
{
# force the levels to be monotonic and uniform if scaling is FALSE
# i.e. ensure that j=1, ..., J and that no levels are skipped
if (scaling && length(levels) != max(levels))
stop("If scaling is TRUE, the levels vector must contain ",
"indices j=1, 2, ..., max(levels)")
method <- match.arg(method, c("bandpass", "integration"))
if (method == "integration"){
# develop band limits of integration
lower.limit <- 2^(-(levels + 1))
upper.limit <- 2^(-levels)
mult.factor <- 2^(levels + 1)
# set the process variance to unity for now until we establish
# an MLE for delta
variance <- 1
# calculate the average SDF value at the center frequency for
# each band, i.e. C'_j for j=1:J
Cprime <- unlist(lapply(levels,
function(j, delta, lim.low, lim.high, mult.factor){
y <- integrate(f=wavFDPSDF, lower=lim.low[j], upper=lim.high[j],
delta=delta, variance=1, response=NULL)
mult.factor[j] * y$value
},
delta = delta,
lim.low = lower.limit,
lim.high = upper.limit,
mult.factor = mult.factor)
)
"fdpsdf.bandpass.scaling" <- function(n.sample, delta, variance, Cprime)
{
J <- length(Cprime)
levels <- seq_len(J)
n.sample * ( variance * gamma(1 - 2 * delta) / gamma(1 - delta)^2 - sum(Cprime / 2^levels) )
}
# now calculate the average SDF value
# for the scaling coefficient(s) band
if (scaling){
CprimeJP1 <- fdpsdf.bandpass.scaling(n.sample, delta, variance, Cprime)
Cprime <- c(Cprime, CprimeJP1 )
}
}
else if (method == "bandpass"){
if (!scaling)
n.sample <- -1
Cprime <- as.vector(itCall("RS_wavelets_fdp_bandpass_variance",
as.integer(levels),
as.numeric(delta),
as.integer(n.sample))
#,
#COPY=rep(FALSE,3),
#CLASSES=c("integer","numeric","integer"),
#PACKAGE="ifultools")
)
}
else{
stop("Unsupported method specified for mid-octave SDF value approximation")
}
# return the vector of average SDF bands
crystals <- paste("d", levels, sep="")
if (scaling)
crystals <- c(crystals, paste("s",max(levels),sep=""))
names(Cprime) <- crystals
Cprime
}
###
# wavFDPBlock
##
"wavFDPBlock" <- function(x, wavelet="s8", levels=NULL, sdf=NULL,
boundary=NULL, edof.mode=1,
estimator="wlse", delta.range=c(-10.0,10.0),
position=list(from=1,by=1,units=character()), units=character(),
title.data=character(), documentation=character(), keep.series=FALSE )
{
# get the series name from the parent
# (or from this function if this is called as a top level function)
series.title <- deparse(substitute(x))
# map estimator
estimator <- match.arg(estimator, c("wlse", "mle"))
# map filter number and obtain length
filter <- mutilsFilterType(wavelet=wavelet)
# obtain filters (needed for storage in dictionary)
filters <- wavDaubechies(wavelet=wavelet, normalized=TRUE)
# store series
series <- ifelse1(keep.series,
create.signalSeries(x, position=position, units=units,
title.data=title.data, documentation=documentation),
create.signalSeries(NULL, position=position, units=units,
title.data=series.title, documentation=documentation))
# map estimator
estimator <- lowerCase(estimator)
full.decomposition <- FALSE
# set default for boundary variable and make
# sure that it coordinates with the estimator
if (is.null(boundary))
boundary <- ifelse1(estimator == "mle", "nonstationary", "unbiased")
if (estimator == "wlse"){
estim <- 1
imatch <- charmatch(boundary, c("biased", "unbiased"))
if (is.na(imatch))
stop("\nThe weighted least squares estimator requires that\n",
"the boundary mode be either \"biased\" or \"unbiased\"")
if (imatch == 1)
bound <- list(mode=TRUE, description="biased")
else if (imatch == 2)
bound <- list(mode=FALSE, description="unbiased")
else
stop("Boundary method is unsupported")
}
else if (estimator == "mle"){
estim <- 0
imatch <- charmatch(boundary, c("stationary", "nonstationary"))
if (is.na(imatch))
stop("\nThe maximum likelihood estimator requires that\n",
"the boundary mode be either \"stationary\" or \"nonstationary\"")
if (imatch == 1){
bound <- list(mode=TRUE, description="stationary FD process model")
full.decomposition <- TRUE
}
else if (imatch == 2)
bound <- list(mode=FALSE,
description="stationary-nonstationary FD process model")
else
stop("Boundary method is unsupported")
}
else
stop("Unsupported estimator method")
# if not supplied, define the LEVELS vector
max.level <- wavMaxLevel(n.taps=filter$length, n.sample=length(x),
xform=switch(estimator, mle="modwt", wlse="dwt"))
if (is.null(levels))
levels <- ifelse1(length(x) > filter$length,
ifelse1(max.level > 2, seq(from=2, length=max.level - 1), seq(1, max.level)), 1)
if (full.decomposition)
levels <- seq(1, max.level)
if (max(levels) > max.level)
stop(paste("\nNumber of levels exceeds maximum:", max.level,
"\nEither reduce the filter length or maximum decomposition level."))
names(levels) <- paste("d",levels,sep="")
# calculate the block-averaged FD model
# parameter estimates
if (is.null(sdf) || !length(sdf)){
sdf <- (1:2)
sdf.empty <- TRUE
}
else
sdf.empty <- FALSE
data <- itCall("RS_wavelets_fdp_estimator_block",
as.double(x),
as.integer(levels),
as.integer(filter$type),
as.integer(filter$length),
as.integer(estim),
as.logical(bound$mode),
as.integer(edof.mode),
as.double(sdf),
as.double(delta.range))
#,
#COPY=rep(FALSE,9),
#CLASSES=c("numeric",rep("integer",4),"logical","integer","numeric","numeric"),
#PACKAGE="ifultools")
# remove values outside of delta range. these values
# are assigned by the C code as a flag for the case
# where there are not enough interior coefficients
# in which to form a maximum likelihood estimation
na.indices <- which(data[[1]] < min(delta.range) | data[[1]] > max(delta.range))
data <- lapply(data, function(x, na.indices){
x[na.indices] <- NA
x}, na.indices=na.indices)
# create a dictionary
dictionary <- wavDictionary(wavelet=wavelet,
dual=FALSE, decimate=FALSE, n.sample=length(x),
attr.x=NULL, n.levels=length(levels),
boundary="periodic", conv=TRUE, filters=filters,
fast=TRUE, is.complex=FALSE)
# pack results into a list and return
wavFDP(estimator=estimator,
delta=as.vector(data[[1]]),
variance.delta=ifelse1(estimator == "mle", NA, as.vector(data[[2]])),
innovations.variance=ifelse1(estimator == "wlse", NA, as.vector(data[[3]])),
delta.range=delta.range,
dictionary=dictionary,
levels=as.integer(levels),
edof.mode=as.integer(edof.mode),
boundary=bound,
series=series,
sdf.method="Integration lookup table",
type="block")
}
###
# wavFDPSDF
##
"wavFDPSDF" <- function(f, delta=0.45, variance=1, response=NULL)
{
mult.fact <- ifelse1(!is.null(response),
approx(response$frequency, response$sqrgain, f, rule=2)$y, 1)
mult.fact * variance / (2 * sin(pi * f))^(2 * delta)
}
###
# wavFDPTime
##
"wavFDPTime" <- function(x, wavelet="s8", levels=NULL,
biased=FALSE, estimator="mle",
dof.order=0, delta.range=c(-10.0,10.0),
position=list(from=1,by=1,units=character()), units=character(),
title.data=character(), documentation=character(), keep.series=FALSE)
{
# get the series name from the parent
# (or from this function if this is called as a top level function)
series.title <- deparse(substitute(x))
# map filter number and obtain length
filter <- mutilsFilterType(wavelet=wavelet)
# obtain filters (needed for storage in dictionary)
filters <- wavDaubechies(wavelet=wavelet, normalized=TRUE)
# set the default number of decomposition levels
# to be the largest scale at which there exists
# at least one interior (non-boundary) wavelet
# coefficient
if (is.null(levels)){
L <- filter$length
N <- length(x)
if (!biased){
if (N > L)
n.levels <- ilogb(((N - 1) / (L - 1) + 1), base=2)
else if (N > 0)
n.levels <- 1
else stop("length of time series must be positive")
}
else
n.levels <- 2
levels <- seq(n.levels)
}
names(levels) <- paste("d",levels,sep="")
# store series
series <- ifelse1(keep.series,
create.signalSeries(x, position=position, units=units,
title.data=title.data, documentation=documentation),
create.signalSeries(NULL, position=position, units=units,
title.data=series.title, documentation=documentation))
# test estimator
estimator <- match.arg(lowerCase(estimator), c("mle","lse"))
method <- as.integer(estimator == "lse")
# obtain MLEs of FDP paramter delta
# and corresponding innovation variance
data <- itCall("RS_wavelets_fdp_estimator_instantaneous",
as.numeric(x),
as.integer(levels),
as.integer(filter$type),
as.integer(filter$length),
as.integer(method),
as.logical(biased),
as.integer(dof.order),
as.numeric(delta.range))
#,
#COPY=rep(FALSE,8),
#CLASSES=c("numeric",rep("integer",4),"logical","integer","numeric"),
#PACKAGE="ifultools")
# remove values outside of delta range. these values
# are assigned by the C code as a flag for the case
# where there are not enough interior coefficients
# in which to form a maximum likelihood estimation
na.indices <- which(data[[1]] < min(delta.range) | data[[1]] > max(delta.range))
data <- lapply(data, function(x, na.indices){
x[na.indices] <- NA
x}, na.indices=na.indices)
# create a dictionary
dictionary <- wavDictionary(wavelet=wavelet,
dual=FALSE, decimate=FALSE, n.sample=length(x),
attr.x=NULL, n.levels=length(levels),
boundary="periodic", conv=TRUE, filters=filters,
fast=TRUE, is.complex=FALSE)
# pack results into a list and return
boundary <- ifelse1(biased, list(mode=TRUE, description="biased"),
list(mode=FALSE, description="unbiased"))
# pack results into a list and return
wavFDP(estimator=estimator,
delta=as.vector(data[[1]]),
variance.delta=as.vector(data[[2]]),
innovations.variance=as.vector(data[[3]]),
delta.range=delta.range,
dictionary=dictionary,
levels=as.integer(levels),
edof.mode=as.integer(dof.order),
boundary=boundary,
series=series,
sdf.method="Integration lookup table",
type="instantaneous")
}
################################################
##
## Constructor Functions and methods:
##
## wavFDP
##
## eda.plot.wavFDP
## plot.wavFDP
## print.wavFDP
## print.summary.wavFDP
## summary.wavFDP
##
################################################
###
# wavFDP (constructor)
##
"wavFDP" <- function(estimator, delta, variance.delta,
innovations.variance, delta.range, dictionary, levels,
edof.mode, boundary, series, sdf.method, type)
{
# check input arguments
checkScalarType(estimator, "character")
checkVectorType(delta)
if (all(!is.na(variance.delta)))
checkVectorType(variance.delta)
if (all(!is.na(innovations.variance)))
checkVectorType(innovations.variance)
checkVectorType(delta.range)
checkVectorType(delta)
if (!is(dictionary, "wavDictionary"))
stop("dictionary must be an object of class \"wavDictionary\"")
checkVectorType(levels, "integer")
checkScalarType(edof.mode, "integer")
if (!is.list(boundary) || !all(is.element(names(boundary), c("mode","description"))))
stop("boundary must be a list with named objects ",
"\"mode\" and \"description\"")
checkScalarType(boundary$mode, "logical")
checkScalarType(boundary$description, "character")
if (!is(series, "signalSeries"))
stop("series must be an object of class \"signalSeries\"")
checkScalarType(sdf.method, "character")
checkScalarType(type, "character")
type <- match.arg(type, c("block","instantaneous"))
z <- list(estimator=estimator,
delta=delta,
variance.delta=variance.delta,
innovations.variance=innovations.variance,
delta.range=delta.range,
dictionary=dictionary,
levels=levels,
edof.mode=edof.mode,
boundary=boundary,
series=series,
sdf.method=sdf.method,
type=type)
oldClass(z) <- "wavFDP"
return(z)
}
###
# print.wavFDP
##
"print.wavFDP" <- function(x, digits=5, justify="left", sep=":", ...)
{
# define series string
name <- ifelse1(!length(x$series@title), "", paste("of", x$series@title))
type <- lowerCase(x$type)
est <- lowerCase(x$estimator)
if (type == "block"){
main <- paste("Block-dependent FD parameter estimation", name)
# display block dependent results
delta <- x$delta
var.delta <- x$variance.delta
innov <- x$innovations.variance
# remove certain figures until resolved
if (est == "mle")
var.delta <- NA
else if (est == "wlse")
innov <- NA
z <- list(
"FD parameter estimate (delta)"=round(delta, digits),
"var{delta} estimate"=round(var.delta, digits),
"Innovations variance estimate"=round(innov, digits))
}
else if (type == "instantaneous"){
#cat("\nInstantaneous FD parameter estimation", name, "\n\n")
main <- paste("Instantaneous FD parameter estimation", name)
# display block dependent results
mean.delta <- mean(x$delta, na.rm=TRUE)
mean.var.delta <- mean(x$variance.delta, na.rm=TRUE)
mean.innovations <- mean(x$innovations.variance, na.rm=TRUE)
# remove certain figures until resolved
if (est == "mle")
mean.var.delta <- NA
else if (est == "lse")
mean.innovations <- NA
z <- list(
"Mean of FD parameter estimates (deltas)"=round(mean.delta, digits),
"Mean of var{delta} estimates"=round(mean.var.delta, digits),
"Mean of innovations variance estimates"=round(mean.innovations, digits))
}
else{
main <- "Unknown method."
}
dict <- x$dictionary
z <- c(z, list(
"Estimator"=upperCase(x$estimator),
"Levels"=x$levels,
"Boundary mode"=x$boundary$description,
"EDOF mode"= ifelse1(est == "wlse" && type == "block", x$edof, NULL),
"Chi-squared DOF"=ifelse1(type == "instantaneous", x$edof * 2 + 1, NULL),
"Delta range"=x$delta.range,
"Wavelet"=dict$wavelet,
"Length of series"=dict$n.sample,
"Number of levels"=dict$n.levels,
"Boundary correction rule"=dict$boundary,
"Filtering technique"=ifelse1(dict$conv,"Convolution","Correlation")))
prettyPrintList(z, header=main, sep=sep, justify=justify, ...)
}
###
# summary.wavFDP
##
"summary.wavFDP" <- function(object, ...)
{
obj <- object
oldClass(obj) <- "summary.wavFDP"
return(obj)
}
###
# print.summary.wavFDP
##
"print.summary.wavFDP" <- function(x, ...)
{
print.wavFDP(x, ...)
}
###
# plot.wavFDP
##
"plot.wavFDP" <- function(x, mean.delta=NULL, xlab="Time",
ylab=paste(ifelse1(x$boundary$mode, "Biased ", "Unbiased "), upperCase(x$estimator),"(delta)",sep=""),
title.str=NULL, type="l", show.key=TRUE, conf.color="lightblue", ...)
{
# save plot parameters and restore upon exit
old.plt <- splitplot(1,1,1)
on.exit(par(old.plt))
# local functions
# shade the area between two discrete curves
"patchLimits" <- function(x, y, t=seq(x), color=13, border=FALSE, ...)
{
# perform argument checking
if (length(x) != length(y) | length(x) != length(t))
stop("All input vectors must be the same length")
if (any(is.na(c(x,y,t))))
stop("NA value not allowed")
# concatenate the vertices in the proper order
xdata <- c(t, rev(t))
ydata <- c(x, rev(y))
# plot the points without actually plotting the points
plot(xdata, ydata,type="n", ...)
# apply polygon shaded area
polygon(xdata, ydata, col=color, border=border)
}
# check for instantaneous
if (x$type != "instantaneous")
stop("The plot method for FDP class objects is only supported for instantaneous estimators")
# obtain the FDP parameters
delta <- x$delta
stdev <- sqrt(x$variance.delta);
bad <- which (is.na(delta));
good <- which (!is.na(delta));
# define local plot variables
pch.bad <- 5
refstr <- ifelse1(is.numeric(mean.delta), paste(", ref=", mean.delta, sep=""), "")
# for now this routine supports LSE biased/unbiased and MLE unbiased
use.confidence <- !x$boundary$mode & (lowerCase(x$estimator) == "lse")
# plot the error bars
if (use.confidence){
if (is.numeric(mean.delta)){
high <- mean.delta + stdev[good]
low <- mean.delta - stdev[good]
}
else{
low <- delta[good] - stdev[good]
high <- delta[good] + stdev[good]
}
ylim.low <- min(low, delta[good])
ylim.high <- max(high, delta[good])
ylim <- c(ylim.low, ylim.high)
patchLimits(low, high, t=good, xlab=xlab, ylab=ylab, ylim=ylim, color=conf.color)
# plot delta
lines(good, delta[good], type=type, col=1)
# plot boundary coefficient points
points(bad, rep(0, length(bad)), pch=pch.bad, col="red")
# add a legend
if (show.key){
legend("bottom",
c(paste("Delta(t)", refstr, sep=""),"NA Values","95% Confidence "),
lty=c(1,-1,-1),
pch=c(-1,pch.bad,15),
col=c(1,"red",conf.color),
cex=1.5)
}
}
else{
# plot delta
plot(good, delta[good], type=type, col=1, xlab=xlab, ylab=ylab)
# plot boundary coefficient points
points(bad, rep(0, length(bad)), pch=pch.bad, col="red")
# add a mean reference line
if (is.numeric(mean.delta)){
# plot mean reference line
lines(range(good), rep(mean.delta, 2), lty= 2, col=3)
# add a legend
if (show.key){
legend("bottom",
c("Delta(t)","NA Values", paste("Reference delta:", mean.delta," ")),
lty=c(1,1,2),
pch=c(-1,pch.bad,-1),
col=c(1,"red",conf.color),
cex=1.5)
}
}
else{
if (show.key){
legend("bottom",
c("Delta(t)","NA Values"),
lty=c(1,-1),
pch=c(-1,pch.bad),
col=c(1,"red"),
cex=1.5)
}
}
}
# add a title
if (is.null(title.str) && (x$type == "instantaneous"))
title.str <- paste("Instantaneous FD parameter estimation of", wavTitle(x))
if (!is.null(title.str))
title(title.str, cex=0.7)
invisible(NULL)
}
###
# eda.plot.wavFDP
##
"eda.plot.wavFDP" <- function(x, mean.delta=NULL, xlab="Time",
ylab=paste(ifelse1(x$boundary$mode, "Biased ", "Unbiased "), x$estimator,"(delta)",sep=""),
title.str=NULL, type="l", ...)
{
# local functions
"edatitle" <- function(x) title(paste("\n",x,sep=""), adj=1, cex=0.6)
# obtain the FDP parameters
delta <- x$delta
stdev <- sqrt(x$variance.delta);
bad <- which (is.na(delta));
good <- which (!is.na(delta));
if (is.null(mean.delta))
mean.delta <- mean(delta[good])
# set up a 2x2 grid of plots
old.plt <- splitplot(2,2,1)
on.exit(par(old.plt))
# create lines for the histogram ...
delta.std <- unique(stdev)
yloc <- c(mean.delta - delta.std, mean.delta + delta.std)
# plot the delta results with confidence intervals
plot(x, show.key=FALSE, mean.delta=mean.delta)
edatitle(paste("Delta(t) with 95% Confidence, ref=", round(mean.delta, 5), sep=""))
ylim <- par("yaxp")
# plot the probability density ...
splitplot(2,2,2)
hist(delta, prob=TRUE, xlab=ylab, horiz=TRUE, xpd=TRUE, ylim=ylim, col=2)
edatitle("Probability Density")
box(1)
xrange <- par("usr")[1:2]
# ... draw mean line
segments(xrange[1], mean.delta, xrange[2], mean.delta, lwd=3, col=8)
# ... draw confidence interval limits
xlocmin <- rep(xrange[1], length(yloc))
xlocmax <- rep(xrange[2], length(yloc))
segments(xlocmin, yloc, xlocmax, yloc, lty=2)
# plot the delta results with confidence intervals
splitplot(2,2,3)
plot(x, show.key=FALSE, mean.delta=NULL)
edatitle("Delta(t) with 95% Confidence")
ylim <- par("yaxp")
# plot qqnorm and qqline
splitplot(2,2,4)
ylab <- paste(ifelse1(x$boundary$mode, "Biased ", "Unbiased "), upperCase(x$estimator),"(delta)",sep="")
qqnorm(delta[good], ylab=ylab, mkh=0.02, ylim=ylim)
qqline(delta[good], col=8)
edatitle("QQ Normal")
invisible(NULL)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.