# tests/forecastTests.R In glarma: Generalized Linear Autoregressive Moving Average Models

require(glarma)
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 = 1e-6)
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 = 1e-6)
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)

### Polio:
data(Polio)
y  <-  Polio[, 2]
X  <-  as.matrix(Polio[, 3:8])
glarmamod <- glarma(y, X, thetaLags = c(1,2,5),
type = "NegBin", method = "FS",
residuals = "Score", maxit = 100, grad = 1e-6)

summary(glarmamod)

XT1 <- matrix(X[168, ], nrow = 1)
mu <- forecast(glarmamod, 1, XT1, offsetT1 = 0)\$mu

plot(0:11, dnbinom(0:11, mu = mu, size =  coef(glarmamod)\$NB))

### Save some values
allX <- X
allFits <- fitted(glarmamod)
ally <- y

### Look at a succession of forecasts
### Using actual values in forecasts
forecasts <- numeric(168)
for (i in (156:167)){
y  <-  Polio[1:i, 2]
X  <-  as.matrix(Polio[1:i, 3:8])
glarmamod <- glarma(y, X, thetaLags = c(1,2,5),
type = "NegBin", method = "FS",
residuals = "Score", maxit = 100, grad = 1e-6)
XT1 <- matrix(c(1, (i + 1 - 72)/1000, cos(2*pi*(i + 1)/12),
sin(2*pi*(i + 1)/12), cos(4*pi*(i + 1)/12),
sin(4*pi*(i + 1)/12)), nrow = 1)
offsetT1 <- 0
mu <- forecast(glarmamod, 1, XT1, offsetT1)\$mu
if (i == 156){
forecasts[1:156] <- fitted(glarmamod)
}
forecasts[i+1] <- mu
}
par(mfrow = c(1,1))
forecasts <- ts(forecasts[157:168], start = c(1983, 1), deltat = 1/12)
fitted <- ts(allFits, start = c(1970, 1), deltat = 1/12)
data <- ts(Polio\$Cases, start = c(1970, 1), deltat = 1/12)
plot(data, ylab = "Polio Cases",
main = "Monthly Numbers of Polio Cases in the USA")
points(data)
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])),
lty = c(2,1,1),
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)

### Boat Race:
data(OxBoatRace)

y1 <- OxBoatRace\$Camwin
n1 <- rep(1, length(OxBoatRace\$Year))
Y <- cbind(y1, n1 - y1)
X <- cbind(OxBoatRace\$Intercept, OxBoatRace\$Diff)
colnames(X) <- c("Intercept", "Weight Diff")

oxcamglm <- glm(Y ~ Diff + I(Diff^2),
data = OxBoatRace,
family = binomial(link = "logit"), x = TRUE)

X <- oxcamglm\$x

glarmamod <- glarma(Y, X, thetaLags = c(1, 2), type = "Bin", method = "NR",
residuals = "Pearson", maxit = 100, grad = 1e-6)

summary(glarmamod)

par(mfrow=c(3,2))
plot(glarmamod)

### Save some values
allX <- X
allFits <- fitted(glarmamod)
ally1 <- y1

### Look at a succession of forecasts
### Using actual values in forecasts
forecasts <- numeric(156)
for (i in (146:155)){
y1 <- OxBoatRace\$Camwin[1:i]
n1 <- rep(1, i)
Y <- cbind(y1, n1 - y1)
X <- cbind(OxBoatRace\$Intercept[1:i], OxBoatRace\$Diff[1:i])
colnames(X) <- c("Intercept", "Weight Diff")

oxcamglm <- glm(Y ~ Diff + I(Diff^2),
data = OxBoatRace[1:i, ],
family = binomial(link = "logit"), x = TRUE)
X <- oxcamglm\$x
glarmamod <- glarma(Y, X, thetaLags = c(1, 2), type = "Bin", method = "NR",
residuals = "Pearson", maxit = 100, grad = 1e-6)
XT1 <- matrix(allX[i + 1, ], nrow = 1)
offsetT1 <- 0
mu <- forecast(glarmamod, 1, XT1, offsetT1 = 0)\$mu
if (i == 146){
forecasts[1:146] <- fitted(glarmamod)
}
forecasts[i+1] <- mu
}

par(mfrow = c(1,1))
forecasts <- zoo(forecasts, order.by = OxBoatRace\$Year)
fitted <- zoo(allFits, order.by = OxBoatRace\$Year)
plot(fitted, lwd = 2, xlab = "Year", ylab = "Cambridge Win")
lines(OxBoatRace\$Year[147:156], forecasts[147:156], col = "red")
plot(fitted[137:156], lwd = 2, ylim = c(0,1),
xlab = "Year", ylab = "Cambridge Win",
main = "Oxford versus Cambridge Boat Race")
lines(OxBoatRace\$Year[147:156], forecasts[147:156], col = "red")
points(OxBoatRace\$Year[137:156],OxBoatRace\$Camwin[137:156])
abline(h=0.5)
par(xpd = NA)
graph.param <-
legend("top",
legend = c("observations",expression(estimated~mu[t]),
expression(predicted~mu[t])),
lty = c(2,1,1),
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(NA,2,1), lty = c(NA,1,1),
pch = c(1, NA_integer_, NA_integer_),
ncol = 3,
cex = 0.7,
bty = "n",
text.font = 4)
par(xpd = FALSE)

### is there skill in the GLARMA model?
results <- cbind(round(forecasts[147:156]),OxBoatRace\$Camwin[147:156])
correct <- (results[,1] - results[,2] == 0)
mean(correct)

###############################
###############################
### 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 = 1e-6)
XT1 <- as.matrix(X[61:72, ])
offsetT1 <- log(Population/100000)[61:72]

nObs <- NROW(X)
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)

for(i in 1:nSims){
temp <-  forecast(glarmamod, n.ahead, XT1, offsetT1)
forecastY[i, ] <- temp\$Y
forecastMu[i, ] <- temp\$mu
}
par(mfrow = c(4, 3))
hist(forecastY[, lead], freq = FALSE, breaks = 0:15,
}

##########################
### Multiple steps ahead example for man page
### 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 = 1e-6)

nObs <- NROW(X)
### 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",
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")

## Try the glarma package in your browser

Any scripts or data that you put into this service are public.

glarma documentation built on Feb. 9, 2018, 6:08 a.m.