AR1 | R Documentation |

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

AR1(ldrift = "identitylink", lsd = "loglink", lvar = "loglink", lrho = "rhobitlink", idrift = NULL, isd = NULL, ivar = NULL, irho = NULL, imethod = 1, ishrinkage = 0.95, type.likelihood = c("exact", "conditional"), type.EIM = c("exact", "approximate"), var.arg = FALSE, nodrift = FALSE, print.EIM = FALSE, zero = c(if (var.arg) "var" else "sd", "rho"))

`ldrift, lsd, lvar, lrho` |
Link functions applied to the scaled mean, standard deviation
or variance, and correlation parameters.
The parameter |

`idrift, isd, ivar, irho` |
Optional initial values for the parameters.
If failure to converge occurs then try different values
and monitor convergence by using |

`ishrinkage, imethod, zero` |
See |

`var.arg` |
Same meaning as |

`nodrift` |
Logical, for determining whether to estimate the drift parameter.
The default is to estimate it.
If |

`type.EIM` |
What type of expected information matrix (EIM) is used in
Fisher scoring. By default, this family function calls
If |

`print.EIM` |
Logical. If |

`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., |

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 weak stationarity might be to verify whether
*1/rho* lies outside the unit circle.
(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*.

An object of class `"vglmff"`

(see `vglmff-class`

).
The object is used by modelling functions such as `vglm`

,
and `vgam`

.

Monitoring convergence is urged, i.e., set `trace = TRUE`

.

Moreover, if the exact EIMs are used, set `print.EIM = TRUE`

to compare the computed exact to the approximate EIM.

Under the VGLM/VGAM approach, parameters can be modelled in terms
of covariates. Particularly, if the standard deviation of
the white noise is modelled in this way, then
`type.EIM = "exact"`

may certainly lead to unstable
results. The reason is that white noise is a stationary
process, and consequently, its variance must remain as a constant.
Consequently, the use of variates to model
this parameter contradicts the assumption of
stationary random components to compute the exact EIMs proposed
by Porat and Friedlander (1987).

To prevent convergence issues in such cases, this family function
internally verifies whether the variance of the white noise remains
as a constant at each Fisher scoring iteration.
If this assumption is violated and `type.EIM = "exact"`

is set,
then `AR1`

automatically shifts to
`type.EIM = "approximate"`

.
Also, a warning is accordingly displayed.

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

Victor Miranda (exact method) and Thomas W. Yee (approximate method).

Porat, B. and Friedlander, B. (1987).
The Exact Cramer-Rao Bond for Gaussian Autoregressive Processes.
*IEEE Transactions on Aerospace and Electronic Systems*,
**AES-23(4)**, 537–542.

`AR1EIM`

,
`vglm.control`

,
`dAR1`

,
`arima.sim`

.

### Example 1: using arima.sim() to generate a 0-mean stationary time series. nn <- 500 tsdata <- data.frame(x2 = runif(nn)) ar.coef.1 <- rhobitlink(-1.55, inverse = TRUE) # Approx -0.65 ar.coef.2 <- rhobitlink( 1.0, inverse = TRUE) # Approx 0.50 set.seed(1) tsdata <- transform(tsdata, index = 1:nn, TS1 = arima.sim(nn, model = list(ar = ar.coef.1), sd = exp(1.5)), TS2 = arima.sim(nn, model = list(ar = ar.coef.2), sd = exp(1.0 + 1.5 * x2))) ### An autoregressive intercept--only model. ### ### Using the exact EIM, and "nodrift = TRUE" ### fit1a <- vglm(TS1 ~ 1, data = tsdata, trace = TRUE, AR1(var.arg = FALSE, nodrift = TRUE, type.EIM = "exact", print.EIM = FALSE), crit = "coefficients") Coef(fit1a) summary(fit1a) ## Not run: ### Two responses. Here, the white noise standard deviation of TS2 ### ### is modelled in terms of 'x2'. Also, 'type.EIM = exact'. ### fit1b <- vglm(cbind(TS1, TS2) ~ x2, AR1(zero = NULL, nodrift = TRUE, var.arg = FALSE, type.EIM = "exact"), constraints = list("(Intercept)" = diag(4), "x2" = rbind(0, 0, 1, 0)), data = tsdata, trace = TRUE, crit = "coefficients") coef(fit1b, matrix = TRUE) summary(fit1b) ### Example 2: another stationary time series nn <- 500 my.rho <- rhobitlink(1.0, inverse = TRUE) my.mu <- 1.0 my.sd <- exp(1) tsdata <- data.frame(index = 1:nn, TS3 = runif(nn)) set.seed(2) for (ii in 2:nn) tsdata$TS3[ii] <- my.mu/(1 - my.rho) + my.rho * tsdata$TS3[ii-1] + rnorm(1, sd = my.sd) tsdata <- tsdata[-(1:ceiling(nn/5)), ] # Remove the burn-in data: ### Fitting an AR(1). The exact EIMs are used. fit2a <- vglm(TS3 ~ 1, AR1(type.likelihood = "exact", # "conditional", type.EIM = "exact"), data = tsdata, trace = TRUE, crit = "coefficients") Coef(fit2a) summary(fit2a) # SEs are useful to know Coef(fit2a)["rho"] # Estimate of rho, for intercept-only models my.rho # The 'truth' (rho) Coef(fit2a)["drift"] # Estimate of drift, for intercept-only models my.mu /(1 - my.rho) # The 'truth' (drift) ## End(Not run)

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.