fixes for CRAN submission

Author: alexiosg

NAMESPACE

@@ -21,6 +21,7 @@ S3method(plot,tsissm.simulate)
 S3method(predict,tsissm.estimate)
 S3method(predict,tsissm.selection)
 S3method(print,summary.tsissm.estimate)
+S3method(print,tsissm.diagnose)
 S3method(residuals,tsissm.estimate)
 S3method(sigma,tsissm.estimate)
 S3method(simulate,tsissm.estimate)
@@ -171,7 +172,6 @@ importFrom(utils,head)
 importFrom(utils,setTxtProgressBar)
 importFrom(utils,tail)
 importFrom(utils,txtProgressBar)
-importFrom(utils,write.table)
 importFrom(viridisLite,viridis)
 importFrom(xts,as.xts)
 importFrom(xts,is.xts)

R/estimation.R

@@ -113,8 +113,8 @@ estimate.tsissm.autospec <- function(object, control = NULL, trace = FALSE, ...)
         tmp_toc <- Sys.time() - tmp_tic    
         est_time_one <- tmp_toc
         estimated_time <- round(as.numeric(est_time_one/n_cores) * NROW(args_grid), 2) %/% 60
-        print(paste0("no. of models to evaluate: ", NROW(args_grid)))
-        print(paste0("estimated evaluation time (mins): ", estimated_time))
+        cat(paste0("no. of models to evaluate: ", NROW(args_grid)))
+        cat(paste0("estimated evaluation time (mins): ", estimated_time))
         prog_trace <- progressor(n)
     }
     tic <- Sys.time()

R/methods.R

@@ -5,6 +5,7 @@
 #' @description Extract the fitted values from an estimated model.
 #' @param object an object of class \dQuote{tsissm.estimate}.
 #' @param ... not currently used.
+#' @returns an xts object of the fitted values.
 #' @aliases fitted
 #' @method fitted tsissm.estimate
 #' @rdname fitted
@@ -941,6 +942,8 @@ BIC.tsissm.selection <- function(object, ...)
 #' @param actual the actual data matched to the dates of the forecasts.
 #' @param alpha the coverage level for distributional forecast metrics.
 #' @param ... not currently used.
+#' @returns a data.frame of performance metrics including MAPE, MSLRE, BIAS and AIC
+#' for the estimate object and MAPE, MSLRE, BIAS, MASE, MIS and CRPS for predict object.
 #' @aliases tsmetrics
 #' @method tsmetrics tsissm.predict
 #' @rdname tsmetrics
@@ -979,21 +982,8 @@ tsmetrics.tsissm.estimate <- function(object, ...)
     MAPE <- mape(object$spec$target$y_orig, fitted(object))
     BIAS <- bias(object$spec$target$y_orig, fitted(object))
     MSLRE <- mslre(object$spec$target$y_orig, fitted(object))
-    # yt <- object$spec$transform$transform(object$spec$target$y_orig, lambda = lambda)
-    # ft <- object$spec$transform$transform(as.numeric(fitted(object)), lambda = lambda)
-    # r <- yt - ft
-    cat("\ntsissm: Performance Metrics")
-    cat("\n----------------------------------\n")
-    cat(paste0("AIC\t: ", round(AIC,2), " (n = ", nr,")"))
-    cat("\n")
-    cat(paste0("MAPE\t: ", round(MAPE,5)))
-    cat("\n")
-    cat(paste0("BIAS\t: ", round(BIAS,5)))
-    cat("\n")
-    cat(paste0("MSLRE\t: ", round(MSLRE, 5)))
-    metrics = c(AIC, MAPE, BIAS, MSLRE)
-    names(metrics) <- c("AIC","MAPE","BIAS","MSLRE")
-    return(invisible(metrics))
+    metrics <- data.frame("MAPE" = MAPE, "MSLRE" = MSLRE, "BIAS" = BIAS, "AIC" = AIC)
+    return(metrics)
 }
 
 # vcov ---------------------------------------------------
