AR1 {VGAM}R Documentation

Autoregressive Process with Order-1 Family Function

Description

Maximum likelihood estimation of the three-parameter AR-1 model

Usage

AR1(ldrift = "identitylink", lsd  = "loge", lvar = "loge",
    lrho = "rhobit", idrift  = NULL,
    isd  = NULL, ivar = NULL, irho = NULL, imethod = 1,
    ishrinkage = 1, type.likelihood = c("exact", "conditional"),
    var.arg = FALSE, nodrift = FALSE, almost1 = 0.99,
    zero = c(if (var.arg) "var" else "sd", "rho"))
 AR1.control(half.stepsizing = FALSE, ...)

Arguments

ldrift, lsd, lvar, lrho

Link functions applied to the scaled mean, standard deviation or variance, and correlation parameters. The parameter drift is known as the drift, and it is a scaled mean. See Links for more choices.

idrift, isd, ivar, irho

Optional initial values for the parameters. If failure to converge occurs then try different values and monitor convergence by using trace = TRUE. For a S-column response, these arguments can be of length S, and they are recycled by the columns first. A value NULL means an initial value for each response is computed internally.

ishrinkage, imethod, zero

See CommonVGAMffArguments for more information. The default for zero assumes there is a drift parameter to be estimated (the default for that argument), so if a drift parameter is suppressed and there are covariates, then zero will need to be assigned the value 1 or 2 or NULL.

var.arg

Same meaning as uninormal.

nodrift

Logical, for determining whether to estimate the drift parameter. The default is to estimate it. If TRUE, the drift parameter is set to 0 and not estimated.

type.likelihood

What type of likelihood function is maximized. The first choice (default) is the sum of the marginal likelihood and the conditional likelihood. Choosing the conditional likelihood means that the first observation is effectively ignored (this is handled internally by setting the value of the first prior weight to be some small positive number, e.g., 1.0e-6). See the note below.

almost1

A value close to 1 but slightly smaller. One of the off-diagonal EIM elements is multiplied by this, to ensure that the EIM is positive-definite.

half.stepsizing, ...

A logical value, overwriting that of vglm.control. Currently this setting is potentially dangerous, and is used for aesthetics at the solution—no jittering occurs. This can often be seen by setting trace = TRUE when the value is set to TRUE. The jittering is due to some heuristics applied to handle the first observation—either by setting its prior weight to a value very close to 0, else adjust for its EIM which is not of full rank.

Details

The AR-1 model implemented here has

Y(1) ~ N(mu, sigma^2 / (1-rho^2),

and

Y(i) = mu^* + rho * Y(i-1) + e(i)

where the e(i) are i.i.d. Normal(0, sd = sigma) random variates.

Here are a few notes: 1. A test for stationarity might be to test whether mu^* is intercept-only. 2. The mean of all the Y(i) is mu^* / (1-rho) and these are returned as the fitted values. 3. The correlation of all the Y(i) with Y(i-1) is rho. 4. The default link function ensures that -1 < rho < 1.

Value

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

Warning

Monitoring convergence is urged: set trace = TRUE.

Yet to do: add an argument that allows the scaled mean parameter to be deleted, i.e, a 2-parameter model is fitted. Yet to do: ARff(p.lag = 1) should hopefully be written soon.

Note

For type.likelihood = "conditional", the prior weight for the first observation is set to some small positive number, which has the effect of deleting that observation. However, n is still the original n so that statistics such as the residual degrees of freedom are unchanged (uncorrected possibly).

Multiple responses are handled. The mean is returned as the fitted values.

Practical experience has shown that half-stepping is a very good idea. The default options use step sizes that are about one third the usual step size. Consequently, maxit is increased to about 100, by default.

Author(s)

Thomas W. Yee and Victor Miranda

See Also

vglm.control, dAR1, uninormal, arima.sim,

Examples

# Example 1: using  arima.sim() to generate a stationary time series
nn <- 1000; set.seed(1)
tsdata <- data.frame(x2 =  runif(nn))
ar.coef.1 <- rhobit(-2, inverse = TRUE)  # Approx -0.8
ar.coef.2 <- rhobit( 1, inverse = TRUE)  # Approx  0.5
tsdata  <- transform(tsdata,
              index = 1:nn,
              TS1 = arima.sim(nn, model = list(ar = ar.coef.1),
                              sd = exp(1.0)),
              TS2 = arima.sim(nn, model = list(ar = ar.coef.2),
                              sd = exp(1.0 + 2 * x2)))
fit1a <- vglm(cbind(TS1, TS2) ~ x2, AR1(zero = "rho", nodrift = TRUE),
              data = tsdata, trace = TRUE)
coef(fit1a, matrix = TRUE)
summary(fit1a)  # SEs are useful to know

# Example 2: another stationary time series
nn <- 1000
my.rho <- rhobit(-1.0, inverse = TRUE)
my.mu  <- 2.5
my.sd  <- exp(1)
tsdata  <- data.frame(index = 1:nn, TS3 = runif(nn))
for (ii in 2:nn)
  tsdata$TS3[ii] <- my.mu + my.rho * tsdata$TS3[ii-1] + rnorm(1, sd = my.sd)
tsdata <- tsdata[-(1:ceiling(nn/5)), ]  # Remove the burn-in data:
fit2a <- vglm(TS3 ~ 1, AR1(type.likelihood = "conditional"),
             data = tsdata, trace = TRUE)
coef(fit2a, matrix = TRUE)
summary(fit2a)  # SEs are useful to know
Coef(fit2a)["rho"]  # Estimate of rho for intercept-only models
my.rho
coef(fit2a)[1]  # drift
my.mu           # Should be the same
head(weights(fit2a, type= "prior"))    # First one is effectively deleted
head(weights(fit2a, type= "working"))  # Ditto

[Package VGAM version 1.0-1 Index]