#' Model Diagnostics
#'
#' @description Creates a short summary of model based diagnostics.
#' @param object an object of class \dQuote{tsissm.estimate}.
#' @param plot whether to generate diagnostic plots to accompany summary.
#' @param ... not currently used.
#' @return A list of tables (printed out and returned invisibly) with
#' Ljung-Box test for residual autocorrelation, parameter and model bounds
#' diagnostics and outlier dates using the Rosner test from the
#' \code{EnvStats} package.
#' @aliases tsdiagnose
#' @method tsdiagnose tsissm.estimate
#' @rdname tsdiagnose
#' @export
#'
#'
tsdiagnose.tsissm.estimate <- function(object, plot = FALSE, ...)
{
if (sum(object$spec$arma$order) > 0) {
cat("\nARMA roots")
cat("\n------------------------------------------\n")
armav <- coef(object)
armav <- armav[grepl("theta|psi", names(armav))]
rt <- .armaroots(armav)
if (object$spec$arma$order[1] > 0) {
cat("Inverse AR roots:", 1/abs(rt$root$ar))
cat("\n")
}
if (object$spec$arma$order[2] > 0) {
cat("Inverse MA roots:", 1/abs(rt$root$ma))
cat("\n")
}
} else {
rt <- NULL
}
cat("\nForecastability (D roots)")
cat("\n------------------------------------------\n")
e <- abs(Re(eigen(object$model$D)$values))
cat("Real Eigenvalues (D):", round(e,3))
cat("\n")
cat("\nWeighted Ljung-Box Test [scaled residuals]")
cat("\n------------------------------------------\n")
df <- sum(object$spec$arma$order)
r <- as.numeric(na.omit(scale(residuals(object))))
if (sum(object$spec$arma$order) > 0 ) j = 0 else j = 1
b1 <- weighted_box_test(r, lag = 1, type = "Ljung-Box", fitdf = 0)
b2j <- pmax(2 * df + df - 1, 1 + df + j)
b2 <- weighted_box_test(r, lag = b2j, type = "Ljung-Box", fitdf = df)
b3j <- pmax(2 * df + df - 1, 2 + df + j)
b3 <- weighted_box_test(r, lag = b3j, type = "Ljung-Box", fitdf = df)
b4j <- pmax(2 * df + df - 1, 3 + df + j)
b4 <- weighted_box_test(r, lag = b4j, type = "Ljung-Box", fitdf = df)
lbsr <- data.table(Lag = c("Lag[1]", paste0("Lag[",b2j,"]"), paste0("Lag[",b3j,"]"), paste0("Lag[",b4j,"]")),
statistic = c(b1$statistic[[1]], b2$statistic[[1]], b3$statistic[[1]],b4$statistic[[1]]),
pvalue = c(b1$p.value[[1]], b2$p.value[[1]],b3$p.value[[1]], b4$p.value[[1]]))
print(lbsr, row.names = FALSE, digits = 3)
rtest <- rosnerTest(as.numeric(na.omit(residuals(object))), k = 10)
if (any(rtest$all.stats$Outlier)) {
out.index <- object$spec$target$index[which(object$spec$good == 1)][rtest$all.stats$Obs.Num[rtest$all.stats$Outlier]]
cat("\nOutlier Diagnostics (based on Rosner Test)")
cat("\n------------------------------------------")
cat("\nOutliers:", as.character(out.index))
} else {
out.index <- NULL
}
if (plot) {
par(mfrow = c(2,2), mar = c(3,3,3,3))
if (df > 0) .plotarmaroots(.armaroots(armav))
acf(as.numeric(r), type = "correlation", main = "Scaled Residuals Autocorrelation")
hist(r, breaks = "fd", main = "Scaled Residuals Histogram", probability = T)
box()
qqnorm(r)
qqline(r, col = 2)
}
L <- list(armaroots = rt, D.eigenvalues = e, lb_test = lbsr, outliers = rtest$all.stats, outlier_index = out.index)
return(invisible(L))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.