ARXff: VGLTSMs family functions: Order--p Autoregressive Model with covariates

Description

Maximum likelihood estimation of the order--p autoregressive model (AR(p)) with covariates. Estimates the drift, standard deviation (or variance) of the random noise (not necessarily constant), and coefficients of the conditional--mean model.

Usage

ARXff(order    = 1,
            zero     = c(if (nodrift) NULL else "ARdrift", "ARcoeff"),
            xLag     = 0,
            type.EIM = c("exact", "approximate")[1],
            var.arg  = TRUE, 
            nodrift  = FALSE,
            noChecks = FALSE,
            ldrift   = "identitylink", 
            lsd      = "loglink",
            lvar     = "loglink",
            lARcoeff = "identitylink",
            idrift   = NULL,
            isd      = NULL,
            ivar     = NULL,
            iARcoeff = NULL)

Value

An object of class "vglmff"

(see vglmff-class). The object is used by VGLM/VGAM modelling functions, such as vglm or vgam.

Arguments

order

The order (i.e., 'p') of the AR model, which is recycled if needed. See below for further details. By default, an autoregressive model of order-$1$ is fitted.

zero

Integer or character--strings vector. Name(s) or position(s) of the parameters/linear predictors to be modeled as intercept-only. Details at zero.

xLag

Same as ARIMAXff.

type.EIM, var.arg, nodrift, noChecks

Same as ARIMAXff.

ldrift, lsd, lvar, lARcoeff

Link functions applied to the drift, the standar deviation (or variance) of the noise, and the AR coefficients. Same as ARIMAXff.

Further details on CommonVGAMffArguments.

idrift, isd, ivar, iARcoeff

Same as ARIMAXff.

Author

Victor Miranda and T. W. Yee

Warning

Values of the estimates may not correspond to stationary ARs, leading to low accuracy in the MLE estimates, e.g., values very close to 1.0. Stationarity is then examined, via checkTS.VGAMextra, if noChecks = FALSE (default) and no constraint matrices are set (See constraints for further details on this). If the estimated model very close to be non-stationary, then a warning will be outlined. Set noChecks = TRUE to completely ignore this.

NOTE: Full details on these 'checks' are shown within the summary() output.

Details

This family function describes an autoregressive model of order-$p$ with covariates (ARX(p)). It is a special case of the subclass VGLM--ARIMA (Miranda and Yee, 2018):

$$ Y_t | \Phi_{t - 1} = \mu_t + \theta_{1} Y_{t - 1} + \ldots + \theta_p Y_{t - p} + \varepsilon_t,$$

where $\boldsymbol{x}_t$ a (possibly time--varying) covariate vector and $\mu_t = \mu^{\star} + \boldsymbol{\beta}^T \boldsymbol{x}_t$ is a (time--dependent) scaled--mean, known as drift.

At this stage, conditional Gaussian white noise, $\varepsilon_t| \Phi_{t - 1}$ is handled, in the form $$\varepsilon_t | \Phi_{t - 1} \sim N(0, \sigma^2_{\varepsilon_t | \Phi_{t - 1}}).$$ The distributional assumptions on the observations are then

$$Y_t | \Phi_{t - 1} \sim N(\mu_{t | \Phi_{t - 1}}, \sigma^2_{\varepsilon_t | \Phi_{t - 1}}), $$ involving the conditional mean equation for the ARX(p) model: $$\mu_{t | \Phi_{t - 1}} = \mu_t + \boldsymbol{\beta}^T * \boldsymbol{x}_t \theta_{1} Y_{t - 1} + \ldots + \theta_p Y_{t - p}.$$

$\Phi_{t}$ denotes the information of the joint process $\left(Y_{t}, \boldsymbol{x}_{t + 1}^T \right)$, at time $t$.

The loglikelihood is computed by dARp, at each Fisher scoring iteration.

The linear predictor is $$\boldsymbol{\eta} = \left( \mu_t, \log \sigma^{2}_{\varepsilon_{t | \Phi_{t - 1}}}, \theta_1, \ldots, \theta_p \right)^T.$$

Note, the covariates may also intervene in the conditional variance model $\log \sigma^{2}_{\varepsilon_{t | \Phi_{t - 1}}}.$ Hence, this family function does not restrict the noise to be strictly white noise (in the sense of constant variance).

The unconditional mean, $ E(Y_{t}) = \mu$, satisfies $$\mu \rightarrow \frac{\mu^{\star}}{1 - (\theta_1 + \ldots + \theta_p)} $$ when the process is stationary, and no covariates are involved.

This family function currently handles multiple responses so that a matrix can be used as the response. Also, for further details on VGLM/VGAM--link functions refer to Links.

Further choices for the random noise, besides Gaussian, will be implemented over time.

