VGAM (version 0.8-7)

garma: GARMA (Generalized Autoregressive Moving-Average) Models

Description

Fits GARMA models to time series data.

Usage

garma(link = "identity", earg=list(), p.ar.lag = 1, q.ma.lag = 0,
      coefstart = NULL, step = 1)

Arguments

link

Link function applied to the mean response. The default is suitable for continuous responses. The link loge should be chosen if the data are counts. Links such as logi

earg

List. Extra argument for the link. See earg in Links for general information. In particular, this argument is useful when the log or logit link is chosen: for log and logit, zero values

p.ar.lag

A positive integer, the lag for the autoregressive component. Called $p$ below.

q.ma.lag

A non-negative integer, the lag for the moving-average component. Called $q$ below.

coefstart

Starting values for the coefficients. For technical reasons, the argument coefstart in vglm cannot be used.

step

Numeric. Step length, e.g., 0.5 means half-stepsizing.

Value

• An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm.

Warning

This VGAM family function is 'non-standard' in that the model does need some coercing to get it into the VGLM framework. Special code is required to get it running. A consequence is that some methods functions may give wrong results when applied to the fitted object.

Details

This function draws heavily on Benjamin et al. (1998). See also Benjamin et al. (2003). GARMA models extend the ARMA time series model to generalized responses in the exponential family, e.g., Poisson counts, binary responses. Currently, this function can handle continuous, count and binary responses only. The possible link functions given in the link argument reflect this, and the user must choose an appropriate link.

The GARMA($p, q$) model is defined by firstly having a response belonging to the exponential family $$f(y_t|D_t) = \exp \left{ \frac{y_t \theta_t - b(\theta_t)}{\phi / A_t} + c(y_t, \phi / A_t) \right}$$ where $\theta_t$ and $\phi$ are the canonical and scale parameters respectively, and $A_t$ are known prior weights. The mean $\mu_t=E(Y_t|D_t)=b'(\theta_t)$ is related to the linear predictor $\eta_t$ by the link function $g$. Here, $D_t={x_t,\ldots,x_1,y_{t-1},\ldots,y_1,\mu_{t-1},\ldots,\mu_1}$ is the previous information set. Secondly, the GARMA($p, q$) model is defined by $$g(\mu_t) = \eta_t = x_t^T \beta + \sum_{k=1}^p \phi_k (g(y_{t-k}) - x_{t-k}^T \beta) + \sum_{k=1}^q \theta_k (g(y_{t-k}) - \eta_{t-k}).$$ Parameter vectors $\beta$, $\phi$ and $\theta$ are estimated by maximum likelihood.

References

Benjamin, M. A., Rigby, R. A. and Stasinopoulos, M. D. (1998) Fitting Non-Gaussian Time Series Models. Pages 191--196 in: Proceedings in Computational Statistics COMPSTAT 1998 by Payne, R. and P. J. Green. Physica-Verlag.

Benjamin, M. A., Rigby, R. A. and Stasinopoulos, M. D. (2003) Generalized Autoregressive Moving Average Models. Journal of the American Statistical Association, 98: 214--223.

Zeger, S. L. and Qaqish, B. (1988) Markov regression models for time series: a quasi-likelihood approach. Biometrics, 44: 1019--1031.

See Also

The site http://www.stat.auckland.ac.nz/~yee contains more documentation about this family function.

Examples

gdata = data.frame(interspike = c(68, 41, 82, 66, 101, 66, 57,  41,  27, 78,
59, 73,  6, 44,  72, 66, 59,  60,  39, 52,
50, 29, 30, 56,  76, 55, 73, 104, 104, 52,
25, 33, 20, 60,  47,  6, 47,  22,  35, 30,
29, 58, 24, 34,  36, 34,  6,  19,  28, 16,
36, 33, 12, 26,  36, 39, 24,  14,  28, 13,
 2, 30, 18, 17,  28,  9, 28,  20,  17, 12,
19, 18, 14, 23,  18, 22, 18,  19,  26, 27,
23, 24, 35, 22,  29, 28, 17,  30,  34, 17,
20, 49, 29, 35,  49, 25, 55,  42,  29, 16)) # See Zeger and Qaqish (1988)
gdata = transform(gdata, spikenum = seq(interspike))
bvalue = 0.1  # .Machine$double.xmin # Boundary value
fit = vglm(interspike ~ 1, trace = TRUE, data = gdata,
           garma("loge", earg = list(bvalue = bvalue),
                 p = 2, coef = c(4, 0.3, 0.4)))
summary(fit)
coef(fit, matrix = TRUE)
Coef(fit)  # A bug here
with(gdata, plot(interspike, ylim = c(0, 120), las = 1,
     xlab = "Spike Number", ylab = "Inter-Spike Time (ms)", col = "blue"))
with(gdata, lines(spikenum[-(1:fit@misc$plag)], fitted(fit), col = "orange"))
abline(h = mean(with(gdata, interspike)), lty = "dashed", col = "gray")