R version 4.5.0 beta (2025-03-27 r88065 ucrt) -- "How About a Twenty-Six"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require(glarma)
Loading required package: glarma
> require(zoo)
Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

> ### 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))

Call: glarma(y = y, X = X, offset = log(Population/1e+05), type = "Poi", 
    method = "FS", residuals = "Pearson", phiLags = c(12), thetaLags = c(1), 
    maxit = 100, grad = 1e-06)

Pearson Residuals:
    Min       1Q   Median       3Q      Max  
-1.7016  -0.7938   0.1189   0.5897   3.7199  

GLARMA Coefficients:
          Estimate  Std.Error z-ratio Pr(>|z|)
phi_12  -0.0338742  0.0899514  -0.377    0.706
theta_1 -0.0005256  0.0827289  -0.006    0.995

Linear Model Coefficients:
           Estimate Std.Error z-ratio Pr(>|z|)   
Intercept  -2.13559   0.80711  -2.646  0.00815 **
ReducedBAC -0.38582   0.15598  -2.474  0.01338 * 
FriSat      0.08663   0.09091   0.953  0.34060   
lnOMVDRate  0.55792   0.23779   2.346  0.01896 * 

    Null deviance: 95.499  on 71  degrees of freedom
Residual deviance: 67.560  on 66  degrees of freedom
AIC: 264.1927 

Number of Fisher Scoring iterations: 17

LRT and Wald Test:
Alternative hypothesis: model is a GLARMA process
Null hypothesis: model is a GLM with the same regression structure
          Statistic p-value
LR Test       0.160   0.923
Wald Test     0.142   0.932

> 
> XT1 <- matrix(X[72,], nrow = 1)
> offsetT1 <- log(Population/100000)[72]
> 
> mu <- forecast(glarmamod, 1, XT1, offsetT1)$mu
> print(mu)
[1] 1.933074
> 
> 
> ### 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)

Call: glarma(y = y, X = X, type = "NegBin", method = "FS", residuals = "Score", 
    thetaLags = c(1, 2, 5), maxit = 100, grad = 1e-06)

Score Residuals:
    Min       1Q   Median       3Q      Max  
-1.0000  -1.0000  -0.2070   0.6889   7.8095  

Negative Binomial Parameter:
      Estimate Std.Error z-ratio Pr(>|z|)   
alpha    2.856     1.039   2.748    0.006 **

GLARMA Coefficients:
        Estimate Std.Error z-ratio Pr(>|z|)    
theta_1  0.27621   0.05708   4.840 1.30e-06 ***
theta_2  0.23545   0.05812   4.051 5.09e-05 ***
theta_5 -0.00880   0.06265  -0.140    0.888    

Linear Model Coefficients:
              Estimate Std.Error z-ratio Pr(>|z|)    
Intcpt         0.09423   0.13827   0.681 0.495572    
Trend         -4.73820   2.65014  -1.788 0.073791 .  
CosAnnual     -0.04591   0.15980  -0.287 0.773869    
SinAnnual     -0.58596   0.17514  -3.346 0.000821 ***
CosSemiAnnual  0.29616   0.12207   2.426 0.015263 *  
SinSemiAnnual -0.26011   0.13196  -1.971 0.048705 *  

    Null deviance: 201.22  on 167  degrees of freedom
Residual deviance: 255.41  on 158  degrees of freedom
AIC: 506.134 

Number of Fisher Scoring iterations: 26

LRT and Wald Test:
Alternative hypothesis: model is a GLARMA process
Null hypothesis: model is a GLM with the same regression structure
          Statistic  p-value    
LR Test       21.52 8.20e-05 ***
Wald Test     31.07 8.21e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> 
> 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)

Call: glarma(y = Y, X = X, type = "Bin", method = "NR", residuals = "Pearson", 
    thetaLags = c(1, 2), maxit = 100, grad = 1e-06)

Pearson Residuals:
    Min       1Q   Median       3Q      Max  
-1.8865  -0.7968   0.4395   0.8436   2.9658  

GLARMA Coefficients:
        Estimate Std.Error z-ratio Pr(>|z|)    
theta_1   0.3396    0.1709   1.987 0.046975 *  
theta_2   0.5552    0.1459   3.804 0.000142 ***

Linear Model Coefficients:
             Estimate Std.Error z-ratio Pr(>|z|)   
(Intercept)  0.349954  0.267249   1.309  0.19038   
Diff         0.114755  0.038238   3.001  0.00269 **
I(Diff^2)   -0.011333  0.004916  -2.305  0.02114 * 

    Null deviance: 216.16  on 155  degrees of freedom
Residual deviance: 148.53  on 151  degrees of freedom
AIC: 193.0553 

Number of Newton Raphson iterations: 5

LRT and Wald Test:
Alternative hypothesis: model is a GLARMA process
Null hypothesis: model is a GLM with the same regression structure
          Statistic  p-value    
LR Test       15.56 0.000417 ***
Wald Test     17.51 0.000158 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> 
> 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)
[1] 0.3
> 
> ###############################
> ### Test multiple steps ahead
> ###############################
> ### 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)
> n.ahead <- 10
> 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))
> for (lead in 1:10) {
+     hist(forecastY[, lead], freq = FALSE, breaks = 0:15,
+          main = as.character(lead))
+     print(table(forecastY[, lead]))
+ }

  0   1   2   3   4   5   6 
112 172 134  56  19   4   3 

  0   1   2   3   4   5   6   7   8   9  10 
 34  66 125  99  82  56  21  10   4   2   1 

  0   1   2   3   4   5   6 
 73 149 146  77  36  15   4 

  0   1   2   3   4   5 
123 160 109  76  19  13 

  0   1   2   3   4 
189 181  86  39   5 

  0   1   2   3   4   5   6   7 
 72 163 130  78  33  16   7   1 

  0   1   2   3   4   5   6   8 
112 152 127  64  31  10   3   1 

  0   1   2   3   4   5   6   7 
 56 114 130 106  61  24   6   3 

  0   1   2   3   4   5   6   7 
 57 117 114 105  73  24   5   5 

  0   1   2   3   4   5   6   7 
 87 153 124  76  33  21   4   2 
> 
> ##########################
> ### 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)
> 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])

  0   1   2   3   4   5   6   7   8 
 77 139 133  87  45  12   5   1   1 
> 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")
> 
> proc.time()
   user  system elapsed 
   9.20    0.56   10.51