References

Madsen, H. (2008) Time Series Analysis. Florida, USA: Chapman & Hall(Sections 5.3 and 5.5).

Porat, B., and Friedlander, B. (1986) Computation of the Exact Information Matrix of Gaussian Time Series with Stationary Random Components. IEEE Transactions on Acoustics, Speech and Signal Processing. ASSp-34(1), 118--130.

Examples

Run this code

set.seed(1)
nn     <- 150
tsdata <- data.frame(x2 =  runif(nn))             # A single covariate.
theta1 <- 0.45; theta2 <- 0.31; theta3 <- 0.10     # Coefficients
drift  <- c(1.3, -1.1)                             # Two responses.
sdAR   <- c(sqrt(4.5), sqrt(6.0))                  # Two responses.

# Generate AR sequences of order 2 and 3, under Gaussian noise.
# Note, the drift for 'TS2' depends on x2 !
tsdata  <-  data.frame(tsdata, TS1 = arima.sim(nn, 
              model = list(ar = c(theta1, theta1^2)),  rand.gen = rnorm, 
              mean = drift[1], sd = sdAR[1]),
                                TS2 = arima.sim(nn,  
              model = list(ar = c(theta1, theta2, theta3)), rand.gen = rnorm, 
              mean = drift[2] + tsdata$x2 , sd = sdAR[2]))

# EXAMPLE 1. A simple AR(2), maximizing the exact log-likelihood
# Note that parameter constraints are involved for TS1, but not 
# considered in this fit. "rhobitlink" is used as link for AR coeffs.

fit.Ex1 <- vglm(TS1 ~ 1, ARXff(order = 2, type.EIM = "exact",
                      #iARcoeff = c(0.3, 0.3, 0.3), # OPTIONAL INITIAL VALUES
                      # idrift = 1, ivar = 1.5, isd = sqrt(1.5),
                      lARcoeff = "rhobitlink"), 
              data = tsdata,  trace = TRUE, crit = "loglikelihood")
Coef(fit.Ex1)
summary(fit.Ex1)
vcov(fit.Ex1, untransform = TRUE)       # Conformable with this fit.
AIC(fit.Ex1)
#------------------------------------------------------------------------#
# Fitting same model using arima(). 
#------------------------------------------------------------------------#
(fitArima <- arima(tsdata$TS1, order = c(2, 0, 0)))
# Compare with 'fit.AR'. True are theta1 = 0.45; theta1^2 = 0.2025
Coef(fit.Ex1)[c(3, 4, 2)]    # Coefficients estimated in 'fit.AR'


# EXAMPLE 2. An AR(3) over TS2, with one covariate affecting the drift only.
# This analysis makes sense as the TS2's drift is a function ox 'x2', 
# i.e., 'x2' affects the 'drift' parameter only. The noise variance 
# (var.arg = TRUE) is estimated, as intercept-only. See the 'zero' argument.


#------------------------------------------------------------------------#
# This model CANNOT be fitted using arima()
#------------------------------------------------------------------------#
fit.Ex2 <- vglm(TS2 ~ x2,  ARXff(order = 3, zero = c("noiseVar", "ARcoeff"), 
                             var.arg = TRUE), 
                  ## constraints = cm.ARMA(Model = ~ 1, lags.cm = 3, Resp = 1),
              data = tsdata,  trace = TRUE, crit = "log")

# True are theta1 <- 0.45; theta2 <- 0.31; theta3 <- 0.10
coef(fit.Ex2, matrix = TRUE)
summary(fit.Ex2)    
vcov(fit.Ex2)
BIC(fit.Ex2)
constraints(fit.Ex2)


# EXAMPLE 3. Fitting an ARX(3) on two responses TS1, TS2; intercept-only model with
#  constraints over the drifts. Here, 
# a) No checks on invertibility performed given the use of cm.ARMA().
# b) Only the drifts are modeled in terms of 'x2'. Then,  'zero' is
# set correspondingly.
#------------------------------------------------------------------------#
# arima() does not handle this model.
#------------------------------------------------------------------------#
fit.Ex3 <- vglm(cbind(TS1, TS2) ~ x2, ARXff(order = c(3, 3), 
                            zero = c("noiseVar", "ARcoeff"), var.arg = TRUE), 
             constraints = cm.ARMA(Model = ~ 1 + x2, lags.cm = c(3, 3), Resp = 2),
             trace = TRUE, data = tsdata, crit = "log")

coef(fit.Ex3, matrix = TRUE)
summary(fit.Ex3)
vcov(fit.Ex3)
constraints(fit.Ex3)

Run the code above in your browser using DataLab

Description

Usage

Value

Arguments

Author

Warning

Details

References

See Also

Examples