Description Usage Arguments Details Value Author(s) Examples
forecast
is a generic function for forecasting time series or
time series models. The function invokes particular methods which
depend on the class of the first argument.
Currently the only method provided by the package is for objects of
class "glarma"
.
1 2 3 4 
object 
An object of class "glarma" obtained from a call to

n.ahead 
The number of periods ahead to be forecast. 
newdata 
The model matrix X comprising the values of the predictors for the times for which the series is to be predicted. 
newoffset 
A vector containing the values of the offset for the times for which the series is to be predicted. 
newm 
A vector containing the number of trials when forecasting binomial or binary time series. Defaults to the binary case. 
... 
Further arguments for the call, currently unused. 
Only one forecasting method is currently provided, for objects of class "glarma". This produces an object of class "glarmaForecast".
When forecasting one step ahead, the values in the matrix
newdata
(and offset
if there is an offset) in the GLARMA
model are used along with the regression coefficients in the model to
obtain the predicted value of eta, the regression
component of the state variable W. The predicted value of the
ARMA component of the state variable is then added to this value to
give the predicted value of W.
When further predictions are required, since no data is available to
calculate the predicted value of the state variable, an observation is
generated from the predicted distribution and the methodology for one
step ahead is then used on this generated data. This process is
repeated until predictions are obtained for the required number of
time periods (specified by n.ahead
). Note that the value of
n.ahead
must equal the row dimension of newdata
and if
they are specified, of newoffset
and newm
.
For completeness a randomly generated value of the time series is produced even for one stepahead prediction.
Note that the forecasted time series returned as the component
fitted
is then a randomly generated sample path for the
predicted time series. If a sample of such paths is produced by
repeated calls to forecast
then sample predicted distributions
can be obtained for the forecast series.
In the case of binary or binomial time series in addition to values of
the predictors in the regression component of the state variable and
the values of any offset, the numbers of trials for the binomially
distributed future observations are required. This information should
be provided in the argument newm
. If not, the number of trials
defaults to 1, which is the case of binary responses.
forecast
currently has no default method.
When object
is of class "glarma"
, forecast
returns an object of class "glarmaForecast"
with components:
eta 
the forecast values of the regression component of the state variable 
W 
the forecast values of the state variable W 
mu 
the conditional mean mu_t 
Y 
the simulated series based on the fitted model 
n.ahead 
the number of steps ahead for which the forecasts were
requested in the call to 
newdata 
the model matrix X comprising the values of the predictors for the times for which the series is to be predicted 
newoffset 
the vector containing the values of the offset for the times for which the series is to be predicted 
newm 
the vector giving the number of trials when forecasting binomial or binary time series 
model 
the 
"William T.M. Dunsmuir" <[email protected]> and "David J Scott" <[email protected]>
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118  require(zoo)
### Model number of deaths
data(DriverDeaths)
y < DriverDeaths[, "Deaths"]
X < as.matrix(DriverDeaths[, 2:5])
Population < DriverDeaths[, "Population"]
### Offset included
glarmamod < glarma(y, X, offset = log(Population/100000),
phiLags = c(12),
thetaLags = c(1),
type = "Poi", method = "FS",
residuals = "Pearson", maxit = 100, grad = 1e6)
print(summary(glarmamod))
XT1 < matrix(X[72,], nrow = 1)
offsetT1 < log(Population/100000)[72]
mu < forecast(glarmamod, 1, XT1, offsetT1)$mu
print(mu)
### Save some values
allX < X
allFits < fitted(glarmamod)
ally < y
### Look at a succession of forecasts
### Using actual values in forecasts
forecasts < numeric(72)
for (i in (62:71)){
y < DriverDeaths[1:i, "Deaths"]
X < as.matrix(DriverDeaths[1:i, 2:5])
Population < DriverDeaths[1:i, "Population"]
## Offset included
glarmamod < glarma(y, X, offset = log(Population/100000),
phiLags = c(12),
thetaLags = c(1),
type = "Poi", method = "FS",
residuals = "Pearson", maxit = 100, grad = 1e6)
XT1 < matrix(allX[i + 1, ], nrow = 1)
offsetT1 < log(DriverDeaths$Population[i + 1]/100000)
mu < forecast(glarmamod, 1, XT1, offsetT1)$mu
if (i == 62){
forecasts[1:62] < fitted(glarmamod)
}
forecasts[i+1] < mu
}
par(mfrow = c(1,1))
forecasts < ts(forecasts[63:72], start = c(1985, 10), deltat = 1/12)
fitted < ts(allFits, start = c(1980, 8), deltat = 1/12)
obs < ts(DriverDeaths$Deaths, start = c(1980, 8), deltat = 1/12)
plot(obs, ylab = "Driver Deaths", lty = 2,
main = "Single Vehicle Nighttime Driver Deaths in Utah")
points(obs)
lines(fitted, lwd = 2)
lines(forecasts, col = "red")
par(xpd = NA)
graph.param <
legend("top",
legend = c("observations",expression(estimated~mu[t]),
expression(predicted~mu[t])),
ncol = 3,
cex = 0.7,
bty = "n", plot = FALSE)
legend(graph.param$rect$left,
graph.param$rect$top + graph.param$rect$h,
legend = c("observations", expression(estimated~mu[t]),
expression(predicted~mu[t])),
col = c("black","black","red"),
lwd = c(1,2,1), lty = c(2,1,1),
pch = c(1, NA_integer_, NA_integer_),
ncol = 3,
cex = 0.7,
bty = "n",
text.font = 4)
par(xpd = FALSE)
### Generate a sample of Y values 2 steps ahead and examine the distribution
data(DriverDeaths)
y < DriverDeaths[, "Deaths"]
X < as.matrix(DriverDeaths[, 2:5])
Population < DriverDeaths[, "Population"]
### Fit the glarma model to the first 70 observations
glarmamod < glarma(y[1:70], X[1:70, ],
offset = log(Population/100000)[1:70],
phiLags = c(12),
thetaLags = c(1),
type = "Poi", method = "FS",
residuals = "Pearson", maxit = 100, grad = 1e6)
nObs < NROW(X)
n.ahead < 2
### Specify the X matrix and offset for the times where predictions
### are required
XT1 < as.matrix(X[(nObs  n.ahead + 1):nObs, ])
offsetT1 < log(Population/100000)[(nObs  n.ahead + 1):nObs]
nSims < 500
forecastY < matrix(ncol = n.ahead, nrow = nSims)
forecastMu < matrix(ncol = n.ahead, nrow = nSims)
### Generate sample predicted values
for(i in 1:nSims){
temp < forecast(glarmamod, n.ahead, XT1, offsetT1)
forecastY[i, ] < temp$Y
forecastMu[i, ] < temp$mu
}
### Examine distribution of sample of Y values n.ahead
table(forecastY[, 2])
par(mfrow = c(2,1))
barplot(table(forecastY[, 2]),
main = "Barplot of Sample Y Values 2 Steps Ahead")
hist(forecastY[, 2], xlab = "Sample Y values",
breaks=seq(0,max(forecastY[, 2])),
main = "Histogram of Sample Y Values 2 Steps Ahead\nwith 0.025 and 0.975 Quantiles")
abline(v = quantile(forecastY[, 2], c(0.025, 0.975)), col = "red")

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.