VGAM (version 1.0-1)

AR1: 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

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

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 paramete

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 (th

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 set

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.

Details

The AR-1 model implemented here has $$Y_1 \sim 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$.

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