@@ -1070,7 +1060,7 @@ sigma.tsissm.estimate <- function(object, ...) {
 #' @description Generates a list of model equations in LaTeX.
 #' @param object an object of class \dQuote{tsissm.estimate}.
 #' @param ... not currently used.
-#' @return A list of equations in LaTeX which can be used in documents. This is
+#' @returns A list of equations in LaTeX which can be used in documents. This is
 #' a list with 3 slots for the observation, state and variance equations,
 #' @details This method is called in the summary when the format output option
 #' chosen is \dQuote{flextable}.

R/plots.R

@@ -15,7 +15,8 @@
 #'
 plot.tsissm.estimate <- function(x, y = NULL, ...)
 {
-    opar <- par(mfrow = c(1,1))
+    opar <- par(no.readonly = TRUE)
+    on.exit(par(opar))
     tsd <- tsdecompose(x)
     # fitted+actual and then states
     a <- x$spec$target$y_orig
@@ -35,15 +36,15 @@ plot.tsissm.estimate <- function(x, y = NULL, ...)
         mtext(colnames(tsd)[i], side = 4, adj = 0.5, padj = 0.5, cex = 0.7, font = 2, family = "mono")
         grid()
     }
-    suppressWarnings(par(opar))
 }
 
 #' @method plot tsissm.simulate
 #' @rdname plot
 #' @export
 plot.tsissm.simulate <- function(x, y = NULL, ...)
 {
-    opar <- par(mfrow = c(1,1))
+    opar <- par(no.readonly = TRUE)
+    on.exit(par(opar))
     components <- tsdecompose(x)
     n <- length(components) + 1
     colx <- (.viridis_fun(option = "H", begin = 0.4, end = 0.9, alpha = 0.5)(n))
@@ -57,7 +58,6 @@ plot.tsissm.simulate <- function(x, y = NULL, ...)
         plot(components[[component_names[i]]], gradient_color = colx[i], main = "", ylab = "", xlab = "", x_axes = FALSE, cex.axis = 0.8, interval_color = "steelblue", median_width = 1, interval_type = 1, interval_width = 1)
         mtext(component_names[i], side = 4, adj = 0.5, padj = 0.5, cex = 0.7, font = 2, family = "mono")
     }
-    suppressWarnings(par(opar))
 }
 
 #' @method plot tsissm.profile
