R/arfima_function.R
In JQVeenstra/arfima: Fractional ARIMA (and Other Long Memory) Time Series Modeling
Defines functions
bestModes
distance
weed
pcircle
Documented in
bestModes
distance
weed
#' Fit ARFIMA, ARIMA-FGN, and ARIMA-PLA (multi-start) models
#'
#' Fits ARFIMA/ARIMA-FGN/ARIMA-PLA multi-start models to times series data.
#' Options include fixing parameters, whether or not to fit fractional noise,
#' what type of fractional noise (fractional Gaussian noise (FGN), fractionally
#' differenced white noise (FDWN), or the newly introduced power-law
#' autocovariance noise (PLA)), etc. This function can fit regressions with
#' ARFIMA/ARIMA-FGN/ARIMA-PLA errors via the xreg argument, including dynamic
#' regression (transfer functions).
#'
#' Fits by direct optimization using optim. The optimizer choices are: 0 -
#' BFGS; 1 - Nealder-Mead; 2 - SANN; otherwise CG.
#'
#' A word of warning: it is generally better to use the default, and only use
#' Nelder-Mead to check for spurious modes. SANN takes a long time (and may
#' only find one mode), and CG may not be stable.
#'
#' If using Nelder-Mead, it must be stressed that Nelder-Mead can take out
#' non-spurious modes or add spurious modes: we have checked visually where we
#' could. Therefore it is wise to use BFGS as the default and if there are
#' modes close to the boundaries, check using Nelder-Mead.
#'
#' The moving average parameters are in the Box-Jenkins convention: they are
#' the negative of the parameters given by \code{\link{arima}}. That is, the
#' model to be fit is, in the case of a non-seasonal ARIMA model, phi(B)
#' (1-B)^d z[t] = theta(B) a[t], where phi(B) = 1 - phi(1) B - ... - phi(p) B^p
#' and theta(B) = 1 - theta(1) B - ... - theta(q) B^q.
#'
#' For the useC parameter, a "0" means no C is used; a "1" means C is only used
#' to compute the log-likelihood, but not the theoretical autocovariance
#' function (tacvf); a "2" means that C is used to compute the tacvf and not
#' the log-likelihood; and a "3" means C is used to compute everything.
#'
#' @param z The data set (time series)
#' @param order The order of the ARIMA model to be fit: c(p, d, q). We have
#' that p is the number of AR parameters (phi), d is the amount of integer
#' differencing, and q is the number of MA parameters (theta). Note we use the
#' Box-Jenkins convention for the MA parameters, in that they are the negative
#' of \code{\link{arima}}: see "Details".
#' @param numeach The number of starts to fit for each parameter. The first
#' argument in the vector is the number of starts for each AR/MA parameter,
#' while the second is the number of starts for the fractional parameter. When
#' this is set to 0, no fractional noise is fit. Note that the number of
#' starts in total is multiplicative: if we are fitting an ARFIMA(2, d, 2), and
#' use the older number of starts (c(2, 2)), we will have 2^2 * 2 * 2^2 = 32
#' starting values for the fits. \strong{Note that the default has changed
#' from c(2, 2) to c(1, 1) since package version 1.4-0}
#' @param back Setting this to true will restore the defaults in numeach.
#' @param dmean Whether the mean should be fit dynamically with the optimizer.
#' Note that the likelihood surface will change if this is TRUE, but this
#' is usually not worrisome.
#' See the referenced thesis for details.
#' @param whichopt Which optimizer to use in the optimization: see "Details".
#' @param itmean This option is under investigation, and will be set to FALSE
#' automatically until it has been decided what to do.
#'
#' Whether the mean should be fit iteratively using the function
#' \code{\link{TrenchMean}}. Currently itmean, if set to TRUE, has higher
#' priority that dmean: if both are TRUE, dmean will be set to FALSE, with a
#' warning.
#' @param fixed A list of parameters to be fixed. If we are to fix certain
#' elements of the AR process, for example, fixed$phi must have length equal to
#' p. Any numeric value will fix the parameter at that value; for example, if
#' we are modelling an AR(2) process, and we wish to fix only the first
#' autoregressive parameter to 0, we would have fixed = list(phi = c(0, NA)).
#' NA corresponds to that parameter being allowed to change in the optimization
#' process. We can fix the fractional parameters, and unlike
#' \code{\link{arima}}, can fix the seasonal parameters as well. Currently,
#' fixing regression/transfer function parameters is disabled.
#' @param lmodel The long memory model (noise type) to be used: "d" for FDWN,
#' "g" for FGN, "h" for PLA, and "n" for none (i.e. ARMA short memory models).
#' Default is "d".
#' @param seasonal The seasonal components of the model we wish to fit, with
#' the same components as above. The period must be supplied.
#' @param useC How much interfaced C code to use: an integer between 0 and 3.
#' The value 3 is strongly recommended. See "Details".
#' @param cpus The number of CPUs used to perform the multi-start fits. A
#' small number of fits and a high number of cpus (say both equal 4) with n not
#' large can actually be slower than when cpus = 1. The number of CPUs should
#' not exceed the number of threads available to R.
#' @param rand Whether random starts are used in the multistart method.
#' Defaults to FALSE.
#' @param numrand The number of random starts to use.
#' @param seed The seed for the random starts.
#' @param eps3 How far to start from the boundaries when using a grid for the
#' multi-starts (i.e. when rand is FALSE.)
#' @param xreg A matrix, data frame, or vector of regressors for regression or
#' transfer functions.
#' @param reglist A list with the following elements:
#' \itemize{
#' \item regpar -
#' either NA or a list, matrix, data frame, or vector with 3 columns. If
#' regpar is a vector, the matrix xreg must have one row or column only. In
#' order, the elements of regpar are: r, s, and b. The values of r are the the
#' orders of the delta parameters as in Box, Jenkins and Reinsel, the values of
#' s are the orders of omega parameters, and the values of b are the
#' backshifting to be done.
#'
#' \item minn - the minimum value for the starting value of the search, if
#' reglist$numeach > 1.
#' \item maxx - the maximum value for the starting value
#' of the search, if reglist$numeach > 1.
#' \item numeach - the number of starts
#' to try for each regression parameter.
#'
#' }
#' @param check If TRUE, checks at each optim iteration whether the model is
#' identifiable. This makes the optimization much slower.
#' @param autoweed Whether to automatically (before the fit is returned) weed
#' out modes found that are found that are close together (usually the same
#' point.)
#' @param weedeps The maximum distance between modes that are close together
#' for the mode with the lower log-likelihood to be weeded out. If adapt is
#' TRUE (default) this value changes.
#' @param adapt If TRUE, if dim is the dimensionality of the search, weedeps is
#' changed to \eqn{(1 + weedeps)^{dim} - 1}.
#' @param weedtype The type of weeding to be done. See \code{\link{weed}}.
#' @param weedp The p in the p-norm to be used in the weeding. p = 2 (default)
#' is Euclidean distance.
#' @param quiet If TRUE, no auxiliary output is generated. The default (FALSE)
#' has information of fits being proformed.
#' @param startfit Meant primarily for debugging (for now), allows starting places
#' for the fitting process. Overrides \code{numeach}.
#' @return An object of class "arfima". In it, full information on the fit is
#' given, though not printed under the print.arfima method. The phis are the
#' AR parameters, and the thetas are the MA parameters. Residuals, regression
#' residuals, etc., are all available, along with the parameter values and
#' standard errors. Note that the muHat returned in the arfima object is
#' of the \strong{differenced} series, if differencing is applied.
#'
#' Note that if multiple modes are found, they are listed in order of
#' log-likelihood value.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima.sim}}, \code{\link{SeriesJ}},
#' \code{\link{arfima-package}}
#' @references McLeod, A. I., Yu, H. and Krougly, Z. L. (2007) Algorithms for
#' Linear Time Series Analysis: With R Package Journal of Statistical Software,
#' Vol. 23, Issue 5
#'
#' 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
#'
#'\donttest{
#' set.seed(8564)
#' sim <- arfima.sim(1000, model = list(phi = c(0.2, 0.1),
#' dfrac = 0.4, theta = 0.9))
#' fit <- arfima(sim, order = c(2, 0, 1), back=TRUE)
#'
#' fit
#'
#' data(tmpyr)
#'
#' fit <- arfima(tmpyr, order = c(1, 0, 1), numeach = c(3, 3))
#' fit
#'
#' plot(tacvf(fit), maxlag = 30, tacf = TRUE)
#'
#' data(SeriesJ)
#' attach(SeriesJ)
#'
#' fitTF <- arfima(YJ, order= c(2, 0, 0), xreg = XJ, reglist =
#' list(regpar = c(2, 2, 3)), lmodel = "n")
#' fitTF
#'
#' detach(SeriesJ)
#'}
#'
#' @export arfima
"arfima" <- function(z, order = c(0, 0, 0), numeach = c(1, 1), dmean = TRUE, whichopt = 0,
itmean = FALSE, fixed = list(phi = NA, theta = NA, frac = NA, seasonal = list(phi = NA,
theta = NA, frac = NA), reg = NA), lmodel = c("d", "g", "h", "n"), seasonal = list(order = c(0,
0, 0), period = NA, lmodel = c("d", "g", "h", "n"), numeach = c(1, 1)), useC = 3,
cpus = 1, rand = FALSE, numrand = NULL, seed = NA, eps3 = 0.01, xreg = NULL, reglist = list(regpar = NA,
minn = -10, maxx = 10, numeach = 1), check = F, autoweed = TRUE, weedeps = 0.01,
adapt = TRUE, weedtype = c("A", "P", "B"), weedp = 2, quiet = FALSE, startfit = NULL, back = FALSE) {
if (itmean) {
warning("itmean is currently under examination, and is not useable at the moment.\nSetting itmean to FALSE")
itmean <- FALSE
}
if(back) {
numeach <- c(2, 2)
seasonal$numeach <- c(2, 2)
dmean <- TRUE
}
if(!is.null(startfit)) {
stopifnot(is.matrix(startfit))
starts <- startfit
num <- nrow(starts)
startflag <- TRUE
}
else startflag <- FALSE
if (any(is.na(z))) {
stop("Missing data not supported at the moment.")
}
cpus <- abs(round(cpus))
if (cpus < 1)
cpus <- 1
getHess <- TRUE ## to save on options but not change code
inclbdries <- FALSE ## to save on options but not change code
tse <- FALSE ## to save on options but not change code
regpar <- reglist$regpar
straightRegress <- FALSE
namexreg <- NULL
intindex <- 0
intname <- NULL
if (!is.null(regpar) && all(is.na(regpar)))
regpar <- NULL
if (!is.null(xreg)) {
if (!(is.data.frame(xreg) || is.matrix(xreg) || is.vector(xreg) || is.ts(xreg)))
stop("Invalid xreg type")
if (any(is.na(xreg)))
stop("no missing values allowed in xreg: if dynamic regression needs padding, please pad with zeroes")
constantInds <- numeric(0)
if (is.null(regpar)) {
if (!is.data.frame(xreg) && !is.matrix(xreg)) {
xreg <- data.frame(xreg)
colnames(xreg) <- as.character(as.list((match.call()))$xreg)
}
namexreg <- colnames(xreg)
if (is.null(namexreg))
namexreg <- paste(rep("X", ncol(xreg)), 1:ncol(xreg), sep = "")
regpar <- matrix(0, nrow = 3, ncol = dim(xreg)[2])
regpar <- as.data.frame(t(regpar))
names(regpar) <- c("r", "s", "b")
constantInds <- c(1:dim(xreg)[2])[apply(xreg, 2, var) == 0]
straightRegress <- ncol(xreg)
if(length(constantInds)) {
for(i in constantInds) {
if(all(xreg[, i] == 1)) {
constantInds <- setdiff(constantInds, i)
if(i > 1) {
xreg <- xreg[, c(i, setdiff(1:ncol(xreg), i))]
namexreg <- namexreg[c(i, setdiff(1:ncol(xreg), i))]
intindex <- 1i
if(length(constantInds)>0)
constantInds <- c(1:dim(xreg)[2])[apply(xreg, 2, var) == 0][-1]
}
intindex <- 1
break
}
}
straightRegress <- straightRegress - length(constantInds) - 1
}
straightRegress <- (straightRegress>0)
if(!straightRegress) {
warning("All regression terms constant. Setting xreg to NULL")
xreg <- NULL
}
else if(straightRegress && dmean != TRUE) {
warning("Regression with time series errors requires dmean=T, setting")
dmean <- TRUE
}
else if(straightRegress && intindex && is.numeric(dmean)) {
warning("with an intercept term, dmean must be TRUE. Setting it so.")
dmean <- TRUE
}
}
else {
if (is.matrix(regpar) || is.data.frame(regpar)) {
if (dim(regpar)[2] != dim(xreg)[2]) {
stop("xreg and regpar sizes do not match")
}
namexreg <- colnames(xreg)
if (is.null(namexreg))
namexreg <- paste(rep("X", ncol(xreg)), 1:ncol(xreg), sep = "")
if (any(is.na(regpar)))
stop("no missing values allowed in regpar")
if (dim(regpar)[1] != 3)
stop("regpar does not have 3 rows")
if (!is.null(names(regpar)) && names(regpar) != c("r", "s", "b"))
stop("names(regpar), if not NULL, must be r, s, and b")
regpar <- as.data.frame(t(regpar))
names(regpar) <- c("r", "s", "b")
} else if (is.vector(xreg) || is.ts(xreg)) {
pars <- data.frame(r = regpar[1], s = regpar[2], b = regpar[3])
regpar <- pars
xreg <- data.frame(xreg)
colnames(xreg) <- as.character(as.list((match.call()))$xreg)
namexreg <- colnames(xreg)
if(is.null(namexreg))
namexreg <- "X1"
} else stop("invalid regpar type")
}
xreg <- as.matrix(xreg)
if (length(constantInds) > 0) {
warning(paste("columns", paste(constantInds, collapse = ", "), "are constant and will be removed"))
xreg <- xreg[, -constantInds]
namexreg <- namexreg[-constantInds]
}
if(intindex && is.logical(dmean) && dmean) {
regpar <- regpar[-1, ]
xreg <- xreg[, -1]
intname <- namexreg[1]
namexreg <- namexreg[-1]
}
else if(intindex) {
warning(paste(paste("there is a constant column named", namexreg[1]),
"take to be the intercept, but dmean is not TRUE. This column",
"will be removed and ignored.", sep = "\n"))
regpar <- regpar[-1, ]
xreg <- xreg[, -1]
intindex <- 0
}
else if(is.logical(dmean) && dmean) {
intname <- "Intercept"
}
s <- regpar$s <- regpar$s + 1
r <- regpar$r
b <- regpar$b
if (any(r > 0) || any(s>1)) {
#cat("note: transfer functions do not work with dynamic mean: setting dmean to FALSE\n")
dmean <- FALSE
}
numvarreg <- sum(s) + sum(r)
regeach <- if(straightRegress) 1 else reglist$numeach
if (length(regeach) == 0)
regeach <- 1
if (length(regeach) != 1 || regeach != round(regeach) || regeach <= 0) {
warning("invalid reglist$numeach: resetting to 1")
regeach <- 1
}
regmin <- reglist$minn
regmax <- reglist$maxx
if (length(regmin) == 0 || is.na(regmin))
regmin <- if(regeach>1) -10 else 0
if (length(regmax) == 0 || is.na(regmax))
regmax <- if(regeach>1) 10 else 0
if (itmean)
stop("Please note that for regression problems, the mean is calculated separately\n via optimization or on the regression parameters or by the mean of the series;\n therefore itmean is not valid")
}
else {
r <- s <- b <- 0
numvarreg <- regeach <- 0
namexreg <- reglist <- regpar <- NULL
}
if (is.numeric(dmean) && !is.double(dmean))
stop("dmean must be logical or a double")
if (itmean && ((is.logical(dmean) && dmean) || is.numeric(dmean))) {
warning("itmean trumps dmean: setting dmean to FALSE")
dmean <- FALSE
}
mean.dof = if(itmean || (is.logical(dmean) && dmean)) 1 else 0
if(!is.null(xreg)) mean.dof <- mean.dof + dim(xreg)[2]
info <- NULL
dint <- order[2]
if (length(dint) == 0)
dint <- 0
lmodel <- lmodel[1]
if (is.null(numeach) || (length(numeach) != 2 && (lmodel != "n" && length(numeach) !=
1)) || !is.vector(numeach))
stop("numeach must be a vector of length 1 or 2:\n \n\t\tposition 1 is the values for ARMA\n parameters, position 2 is the value for the fractional parameter\n\n\t\tif a vector of length 1, we must have lmodel = \"n\"")
if (numeach[1] != round(numeach[1]) || numeach[1] <= 0) {
warning("incorrect numeach ARMA: resetting to 1")
numeach[1] <- 1
}
if (numeach[2] != round(numeach[2]) || numeach[2] < 0) {
warning("incorrect numeach fractional: resetting to 1")
numeach[2] <- 1
}
if (numeach[2] == 0 && lmodel != "n") {
cat("lmodel = \"", lmodel, "\" and numeach[2] = 0; setting lmodel = \"n\"\n", sep = "")
lmodel <- "n"
}
if (lmodel == "n") {
if (length(numeach) == 1)
numeach <- c(numeach, 0) else numeach[2] <- 0
}
if (!is.null(seasonal) && !is.na(seasonal$period) && !is.null(seasonal$period) && (seasonal$period !=
0) && ((!is.null(seasonal$order) || !all(seasonal$order == 0)) || (!is.null(seasonal$lmodel) &&
seasonal$lmodel != "n"))) {
dseas <- seasonal$order[2]
if (length(dseas) == 0)
dseas <- 0
pseas <- seasonal$order[1]
qseas <- seasonal$order[3]
if (length(pseas) == 0)
pseas <- 0
if (length(qseas) == 0)
qseas <- 0
if (!is.null(seasonal$lmodel))
slmodel <- seasonal$lmodel[1] else slmodel <- "d"
if (!is.null(seasonal$numeach))
snumeach <- seasonal$numeach else {
warning("no seasonal$numeach: setting to numeach")
snumeach <- numeach
}
if (is.null(snumeach) || (length(snumeach) != 2 && (slmodel != "n" && length(snumeach) !=
1)) || !is.vector(snumeach))
stop("seasonal$numeach ust be a vector of length 1 or 2:\n \n\t\tposition 1 is the values for ARMA\n parameters, position 2 is the value for the fractional parameter\n\n\t\tif a vector of length 1, we must have seasonal$lmodel = \"n\"")
period <- seasonal$period
if (snumeach[1] != round(snumeach[1]) || snumeach[1] <= 0) {
warning("incorrect seasonal$numeach ARMA: resetting to 1")
snumeach[1] <- 1
}
if (snumeach[2] != round(snumeach[2]) || snumeach[2] < 0) {
warning("incorrect seasonal$numeach fractional: resetting to 1")
snumeach[2] <- 1
}
if (snumeach[2] == 0 && slmodel != "n") {
cat("seasonal lmodel =", slmodel, "and seasonal numeach[2] = 0; setting seasonal lmodel = \"n\"\n")
slmodel <- "n"
}
if (slmodel == "n") {
if (length(snumeach) == 1)
snumeach <- c(snumeach, 0) else snumeach[2] <- 0
}
slmr <- if (slmodel != "n")
1 else 0
} else {
dseas <- pseas <- qseas <- period <- slmr <- 0
snumeach <- c(0, 0)
slmodel <- "n"
}
p <- order[1]
q <- order[3]
lmr <- if (lmodel != "n")
1 else 0
numvars <- p + q + pseas + qseas + lmr + slmr + numvarreg +
if (dmean == TRUE)
1 else 0
if (length(fixed) > 0) {
if (!(is.null(fixed$reg)) && !(all(is.na(fixed$reg))) && length(fixed$reg) > 0) {
warning("no fixed regression parameters for now.")
fixed$reg <- NA
}
arfixed <- fixed$phi
mafixed <- fixed$theta
fracfix <- fixed$frac
if((length(arfixed)==1 && is.na(arfixed))||length(arfixed)==0)
arfixed <- rep(NA, p)
if((length(mafixed)==1 && is.na(mafixed))||length(mafixed)==0)
mafixed <- rep(NA, q)
if(length(fracfix)==0 && lmodel !="n") {
fracfix <- NA
}
warnstring <- "In the fixed argument:"
if (length(arfixed) != p) {
warnstring <- paste0(warnstring, "\n length phi is ", length(arfixed),
" but p = ", p, " unfixing phi")
arfixed <- if (p > 0)
rep(NA, p) else numeric(0)
}
else if(p>0 && length(arfixed)==p&&any(!is.na(arfixed))) {
if(sum(is.na(arfixed))==0) {
if(!InvertibleQ(arfixed)) {
stop("phi parameters fixed at nonstationary parameters")
}
}
}
if (length(mafixed) != q) {
warnstring <- paste0(warnstring, "\n length theta is ", length(mafixed),
" but q = ", q, " unfixing theta")
mafixed <- if (q > 0)
rep(NA, q) else numeric(0)
}
else if(any(!is.na(mafixed))) {
if(sum(is.na(mafixed))==0) {
if(!InvertibleQ(mafixed)) {
stop("theta parameters fixed at noninvertible parameters")
}
}
}
if (!is.null(fracfix) && length(fracfix) != 1 && !is.na(fracfix)) {
warnstring <- paste0(warnstring, "\n fractional parameter is fixed greater",
"\n than length 1, resetting fractional prameter ")
if(lmodel=="n") {
warnstring<- paste0(warnstring,"\n lmodel == 'n' and fractional parameter fixed, resetting")
fracfix <- numeric(0)
}
else fracfix <- NA
}
if (lmodel == "n")
fracfix <- numeric(0)
if(lmodel!="n"&&(!is.null(fracfix)&&!is.na(fracfix))) {
fix_stop <- F
if(lmodel=="d") {
if(!InvertibleD(fracfix)) fix_stop <- T
}
else if(lmodel=="g") {
if(!InvertibleH(fracfix)) fix_stop <- T
}
else {
if(!InvertibleAlpha(fracfix)) fix_stop <- T
}
if(fix_stop){
stop("fractional parameter fixed at nonstationary/noninvertible value")
}
}
fixx <- c(arfixed, mafixed)
if (period > 0) {
sarfixed <- fixed$seasonal$phi
smafixed <- fixed$seasonal$theta
sfracfix <- fixed$seasonal$frac
if((length(sarfixed)==1 && is.na(sarfixed))||length(sarfixed)==0)
sarfixed <- rep(NA, pseas)
if((length(smafixed)==1 && is.na(smafixed))||length(smafixed)==0)
smafixed <- rep(NA, qseas)
if(length(sfracfix)==0 && slmodel !="n") {
sfracfix <- NA
}
if (length(sarfixed) != pseas) {
warnstring <- paste0(warnstring, "\n length seasonal phi is ", length(sarfixed),
" but pseas = ", pseas, " unfixing s_phi")
sarfixed <- if (pseas > 0)
rep(NA, pseas) else numeric(0)
}
else if(pseas>0 && length(sarfixed)==pseas&&any(!is.na(sarfixed))) {
if(sum(is.na(sarfixed))==0) {
if(!InvertibleQ(sarfixed)) {
stop("seasonal phi parameters fixed at nonstationary parameters")
}
}
}
if (length(smafixed) != qseas) {
warnstring <- paste0(warnstring, "\n length seasonal theta is ", length(smafixed),
" but qseas = ", qseas, " unfixing s_theta")
smafixed <- if (qseas > 0)
rep(NA, qseas) else numeric(0)
}
else if(qseas>0 && length(smafixed)==qseas&&any(!is.na(smafixed))) {
if(sum(is.na(smafixed))==0) {
if(!InvertibleQ(smafixed)) {
stop("seasonal theta parameters fixed at noninvertible parameters")
}
}
}
if (!is.null(sfracfix) && length(sfracfix) != 1 && !is.na(sfracfix)) {
warnstring <- paste0(warnstring, "\n seasonal fractional parameter is",
"\n fixed greater than length 1, resetting seasonal",
"\n fractional prameter ")
if(slmodel=="n") {
warnstring("\n slmodel == 'n' and seasonal fractional parameter fixed",
"\n resetting")
sfracfix <- numeric(0)
}
else sfracfix <- NA
}
if(slmodel!="n"&&(!is.null(sfracfix)&&!is.na(sfracfix))) {
fix_stop <- F
if(slmodel=="d") {
if(!InvertibleD(sfracfix)) fix_stop <- T
}
else if(slmodel=="g") {
if(!InvertibleH(sfracfix)) fix_stop <- T
}
else {
if(!InvertibleAlpha(sfracfix)) fix_stop <- T
}
if(fix_stop) {
stop("seasonal fractional parameter fixed at nonstationary/noninvertible value")
}
}
sfixx <- c(sarfixed, smafixed)
} else {
sfixx <- NULL
sfracfix <- NULL
}
if(nchar(warnstring)>31) warning(warnstring)
fixx <- c(fixx, sfixx, fracfix, sfracfix, rep(NA, sum(r)+sum(s)))
indfixx <- !is.na(fixx)
if (is.logical(dmean) && dmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
} else if (is.logical(itmean) && itmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
}
} else {
fixx <- rep(NA, numvars)
indfixx <- rep(FALSE, numvars)
if (is.logical(dmean) && dmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
} else if (is.logical(itmean) && itmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
}
}
numtries <- 1
if (p + q > 0)
numtries <- numtries * numeach[1]^(p + q)
if (numeach[2] > 1)
numtries <- numtries * numeach[2]
if (pseas + qseas > 0)
numtries <- numtries * snumeach[1]^(pseas + qseas)
if (snumeach[2] > 1)
numtries <- numtries * snumeach[2]
if (regeach > 1)
numtries <- numtries * regeach
if (sum(as.integer(indfixx)) > 0) {
if (p > 0 && any(!is.na(arfixed)))
numtries <- numtries/(numeach[1]^(sum(!is.na(arfixed))))
if (q > 0 && any(!is.na(mafixed)))
numtries <- numtries/(numeach[1]^(sum(!is.na(mafixed))))
if (lmodel != "n" && (is.null(fracfix)||!is.na(fracfix)))
numtries <- numtries/numeach[2]
if (pseas > 0 && any(!is.na(sarfixed)))
numtries <- numtries/(snumeach[1]^(sum(!is.na(sarfixed))))
if (qseas > 0 && any(!is.na(smafixed)))
numtries <- numtries/(snumeach[1]^(sum(!is.na(smafixed))))
if (slmodel != "n" && (is.null(sfracfix)||!is.na(sfracfix)))
numtries <- numtries/snumeach[2]
}
# if(numvars == 0 && numtries > 1) stop('invalid specification of variables')
if (round(dint) != dint || dint < 0)
stop("d must be an integer >= 0")
if (round(dseas) != dseas || dseas < 0)
stop("dseas must be an integer >= 0")
if (!quiet) {
if(!startflag) {
if (autoweed && numtries > 1 && !startflag)
cat("Note: autoweed is ON. It is possible, but not likely,\n that unique modes may be lost.\n") else if (!autoweed && numtries > 1)
cat("Note: autoweed is OFF. Running weed after the\n fit is highly recommended.\n") else {
cat("Note: only one starting point. Only one mode can be found -- this is now the default behavior.\n")
}
}
else {
if(autoweed)
cat("Note: autoweed is ON. It is possible, but not likely,\n that unique modes may be lost.\n")
else
cat("Note: autoweed is OFF. Running weed after the\n fit may be required.\n")
cat("Starting fit with given start matrix.\n")
}
}
y <- z
xr <- xreg
if (dint > 0) {
y <- diff(y, differences = dint)
if (!is.null(xreg))
xreg <- diff(xreg, differences = dint)
}
if (dseas > 0) {
y <- diff(y, differences = dseas, lag = period)
if (!is.null(xreg))
xreg <- diff(xreg, differences = dseas, lag = period)
}
D <- dint + dseas * period
if (any(r > 0) && any(abs(colMeans(xreg)) > 1e-14)) {
cat("Please note for transfer functions the means of each X variable must be 0: subtracting mean from each X\n")
for (i in 1:ncol(xreg)) {
xreg[, i] <- xreg[, i] - colMeans(xreg)[i]
}
}
flag0 <- regOnly <- FALSE
if (p + q + numeach[2] + pseas + qseas + snumeach[2] == 0) {
if(is.null(xreg))
flag0 <- TRUE
else
regOnly <- TRUE
}
differencing <- if (D > 0) TRUE else FALSE
if (is.logical(dmean) && !dmean && !straightRegress)
meanval <- mean(y)
else if (is.numeric(dmean))
meanval <- dmean
else
meanval <- 0
sumfixed <- sum(indfixx)
numvars1 <- numvars - sumfixed
if (numvars1 == 0)
getHess <- FALSE
if (flag0) {
if(sum(numeach)>2||sum(seasonal$numeach)>2||
(!is.null(numrand)) && numrand>0)
warning('With only the mean to be fit, since the problem is convex,
it has a unique minimum and thus only one fit will be performed.')
if ((is.logical(dmean) && dmean) || itmean) {
fun <- function(m) {
lARFIMA(y - m)
}
fit <- optim(mean(y), fn = fun, method = "Brent", lower = mean(y) - 10, upper = mean(y) +
10, control = list(trace = 0, fnscale = -1), hessian = getHess)
hess <- fit$hessian[1, 1]
se <- tryCatch(sqrt(diag(solve(-hess))), error = function(err) rep(NA, dim(hess)[1]))
res <- list(maximum = fit[[1]][1], objective = fit[[2]][1], hess = hess, se = se)
} else if (is.logical(dmean)) {
res <- list(maximum = mean(y), objective = lARFIMA(y - mean(y)), hess = NA,
se = -Inf)
} else {
res <- list(maximum = dmean, objective = lARFIMA(y - dmean), hess = NA, se = -Inf)
}
allpars <- vector("list", 1)
allpars[[1]]$muHat <- res$maximum
allpars[[1]]$loglik <- res$objective
allpars[[1]]$hessian <- res$hess
allpars[[1]]$se <- res$se
allpars[[1]]$phi <- allpars[[1]]$theta <- allpars[[1]]$phiseas <-
allpars[[1]]$thetaseas <- allpars[[1]]$dfrac <- allpars[[1]]$dfs <-
allpars[[1]]$phip <- allpars[[1]]$thetap <- allpars[[1]]$H <-
allpars[[1]]$alpha <- allpars[[1]]$alphas <- allpars[[1]]$delta <-
allpars[[1]]$omega <- allpars[[1]]$Hs <- numeric(0)
allpars[[1]]$pars <-
c(allpars[[1]]$phi, allpars[[1]]$theta, allpars[[1]]$phiseas,
allpars[[1]]$thetaseas, allpars[[1]]$dfrac, allpars[[1]]$H,
allpars[[1]]$alpha, allpars[[1]]$dfs, allpars[[1]]$Hs,
allpars[[1]]$alphas, allpars[[1]]$delta, allpars[[1]]$omega)
sigma2muse <- getsigma2muhat(y - allpars[[1]]$muHat, mean.dof = mean.dof)
allpars[[1]]$sigma2 <- sigma2muse
rr <- tacvfARFIMA(maxlag = length(y) - 1)
res <- DLResiduals(rr, y - allpars[[1]]$muHat)
allpars[[1]]$residuals <- res
allpars[[1]]$fitted <- y - allpars[[1]]$muHat - res
class(allpars[[1]]) <- "ARFIMA"
ans <- list(z = z, differencing = differencing, dint = dint, dseas = dseas, period = period,
modes = allpars, dmean = dmean, itmean = itmean, p = p, q = q, pseas = pseas,
qseas = qseas, lmodel = lmodel, slmodel = slmodel, weeded = TRUE, getHess = getHess,
numvars = numvars, numcut = 0, n = length(z), xreg = xr, r = r, s = s, b = b,
call = match.call(), flag0 = flag0, numeach = numeach, strReg = F,
regOnly = regOnly, intindex = intindex, intname = intname, namexreg=namexreg,
fixx = fixx, indfixx = indfixx)
} else {
num <- if (!rand)
numtries else numrand
if(startflag) num <- nrow(starts)
if (rand && !startflag) {
if (length(num) == 0) {
warning("random starts selected and no number of starts selected: setting to 8")
num <- 8
}
starts <- matrix(0, nrow = num, ncol = numvars)
pl1 <- 0
if (!is.na(seed) && !is.null(seed))
set.seed(seed) else warning("Random start seed not selected: results may not be reproducible")
if (p + q + pseas + qseas > 0) {
starts[, 1:(p + q + pseas + qseas)] <- matrix(runif(num * (p + q + pseas +
qseas), min = -1, max = 1), ncol = (p + q + pseas + qseas))
}
if (lmodel != "n") {
pl1 <- 1
if (lmodel == "d")
starts[, (1 + p + q + pseas + qseas)] <- runif(1, min = -1, max = 0.5)
else if (lmodel == "g")
starts[, (1 + p + q + pseas + qseas)] <- runif(1, min = 0, max = 1)
else
starts[, (1 + p + q + pseas + qseas)] <- runif(1, min = 0, max = 3)
}
if (slmodel != "n") {
if (slmodel == "d")
starts[, (1 + p + q + pl1 + pseas + qseas)] <- runif(1, min = -1, max = 0.5)
else if (slmodel == "g")
starts[, (1 + p + q + pl1 + pseas + qseas)] <- runif(1, min = 0, max = 1)
else
starts[, (1 + p + q + pl1 + pseas + qseas)] <- runif(1, min = 0, max = 3)
}
if (!is.null(xreg)) {
if(!straightRegress) {
starts[, (1 + p + q + pl1 + pseas + qseas):(1 + p + q + pl1 + pseas + qseas +
numvarreg)] <- matrix(runif(num * (numvarreg), min = regmin,
max = regmax), nrow = num)
}
else if(is.logical(dmean)&&dmean) {
coeffs <- lsfit(xreg, y)$coef
for(ii in 1:numvarreg) starts[, (1 + p + q + pl1 + pseas + qseas+ii)] <- coeffs
}
else {
coeffs <- lsfit(xreg, y, intercept = FALSE)$coef
for(ii in 1:numvarreg) starts[, (1 + p + q + pl1 + pseas + qseas+ii)] <- coeffs
}
}
} else if (num > 1 && !startflag) {
pl1 <- 0
li <- vector("list", numvars)
if (p + q > 0) {
if (inclbdries)
seqq <- seq(-1 + eps3, 1 - eps3, length.out = numeach[1]) else {
seqq <- seq(-1 + eps3, 1 - eps3, length.out = numeach[1] + 2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
}
if (pseas + qseas > 0) {
if (inclbdries)
seqq <- seq(-1 + eps3, 1 - eps3, length.out = snumeach[1]) else {
seqq <- seq(-1 + eps3, 1 - eps3, length.out = snumeach[1] + 2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
}
for (i in 1:(p + q + pseas + qseas))
li[[i]] <- seqq
if (lmodel != "n") {
pl1 <- 1
if (lmodel == "d") {
minner <- -1
maxxer <- 0.5
} else if (lmodel == "g") {
minner <- 0
maxxer <- 1
} else {
minner <- 0
maxxer <- 3
}
if (inclbdries)
seqq <- seq(minner + eps3, maxxer - eps3, length.out = numeach[2]) else {
seqq <- seq(minner + eps3, maxxer - eps3, length.out = numeach[2] +
2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
li[[(p + q + pl1 + pseas + qseas)]] <- seqq
}
pl1_1 <- 0
if (slmodel != "n") {
pl1_1 <- 1
if (slmodel == "d") {
minner <- -1
maxxer <- 0.5
} else if (slmodel == "g") {
minner <- 0
maxxer <- 1
} else {
minner <- 0
maxxer <- 3
}
if (inclbdries)
seqq <- seq(minner + eps3, maxxer - eps3, length.out = snumeach[2]) else {
seqq <- seq(minner + eps3, maxxer - eps3, length.out = snumeach[2] +
2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
li[[(p + q + 1 + pseas + qseas + pl1)]] <- seqq
}
if (!is.null(xreg)) {
if(!straightRegress) {
if (inclbdries)
seqq <- seq(regmin + eps3, regmax - eps3, length.out = regeach)
else {
seqq <- seq(regmin + eps3, regmax - eps3, length.out = regeach + 2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
for (i in 1:(numvarreg)) li[[(p + q + pseas + qseas + pl1 + pl1_1) + i]] <- seqq
}
else if(is.logical(dmean)&&dmean){
coeffs <- lsfit(xreg, y)$coef[c(2:(ncol(xreg)+1), 1)]
for (i in 1:(numvarreg+1)) li[[(p + q + pseas + qseas + pl1 + pl1_1) + i]] <- coeffs[i]
}
else {
coeffs <- lsfit(xreg, y, intercept = FALSE)$coef
for (i in 1:(numvarreg)) li[[(p + q + pseas + qseas + pl1 + pl1_1) + i]] <- coeffs[i]
}
}
else if(is.logical(dmean) && dmean)
li[[length(li)]] <- mean(y)
starts <- expand.grid(li, KEEP.OUT.ATTRS = FALSE)
if(sumfixed>0) {
for(i in 1:numvars) {
if(indfixx[i]) starts[, i] <- rep(fixx[i], nrow(starts))
}
starts <- unique(starts)
}
} else if(!startflag) {
starts <- if (numvars > 0)
matrix(0, nrow=1,ncol = numvars) else matrix(0)
if(straightRegress && dmean == TRUE) {
coeffs <- lsfit(xreg, y)$coef[c(2:(ncol(xreg)+1), 1)]
starts[, (numvars-length(coeffs)+1):numvars] <- coeffs
}
else if(straightRegress) {
coeffs <- lsfit(xreg, y, intercept = FALSE)$coef
starts[, (numvars-length(coeffs)+1):numvars] <- coeffs
}
else if(dmean == TRUE)
starts[, numvars] <- mean(y)
if(sumfixed>0) {
for(i in 1:numvars) {
if(indfixx[i]) starts[, i] <- fixx[i]
}
}
}
if(dmean==F||is.numeric(dmean)) {
indfixxM <- c(indfixx, TRUE)
numvarsM <- numvars + 1
}
else {
indfixxM <- indfixx
numvarsM <- numvars
}
if (nrow(starts) != num)
stop("error in starts")
if (rand)
strg <- " the random start " else strg <- " the "
if (!quiet) {
cat("Beginning", strg, "fits with ", num, " starting values.\n", sep = "")
cat("\n")
}
fit <- vector("list", num)
if(!straightRegress)
xreg <- as.vector(xreg)
if (cpus == 1) {
for (i in 1:num) {
fit[[i]] <- arfimaFit(y = y, p = p, q = q, pseas = pseas, qseas = qseas,
period = period, lmodel = lmodel, slmodel = slmodel, whichopt = whichopt,
dmean = dmean, parstart = as.double(starts[i, ]), getHess = getHess, itmean = itmean,
indfixx = indfixx, fixx = fixx, xreg = xreg, r = r, s = s, b = b, straightRegress = straightRegress)
}
} else {
cl <- makePSOCKcluster(cpus, methods = FALSE)
numms <- 1:num
clusterExport(cl, c("lARFIMA", "PacfToAR", "DLLoglikelihood", "lARFIMAwTF"),
envir = environment(arfima))
fit <- clusterApplyLB(cl, numms, arfimaFitpar, parstart = starts, y = y, p = p,
q = q, pseas = pseas, qseas = qseas, period = period, lmodel = lmodel, slmodel = slmodel,
whichopt = whichopt, useC = useC, dmean = dmean, getHess = getHess, itmean = itmean,
indfixx = indfixx, fixx = fixx, xreg = xreg, r = r, s = s, b = b, straightRegress = straightRegress)
stopCluster(cl)
setDefaultCluster(NULL)
}
allpars <- vector("list", num)
logl <- rep(0, num)
for (i in 1:num) {
pars <- fit[[i]][[1]]
pars[indfixx] <- fixx[indfixx]
allpars[[i]]$phi <- if (p > 0)
pars[1:p] else numeric(0)
allpars[[i]]$theta <- if (q > 0)
pars[p + (1:q)] else numeric(0)
allpars[[i]]$phiseas <- if (pseas > 0)
pars[p + q + (1:pseas)] else numeric(0)
allpars[[i]]$thetaseas <- if (qseas > 0)
pars[p + q + pseas + (1:qseas)] else numeric(0)
allpars[[i]]$phip <- if (p > 0)
ARToPacf(pars[1:p]) else numeric(0)
allpars[[i]]$thetap <- if (q > 0)
ARToPacf(pars[p + (1:q)]) else numeric(0)
allpars[[i]]$phiseasp <- if (pseas > 0)
ARToPacf(pars[p + q + (1:pseas)]) else numeric(0)
allpars[[i]]$thetaseasp <- if (qseas > 0)
ARToPacf(pars[p + q + pseas + (1:qseas)]) else numeric(0)
allpars[[i]]$dfrac <- if (lmodel == "d")
pars[p + q + pseas + qseas + 1] else numeric(0)
allpars[[i]]$H <- if (lmodel == "g")
pars[p + q + pseas + qseas + 1] else numeric(0)
allpars[[i]]$alpha <- if (lmodel == "h")
pars[p + q + pseas + qseas + 1] else numeric(0)
allpars[[i]]$dfs <- if (slmodel == "d")
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$alpha) + 1] else numeric(0)
allpars[[i]]$Hs <- if (slmodel == "g")
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$alpha) + 1] else numeric(0)
allpars[[i]]$alphas <- if (slmodel == "h")
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$alpha) + 1] else numeric(0)
allpars[[i]]$delta <- if (!is.null(xreg) && sum(r) > 0)
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + 1:(sum(r))] else numeric(0)
allpars[[i]]$omega <- if (!is.null(xreg) && sum(s) > 0)
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + sum(r) + 1:sum(s)] else numeric(0)
allpars[[i]]$muHat <- if (is.logical(dmean) && dmean)
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + sum(r) + sum(s) + 1] else meanval
rr <- tacvfARFIMA(phi = allpars[[i]]$phi, theta = allpars[[i]]$theta, phiseas = allpars[[i]]$phiseas,
thetaseas = allpars[[i]]$thetaseas, dfrac = allpars[[i]]$dfrac, dfs = allpars[[i]]$dfs,
H = allpars[[i]]$H, Hs = allpars[[i]]$Hs, alpha = allpars[[i]]$alpha, alphas = allpars[[i]]$alphas,
period = period, maxlag = length(y) - 1)
if (itmean) {
allpars[[i]]$muHat <- TrenchMean(rr, y)
}
if (is.null(fit[[i]][2])) {
cat("mode ", i, " has no loglikelihood value", sep = "")
allpars[[i]]$loglik <- logl[i] <- -3e+10
} else {
allpars[[i]]$loglik <- logl[i] <- fit[[i]][[2]]
}
allpars[[i]]$pars <- c(allpars[[i]]$phi, allpars[[i]]$theta, allpars[[i]]$phiseas,
allpars[[i]]$thetaseas, allpars[[i]]$dfrac, allpars[[i]]$H, allpars[[i]]$alpha,
allpars[[i]]$dfs, allpars[[i]]$Hs, allpars[[i]]$alphas, allpars[[i]]$delta,
allpars[[i]]$omega)
allpars[[i]]$parsp <- c(allpars[[i]]$phip, allpars[[i]]$thetap, allpars[[i]]$phiseasp,
allpars[[i]]$thetaseasp, allpars[[i]]$dfrac, allpars[[i]]$H, allpars[[i]]$alpha,
allpars[[i]]$dfs, allpars[[i]]$Hs, allpars[[i]]$alphas, allpars[[i]]$delta,
allpars[[i]]$omega)
if (getHess) {
hess <- fit[[i]]$hessian[!indfixx, !indfixx]
se <- rep(-Inf, numvarsM)
se[!indfixxM] <- tryCatch(sqrt(diag(solve(-hess))),
error = function(err) rep(NA, max(dim(hess)[1], 1)))
allpars[[i]]$se <- se
allpars[[i]]$hess <- hess
} else hess <- numeric(0)
if (!is.null(xreg)) {
omega <- allpars[[i]]$omega
delta <- allpars[[i]]$delta
xreglist <- list(s = s, r = r, b = b, omega = vector("list", length(r)),
delta = vector("list", length(r)))
if (!straightRegress) {
for (j in 1:length(r)) {
if (r[j] == 0) {
xreglist$delta[[j]] <- NA
if (s[j] == 1) {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- namexreg[j]
} else {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- paste("omega(", 0:(s[j] - 1), ").",
namexreg[j], sep = "")
}
} else {
if (j == 1)
xreglist$delta[[j]] <- delta[1:r[1]] else {
rrr <- cumsum(r)[j - 1]
xreglist$delta[[j]] <- delta[(rrr + 1):(rrr + r[j])]
}
names(xreglist$delta[[j]]) <- paste("delta(", 1:r[j], ").", namexreg[j],
sep = "")
if (s[j] == 1) {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- paste("omega(0).", namexreg[j], sep = "")
} else {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- paste("omega(", 0:(s[j] - 1), ").",
namexreg[j], sep = "")
}
}
}
}
else if (straightRegress) {
xreglist$omega <- as.list(omega)
names(xreglist$omega) <- namexreg
}
allpars[[i]]$xreglist <- xreglist
mmu <- allpars[[i]]$muHat
if(straightRegress) {
resR <- z - xr %*% omega
resRdiff <- y - xreg %*% omega - mmu
}
else {
resR <- funcTF(y = z, x = as.vector(xr), delta = delta, omega = omega, b = b,
rx = r, sx = s, nx = length(z), meanval = 0)$y ##regression residuals
resRdiff <- funcTF(y = y, x = as.vector(xreg), delta = delta, omega = omega,
b = b, rx = r, sx = s, nx = length(y), meanval = 0)$y
}
res <- DLResiduals(rr, resRdiff)
sigma2muse <- getsigma2muhat(resRdiff, phi = allpars[[i]]$phi, theta = allpars[[i]]$theta,
thetaseas = allpars[[i]]$thetaseas, phiseas = allpars[[i]]$phiseas, period = period,
dfrac = allpars[[i]]$dfrac, dfs = allpars[[i]]$dfs, H = allpars[[i]]$H,
Hs = allpars[[i]]$Hs, alpha = allpars[[i]]$alpha, alphas = allpars[[i]]$alphas,
mean.dof = mean.dof)
allpars[[i]]$fitted <- resRdiff - res
} else {
yy <- y - allpars[[i]]$muHat
resR <- NULL
sigma2muse <- getsigma2muhat(yy, phi = allpars[[i]]$phi, theta = allpars[[i]]$theta,
thetaseas = allpars[[i]]$thetaseas, phiseas = allpars[[i]]$phiseas, period = period,
dfrac = allpars[[i]]$dfrac, dfs = allpars[[i]]$dfs, H = allpars[[i]]$H,
Hs = allpars[[i]]$Hs, alpha = allpars[[i]]$alpha, alphas = allpars[[i]]$alphas,
mean.dof = mean.dof)
res <- DLResiduals(rr, yy)
allpars[[i]]$fitted <- yy - res
}
allpars[[i]]$regResiduals <- resR
allpars[[i]]$residuals <- res
allpars[[i]]$sigma2 <- sigma2muse
if (getHess && !is.null(xreg) && sum(s) > 0) {
inds <- p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + sum(r) + 1:(sum(s) + 1)
XX <- cbind(xreg, rep(1, length(z)))
XXi <- solve(t(XX)%*% XX)
V <- toeplitz(rr)
covTS <- XXi %*% t(XX) %*% V %*% XX %*% XXi
allpars[[i]]$se[inds] <- sqrt(diag(sigma2muse*covTS))
}
class(allpars[[i]]) <- "ARFIMA"
}
ord <- order(logl, decreasing = TRUE)
allpars <- allpars[ord]
index <- 1:num
weeded <- (num == 1)
ans <- list(z = z, differencing = differencing, dint = dint, dseas = dseas, period = period,
modes = allpars, dmean = dmean, itmean = itmean, p = p, q = q, pseas = pseas,
qseas = qseas, lmodel = lmodel, slmodel = slmodel, weeded = weeded, getHess = getHess,
numvars = numvars, numcut = 0, n = length(z), xreg = xr, r = r, s = s, b = b,
call = match.call(), flag0 = flag0, numeach = numeach, strReg = straightRegress,
regOnly = regOnly, intindex = intindex, intname = intname, namexreg=namexreg,
fixx = fixx, indfixx = indfixx)
}
class(ans) <- "arfima"
if (autoweed && !ans$weeded) {
ans <- weed(ans, type = weedtype[1], eps2 = weedeps, adapt = adapt, pn = weedp)
}
ans
}
# Tells whether points share a circle in the p-norm
pcircle <- function(x1, x2, rad, pn = 2) {
zz <- x1 - x2
z <- (sum(abs(zz)^pn))^(1/pn)
return(z > rad)
}
# Weeds out similar (or the same) modes
#' Weeds out fits from a call to arfima that are too close to each other.
#'
#' Weeds out fits from a call to arfima that are too close to each other.
#'
#'
#' @param ans The result from a call to arfima.
#' @param type The space to perform the weeding in. "A" is for operating
#' parameters. "P" is in the PACF space. "B" performs weeding in both. "N"
#' performs no weeding and is only used internally.
#' @param walls If more than one mode is on a wall in the PACF space, all modes
#' but the one with the highest log-likelihood on the same wall are deleted.
#' @param eps2 The maximum distance between modes that are close together for
#' the mode with the lower log-likelihood to be weeded out. If adapt is TRUE
#' (default) this value changes.
#' @param eps3 The minimum distance from a wall for a secondary mode to be
#' weeded out, if walls are TRUE.
#' @param adapt If TRUE, if dim is the dimensionality of the search, eps2 is
#' changed to \eqn{(1 + eps2)^{dim} - 1}.
#' @param pn The p in the p-norm to be used in the weeding. p = 2 (default) is
#' Euclidean distance.
#' @return An object of class "arfima" with modes possibly weeded out.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}, \code{\link{distance}}
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(1234)
#' sim <- arfima.sim(1000, model = list(theta = 0.9, dfrac = 0.4))
#' fit <- arfima(sim, order = c(0, 0, 1), autoweed = FALSE, back=TRUE)
#' fit
#' distance(fit)
#' fit1 <- weed(fit)
#' fit1
#' distance(fit1)
#'}
#' @export weed
weed <- function(ans, type = c("A", "P", "B", "N"), walls = FALSE, eps2 = 0.025, eps3 = 0.01,
adapt = TRUE, pn = 2) {
if (!inherits(ans, "arfima"))
stop("weed only defined for arfima objects")
allpars <- ans$modes
if(length(allpars)==1) {
ans$weeded <- TRUE
return(ans)
}
if (ans$flag0){
stopifnot(length(allpars[[1]]$pars)==0)
allpars <- allpars[[1]]
ans$modes <- allpars
ans$weeded <- TRUE
return(ans)
}
if (length(pn) == 0 || pn <= 0 || !is.finite(pn)) {
warning("invalid pn in weed: setting to 2")
pn <- 2
}
if (length(eps2) == 0 || eps2 < 0) {
warning("invalid epsilon in weed: setting to 0.025")
eps2 <- 0.025
}
type <- type[1]
if (!(type %in% c("A", "P", "B", "N"))) {
warning("invalid type in weed: setting to \"A\"")
type <- "A"
}
if (length(eps3) == 0 || eps3 <= 0) {
warning("invalid wall epsilon in weed: setting to 0.01")
eps3 <- 0.01
}
if (!is.logical(walls)) {
warning("invalid walls selection in weed: setting to FALSE")
walls <- FALSE
}
modes <- vector("list", 1)
p <- ans$p
q <- ans$q
pseas <- ans$pseas
qseas <- ans$qseas
lmodel <- ans$lmodel
slmodel <- ans$slmodel
numvars <- p + q + qseas + pseas + if (lmodel != "n")
1 else 0
numvars <- numvars + if (slmodel != "n")
1 else 0
numvars <- numvars + sum(ans$r) + sum(ans$s)
if(numvars==0) {
stop("Should not get here. Please contact the package maintainer.")
}
if (adapt)
eps2 <- (1 + eps2)^numvars - 1
num <- length(allpars)
logl <- rep(0, num)
for (i in 1:num) logl[i] <- allpars[[i]]$loglik
ord <- order(logl, decreasing = TRUE)
allpars <- allpars[ord]
if (type %in% c("A", "B")) {
modes <- vector("list", 1)
modes[[1]] <- allpars[[1]]
flag <- NULL
if (num > 1)
for (i in 2:num) {
flag <- NULL
for (j in 1:length(modes)) {
flag <- c(flag, pcircle(x1 = allpars[[i]]$pars, x2 = modes[[j]]$pars,
rad = eps2, pn = pn))
}
if (all(flag)) {
modes[[length(modes) + 1]] <- allpars[[i]]
}
}
} else modes <- allpars
num <- length(modes)
logl <- rep(0, num)
for (i in 1:num) logl[i] <- modes[[i]]$loglik
ord <- order(logl, decreasing = TRUE)
modes <- modes[ord]
if (type %in% c("P", "B")) {
modes1 <- modes
modes <- vector("list", 1)
modes[[1]] <- modes1[[1]]
flag <- NULL
if (!ans$regOnly && num > 1) {
for (i in 2:num) {
flag <- NULL
for (j in 1:length(modes)) {
flag <- c(flag, pcircle(x1 = modes1[[i]]$parsp, x2 = modes[[j]]$parsp,
rad = eps2, pn = pn))
}
if (all(flag)) {
modes[[length(modes) + 1]] <- modes1[[i]]
}
}
}
else if(num>1) {
for (i in 2:num) {
flag <- NULL
for (j in 1:length(modes)) {
flag <- c(flag, pcircle(x1 = modes1[[i]]$omega, x2 = modes[[j]]$omega,
rad = eps2, pn = pn))
}
if (all(flag)) {
modes[[length(modes) + 1]] <- modes1[[i]]
}
}
}
}
wall <- inds <- vector("list", numvars * 2)
pl1 <- p2 <- 0
arma_varnum <- p + q + pseas + qseas
for (i in 1:numvars) {
if(ans$regOnly) break
if (i <= arma_varnum) {
minner = -1
maxxer = 1
} else if (lmodel != "n" && !pl1) {
pl1 <- 1
if (lmodel == "d") {
minner = -1
maxxer = 0.5
} else if (lmodel == "g") {
minner = 0
maxxer = 1
} else if (lmodel == "h") {
minner = 0
maxxer = 3
}
} else {
if (slmodel == "d") {
minner = -1
maxxer = 0.5
} else if (slmodel == "g") {
minner = 0
maxxer = 1
} else if (slmodel == "h") {
minner = 0
maxxer = 3
}
}
for (j in 1:length(modes)) {
wall[[2 * i - 1]] <- c(wall[[2 * i - 1]], abs(maxxer - modes[[j]]$parsp[i]) <
eps3)
wall[[2 * i]] <- c(wall[[2 * i]], abs(minner - modes[[j]]$parsp[i]) < eps3)
}
inds[[2 * i - 1]] <- which(wall[[2 * i - 1]])
inds[[2 * i]] <- which(wall[[2 * i]])
}
index <- numeric(0)
for (i in 1:(2 * numvars)) index <- c(index, inds[[i]])
index <- unique(index)
numcut <- length(ans$index) - length(index)
ans1 <- list(z = ans$z, differencing = ans$differencing, dint = ans$dint, dseas = ans$dseas,
period = ans$period, modes = modes, dmean = ans$dmean, itmean = ans$itmean, getHess = ans$getHess,
p = p, q = q, pseas = pseas, qseas = qseas, lmodel = lmodel, slmodel = slmodel,
numcut = numcut, index = index, weeded = TRUE, n = ans$n, xreg = ans$xreg, r = ans$r,
s = ans$s, b = ans$b, call = ans$call, flag0 = ans$flag0, numeach = ans$numeach,
strReg = ans$strReg, namexreg = ans$namexreg, regOnly = ans$regOnly, intindex = ans$intindex,
intname = ans$intname, fixx = ans$fixx, indfixx = ans$indfixx)
class(ans1) <- "arfima"
ans1
}
# computes and reports the distance between all modes
#' The distance between modes of an \code{arfima} fit.
#'
#' The distance between modes of an \code{arfima} fit.
#'
#'
#' @param ans An object of class "arfima".
#' @param p The p in the p-norm to be used.
#' @param digits The number of digits to print.
#' @return A list of two data frames: one with distances in operator space, the
#' second with distances in the transformed (PACF) space.
#' @author JQ (Justin) Veensta
#' @references Veenstra, J.Q. Persistence and Antipersistence: Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(8564)
#' sim <- arfima.sim(1000, model = list(phi = c(0.2, 0.1), dfrac = 0.4, theta = 0.9))
#' fit <- arfima(sim, order = c(2, 0, 1), back=TRUE)
#'
#' fit
#'
#' distance(fit)
#'}
#' @export distance
distance <- function(ans, p = 2, digits = 4) {
if (!inherits(ans, "arfima"))
stop("distance only defined for arfima objects")
if (p <= 0)
stop("invalid p")
num <- length(ans$modes)
distances1 <- distances2 <- matrix(0, nrow = num, ncol = num)
if (num == 1)
stop("only one mode")
for (i in 1:(num)) {
for (j in 1:num) {
vals <- ans$modes[[i]]$pars - ans$modes[[j]]$pars
distances1[i, j] <- (sum(abs(vals)^p))^(1/p)
vals <- ans$modes[[i]]$parsp - ans$modes[[j]]$parsp
distances2[i, j] <- (sum(abs(vals)^p))^(1/p)
}
}
distances1 <- round(distances1, digits = digits)
distances2 <- round(distances2, digits = digits)
rownames(distances1) <- colnames(distances1) <- rownames(distances2) <- colnames(distances2) <- paste("Mode",
1:num)
distances <- vector("list", 2)
distances[[1]] <- distances1
distances[[2]] <- distances2
names(distances) <- c(paste(p, "-norm distance", sep = ""), paste(p, "-norm distance in transformed space",
sep = ""))
distances
}
#' Finds the best modes of an \code{arfima} fit.
#'
#' Finds the best modes of an \code{arfima} fit with respect to log-likelihood.
#'
#' This is the easiest way to remove modes with lower log-likelihoods.
#'
#' @param object An object of class "arfima".
#' @param bestn The top number of modes to keep with respect to the
#' log-likelihood.
#' @return The \code{bestn} "best" modes.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(8765)
#' sim <- arfima.sim(1000, model = list(phi = 0.4, theta = 0.9, dfrac = 0.4))
#' fit <- arfima(sim, order = c(1, 0, 1), back=TRUE)
#' fit
#' fit <- bestModes(fit, 2)
#' fit
#'}
#' @export bestModes
bestModes <- function(object, bestn) {
if (!inherits(object, "arfima"))
stop("bestModes only defined for arfima objects")
if (!object$weeded)
stop("Please weed the object first before calling bestModes.")
m <- length(object$modes)
if (bestn > m)
stop("bestn is larger than the number of modes")
modes <- object$modes
logls <- sapply(1:m, function(i) modes[[i]]$logl)
ord <- order(logls, decreasing = TRUE)
ord <- ord[1:bestn]
modes <- modes[ord]
object$modes <- modes
object <- weed(object, eps2 = 0)
object
}
JQVeenstra/arfima documentation built on Jan. 31, 2024, 12:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions bestModes distance weed pcircle
Documented in bestModes distance weed
#' Fit ARFIMA, ARIMA-FGN, and ARIMA-PLA (multi-start) models
#'
#' Fits ARFIMA/ARIMA-FGN/ARIMA-PLA multi-start models to times series data.
#' Options include fixing parameters, whether or not to fit fractional noise,
#' what type of fractional noise (fractional Gaussian noise (FGN), fractionally
#' differenced white noise (FDWN), or the newly introduced power-law
#' autocovariance noise (PLA)), etc. This function can fit regressions with
#' ARFIMA/ARIMA-FGN/ARIMA-PLA errors via the xreg argument, including dynamic
#' regression (transfer functions).
#'
#' Fits by direct optimization using optim. The optimizer choices are: 0 -
#' BFGS; 1 - Nealder-Mead; 2 - SANN; otherwise CG.
#'
#' A word of warning: it is generally better to use the default, and only use
#' Nelder-Mead to check for spurious modes. SANN takes a long time (and may
#' only find one mode), and CG may not be stable.
#'
#' If using Nelder-Mead, it must be stressed that Nelder-Mead can take out
#' non-spurious modes or add spurious modes: we have checked visually where we
#' could. Therefore it is wise to use BFGS as the default and if there are
#' modes close to the boundaries, check using Nelder-Mead.
#'
#' The moving average parameters are in the Box-Jenkins convention: they are
#' the negative of the parameters given by \code{\link{arima}}. That is, the
#' model to be fit is, in the case of a non-seasonal ARIMA model, phi(B)
#' (1-B)^d z[t] = theta(B) a[t], where phi(B) = 1 - phi(1) B - ... - phi(p) B^p
#' and theta(B) = 1 - theta(1) B - ... - theta(q) B^q.
#'
#' For the useC parameter, a "0" means no C is used; a "1" means C is only used
#' to compute the log-likelihood, but not the theoretical autocovariance
#' function (tacvf); a "2" means that C is used to compute the tacvf and not
#' the log-likelihood; and a "3" means C is used to compute everything.
#'
#' @param z The data set (time series)
#' @param order The order of the ARIMA model to be fit: c(p, d, q). We have
#' that p is the number of AR parameters (phi), d is the amount of integer
#' differencing, and q is the number of MA parameters (theta). Note we use the
#' Box-Jenkins convention for the MA parameters, in that they are the negative
#' of \code{\link{arima}}: see "Details".
#' @param numeach The number of starts to fit for each parameter. The first
#' argument in the vector is the number of starts for each AR/MA parameter,
#' while the second is the number of starts for the fractional parameter. When
#' this is set to 0, no fractional noise is fit. Note that the number of
#' starts in total is multiplicative: if we are fitting an ARFIMA(2, d, 2), and
#' use the older number of starts (c(2, 2)), we will have 2^2 * 2 * 2^2 = 32
#' starting values for the fits. \strong{Note that the default has changed
#' from c(2, 2) to c(1, 1) since package version 1.4-0}
#' @param back Setting this to true will restore the defaults in numeach.
#' @param dmean Whether the mean should be fit dynamically with the optimizer.
#' Note that the likelihood surface will change if this is TRUE, but this
#' is usually not worrisome.
#' See the referenced thesis for details.
#' @param whichopt Which optimizer to use in the optimization: see "Details".
#' @param itmean This option is under investigation, and will be set to FALSE
#' automatically until it has been decided what to do.
#'
#' Whether the mean should be fit iteratively using the function
#' \code{\link{TrenchMean}}. Currently itmean, if set to TRUE, has higher
#' priority that dmean: if both are TRUE, dmean will be set to FALSE, with a
#' warning.
#' @param fixed A list of parameters to be fixed. If we are to fix certain
#' elements of the AR process, for example, fixed$phi must have length equal to
#' p. Any numeric value will fix the parameter at that value; for example, if
#' we are modelling an AR(2) process, and we wish to fix only the first
#' autoregressive parameter to 0, we would have fixed = list(phi = c(0, NA)).
#' NA corresponds to that parameter being allowed to change in the optimization
#' process. We can fix the fractional parameters, and unlike
#' \code{\link{arima}}, can fix the seasonal parameters as well. Currently,
#' fixing regression/transfer function parameters is disabled.
#' @param lmodel The long memory model (noise type) to be used: "d" for FDWN,
#' "g" for FGN, "h" for PLA, and "n" for none (i.e. ARMA short memory models).
#' Default is "d".
#' @param seasonal The seasonal components of the model we wish to fit, with
#' the same components as above. The period must be supplied.
#' @param useC How much interfaced C code to use: an integer between 0 and 3.
#' The value 3 is strongly recommended. See "Details".
#' @param cpus The number of CPUs used to perform the multi-start fits. A
#' small number of fits and a high number of cpus (say both equal 4) with n not
#' large can actually be slower than when cpus = 1. The number of CPUs should
#' not exceed the number of threads available to R.
#' @param rand Whether random starts are used in the multistart method.
#' Defaults to FALSE.
#' @param numrand The number of random starts to use.
#' @param seed The seed for the random starts.
#' @param eps3 How far to start from the boundaries when using a grid for the
#' multi-starts (i.e. when rand is FALSE.)
#' @param xreg A matrix, data frame, or vector of regressors for regression or
#' transfer functions.
#' @param reglist A list with the following elements:
#' \itemize{
#' \item regpar -
#' either NA or a list, matrix, data frame, or vector with 3 columns. If
#' regpar is a vector, the matrix xreg must have one row or column only. In
#' order, the elements of regpar are: r, s, and b. The values of r are the the
#' orders of the delta parameters as in Box, Jenkins and Reinsel, the values of
#' s are the orders of omega parameters, and the values of b are the
#' backshifting to be done.
#'
#' \item minn - the minimum value for the starting value of the search, if
#' reglist$numeach > 1.
#' \item maxx - the maximum value for the starting value
#' of the search, if reglist$numeach > 1.
#' \item numeach - the number of starts
#' to try for each regression parameter.
#'
#' }
#' @param check If TRUE, checks at each optim iteration whether the model is
#' identifiable. This makes the optimization much slower.
#' @param autoweed Whether to automatically (before the fit is returned) weed
#' out modes found that are found that are close together (usually the same
#' point.)
#' @param weedeps The maximum distance between modes that are close together
#' for the mode with the lower log-likelihood to be weeded out. If adapt is
#' TRUE (default) this value changes.
#' @param adapt If TRUE, if dim is the dimensionality of the search, weedeps is
#' changed to \eqn{(1 + weedeps)^{dim} - 1}.
#' @param weedtype The type of weeding to be done. See \code{\link{weed}}.
#' @param weedp The p in the p-norm to be used in the weeding. p = 2 (default)
#' is Euclidean distance.
#' @param quiet If TRUE, no auxiliary output is generated. The default (FALSE)
#' has information of fits being proformed.
#' @param startfit Meant primarily for debugging (for now), allows starting places
#' for the fitting process. Overrides \code{numeach}.
#' @return An object of class "arfima". In it, full information on the fit is
#' given, though not printed under the print.arfima method. The phis are the
#' AR parameters, and the thetas are the MA parameters. Residuals, regression
#' residuals, etc., are all available, along with the parameter values and
#' standard errors. Note that the muHat returned in the arfima object is
#' of the \strong{differenced} series, if differencing is applied.
#'
#' Note that if multiple modes are found, they are listed in order of
#' log-likelihood value.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima.sim}}, \code{\link{SeriesJ}},
#' \code{\link{arfima-package}}
#' @references McLeod, A. I., Yu, H. and Krougly, Z. L. (2007) Algorithms for
#' Linear Time Series Analysis: With R Package Journal of Statistical Software,
#' Vol. 23, Issue 5
#'
#' 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
#'
#'\donttest{
#' set.seed(8564)
#' sim <- arfima.sim(1000, model = list(phi = c(0.2, 0.1),
#' dfrac = 0.4, theta = 0.9))
#' fit <- arfima(sim, order = c(2, 0, 1), back=TRUE)
#'
#' fit
#'
#' data(tmpyr)
#'
#' fit <- arfima(tmpyr, order = c(1, 0, 1), numeach = c(3, 3))
#' fit
#'
#' plot(tacvf(fit), maxlag = 30, tacf = TRUE)
#'
#' data(SeriesJ)
#' attach(SeriesJ)
#'
#' fitTF <- arfima(YJ, order= c(2, 0, 0), xreg = XJ, reglist =
#' list(regpar = c(2, 2, 3)), lmodel = "n")
#' fitTF
#'
#' detach(SeriesJ)
#'}
#'
#' @export arfima
"arfima" <- function(z, order = c(0, 0, 0), numeach = c(1, 1), dmean = TRUE, whichopt = 0,
itmean = FALSE, fixed = list(phi = NA, theta = NA, frac = NA, seasonal = list(phi = NA,
theta = NA, frac = NA), reg = NA), lmodel = c("d", "g", "h", "n"), seasonal = list(order = c(0,
0, 0), period = NA, lmodel = c("d", "g", "h", "n"), numeach = c(1, 1)), useC = 3,
cpus = 1, rand = FALSE, numrand = NULL, seed = NA, eps3 = 0.01, xreg = NULL, reglist = list(regpar = NA,
minn = -10, maxx = 10, numeach = 1), check = F, autoweed = TRUE, weedeps = 0.01,
adapt = TRUE, weedtype = c("A", "P", "B"), weedp = 2, quiet = FALSE, startfit = NULL, back = FALSE) {
if (itmean) {
warning("itmean is currently under examination, and is not useable at the moment.\nSetting itmean to FALSE")
itmean <- FALSE
}
if(back) {
numeach <- c(2, 2)
seasonal$numeach <- c(2, 2)
dmean <- TRUE
}
if(!is.null(startfit)) {
stopifnot(is.matrix(startfit))
starts <- startfit
num <- nrow(starts)
startflag <- TRUE
}
else startflag <- FALSE
if (any(is.na(z))) {
stop("Missing data not supported at the moment.")
}
cpus <- abs(round(cpus))
if (cpus < 1)
cpus <- 1
getHess <- TRUE ## to save on options but not change code
inclbdries <- FALSE ## to save on options but not change code
tse <- FALSE ## to save on options but not change code
regpar <- reglist$regpar
straightRegress <- FALSE
namexreg <- NULL
intindex <- 0
intname <- NULL
if (!is.null(regpar) && all(is.na(regpar)))
regpar <- NULL
if (!is.null(xreg)) {
if (!(is.data.frame(xreg) || is.matrix(xreg) || is.vector(xreg) || is.ts(xreg)))
stop("Invalid xreg type")
if (any(is.na(xreg)))
stop("no missing values allowed in xreg: if dynamic regression needs padding, please pad with zeroes")
constantInds <- numeric(0)
if (is.null(regpar)) {
if (!is.data.frame(xreg) && !is.matrix(xreg)) {
xreg <- data.frame(xreg)
colnames(xreg) <- as.character(as.list((match.call()))$xreg)
}
namexreg <- colnames(xreg)
if (is.null(namexreg))
namexreg <- paste(rep("X", ncol(xreg)), 1:ncol(xreg), sep = "")
regpar <- matrix(0, nrow = 3, ncol = dim(xreg)[2])
regpar <- as.data.frame(t(regpar))
names(regpar) <- c("r", "s", "b")
constantInds <- c(1:dim(xreg)[2])[apply(xreg, 2, var) == 0]
straightRegress <- ncol(xreg)
if(length(constantInds)) {
for(i in constantInds) {
if(all(xreg[, i] == 1)) {
constantInds <- setdiff(constantInds, i)
if(i > 1) {
xreg <- xreg[, c(i, setdiff(1:ncol(xreg), i))]
namexreg <- namexreg[c(i, setdiff(1:ncol(xreg), i))]
intindex <- 1i
if(length(constantInds)>0)
constantInds <- c(1:dim(xreg)[2])[apply(xreg, 2, var) == 0][-1]
}
intindex <- 1
break
}
}
straightRegress <- straightRegress - length(constantInds) - 1
}
straightRegress <- (straightRegress>0)
if(!straightRegress) {
warning("All regression terms constant. Setting xreg to NULL")
xreg <- NULL
}
else if(straightRegress && dmean != TRUE) {
warning("Regression with time series errors requires dmean=T, setting")
dmean <- TRUE
}
else if(straightRegress && intindex && is.numeric(dmean)) {
warning("with an intercept term, dmean must be TRUE. Setting it so.")
dmean <- TRUE
}
}
else {
if (is.matrix(regpar) || is.data.frame(regpar)) {
if (dim(regpar)[2] != dim(xreg)[2]) {
stop("xreg and regpar sizes do not match")
}
namexreg <- colnames(xreg)
if (is.null(namexreg))
namexreg <- paste(rep("X", ncol(xreg)), 1:ncol(xreg), sep = "")
if (any(is.na(regpar)))
stop("no missing values allowed in regpar")
if (dim(regpar)[1] != 3)
stop("regpar does not have 3 rows")
if (!is.null(names(regpar)) && names(regpar) != c("r", "s", "b"))
stop("names(regpar), if not NULL, must be r, s, and b")
regpar <- as.data.frame(t(regpar))
names(regpar) <- c("r", "s", "b")
} else if (is.vector(xreg) || is.ts(xreg)) {
pars <- data.frame(r = regpar[1], s = regpar[2], b = regpar[3])
regpar <- pars
xreg <- data.frame(xreg)
colnames(xreg) <- as.character(as.list((match.call()))$xreg)
namexreg <- colnames(xreg)
if(is.null(namexreg))
namexreg <- "X1"
} else stop("invalid regpar type")
}
xreg <- as.matrix(xreg)
if (length(constantInds) > 0) {
warning(paste("columns", paste(constantInds, collapse = ", "), "are constant and will be removed"))
xreg <- xreg[, -constantInds]
namexreg <- namexreg[-constantInds]
}
if(intindex && is.logical(dmean) && dmean) {
regpar <- regpar[-1, ]
xreg <- xreg[, -1]
intname <- namexreg[1]
namexreg <- namexreg[-1]
}
else if(intindex) {
warning(paste(paste("there is a constant column named", namexreg[1]),
"take to be the intercept, but dmean is not TRUE. This column",
"will be removed and ignored.", sep = "\n"))
regpar <- regpar[-1, ]
xreg <- xreg[, -1]
intindex <- 0
}
else if(is.logical(dmean) && dmean) {
intname <- "Intercept"
}
s <- regpar$s <- regpar$s + 1
r <- regpar$r
b <- regpar$b
if (any(r > 0) || any(s>1)) {
#cat("note: transfer functions do not work with dynamic mean: setting dmean to FALSE\n")
dmean <- FALSE
}
numvarreg <- sum(s) + sum(r)
regeach <- if(straightRegress) 1 else reglist$numeach
if (length(regeach) == 0)
regeach <- 1
if (length(regeach) != 1 || regeach != round(regeach) || regeach <= 0) {
warning("invalid reglist$numeach: resetting to 1")
regeach <- 1
}
regmin <- reglist$minn
regmax <- reglist$maxx
if (length(regmin) == 0 || is.na(regmin))
regmin <- if(regeach>1) -10 else 0
if (length(regmax) == 0 || is.na(regmax))
regmax <- if(regeach>1) 10 else 0
if (itmean)
stop("Please note that for regression problems, the mean is calculated separately\n via optimization or on the regression parameters or by the mean of the series;\n therefore itmean is not valid")
}
else {
r <- s <- b <- 0
numvarreg <- regeach <- 0
namexreg <- reglist <- regpar <- NULL
}
if (is.numeric(dmean) && !is.double(dmean))
stop("dmean must be logical or a double")
if (itmean && ((is.logical(dmean) && dmean) || is.numeric(dmean))) {
warning("itmean trumps dmean: setting dmean to FALSE")
dmean <- FALSE
}
mean.dof = if(itmean || (is.logical(dmean) && dmean)) 1 else 0
if(!is.null(xreg)) mean.dof <- mean.dof + dim(xreg)[2]
info <- NULL
dint <- order[2]
if (length(dint) == 0)
dint <- 0
lmodel <- lmodel[1]
if (is.null(numeach) || (length(numeach) != 2 && (lmodel != "n" && length(numeach) !=
1)) || !is.vector(numeach))
stop("numeach must be a vector of length 1 or 2:\n \n\t\tposition 1 is the values for ARMA\n parameters, position 2 is the value for the fractional parameter\n\n\t\tif a vector of length 1, we must have lmodel = \"n\"")
if (numeach[1] != round(numeach[1]) || numeach[1] <= 0) {
warning("incorrect numeach ARMA: resetting to 1")
numeach[1] <- 1
}
if (numeach[2] != round(numeach[2]) || numeach[2] < 0) {
warning("incorrect numeach fractional: resetting to 1")
numeach[2] <- 1
}
if (numeach[2] == 0 && lmodel != "n") {
cat("lmodel = \"", lmodel, "\" and numeach[2] = 0; setting lmodel = \"n\"\n", sep = "")
lmodel <- "n"
}
if (lmodel == "n") {
if (length(numeach) == 1)
numeach <- c(numeach, 0) else numeach[2] <- 0
}
if (!is.null(seasonal) && !is.na(seasonal$period) && !is.null(seasonal$period) && (seasonal$period !=
0) && ((!is.null(seasonal$order) || !all(seasonal$order == 0)) || (!is.null(seasonal$lmodel) &&
seasonal$lmodel != "n"))) {
dseas <- seasonal$order[2]
if (length(dseas) == 0)
dseas <- 0
pseas <- seasonal$order[1]
qseas <- seasonal$order[3]
if (length(pseas) == 0)
pseas <- 0
if (length(qseas) == 0)
qseas <- 0
if (!is.null(seasonal$lmodel))
slmodel <- seasonal$lmodel[1] else slmodel <- "d"
if (!is.null(seasonal$numeach))
snumeach <- seasonal$numeach else {
warning("no seasonal$numeach: setting to numeach")
snumeach <- numeach
}
if (is.null(snumeach) || (length(snumeach) != 2 && (slmodel != "n" && length(snumeach) !=
1)) || !is.vector(snumeach))
stop("seasonal$numeach ust be a vector of length 1 or 2:\n \n\t\tposition 1 is the values for ARMA\n parameters, position 2 is the value for the fractional parameter\n\n\t\tif a vector of length 1, we must have seasonal$lmodel = \"n\"")
period <- seasonal$period
if (snumeach[1] != round(snumeach[1]) || snumeach[1] <= 0) {
warning("incorrect seasonal$numeach ARMA: resetting to 1")
snumeach[1] <- 1
}
if (snumeach[2] != round(snumeach[2]) || snumeach[2] < 0) {
warning("incorrect seasonal$numeach fractional: resetting to 1")
snumeach[2] <- 1
}
if (snumeach[2] == 0 && slmodel != "n") {
cat("seasonal lmodel =", slmodel, "and seasonal numeach[2] = 0; setting seasonal lmodel = \"n\"\n")
slmodel <- "n"
}
if (slmodel == "n") {
if (length(snumeach) == 1)
snumeach <- c(snumeach, 0) else snumeach[2] <- 0
}
slmr <- if (slmodel != "n")
1 else 0
} else {
dseas <- pseas <- qseas <- period <- slmr <- 0
snumeach <- c(0, 0)
slmodel <- "n"
}
p <- order[1]
q <- order[3]
lmr <- if (lmodel != "n")
1 else 0
numvars <- p + q + pseas + qseas + lmr + slmr + numvarreg +
if (dmean == TRUE)
1 else 0
if (length(fixed) > 0) {
if (!(is.null(fixed$reg)) && !(all(is.na(fixed$reg))) && length(fixed$reg) > 0) {
warning("no fixed regression parameters for now.")
fixed$reg <- NA
}
arfixed <- fixed$phi
mafixed <- fixed$theta
fracfix <- fixed$frac
if((length(arfixed)==1 && is.na(arfixed))||length(arfixed)==0)
arfixed <- rep(NA, p)
if((length(mafixed)==1 && is.na(mafixed))||length(mafixed)==0)
mafixed <- rep(NA, q)
if(length(fracfix)==0 && lmodel !="n") {
fracfix <- NA
}
warnstring <- "In the fixed argument:"
if (length(arfixed) != p) {
warnstring <- paste0(warnstring, "\n length phi is ", length(arfixed),
" but p = ", p, " unfixing phi")
arfixed <- if (p > 0)
rep(NA, p) else numeric(0)
}
else if(p>0 && length(arfixed)==p&&any(!is.na(arfixed))) {
if(sum(is.na(arfixed))==0) {
if(!InvertibleQ(arfixed)) {
stop("phi parameters fixed at nonstationary parameters")
}
}
}
if (length(mafixed) != q) {
warnstring <- paste0(warnstring, "\n length theta is ", length(mafixed),
" but q = ", q, " unfixing theta")
mafixed <- if (q > 0)
rep(NA, q) else numeric(0)
}
else if(any(!is.na(mafixed))) {
if(sum(is.na(mafixed))==0) {
if(!InvertibleQ(mafixed)) {
stop("theta parameters fixed at noninvertible parameters")
}
}
}
if (!is.null(fracfix) && length(fracfix) != 1 && !is.na(fracfix)) {
warnstring <- paste0(warnstring, "\n fractional parameter is fixed greater",
"\n than length 1, resetting fractional prameter ")
if(lmodel=="n") {
warnstring<- paste0(warnstring,"\n lmodel == 'n' and fractional parameter fixed, resetting")
fracfix <- numeric(0)
}
else fracfix <- NA
}
if (lmodel == "n")
fracfix <- numeric(0)
if(lmodel!="n"&&(!is.null(fracfix)&&!is.na(fracfix))) {
fix_stop <- F
if(lmodel=="d") {
if(!InvertibleD(fracfix)) fix_stop <- T
}
else if(lmodel=="g") {
if(!InvertibleH(fracfix)) fix_stop <- T
}
else {
if(!InvertibleAlpha(fracfix)) fix_stop <- T
}
if(fix_stop){
stop("fractional parameter fixed at nonstationary/noninvertible value")
}
}
fixx <- c(arfixed, mafixed)
if (period > 0) {
sarfixed <- fixed$seasonal$phi
smafixed <- fixed$seasonal$theta
sfracfix <- fixed$seasonal$frac
if((length(sarfixed)==1 && is.na(sarfixed))||length(sarfixed)==0)
sarfixed <- rep(NA, pseas)
if((length(smafixed)==1 && is.na(smafixed))||length(smafixed)==0)
smafixed <- rep(NA, qseas)
if(length(sfracfix)==0 && slmodel !="n") {
sfracfix <- NA
}
if (length(sarfixed) != pseas) {
warnstring <- paste0(warnstring, "\n length seasonal phi is ", length(sarfixed),
" but pseas = ", pseas, " unfixing s_phi")
sarfixed <- if (pseas > 0)
rep(NA, pseas) else numeric(0)
}
else if(pseas>0 && length(sarfixed)==pseas&&any(!is.na(sarfixed))) {
if(sum(is.na(sarfixed))==0) {
if(!InvertibleQ(sarfixed)) {
stop("seasonal phi parameters fixed at nonstationary parameters")
}
}
}
if (length(smafixed) != qseas) {
warnstring <- paste0(warnstring, "\n length seasonal theta is ", length(smafixed),
" but qseas = ", qseas, " unfixing s_theta")
smafixed <- if (qseas > 0)
rep(NA, qseas) else numeric(0)
}
else if(qseas>0 && length(smafixed)==qseas&&any(!is.na(smafixed))) {
if(sum(is.na(smafixed))==0) {
if(!InvertibleQ(smafixed)) {
stop("seasonal theta parameters fixed at noninvertible parameters")
}
}
}
if (!is.null(sfracfix) && length(sfracfix) != 1 && !is.na(sfracfix)) {
warnstring <- paste0(warnstring, "\n seasonal fractional parameter is",
"\n fixed greater than length 1, resetting seasonal",
"\n fractional prameter ")
if(slmodel=="n") {
warnstring("\n slmodel == 'n' and seasonal fractional parameter fixed",
"\n resetting")
sfracfix <- numeric(0)
}
else sfracfix <- NA
}
if(slmodel!="n"&&(!is.null(sfracfix)&&!is.na(sfracfix))) {
fix_stop <- F
if(slmodel=="d") {
if(!InvertibleD(sfracfix)) fix_stop <- T
}
else if(slmodel=="g") {
if(!InvertibleH(sfracfix)) fix_stop <- T
}
else {
if(!InvertibleAlpha(sfracfix)) fix_stop <- T
}
if(fix_stop) {
stop("seasonal fractional parameter fixed at nonstationary/noninvertible value")
}
}
sfixx <- c(sarfixed, smafixed)
} else {
sfixx <- NULL
sfracfix <- NULL
}
if(nchar(warnstring)>31) warning(warnstring)
fixx <- c(fixx, sfixx, fracfix, sfracfix, rep(NA, sum(r)+sum(s)))
indfixx <- !is.na(fixx)
if (is.logical(dmean) && dmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
} else if (is.logical(itmean) && itmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
}
} else {
fixx <- rep(NA, numvars)
indfixx <- rep(FALSE, numvars)
if (is.logical(dmean) && dmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
} else if (is.logical(itmean) && itmean) {
fixx <- c(fixx, NA)
indfixx <- c(indfixx, FALSE)
}
}
numtries <- 1
if (p + q > 0)
numtries <- numtries * numeach[1]^(p + q)
if (numeach[2] > 1)
numtries <- numtries * numeach[2]
if (pseas + qseas > 0)
numtries <- numtries * snumeach[1]^(pseas + qseas)
if (snumeach[2] > 1)
numtries <- numtries * snumeach[2]
if (regeach > 1)
numtries <- numtries * regeach
if (sum(as.integer(indfixx)) > 0) {
if (p > 0 && any(!is.na(arfixed)))
numtries <- numtries/(numeach[1]^(sum(!is.na(arfixed))))
if (q > 0 && any(!is.na(mafixed)))
numtries <- numtries/(numeach[1]^(sum(!is.na(mafixed))))
if (lmodel != "n" && (is.null(fracfix)||!is.na(fracfix)))
numtries <- numtries/numeach[2]
if (pseas > 0 && any(!is.na(sarfixed)))
numtries <- numtries/(snumeach[1]^(sum(!is.na(sarfixed))))
if (qseas > 0 && any(!is.na(smafixed)))
numtries <- numtries/(snumeach[1]^(sum(!is.na(smafixed))))
if (slmodel != "n" && (is.null(sfracfix)||!is.na(sfracfix)))
numtries <- numtries/snumeach[2]
}
# if(numvars == 0 && numtries > 1) stop('invalid specification of variables')
if (round(dint) != dint || dint < 0)
stop("d must be an integer >= 0")
if (round(dseas) != dseas || dseas < 0)
stop("dseas must be an integer >= 0")
if (!quiet) {
if(!startflag) {
if (autoweed && numtries > 1 && !startflag)
cat("Note: autoweed is ON. It is possible, but not likely,\n that unique modes may be lost.\n") else if (!autoweed && numtries > 1)
cat("Note: autoweed is OFF. Running weed after the\n fit is highly recommended.\n") else {
cat("Note: only one starting point. Only one mode can be found -- this is now the default behavior.\n")
}
}
else {
if(autoweed)
cat("Note: autoweed is ON. It is possible, but not likely,\n that unique modes may be lost.\n")
else
cat("Note: autoweed is OFF. Running weed after the\n fit may be required.\n")
cat("Starting fit with given start matrix.\n")
}
}
y <- z
xr <- xreg
if (dint > 0) {
y <- diff(y, differences = dint)
if (!is.null(xreg))
xreg <- diff(xreg, differences = dint)
}
if (dseas > 0) {
y <- diff(y, differences = dseas, lag = period)
if (!is.null(xreg))
xreg <- diff(xreg, differences = dseas, lag = period)
}
D <- dint + dseas * period
if (any(r > 0) && any(abs(colMeans(xreg)) > 1e-14)) {
cat("Please note for transfer functions the means of each X variable must be 0: subtracting mean from each X\n")
for (i in 1:ncol(xreg)) {
xreg[, i] <- xreg[, i] - colMeans(xreg)[i]
}
}
flag0 <- regOnly <- FALSE
if (p + q + numeach[2] + pseas + qseas + snumeach[2] == 0) {
if(is.null(xreg))
flag0 <- TRUE
else
regOnly <- TRUE
}
differencing <- if (D > 0) TRUE else FALSE
if (is.logical(dmean) && !dmean && !straightRegress)
meanval <- mean(y)
else if (is.numeric(dmean))
meanval <- dmean
else
meanval <- 0
sumfixed <- sum(indfixx)
numvars1 <- numvars - sumfixed
if (numvars1 == 0)
getHess <- FALSE
if (flag0) {
if(sum(numeach)>2||sum(seasonal$numeach)>2||
(!is.null(numrand)) && numrand>0)
warning('With only the mean to be fit, since the problem is convex,
it has a unique minimum and thus only one fit will be performed.')
if ((is.logical(dmean) && dmean) || itmean) {
fun <- function(m) {
lARFIMA(y - m)
}
fit <- optim(mean(y), fn = fun, method = "Brent", lower = mean(y) - 10, upper = mean(y) +
10, control = list(trace = 0, fnscale = -1), hessian = getHess)
hess <- fit$hessian[1, 1]
se <- tryCatch(sqrt(diag(solve(-hess))), error = function(err) rep(NA, dim(hess)[1]))
res <- list(maximum = fit[[1]][1], objective = fit[[2]][1], hess = hess, se = se)
} else if (is.logical(dmean)) {
res <- list(maximum = mean(y), objective = lARFIMA(y - mean(y)), hess = NA,
se = -Inf)
} else {
res <- list(maximum = dmean, objective = lARFIMA(y - dmean), hess = NA, se = -Inf)
}
allpars <- vector("list", 1)
allpars[[1]]$muHat <- res$maximum
allpars[[1]]$loglik <- res$objective
allpars[[1]]$hessian <- res$hess
allpars[[1]]$se <- res$se
allpars[[1]]$phi <- allpars[[1]]$theta <- allpars[[1]]$phiseas <-
allpars[[1]]$thetaseas <- allpars[[1]]$dfrac <- allpars[[1]]$dfs <-
allpars[[1]]$phip <- allpars[[1]]$thetap <- allpars[[1]]$H <-
allpars[[1]]$alpha <- allpars[[1]]$alphas <- allpars[[1]]$delta <-
allpars[[1]]$omega <- allpars[[1]]$Hs <- numeric(0)
allpars[[1]]$pars <-
c(allpars[[1]]$phi, allpars[[1]]$theta, allpars[[1]]$phiseas,
allpars[[1]]$thetaseas, allpars[[1]]$dfrac, allpars[[1]]$H,
allpars[[1]]$alpha, allpars[[1]]$dfs, allpars[[1]]$Hs,
allpars[[1]]$alphas, allpars[[1]]$delta, allpars[[1]]$omega)
sigma2muse <- getsigma2muhat(y - allpars[[1]]$muHat, mean.dof = mean.dof)
allpars[[1]]$sigma2 <- sigma2muse
rr <- tacvfARFIMA(maxlag = length(y) - 1)
res <- DLResiduals(rr, y - allpars[[1]]$muHat)
allpars[[1]]$residuals <- res
allpars[[1]]$fitted <- y - allpars[[1]]$muHat - res
class(allpars[[1]]) <- "ARFIMA"
ans <- list(z = z, differencing = differencing, dint = dint, dseas = dseas, period = period,
modes = allpars, dmean = dmean, itmean = itmean, p = p, q = q, pseas = pseas,
qseas = qseas, lmodel = lmodel, slmodel = slmodel, weeded = TRUE, getHess = getHess,
numvars = numvars, numcut = 0, n = length(z), xreg = xr, r = r, s = s, b = b,
call = match.call(), flag0 = flag0, numeach = numeach, strReg = F,
regOnly = regOnly, intindex = intindex, intname = intname, namexreg=namexreg,
fixx = fixx, indfixx = indfixx)
} else {
num <- if (!rand)
numtries else numrand
if(startflag) num <- nrow(starts)
if (rand && !startflag) {
if (length(num) == 0) {
warning("random starts selected and no number of starts selected: setting to 8")
num <- 8
}
starts <- matrix(0, nrow = num, ncol = numvars)
pl1 <- 0
if (!is.na(seed) && !is.null(seed))
set.seed(seed) else warning("Random start seed not selected: results may not be reproducible")
if (p + q + pseas + qseas > 0) {
starts[, 1:(p + q + pseas + qseas)] <- matrix(runif(num * (p + q + pseas +
qseas), min = -1, max = 1), ncol = (p + q + pseas + qseas))
}
if (lmodel != "n") {
pl1 <- 1
if (lmodel == "d")
starts[, (1 + p + q + pseas + qseas)] <- runif(1, min = -1, max = 0.5)
else if (lmodel == "g")
starts[, (1 + p + q + pseas + qseas)] <- runif(1, min = 0, max = 1)
else
starts[, (1 + p + q + pseas + qseas)] <- runif(1, min = 0, max = 3)
}
if (slmodel != "n") {
if (slmodel == "d")
starts[, (1 + p + q + pl1 + pseas + qseas)] <- runif(1, min = -1, max = 0.5)
else if (slmodel == "g")
starts[, (1 + p + q + pl1 + pseas + qseas)] <- runif(1, min = 0, max = 1)
else
starts[, (1 + p + q + pl1 + pseas + qseas)] <- runif(1, min = 0, max = 3)
}
if (!is.null(xreg)) {
if(!straightRegress) {
starts[, (1 + p + q + pl1 + pseas + qseas):(1 + p + q + pl1 + pseas + qseas +
numvarreg)] <- matrix(runif(num * (numvarreg), min = regmin,
max = regmax), nrow = num)
}
else if(is.logical(dmean)&&dmean) {
coeffs <- lsfit(xreg, y)$coef
for(ii in 1:numvarreg) starts[, (1 + p + q + pl1 + pseas + qseas+ii)] <- coeffs
}
else {
coeffs <- lsfit(xreg, y, intercept = FALSE)$coef
for(ii in 1:numvarreg) starts[, (1 + p + q + pl1 + pseas + qseas+ii)] <- coeffs
}
}
} else if (num > 1 && !startflag) {
pl1 <- 0
li <- vector("list", numvars)
if (p + q > 0) {
if (inclbdries)
seqq <- seq(-1 + eps3, 1 - eps3, length.out = numeach[1]) else {
seqq <- seq(-1 + eps3, 1 - eps3, length.out = numeach[1] + 2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
}
if (pseas + qseas > 0) {
if (inclbdries)
seqq <- seq(-1 + eps3, 1 - eps3, length.out = snumeach[1]) else {
seqq <- seq(-1 + eps3, 1 - eps3, length.out = snumeach[1] + 2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
}
for (i in 1:(p + q + pseas + qseas))
li[[i]] <- seqq
if (lmodel != "n") {
pl1 <- 1
if (lmodel == "d") {
minner <- -1
maxxer <- 0.5
} else if (lmodel == "g") {
minner <- 0
maxxer <- 1
} else {
minner <- 0
maxxer <- 3
}
if (inclbdries)
seqq <- seq(minner + eps3, maxxer - eps3, length.out = numeach[2]) else {
seqq <- seq(minner + eps3, maxxer - eps3, length.out = numeach[2] +
2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
li[[(p + q + pl1 + pseas + qseas)]] <- seqq
}
pl1_1 <- 0
if (slmodel != "n") {
pl1_1 <- 1
if (slmodel == "d") {
minner <- -1
maxxer <- 0.5
} else if (slmodel == "g") {
minner <- 0
maxxer <- 1
} else {
minner <- 0
maxxer <- 3
}
if (inclbdries)
seqq <- seq(minner + eps3, maxxer - eps3, length.out = snumeach[2]) else {
seqq <- seq(minner + eps3, maxxer - eps3, length.out = snumeach[2] +
2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
li[[(p + q + 1 + pseas + qseas + pl1)]] <- seqq
}
if (!is.null(xreg)) {
if(!straightRegress) {
if (inclbdries)
seqq <- seq(regmin + eps3, regmax - eps3, length.out = regeach)
else {
seqq <- seq(regmin + eps3, regmax - eps3, length.out = regeach + 2)
seqq <- seqq[-1]
seqq <- seqq[-length(seqq)]
}
for (i in 1:(numvarreg)) li[[(p + q + pseas + qseas + pl1 + pl1_1) + i]] <- seqq
}
else if(is.logical(dmean)&&dmean){
coeffs <- lsfit(xreg, y)$coef[c(2:(ncol(xreg)+1), 1)]
for (i in 1:(numvarreg+1)) li[[(p + q + pseas + qseas + pl1 + pl1_1) + i]] <- coeffs[i]
}
else {
coeffs <- lsfit(xreg, y, intercept = FALSE)$coef
for (i in 1:(numvarreg)) li[[(p + q + pseas + qseas + pl1 + pl1_1) + i]] <- coeffs[i]
}
}
else if(is.logical(dmean) && dmean)
li[[length(li)]] <- mean(y)
starts <- expand.grid(li, KEEP.OUT.ATTRS = FALSE)
if(sumfixed>0) {
for(i in 1:numvars) {
if(indfixx[i]) starts[, i] <- rep(fixx[i], nrow(starts))
}
starts <- unique(starts)
}
} else if(!startflag) {
starts <- if (numvars > 0)
matrix(0, nrow=1,ncol = numvars) else matrix(0)
if(straightRegress && dmean == TRUE) {
coeffs <- lsfit(xreg, y)$coef[c(2:(ncol(xreg)+1), 1)]
starts[, (numvars-length(coeffs)+1):numvars] <- coeffs
}
else if(straightRegress) {
coeffs <- lsfit(xreg, y, intercept = FALSE)$coef
starts[, (numvars-length(coeffs)+1):numvars] <- coeffs
}
else if(dmean == TRUE)
starts[, numvars] <- mean(y)
if(sumfixed>0) {
for(i in 1:numvars) {
if(indfixx[i]) starts[, i] <- fixx[i]
}
}
}
if(dmean==F||is.numeric(dmean)) {
indfixxM <- c(indfixx, TRUE)
numvarsM <- numvars + 1
}
else {
indfixxM <- indfixx
numvarsM <- numvars
}
if (nrow(starts) != num)
stop("error in starts")
if (rand)
strg <- " the random start " else strg <- " the "
if (!quiet) {
cat("Beginning", strg, "fits with ", num, " starting values.\n", sep = "")
cat("\n")
}
fit <- vector("list", num)
if(!straightRegress)
xreg <- as.vector(xreg)
if (cpus == 1) {
for (i in 1:num) {
fit[[i]] <- arfimaFit(y = y, p = p, q = q, pseas = pseas, qseas = qseas,
period = period, lmodel = lmodel, slmodel = slmodel, whichopt = whichopt,
dmean = dmean, parstart = as.double(starts[i, ]), getHess = getHess, itmean = itmean,
indfixx = indfixx, fixx = fixx, xreg = xreg, r = r, s = s, b = b, straightRegress = straightRegress)
}
} else {
cl <- makePSOCKcluster(cpus, methods = FALSE)
numms <- 1:num
clusterExport(cl, c("lARFIMA", "PacfToAR", "DLLoglikelihood", "lARFIMAwTF"),
envir = environment(arfima))
fit <- clusterApplyLB(cl, numms, arfimaFitpar, parstart = starts, y = y, p = p,
q = q, pseas = pseas, qseas = qseas, period = period, lmodel = lmodel, slmodel = slmodel,
whichopt = whichopt, useC = useC, dmean = dmean, getHess = getHess, itmean = itmean,
indfixx = indfixx, fixx = fixx, xreg = xreg, r = r, s = s, b = b, straightRegress = straightRegress)
stopCluster(cl)
setDefaultCluster(NULL)
}
allpars <- vector("list", num)
logl <- rep(0, num)
for (i in 1:num) {
pars <- fit[[i]][[1]]
pars[indfixx] <- fixx[indfixx]
allpars[[i]]$phi <- if (p > 0)
pars[1:p] else numeric(0)
allpars[[i]]$theta <- if (q > 0)
pars[p + (1:q)] else numeric(0)
allpars[[i]]$phiseas <- if (pseas > 0)
pars[p + q + (1:pseas)] else numeric(0)
allpars[[i]]$thetaseas <- if (qseas > 0)
pars[p + q + pseas + (1:qseas)] else numeric(0)
allpars[[i]]$phip <- if (p > 0)
ARToPacf(pars[1:p]) else numeric(0)
allpars[[i]]$thetap <- if (q > 0)
ARToPacf(pars[p + (1:q)]) else numeric(0)
allpars[[i]]$phiseasp <- if (pseas > 0)
ARToPacf(pars[p + q + (1:pseas)]) else numeric(0)
allpars[[i]]$thetaseasp <- if (qseas > 0)
ARToPacf(pars[p + q + pseas + (1:qseas)]) else numeric(0)
allpars[[i]]$dfrac <- if (lmodel == "d")
pars[p + q + pseas + qseas + 1] else numeric(0)
allpars[[i]]$H <- if (lmodel == "g")
pars[p + q + pseas + qseas + 1] else numeric(0)
allpars[[i]]$alpha <- if (lmodel == "h")
pars[p + q + pseas + qseas + 1] else numeric(0)
allpars[[i]]$dfs <- if (slmodel == "d")
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$alpha) + 1] else numeric(0)
allpars[[i]]$Hs <- if (slmodel == "g")
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$alpha) + 1] else numeric(0)
allpars[[i]]$alphas <- if (slmodel == "h")
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$alpha) + 1] else numeric(0)
allpars[[i]]$delta <- if (!is.null(xreg) && sum(r) > 0)
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + 1:(sum(r))] else numeric(0)
allpars[[i]]$omega <- if (!is.null(xreg) && sum(s) > 0)
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + sum(r) + 1:sum(s)] else numeric(0)
allpars[[i]]$muHat <- if (is.logical(dmean) && dmean)
pars[p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + sum(r) + sum(s) + 1] else meanval
rr <- tacvfARFIMA(phi = allpars[[i]]$phi, theta = allpars[[i]]$theta, phiseas = allpars[[i]]$phiseas,
thetaseas = allpars[[i]]$thetaseas, dfrac = allpars[[i]]$dfrac, dfs = allpars[[i]]$dfs,
H = allpars[[i]]$H, Hs = allpars[[i]]$Hs, alpha = allpars[[i]]$alpha, alphas = allpars[[i]]$alphas,
period = period, maxlag = length(y) - 1)
if (itmean) {
allpars[[i]]$muHat <- TrenchMean(rr, y)
}
if (is.null(fit[[i]][2])) {
cat("mode ", i, " has no loglikelihood value", sep = "")
allpars[[i]]$loglik <- logl[i] <- -3e+10
} else {
allpars[[i]]$loglik <- logl[i] <- fit[[i]][[2]]
}
allpars[[i]]$pars <- c(allpars[[i]]$phi, allpars[[i]]$theta, allpars[[i]]$phiseas,
allpars[[i]]$thetaseas, allpars[[i]]$dfrac, allpars[[i]]$H, allpars[[i]]$alpha,
allpars[[i]]$dfs, allpars[[i]]$Hs, allpars[[i]]$alphas, allpars[[i]]$delta,
allpars[[i]]$omega)
allpars[[i]]$parsp <- c(allpars[[i]]$phip, allpars[[i]]$thetap, allpars[[i]]$phiseasp,
allpars[[i]]$thetaseasp, allpars[[i]]$dfrac, allpars[[i]]$H, allpars[[i]]$alpha,
allpars[[i]]$dfs, allpars[[i]]$Hs, allpars[[i]]$alphas, allpars[[i]]$delta,
allpars[[i]]$omega)
if (getHess) {
hess <- fit[[i]]$hessian[!indfixx, !indfixx]
se <- rep(-Inf, numvarsM)
se[!indfixxM] <- tryCatch(sqrt(diag(solve(-hess))),
error = function(err) rep(NA, max(dim(hess)[1], 1)))
allpars[[i]]$se <- se
allpars[[i]]$hess <- hess
} else hess <- numeric(0)
if (!is.null(xreg)) {
omega <- allpars[[i]]$omega
delta <- allpars[[i]]$delta
xreglist <- list(s = s, r = r, b = b, omega = vector("list", length(r)),
delta = vector("list", length(r)))
if (!straightRegress) {
for (j in 1:length(r)) {
if (r[j] == 0) {
xreglist$delta[[j]] <- NA
if (s[j] == 1) {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- namexreg[j]
} else {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- paste("omega(", 0:(s[j] - 1), ").",
namexreg[j], sep = "")
}
} else {
if (j == 1)
xreglist$delta[[j]] <- delta[1:r[1]] else {
rrr <- cumsum(r)[j - 1]
xreglist$delta[[j]] <- delta[(rrr + 1):(rrr + r[j])]
}
names(xreglist$delta[[j]]) <- paste("delta(", 1:r[j], ").", namexreg[j],
sep = "")
if (s[j] == 1) {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- paste("omega(0).", namexreg[j], sep = "")
} else {
if (j == 1)
xreglist$omega[[j]] <- omega[1:s[1]] else {
ss <- cumsum(s)[j - 1]
xreglist$omega[[j]] <- omega[(ss + 1):(ss + s[j])]
}
names(xreglist$omega[[j]]) <- paste("omega(", 0:(s[j] - 1), ").",
namexreg[j], sep = "")
}
}
}
}
else if (straightRegress) {
xreglist$omega <- as.list(omega)
names(xreglist$omega) <- namexreg
}
allpars[[i]]$xreglist <- xreglist
mmu <- allpars[[i]]$muHat
if(straightRegress) {
resR <- z - xr %*% omega
resRdiff <- y - xreg %*% omega - mmu
}
else {
resR <- funcTF(y = z, x = as.vector(xr), delta = delta, omega = omega, b = b,
rx = r, sx = s, nx = length(z), meanval = 0)$y ##regression residuals
resRdiff <- funcTF(y = y, x = as.vector(xreg), delta = delta, omega = omega,
b = b, rx = r, sx = s, nx = length(y), meanval = 0)$y
}
res <- DLResiduals(rr, resRdiff)
sigma2muse <- getsigma2muhat(resRdiff, phi = allpars[[i]]$phi, theta = allpars[[i]]$theta,
thetaseas = allpars[[i]]$thetaseas, phiseas = allpars[[i]]$phiseas, period = period,
dfrac = allpars[[i]]$dfrac, dfs = allpars[[i]]$dfs, H = allpars[[i]]$H,
Hs = allpars[[i]]$Hs, alpha = allpars[[i]]$alpha, alphas = allpars[[i]]$alphas,
mean.dof = mean.dof)
allpars[[i]]$fitted <- resRdiff - res
} else {
yy <- y - allpars[[i]]$muHat
resR <- NULL
sigma2muse <- getsigma2muhat(yy, phi = allpars[[i]]$phi, theta = allpars[[i]]$theta,
thetaseas = allpars[[i]]$thetaseas, phiseas = allpars[[i]]$phiseas, period = period,
dfrac = allpars[[i]]$dfrac, dfs = allpars[[i]]$dfs, H = allpars[[i]]$H,
Hs = allpars[[i]]$Hs, alpha = allpars[[i]]$alpha, alphas = allpars[[i]]$alphas,
mean.dof = mean.dof)
res <- DLResiduals(rr, yy)
allpars[[i]]$fitted <- yy - res
}
allpars[[i]]$regResiduals <- resR
allpars[[i]]$residuals <- res
allpars[[i]]$sigma2 <- sigma2muse
if (getHess && !is.null(xreg) && sum(s) > 0) {
inds <- p + q + pseas + qseas + length(allpars[[i]]$dfrac) + length(allpars[[i]]$H) +
length(allpars[[i]]$dfs) + length(allpars[[i]]$Hs) + length(allpars[[i]]$alpha) +
length(allpars[[i]]$alphas) + sum(r) + 1:(sum(s) + 1)
XX <- cbind(xreg, rep(1, length(z)))
XXi <- solve(t(XX)%*% XX)
V <- toeplitz(rr)
covTS <- XXi %*% t(XX) %*% V %*% XX %*% XXi
allpars[[i]]$se[inds] <- sqrt(diag(sigma2muse*covTS))
}
class(allpars[[i]]) <- "ARFIMA"
}
ord <- order(logl, decreasing = TRUE)
allpars <- allpars[ord]
index <- 1:num
weeded <- (num == 1)
ans <- list(z = z, differencing = differencing, dint = dint, dseas = dseas, period = period,
modes = allpars, dmean = dmean, itmean = itmean, p = p, q = q, pseas = pseas,
qseas = qseas, lmodel = lmodel, slmodel = slmodel, weeded = weeded, getHess = getHess,
numvars = numvars, numcut = 0, n = length(z), xreg = xr, r = r, s = s, b = b,
call = match.call(), flag0 = flag0, numeach = numeach, strReg = straightRegress,
regOnly = regOnly, intindex = intindex, intname = intname, namexreg=namexreg,
fixx = fixx, indfixx = indfixx)
}
class(ans) <- "arfima"
if (autoweed && !ans$weeded) {
ans <- weed(ans, type = weedtype[1], eps2 = weedeps, adapt = adapt, pn = weedp)
}
ans
}
# Tells whether points share a circle in the p-norm
pcircle <- function(x1, x2, rad, pn = 2) {
zz <- x1 - x2
z <- (sum(abs(zz)^pn))^(1/pn)
return(z > rad)
}
# Weeds out similar (or the same) modes
#' Weeds out fits from a call to arfima that are too close to each other.
#'
#' Weeds out fits from a call to arfima that are too close to each other.
#'
#'
#' @param ans The result from a call to arfima.
#' @param type The space to perform the weeding in. "A" is for operating
#' parameters. "P" is in the PACF space. "B" performs weeding in both. "N"
#' performs no weeding and is only used internally.
#' @param walls If more than one mode is on a wall in the PACF space, all modes
#' but the one with the highest log-likelihood on the same wall are deleted.
#' @param eps2 The maximum distance between modes that are close together for
#' the mode with the lower log-likelihood to be weeded out. If adapt is TRUE
#' (default) this value changes.
#' @param eps3 The minimum distance from a wall for a secondary mode to be
#' weeded out, if walls are TRUE.
#' @param adapt If TRUE, if dim is the dimensionality of the search, eps2 is
#' changed to \eqn{(1 + eps2)^{dim} - 1}.
#' @param pn The p in the p-norm to be used in the weeding. p = 2 (default) is
#' Euclidean distance.
#' @return An object of class "arfima" with modes possibly weeded out.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}, \code{\link{distance}}
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(1234)
#' sim <- arfima.sim(1000, model = list(theta = 0.9, dfrac = 0.4))
#' fit <- arfima(sim, order = c(0, 0, 1), autoweed = FALSE, back=TRUE)
#' fit
#' distance(fit)
#' fit1 <- weed(fit)
#' fit1
#' distance(fit1)
#'}
#' @export weed
weed <- function(ans, type = c("A", "P", "B", "N"), walls = FALSE, eps2 = 0.025, eps3 = 0.01,
adapt = TRUE, pn = 2) {
if (!inherits(ans, "arfima"))
stop("weed only defined for arfima objects")
allpars <- ans$modes
if(length(allpars)==1) {
ans$weeded <- TRUE
return(ans)
}
if (ans$flag0){
stopifnot(length(allpars[[1]]$pars)==0)
allpars <- allpars[[1]]
ans$modes <- allpars
ans$weeded <- TRUE
return(ans)
}
if (length(pn) == 0 || pn <= 0 || !is.finite(pn)) {
warning("invalid pn in weed: setting to 2")
pn <- 2
}
if (length(eps2) == 0 || eps2 < 0) {
warning("invalid epsilon in weed: setting to 0.025")
eps2 <- 0.025
}
type <- type[1]
if (!(type %in% c("A", "P", "B", "N"))) {
warning("invalid type in weed: setting to \"A\"")
type <- "A"
}
if (length(eps3) == 0 || eps3 <= 0) {
warning("invalid wall epsilon in weed: setting to 0.01")
eps3 <- 0.01
}
if (!is.logical(walls)) {
warning("invalid walls selection in weed: setting to FALSE")
walls <- FALSE
}
modes <- vector("list", 1)
p <- ans$p
q <- ans$q
pseas <- ans$pseas
qseas <- ans$qseas
lmodel <- ans$lmodel
slmodel <- ans$slmodel
numvars <- p + q + qseas + pseas + if (lmodel != "n")
1 else 0
numvars <- numvars + if (slmodel != "n")
1 else 0
numvars <- numvars + sum(ans$r) + sum(ans$s)
if(numvars==0) {
stop("Should not get here. Please contact the package maintainer.")
}
if (adapt)
eps2 <- (1 + eps2)^numvars - 1
num <- length(allpars)
logl <- rep(0, num)
for (i in 1:num) logl[i] <- allpars[[i]]$loglik
ord <- order(logl, decreasing = TRUE)
allpars <- allpars[ord]
if (type %in% c("A", "B")) {
modes <- vector("list", 1)
modes[[1]] <- allpars[[1]]
flag <- NULL
if (num > 1)
for (i in 2:num) {
flag <- NULL
for (j in 1:length(modes)) {
flag <- c(flag, pcircle(x1 = allpars[[i]]$pars, x2 = modes[[j]]$pars,
rad = eps2, pn = pn))
}
if (all(flag)) {
modes[[length(modes) + 1]] <- allpars[[i]]
}
}
} else modes <- allpars
num <- length(modes)
logl <- rep(0, num)
for (i in 1:num) logl[i] <- modes[[i]]$loglik
ord <- order(logl, decreasing = TRUE)
modes <- modes[ord]
if (type %in% c("P", "B")) {
modes1 <- modes
modes <- vector("list", 1)
modes[[1]] <- modes1[[1]]
flag <- NULL
if (!ans$regOnly && num > 1) {
for (i in 2:num) {
flag <- NULL
for (j in 1:length(modes)) {
flag <- c(flag, pcircle(x1 = modes1[[i]]$parsp, x2 = modes[[j]]$parsp,
rad = eps2, pn = pn))
}
if (all(flag)) {
modes[[length(modes) + 1]] <- modes1[[i]]
}
}
}
else if(num>1) {
for (i in 2:num) {
flag <- NULL
for (j in 1:length(modes)) {
flag <- c(flag, pcircle(x1 = modes1[[i]]$omega, x2 = modes[[j]]$omega,
rad = eps2, pn = pn))
}
if (all(flag)) {
modes[[length(modes) + 1]] <- modes1[[i]]
}
}
}
}
wall <- inds <- vector("list", numvars * 2)
pl1 <- p2 <- 0
arma_varnum <- p + q + pseas + qseas
for (i in 1:numvars) {
if(ans$regOnly) break
if (i <= arma_varnum) {
minner = -1
maxxer = 1
} else if (lmodel != "n" && !pl1) {
pl1 <- 1
if (lmodel == "d") {
minner = -1
maxxer = 0.5
} else if (lmodel == "g") {
minner = 0
maxxer = 1
} else if (lmodel == "h") {
minner = 0
maxxer = 3
}
} else {
if (slmodel == "d") {
minner = -1
maxxer = 0.5
} else if (slmodel == "g") {
minner = 0
maxxer = 1
} else if (slmodel == "h") {
minner = 0
maxxer = 3
}
}
for (j in 1:length(modes)) {
wall[[2 * i - 1]] <- c(wall[[2 * i - 1]], abs(maxxer - modes[[j]]$parsp[i]) <
eps3)
wall[[2 * i]] <- c(wall[[2 * i]], abs(minner - modes[[j]]$parsp[i]) < eps3)
}
inds[[2 * i - 1]] <- which(wall[[2 * i - 1]])
inds[[2 * i]] <- which(wall[[2 * i]])
}
index <- numeric(0)
for (i in 1:(2 * numvars)) index <- c(index, inds[[i]])
index <- unique(index)
numcut <- length(ans$index) - length(index)
ans1 <- list(z = ans$z, differencing = ans$differencing, dint = ans$dint, dseas = ans$dseas,
period = ans$period, modes = modes, dmean = ans$dmean, itmean = ans$itmean, getHess = ans$getHess,
p = p, q = q, pseas = pseas, qseas = qseas, lmodel = lmodel, slmodel = slmodel,
numcut = numcut, index = index, weeded = TRUE, n = ans$n, xreg = ans$xreg, r = ans$r,
s = ans$s, b = ans$b, call = ans$call, flag0 = ans$flag0, numeach = ans$numeach,
strReg = ans$strReg, namexreg = ans$namexreg, regOnly = ans$regOnly, intindex = ans$intindex,
intname = ans$intname, fixx = ans$fixx, indfixx = ans$indfixx)
class(ans1) <- "arfima"
ans1
}
# computes and reports the distance between all modes
#' The distance between modes of an \code{arfima} fit.
#'
#' The distance between modes of an \code{arfima} fit.
#'
#'
#' @param ans An object of class "arfima".
#' @param p The p in the p-norm to be used.
#' @param digits The number of digits to print.
#' @return A list of two data frames: one with distances in operator space, the
#' second with distances in the transformed (PACF) space.
#' @author JQ (Justin) Veensta
#' @references Veenstra, J.Q. Persistence and Antipersistence: Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(8564)
#' sim <- arfima.sim(1000, model = list(phi = c(0.2, 0.1), dfrac = 0.4, theta = 0.9))
#' fit <- arfima(sim, order = c(2, 0, 1), back=TRUE)
#'
#' fit
#'
#' distance(fit)
#'}
#' @export distance
distance <- function(ans, p = 2, digits = 4) {
if (!inherits(ans, "arfima"))
stop("distance only defined for arfima objects")
if (p <= 0)
stop("invalid p")
num <- length(ans$modes)
distances1 <- distances2 <- matrix(0, nrow = num, ncol = num)
if (num == 1)
stop("only one mode")
for (i in 1:(num)) {
for (j in 1:num) {
vals <- ans$modes[[i]]$pars - ans$modes[[j]]$pars
distances1[i, j] <- (sum(abs(vals)^p))^(1/p)
vals <- ans$modes[[i]]$parsp - ans$modes[[j]]$parsp
distances2[i, j] <- (sum(abs(vals)^p))^(1/p)
}
}
distances1 <- round(distances1, digits = digits)
distances2 <- round(distances2, digits = digits)
rownames(distances1) <- colnames(distances1) <- rownames(distances2) <- colnames(distances2) <- paste("Mode",
1:num)
distances <- vector("list", 2)
distances[[1]] <- distances1
distances[[2]] <- distances2
names(distances) <- c(paste(p, "-norm distance", sep = ""), paste(p, "-norm distance in transformed space",
sep = ""))
distances
}
#' Finds the best modes of an \code{arfima} fit.
#'
#' Finds the best modes of an \code{arfima} fit with respect to log-likelihood.
#'
#' This is the easiest way to remove modes with lower log-likelihoods.
#'
#' @param object An object of class "arfima".
#' @param bestn The top number of modes to keep with respect to the
#' log-likelihood.
#' @return The \code{bestn} "best" modes.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(8765)
#' sim <- arfima.sim(1000, model = list(phi = 0.4, theta = 0.9, dfrac = 0.4))
#' fit <- arfima(sim, order = c(1, 0, 1), back=TRUE)
#' fit
#' fit <- bestModes(fit, 2)
#' fit
#'}
#' @export bestModes
bestModes <- function(object, bestn) {
if (!inherits(object, "arfima"))
stop("bestModes only defined for arfima objects")
if (!object$weeded)
stop("Please weed the object first before calling bestModes.")
m <- length(object$modes)
if (bestn > m)
stop("bestn is larger than the number of modes")
modes <- object$modes
logls <- sapply(1:m, function(i) modes[[i]]$logl)
ord <- order(logls, decreasing = TRUE)
ord <- ord[1:bestn]
modes <- modes[ord]
object$modes <- modes
object <- weed(object, eps2 = 0)
object
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.