The function glarma
is used to fit generalized linear
autoregressive moving average models with various distributions
(Poisson, binomial, negative binomial) using either Pearson residuals
or score residuals, and for the binomial distribution, identity
residuals. It also estimates the parameters of the GLARMA model with
various distributions by using Fisher scoring or NewtonRaphson
iteration.
For Poisson and negative binomial response distributions the log link is currently used. For binomial responses the logit link is currently used.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  glarma(y, X, offset = NULL, type = "Poi", method = "FS", residuals = "Pearson",
phiLags, thetaLags, phiInit, thetaInit, beta, alphaInit,
alpha = 1, maxit = 30, grad = 2.22e16)
glarmaPoissonPearson(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS")
glarmaPoissonScore(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS")
glarmaBinomialIdentity(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS")
glarmaBinomialPearson(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS")
glarmaBinomialScore(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS")
glarmaNegBinPearson(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS")
glarmaNegBinScore(y, X, offset = NULL, delta, phiLags, thetaLags,
method = "FS")

y 
Numeric vector; the response variable. If the response variable is for the model with the binomial distribution, it should be a n by 2 matrix, one column is the number of successes and another is the number of failures. 
X 
Matrix; the explanatory variables. A vector of ones should be
added to the data matrix as the first column for the 
offset 
Either 
beta 
Numeric vector; initial values of the regression coefficients. 
phiLags 
Numeric vector; AR orders. 
phiInit 
Numeric vector; initial values for the corresponding AR orders. 
thetaLags 
Numeric vector; MA orders. 
thetaInit 
Numeric vector; initial values for the corresponding MA orders. 
delta 
Numeric vector; initial values of the parameters for the
GLARMA estimation procedure. It is a combination of the parameters of

alpha 
Numeric; an optional initial shape parameter for

alphaInit 
Numeric; an initial shape parameter for

type 
Character; the count distribution. Possible values are

method 
Character; method of iteration to be used. Possible
values are 
residuals 
Character; the type of residuals to be used. Possible
values are 
maxit 
Numeric; the maximum number of iterations allowed. 
grad 
Numeric; the tolerance for recognizing numbers, which are smaller than the specified tolerance, as zero. 
Models for glarma
are specified symbolically. A typical model
has the form y
(response), X
(terms) where y
is
the count or factor reponse vector, X
is a series of terms
which specifies a linear predictor for the response. It should be noted
that the first column of X
should be a vector of 1s as the
intercept in the model. Four initial parameters that need to be
estimated are combined into delta = (beta, phi, theta, alpha), where alpha
is an optional parameter to accomodate the negative binomial
model. Note that in the function glm.nb
from the
package MASS, this parameter is called theta
.
For Poisson and negative binomial response distributions the log link is currently used. For binomial responses the logit link is currently used.
The generalized linear autoregressive moving average models are computed as follows.
The linear predictor for the response is
log(mu_t) = W_t = transpose(X_t) * beta + offset + Z_t.
The infinite moving average from the linear predictor is
Z_t = sum(gamma_i * residuals_(ti)).
This infinite moving average, is computed using the autoregressive moving average recursions
Z_t = phi_1 * (Z_(t1) + e_(t1)) + ... + phi_p * (Z_{tp} + e_(tp)) + theta_1 * e_{1} + ... + theta_q * e_{tq}
where p and q are the orders of phi
and theta respectively and the nonzero lags of the vectors
phi
and theta
may be specified by the user via the
arguments phiLag
and thetaLag
.
There are two types of residuals which may be used in each recursion, Pearson residuals or score residuals, and in addition, for the binomial distribution, identity residuals may be used. The infinite moving average, Z_t, depends on the type of residuals used, as do the final parameters obtained from the filter. Standardisation of past observed counts is necessary to avoid instability, therefore the user should choose the appropriate type of residuals depending on the situation.
The method of estimation for parameters implemented in the function aims to maximise the log likelihood by an iterative method commencing from suitably chosen initial values for the parameters. Starting from initial values delta hat^(0) for the vector of parameters updates are obtained using the iterations
delta_(k+1) = delta_(k) + Omega(delta_k) * first derivative of log(delta_k)
where Omega(delta hat ^(k)) is some suitably chosen matrix.
Iterations continue for k >= 1 until convergence is reached or the number of iterations k reaches a user specified upper limit on maximum iterations in which case they will stop. The convergence criterion used in our implementation is that based on eta, the maximum of absolute values of the first derivatives.
When eta is less than a user specified value grad
the iterations stop. There are two methods of optimization of the
likelihood, NewtonRaphson and Fisher scoring. The method used is
specified by the argument method
. It should be noticed that if
the initial value for parameters are not chosen well, the
optimization of the likelihood might fail to converge. Care is needed
when fitting mixed ARMA specifications because there is potential for
the AR and MA parameters to be nonidentifiable if the orders p and
q are too large. Lack of identifiability manifests itself in the
algorithm to optimize the likelihood failing to converge and/or the
hessian being singularâ€”check the warning messages and convergence
error codes.
The function summary
(i.e., summary.glarma
)
can be used to obtain or print a summary of the results.
The generic accessor functions coef
(i.e.,
coef.glarma
), logLik
(i.e.,
logLik.glarma
), fitted
(i.e.,
fitted.glarma
), residuals
(i.e.,
residuals.glarma
), nobs
(i.e.,
nobs.glarma
), model.frame
(i.e.,
model.frame.glarma
) and extractAIC
(i.e.,
extractAIC.glarma
) can be used to extract various useful
features of the value returned by glarma
.
glarma
returns an object of class "glarma" with components:
delta 
a vector of coefficients for 
logLik 
the loglikelihood of the specific distribution. 
logLikDeriv 
the derivative of the loglikelhood of the specified distribution. 
logLikDeriv2 
the second derivative of the loglikelihood of the specified distribution. 
eta 
the estimated linear predictor. 
mu 
the GLARMA estimated mean. 
fitted.values 
the GLARMA fitted values. 
residuals 
the residuals of the type specified. 
cov 
the estimated covariance matrix of the maximum likelihood estimators. 
phiLags 
vector of AR orders. 
thetaLags 
vector of MA orders. 
r 
the number of columns in the model matrix. 
pq 
the number of 
null.deviance 
the deviance from the initial GLM fit. 
df.null 
the degrees of freedom from the initial GLM fit. 
y 
the y vector used in the GLARMA model. 
X 
the model matrix. 
offset 
the offset, 
type 
the distribution of the counts. 
method 
the method of iteration used. 
residType 
the type of the residuals returned. 
call 
the matched call. 
iter 
the number of iterations. 
errCode 
the error code; 0 indicating successful convergence of the iteration method, 1 indicating failure. 
WError 
error code for finiteness of W; 0 indicating all values of W are finite, 1 indicating at least one infinite value. 
min 
the minimum of the absolute value of the gradient. 
aic 
A version of Akaike's An Information Criterion, minus twice the maximized loglikelihood plus twice the number of parameters, computed by the aic component of the family. For binomial and Poisson families the dispersion is fixed at one and the number of parameters is the number of coefficients. 
The original GLARMA routine for Poisson responses was developed in collaboration with Richard A. Davis and Ying Wang. The binomial response version was developed with the assistance of Haolan Lu. The extension to negative binomial response was carried out by Bo Wang. Daniel Drescher contributed to the initial structure of the software used as the basis of the package.
The main author of the package is "William T.M. Dunsmuir" <w.dunsmuir@unsw.edu.au>. Package development was carried out by Cenanning Li supervised by David J. Scott.
Dunsmuir, William T. M. and Scott, David J. (2015) The glarma Package for ObservationDriven Time Series Regression of Counts. Journal of Statistical Software, 67(7), 1â€“36. http://dx.doi.org/10.18637/jss.v067.i07
Additional examples may be found in Asthma
,
OxBoatRace
, RobberyConvict
, and
DriverDeaths
.
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  ### Example from Davis, Dunsmuir Wang (1999)
## MA(1,2,5), Pearson Residuals, Fisher Scoring
data(Polio)
y < Polio[, 2]
X < as.matrix(Polio[, 3:8])
glarmamod < glarma(y, X, thetaLags = c(1,2,5), type = "Poi", method = "FS",
residuals = "Pearson", maxit = 100, grad = 1e6)
glarmamod
summary(glarmamod)
## Score Type (GAS) Residuals, Fisher Scoring
glarmamod < glarma(y, X, thetaLags = c(1,2,5), type = "Poi", method = "FS",
residuals = "Score", maxit = 100, grad = 1e6)
glarmamod
summary(glarmamod)
## Score Type (GAS) Residuals, Newton Raphson
## Note: Newton Raphson fails to converge from GLM initial estimates.
## Setting up the initial estimates by ourselves
init.delta < glarmamod$delta
beta < init.delta[1:6]
thetaInit < init.delta[7:9]
glarmamod < glarma(y, X, beta = beta, thetaLags = c(1, 2, 5),
thetaInit = thetaInit, type ="Poi", method = "NR",
residuals = "Score", maxit = 100, grad = 1e6)
glarmamod
summary(glarmamod)
## AR(1,5), Pearson Residuals, Fisher Scoring
glarmamod < glarma(y, X, phiLags = c(1, 5), type = "Poi", method = "FS",
residuals = "Pearson", maxit = 100, grad = 1e6)
glarmamod
summary(glarmamod)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.