mix <- function(x, y) {
n <- 2 * length(x) - 2
rev(Re(fft(fft(symtacvf(x)) * fft(symtacvf(y)), inverse = TRUE)/n)[(n/2 - 1):(n - 1)])
}
#' Extracts the tacvfs of a fitted object
#'
#' Extracts the theoretical autocovariance functions (tacvfs) from a fitted
#' \code{arfima} or one of its modes (an \code{ARFIMA}) object.
#'
#'
#' @param obj An object of class "arfima" or "ARFIMA". The latter class is a
#' mode of the former.
#' @param xmaxlag The number of extra points to be added on to the end. That
#' is, if the original series has length 300, and xmaxlag = 5, the tacvfs will
#' go from lag 0 to lag 304.
#' @param forPred Should only be \code{TRUE} from a call to
#' \code{predict.arfima}.
#' @param n.ahead Only used internally.
#' @param nuse Only used internally.
#' @param \dots Optional arguments, currently not used.
#' @return A list of tacvfs, one for each mode, the length of the time series.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{plot.tacvf}}, \code{\link{print.tacvf}},
#' \code{\link{tacfplot}}, \code{\link{arfima}}
#' @references Veenstra, J.Q. Persistence and Antipersistence: Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @export tacvf
tacvf <- function(obj, xmaxlag = 0, forPred = FALSE, n.ahead = 0, nuse = -1,...) {
if (xmaxlag < 0)
stop("xmaxlag must be >= 0")
if (inherits(obj, "arfima")) {
if (!obj$weeded) {
warning("The object has not been weeded: there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
}
m <- length(obj$modes)
res <- vector("list", m + 1)
name <- deparse(substitute(obj))
res[[1]] <- name
period <- obj$period
dint <- obj$dint
dseas <- obj$dseas
if(nuse<0)
maxlag <- obj$n - dint - dseas * period + xmaxlag - 1
else
maxlag <- nuse - dint - dseas * period + xmaxlag - 1
if(maxlag < 1) stop('proportion of series used is not sufficient considering differencing!')
for (i in 1:m) {
phi <- obj$modes[[i]]$phi
theta <- obj$modes[[i]]$theta
dfrac <- obj$modes[[i]]$dfrac
phiseas <- obj$modes[[i]]$phiseas
thetaseas <- obj$modes[[i]]$thetaseas
dfs <- obj$modes[[i]]$dfs
H <- obj$modes[[i]]$H
Hs <- obj$modes[[i]]$Hs
alpha <- obj$modes[[i]]$alpha
alphas <- obj$modes[[i]]$alphas
sigma2 <- obj$modes[[i]]$sigma2
rr <- tacvfARFIMA(phi = phi, theta = theta, dfrac = dfrac, phiseas = phiseas,
thetaseas = thetaseas, dfs = dfs, H = H, Hs = Hs, alpha = alpha, alphas = alphas,
period = period, maxlag = maxlag, useCt = TRUE, sigma2 = sigma2)
if (!forPred)
res[[i + 1]] <- rr else {
res[[i + 1]]$muHat <- obj$modes[[i]]$muHat
res[[i + 1]]$tacvf <- rr
res[[i + 1]]$sigma2 <- obj$modes[[i]]$sigma2
if (length(H) == 0 && length(Hs) == 0) {
## need to do alpha, too.
if (n.ahead == 0)
stop("invalid n.ahead")
res[[i + 1]]$psis <- psiwts(phi = phi, theta = theta, phiseas = phiseas,
thetaseas = thetaseas, dfrac = dfrac, dfs = dfs, dint = dint, dseas = dseas,
period = period, len = n.ahead)
}
}
}
} else if (inherits(obj, "ARFIMA")) {
name <- deparse(substitute(obj))
phi <- obj$phi
theta <- obj$theta
dfrac <- obj$dfrac
phiseas <- obj$phiseas
thetaseas <- obj$thetaseas
dfs <- obj$dfs
H <- obj$H
Hs <- obj$Hs
alpha <- obj$alpha
alphas <- obj$alphas
period <- obj$period
maxlag <- length(obj$res) + xmaxlag - 1
sigma2 <- obj$sigma2
res <- vector("list", 2)
res[[1]] <- name
res[[2]] <- tacvfARFIMA(phi = phi, theta = theta, dfrac = dfrac, phiseas = phiseas,
thetaseas = thetaseas, H = H, Hs = Hs, alpha = alpha, alphas = alphas, period = period,
maxlag = maxlag, useCt = TRUE, sigma2 = sigma2)
} else stop("The tacvf function is only defined for arfima or ARFIMA objects.")
class(res) <- "tacvf"
res
}
#' Prints a tacvf object.
#'
#' Prints the output of a call to \code{tacvf} on an \code{arfima} object
#'
#'
#' @param x The \code{tacvf} object.
#' @param \dots Optional arguments. See \code{\link{print}}.
#' @return The object is returned invisibly
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{tacvf}}, \code{\link{plot.tacvf}}
#' @keywords ts
print.tacvf <- function(x, ...) {
m <- length(x) - 1
if (m > 1) {
cat("Printing the tacvfs of", x[[1]], "\n")
for (i in 2:(m + 1)) {
cat("Mode", i - 1, "\n")
if (is.vector(x[[i]], mode = "numeric")) {
print(x[[i]])
} else print(x[[i]]$tacvf)
}
} else {
cat("Printing the tacvf of", x[[1]], "\n")
if (is.vector(x[[2]], mode = "numeric")) {
print(x[[2]])
} else print(x[[2]]$tacvf)
}
invisible(x)
}
#' Plots the output from a call to \code{tacvf}
#'
#' Plots the theoretical autocovariance functions of the modes for a fitted
#' \code{arfima} object
#'
#' Only plots up to nine tacvfs. It is highly recommended that the
#' \code{arfima} object be weeded before calling \code{tacvf}
#'
#' @param x A \code{tacvf} object from a call to said function
#' @param type See \code{\link{plot}}. The default is recommended for short
#' \code{maxlag}
#' @param pch See \code{\link{plot}}
#' @param xlab See \code{\link{plot}}
#' @param ylab See \code{\link{plot}}
#' @param main See \code{\link{plot}}
#' @param xlim See \code{\link{plot}}
#' @param ylim See \code{\link{plot}}
#' @param tacf If \code{TRUE}, plots the theoretical autocorellations instead
#' @param maxlag The maximum lag for the plot
#' @param lag0 Whether or not to plot lag 0 of the tacvfs/tacfs. Default
#' \code{!tacf}. Used by \code{\link{tacfplot}}.
#' @param \dots Currently not used
#' @return None. There is a plot as output.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{tacvf}}
#' @references Veenstra, J.Q. Persistence and Antipersistence: Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#' set.seed(1234)
#' sim <- arfima.sim(1000, model = list(theta = 0.99, dfrac = 0.49))
#' fit <- arfima(sim, order = c(0, 0, 1))
#' plot(tacvf(fit))
#' plot(tacvf(fit), tacf = TRUE)
#'
plot.tacvf <- function(x, type = "o", pch = 20, xlab = NULL, ylab = NULL, main = NULL, xlim = NULL,
ylim = NULL, tacf = FALSE, maxlag = NULL, lag0 = !tacf, ...) {
if (is.null(ylab))
ylab <- if (!tacf)
"tacvf" else "tacf"
if (is.null(xlab))
xlab <- "lag"
m <- length(x) - 1
for (i in 2:(m + 1)) {
if (!is.vector(x[[i]], mode = "numeric"))
x[[i]] <- x[[i]]$tacvf
if (tacf)
x[[i]] <- x[[i]]/x[[i]][1]
}
if (m > 9) {
stop("more than nine modes currently not supported: please use the weed() function")
}
if (is.null(main)) {
if (!tacf) {
if (m > 1)
main <- paste("The tacvfs of ", x[[1]], sep = "") else main <- paste("The tacvf of ", x[[1]], sep = "")
} else {
if (m > 1)
main <- paste("The tacfs of ", x[[1]], sep = "") else main <- paste("The tacf of ", x[[1]], sep = "")
}
}
if (!lag0 && is.null(xlim)) {
xbds <- if (is.null(maxlag))
1:(length(x[[2]]) - 1) else 1:maxlag
xlim <- range(xbds)
} else if (is.null(xlim)) {
xbds <- if (is.null(maxlag))
0:(length(x[[2]]) - 1) else 0:maxlag
xlim <- range(xbds)
} else {
if (xlim[1] > xlim[2]) {
warning("upper x limit lower than lower x limit: resetting")
temp <- xlim[1]
xlim[1] <- xlim[2]
xlim[2] <- temp
}
if (xlim[2] > length(x[[2]]) - 1) {
warning("the bounds of the x limits (lag) end at n - 1 + xmaxlag; subtracting one from each of the limits")
xlim <- xlim - 1
}
if (xlim[1] < 0) {
warning("lower x limit (lag) less than 0: resetting to 0")
xlim[1] <- 0
}
xbds <- seq(xlim[1], xlim[2])
}
if (is.null(ylim)) {
minner <- Inf
maxxer <- -Inf
for (i in 1:m) {
minner <- min(minner, x[[i + 1]][xbds + 1])
maxxer <- max(maxxer, x[[i + 1]][xbds + 1])
}
ylim <- c(minner, maxxer)
} else if (ylim[1] > ylim[2]) {
warning("upper y limit lower than lower y limit: resetting")
temp <- ylim[1]
ylim[1] <- ylim[2]
ylim[2] <- temp
}
cols <- c("black", "red", "blue", "orange", "gray", "purple", "navy", "gold", "cyan")
if (m > 1) {
plot(x = xbds, y = x[[2]][xbds + 1], main = main, ylab = ylab, xlab = xlab, ylim = ylim,
xlim = xlim, type = type, pch = pch, col = cols[1])
for (i in 1:(m - 1)) lines(x = xbds, y = x[[i + 2]][xbds + 1], type = type, pch = pch,
col = cols[i + 1])
nncol <- floor(m/3)
if (nncol == 0)
nncol <- 1
legend(x = "topright", ncol = nncol, legend = paste("Mode", 1:m), pch = pch, col = cols[1:m])
} else plot(x = xbds, y = x[[2]][xbds + 1], main = main, ylab = ylab, xlab = xlab, ylim = ylim,
xlim = xlim, type = type, pch = pch, col = cols[1])
}
#' Plots the theoretical autocorralation functions (tacfs) of one or more fits.
#'
#' Plots the theoretical autocorralation functions (tacfs) of one or more fits.
#'
#'
#' @param fits A list of objects of class "arfima".
#' @param modes Either "all" or a vector of the same length as fits for which
#' the tacfs will be ploted.
#' @param xlab Optional. Usually better to be generated by the function.
#' @param ylab Optional. Usually better to be generated by the function.
#' @param main Optional. Usually better to be generated by the function.
#' @param xlim Optional. Usually better to be generated by the function.
#' @param ylim Optional. Usually better to be generated by the function.
#' @param maxlag Optional. Used to limit the length of tacfs. Highly
#' recommended to be a value from 20 - 50.
#' @param lag0 Whether or not the lag 0 tacf should be printed. Since this is
#' always 1 for all tacfs, recommended to be \code{TRUE}. It is easier to see
#' the shape of the tacfs.
#' @param \dots Optional. Currently not used.
#' @return NULL. However, there is a plot output.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{tacvf}}, \code{\link{plot.tacvf}}
#' @references Veenstra, J.Q. Persistence and Antipersistence: Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#' set.seed(34577)
#' sim <- arfima.sim(500, model = list(theta = 0.9, phi = 0.5, dfrac = 0.4))
#' fit1 <- arfima(sim, order = c(1, 0, 1), cpus = 2, back=TRUE)
#' fit2 <- arfima(sim, order = c(1, 0, 1), cpus = 2, lmodel = "g", back=TRUE)
#' fit3 <- arfima(sim, order = c(1, 0, 1), cpus = 2, lmodel = "h", back=TRUE)
#' fit1
#' fit2
#' fit3
#' tacfplot(fits = list(fit1, fit2, fit3), maxlag = 30)
#'
#' @export tacfplot
tacfplot <- function(fits = list(), modes = "all", xlab = NULL, ylab = NULL, main = NULL,
xlim = NULL, ylim = NULL, maxlag = 20, lag0 = FALSE, ...) {
if (length(fits) == 0)
stop("need at least one tacvf to compare")
if (length(fits) == 1) {
stop("for a single fit, only the fit should be provided: recommended to use plot.tacvf with tacf = TRUE")
}
if (inherits(fits, "arfima")) {
if (modes != "all") {
warning("only one fit provided and modes != 'all': using all modes")
}
warning("tacfplot is intended to compare two fits on the same data: recommended to use plot.tacvf with tacf = TRUE")
tacvf <- tacvf(fits)
if (is.null(main))
main <- "The tacf plot of the fit"
plot(tacvf, type = "o", xlab = xlab, ylab = ylab, main = main, xlim = xlim, ylim = ylim,
tacf = TRUE, maxlag = maxlag, lag0 = lag0, ...)
} else {
if (modes != "all" && length(modes) != length(fits))
stop("length of 'modes' not equal to length of 'fits', and modes != 'all'")
numfits <- length(fits)
if (numfits > 9)
stop("comparing more than 9 fits: stopping")
z1 <- fits[[1]]$z
num <- length(fits[[1]]$modes)
for (i in 2:numfits) {
z2 <- fits[[i]]$z
if (length(z1) != length(z2) || !all(z1 == z2))
stop("fits are not on the same series")
}
lenner <- length(tacvf(fits[[1]])[[1]])
if (!lag0 && is.null(xlim)) {
xbds <- if (is.null(maxlag))
1:(lenner - 1) else 1:maxlag
xlim <- range(xbds)
} else if (is.null(xlim)) {
xbds <- if (is.null(maxlag))
0:(lenner - 1) else 0:maxlag
xlim <- range(xbds)
} else {
if (xlim[1] > xlim[2]) {
warning("upper x limit lower than lower x limit: resetting")
temp <- xlim[1]
xlim[1] <- xlim[2]
xlim[2] <- temp
}
if (xlim[2] > length(z1) - 1) {
warning("the bounds of the x limits (lag) end at n - 1; subtracting one from each of the limits")
xlim <- xlim - 1
}
if (xlim[1] < 0) {
warning("lower x limit (lag) less than 0: resetting to 0")
xlim[1] <- 0
}
xbds <- seq(xlim[1], xlim[2])
}
tacfs <- vector("list", numfits)
maxnum <- 1
for (i in 1:numfits) {
if (modes == "all") {
maxnum <- max(maxnum, length(fits[[i]]$modes))
tacfs[[i]] <- tacvf(fits[[i]])
tacfs[[i]] <- tacfs[[i]][-1]
for (j in 1:length(tacfs[[i]])) {
num <- num + 1
tacfs[[i]][[j]] <- tacfs[[i]][[j]]/tacfs[[i]][[j]][1]
}
} else {
tacfs[[i]] <- tacfs[[i]][[modes[i]]]/tacfs[[i]][[modes[i]]][1]
}
}
if (is.null(main)) {
main <- paste("Comparing tacfs of ", numfits, " fits", sep = "")
}
namer <- paste("Fit ", 1:numfits)
if (num > 9)
warning("more than 9 tacfs to plot: plot will be very hard to read")
if (is.null(xlab))
xlab <- "lag"
if (is.null(ylab))
ylab <- "tacf"
if (is.null(ylim)) {
minner <- Inf
maxxer <- -Inf
for (i in 1:numfits) {
if (modes == "all") {
for (j in 1:length(tacfs[[i]])) {
minner <- min(minner, tacfs[[i]][[j]][xbds + 1])
maxxer <- max(maxxer, tacfs[[i]][[j]][xbds + 1])
}
} else {
minner <- min(minner, tacfs[[i]][xbds + 1])
maxxer <- max(maxxer, tacfs[[i]][xbds + 1])
}
}
ylim <- c(minner, maxxer)
} else if (ylim[1] > ylim[2]) {
warning("upper y limit lower than lower y limit: resetting")
temp <- ylim[1]
ylim[1] <- ylim[2]
ylim[2] <- temp
}
cols <- c("black", "red", "blue", "orange", "gray", "purple", "navy", "gold", "cyan")
pchlist <- c(16, 17, 18, 0, 1, 2, 3, 4, 5)
pchs <- pchlist[1:maxnum]
if (modes != "all") {
plot(x = xbds, y = tacfs[[1]][xbds + 1], main = main, ylab = ylab, xlab = xlab,
ylim = ylim, xlim = xlim, type = "o", pch = pchs[modes[1]], col = cols[1])
for (i in 2:(numfits)) lines(x = xbds, y = tacfs[[i]][xbds + 1], type = "o",
pch = pchs[modes[i]], col = cols[i])
legend(x = "topright", ncol = floor(numfits/3), legend = paste(namer, "mode",
modes), pch = pchs[modes], col = cols[1:numfits])
} else {
plot(x = xbds, y = tacfs[[1]][[1]][xbds + 1], main = main, ylab = ylab, xlab = xlab,
ylim = ylim, xlim = xlim, type = "o", pch = pchs[1], col = cols[1])
if (length(tacfs[[1]]) > 1) {
for (j in 2:length(tacfs[[1]])) lines(x = xbds, y = tacfs[[1]][[j]][xbds +
1], type = "o", pch = pchs[j], col = cols[1])
}
for (i in 2:numfits) {
for (j in 1:length(tacfs[[i]])) {
lines(x = xbds, y = tacfs[[i]][[j]][xbds + 1], type = "o", pch = pchs[j],
col = cols[i])
}
}
legend(x = "topright", ncol = floor(num/3), legend = c(namer, paste("Mode",
1:maxnum)), pch = c(rep(19, numfits), pchs), col = c(cols[1:numfits], rep("black",
maxnum)))
}
}
}
#' The theoretical autocovariance function of a long memory process.
#'
#' Calculates the tacvf of a mixed long memory-ARMA (with posible seasonal
#' components). Combines long memory and ARMA (and non-seasonal and seasonal)
#' parts via convolution.
#'
#' The log-likelihood is computed for the given series z and the parameters.
#' If two or more of \code{dfrac}, \code{H} or \code{alpha} are present and/or
#' two or more of \code{dfs}, \code{Hs} or \code{alphas} are present, an error
#' will be thrown, as otherwise there is redundancy in the model. Note that
#' non-seasonal and seasonal components can be of different types: for example,
#' there can be seasonal FGN with FDWN at the non-seasonal level.
#'
#' The moving average parameters are in the Box-Jenkins convention: they are
#' the negative of the parameters given by \code{\link{arima}}.
#'
#' @param phi The autoregressive parameters in vector form.
#' @param theta The moving average parameters in vector form. See Details for
#' differences from \code{\link{arima}}.
#' @param dfrac The fractional differencing parameter.
#' @param phiseas The seasonal autoregressive parameters in vector form.
#' @param thetaseas The seasonal moving average parameters in vector form. See
#' Details for differences from \code{\link{arima}}.
#' @param dfs The seasonal fractional differencing parameter.
#' @param H The Hurst parameter for fractional Gaussian noise (FGN). Should
#' not be mixed with \code{dfrac} or \code{alpha}: see "Details".
#' @param Hs The Hurst parameter for seasonal fractional Gaussian noise (FGN).
#' Should not be mixed with \code{dfs} or \code{alphas}: see "Details".
#' @param alpha The decay parameter for power-law autocovariance (PLA) noise.
#' Should not be mixed with \code{dfrac} or \code{H}: see "Details".
#' @param alphas The decay parameter for seasonal power-law autocovariance
#' (PLA) noise. Should not be mixed with \code{dfs} or \code{Hs}: see
#' "Details".
#' @param period The periodicity of the seasonal components. Must be >= 2.
#' @param maxlag The number of terms to compute: technically the output
#' sequence is from lags 0 to maxlag, so there are maxlag + 1 terms.
#' @param useCt Whether or not to use C to compute the (parts of the) tacvf.
#' @param sigma2 Used in \code{\link{arfima.sim}}: determines the value of the
#' innovation variance. The tacvf sequence is multiplied by this value.
#' @return A sequence of length maxlag + 1 (lags 0 to maxlag) of the tacvf of
#' the given process.
#' @author JQ (Justin) Veenstra and A. I. McLeod
#' @references Veenstra, J.Q. Persistence and Antipersistence: Theory and
#' Software (PhD Thesis)
#'
#' P. Borwein (1995) An efficient algorithm for Riemann Zeta function Canadian
#' Math. Soc. Conf. Proc., 27, pp. 29-34.
#' @keywords ts
#' @examples
#'
#' t1 <- tacvfARFIMA(phi = c(0.2, 0.1), theta = 0.4, dfrac = 0.3, maxlag = 30)
#' t2 <- tacvfARFIMA(phi = c(0.2, 0.1), theta = 0.4, H = 0.8, maxlag = 30)
#' t3 <- tacvfARFIMA(phi = c(0.2, 0.1), theta = 0.4, alpha = 0.4, maxlag = 30)
#' plot(t1, type = "o", col = "blue", pch = 20)
#' lines(t2, type = "o", col = "red", pch = 20)
#' lines(t3, type = "o", col = "purple", pch = 20) #they decay at about the same rate
#'
#'
#' @export tacvfARFIMA
"tacvfARFIMA" <- function(phi = numeric(0), theta = numeric(0), dfrac = numeric(0), phiseas = numeric(0),
thetaseas = numeric(0), dfs = numeric(0), H = numeric(0), Hs = numeric(0), alpha = numeric(0),
alphas = numeric(0), period = 0, maxlag, useCt = T, sigma2 = 1) {
if (length(maxlag) == 0)
stop("maxlag must be defined")
if (!IdentInvertQ(phi = phi, theta = theta, dfrac = dfrac, phiseas = phiseas, thetaseas = thetaseas,
dfs = dfs, H = H, Hs = Hs, alpha = alpha, alphas = alphas, period = period, ident = F)) {
#print("oops")
return(NULL)
}
if (length(period) == 0 || is.na(period) || period == 0) {
if (length(phiseas) > 0 && phiseas != 0)
stop("Period = 0, but seasonal phi != 0")
if (length(thetaseas) > 0 && thetaseas != 0)
stop("Period = 0, but seasonal theta != 0")
if ((length(dfs) > 0 && dfs != 0) || (length(Hs) > 0 && Hs != 0) || (length(alphas) >
0 && alphas != 0))
stop("Period = 0, but fractional seasonal parameter != 0")
return(tacvfFARMA(phi = phi, theta = theta, dfrac = dfrac, H = H, alpha = alpha,
maxlag = maxlag, useCt = useCt, sigma2 = sigma2))
}
if ((period != round(period)) || (period < 2))
stop("Period must be an integer >= 2")
lagTrunc <- 2 * max(128, nextn(maxlag, factors = 2))
iseas <- nextn(period, factors = 2)
model <- tacvfFARMA(phi = phi, theta = theta, dfrac = dfrac, H = H, alpha = alpha, maxlag = (lagTrunc *
iseas), nolagtrunc = T, useCt = useCt, sigma2 = 1)
modelseas <- tacvfFARMA(phi = phiseas, theta = thetaseas, dfrac = dfs, H = Hs, alpha = alphas,
maxlag = (lagTrunc * 2), nolagtrunc = T, useCt = useCt, sigma2 = 1)
modelseas <- shift(modelseas, period, useCt = T)
if (length(modelseas) < length(model))
stop("oops")
modelseas <- modelseas[1:(lagTrunc * iseas + 1)]
z <- sigma2 * mix(model, modelseas)
return(z[1:(maxlag + 1)])
}
"tacvfFARMA" <- function(phi = numeric(0), theta = numeric(0), dfrac = numeric(0), H = numeric(0),
alpha = numeric(0), maxlag, nolagtrunc = F, useCt = T, sigma2 = 1) {
if (length(H) + length(dfrac) + length(alpha) > 1)
stop("more than one GHD process specified: stopping")
if ((length(H) == 0) && (length(dfrac) == 0) && (length(alpha) == 0) && (length(phi) ==
0) && (length(theta) == 0)) {
return(c(sigma2, rep(0, maxlag)))
}
if ((length(dfrac) == 0) && (length(H) == 0) && (length(alpha) == 0)) {
return(tacvfARMA(phi = phi, theta = theta, maxlag = maxlag, useCt = useCt, sigma2 = sigma2))
}
onlyHD <- (length(phi) == 0 && length(theta) == 0)
if (nolagtrunc || onlyHD)
lagTrunc <- maxlag else {
lagTrunc <- 2 * max(128, nextn(maxlag, factors = 2))
}
if (length(dfrac) > 0)
x <- tacvfFDWN(dfrac = dfrac, maxlag = lagTrunc, useCt = useCt) else if (length(H) > 0)
x <- tacvfFGN(H = H, maxlag = lagTrunc, useCt = useCt) else x <- tacvfHD(alpha = alpha, maxlag = lagTrunc, useCt = useCt)
if (onlyHD)
return(sigma2 * x[1:(maxlag + 1)])
y <- tacvfARMA(phi = phi, theta = theta, maxlag = lagTrunc, useCt = useCt, sigma2 = 1)
z <- sigma2 * mix(x, y)
return(z[1:(maxlag + 1)])
}
"symtacvf" <- function(x) {
c(rev(x[-1])[-1], x)
}
"tacvfARMA" <- function(phi = numeric(0), theta = numeric(0), maxlag = 20, useCt = T, sigma2 = 1) {
if (useCt) {
dd <- function(x, y, ml) .C("tacvfARMA_C", as.double(x), as.integer(length(x)),
as.double(y), as.integer(length(y)), as.integer(ml), res = double(ml + 1))$res
return(sigma2 * dd(phi, theta, maxlag))
}
p <- length(phi)
q <- length(theta)
maxlagp1 <- maxlag + 1
if (max(p, q) == 0) {
return(c(sigma2, numeric(maxlag)))
}
r <- max(p, q) + 1
b <- numeric(r)
C <- numeric(q + 1)
C[1] <- 1
theta2 <- c(-1, theta)
phi2 <- numeric(3 * r)
phi2[r] <- -1
if (p > 0) {
phi2[r + 1:p] <- phi
}
if (q > 0) {
for (k in 1:q) {
C[k + 1] <- -theta[k]
if (p > 0) {
for (i in 1:min(p, k)) {
C[k + 1] <- C[k + 1] + phi[i] * C[k + 1 - i]
}
}
}
}
for (k in 0:q) {
for (i in k:q) {
b[k + 1] <- b[k + 1] - theta2[i + 1] * C[i - k + 1]
}
}
if (p == 0) {
g <- c(b, numeric(maxlagp1))[1:maxlagp1]
return(g)
} else if (p > 0) {
a <- matrix(numeric(r^2), ncol = r)
for (i in 1:r) {
for (j in 1:r) {
if (j == 1) {
a[i, j] <- phi2[r + i - 1]
} else if (j != 1) {
a[i, j] <- phi2[r + i - j] + phi2[r + i + j - 2]
}
}
}
g <- solve(a, -b)
if (length(g) <= maxlag) {
g <- c(g, numeric(maxlagp1 - r))
for (i in (r + 1):maxlagp1) {
g[i] <- phi %*% g[i - 1:p]
}
return(sigma2 * g[1:maxlagp1])
} else if (length(g) >= maxlagp1) {
return(sigma2 * g[1:maxlagp1])
}
}
}
"tacvfFDWN" <- function(dfrac, maxlag, useCt = T) {
if (!InvertibleD(dfrac)) {
warning("Model is non-causal or non-invertible\n")
return(NULL)
}
if (useCt) {
ta <- function(ds, ma) .C("tacvfFDWN_C", as.double(ds), as.integer(ma), x = double(ma +
1))$x
return(ta(dfrac, maxlag))
}
x <- numeric(maxlag + 1)
x[1] <- gamma(1 - 2 * dfrac)/gamma(1 - dfrac)^2
for (i in 1:maxlag) {
x[i + 1] <- ((i - 1 + dfrac)/(i - dfrac)) * x[i]
}
return(x)
}
tacvfFGN <- function(H, maxlag, useCt = TRUE) {
if (!InvertibleH(H))
return(NULL)
if (useCt) {
tg <- function(H, ma) .C("tacfFGN_C", as.double(H), as.integer(ma), x = double(ma +
1))$x
return(tg(H, maxlag))
}
h2 <- 2 * H
r <- sapply(0:maxlag, function(k) 0.5 * (abs(k + 1)^h2 - 2 * abs(k)^h2 + abs(k - 1)^h2))
return(r)
}
## will lgamma make it faster??
Zeta <- function(s, n = 20) {
d <- rep(0, n + 1)
d[1] <- n * (factorial(n - 1))/(factorial(n))
for (i in 1:(n)) {
d[i + 1] <- d[i] + n * factorial(n + i - 1) * 4^i/(factorial(n - i) * factorial(2 *
i))
}
zeta <- 0
for (i in 0:(n - 1)) {
zeta <- zeta + (-1)^i * (d[i + 1] - d[n + 1])/(i + 1)^s
}
zeta <- zeta * (-1)/(d[n + 1] * (1 - 2^(1 - s)))
zeta
}
## Reference paper. (zetaalgorithm.)
tacvfHD <- function(alpha, maxlag, useCt = TRUE) {
if (!InvertibleAlpha(alpha))
return(NULL)
if (useCt) {
th <- function(alpha, ma) .C("tacfHD_C", as.double(alpha), as.integer(ma), x = double(ma +
1))$x
return(th(alpha, maxlag))
}
val <- (-2 * Zeta(alpha))^(-1)
r <- c(1, val * sapply(1:maxlag, function(k) k^(-alpha)))
return(r)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.