Nothing
#' Wavelet-based Maximum Likelihood Estimation for Seasonal Persistent
#' Processes
#'
#' Parameter estimation for a seasonal persistent (seasonal long-memory)
#' process is performed via maximum likelihood on the wavelet coefficients.
#'
#' The variance-covariance matrix of the original time series is approximated
#' by its wavelet-based equivalent. A Whittle-type likelihood is then
#' constructed where the sums of squared wavelet coefficients are compared to
#' bandpass filtered version of the true spectral density function.
#' Minimization occurs for the fractional difference parameter \eqn{d} and the
#' Gegenbauer frequency \eqn{f_G}, while the innovations variance is
#' subsequently estimated.
#'
#' @usage spp.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, frac = 1)
#' @usage spp2.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, dyadic = TRUE, frac = 1)
#' @aliases spp.mle spp2.mle
#' @param y Not necessarily dyadic length time series.
#' @param wf Name of the wavelet filter to use in the decomposition. See
#' \code{\link{wave.filter}} for those wavelet filters available.
#' @param J Depth of the discrete wavelet packet transform.
#' @param p Level of significance for the white noise testing procedure.
#' @param dyadic Logical parameter indicating whether or not the original time
#' series is dyadic in length.
#' @param frac Fraction of the time series that should be used in constructing
#' the likelihood function.
#' @return List containing the maximum likelihood estimates (MLEs) of
#' \eqn{\delta}, \eqn{f_G} and \eqn{\sigma^2}, along with the value of the
#' likelihood for those estimates.
#' @author B. Whitcher
#' @seealso \code{\link{fdp.mle}}
#' @references Whitcher, B. (2004) Wavelet-based estimation for seasonal
#' long-memory processes, \emph{Technometrics}, \bold{46}, No. 2, 225-238.
#' @keywords ts
#' @export spp.mle
spp.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01, frac=1)
{
sppLL <- function(x, y) {
delta <- x[1]
fG <- x[2]
## cat("Parameters are: d =", delta, ", and f =", fG, fill=TRUE)
y.dwpt <- y[[1]]
y.basis <- y[[2]]
n <- y[[3]]
J <- y[[4]]
## Establish the limits of integration for the band-pass variances
a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
## Define some useful parameters for the wavelet packet tree
# n <- length(y)
length.jn <- n / rep(2^(1:J), 2^(1:J))
scale.jn <- rep(2^(1:J+1), 2^(1:J))
## Initialize various parameters for the reduced LL
Basis <- (1:length(y.basis))[y.basis]
bp.var <- numeric(length(Basis))
delta.n <- 100
## Compute the band-pass variances according to \delta and f_G
omega.diag <- NULL
for(i in 1:sum(y.basis)) {
jn <- Basis[i]
bp.var[i] <- bandpass.spp(a[jn], b[jn], delta, fG)
omega.diag <- c(omega.diag,
scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
}
## Compute reduced log-likelihood
rLL <- n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var))
rLL
}
n <- length(y)
x0 <- numeric(2)
## Perform discrete wavelet packet transform (DWPT) on Y
y.dwpt <- dwpt(y, wf, n.levels=J)
n <- length(y)
if(frac < 1) {
for(i in 1:length(y.dwpt)) {
vec <- y.dwpt[[i]]
ni <- length(vec)
j <- rep(1:J, 2^(1:J))[i]
vec[trunc(frac * n/2^j):ni] <- NA
y.dwpt[[i]] <- vec
}
}
y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p)))
y.dwpt <- as.matrix(unlist(y.dwpt[y.basis]))
## Compute initial estimate of the Gegenbauer frequency
y.per <- per(y - mean(y))
x0[2] <- (0:(n/2)/n)[max(y.per) == y.per]
## Compute initial estimate of the fractional difference parameter
muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) +
unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J)))) / 4
y.modwpt <- modwpt(y, wf=wf, n.levels=J)
y.varJ <- rep(2^(1:J), 2^(1:J)) *
unlist(lapply(y.modwpt,
FUN=function(x)sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
x0[1] <- min(-0.5 * lsfit(log(abs(muJ[y.basis] - x0[2])),
log(y.varJ[y.basis]))$coef[2], 0.49)
cat(paste("Initial parameters are: delta =", round(x0[1],4),
"freqG =", round(x0[2],4), "\n"))
result <- optim(par=x0, fn=sppLL, method="L-BFGS-B",
lower=c(0.001,0.001), upper=c(0.499,0.499),
control=list(trace=0, fnscale=2),
y=list(y.dwpt, y.basis, n, J))
return(result)
}
spp2.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01,
dyadic=TRUE, frac=1)
{
spp2LL <- function(x, y) {
d1 <- x[1]
f1 <- x[2]
d2 <- x[3]
f2 <- x[4]
## cat("Parameters are: d1 =", round(d1,6), ", and f1 =", round(f1,6),
## ", d2 =", round(d2,6), ", and f2 =", round(f2,6), fill=TRUE)
y.dwpt <- y[[1]]
y.basis <- y[[2]]
n <- y[[3]]
J <- y[[4]]
## Establish the limits of integration for the band-pass variances
a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
## Define some useful parameters for the wavelet packet tree
length.jn <- n / rep(2^(1:J), 2^(1:J))
scale.jn <- rep(2^(1:J+1), 2^(1:J))
## Initialize various parameters for the reduced LL
Basis <- (1:length(y.basis))[y.basis]
bp.var <- numeric(length(Basis))
delta.n <- 100
## Compute the band-pass variances according to \delta and f_G
omega.diag <- NULL
for(i in 1:sum(y.basis)) {
jn <- Basis[i]
bp.var[i] <- bandpass.spp2(a[jn], b[jn], d1, f1, d2, f2)
omega.diag <- c(omega.diag,
scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
}
## Compute reduced log-likelihood
n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var), na.rm=TRUE)
}
n <- length(y)
x0 <- numeric(4)
## Perform discrete wavelet packet transform (DWPT) on Y
y.dwpt <- dwpt(y, wf, n.levels=J)
if(!dyadic) {
for(i in 1:length(y.dwpt)) {
vec <- y.dwpt[[i]]
ni <- length(vec)
j <- rep(1:J, 2^(1:J))[i]
vec[trunc(frac * n/2^j):ni] <- NA
y.dwpt[[i]] <- vec
}
}
y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p, type="other")))
y.dwpt <- as.vector(unlist(y.dwpt[y.basis]))
## Compute initial estimate of the Gegenbauer frequencies
if(dyadic)
y.per <- per(y - mean(y))
else
y.per <- per(y[1:(frac*n)] - mean(y[1:(frac*n)]))
freq.y <- (0:(frac*n %/% 2))/(frac*n)
x0[2] <- freq.y[max(y.per) == y.per]
x0[4] <- freq.y[max(y.per[freq.y > x0[2] + freq.y[10] |
freq.y < x0[2] - freq.y[10]]) == y.per]
if(x0[2] > x0[4]) {
xx <- x0[2]
x0[2] <- x0[4]
x0[4] <- xx
rm(xx)
}
## Compute initial estimate of the fractional difference parameters
muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) +
unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J)))) / 4
y.modwpt <- modwpt(y, wf=wf, n.levels=J)
y.varJ <- rep(2^(1:J), 2^(1:J)) *
unlist(lapply(y.modwpt,
FUN = function(x) sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
x0.mid <- (x0[2] + x0[4]) / 2
muJ <- muJ[y.basis]
y.varJ <- y.varJ[y.basis]
x0[1] <- min(-0.5 * lsfit(log(abs(muJ[muJ < x0.mid] - x0[2])),
log(y.varJ[muJ < x0.mid]))$coef[2], 0.49)
x0[3] <- min(-0.5 * lsfit(log(abs(muJ[muJ > x0.mid] - x0[4])),
log(y.varJ[muJ > x0.mid]))$coef[2], 0.49)
cat(paste("Initial parameters: d1 = ", round(x0[1],4),
", f1 = ", round(x0[2],4), ", d2 = ", round(x0[3],4),
", f2 = ", round(x0[4],4), sep=""), fill=TRUE)
result <- optim(par=x0, fn=spp2LL, method="L-BFGS-B",
lower=rep(0.001,4), upper=rep(0.499,4),
control=list(trace=1, fnscale=2),
y=list(y.dwpt, y.basis, n, J))
return(result)
}
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.