### univariate example
data("salmAllOnset")
## fit a hhh4 model to the first 13 years
salmModel <- list(end = list(f = addSeason2formula(~1 + t)),
ar = list(f = ~1), family = "NegBin1", subset = 2:678)
salmFit <- hhh4(salmAllOnset, salmModel)
## simulate the next 20 weeks ahead
salmSims <- simulate(salmFit, nsim = 300, seed = 3, subset = 678 + seq_len(20),
y.start = observed(salmAllOnset)[678,])
## compare final size distribution to observed value
summary(aggregate(salmSims, time = TRUE)) # summary of simulated values
plot(salmSims, type = "size")
## individual and average simulated time series with a confidence interval
plot(salmSims, type = "time", main = "20-weeks-ahead simulation")
## fan chart based on the quantiles of the simulated counts at each time point
## point forecasts are represented by a white line within the fan
if (requireNamespace("fanplot")) {
plot(salmSims, type = "fan", main = "20-weeks-ahead simulation",
fan.args = list(ln = 1:9/10), means.args = list())
}
### multivariate example
data("measlesWeserEms")
## fit a hhh4 model to the first year
measlesModel <- list(
end = list(f = addSeason2formula(~1), offset = population(measlesWeserEms)),
ar = list(f = ~1),
ne = list(f = ~1 + log(pop),
weights = W_powerlaw(maxlag = 5, normalize = TRUE)),
family = "NegBin1", subset = 2:52,
data = list(pop = population(measlesWeserEms)))
measlesFit1 <- hhh4(measlesWeserEms, control = measlesModel)
## use a Poisson distribution instead (just for comparison)
measlesFit2 <- update(measlesFit1, family = "Poisson")
## simulate realizations from these models during the second year
measlesSims <- lapply(X = list(NegBin = measlesFit1, Poisson = measlesFit2),
FUN = simulate, nsim = 50, seed = 1, subset = 53:104,
y.start = observed(measlesWeserEms)[52,])
## final size of the first model
plot(measlesSims[[1]])
## stratified by groups of districts
mygroups <- factor(substr(colnames(measlesWeserEms), 4, 4))
apply(aggregate(measlesSims[[1]], time = TRUE, units = mygroups), 1, summary)
plot(measlesSims[[1]], groups = mygroups)
## a class and plot-method for a list of simulations from different models
measlesSims <- as.hhh4simslist(measlesSims)
plot(measlesSims)
## simulated time series
plot(measlesSims, type = "time", individual = TRUE, ylim = c(0, 80))
## fan charts
if (requireNamespace("fanplot")) {
opar <- par(mfrow = c(2,1))
plot(measlesSims, type = "fan", which = 1, ylim = c(0, 80), main = "NegBin",
key.args = list())
plot(measlesSims, type = "fan", which = 2, ylim = c(0, 80), main = "Poisson")
par(opar)
}
Run the code above in your browser using DataLab