"iARMAapprox" <- function(phi = numeric(0), theta = numeric(0), seas = F) {
maxlag <- 2048
if (!(InvertibleQ(phi) & InvertibleQ(theta))) {
return(NULL)
}
p <- length(phi)
q <- length(theta)
unames <- character(0)
vnames <- character(0)
if (p > 0) {
psisphi <- psiwtsAR(phi, maxlag = maxlag)
vv <- numeric(p)
vv[1] <- crossprod(psisphi)
if (p > 1)
vv[2:p] <- sapply(1:(p - 1), function(i) crossprod(psisphi[-(1:i)], psisphi[-((maxlag -
i + 1):maxlag)]))
vvmatrix <- toeplitz(vv)
imatrix <- vvmatrix
if (!seas)
vnames <- paste("phi(", 1:p, ")", sep = "") else vnames <- paste("phiseas(", 1:p, ")", sep = "")
}
if (q > 0) {
psisth <- -psiwtsAR(theta, maxlag = maxlag)
uu <- numeric(q)
uu[1] <- crossprod(psisth)
if (q > 1)
uu[2:q] <- sapply(1:(q - 1), function(i) crossprod(psisth[-(1:i)], psisth[-((maxlag -
i + 1):maxlag)]))
uumatrix <- toeplitz(uu)
imatrix <- uumatrix
if (!seas)
unames <- paste("theta(", 1:q, ")", sep = "") else unames <- paste("thetaseas(", 1:q, ")", sep = "")
}
if (p > 0 && q > 0) {
uvmatrix <- matrix(numeric(1), nrow = p, ncol = q)
uvmatrix[1, 1] <- crossprod(psisphi, psisth)
if (q > 1)
uvmatrix[1, -1] <- sapply(1:(q - 1), function(i) crossprod(psisphi[-(1:i)],
psisth[-((maxlag - i + 1):maxlag)]))
if (p > 1)
uvmatrix[-1, 1] <- sapply(1:(p - 1), function(i) crossprod(psisth[-(1:i)], psisphi[-((maxlag -
i + 1):maxlag)]))
if (p > 1 && q > 1) {
for (i in 2:p) {
for (j in 2:q) {
if (i == j)
uvmatrix[i, j] <- uvmatrix[1, 1] else if ((k = i - j) < 0)
uvmatrix[i, j] <- crossprod(psisphi[-(1:(-k))], psisth[-((maxlag + k +
1):maxlag)]) else uvmatrix[i, j] <- crossprod(psisth[-(1:k)], psisphi[-((maxlag - k +
1):maxlag)])
}
}
}
imatrix <- cbind(rbind(vvmatrix, t(uvmatrix)), rbind(uvmatrix, uumatrix))
}
inames <- c(vnames, unames)
dimnames(imatrix) <- list(inames, inames)
return(imatrix)
}
"iARMASeasMultapprox" <- function(phi = numeric(0), theta = numeric(0), phiseas = numeric(0),
thetaseas = numeric(0), period = 0) {
if (period == 0) {
warning("Seasonal call, but period = 0; reverting to nonseasonal model")
return(iARMAapprox(phi = phi, theta = theta))
}
if (!(InvertibleQ(phi) & InvertibleQ(theta) & InvertibleQ(phiseas) & InvertibleQ(thetaseas))) {
return(NULL)
}
maxlag <- 2048
maxlag1 <- maxlag * period
imatrix <- iarma <- isarma <- block <- ppsmatrix <- qpsmatrix <- pqsmatrix <- qqsmatrix <- NULL
p <- length(phi)
q <- length(theta)
ps <- length(phiseas)
qs <- length(thetaseas)
if (p + q > 0)
iarma <- iARMAapprox(phi = phi, theta = theta)
if (ps + qs > 0)
isarma <- iARMAapprox(phi = phiseas, theta = thetaseas, seas = T)
inames <- c(dimnames(iarma)[[1]], dimnames(isarma)[[1]])
if (p > 0)
psiph <- psiwtsAR(phi, maxlag1)
if (q > 0)
psith <- -psiwtsAR(theta, maxlag1)
if (ps > 0)
psips <- shift(psiwtsAR(phiseas, maxlag), period)
if (qs > 0)
psits <- -shift(psiwtsAR(thetaseas, maxlag), period)
if (p > 0 && ps > 0) {
ppsmatrix <- matrix(0, nrow = p, ncol = ps)
ppsmatrix[1, 1] <- crossprod(psiph[-(1:(period - 1))], psips[-((maxlag1 - period +
2):maxlag1)])
if (ps > 1)
ppsmatrix[1, -1] <- sapply(2:ps, function(i) crossprod(psiph[-(1:(period * i -
1))], psips[-((maxlag1 - i * period + 2):maxlag1)]))
if (p > 1)
ppsmatrix[-1, 1] <- sapply(2:p, function(i) crossprod(psiph[-(1:(period + i -
1))], psips[-((maxlag1 - i - period + 2):maxlag1)]))
if (p > 1 && ps > 1) {
for (i in 2:p) {
for (j in 2:ps) {
ppsmatrix[i, j] <- crossprod(psiph[-(1:(period * j + i - 1))], psips[-((maxlag1 -
i - period * j + 2):maxlag1)])
}
}
}
}
if (q > 0 && ps > 0) {
qpsmatrix <- matrix(0, nrow = q, ncol = ps)
qpsmatrix[1, 1] <- crossprod(psith[-(1:(period - 1))], psips[-((maxlag1 - period +
2):maxlag1)])
if (ps > 1)
qpsmatrix[1, -1] <- sapply(2:ps, function(i) crossprod(psith[-(1:(period * i -
1))], psips[-((maxlag1 - i * period + 2):maxlag1)]))
if (q > 1)
qpsmatrix[-1, 1] <- sapply(2:q, function(i) crossprod(psith[-(1:(period + i -
1))], psips[-((maxlag1 - i - period + 2):maxlag1)]))
if (q > 1 && ps > 1) {
for (i in 2:q) {
for (j in 2:ps) {
qpsmatrix[i, j] <- crossprod(psith[-(1:(period * j + i - 1))], psips[-((maxlag1 -
i - period * j + 2):maxlag1)])
}
}
}
}
if (p > 0 && qs > 0) {
pqsmatrix <- matrix(0, nrow = p, ncol = qs)
pqsmatrix[1, 1] <- crossprod(psiph[-(1:(period - 1))], psits[-((maxlag1 - period +
2):maxlag1)])
if (qs > 1)
pqsmatrix[1, -1] <- sapply(2:qs, function(i) crossprod(psiph[-(1:(period * i -
1))], psits[-((maxlag1 - i * period + 2):maxlag1)]))
if (p > 1)
pqsmatrix[-1, 1] <- sapply(2:p, function(i) crossprod(psiph[-(1:(period + i -
1))], psits[-((maxlag1 - i - period + 2):maxlag1)]))
if (p > 1 && qs > 1) {
for (i in 2:p) {
for (j in 2:qs) {
pqsmatrix[i, j] <- crossprod(psiph[-(1:(period * j + i - 1))], psits[-((maxlag1 -
i - period * j + 2):maxlag1)])
}
}
}
}
if (q > 0 && qs > 0) {
qqsmatrix <- matrix(0, nrow = q, ncol = qs)
qqsmatrix[1, 1] <- crossprod(psith[-(1:(period - 1))], psits[-((maxlag1 - period +
2):maxlag1)])
if (qs > 1)
pqsmatrix[1, -1] <- sapply(2:qs, function(i) crossprod(psith[-(1:(period * i -
1))], psits[-((maxlag1 - i * period + 2):maxlag1)]))
if (q > 1)
qqsmatrix[-1, 1] <- sapply(2:q, function(i) crossprod(psith[-(1:(period + i -
1))], psits[-((maxlag1 - i - period + 2):maxlag1)]))
if (q > 1 && qs > 1) {
for (i in 2:q) {
for (j in 2:qs) {
qqsmatrix[i, j] <- crossprod(psith[-(1:(period * j + i - 1))], psits[-((maxlag1 -
i - period * j + 2):maxlag1)])
}
}
}
}
block <- rbind(cbind(ppsmatrix, pqsmatrix), cbind(qpsmatrix, qqsmatrix))
imatrix <- cbind(iarma, block)
imatrix <- rbind(imatrix, cbind(if (!is.null(block))
t(block) else NULL, isarma))
dimnames(imatrix) <- list(inames, inames)
return(imatrix)
}
"iARMASeasMult" <- function(phi = numeric(0), theta = numeric(0), phiseas = numeric(0),
thetaseas = numeric(0), period = 0, warn = T) {
if (period == 0) {
if (warn)
warning("Seasonal call, but period = 0; reverting to nonseasonal model")
return(iARMA(phi = phi, theta = theta))
}
if (!(InvertibleQ(phi) & InvertibleQ(theta) & InvertibleQ(phiseas) & InvertibleQ(thetaseas))) {
warning("Model is non-causal or non-invertible\n")
return(NULL)
}
imatrix <- iarma <- isarma <- block <- ppsmatrix <- qpsmatrix <- pqsmatrix <- qqsmatrix <- NULL
p <- length(phi)
q <- length(theta)
ps <- length(phiseas)
qs <- length(thetaseas)
if (p + q > 0)
iarma <- iARMA(phi = phi, theta = theta)
if (ps + qs > 0)
isarma <- iARMA(phi = phiseas, theta = thetaseas, period = period)
inames <- c(dimnames(iarma)[[1]], dimnames(isarma)[[1]])
phi.s <- shift(c(1, phiseas), period)[-1]
theta.s <- shift(c(1, thetaseas), period)[-1]
if (p > 0 && ps > 0) {
ppsmatrix <- matrix(0, nrow = p, ncol = ps)
tt <- tccfAR(phi, phi.s)
for (i in 1:p) for (j in 1:ps) ppsmatrix[i, j] <- tt[ps * period + (i - 0) - (j -
0) * period]
}
if (q > 0 && ps > 0) {
qpsmatrix <- matrix(0, nrow = q, ncol = ps)
tt <- -tccfAR(theta, phi.s)
for (i in 1:q) for (j in 1:ps) qpsmatrix[i, j] <- tt[ps * period + (i - 0) - (j -
0) * period]
}
if (p > 0 && qs > 0) {
pqsmatrix <- matrix(0, nrow = p, ncol = qs)
tt <- -tccfAR(phi, theta.s)
for (i in 1:p) for (j in 1:qs) pqsmatrix[i, j] <- tt[qs * period + (i - 0) - (j -
0) * period]
}
if (q > 0 && qs > 0) {
qqsmatrix <- matrix(0, nrow = q, ncol = qs)
tt <- tccfAR(theta, theta.s)
for (i in 1:q) for (j in 1:qs) qqsmatrix[i, j] <- tt[qs * period + (i - 0) - (j -
0) * period]
}
block <- rbind(cbind(ppsmatrix, pqsmatrix), cbind(qpsmatrix, qqsmatrix))
imatrix <- cbind(iarma, block)
imatrix <- rbind(imatrix, cbind(if (!is.null(block))
t(block) else NULL, isarma))
dimnames(imatrix) <- list(inames, inames)
return(imatrix)
}
"iFARMA" <- function(phi = numeric(0), theta = numeric(0)) {
if (!(InvertibleQ(phi) & InvertibleQ(theta))) {
warning("Model is non-causal or non-invertible\n")
return(NULL)
}
I22 <- (pi^2)/6
if ((length(phi) == 0) && (length(theta) == 0)) {
I22 <- as.matrix(I22)
dimnames(I22) <- list("d.f", "d.f")
return(I22)
}
I11 <- iARMA(phi = phi, theta = theta)
J11 <- numeric(0)
J12 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi)
if (length(theta) > 0)
J12 <- -jFARMA(theta)
J <- c(J11, J12)
I <- rbind(I11, J)
J <- c(J, I22)
I <- cbind(I, J)
inames <- c(dimnames(I11)[[1]], "d.f")
dimnames(I) <- list(inames, inames)
return(I)
}
"iFARMAapprox" <- function(phi = numeric(0), theta = numeric(0)) {
if (!(InvertibleQ(phi) & InvertibleQ(theta))) {
warning("Model is non-causal or non-invertible\n")
return(NULL)
}
I22 <- (pi^2)/6
if ((length(phi) == 0) && (length(theta) == 0)) {
I22 <- as.matrix(I22)
dimnames(I22) <- list("d.f", "d.f")
return(I22)
}
I11 <- iARMAapprox(phi = phi, theta = theta)
J11 <- numeric(0)
J12 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi)
if (length(theta) > 0)
J12 <- -jFARMA(theta)
J <- c(J11, J12)
I <- rbind(I11, J)
J <- c(J, I22)
I <- cbind(I, J)
inames <- c(dimnames(I11)[[1]], "d.f")
dimnames(I) <- list(inames, inames)
return(I)
}
# write about below
#' The Fisher information matrix of an ARFIMA process
#'
#' Computes the approximate or (almost) exact Fisher information matrix of an
#' ARFIMA process
#'
#' The matrices are calculated as outlined in Veenstra and McLeod (2012), which
#' draws on many references. The psi-weights approximation has a fixed maximum
#' lag for the weights as 2048 (to be changed to be adaptable.) The fractional
#' difference(s) by AR/MA components have a fixed maximum lag of 256, also to
#' be changed. Thus the exact matrix has some approximation to it. Also note
#' that the approximate method takes much longer than the "exact" one.
#'
#' 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 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 period The periodicity of the seasonal components. Must be >= 2.
#' @param dfrac TRUE if we include the fractional d parameter, FALSE otherwise
#' @param dfs TRUE if we include the seasonal fractional d parameter, FALSE
#' otherwise
#' @param exact If FALSE, calculate the approximate information matrix via
#' psi-weights. Otherwise the (almost) exact information matrix will be
#' calculated. See "Details".
#' @return The information matrix of the model.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{IdentInvertQ}}
#' @references Veenstra, J.Q. Persistence and Antipersistence: Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#' tick <- proc.time()
#' exactI <- iARFIMA(phi = c(.4, -.2), theta = c(.7), phiseas = c(.8, -.4),
#' d = TRUE, dfs = TRUE, period = 12)
#' proc.time() - tick
#' tick <- proc.time()
#' approxI <- iARFIMA(phi = c(.4, -.2), theta = c(.7), phiseas = c(.8, -.4),
#' d = TRUE, dfs = TRUE, period = 12, exact = FALSE)
#' proc.time() - tick
#' exactI
#' max(abs(exactI - approxI))
#'
#' @export iARFIMA
"iARFIMA" <- function(phi = numeric(0), theta = numeric(0), phiseas = numeric(0), thetaseas = numeric(0),
period = 0, dfrac = TRUE, dfs = FALSE, exact = TRUE) {
if (exact)
ans <- iARFIMAexact(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period, d = dfrac, dfs = dfs) else ans <- iARFIMAapprox(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period, d = dfrac, dfs = dfs)
ans
}
"iARFIMAexact" <- function(phi = numeric(0), theta = numeric(0), phiseas = numeric(0), thetaseas = numeric(0),
period = 0, d = TRUE, dfs = T) {
if (period == 0) {
if (dfs || length(phiseas) > 0 || length(thetaseas) > 0)
stop("0 period, but seasonal parts")
if (d)
return(iFARMA(phi, theta)) else {
if ((length(phi) == 0) && (length(theta) == 0))
return(numeric(0)) else return(iARMA(phi, theta))
}
}
if (d && dfs) {
I22 <- (pi^2)/6
I33 <- (pi^2)/(6 * period)
I4 <- cbind(c(I22, I33), c(I33, I22))
if ((length(phi) == 0) && (length(theta) == 0) && (length(phiseas) == 0) && (length(thetaseas) ==
0)) {
dimnames(I4) <- list(c("d.f", paste("d.f.", period, sep = "")), c("d.f", paste("d.f.",
period, sep = "")))
return(I4)
}
I11 <- iARMASeasMult(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
J11 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi) * -1
J12 <- numeric(0)
if (length(phiseas) > 0)
J12 <- jFARMA(phiseas, period = period, fs = F) * -1
if (length(theta) > 0)
J11 <- c(J11, jFARMA(theta))
if (length(thetaseas) > 0)
J12 <- c(J12, jFARMA(thetaseas, period = period, fs = F))
J <- c(J11, J12)
K11 <- numeric(0)
if (length(phi) > 0)
K11 <- jFARMA(phi, period = period, fs = T) * -1
K12 <- numeric(0)
if (length(phiseas) > 0)
K12 <- jFARMA(phiseas) * -1
if (length(theta) > 0)
K11 <- c(K11, jFARMA(theta, period = period, fs = T))
if (length(thetaseas) > 0)
K12 <- c(K12, jFARMA(thetaseas))
K <- c(K11, K12)
I <- rbind(I11, J, K)
I <- cbind(I, rbind(cbind(J, K), I4))
inames <- c(dimnames(I11)[[1]], "d.f", paste("d.f.", period, sep = ""))
dimnames(I) <- list(inames, inames)
return(I)
} else if (d) {
I22 <- (pi^2)/6
if ((length(phi) == 0) && (length(theta) == 0) && (length(phiseas) == 0) && (length(thetaseas) ==
0))
return(I22)
I11 <- iARMASeasMult(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
J11 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi) * -1
J12 <- numeric(0)
if (length(phiseas) > 0)
J12 <- jFARMA(phiseas, period = period, fs = F) * -1
if (length(theta) > 0)
J11 <- c(J11, jFARMA(theta))
if (length(thetaseas) > 0)
J12 <- c(J12, jFARMA(thetaseas, period = period, fs = F))
J <- c(J11, J12)
I <- rbind(I11, J)
I <- cbind(I, rbind(cbind(J), I22))
inames <- c(dimnames(I11)[[1]], "d.f")
dimnames(I) <- list(inames, inames)
return(I)
} else if (dfs) {
I22 <- (pi^2)/6
if ((length(phi) == 0) && (length(theta) == 0) && (length(phiseas) == 0) && (length(thetaseas) ==
0)) {
I22 <- as.matrix(I22)
dimnames(I22) <- list(paste("d.f.", period, sep = ""), paste("d.f.", period,
sep = ""))
return(I22)
}
I11 <- iARMASeasMult(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
J11 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi, period = period, fs = T) * -1
J12 <- numeric(0)
if (length(phiseas) > 0)
J12 <- jFARMA(phiseas) * -1
if (length(theta) > 0)
J11 <- c(J11, jFARMA(theta, period = period, fs = T))
if (length(thetaseas) > 0)
J12 <- c(J12, jFARMA(thetaseas))
J <- c(J11, J12)
I <- rbind(I11, J)
I <- cbind(I, rbind(cbind(J), I22))
inames <- c(dimnames(I11)[[1]], paste("d.f.", period, sep = ""))
dimnames(I) <- list(inames, inames)
return(I)
} else {
I <- iARMASeasMult(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
return(I)
}
}
"iARFIMAapprox" <- function(phi = numeric(0), theta = numeric(0), phiseas = numeric(0),
thetaseas = numeric(0), period = 0, d = T, dfs = T) {
if (period == 0) {
if (dfs || length(phiseas) > 0 || length(thetaseas) > 0)
stop("0 period, but seasonal parts")
if (d)
return(iFARMAapprox(phi, theta)) else {
if ((length(phi) == 0) && (length(theta) == 0))
return(numeric(0)) else return(iARMAapprox(phi, theta))
}
}
if (d && dfs) {
I22 <- (pi^2)/6
I33 <- (pi^2)/(6 * period)
I4 <- cbind(c(I22, I33), c(I33, I22))
if ((length(phi) == 0) && (length(theta) == 0) && (length(phiseas) == 0) && (length(thetaseas) ==
0)) {
dimnames(I4) <- list(c("d.f", paste("d.f.", period, sep = "")), c("d.f", paste("d.f.",
period, sep = "")))
return(I4)
}
I11 <- iARMASeasMultapprox(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
J11 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi) * -1
J12 <- numeric(0)
if (length(phiseas) > 0)
J12 <- jFARMA(phiseas, period = period, fs = F) * -1
if (length(theta) > 0)
J11 <- c(J11, jFARMA(theta))
if (length(thetaseas) > 0)
J12 <- c(J12, jFARMA(thetaseas, period = period, fs = F))
J <- c(J11, J12)
K11 <- numeric(0)
if (length(phi) > 0)
K11 <- jFARMA(phi, period = period, fs = T) * -1
K12 <- numeric(0)
if (length(phiseas) > 0)
K12 <- jFARMA(phiseas) * -1
if (length(theta) > 0)
K11 <- c(K11, jFARMA(theta, period = period, fs = T))
if (length(thetaseas) > 0)
K12 <- c(K12, jFARMA(thetaseas))
K <- c(K11, K12)
I <- rbind(I11, J, K)
I <- cbind(I, rbind(cbind(J, K), I4))
inames <- c(dimnames(I11)[[1]], "d.f", paste("d.f.", period, sep = ""))
dimnames(I) <- list(inames, inames)
return(I)
} else if (d) {
I22 <- (pi^2)/6
if ((length(phi) == 0) && (length(theta) == 0) && (length(phiseas) == 0) && (length(thetaseas) ==
0)) {
I22 <- as.matrix(I22)
dimnames(I22) <- list(paste("d.f.", period, sep = ""), paste("d.f.", period,
sep = ""))
return(I22)
}
I11 <- iARMASeasMultapprox(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
J11 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi) * -1
J12 <- numeric(0)
if (length(phiseas) > 0)
J12 <- jFARMA(phiseas, period = period, fs = F) * -1
if (length(theta) > 0)
J11 <- c(J11, jFARMA(theta))
if (length(thetaseas) > 0)
J12 <- c(J12, jFARMA(thetaseas, period = period, fs = F))
J <- c(J11, J12)
I <- rbind(I11, J)
I <- cbind(I, rbind(cbind(J), I22))
inames <- c(dimnames(I11)[[1]], "d.f")
dimnames(I) <- list(inames, inames)
return(I)
} else if (dfs) {
I22 <- (pi^2)/6
if ((length(phi) == 0) && (length(theta) == 0) && (length(phiseas) == 0) && (length(thetaseas) ==
0))
return(I22)
I11 <- iARMASeasMultapprox(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
J11 <- numeric(0)
if (length(phi) > 0)
J11 <- jFARMA(phi, period = period, fs = T) * -1
J12 <- numeric(0)
if (length(phiseas) > 0)
J12 <- jFARMA(phiseas) * -1
if (length(theta) > 0)
J11 <- c(J11, jFARMA(theta, period = period, fs = T))
if (length(thetaseas) > 0)
J12 <- c(J12, jFARMA(thetaseas))
J <- c(J11, J12)
I <- rbind(I11, J)
I <- cbind(I, rbind(cbind(J), I22))
inames <- c(dimnames(I11)[[1]], paste("d.f.", period, sep = ""))
dimnames(I) <- list(inames, inames)
return(I)
} else {
I <- iARMASeasMultapprox(phi = phi, theta = theta, phiseas = phiseas, thetaseas = thetaseas,
period = period)
return(I)
}
}
"jFARMA" <- function(theta, period = 1, fs = F, maxlag = 256) {
q <- length(theta)
J <- numeric(q)
maxlag <- maxlag * period
maxlagger <- 0:maxlag
psis <- psiwtsAR(theta, maxlag = maxlag)
if (fs) {
for (k in 1:q) {
a <- shift(psis[-(1:(period - k))], -period)
J[k] <- sum(a/(1 + maxlagger[1:length(a)]))
}
} else {
for (k in 1:q) {
J[k] <- sum(psis/((k + maxlagger) * period))
}
}
return(J)
}
"iARMA" <- function(phi = numeric(0), theta = numeric(0), period = 0) {
######## information matrix of arma#####################
if (!(InvertibleQ(phi) & InvertibleQ(theta))) {
cat("Model is non-causal or non-invertible\n")
return(NULL)
}
p <- length(phi)
q <- length(theta)
unames <- character(0)
vnames <- character(0)
if (p > 0) {
if (p > 1) {
vvmatrix <- (tccfAR(phi, phi)[-(1:(p - 1))])[-(p + 1)]
} else if (p == 1) {
vvmatrix <- tccfAR(phi, phi)[-(p + 1)]
}
vvmatrix <- toeplitz(vvmatrix)
imatrix <- vvmatrix
if (!period)
vnames <- paste("phi(", 1:p, ")", sep = "") else vnames <- paste("phi.", period, "(", 1:p, ")", sep = "")
}
if (q > 0) {
if (q > 1) {
uumatrix <- (tccfAR(theta, theta)[-(1:(q - 1))])[-(q + 1)]
} else if (q == 1) {
uumatrix <- tccfAR(theta, theta)[-(q + 1)]
}
uumatrix <- toeplitz(uumatrix)
imatrix <- uumatrix
if (!period)
unames <- paste("theta(", 1:q, ")", sep = "") else unames <- paste("theta.", period, "(", 1:q, ")", sep = "")
}
if (p > 0 && q > 0) {
uvmatrix <- matrix(numeric(1), nrow = p, ncol = q)
tuv <- -tccfAR(phi, theta)
for (i in 1:p) {
for (j in 1:q) {
uvmatrix[i, j] <- tuv[q + i - j]
}
}
imatrix <- cbind(rbind(vvmatrix, t(uvmatrix)), rbind(uvmatrix, uumatrix))
}
inames <- c(vnames, unames)
dimnames(imatrix) <- list(inames, inames)
return(imatrix)
}
"tccfAR" <- function(phi, theta) {
# auxilary function used with iarma######### computes the theoretical cross-covariance
# function of two autoregressions z[t]-phi[1] z_[t-1] --- phi[p] z[t-p] = a[t]
# z[t]-theta[1] z_[t-1] --- theta[q] z[t-q] = a[t] where p, q are length(phi),
# length(theta)
p <- length(phi)
q <- length(theta)
if (p == 0 || q == 0)
return(numeric(0))
k <- p + q
rhs <- c(-1, rep(0, k - 1))
A <- matrix(numeric(k^2), nrow = k, ncol = k)
for (i in 1:k) {
for (j in 1:k) {
imj <- i - j
ijq <- i + j - q - 1
if (i > q) {
if (i > j && imj <= q)
A[i, j] <- theta[imj] else if (i > q && imj == 0)
A[i, j] <- -1
} else {
if (ijq > 0 && ijq <= p)
A[i, j] <- phi[ijq] else if (ijq == 0)
A[i, j] <- -1
}
}
}
return(solve(A, rhs))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.