Nothing
#' @import graphics
#'
#' @title Plot method for class 'gsmarpred' objects
#'
#' @description \code{plot.gsmarpred} is plot method for class 'gsmarpred' objects
#'
#' @inheritParams predict.gsmar
#' @param x object of class \code{'gsmarpred'} created with \code{predict.gsmar}.
#' @param add_grid should grid a be added to the plots?
#' @param ... arguments passed to function \code{grid}.
#' @details This method is intended for plotting forecasts of GSMAR processes.
#' @export
plot.gsmarpred <- function(x, ..., nt, mix_weights=TRUE, add_grid=TRUE) {
# Pick statistics etc
gsmarpred <- x
data <- as.ts(gsmarpred$gsmar$data)
n_obs <- length(data)
q <- gsmarpred$q
M <- sum(gsmarpred$gsmar$model$M)
mix_weights <- gsmarpred$mix_weights & M > 1
mixing_weights <- gsmarpred$gsmar$mixing_weights
n_mix <- nrow(mixing_weights)
# Graphical parameters
old_par <- par(no.readonly=TRUE) # Save old settings
on.exit(par(old_par)) # Restore the settings before quitting
if(mix_weights) {
par(mfrow=c(2, 1), mar=c(2.1, 2.6, 2.1, 1.1))
} else {
par(mfrow=c(1, 1), mar=c(2.6, 2.6, 2.1, 1.1))
}
# How many observations to plot preceding the prediction
if(missing(nt)) {
nt <- round(length(data)*0.15)
} else {
stopifnot(nt > 0 & nt %% 1 == 0)
if(nt > length(data)) {
warning("nt > length(data); using nt = length(data)")
nt <- length(data)
}
}
# Function that creates ts objects to plot
make_ts <- function(a, mix=FALSE, m) { # a is vector of values; if mix == TRUE, use component m mixing weight last observation as the first value
last_obs <- ifelse(mix, mixing_weights[n_mix, m], data[n_obs])
ts(c(last_obs, a), start=time(data)[n_obs], frequency=frequency(data))
}
# Prediction time series and all values
ts_pred <- make_ts(gsmarpred$pred)
ts_dat <- ts(data[(n_obs - nt + 1):n_obs], start=time(data)[n_obs - nt + 1], frequency=frequency(data))
t0 <- time(ts_pred)
all_val <- c(ts_dat, ts_pred, gsmarpred$pred_ints)
# Prediction invervals
if(gsmarpred$pi_type == "two-sided") {
# Functions that create prediction interval time series
ts1_lapply <- function(pred_ints, mix=FALSE, m) lapply(1:(length(q)/2), function(i1) make_ts(pred_ints[,i1], mix, m)) # Lower bounds
ts2_lapply <- function(pred_ints, mix=FALSE, m) lapply((length(q)/2 + 1):length(q), function(i1) make_ts(pred_ints[,i1], mix, m)) # Upper bounds
ints1 <- gsmarpred$pred_ints # Process pi
ints1_mix <- gsmarpred$mix_pred_ints # Mixing weight pi
} else { # pi_type = upper, lower, or none
# If upper/lower/none pi, at least one of the "bounds" is just constants
what_to_rep <- ifelse(gsmarpred$pi_type == "upper", round(min(all_val)) - 3, round(max(all_val)) + 3) # Otherwise "lower" or "none"
what_to_rep_mix <- ifelse(gsmarpred$pi_type == "upper", 0, 1)
# Functions that create prediction interval time series
ts1_lapply <- function(pred_ints, mix=FALSE, m) lapply(seq_along(q), function(i1) make_ts(pred_ints[,i1], mix, m)[-1]) # Upper or lower graphical device box bound (redundant argument)
ts2_lapply <- function(pred_ints, mix=FALSE, m) lapply(seq_along(q), function(i1) make_ts(pred_ints[,i1], mix, m)) # Make upper or lower prediction bound
# Constant lower or upper "bound" for one-sided prediction intervals
ints1 <- matrix(rep(what_to_rep, times=length(q)*length(ts_pred)), ncol=length(q))
ints1_mix <- array(rep(what_to_rep_mix, times=length(ts_pred)*length(q)*M), dim=c(length(ts_pred), length(q), M))
}
# Create the ts objects for confidence intervals
ts1 <- ts1_lapply(ints1)
ts2 <- ts2_lapply(gsmarpred$pred_ints)
# Mixing weight prediction intervals
if(mix_weights) {
ts1_mw <- vector(mode="list", length=M) # Sublist for each regime
ts2_mw <- vector(mode="list", length=M)
for(m in 1:M) {
ts1_mw[[m]] <- ts1_lapply(matrix(ints1_mix[, , m], ncol=length(q)), mix=TRUE, m=m)
ts2_mw[[m]] <- ts2_lapply(matrix(gsmarpred$mix_pred_ints[, , m], ncol=length(q)), mix=TRUE, m=m)
}
}
# Plot the point predictions
ts.plot(ts_dat, ts_pred, gpars=list(col=c("black", "blue"), lty=1:2,
ylim=c(round(min(all_val)) - 1,
round(max(all_val)) + 1),
main=paste("Forecast", gsmarpred$n_ahead, "steps ahead")))
if(add_grid) grid(...)
# Plot the prediction intervals
draw_poly <- function(ts1_or_ts2, pred_ts, col) polygon(x=c(t0, rev(t0)), y=c(ts1_or_ts2, rev(pred_ts)), col=col, border=NA)
col_pred <- grDevices::rgb(0, 0, 1, 0.2)
if(gsmarpred$pi_type %in% c("two-sided", "upper", "lower")) {
for(i1 in 1:length(gsmarpred$pi)) { # Go through the prediction intervals
draw_poly(ts1[[i1]], ts_pred, col=col_pred)
draw_poly(ts2[[i1]], ts_pred, col=col_pred)
}
}
# Plot mixing weight predictions
if(mix_weights) {
# Plot point predictions
colpal_mw <- grDevices::colorRampPalette(c("blue", "turquoise1", "green", "red"))(M)
colpal_mw2 <- grDevices::adjustcolor(colpal_mw, alpha.f=0.5)
mix_ts <- ts(mixing_weights[(n_mix - nt + 1):n_mix,], start=time(data)[n_obs - nt + 1],
frequency=frequency(data))
mix_pred_ts <- ts(rbind(mix_ts[nrow(mix_ts),], gsmarpred$mix_pred), start=time(data)[n_obs], frequency=frequency(data))
ts.plot(mix_ts, mix_pred_ts, gpars=list(col=c(colpal_mw2, colpal_mw), ylim=c(0, 1), lty=c(rep(1, M), rep(2, M))),
main="Mixing weights")
legend("topleft", legend=paste0("regime ", 1:M), bty="n", col=colpal_mw, lty=1, lwd=2,
text.font=2, cex=0.70, x.intersp=0.3, y.intersp=1, seg.len=1, inset=c(-0.01, -0.035))
if(add_grid) grid(...)
# Plot mixing weight prediction intervals
colpal_mw3 <- grDevices::adjustcolor(colpal_mw, alpha.f=0.2)
if(gsmarpred$pi_type %in% c("two-sided", "upper", "lower")) { # Don't plot if pi_type == "none"
for(m in 1:M) { # Go through regimes
for(i1 in 1:length(gsmarpred$pi)) { # Go through the prediction intervals
draw_poly(ts1_mw[[m]][[i1]], mix_pred_ts[, m], col=colpal_mw3[m])
draw_poly(ts2_mw[[m]][[i1]], mix_pred_ts[, m], col=colpal_mw3[m])
}
}
}
}
invisible(gsmarpred)
}
#' @import graphics
#' @describeIn quantile_residual_tests Plot p-values of the autocorrelation and conditional
#' heteroskedasticity tests.
#' @inheritParams print.qrtest
#' @export
plot.qrtest <- function(x, ...) {
# Graphical settings
old_par <- par(no.readonly = TRUE)
on.exit(par(old_par))
qrtest <- x
par(mfrow=c(1, 2), mar=c(5.1, 3.1, 3.1, 1.1))
# Function to plot p-values
plot_pvalues <- function(which_ones) { # ac_res, ch_res
res <- qrtest[[which(names(qrtest) == which_ones)]]
pvals <- res$pvalue
seq_pvals <- seq_along(pvals)
plot(pvals, ylim=c(0, 1), xlim=c(min(seq_pvals) - 0.2, max(seq_pvals) + 0.2), ylab="", xlab="lags",
main=ifelse(which_ones == "ac_res", "Autocorrelation", "Cond. h.skedasticity"),
xaxt="n", yaxt="n", pch=16, col="blue")
axis(side=1, at=seq_pvals, labels=res$lags)
levels <- c(0.01, 0.05, 0.10, seq(from=0.20, to=1.00, by=0.20))
axis(side=2, at=levels, las=1, cex.axis=0.8)
abline(h=0, lwd=2)
abline(h=c(0.01, 0.05, 0.10, 1.00), lty=2, col=c("red", "red", "red", "darkgreen"))
segments(x0=seq_pvals, y0=0, y1=pvals, x1=seq_pvals, ...)
points(pvals)
}
# Plot the p-values
plot_pvalues("ac_res") # Autocorrelation tests
plot_pvalues("ch_res") # Conditional heteroskedasticity tests
invisible(qrtest)
}
#' @import graphics
#' @describeIn GSMAR Plot method for class 'gsmar'
#' @param x object of class \code{'gsmar'} created with \code{fitGSMAR} or \code{GSMAR}.
#' @param ... in the plot method: arguments passed to the function \code{density} which
#' calculates the kernel density estimate of the data.
#' @param include_dens Plot also kernel density estimate of the data and model implied stationary
#' density with regimewise densities? See the details.
#' @details If \code{include_dens == TRUE}, the kernel density estimate of the data is calculated
#' with the function \code{density} from the package \code{stats}. By default, the default
#' settings of that function are used, including the usage of Gaussian kernel. Use the dot
#' parameters to adjust the settings if desired.
#'
#' By the model implied stationary density we mean the stationary one-dimensional mixture
#' density of M regimes (see KMS 2015, Theorem 1 and the discussion following it for the Gaussian
#' case and Theorem 2 in PMS 2018 for the Student's t case). The regimewise densities (i.e. each
#' density 1,...,M in the stationary mixture density) are multiplied with the mixing weight
#' parameters accordingly.
#'
#' In the density plot black represents the kernel density estimate of the data, grey dashed line the
#' model implied density, and the colored dotted lines the regime wise densities.
#' @export
plot.gsmar <- function(x, ..., include_dens=TRUE) {
# Pick the relevant statistics etc
gsmar <- x
check_data(gsmar)
data <- as.ts(gsmar$data)
n_obs <- length(data)
p <- gsmar$model$p
M <- sum(gsmar$model$M)
ts_mw <- ts(rbind(matrix(NA, nrow=p, ncol=M), as.matrix(gsmar$mixing_weights)),
start=start(data), frequency=frequency(data)) # First p observations are starting values
# Graphical settings
old_par <- par(no.readonly=TRUE)
on.exit(par(old_par))
if(include_dens) {
layout(mat=matrix(c(1, 3, 2, 3), nrow=2, byrow=TRUE), widths=c(2, 1))
par(mar=c(2.5, 2.3, 1.5, 0.3))
} else {
par(mfrow=c(2, 1), mar=c(2.5, 2.5, 2.1, 1))
}
colpal_mw <- grDevices::colorRampPalette(c("blue", "turquoise1", "green", "red"))(M)
# Plot the time series and mixing weights
ts.plot(data, gpars=list(main="Time series"))
ts.plot(ts_mw, gpars=list(main="Mixing weights", ylim=c(0, 1), col=colpal_mw, lty=2))
legend("topleft", legend=paste0("regime ", 1:M), bty="n", col=colpal_mw, lty=1, lwd=2,
text.font=2, cex=0.75, x.intersp=0.2, y.intersp=1, seg.len=1, inset=c(-0.03, -0.035))
# Plot kernel density estimate of the data with the model implied density
if(include_dens) {
# Collect the required statistics
params <- gsmar$params
M_orig <- gsmar$model$M
model <- gsmar$model$model
restricted <- gsmar$model$restricted
constraints <- gsmar$model$constraints
alphas <- pick_alphas(p=p, M=M_orig, params=params, model=model, restricted=restricted, constraints=constraints)
means <- get_regime_means(gsmar)
vars <- get_regime_vars(gsmar)
if(model == "GMAR") {
M1 <- M
M2 <- 0
} else if(model == "StMAR") {
M1 <- 0
M2 <- M
} else { # model == "G-StMAR"
M1 <- M_orig[1]
M2 <- M_orig[2]
}
# Degrees of freedom and Student's t density
if(model %in% c("StMAR", "G-StMAR")) {
dfs <- pick_dfs(p=p, M=M_orig, params=params, model=model)
# The student's density function as in the model definition
my_dt <- function(y, mean, var, df) {
tmp <- lgamma(0.5*(1 + df)) - lgamma(0.5*df) - 0.5*log(pi*(df - 2)) - 0.5*log(var) -
0.5*(1 + df)*log(1 + (y - mean)^2/(var*(df - 2)))
exp(tmp)
}
}
# Regime density multiplied with the corresponding mixing weight parameter
reg_dens <- function(m, xx) {
if(m <= M1) {
return(alphas[m]*dnorm(xx, mean=means[m], sd=sqrt(vars[m])))
} else {
return(alphas[m]*my_dt(xx, mean=means[m], var=vars[m], df=dfs[m - M1]))
}
}
# The model density
mod_dens_f <- function(xx) {
rowSums(vapply(1:M, function(m) reg_dens(m, xx), numeric(length(xx))))
}
# Densities and figure bounds
data_dens <- density(data)
mod_mean <- gsmar$uncond_moments$uncond_mean
mod_sd <- sqrt(gsmar$uncond_moments$uncond_var)
x0 <- min(mod_mean - 3*mod_sd, min(data_dens$x))
x1 <- max(mod_mean + 3*mod_sd, max(data_dens$x))
xpp <- seq(from=x0, to=x1, length.out=500)
mod_dens <- mod_dens_f(xpp)
y0 <- 0
y1 <- max(c(data_dens$y, mod_dens))
# Plot the densities
col_seriesdens <- "black"
col_modeldens <- "darkgrey"
plot(x=data_dens$x, y=data_dens$y, xlim=c(x0, x1), ylim=c(y0, y1), main="Density",
ylab="", xlab="", cex.axis=0.8, font.axis=2, type="l", col=col_seriesdens)
lines(x=xpp, y=mod_dens, type="l", lty=2, lwd=2, col=col_modeldens)
for(m in 1:M) {
lines(x=xpp, y=reg_dens(m, xx=xpp), type="l", lty=3, col=colpal_mw[m])
}
legend("topleft", legend=c("series", "model", paste0("regime ", 1:M)),
bty="n", col=c(col_seriesdens, col_modeldens, colpal_mw), lty=c(1, 2, 3, 3), lwd=2,
text.font=2, cex=0.75, x.intersp=0.2, y.intersp=1, seg.len=1, inset=c(-0.075, -0.01))
}
invisible(gsmar)
}
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.