@@ -66,7 +66,8 @@ plot.tsissm.simulate <- function(x, y = NULL, ...)
 plot.tsissm.profile <- function(x, y = NULL, type = c("coef","mape","mase","crps"), ...)
 {
     Simulation <- Variable <- NULL
-    opar <- par(mfrow = c(1,1))
+    opar <- par(no.readonly = TRUE)
+    on.exit(par(opar))
     type <- match.arg(type[1], c("coef","mape","mase","crps"))
     if (type == "coef") {
         true_values <- data.table(Variable = names(x$true_coef),
@@ -99,6 +100,5 @@ plot.tsissm.profile <- function(x, y = NULL, type = c("coef","mape","mase","crps
         colx <- rev(.viridis_fun(alpha = 0.5, begin = 0, end = 0.7, option = "H")(ncol(tmp)))
         boxplot(round(tmp, 2), xlab = "Horizon", ylab = "CRPS", col = colx, outline = FALSE, main = "CRPS by Horizon")
     }
-    suppressWarnings(par(opar))
     return(invisible(x))
 }

R/tsdiagnose.R

@@ -4,9 +4,10 @@
 #' @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 (\code{\link[EnvStats]{rosnerTest}}).
+#' @returns A list of tests including the weighted Ljung-Box test for residual 
+#' autocorrelation, system forecastability test, and outlier dates using the 
+#' Rosner test (\code{\link[EnvStats]{rosnerTest}}). Optionally generates a plot
+#' of relevant diagnostics.
 #' @aliases tsdiagnose
 #' @method tsdiagnose tsissm.estimate
 #' @rdname tsdiagnose
@@ -16,29 +17,13 @@
 tsdiagnose.tsissm.estimate <- function(object, plot = FALSE, ...)
 {
     if (sum(object$spec$arma$order) > 0) {
-        cat("\nARMA roots (<1)")
-        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")
-    cat("\n------------------------------------------\n")
     e <- abs(eigen(object$model$D, symmetric = FALSE)$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)
     sigma <- sigma(object)
     r <- as.numeric(na.omit(residuals(object, transformed = TRUE)/sigma))
@@ -53,25 +38,72 @@ tsdiagnose.tsissm.estimate <- function(object, plot = FALSE, ...)
     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 <- .rosner_test(as.numeric(na.omit(residuals(object, transformed = TRUE))), k = 10)
     if (any(rtest$Outlier)) {
-        out.index <- object$spec$target$index[which(object$spec$good == 1)][rtest$Obs.Num[rtest$Outlier]]
-        cat("\nOutlier Diagnostics (based on Rosner Test)")
-        cat("\n------------------------------------------")
-        cat("\nOutliers:", as.character(out.index))
+        outliers_index <- object$spec$target$index[which(object$spec$good == 1)][rtest$Obs.Num[rtest$Outlier]]
     } else {
-        out.index <- NULL
+        outliers_index <- NULL
     }
     if (plot) {
+        opar <- par(no.readonly = TRUE)
+        on.exit(par(opar))
         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)
+        hist(r, breaks = "fd", main = "Scaled Residuals Histogram", probability = TRUE)
         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))
+    L <- list(arma_test = rt, 
+              stability_test = e, 
+              weighted_box_test = lbsr, 
+              rosner_test = rtest, 
+              outliers_index = outliers_index)
+    class(L) <- "tsissm.diagnose"
+    return(L)
+}
+
+
+#' Model Diagnostics Print method
+#'
+#' @description Print method for class \dQuote{tsissm.diagnose}
+#' @param x an object of class \dQuote{tsissm.duagnose} generated from
+#' calling \code{\link[tsissm]{tsdiagnose}}.
+#' @param ... not currently used.
+#' @returns Invisibly returns the original object and prints the output to console.
+#' @aliases print.tsissm.diagnose
+#' @method print tsissm.diagnose
+#' @rdname print.tsissm.diagnose
+#' @export
+#'
+#'
+print.tsissm.diagnose <- function(x, ...)
+{
+    arma_test <- x$arma_test
+    if (!is.null(arma_test)) {
+        cat("\nARMA roots (<1)")
+        cat("\n------------------------------------------\n")
+        if (!is.null(arma_test$root$ar)) {
+            cat("Inverse AR roots:", 1/abs(arma_test$root$ar))
+            cat("\n")
+        }
+        if (!is.null(arma_test$root$ma)) {
+            cat("Inverse MA roots:", 1/abs(arma_test$root$ma))
+            cat("\n")
+        }
+    }
+    cat("\nForecastability")
+    cat("\n------------------------------------------\n")
+    cat("Real Eigenvalues (D):", round(x$stability_test,3))
+    cat("\n")
+    cat("\nWeighted Ljung-Box Test [scaled residuals]")
+    cat("\n------------------------------------------\n")
+    print(x$weighted_box_test, row.names = FALSE, digits = 3)
+    if (any(x$rosner_test$Outlier)) {
+        cat("\nOutlier Diagnostics (based on Rosner Test)")
+        cat("\n------------------------------------------")
+        cat("\nOutliers:", as.character(x$outliers_index))
+    }
+    return(invisible(x))
 }

R/tsissm-package.R

@@ -5,7 +5,7 @@
 #' @import data.table
 #' @importFrom RTMB AD ADoverload ADREPORT REPORT matrix eigen as.vector.advector
 #' @importFrom TMB MakeADFun
-#' @importFrom utils head tail data txtProgressBar setTxtProgressBar write.table
+#' @importFrom utils head tail data txtProgressBar setTxtProgressBar
 #' @importFrom stats sd var acf na.pass pchisq pgamma qqline qqnorm simulate na.omit median fitted coef quantile residuals predict logLik cov cor qt pnorm AIC BIC nobs sigma ar arima dnorm printCoefmat vcov density
 #' @importFrom graphics grid legend lines par plot points abline axis axis.Date axis.POSIXct box polygon layout mtext title hist boxplot
 #' @importFrom grDevices gray n2mfrow

R/tsprofile.R

@@ -74,8 +74,8 @@ profile_fun <- function(sim, object, h, control, trace)
     toc <- Sys.time()
     dtime <- difftime(toc, tic, units = "secs")
     if (trace) {
-        print(paste0("\nCompleted Profiling in ", round(as.numeric(dtime),2)," secs."))
-        print(paste0("\nCompiling Performance Metrics..."))
+        cat(paste0("\nCompleted Profiling in ", round(as.numeric(dtime),2)," secs."))
+        cat(paste0("\nCompiling Performance Metrics..."))
     }
     C <- rbindlist(lapply(1:length(prof), function(i) prof[[i]]$L1))
     M <- rbindlist(lapply(1:length(prof), function(i) prof[[i]]$L2))