## Do not edit this file manually.
## It has been automatically generated from *.org sources.
## todo: expand_ar() seems unused, maybe it is obsolete(d)?
expand_ar <- function(ar, sar, s){
x <- polynom(c(0,1))
phi <- polynom(c(1, -ar))
sphi <- polynom(c(1, -sar))(x^s)
coef(phi * sphi)
}
poly_x <- polynom(c(0,1))
poly_delta <- polynom(c(1, -1)) # = 1 - poly_x
.sarima_env <- function(){
.process_fixed <- function(x, fixed = NULL){
if(is.null(x))
stop("'x' must have a value here")
xlen <- length(x)
## (2018-08-23) TODO: this if/else sequence needs cleaning up!
if(is.null(fixed))
rep(FALSE, xlen)
else if(isTRUE(fixed))
rep(TRUE, xlen)
else if(is.logical(fixed)){
## 2018-08-23 fiixng missing case length=1 and FALSE
## TODO: check nevertheless if this was omitted on purpose?
if(length(fixed) == 1)
rep(fixed, xlen)
else if(length(fixed) == xlen)
fixed
else
stop("if 'fixed' is logical, its length must be equal to one or the order.")
}else if(is.numeric(fixed)){
## indices of fixed elements, must be whole numbers
if(!is.integer(fixed) && any(round(fixed) != fixed))
stop("if numeric, 'fixed' must contain integer numbers")
if(any(fixed <= 0) || any(fixed > xlen))
stop("if numeric, 'fixed' must contain integers between 1 and 'order'")
res <- rep(FALSE, xlen)
res[fixed] <- TRUE
res # needed since 'res[fixed] <- value' returns the rhs
}else{
stop("unsupported type of 'fixed'")
}
}
## currently essentially identical to .process_fixed()
.process_atanh.tr <- function(x, transf = NULL){
if(is.null(x))
stop("'x' must have a value here")
xlen <- length(x)
if(is.null(transf))
rep(FALSE, xlen)
else if(isTRUE(transf))
rep(TRUE, xlen)
else if(is.logical(transf)){
if(length(transf) == 1)
rep(transf, xlen)
else if(length(transf) == xlen)
transf
else
stop("if 'transf' is logical, its length must be equal to one or the order.")
}else if(is.numeric(transf)){
## indices of transf elements, must be whole numbers
if(!is.integer(transf) && any(round(transf) != transf))
stop("if numeric, 'transf' must contain integer numbers")
if(any(transf <= 0) || any(transf > xlen))
stop("if numeric, 'transf' must contain integers between 1 and 'order'")
res <- rep(FALSE, xlen)
res[transf] <- TRUE
res
}else{
stop("unsupported type of 'transf'")
}
}
## parcor
## roots ?
## fixed
## nonfixed
## convention ?
## dispname ?
##
## based on co_ar() in sarima() and ar() in .sarima_env
.set_coef <- function(order, coef = NULL, fixed = NULL, nonfixed = NULL,
atanh.tr = NULL, ...){
wrk <- NULL
coef.type <- "coef" # todo: any subtleties?
if(!is.null(order))
wrk <- rep(NA_real_, order)
if(!missing(coef) && !is.null(coef)){
if(is.null(wrk))
wrk <- coef
else{
# todo: is it reasonable to allow a smaller length, so that coef's from
# lower order fits can be used as initial values more easily?
## Note (JH, 11/10/18): I like that length(coef) == length(wrk)
## but understand it is not necessary.
## If this is true, 'p' is not required (see note below)
if(length(coef) != length(wrk))
stop("if given, 'coef' must have length equal to 'order'")
## 2018-08-23 rewriting using seq_along:
wrk[seq_along(coef)] <- coef # wrk[1:length(coef)] <- coef
}
}
a <- .process_fixed(wrk, fixed)
b <- .process_atanh.tr(wrk, atanh.tr)
if(!missing(nonfixed) || !is.null(nonfixed)){
# declare non-fixed some coef's previously marked as fixed
## !!! 2018-08-23 - this 'b' destroys the 'b' from atanh.tr, seems oversight!
## Changing this 'b' to 'bb'.
## TODO: verify nevertheless!
## b <- .process_fixed(wrk, nonfixed)
bb <- .process_fixed(wrk, nonfixed)
# !bb is FALSE only for the elements declared nonfixed by 'nonfixed', so this
# ! will transfer this to 'a', leaving others as they are (fixed or nonfixed)
a <- a & !bb
}
res <- wrk
attr(res, "fixed") <- a
attr(res, "atanh.tr") <- b
## store addtional info, e.g. parcor, negate
dots <- list(...)
for(at in names(dots))
attr(res, at) <- dots[[at]]
res
}
ar <- function(p = 0, ar, sign = "-", atanh.tr = TRUE, ...){
## NOTE: do not name p and ar in the call for flexibility:
coef <- .set_coef(p, ar, sign = sign, dispname = "ar", atanh.tr = atanh.tr, ...)
list(name = "ar", p = p, coef = coef) # TODO: remove 'p' from the return value?
## Note (JH, 11/10/18): 'p' is never used, I would agree on removing 'p'.
## If one needed, p could be obtained from the length of coef
}
ma <- function(q = 0, ma, sign = "+", atanh.tr = TRUE, ...){
coef <- .set_coef(q, ma, sign = sign, dispname = "ma", atanh.tr = atanh.tr, ...)
list(name = "ma", q = q, coef = coef)
}
s <- function(...){
## TODO: for now 'm' of length 1
## Note (JH, 11/10/18): what else would be an argument here?
dots <- c(...)
if(length(dots) != 1)
stop("currently only one argument is supported")
p <- dots[1] - 1
coef <- .set_coef(p, rep(1, p), sign = "+", fixed = TRUE, operator = TRUE, dispname = "s")
list(name = "s", p = p, coef = coef)
}
## unit roots
i <- function(d, ...){
## TODO: this needs special method later!
coef <- .set_coef(1, 1, sign = "-", fixed = TRUE, d = d, operator = TRUE,
dispname = "i", ...)
list(name = "i", d = d, p = 1, coef = coef)
}
u <- function(u, fixed = TRUE, operator = all(fixed), ...){
# u may be a list to accommodate mixed complex and real values, probably
# not a good idea. Do not convert to use Recall() if u is 'list'
f <- function(x){
if(is.complex(x))
-2*cos(Arg(x)) # equivalently: -2*Re(x), but this assumes that Mod(x)=1
else -2 * cospi(2*x) # x is a fraction of 2pi
}
co <- sapply(u, f)
coef <- lapply(co, function(x) .set_coef(2, c(x, 1), sign = "+", fixed = fixed,
operator = operator, dispname = "su") )
list(name = "u", u = u, coef = coef)
}
## spec for unit roots to be estimated (unless fixed)
## TODO: require p >= 2 in "uar"?
## Note (JH 11/10/18): We should require this.
uar <- function(p = 2, parcor, sign = "-", atanh.tr = TRUE, fixed = NULL, ...){
## Note (JH 11/10/18): adding checks
## check p is at least 2
if(p < 2)
stop("'p' must be greater than 2 for unit AR polynomials")
## the last value should be +/-1
if(abs(parcor[p]) != 1)
stop("the last coefficient supplied should be +1 or -1")
if(is.null(fixed)){
fixed <- rep(FALSE, p)
if(p > 0)
fixed[p] <- TRUE
}
# NOTE: do not name p and parcor in the call for flexibility
coef <- .set_coef(p, parcor, sign = sign, dispname = "uar", atanh.tr = atanh.tr,
fixed = fixed, ...)
list(name = "uar", p = p, coef = coef) # TODO: rename 'coef' to 'parcor' here?
}
## seasonal
sar <- function(s, p = 0, ar, sign = "-", atanh.tr = TRUE, ...){
# NOTE: do not name p and ar in the call for flexibility
coef <- .set_coef(p, ar, sign = sign, nseasons = s, dispname = "sar",
atanh.tr = atanh.tr, ...)
list(name = "sar", s = s, p = p, coef = coef) # TODO: remove 's' and 'p'?
}
## Note (JH, 11/10/18): 's' and 'p' aren't used, so removing them won't be a problem.
sma <- function(s, q = 0, ma, sign = "+", atanh.tr = TRUE, ...){
coef <- .set_coef(q, ma, sign = sign, nseasons = s, dispname = "sma",
atanh.tr = atanh.tr, ...)
list(name = "sma", s = s, q = q, coef = coef)
}
ss <- function(s, ...){
## TODO: for now 'm' of length 1
dots <- c(...)
if(length(dots) != 1)
stop("currently only one argument is supported")
p <- dots[1] - 1
coef <- .set_coef(p, rep(1, p), sign = "+", nseasons = s, fixed = TRUE,
operator = TRUE, dispname = "ss")
list(name = "ss", s = s, p = p, coef = coef)
}
## seasonal unit roots
si <- function(s, d, ...){
coef <- .set_coef(1, 1, sign = "-", nseasons = s, fixed = TRUE,
operator = TRUE, dispname = "si", ...)
list(name = "si", s = s, d = d, coef = coef)
}
su <- function(s, h){
# TODO: check that elements of h are among 1,2, ..., s/2 (excluding s/s
# if s is even) coef is a list here!
## (JH, 11/10/18): Checking h
int_cond <- any(h %% 1 > 0)
if(int_cond)
stop("'h' must contain only (positive) integer values")
pos_cond <- any(h < 1)
if(int_cond)
stop("'h' must contain only positive (integer) values")
max_cond <- !all(h < s/2)
if(max_cond)
stop("'h' must be a positive integer less than 's/2'")
coef <- lapply(h, function(x)
.set_coef(2, c(- 2 * cospi(2*x/s), 1),
sign = "+", fixed = TRUE, operator = TRUE,
dispname = "su", su.nseasons = s, su.harmonic = x) )
list(name = "su", s = s, u = h, coef = coef)
}
.specials <- ls() # assumes that user facing objects do not start with dot!
.sarima_descr <- list()
class(.sarima_descr) <- "sarimadescr"
environment()
}# end of .sarima_env()
SARIMA <- function(formula){
e <- .sarima_env()
parent.env(e) <- environment(formula)
# keep.order keeps the order for the data frame.
# The order in attribute "specials" is not affected by this.
te <- terms(formula, specials = e$.specials, keep.order = TRUE)
# tmp <- (attr(te, "variables")[[3]]); eval(tmp, envir = e)
termvars <- attr(te, "variables")
sp <- attr(te, "specials")
res <- list()
indices <- lapply(names(sp),
function(nam){
ind <- sp[[nam]]
if(is.null(ind))
return(integer(0))
names(ind) <- rep(nam, length(ind))
ind
}
)
indices <- sort(unlist(indices))
nams <- names(indices)
indices <- indices + 1
res <- vector("list", length = length(indices))
names(res) <- nams
for(i in seq_along(indices)){
## evaluate the terms (like ar(), ma(), etc.)
res[[i]] <- eval(termvars[[indices[i]]], envir = e)
}
res
}
trendMaker <- function(formula, data = NULL, time = NULL){
if(is.null(data))
data.env <- new.env(hash = FALSE, parent = environment(formula))
else{
stopifnot(is.environment(data))
data.env <- data
}
environment(formula) <- data.env
data.env$formula <- formula
# data.env$.t.orig can be used to detect if 't' is not the original one
data.env$.t.orig <- data.env$t <- time
.t.orig <- "dummy for R CMD check; .t.orig is needed in dataenv, as set above"
data.env$.p <- function(degree){
if(identical(t, .t.orig))
res <- poly(t, degree = degree, simple = TRUE)
else{
wrk <- poly(.t.orig, degree = degree, simple = FALSE)
res <- predict(wrk, newdata = t)
}
colnames(res) <- paste0("ortht", seq_len(degree))
res
}
environment(data.env$.p) <- data.env
# cs <- function(s, freq){
# ## s - scalar; k in 1,2,...,[s/2]
# wrk <- lapply(freq, function(k) cbind(cospi(2*time*k/s), sinpi(2*time*k/s)))
# res <- do.call("cbind", wrk)
#
# # colnames(res) <-
# # unlist(lapply(freq, function(x) paste0(c("cos", "sin"), x, ".", s)))
# res
# }
data.env$.cs <- function(s, k){
## s - scalar; k in 1,2,...,[s/2]
## TODO: for now assume scalar k too:
res <- if(length(k) == 1){ # cbind(cospi(2*t*k/s), sinpi(2*t*k/s))
structure(c(cospi(2*t*k/s), sinpi(2*t*k/s)),
dim = c(length(t), 2),
.Dimnames = list(NULL, paste0(c("c", "s"), s, ".", k) ))
}else if(length(k) > 1){
wrk <- lapply(k, function(x) cbind(cospi(2*t*x/s), sinpi(2*t*x/s)) )
wrk <- do.call("cbind", wrk)
colnames(wrk) <- paste0(c("c", "s"), s, ".", rep(k, each = 2) )
wrk
}else
stop("'k' must have positive length")
# TODO: put names, keeping them short (Hm... 7*48 ...)
# TODO: colnames(res) <- paste0(c("cos", "sin"), k, ".", s)
res
}
environment(data.env$.cs) <- data.env
## TODO: test thoroughly
data.env$.B <- function(x, lags){
flags <- t > 0
t.pos <- t[flags]
t.neg <- t[!flags]
if(any(flags)){
## 15/01/2018 - krapka to accept matrix x
if(is.matrix(x)){ ## TODO: maybe check also if ncol(1) ?
res <- matrix(NA_real_, nrow = length(t), ncol = ncol(x) * length(lags))
nc.x <- ncol(x)
curcols <- seq_len(nc.x)
for(i in seq_along(lags)){
lag <- lags[i]
res[t - lag > 0, curcols] <- x[pmax(t - lag, 0), ]
curcols <- curcols + nc.x
}
## TODO: don't know why there are no names in this case.
xnam <- deparse(substitute(x))
## TODO: check for "[" or "$" in xnam and adjust the names accordingly
colnames(res) <- paste0("L(", xnam, rep(seq_len(ncol(x)), length(lags)),
", ", rep(lags, each = ncol(x)), ")")
}else{
res <- matrix(NA_real_, nrow = length(t), ncol = length(lags))
for(i in seq_along(lags)){
lag <- lags[i]
res[t - lag > 0, i] <- x[pmax(t - lag, 0)]
}
}
}
res
}
environment(data.env$.B) <- data.env
res <- function(index){
if(missing(index)){# TODO: napravi tova s na.action to avoid losing NA's
# 15/01/2018 was: model.matrix(formula)
# changing to use model.frame since we need to keep rows with NA's
# and it seems that model.matrix() doesn't have this argument.
# TODO: !!! However, model.frame doesn't contain the intercept,
# if present! So, it gets lost here.
#
# OK, sorted, this is a reply by Brian Ripley to a similar question
# (TODO: take note also of the user defined function for na.action.):
#
# Yes, but it uses options("na.action"), so you could set that.
# Better, call model.frame then model.matrix, something like
#
# stratmat <- model.matrix(myformula,
# model.frame(myformula,mydata, na.action = function(x) x))
#
# ?model.matrix does tell you the second argument should be the result
# of model.frame, which is a pretty strong hint.
res <- model.matrix(formula, model.frame(formula, na.action = NULL))
}else{
t.old <- t
t <<- index
## TODO: krapka!!!! - drop the lhs since the length of t is now diferent.
## TODO: this will do if there is only time trend but exogeneous vars need
## further work also, will this work if there is only intercept? Maybe
## should put a dummy variable on the left?
formula <- formula(Formula::Formula(formula), lhs = FALSE, rhs = 1)
txt <- paste0(".dummy", paste0(as.character(formula), collapse = ""))
formula <- as.formula(txt, env = environment(formula))
# 15/01/2018 - changing to use model.frame since we need to keep rows with
# NA's and it seems that model.matrix() doesn't have this argument.
# (see also Ripley's remark in the 'if' clause above)
# res <- model.matrix(formula, data = data.frame(.dummy = t))
# res <- as.matrix(res[ , -1])
res <- model.frame(formula, data = data.frame(.dummy = t), na.action = NULL)
res <- model.matrix(formula, res)
t <<- t.old
res
}
}
environment(res) <- data.env
res
}
.cat_u_delta <- function(x){
f <- function(poly){
s <- environment(poly)$nseasons
if(!is.null(s))
poly <- poly(poly_x^s)
if(all(coef(poly) == 1) && (deg <- length(coef(poly)) - 1) > 3)
res <- paste0("1 + B + ... + B^", deg)
else
res <- as.character(poly, "B")
d <- environment(poly)$d
if(is.null(d) || d == 1)
paste0("(", res, ")")
else
paste0("(", res, ")^", d)
}
paste0(sapply(x, f), collapse = "")
}
.Sarima_fixed <- function(x) x$internal$fixed
.capture_fo <- function(x) if(is.numeric(x)) x else capture.output(print(x, FALSE))
.Fo.xreg <- function(x) x$internal$Fo.xreg
.Fo.sarima <- function(x) x$internal$Fo.sarima
.Fo.regx <- function(x) if(is.null(res <- x$internal$Fo.regx)) 0 else res
print.Sarima <-
function (x, digits = max(3L, getOption("digits") - 3L), se = TRUE, ...){
cat("*Sarima model*")
cat("\nCall:", deparse(x$call, width.cutoff = 75L), "", sep = "\n")
# cat("\nTrend: \n\t")
# print(Fo.trend, FALSE)
# cat("\nSARIMA specification: \n\t")
# print(Fo.sarima, FALSE)
delta <- x$internal$delta
if(length(delta) > 0){
cat("Unit root terms:\n", " ")
cat(.cat_u_delta(x$internal$delta_poly), sep = "")
cat("\n\n")
}
if (length(x$coef)) {
cat("Coefficients:\n")
coef <- round(x$coef, digits = digits)
## use NROW as if all coefs are fixed there are no var.coef's
if (se && NROW(x$var.coef)) {
ses <- rep.int(0, length(coef))
ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
coef <- matrix(coef, 1L, dimnames = list(NULL, names(coef)))
coef <- rbind(coef, s.e. = ses)
}
print.default(coef, print.gap = 2)
}
cm <- x$call$method
if(is.null(cm) || cm != "css")
cat("\nsigma^2 estimated as ", format(x$sigma2, digits = digits),
": log likelihood = ", format(round(x$loglik, 2L)),
", aic = ", format(round(x$aic, 2L)), "\n", sep = "")
else
cat("\nsigma^2 estimated as ",
format(x$sigma2, digits = digits),
": part log likelihood = ", format(round(x$loglik,2)),
"\n", sep = "")
invisible(x)
}
summary.Sarima <- function(object, ...){
## follow the 'lm' tradition
cat("\nCall:\n", paste(deparse(object$call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
## TODO: these need more stable implementation, probably should be in the object.
# object$call[[2]][[3]][[1]] is '|'
# Fo.xreg <- object$formula[[3]][[2]] # or object$call[[2]][[3]][[2]]
# Fo.sarima <- object$formula[[3]][[3]]
# Fo.regx <- if(length(object$formula[[3]]) > 3) object$formula[[3]][[4]] else 0
Fo.xreg <- .Fo.xreg (object)
Fo.sarima <- .Fo.sarima(object)
Fo.regx <- .Fo.regx (object)
xfo <- .capture_fo(Fo.xreg)
fosarima <- .capture_fo(Fo.sarima)
fox <- .capture_fo(Fo.regx)
cat("Model: Y_t - xreg_t is SarimaX\n")
cat(" xreg: ", xfo, "\n")
cat(" sarima: ", fosarima, "\n")
cat(" regx: ", fox, "\n")
cat("\n")
delta <- object$internal$delta
if(length(delta) > 0){
cat("Unit root terms:\n", " ")
cat(.cat_u_delta(object$internal$delta_poly), sep = "")
cat("\n\n")
}
fixed <- .Sarima_fixed(object)
est <- object$coef[!fixed]
se <- sqrt(diag(object$var.coef))[!fixed]
zval <- est / se
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
coefs <- cbind(est, se, zval, pval)
dimnames(coefs) <-
list(names(est),
c("Estimate", "Std. Error", "Z value", "Pr(>|z|)") )
cat("Coefficients:\n")
printCoefmat(coefs
# , digits = digits, signif.stars = signif.stars, na.print = "NA", ...
)
cat("\nestimated sigma^2 = ", format(object$sigma2, digits = 3),
", log-likelihood = ", format(round(object$loglik, 2)),
", aic = ", format(round(object$aic, 2L)),
sep = "")
#browser()
cat("\n")
invisible(object)
}
.tsdiag_choices <- c(
## "classic (std. residuals, acf, portmanteau p-values)",
"residuals",
"acf of residuals",
"p values for Ljung-Box statistic",
"p values for Li-McLeod statistic",
"p values for Box-Pierce statistic",
"pacf of residuals"
)
tsdiag.Sarima <- function(object, gof.lag = NULL, ask = FALSE, ..., plot = 1:3, layout = NULL)
{
if(is.null(gof.lag))
gof.lag <- 20 # NOTE: arbitrary value
else if(!is.numeric(gof.lag))
stop("'gof.lag' must be numeric and contain positive integers")
lag.max <- max(gof.lag)
## plot standardized residuals, acf of residuals, Ljung-Box p-values
err <- object$residuals
stdres <- err / sqrt(object$sigma2)
choices <- .tsdiag_choices
chnum <- 1:length(choices)
if(!isTRUE(plot)){ # plot is typically numeric index here;
choices <- choices[plot] # FALSE or NULL give zero length result, so no plots
chnum <- chnum[plot]
if(anyNA(choices)){
warning("'plot' should be TRUE/FALSE or vector of positive integers <= ",
length(.tsdiag_choices), ",\n", "ignoring non-existent values")
chnum <- chnum[!is.na(choices)]
choices <- choices[!is.na(choices)]
}
}
if(length(choices) > 0){
old.par <- par(no.readonly = TRUE)
on.exit(par(old.par)) # restore graphics parameters before exiting.
# n_per_page <- 3
n_per_page <- if(is.null(layout))
layout(matrix(1:min(3, length(choices)), ncol = 1))
else
## TODO: this needs further thought!
do.call("layout", layout) # for the time being
ask_user <- interactive() && (ask || length(choices) > n_per_page)
choice_title <- "Select a plot number or 0 to exit"
ch_index <- if(length(choices) == 1)
1
else if(ask)
menu(choices, title = choice_title)
else if(!identical(plot, FALSE))
1
else
integer(0)
choice <- chnum[ch_index]
## precompute common stuff for portmanteau tests
nlag <- gof.lag
pval <- numeric(nlag)
fitdf <- if(inherits(object, "Sarima"))
length(object$internal$nonfixed)
else if(inherits(object, "Arima"))
sum(object$arma[1:4]) # object$arma is: p, q, p_s. q_s, s, d, d_s
else
0
# for(i in 1L:nlag)
# pval[i] <- Box.test(err, i, type="Ljung-Box",
# fitdf = ifelse(i > fitdf, fitdf, i - 1))$p.value
sacf <- autocorrelations(err, maxlag = nlag) # deal with NA's?
res <- list(residuals = err,
"LjungBox" = NULL,
"LiMcLeod" = NULL,
"BoxPierce" = NULL)
while(length(choice) != 0){
switch(choice,
{ # 1: "residuals",
plot(stdres, type = "h", main = "Standardized Residuals", ylab = "")
abline(h = 0)
#acf(err, main = "ACF of residuals from model", lag.max = lag.max)
#pacf(err, main = "PACF of residuals from model", lag.max = lag.max)
},
{ # 2: "ACF of residuals"
## acf
acf(err, plot = TRUE, main = "ACF of Residuals", lag.max = lag.max, na.action = na.pass)
# acf(cdf, main = "", lag.max = lag.max)
# title("ACF of" ~U[t])
# pacf(cdf, main = "", lag.max = lag.max)
# title("PACF of" ~U[t])
},
{ # 3: "Ljung-Box p-values"
acftest <- acfIidTest(sacf, npar = fitdf, nlags = 1:nlag, method = "LjungBox",
interval = NULL)
res[["LjungBox"]] <- acftest
},
{ # 4: "Li-McLeod p-values"
acftest <- acfIidTest(sacf, npar = fitdf, nlags = 1:nlag, method = "LiMcLeod",
interval = NULL)
res[["LiMcLeod"]] <- acftest
},
{ # 5: "Box-Pierce p-values"
acftest <- acfIidTest(sacf, npar = fitdf, nlags = 1:nlag, method = "BoxPierce",
interval = NULL)
res[["BoxPierce"]] <- acftest
},
{ # 6: "PACF of residuals"
## acf
pacf(err, plot = TRUE, main = "PACF of Residuals", lag.max = lag.max, na.action = na.pass)
},
# { # 4: "ACF/Histogram of tau_residuals"
# acf(err2, main = "ACF of tau_residuals", lag.max = lag.max)
# hist(err2, freq = FALSE, main = "Histogram of tau_residuals", xlab = "",
# ylim = c(0, 0.5))
# lines(seq(-5, 5, .01), dnorm(seq(-5, 5, .01)), col = "red")
# }
)
if(choice %in% 3:5){ # common code for portmanteau tests
pval <- acftest$test[ , "pvalue"]
plot(1L:nlag, pval, xlab = "lag", ylab = "p value", ylim = c(0,1),
main = .tsdiag_choices[choice])
abline(h = 0.05, lty = 2, col = "blue")
}
if(length(chnum) == 1) # length(choices) == 1
break
if(ask_user) {
ch_index <- menu(choices, title = choice_title)
choice <- chnum[ch_index]
} else{
## just plot the next one
## Note: this doesn't update ch_index
chnum <- chnum[-1]
choice <- chnum[1]
}
}
}
class(res) <- "tsdiagSarima"
invisible(res)
}
## fkf: y ~ d x 1
## state ~ m x 1
armapqss <- function(ar, ma, sigma) {
p <- length(ar)
q <- length(ma)
r <- max(p, q + 1) # fkf: m
ear <- c(ar, numeric(max(r - p, 0)))
ema <- c(ma, numeric(max(r - q - 1, 0)))
Tt <- cbind(ear, rbind(diag(r - 1), 0))
Rt <- matrix(c(1, ema), ncol = 1)
Zt <- matrix(c(1, numeric(r - 1)), nrow = 1)
ct <- matrix(0) # obs. regression
dt <- matrix(0, nrow = r, ncol = 1) # state regression
GGt <- matrix(0) # d x d
H <- Rt * sigma
HHt <- H %*% t(H)
a0 <- numeric(r)
P0 <- diag(1e6, nrow = r)
list(a0 = a0, P0 = P0, ct = ct, dt = dt,
Zt = Zt, Tt = Tt, GGt = GGt, HHt = HHt)
}
## extention of armapqss
## fkf: y ~ d x 1
## state ~ m x 1
xarmaxss <- function(ar, ma, sigma, xreg, regx) {
d <- 1 # dim of y, currently univariate
p <- length(ar)
q <- length(ma)
r <- max(p, q + 1) # fkf: m = r
ear <- c(ar, numeric(max(r - p, 0)))
ema <- c(ma, numeric(max(r - q - 1, 0)))
Tt <- cbind(ear, rbind(diag(r - 1), 0))
Rt <- matrix(c(1, ema), ncol = 1)
Zt <- matrix(c(1, numeric(r - 1)), nrow = d)
## TODO: more flexible options for ct i dt ?
## obs. regression
ct <- if(missing(xreg) || length(xreg) == 0)
matrix(0, nrow = d)
else
xreg
## state regression
dt <- if(missing(regx) || length(regx) == 0)
matrix(0, nrow = r, ncol = 1)
else if(nrow(regx) == r)
regx
else # assume < r
rbind(regx, matrix(0, nrow = r - nrow(regx), ncol = ncol(regx)))
GGt <- matrix(0, nrow = d) # d x d
H <- Rt * sigma
HHt <- H %*% t(H)
a0 <- numeric(r)
P0 <- diag(1e6, nrow = r)
list(a0 = a0, P0 = P0, ct = ct, dt = dt,
Zt = Zt, Tt = Tt, GGt = GGt, HHt = HHt)
}
## extention of xarmaxss
## fkf: y ~ d x 1
## state ~ m x 1
xarimaxss <- function(ar, ma, sigma, xreg, regx, delta = numeric(0)) {
model <- xarmaxss(ar, ma, sigma, xreg, regx)
dd <- length(delta)
if(dd == 0)
return(model)
model$dt <- c(model$dt, numeric(dd)) # extend dt with zeroes ## TODO: dt is a matrix!
model$Tt <- rbind(c(delta, 1, numeric(ncol(model$Tt) - 1)),
dbind(diag(1, nrow = dd - 1, ncol = dd),
model$Tt ))
model$Ht <- dbind(matrix(0, dd, dd), model$Ht)
model$Zt <- c(model$Zt, numeric(dd)) # [1, 0, ...., 0]
model$a0 <- c(numeric(dd), model$a0)
model$P0 <- dbind(diag(1e6, nrow = dd), model$P0)
# list(a0 = a0, P0 = P0, ct = ct, dt = dt,
# Zt = Zt, Tt = Tt, GGt = GGt, HHt = HHt)
model
}
## modification of makeARIMA() from ~R/src/base/R-3.3.2/src/library/stats/R/~
makeArimaGnb <- function(phi, theta, Delta, kappa = 1e6,
SSinit = "gnb", # c("Rossignol2011", "gnb", "Gardner1980")
tol = .Machine$double.eps){
if(anyNA(phi)) warning(gettextf("NAs in '%s'", "phi"), domain=NA)
if(anyNA(theta)) warning(gettextf("NAs in '%s'", "theta"), domain=NA)
p <- length(phi)
q <- length(theta)
d <- length(Delta)
r <- max(p, q + 1L)
rd <- r + d
V <- Pn <- P <- T <- matrix(0., rd, rd)
h <- 0.
a <- rep(0., rd)
Z <- c(1., rep.int(0, r-1L), Delta)
## if(q < r - 1L)
## theta <- c(theta, rep.int(0, r-1L-q))
## R <- c(1, theta, rep.int(0, d))
## V <- R %o% R
if(q > 0)
V[1:(q+1), 1:(q+1)] <- c(1, theta) %o% c(1, theta)
else
V[1, 1] <- 1.0
if(q < r - 1L) # extent theta - this is in the returned value!
theta <- c(theta, rep.int(0, r - 1L - q))
## The stationary part of T is left-companion of dim (r,r)
if(p > 0)
T[1L:p, 1L] <- phi
if(r > 1L) {
ind <- 2:r
T[cbind(ind-1L, ind)] <- 1
}
## set the top left corner of Pn to the covariance matrix of the stationary part
if(r > 1L)
Pn[1L:r, 1L:r] <- switch(match.arg(SSinit),
# "Gardner1980" = .Call(stats:::C_getQ0, phi, theta),
# "Rossignol2011" = .Call(stats:::C_getQ0bis, phi, theta, tol),
"gnb" = .Call(`_sarima_arma_Q0gnb`, phi, theta, tol),
stop("invalid 'SSinit'"))
else Pn[1L, 1L] <- if(p > 0) 1/(1 - phi^2) else 1
if(d > 0L) {
## The nonstationary part of T is top-companion of dim (d,d);
## the r+1 row of T is equal to Z: its first element is 1,
## followed by zeroes, followed by Delta from column r + 1.
T[r + 1L, ] <- Z
if(d > 1L) {
ind <- r + 2:d
T[cbind(ind, ind-1)] <- 1
}
## set the diagonal of the bottom right corner of Pn to kappa
## (this corresponds to the nonstationary part)
Pn[cbind(r + 1L:d, r + 1L:d)] <- kappa
}
## note P is left equal to the zero matrix
list(phi = phi, theta = theta, Delta = Delta, Z = Z, a = a,
P = P, T = T, V = V, h = h, Pn = Pn)
}
sarimaReport <- function(o1, o2, ...){
list(
all.equal(coef(o1), coef(o2))
)
}
factorizeMA <- function(x, theta = NULL, tol = 1e-12, maxiter = 1000){
q <- length(x) - 1
if(is.null(theta))
theta <- c(sqrt(x[1] + 2 * sum(x[-1])), rep(0, q))
storage.mode(theta) <- "double"
# r <- ltsa::tacvfARMA(theta = -theta[-1]/theta[1], maxLag = q, sigma2 = theta[1]^2)
r <- .Call("_sarima_MAacvf0", theta)
#cat( "theta = ", theta, ",\t r = ", r, "\n")
crit <- 1e100
iter <- 0
while(crit > tol && iter <= maxiter){
iter <- iter + 1
T <- .Call("_sarima_DAcvfWrtMA", theta)
theta <- solve(T, r + x)
r <- ltsa::tacvfARMA(theta = -theta[-1]/theta[1], maxLag = q, sigma2 = theta[1]^2)
rnew <- .Call("_sarima_MAacvf0", as.double(theta)) # MAacvf0(theta)
#cat( "theta = ", theta, ",\t r = ", r, "\n")
#browser()
crit <- sum((x - r)^2)
}
## output as SQUAREM::fpiter(), maybe should use it anyway.
list(par = theta, value.objfn = crit, fpevals = iter,
objfevals = iter, convergence = if(iter > maxiter) 1 else 0)
}
## arimaSS <- function(y, mod)
## {
## ## next call changes mod components a, P, Pn so beware!
## .Call(stats:::C_ARIMA_Like, y, mod, 0L, TRUE)
## }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.