# VGLTSMs family functions: Order--p Autoregressive Model with covariates 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. ```r 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) ``` ## 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`. ## 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,%Y[t] | \Phi[t - 1] = \mu[t] +\theta[1] Y[t - 1] + \ldots +\theta[p] Y[t - p] + e[t], $$ where $x[t]$ a (possibly time--varying) covariate vector and c("$\n$", "$\\mu[t] = \\mu[*] + \\beta^T * x[t]\n$") is a (time--dependent) scaled--mean, known as **drift**. At this stage, conditional Gaussian white noise, $e[t] | \Phi [t - 1]$ is handled, in the form $$ \varepsilon_t | \Phi_{t - 1} \simN(0, \sigma^2_{\varepsilon_t | \Phi_{t - 1}}).e[t] | \Phi[t - 1] ~ N(0 , \sigma[e]^2 | \Phi [t - 1] ). $$ The distributional assumptions on the observations are then $$ Y_t | \Phi_{t - 1} \simN(\mu_{t | \Phi_{t - 1}},\sigma^2_{\varepsilon_t | \Phi_{t - 1}}),Y[t] | \Phi[t - 1] ~ N(\mu[t] | \Phi [t - 1] ,\sigma[e]^2 | \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}.\mu[t] | \Phi [t - 1] = \mu[t] +\beta^T x_t + \theta[1] Y[t - 1] +\ldots + \theta[p] Y[t - p]. $$ $\Phi[t]$ denotes the information of the joint process c("$\n$", "$(Y[t], x[t + 1])\n$"), 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.\eta = (\mu[t], log \sigma[e]^2 | \Phi[t - 1],\theta[1], \ldots, \theta[p] )^T. $$ Note, the covariates may also intervene in the conditional variance model c("$\n$", "$log \\sigma[e]^2 | \\Phi[t - 1].\n$") Hence, this family function does not restrict the noise to be **strictly** white noise (in the sense of **constant variance**). The unconditional mean, c("$%$", "$\n$", "$E(Y[t]) = \\mu $"), satisfies $$ \mu \rightarrow \frac{\mu^{\star}}{1 - (\theta_1 +\ldots + \theta_p)}%\mu -> \mu[*] / ( 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. ## Returns An object of class `"vglmff"` (see `vglmff-class`). The object is used by VGLM/VGAM modelling functions, such as `vglm` or `vgam`. ## 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. ## Author(s) Victor Miranda and T. W. Yee ## Note `zero` can be either an **integer** vector or a vector of character strings specifying either the position(s) or name(s) (partially or not) of the parameter(s) modeled as intercept-only. Numeric values can be set as usual (See `CommonVGAMffArguments`). Character strings can be entered as per parameter names in this family function, given by: `c("drift", "noiseVar" or "noiseSD", "ARcoeff")`. Users can modify the `zero` argument according to their needs. By default, $\mu[t]$ and the coefficients $\theta[1], \ldots, \theta[p]$ are intercept--only. That is, c("$\n$", "$log \\sigma[e]^2 | \\Phi[t - 1]$") is modelled in terms of any explanatories entered in the `formula`. Users, however, can modify this according to their needs via `zero`. For example, set the covariates in the `drift` model, $\mu[t]$. In addition, specific constraints for parameters are handled through the function `cm.ARMA`. If `var.arg = TRUE`, this family function estimates c("$ \\sigma[e]^2 | \n$", "$\\Phi[t - 1]$"). Else, the $ \sigma[e] |\Phi[ t- 1]$ estimate is returned. For this family function the `order` is recycled. That is, `order` will be replicated up to the number of responses given in the `vglm` call is matched. ## 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. ## See Also `ARIMAXff`, `ARMAXff`, `MAXff`, `checkTS.VGAMextra`, `CommonVGAMffArguments`, `Links`, `vglm`, ## Examples ```r 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) ```

ARXff function