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 Newton-Raphson
iteration.
For Poisson and negative binomial response distributions the log link is currently used. For binomial responses the logit link is currently used.
glarma(y, X, offset = NULL, type = "Poi", method = "FS", residuals = "Pearson", phiLags, thetaLags, phiInit, thetaInit, beta, alphaInit, alpha = 1, maxit = 30, grad = 2.22e-16)
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")beta of
the intercept.NULL or a numeric vector of length equal
to the number of cases. Used to specify an a priori known
component to be included in the linear predictor during fitting.beta, the AR terms and the MA terms.glm.nb.glarma for negative binomial counts."Poi" (Poisson), "Bin" (binomial) and "NegBin"
(negative binomial). The default is the
Poisson distribution."FS" (Fisher scoring), and "NR"
(Newton-Raphson). The default is to use Fisher scoring to estimate
the parameters of a GLARMA model."Pearson" and "Score", and for the binomial
distribution "Identity" is also allowed. The default is to
use Pearson residuals.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:beta, AR and
MA.phiLags plus the number of
thetaLags.NULL if there is no offset. 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 = X_t^T\beta + \mbox{offset} + Z_t.$$
The infinite moving average from the linear predictor is $$Z_t = \sum_{i=1}^\infty \gamma_i e_{t-i}.$$
This infinite moving average, is computed using the autoregressive moving average recursions
$$Z_t = \phi_1 (Z_{t-1} + e_{t-1}) + ... +% \phi_p (Z_{t-p} + e_{t-p}) + \theta_1 e_{t-1}% + ... + \theta_q e_{t-q}$$
where $p$ and $q$ are the orders of $phi$
and $theta$ respectively and the non-zero 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
$$\hat{\delta}^{(k+1)}=\hat{\delta}^{(k)}+\Omega(\hat{\delta}^{(k)}% )\frac{\partial l(\hat{\delta}^{(k)})}{\partial\delta}$$
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, Newton-Raphson 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 non-identifiable 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.
Asthma,
OxBoatRace, RobberyConvict, and
DriverDeaths.
### 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 = 1e-6)
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 = 1e-6)
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 = 1e-6)
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 = 1e-6)
glarmamod
summary(glarmamod)
Run the code above in your browser using DataLab