Nothing
tsdiag.MixAR <- function(object, gof.lag = NULL, y, ask = interactive(), ...,
plot = interactive(), std.resid = FALSE){
old.par <- par(no.readonly = TRUE)
on.exit(par(old.par)) # restore graphics parameters before exiting.
## changing 'f' to 'model'. 'model' is used below only in the call to
## mixAR_BIC but the loglik is calculated by the latter anyway. Also, 'f' is
## definitely not suitable name for this variable.
##
## changing 'model' to 'object' for consistency with 'tsdiag'
if(is.list(object))
object <- object$model
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)
n <- length(y)
p <- max(object@order)
g <- length(object@prob)
mu <- mix_hatk(object, y, index = (p+1):n)@m
prob <- object@prob
sigma <- object@scale
ar <- object@arcoef@a
shift <- object@shift
err <- c(y[-c(1:p)] - (mu %*% prob)) #rowSums(t(t(mu) * prob))
if(std.resid)
err <- err / sqrt(mix_variance(object, y, index = (p+1):n))
LBtest <-
function(x, lags){
m <- sapply(lags, function(lag)
c(lag, Box.test(x, lag = lag, type = "Ljung-Box")$p.value))
res <- t(m)
colnames(res) <- c("df", "p.value")
rownames(res) <- paste0("Lag_", lags)
res
}
# lb <- Box.test(err, lag = gof.lag, type = "Ljung-Box")
#or:
# lb <- sapply(gof.lag,
# function(lag) c(lag, Box.test(err, lag = lag, type = "Ljung-Box")$p.value))
# lb <- t(lb)
# colnames(lb) <- c("df", "p.value")
# rownames(lb) <- paste0("Lag_", gof.lag)
lb <- LBtest(err, gof.lag)
cdf <- mix_cdf(object, y, index = (p + 1):n)
lb_cdf <- LBtest(cdf, gof.lag)
# or: ek_new <- mix_ek(object, y, index = (p+1):n, scale = TRUE)
ek_new <- (y[-c(1:p)] - mix_hatk(object, y, index = (p+1):n)) / sigma
pdf_k_new <- noise_dist(object, "pdf") %of% ek_new
tau <- (prob * (pdf_k_new /sigma))@m
## tau <- mu
## for(k in 1:g){
## tau[,k] <- prob[k] * dnorm(y[-c(1:p)], mu[,k], sigma[k])
## }
##
## The calculation for tau can be replaced by any of the following. At the same time it
## becomes valid for any noise distributions.
##
## ## the following are equivalent - if uncommented the stopifnot() conditions pass 'R CMD
## ## check'
## tau_new <- (prob * (noise_dist(object, "pdf") %of%
## ((y[-c(1:p)] - mix_hatk(object, y, index = (p+1):n)) / sigma) /sigma))@m
## ## tau_new split on several lines for clarity
## ##
## ## can be replaced by
## ## ek_new <- mix_ek(object, y, index = (p+1):n, scale = TRUE)
## ## but check
## ek_new <- (y[-c(1:p)] - mix_hatk(object, y, index = (p+1):n)) / sigma
## pdf_k_new <- noise_dist(object, "pdf") %of% ek_new
## tau_new2 <- (prob * (pdf_k_new /sigma))@m
##
## stopifnot(max(abs(tau - tau_new)) < 1e-20)
## stopifnot(max(abs(tau_new2 - tau_new)) < 1e-20)
ks <- ks.test(cdf, punif)
v <- qnorm(cdf)
sw <- shapiro.test(v)
lb_v <- LBtest(v, gof.lag)
index <- apply(tau, 1, which.max) # function(x) which.max(x)
err2 <- numeric(n-p)
for(t in seq_along(index)){
err2[t] <- (y[t+p] - mu[t, index[t]])/sigma[index[t]]
}
sw2 <- shapiro.test(err2)
lb_err2 <- LBtest(err2, gof.lag)
choices <- c(
"ACF/PACF of residuals",
"ACF/PACF of U_residuals",
"ACF/PACF of tau_residuals",
"ACF/Histogram of tau_residuals"
)
chnum <- 1:length(choices)
if(!isTRUE(plot)){ # plot is typically numeric index;
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(choices), ",\n", "ignoring non-existent values")
chnum <- chnum[!is.na(choices)]
choices <- choices[!is.na(choices)]
}
}
if(length(choices) > 0){
par(mfrow = c(2,1))
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]
while(length(choice) != 0){
switch(choice,
{ # 1: "ACF/PACF of residuals",
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/PACF of U_residuals"
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: "ACF/PACF of tau_residuals"
acf(err2, main = "ACF of tau_residuals", lag.max = lag.max)
pacf(err2, main = "PACF of tau_residuals", lag.max = lag.max)
},
{ # 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(length(chnum) == 1) # length(choices) == 1
break
## TODO: argument 'ask' could be used here to present a menu or just
## plot the next plot in choices.
if(ask)
choice <- chnum[ menu(choices, title = choice_title) ]
else{
## just plot the next one
chnum <- chnum[-1]
choice <- chnum[1]
}
}
}
BIC <- mixAR_BIC(y, object)
# res <- list("residuals" = err, "U" = cdf, "V" = v, "tau_residuals" = err2,
# "tests" = list("Ljung-Box" = lb, "K-S" = ks, "Shapiro-Wilk" = sw,
# "Shapiro-Wilk_tau" = sw2),
# lb_cdf = lb_cdf, lb_err2 = lb_err2,
# BIC = BIC)
# res <- list(residuals = list(values = err, "Ljung-Box" = lb),
# U_residuals = list(values = cdf, , "Ljung-Box" = lb_cdf, "KS" = ks),
# V_residuals = list(values = v, "Shapiro-Wilk" = sw),
# tau_residuals = list(values = err2, "Ljung-Box" = lb_err2,
# "Shapiro-Wilk" = sw2),
# BIC = BIC)
# class(res) <- "diagMixAR"
res <- structure(
list(
residuals = list(values = err, "Ljung-Box" = lb),
U_residuals = list(values = cdf, "Ljung-Box" = lb_cdf, "KS" = ks),
V_residuals = list(values = v, "Ljung-Box" = lb_v, "Shapiro-Wilk" = sw),
tau_residuals = list(values = err2, "Ljung-Box" = lb_err2, "Shapiro-Wilk" = sw2),
BIC = BIC),
class = "tsdiagMixAR")
invisible(res)
}
print.tsdiagMixAR <- function(x, ...){
## for now just drop the values of the residuals
dashes <- "--------------------------------------------------"
cat("\n", dashes, "\nTests on the ordinary residuals\n", dashes, "\n")
print(x$residuals[-1])
cat("\n", dashes, "\nTests on the U_residuals\n", dashes, "\n")
print(x$U_residuals[-1])
cat("\n", dashes, "\nTests on the V_residuals\n", dashes, "\n")
print(x$V_residuals[-1])
cat("\n", dashes, "\nTests on the tau_residuals\n", dashes, "\n")
print(x$tau_residuals[-1])
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.