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,
    ishrinkage = 0.9, type.likelihood = c("exact", "conditional"),
    var.arg = FALSE, almost1 = 0.99, zero = c(-2, -3))

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, zero

See CommonVGAMffArguments for more information.

var.arg

Same meaning as uninormal.

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.

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  arimia.sim() to generate a stationary time series
nn <- 100; set.seed(1)
tsdata <- data.frame(x2 =  runif(nn))
tsdata  <- transform(tsdata,
              index = 1:nn,
              TS1 = arima.sim(nn, model = list(ar = -0.80),
                              sd = exp(1.0)),
              TS2 = arima.sim(nn, model = list(ar =  0.50),
                              sd = exp(1.0 + 2 * x2)))
fit1a <- vglm(cbind(TS1, TS2) ~ x2, AR1(zero = c(1:4, 6)),
             data = tsdata, trace = TRUE)
rhobit(-0.8)
rhobit( 0.5)
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 0.9-8 Index]