Package 'VGAMextra'

Additions and extensions of the VGAM package.

Description

VGAMextra supplies additional functions and methods to the package VGAM, under three main topics:

** Time series modelling. A novel class of VGLMs to model univariate time series, called vector generalized linear time series models (VGLTSMs). It is characterized by incorporating past information in the VGLM/VGAM loglikelihood.

** 1–parameter distribution mean modelling. We return full circle by developing new link functions for the mean of 1–parameter distributions. VGAMs, VGLMs and GAMLSSs are restricted to location, scale and shape, however, the VGLM/VGAM framework has infrastructure to accommodate new links as a function of the parameters.

** Quantile modelling of 1–parameter distributions. Similarly, we have implemented link functions to model the quantiles of several 1–parameter distributions, for quantile regression.

Details

The inference infrastructure of VGAMextra relies on the VGLM/VGAM framework. Particularly, estimation is carried out via IRLS using Fisher scoring, whilst additive models and reduced rank regression are also accommodated by all VGAMextra family functions.

At present, this package allows the extent of VGLMs/VGAMs to operate popular time series models as special cases, as well as cointegrated time series (bivariate case), and modelling choices for volatility models incorporating explanatories in the variance equation. The central family functions in this respect are ARXff, MAXff, ARMAX.GARCHff, and VGLM.INGARCHff.

Regarding modelling the mean/quantile-functions, VGAMextra affords links for several 1-parameter distributions, e.g., expMlink, benini1Qlink, or inv.chisqMlink. Collectively, the quantile-links represent an alternative to quantile regression by directly modelling the quantile function for distributions beyond the exponential family (See Example 3 below).

The VGLM/VGAM framework is very large and encompasses a wide range of multivariate response types and models, e.g., it includes univariate and multivariate distributions, categorical data analysis, and extreme values. See VGAM-package for a broad description and further information on VGAM.

Future work

* Implement VGLM time series family functions to handle error correction models (ECMs) for cointegrated time series. Upgrade this framework beyond the bivariate case, e.g., the the Vector ECM (VECMs).

* Upgrade VGLMs time series family functions to handle multivariate time series, e.g., the VAR model (Coming shortly).

* Incorporate VGLM/VGAM-links to model the mean and quantile functions of distributions with > 1 parameters.

* Develop the class of multiple reduced–rank VGLMs towards time series, to handle, e.g., vector error correction models (VECMs), for multiple cointegrated time series.

Warning

VGAMextra is revised, altered, and/or upgraded on a regular basis. Hence, be aware that any feature, e.g., function names, arguments, or methods, may be modified without prior notice. Check the NEWS for the latest changes and additions across the different versions.

Author(s)

Victor Miranda, [email protected].

Maintainer: Victor Miranda, [email protected].

Contributor: Thomas Yee, [email protected].

References

Pfaff, B. (2011) Analysis of Integrated and Cointegrated Time Series with R. Seattle, Washington, USA: Springer.

Chan, N., Li, D., Peng, L. and Zhang, R. (2013) Tail index of an AR(1) model with ARCH(1) errors. Econometric Theory, 29(5):920-940.

Yee, T. W. (2015) Vector Generalized Linear and Additive Models: With an Implementation in R. New York, USA: Springer.

Miranda, V. and Yee, T. W. Vector generalized linear time series models. In preparation.

Miranda, V. and Yee, T. W. On mean modelling of 1-parameter distributions using vector generalized linear models. In preparation

Miranda, V. and Yee, T. W. Two–Parameter Link Functions with Applications to Negative Binomail, Weibull, and Quantile Regression. To be submitted to the Scandinavian Journal of Statistics.

Miranda, V. and Yee, T.W. New Link Functions for Distribution-Specific Quantile Regression Based on Vector Generalized Linear and Additive Models. Journal of Probability and Statistics, Volume 2019, Article ID 3493628.

Yee, T. W. (2008) The VGAM Package. R News, 8, 28–39.

See Also

vglm, vgam, rrvglm, Links, CommonVGAMffArguments, https://www.stat.auckland.ac.nz/~vmir178/.

Examples

##### EXAMPLE 1. An AR(1) model with ARCH(1) errors.
# Chan et.al. (2013) proposed a long and technical methodology to
# estimate the tail index of an AR(1) with ARCH(1) errors involving 
# its estimation by QMLE. I fit this model straightforwardly by MLE
# using the family function ARXff() for time series, and constraining
# the effect of Y^2_{t - 1} to the conditional variance using
# constraint matrices.


# Generate some data
set.seed(1)
nn       <- ceiling(runif(1, 150, 160))
my.rho   <- rhobitlink(-1.0, inverse = TRUE)  # -0.46212
my.mu    <- 0.0
my.omega <- 1
my.b     <- 0.5
tsdata   <- data.frame(x2 =  sort(runif(n = nn)))
tsdata   <- transform(tsdata, index = 1:nn, TS1   = runif(nn))

for (ii in 2:nn)
  tsdata$TS1[ii] <- my.mu + my.rho * tsdata$TS1[ii-1] +
  sqrt(my.omega + my.b * (tsdata$TS1[ii-1])^2) * rnorm(1)
  
# Remove the burn-in data:
nnr <- ceiling(nn/5)
tsdata <- tsdata[-(1:nnr), ]
tsdata["index"] <- 1:(nn - nnr)


# The timeplot.
with(tsdata, plot(ts(TS1), lty = "solid", col = "blue", xlab ="Time", ylab = "Series"))
abline(h = mean(tsdata[, "TS1"]), lty = "dotted")

# The constraint matrices, to inhibit the effect of Y^2_{t - 1}
# over sigma^2 only.
const.mat <- list('(Intercept)' = diag(3), 'TS1l1sq' = cbind(c(0, 1, 0)))

# Set up the data using function WN.lags() from VGAMextra to generate
# our 'explanatory'
tsdata <- transform(tsdata, TS1l1sq = WN.lags(y = cbind(tsdata[, "TS1"])^2, lags = 1))

# Fitting the model
fit.Chan.etal <- vglm(TS1 ~ TS1l1sq, ARXff(order = 1,   # AR order 
                                          zero = NULL, noChecks = FALSE,
                                          var.arg = TRUE, lvar = "identitylink"),
             crit = "loglikelihood", trace = TRUE,
             constraints = const.mat, data = tsdata)  ## Constraints...
summary(fit.Chan.etal, lrt0 = TRUE, score0 = TRUE, wald0 = TRUE)
constraints(fit.Chan.etal)


###### EXAMPLE 2. VGLMs handling cointegrated (bivariate) time series.
# In this example, vglm() accommodates an error correction model
# of order (2, 2) to fit two (non-stationary) cointegrated time series.

# Simulating some data.
set.seed(2017081901)
nn <- 280
rho <- 0.75
s2u <- exp(log(1.5))  # Gaussian noise1
s2w <- exp(0)         # Gaussian noise2
my.errors <- rbinorm(nn, mean1 = 0, mean2 = 0, var1 = s2u, var2 = s2w, cov12 = rho)
ut <- my.errors[, 1]
wt <- my.errors[, 2]
yt <- xt <- numeric(0)

xt[1] <- ut[1]     # Initial value: error.u[0]
yt[1] <- wt[1]     # Initial value: error.w[0]
beta <- c(0.0, 2.5, -0.32)  # Coefficients true values.

for (ii in 2:nn) {
  xt[ii] <-  xt[ii - 1] + ut[ii]
  yt[ii] <- beta[1] + beta[2] * xt[ii]  + beta[3] * yt[ii - 1] + wt[ii]
}

# Regression of yt on xt, save residuals. Compute Order--1 differences.
errors.coint <- residuals(lm(yt ~ xt)) # Residuals from the static regression yt ~ xt
difx1 <- diff(ts(xt), lag = 1, differences = 1)  # First difference for xt
dify1 <- diff(ts(yt), lag = 1, differences = 1)  # First difference for yt

# Set up the dataset (coint.data), including Order-2 lagged differences.
coint.data <- data.frame(embed(difx1, 3), embed(dify1, 3)) 
colnames(coint.data) <- c("difx1", "difxLag1", "difxLag2",
                          "dify1", "difyLag1", "difyLag2")

# Remove unutilized lagged errors accordingly. Here, use from t = 3.
errors.cointLag1 <- errors.coint[-c(1:2, nn)]
coint.data <- transform(coint.data, errors.cointLag1 = errors.cointLag1)


# Plotting the data
plot(ts(yt[-c(1:3, NULL)]), lty = 4, type = "l", col = "red",
     main = "", xlab = "Time", las = 1, ylim = c(-32, 20),
     ylab = expression("Series"~x[t]~"and"~y[t]))
lines(ts(xt[-c(1:3, NULL)]), lty = 4, type = "l", col = "blue")
legend("bottomleft", c(expression("Series"~x[t]),
                       expression("Series"~y[t])),
       col = c("red", "blue"), lty = c(4, 4))


# Fitting an error correction model (2, 2), aka ECM(2, 2)
fit.coint <- vglm(cbind(dify1, difx1) ~ errors.cointLag1 + difxLag1 + difyLag1 + 
                                           difxLag2 + difyLag2,
          binormal(zero = c("sd", "rho")),  # 'sigma', 'rho' are intercept--only.
          trace = FALSE, data = coint.data)
summary(fit.coint)
coef(fit.coint, matrix = TRUE)


##### EXAMPLE 3. Quantile Modelling (QM).
# Here, the quantile function of the Maxwell distribution is modelled
# for percentiles 25%, 50% and 75%. The resulting quantile-curves
# are plotted. The rate parameter is determined by an artificial covariate.

set.seed(123)
# An artificial covariate.
maxdata <- data.frame(x2 = sort(runif(n <- nn <- 120)))
# The 'rate' function.
mymu <- function(x) exp(2 -  6 * sin(2 * x - 0.2) / (x + 0.5)^2)
# Set up the data.
maxdata <- transform(maxdata, y = rmaxwell(n, rate = mymu(x2)))

# 25%, 50% and 75% quantiles are to be modelled.
mytau <- c(0.25, 0.50, 0.75)
mydof <- 4

### Using B-splines with 'mydof'-degrees of freedom on the predictors
fit.QM <- vglm(Q.reg(y, pvector = mytau) ~ bs(x2, df = mydof),
                 family = maxwell(link = maxwellQlink(p = mytau), zero = NULL),
                 data = maxdata, trace = TRUE)
            
summary(fit.QM, lscore0 = TRUE)
head(predictors(fit.QM))         # The 'fitted values'


## The 25%, 50%, and 75% quantile curves.
mylwd <- 1.5
with(maxdata, plot(x2, y, col = "orange",
                 main = "Example 1; Quantile Modelling",
                 ylab = "y", pch = "o", cex = 0.75))
with(maxdata, matlines(x2, predict(fit.QM)[, 1], col = "blue",
                       lty = "dotted", lwd = mylwd))
with(maxdata, matlines(x2, predict(fit.QM)[, 2], col = "chocolate",
                       lty = "dotted", lwd = mylwd))
with(maxdata, matlines(x2, predict(fit.QM)[, 3], col = "brown",
                       lty = "dotted", lwd = mylwd))
legend("topleft", c("percentile25", "percentile50", "percentile75"),
       lty = rep("dotted", 3), lwd = rep(mylwd, 3))


### Double check: The data (in percentage) below the 25%, 50%, and 75% curves
round(length(predict(fit.QM)[, 1][(maxdata$y
                 <= predict(fit.QM)[, 1] )]) /nn, 5) * 100  ## Should be 25% approx
round(length(predict(fit.QM)[, 2][(maxdata$y
                 <= predict(fit.QM)[, 2] )]) /nn, 5) * 100  ## Should be 50% approx
round(length(predict(fit.QM)[, 3][(maxdata$y
                 <= predict(fit.QM)[, 3] )]) /nn, 5) * 100  ## Should be 75% approx

---

Air pollution Data, Mexico City.

Description

Daily air pollution levels in Mexico City, January 2004 – June 2005.

Usage

data(ap.mx)

Format

This data frame stores time series vectors with the following information:

time

Time vector.

PM10

24–hr average concentration of $PM_{10}$, in micrograms per milliliter.

O3

Daily maximum 8–hour moving average of ozone, in micrograms per milliliter.

temp

Daily mean average of temperature, in celsius degrees.

HR

Daily mean average (%) of relative humidity.

Details

These are readings of $PM_{10}$, $O_{3}$, temperature and humidity between 1 January 2004 and 30 June 2005 in Mexico City Metropolitan Area. Each observation is the 24–hr mean average (between 00:00 and 23:59 hrs), except for ozone, where the maximum over all the sliding 8–hour–windows, between 00:00 and 23:59 hrs is reported, viz. the daily maximum 8–hour moving average.

Source

National Institute of Ecology. Gathers and disseminates the data generated by the central air quality monitoring network in Mexico City. Website: https://www.gob.mx/inecc/

Examples

data(ap.mx)
summary(ap.mx[, -1])
class(ap.mx[, "PM10"])

layout(matrix(c(1, 1, 2,3), 2, 2, byrow = TRUE))
plot.ts(ts(ap.mx$PM10), ylab = expression(PM[10]~"Series"), 
        col = "brown", xaxt = "n", las = 1)
xtick <- c(1, 92, 183, 275, 367, 457, 518)
xtext <- c("Jan/04", "April/04", "July/04", "Oct/04", "Jan/05",
           "April/05", "June/05")
axis(side = 1, at = xtick, labels = FALSE)
text(x = xtick, par("usr")[3], labels = xtext,
     pos = 1, xpd = TRUE, col = "black")
pacf(ap.mx$PM10, main = "", ylim= c(-0.5, 1), lag.max = 60, las = 1)
acf(ap.mx$PM10, main = "", ylim= c(-0.5, 1), lag.max = 60, las = 1)

---

VGLTSMs Family Functions: Generalized integrated regression with order–$(p, q)$ ARMA errors

Description

A VLTSMff for dynamic regression. Estimates regression models with order–$(p, d, q)$ ARIMA errors by maximum likelihood.

Usage

ARIMAX.errors.ff(order = c(1, 1, 1),
                       zero = "var",  # optionally, "mean".
                       order.trend = 0,
                       include.int = TRUE,
                       diffCovs  = TRUE,
                       xLag = 0,
                       include.currentX = TRUE,
                       lvar = "loglink",
                       lmean = "identitylink")

Arguments

order	The usual $(p, d, q)$ integer vector as in, e.g., ARIMAXff. By default, an order–$(p, q)$ ARMA model is fitted on the errors, whlist $d$ is the degree of differencing on the response.
zero	What linear predictor is modelled as intercept–only? See zero and CommonVGAMffArguments for further details.
order.trend	Non–negative integer. Allows to incorporate a polynomial trend of order order.trend in the forecast mean function.
include.int	Logical. Should an intercept (int) be included in the model for $y_t$? Default is TRUE. See below for details.
diffCovs	Logical. If TRUE (default) the order–d difference of the covariates is internally computed and then incorporated in the regression model. Else, only the current values are included.
xLag	Integer. If entered, the covariates, say $\boldsymbol{x}_t$ are laggeg (up to order xLag) and then embedded in the regression model. See below for further details.
include.currentX	Logical. If TRUE, the actual values, $\boldsymbol{x}_t$, are included in the regression model. Else, this is ignored and only the lagged $\boldsymbol{x}_{t - 1}, \ldots, \boldsymbol{x}_{t - xLag}$ will be included.
lvar, lmean	Link functions applied to conditional mean and the variance. Same as uninormal.

Details

The generalized linear regression model with ARIMA errors is another subclass of VGLTSMs (Miranda and Yee, 2018).

For a univariate time series, say $y_t$, and a $p$–dimensional vector of covariates $\boldsymbol{x}_t$ covariates, the model described by this VGLTSM family function is

$y_t = \boldsymbol{\beta}^T \boldsymbol{x}_t + u_t,$

$u_t = \theta_1 u_{t - 1} + \cdots + \theta_p u_{t - p} + z_t + \phi_1 z_{t - 1} + \cdots + \phi_1 z_{t - q}.$

The first entry in $x_t$ equals 1, allowing an intercept, for every $t$. Set include.int = FALSE to set this to zero, dimissing the intercept.

Also, if diffCovs = TRUE, then the differences up to order d of the set $\boldsymbol{x}_t$ are embedded in the model for $y_t$. If xLag$> 0$, the lagged values up to order xLag of the covariates are also included.

The random disturbances $z_t$ are by default handled as $N(0, \sigma^2_z)$. Then, denoting $\Phi_{t}$ as the history of the process $(x_{t + 1}, u_t)$ up to time $t$, yields

$E(y_t | \Phi_{t - 1}) = \boldsymbol{\beta}^T \boldsymbol{x}_t + \theta_1 u_{t - 1} + \cdots + \theta_p u_{t - p} + \phi_1 z_{t - 1} + \cdots + \phi_1 z_{t - q}.$

Denoting $\mu_t = E(y_t | \Phi_{t - 1}),$ the default linear predictor for this VGLTSM family function is

$\boldsymbol{\eta} = ( \mu_t, \log \sigma^2_{z})^T.$

Value

An object of class "vglmff" (see vglmff-class) to be used by VGLM/VGAM modelling functions, e.g., vglm or vgam.

Note

If d = 0 in order, then ARIMAX.errors.ff will perform as ARIMAXff.

Author(s)

Victor Miranda

See Also

ARIMAXff, CommonVGAMffArguments, uninormal, vglm.

Examples

### Estimate a regression model with ARMA(1, 1) errors.
## Covariates are included up to lag 1.
set.seed(20171123)
nn <- 250
x2 <- rnorm(nn)                                    # One covariate
sigma2 <- exp(1.15); theta1 <- 0.5; phi1 <- 0.27   # True coefficients
beta0  <- 1.25; beta1 <- 0.25; beta2 <- 0.5

y <- numeric(nn)
u <- numeric(nn)
z <- numeric(nn)

u[1] <- rnorm(1)
z[1] <- rnorm(1, 0, sqrt(sigma2))

for(ii in 2:nn) {
  z[ii] <- rnorm(1, 0, sqrt(sigma2))
  u[ii] <- theta1 * u[ii - 1] + phi1 * z[ii - 1] + z[ii]
  y[ii] <- beta0 + beta1 * x2[ii] + beta2 * x2[ii - 1] + u[ii]
}

# Remove warm-up values.
x2 <- x2[-c(1:100)]
y  <- y[-c(1:100)]


plot(ts(y), lty = 2, col = "blue", type = "b")
abline(h = 0, lty = 2)


## Fit the model.
ARIMAX.reg.fit <- vglm(y ~ x2, ARIMAX.errors.ff(order = c(1, 0, 1), xLag = 1),
             data = data.frame(y = y, x2 = x2), trace = TRUE)
coef(ARIMAX.reg.fit, matrix = TRUE)
summary(ARIMAX.reg.fit, HD = FALSE)


# Compare to arima()
# arima() can't handle lagged values of 'x2' by default, but these 
# may entered at argument 'xreg'.
arima(y, order = c(1, 0, 1), xreg = cbind(x2, c(0, x2[-150])))

---

VGLTSMs Family functions: The Order–$(p, d, q)$ Autoregressive Integrated Moving Average Model (ARIMA(p, d, q)) with covariates

Description

Maximum likelihood estimation of the drift, standard deviation or variance of the random noise, and coefficients of an autoregressive integrated moving average process of order-$(p, d, q)$ with covariates by MLE using Fisher scoring. No seasonal terms handled yet. No seasonal components handled yet.

Usage

ARIMAXff(order   = c(1, 1, 0),
               zero     = c("ARcoeff", "MAcoeff"),
               diffCovs = TRUE,
               xLag     = 0,
               include.current = FALSE,
               type.EIM = c("exact", "approximate")[1], 
               var.arg  = TRUE,
               nodrift  = FALSE,
               noChecks = FALSE,
               ldrift   = "identitylink", 
               lsd      = "loglink",
               lvar     = "loglink",
               lARcoeff = "identitylink",
               lMAcoeff = "identitylink", 
               idrift   = NULL,
               isd      = NULL,
               ivar     = NULL,
               iARcoeff = NULL, 
               iMAcoeff = NULL)

Arguments

order	Integer vector with three components, (p, d, q): The AR order (p), the degree of differencing (d), and the MA order (q).
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.
diffCovs	Logical. The default is diffCovs = TRUE, which means that the order–d differences of the covariates (if entered) are internally computed and then incorporated in the conditional–mean model. Otherwise, only the current actual values of the covariates are included.
xLag	Integer, non–negative. If xLag > 0, the covariates at current time, $\boldsymbol{x}_t$, plus the lagged covariates up to order 'xLag' are embedded into the design matrix (as covariates too). Leave xLag = 0 and only the current value, $\boldsymbol{x}_t$, will be considered. See more details below.
include.current	Logical. Same as ARIMAX.errors.ff.
type.EIM	The type of expected information matrix (EIM) of the ARMA process to be utilized in Fisher scoring. type.EIM = "exact" (default) enables the exact IM (Porat, et.al. 1986), otherwise the approximate version is utilized.
var.arg	Logical. If FALSE (default), then the standard deviation of the random noise is estimated. Else, the variance estimate is returned.
nodrift	Logical. nodrift = TRUE supresses estimation of the intercept (the drift in the ARMA case), which is set to zero internally.
noChecks	Logical. If FALSE (default), this family function internally checks stationarity (AR case) and invertibility (MA case) of the the estimated model. A warning is correspondingly displayed.
ldrift, lsd, lvar, lARcoeff, lMAcoeff	Link functions applied to the intercept, the random noise standard deviation (or optionally, the variance), and the coefficients in the ARMA–type conditional–mean model.
idrift, isd, ivar, iARcoeff, iMAcoeff	Optional initial values for the intercept (drift), noise SD (or variance), and ARMA coeffcients (a vector of length $p + q$). If failure to converge occurs then try different values and monitor convergence by using trace = TRUE in the vglm() call.

Details

Let $\boldsymbol{x}_t$ be a (probably time–varying) vector of suitable covariates. The ARIMAX model handled by ARIMAXff is

$\nabla^d Y_t = \mu^{\star} + \boldsymbol{\beta}^T \nabla^d \boldsymbol{x}_t + \theta_1 \nabla^d Y_{t - 1} + \cdots + \theta_p \nabla^d Y_{t - p} + \phi_1 \varepsilon_{t - 1} + \cdots + \phi_q \varepsilon_{t - q} + \varepsilon,$

with $\nabla^d (\cdot)$ the operator differencing a series d times. If diffCovs = TRUE, this function differencing the covariates d times too.

Similarly, ARMAXff manages

$\nabla^d Y_t = \mu^{\star} + \boldsymbol{\beta}^T \boldsymbol{x}_t + \theta_1 Y_{t - 1} + \cdots + \theta_p Y_{t - p} + \phi_1 \varepsilon_{t - 1} + \cdots + \phi_q \varepsilon_{t - q} + \varepsilon,$

where

$\varepsilon_{t | \Phi_{t - 1}} \sim N(0, \sigma_{\varepsilon_t | \Phi_{t - 1}}^2).$

Note, $\sigma_{\varepsilon | \Phi_{t - 1}}^2$ is conditional on $\Phi_{t - 1}$, the information of the joint process $\left(Y_{t - 1}, \boldsymbol{x}_t \right)$, at time $t$, and hence may be modelled in terms of $\boldsymbol{x}_t$, if required.

ARIMAXff() and ARMAXff() handle multiple responses, thus a matrix can be used as the response. Note, no seasonal terms handled. This feature is to be incorporated shortly.

The default linear predictor is

$\boldsymbol{\eta} = \left( \mu, \log \sigma^{2}_{\varepsilon_{t | \Phi_{t - 1}}}, \theta_1, \ldots, \theta_p, \phi_1, \ldots, \phi_q \right)^T.$

Other links are also handled. See Links.

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

As with ARXff and MAXff, choices for the EIMs are "exact" and "approximate". Covariates may be incorporated in the fit for any linear predictor above. Hence, ARIMAXff supports non–stationary processes ($\sigma_{\varepsilon_t | \Phi_{t - 1}}^2$) may depend on $\boldsymbol{X}_t$. Also, constraint matrices on the linear predictors may be entered through cm.ARMA or using the argument constraints, from vglm.

Checks on stationarity and invertibility on the esitmated process are performed by default. Set noChecks = TRUE to dismiss this step.

Value

An object of class "vglmff" (see vglmff-class) to be used by VGLM/VGAM modelling functions, e.g., vglm or vgam.

Warning

zero can be a numeric or a character–strings vector specifying the position(s) or the name(s) of the parameter(s) modeled as intercept–only. Numeric values can be set as usual (See CommonVGAMffArguments). If names are entered, the parameter names in this family function are:

c("drift.mean", "noiseVar" || "noiseSD", "ARcoeff", "MAcoeff").

Manually modify this if required. For simplicity, the second choice is recommended.

Note

No seasonal components handled yet.

If no covariates, $\boldsymbol{x}_t$, are incorporated in the analysis, then ARIMAXff fits an ordinary ARIMA model. Ditto with ARMAXff.

If nodrift = TRUE, then the 'drift' is removed from the vector of parameters and is not estimated.

By default, an ARMA model of order–$c(1, 0)$ with order–1 differences is fitted. When initial values are entered (isd, iARcoeff, etc.), they are recycled according to the number of responses.

Also, the ARMA coefficients are intercept–only (note, zero = c("ARcoeff", "MAcoeff")) This may altered via zero, or by constraint matrices (See constraints) using cm.ARMA.

Checks on stationarity and/or invertibility can be manually via checkTS.VGAMextra.

Author(s)

Victor Miranda and T. W. Yee

References

Miranda, V. and Yee, T.W. (2018) Vector Generalized Linear Time Series Models. In preparation.

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.

See Also

ARXff, MAXff, checkTS.VGAMextra, cm.ARMA, CommonVGAMffArguments, constraints, vglm.

Examples

set.seed(3)
nn      <- 90
theta   <- c(0.12, 0.17)  # AR coefficients
phi     <- c(-0.15, 0.20)  # MA coefficients.
sdWNN <- exp(1.0)          # SDs
mu    <- c(1.25, 0.85)     # Mean (not drift) of the process.
covX  <- runif(nn + 1)     # A single covariate.
mux3  <- mu[1] + covX
##
## Simulate ARMA processes. Here, the drift for 'tsd3' depends on covX.
##
tsdata <- data.frame(TS1 = mu[1] + arima.sim(model = list(ar = theta, ma = phi,
                             order = c(2, 1, 2)), n = nn, sd = sdWNN ),
                      TS2 = mu[2] + arima.sim(model = list(ar = theta, ma = phi,
                            order = c(2, 1, 2)), n = nn, sd = exp(2 + covX)),
                      TS3 =  mux3 + arima.sim(model = list(ar = theta, ma = phi,
                            order = c(2, 1, 2)), n = nn, sd = exp(2 + covX) ),
                      x2 = covX)

### EXAMPLE 1. Fitting a simple ARIMA(2, 1, 2) using vglm(). 
# Note that no covariates involved. 
fit.ARIMA1 <-  vglm(TS1 ~ 1, ARIMAXff(order = c(2, 1, 2), var.arg = FALSE,
                       # OPTIONAL INITIAL VALUES
                       # idrift = c(1.5)*(1 - sum(theta)), 
                       # ivar = exp(4), isd = exp(2),
                       # iARcoeff = c(0.20, -0.3, 0.1),
                       # iMAcoeff = c(0.25, 0.35, 0.1),
                        type.EIM = "exact"), 
                 data = tsdata, trace = TRUE, crit = "log")
Coef(fit.ARIMA1)
summary(fit.ARIMA1)
vcov(fit.ARIMA1, untransform = TRUE)
#------------------------------------------------------------------------#
# Fitting same model using arima().
#------------------------------------------------------------------------#
# COMPARE to EXAMPLE1 
( fitArima  <- arima(tsdata$TS1, order = c(2, 1, 2)) ) 


### EXAMPLE 2. Here only the ARMA coefficients and drift are intercept-only.
# The random noise variance is not constant.
fit.ARIMA2 <-  vglm(TS2 ~ x2,  ARIMAXff(order = c(2, 1, 2), var.arg = TRUE,
                     lARcoeff = "rhobitlink", lMAcoeff = "identitylink",
                     type.EIM = c("exact", "approximate")[1], 
                     # NOTE THE ZERO ARGUMENT. 
                     zero = c("drift.mean", "ARcoeff", "MAcoeff")),
              data = tsdata, trace = TRUE)

coef(fit.ARIMA2, matrix = TRUE)
summary(fit.ARIMA2)
constraints(fit.ARIMA2)



### EXAMPLE 3. Here only ARMA coefficients are intercept-only.
# The random noise variance is not constant.
# Note that the "drift" and the "variance" are "generated" in 
# terms of 'x2' above for TS3.

fit.ARIMA3 <- vglm(TS3 ~ x2,  ARIMAXff(order = c(1, 1, 2), var.arg = TRUE,
                     lARcoeff = "identitylink", lMAcoeff = "identitylink",
                     type.EIM = c("exact", "approximate")[1], nodrift = FALSE,
                     zero = c( "ARcoeff", "MAcoeff")), # NOTE THE ZERO ARGUMENT. 
              data = tsdata, trace = TRUE)
              
coef(fit.ARIMA3, matrix = TRUE)
summary(fit.ARIMA3)
constraints(fit.ARIMA3)

---

VGLTSMs Family Functions: Generalized autoregressive moving average model with Student-t errors

Description

For an ARMA model, estimates a 3–parameter Student-$t$ distribution characterizing the errors plus the ARMA coefficients by MLE usign Fisher scoring. Central Student–t handled currently.

Usage

ARMA.studentt.ff(order = c(1, 0),
                             zero = c("scale", "df"),
                             cov.Reg = FALSE,
                             llocation = "identitylink",
                             lscale    = "loglink",
                             ldf       = "logloglink",
                             ilocation = NULL,
                             iscale = NULL,
                             idf = NULL)

Arguments

order	Two–entries vector, non–negative. The order $u$ and $v$ of the ARMA model.
zero	Same as studentt3.
cov.Reg	Logical. If covariates are entered, Should these be included in the ARMA model as a Regressand? Default is FALSE, then only embedded in the linear predictors.
llocation, lscale, ldf, ilocation, iscale, idf	Same as studentt3.

Details

The normality assumption for time series analysis is relaxed to handle heavy–tailed data, giving place to the ARMA model with shift-scaled Student-$t$ errors, another subclass of VGLTSMs.

For a univariate time series, say $y_t$, the model described by this VGLTSM family function is

$\theta(B)y_t = \phi(B) \varepsilon_t,$

where $\varepsilon_t$ are distributed as a shift-scaled Student–$t$ with $\nu$ degrees of freedom, i.e., $\varepsilon_t \sim t(\nu_t, \mu_t, \sigma_t)$. This family functions estimates the location ($\mu_t$), scale ($\sigma_t$) and degrees of freedom ($\nu_t$) parameters, plus the ARMA coefficients by MLE.

Currently only centered Student–t distributions are handled. Hence, the non–centrality parameter is set to zero.

The linear/additive predictors are $\boldsymbol{\eta} = (\mu, \log \sigma, \log \log \nu)^T,$ where $\log \sigma$ and $\nu$ are intercept–only by default.

Value

An object of class "vglmff" (see vglmff-class) to be used by VGLM/VGAM modelling functions, e.g., vglm or vgam.

Note

If order = 0, then AR.studentt.ff fits a usual 3–parameter Student–$t$, as with studentt3.

If covariates are incorporated in the analysis, these are embedded in the location–parameter model. Modify this through zero. See CommonVGAMffArguments for details on zero.

Author(s)

Victor Miranda

See Also

ARIMAXff, studentt, vglm.

Examples

### Estimate the parameters of the errors distribution for an
## AR(1) model. Sample size = 50

set.seed(20180218)
nn <- 250
y  <- numeric(nn)
ncp   <- 0           # Non--centrality parameter
nu    <- 3.5         # Degrees of freedom.
theta <- 0.45        # AR coefficient
res <- numeric(250)  # Vector of residuals.

y[1] <- rt(1, df = nu, ncp = ncp)
for (ii in 2:nn) {
  res[ii] <- rt(1, df = nu, ncp = ncp)
  y[ii] <- theta * y[ii - 1] + res[ii]
}
# Remove warm up values.
y <- y[-c(1:200)]
res <- res[-c(1:200)]

### Fitting an ARMA(1, 0) with Student-t errors.
AR.stut.er.fit <- vglm(y ~ 1, ARMA.studentt.ff(order = c(1, 0)),
                       data = data.frame(y = y), trace = TRUE)

summary(AR.stut.er.fit)
Coef(AR.stut.er.fit)


plot(ts(y), col = "red", lty = 1, ylim = c(-6, 6), main = "Plot of series Y with Student-t errors")
lines(ts(fitted.values(AR.stut.er.fit)), col = "blue", lty = 2)
abline( h = 0, lty = 2)

---

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)

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

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.

Value

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

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, $\log \sigma^{2}_{\varepsilon_{t | \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 $\sigma_{\varepsilon_t | \Phi_{t - 1}}^2$. Else, the $\sigma_{\varepsilon_t | \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.

Author(s)

Victor Miranda and T. W. Yee

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.

See Also

ARIMAXff, ARMAXff, MAXff, checkTS.VGAMextra, CommonVGAMffArguments, Links, vglm,

Examples

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)

---

Link functions for the quantiles of several 1–parameter continuous distributions

Description

Computes the benini1Qlink transformation, its inverse and the first two derivatives.

Usage

benini1Qlink(theta, p = stop("Argument 'p' must be entered."),
               y0 = stop("Argument 'y0' must be specified."),
               bvalue = NULL, inverse = FALSE,
               deriv = 0, short = TRUE, tag = FALSE)

Arguments

theta	Numeric or character. See below for further details.
p	Numeric. A single value between 0.0 and 1.0. It is the $p$–quantile to be modeled by this link function.
y0	Same as benini1.
bvalue, inverse, deriv, short, tag	See Links.

Details

This is a link function to model any $p$–quantile of the 1–parameter Benini distribution. It is called the benini1Qlink transformation defined as

$\log y_0 + \sqrt{\frac{ -\log (1 - p) }{s}}$

where $y_0 > 0$ is a scale parameter and $s$ is a positive shape parameter, as in benini1.

Numerical values of $s$ or $p$ out of range may result in Inf, -Inf, NA or NaN.

In particular, arguments inverse and deriv are disregarded if theta is character.

Value

For deriv = 0, the benini1Qlink transformation of theta, when inverse = FALSE. If inverse = TRUE, then the inverse transformation given by -log(1 - p) / (theta - log y0)^2 is returned.

For deriv = 1, this function returns the derivative $d$ eta / $d$ theta, if inverse = FALSE. Else, the reciprocal $d$ theta / $d$ eta as a function of theta.

If deriv = 2, then the second order derivatives in terms of theta are accordingly returned.

Warning

The horizontal straight line $\log y0$ is a lower asymptote for this link function as $\theta$ increases to $\infty$. Thus, when inverse = TRUE and deriv = 0 entries at theta becoming $\eta$ must be greater than $\log y0$. Else, NaN will be returned. See examples 2 and 3 below.

Note

Numerical instability may occur for values theta too close to zero or lower than $\log y0$. Use argument bvalue to replace them before computing the link.

Author(s)

V. Miranda and Thomas W. Yee.

See Also

benini1, Links.

Examples

## E1. benini1Qlink() and its inverse ##
   p <- 0.50; y0 = 1.25         ## Modeling the median
   my.s <- seq(0, 5, by = 0.1)[-1]
    max(my.s - benini1Qlink(benini1Qlink(my.s, p = p, y0 = y0), 
                            p = p, y0 = y0, inverse  =TRUE))    ## Zero

  ## E2. Plot of the benini1Qlink() transformation and its inverse     ##
  ## Note, inverse = TRUE implies that argument 'theta' becomes 'eta'. ##
  ## which must be greater than log(y0). Else, value less than log(y0) ##
  ## are replaced by NaN.                                              ##
  
   #--- THE LINK
   my.b <- seq(0, 5, by = 0.01)[-1]
   plot(benini1Qlink(theta = my.b, p = p, y0 = y0) ~ my.b,
        type = "l", col = "blue", lty = "dotted", lwd = 3,
        xlim = c(-0.1, 6), ylim = c(-0.1, 5), las = 1,
        main = c("Blue is benini1Qlink(), green is the inverse"),
        ylab = "eta = benini1Qlink", xlab = "theta")
   abline(h = 0, v = 0, lwd = 2)
    
   #--- THE INVERSE
   lines(my.b, benini1Qlink(theta = my.b, p = p, y0 = y0, inv = TRUE),
         col = "green", lwd = 2, lty = "dashed")
   #--- Tracing the identity function for double--check
   lines(my.b, my.b)     
   
    
   ## E3. WARNING! The first two values are less than log(y0)  ##
    benini1Qlink(theta = c(0.10, 0.15, 0.25, 0.35) , p = p, y0 = y0, inverse  = TRUE)

---

Link functions for the mean of 1–parameter discrete distributions: The Borel–Tanner distribution.

Description

Computes the borel.tannerMlink transformation, its inverse and the first two derivatives.

Usage

borel.tannerMlink(theta, Qsize = 1, 
                      bvalue = NULL, inverse = FALSE, 
                      deriv = 0, short = TRUE, tag = FALSE)

Arguments

theta	Numeric or character. See below for further details.
Qsize	A positive integer. It is called $Q$. Same as borel.tanner. Default it 1.
bvalue, inverse, deriv, short, tag	Details at Links

Details

As with zetaffMlink or yulesimonMlink, this link function is part of a set of link functions in VGAM developed under a common methodology: by taking the logarithm of the mean of the corresponding distribution.

In particular, this link function emerges by computing the logarithm of the mean of the Borel–Tanner distribution. It is defined as

$\tt{borel.tannerMlink}(a) = - \log (Q^{-1} - aQ^{-1}),$

where $a$, $0 < a < 1$, is a scale parameter as in borel.tanner.

The domain set of borel.tannerMlink is the open interval $(0, 1)$, except when inverse = TRUE and deriv = 0. See below for further details about this. Moreover, unlike zetaffMlink or posPoiMlink, the inverse of borel.tannerMlink can be written in closed–form.

Values of $a$ (i.e. theta) out of range will result in NaN of NA.

If theta is a character, arguments inverse and deriv are discarded.

Value

For deriv = 0, the borel.tannerMlink transformation of theta, if inverse = FALSE. When inverse = TRUE, theta becomes $\eta$ and the inverse of borel.tannerMlink, given by

$1 - \frac{Q}{ e^{\eta} },$

is returned. Here, the domain set changes to $(0, \infty)$.

For deriv = 1, d eta / d theta as a function of theta if inverse = FALSE, else the reciprocal d theta / d eta.

Similarly, when deriv = 2 the second order derivatives in terms of theta are returned.

References

Haight, F. and Brueuer, M. A. (1960) The Borel–Tanner distribution. Biometrika, 47, 143–150.

Note

The vertical line $a = 1$ is an asymptote for this link function, which sharply grows for values of $a$ too close to 1.0 from the left. For such cases, Inf might result when computing borel.tannerMlink.

This link function is useful to model any parameter in $(0, 1)$. Then, some problems may occur if there are covariates causing out of range values.

Author(s)

V. Miranda and T. W. Yee

See Also

borel.tanner, yulesimonMlink, zetaffMlink, posPoiMlink, Links.

Examples

## Example 1. Special values for theta (or eta, accordingly)  ##
a.par <- c(0, 1:10/10, 20, 1e1, Inf, -Inf, NaN, NA) 

# The borel.tannerMlink transformation and the first two derivatives.
print(rbind(a.par, 
 deriv1 = borel.tannerMlink(theta = a.par, inverse = FALSE, deriv = 1),
 deriv2 = borel.tannerMlink(theta = a.par, inverse = FALSE, deriv = 2)),
 digits = 2)
  
# The inverse of 'borel.tannerMlink()' and the first two derivatives.
# 'theta' turns into 'eta'.
print(rbind(a.par, 
    Invderiv1 = borel.tannerMlink(theta = a.par, inverse = TRUE, deriv = 1),
    Invderiv2 = borel.tannerMlink(theta = a.par, inverse = TRUE, deriv = 2)),
    digits = 2)

  
## Example 2 ##
a.param <- c(0, 1, 5, 10, 1e2, 1e3) 
rbind(a.values = a.param, 
      inv.BT = borel.tannerMlink(theta = a.param, inverse = TRUE))
    
data.inv <- borel.tannerMlink(borel.tannerMlink(a.param, inv = TRUE)) - a.param
summary(data.inv)                      ## Should be zero



## Example 3. Some link functions in VGAM with domain set (0, 1) ##
a.param <- ppoints(100)

par(lwd = 2)
plot(a.param, borel.tannerMlink(a.param), ylim = c(-5, 7), xlim = c(-0.01, 1.01),
     type = "l", col = "gray10", ylab = "transformation",
     las = 1, main = "Some probability link functions")
lines(a.param, logffMlink(a.param), col = "blue")
lines(a.param, logitlink(a.param), col = "limegreen")
lines(a.param, probitlink(a.param), col = "purple")
lines(a.param, clogloglink(a.param), col = "chocolate")
lines(a.param, cauchitlink(a.param), col = "tan")
abline(v = c(0.5, 1), lty = "dashed")
abline(v = 0, h = 0, lty = "dashed")
legend(0.05, 7, c("borel.tanneMlink", "logffMlink", "logitlink", "probitlink",
                  "clogloglink", "cauchitlink"),
       col = c("gray10", "blue", "limegreen", "purple", "chocolate", "tan"),
      lwd = 1)
par(lwd = 1)

---

Names/Value of linear predictors/parameters in time series family functions.

Description

Splitting out the names of linear predictors or Numeric values for parameters in time series family functions in VGAMextra.

Usage

break.VGAMextra(eta      = NULL,
                      M1       = NULL,
                      noInter  = NULL,
                      bOrder   = NULL,
                      NOS      = NULL,
                      lInter   = "identitylink",
                      lvar     = "loglink",
                      lsd      = "loglink",
                      lcoeff1  = "rhobitlink",
                      lcoeff2  = "rhobitlink",
                      typeTS   = "AR",
                      namesLP  = FALSE,
                      Complete = FALSE,
                      varArg   = NULL)

Arguments

eta	A matrix of dimensions c(n, M) storing the linear predictors values coming from the vglm fit. Here, $M$ is the number of parameters. See warning below for further information.
M1	Number of parameters involved in the vglm fit.
noInter	Logical. To determine whether the intercept is estimated. If 'TRUE', the intercept is not estimated and set to 0.
bOrder	A vector. The order of the linear process fitted. Either a single number (if one response), or a vector (if multiple responses).
NOS	Integer. Number of respones set in the vglm call.
lInter, lvar, lsd, lcoeff1, lcoeff2	Link functions applied to parameters. Same as in ARXff, or MAXff.
typeTS	Character. Currently, options "AR" for Autoregressive, and "MA" for Moving Average processes are handled.
namesLP	Logical. This function returns either the names of linear the predictors/parameters ( if namesLP = TRUE ) or parameter values (default) broken down from the eta matrix.
Complete	Logical. If TRUE, columns of zeros are incorporated into the matrix eta. See below for further details.
varArg	Sames as in ARXff or MAXff

Details

Time series family functions in VGAMextra currently recycle the order set in the vglm. Particularly, it occurs when the number of responses is fewer than the specified order. For instance, if the order set in vglm is $c(1, 3)$, and 5 responses are managed, then the new order becomes $c(1, 3, 1, 3, 1)$.

Due to such flexibility, time series family functions require specific functions to unload the amount of code within each one.

Moreover, when the order is recycled, the matrix eta is completed, as if the order was the same for each response. This feature is enabled when Complete = TRUE. This ‘common’ order turns out to be the maximum order established in the vector order. This trick makes the family function to work properly. To return to the riginal ‘order’, eta is reduced in the same number of colums initially added.

break.VGAMextra works in this context. It may return either the names of the linear predictors/parameters, or the parameter values splitted out as a list. Thus, link functions entered in the vglm call must be passed down to this functions. For further details on link functions refer to CommonVGAMffArguments.

Value

A list containing either the names of the linear predictors or the parameters values (not linear predictors) unwrapped from tje eta matrix, as follows:

a) If namesLP = FALSE (default), value of parameters are returned in this order: the intercept (1), standard deviation and variance of the white noise (2, 3), and the coefficients (4).

b) If namesLP = TRUE, names of linear predictors are returned in the first entry, whereas parameter names are allocated to the second entry.

Yee and Wild (1996) provide more detailed information about the relationship between linear predictors and parameters within the VGLM statistical framework.

Warning

Note that library VGAM is definitely required.

Warning

Be aware of the dimensions of matrix eta. It is c(n, M), where $n$ is the sample size, and $M$ is the number of parameters. If multiple responses, then $M$ equals the summation of parameters individually.

Author(s)

Victor Miranda and T. W. Yee

References

Yee, T. W. and Wild, C. J. (1996) Vector Generalized Additive Models. Journal of the Royal Statistical Society, Series B, Methodological, 58(3), 481–493.

See Also

ARXff, MAXff, CommonVGAMffArguments, vglm.

Examples

library(VGAM)

eta     <- matrix(runif(100), nrow = 10, ncol = 10)
M1      <- c(5, 5)
noInter <- FALSE
bOrder  <- c(3, 3)
NOS     <- 2

### ONLY LINEAR PREDICTORS/PARAMETERS NAMES!
### RETURNED OBJECT IS A LIST !

break.VGAMextra(M1      = M1, 
                noInter = noInter, 
                bOrder  = bOrder, 
                NOS     = NOS, 
                typeTS  = "AR", 
                namesLP = TRUE, 
                varArg  = TRUE)

### PARAMETER VALUEs... "UNWRAPPED". Inverse link functions are applied.
###  Note that namesLP must be set to FALSE

break.VGAMextra(eta     = eta,
                M1      = M1, 
                noInter = noInter, 
                bOrder  = bOrder, 
                NOS     = NOS, 
                typeTS  = "AR", 
                namesLP = FALSE, 
                varArg  = TRUE)

---

Polynomial roots based on transfer operators in Vector Generalized Time Series Family Functions

Description

checkTS.VGAMextra computes the polynomial roots as per transfer operator in Vector Generalized Time Series Family Functions in VGAMextra

Usage

checkTS.VGAMextra(thetaEst = NULL, 
                          tsclass  = c("AR", "MA"), 
                          chOrder  = 1,
                          NofS     = 1,
                          retmod   = TRUE,
                          pRoots   = TRUE)

Arguments

thetaEst	A vector of coefficients. Its lenght must be NofS * chOrder. If vglm is called, then thetaEst contains the estimated coefficients of the model specified in formula .
tsclass	Character indicating the model class to be checked. Presently, options "AR" and "MA" are handled.
chOrder	Positive integer. The order of polynomial associated to the underlying procees involved: either $p$ or $q$, which apply for "AR" or "MA" rspectively.
NofS	A positive integer denoting the number of Time Series to verify. In the vglm environment, NofS is the number of responses given in formula.
retmod	Logical. If TRUE (default), the Module of all roots as per transfer operator in the process established in tsclass is returned. Else, essentially the roots are returned.
pRoots	Logical. If TRUE (default), the roots computed from estimated models are displayed along with the time series family function execution.

Details

Stationarity and/or Invertibility of time series (TS) are usually verified via the roots of the polynomial derived from the transfer operators.

In particular, checkTS.VGAMextra computes such roots via the coefficients estimated by vector generalized TS family functions available in VGAMextra ( ARXff, and MAXff).

Specifically, checkTS.VGAMextra verifies whether the TS analyzed via vglm is stationary or invertible, accordingly.

Note that an autoregressive process of order-$p$ [AR($p$)] with coefficients $\theta_{1}, \ldots, \theta_{p}$ can be written in the form

$\theta(B) Y_{t} = \varepsilon_{t},$

where

$\theta(B) = 1 - \sum_{k = 1}^{p} \theta_{k} B^{k}$

Here, $\theta(B)$ is referred to as the transfer operator of the process, and $B^{k} Y_{t} = Y_{t - k},$, for $k = 0, 1, \ldots,p$, is the lagged single-function.

In general, an autoregressive process of order-$p$ is stationary if the roots of

$\theta(z) = 1 - \theta_{1} z - \ldots - \theta_{q} z^q$

lie outside the unit circle, i.e. $|z| > 1$.

Similarly, a moving-average process of order-q can be formulated (without loss of generality $\mu = 0$)

$Y_{t} = \psi(B) \varepsilon_{t},$

where $\psi(B)$ is the transfer operator, given by

$\psi(B) = 1 + \sum_{k = 1}^{q} \psi_{k} B^{k},$

Note that $\psi_{0} = 1$, and $B^{k} \varepsilon_{t} = \varepsilon{t - k}$.

Hence, a moving-average process of order-$q$ [MA($q$)], generally given by (note $\mu = 0$)

$Y_{t} = \phi_1 \varepsilon_{t - 1} + \ldots + \phi_q \varepsilon_{t - q} + \varepsilon_{t},$

is invertible if all the roots of

$\phi(B) = 1 + \phi_{1} B + \ldots + \phi_{q} B^q$

lie outside the unit circle., i.e.m $|z| > 1$.

Parallel arguments can be stated for autoregressive moving aberage processes (ARMA). See Box and Jenkins (1970) for further details.

Value

A vector whose elements are the roots of polynomials inherited from transfer operators according to the process analyzed.

Alternatively, the modules of roots can by returned instead of merely roots via the retmod argument.

Warning

The argument thetaEst manages the coefficients of the TS model(s) in turn. Then, it must be NofS * chOrder length, where NofS is the number of responses established in the vglm call. Here the coefficients for each response must be sequentially groped.

A moving average process is always stationary (See Madsen (2007) for further details). Consequently, the MAXff in VGAMextra verifies (by default) only for invertibility. To enable this option set nowarning = FALSE in the MAXff call.

Similarly, ARXff verifies whether the TS data fitted is stationary, whereas ARMAXff() verifies both properties.

Note

For TS family functions in the VGLM/VGAM framework, checkTS.VGAMextra is called at the final iteration of Fisher scoring. It means that the MLE estimates are actually evaluated to verify whether the process is stationary or invertible.

If any root has module less than $1 + 1e-5$, a warning is displayed for informative purposes.

Argument thetaEst manages the parameters of the TS model in turn.

Author(s)

Victor Miranda and T. W. Yee.

References

Box, G.E.P. and Jenkins, G.M. (1970) Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco, USA.

Madsen, H (2007) Time Series Analysis. Chapman and Hall/CRC, Boca Raton, Florida, USA.

Examples

# A moving average process order-3 with coeffs --> c(2.4, -5.6, 0.83)
#-------------------------#
# This is NOT invertible !
#-------------------------#

MAcoeffs <- c(2.4, -5.6, 0.83)
checkTS.VGAMextra(thetaEst = MAcoeffs, 
                  tsclass = "MA", 
                  chOrder = 3,
                  retmod = FALSE)


# AR process order-3 with coeffs --> c( 0.45, 0.45^2, 0.45^3 )
#-------------------------#
# This is stationary !
#-------------------------#

ARcoeffs <- c( 0.45 , 0.45^2 , 0.45^3 )
checkTS.VGAMextra(thetaEst = ARcoeffs, 
                  tsclass = "AR", 
                  chOrder = 3,
                  retmod = TRUE,
                  pRoots = TRUE)  # DEFAULT for 'pRoots'

---

Constraint matrices for vector generalized time series family functions.

Description

Constraint matrices for coefficients of vector genelized time series family functions in VGAMextra.

Usage

cm.ARMA(Model      = ~ 1, 
               Resp       =  1,
               lags.cm    =  2,
               offset     = -2, 
               whichCoeff =  1,
               factorSeq  =  2)

Arguments

Model	A symbolic description of the model being fitted. Must match that formula specified in the vglm() call.
lags.cm	Vector of POSITIVE integers greater than 1 indicating the order for each response. It must match the orders entered in the vglm call. Its default value is 2, assuming that a TS process of order greater than 1 is being fitted. If lags.cm < 2, then NO constraints are required as only one coefficient (AR, MA or ARMA) is being estimated.
offset	Vector of integers specifying the position of the ARMA coefficient at which constraints initiate FOR each response. If negative, it is recycled and the absolute value is used. The default value is -2, which refers to the fourth position on the vector parameter, right after the drift or mean, the white noise sd, and the first ARMA coefficient. Particularly, if only one coefficient is being estimated, i.e, an AR or MA process of order-1 is being fitted, then NO restrictions over the (unique) coefficient are needed. Consequently, abs(offset < 2) leads to a message error.
whichCoeff	Vector of POSITIVE integers strictly less than 'abs(offset)', each entry aplies to each response in the vglm(...) call. This argument allows the user to specify the unrestricted coefficient to be considered for constraints. For instance, whichCoeff = 2 means that $\theta_2$ is the required coefficient to compute the constraint matrices. By default, whichCoeff = -1 which implies that $\theta_1$ is used for this purpose. If whichCoeff is greater than or equal to abs(offset), an error message is displayed since constraints must be function of unrestricted parameters.
Resp	The number of responses in the model fitted. Must match the number of responses given in formula in the vglm call.
factorSeq	Vector of POSITIVE integers. Thus far, restrictions handled are geometric sequences and arithmetic progressions. Hence, factorSeq specifies either the initial power or factor at restrictions. See below for further details.

Details

NOTE: Except for the Model, all arguments of length 1 are recycled when Resp $\geq 2$.

Time Series family functions in VGAMextra that are derived from AR(p) or MA(q) processes include the drift term (or mean) and the white noise standard deviation as the first two elements in the vector parameter. For an MA(4), for example, it is given by

$(\mu, \sigma_\varepsilon, \phi_1, \phi_2, \phi_3, \phi_4).$

Thus, constraint matrices on coefficients can be stated from the second coefficient, i.e., from $\phi_2$. This feature is specified with offset = -2 by default.

In other words, offset indicates the exact position at which parameter restrictions commence. For example, offset = -3 indicates that $\phi_3$ is the first coefficient over which constraints are applied. Then, in order to successfully utilize this argument, it must be greater than or equal to 2 in absolute value. Otherwise, an error message will be displayed as no single restriction are amenable with $\phi_1$ only.

Furthermore, if lags.cm = 1, i.e, a AR or MA process of order one is being fitted, then NO constraints are required either, as only one coefficient is directly considered.

Hence, the miminum absolute value for argument offset is 2

As for the factorSeq argument, its defaul value is 2. Let factorSeq = 4, lags.cm = 5, offset = -3, and whichCoeff = 1. The coefficient restrictions if a geometric progression is assumed are

$\theta_3 = \theta_1^4,$

$\theta_4 = \theta_1^5,$

$\theta_5 = \theta_1^6,$

If coefficient restrictions are in arithmetic sequence, constraints are given by

$\theta_3 = 4 * \theta_1,$

$\theta_4 = 5 * \theta_1,$

$\theta_5 = 6 * \theta_1,$

The difference lies on thelink function used: loglink for the first case, and identitylink for the latter.

Note that conditions above are equivalent to test the following two Null Hypotheses:

$Ho: \theta_k = \theta_1^k$

or

$Ho: \theta_k = j * \theta_1$

for $k = 3, 4, 5$.

Simpler hypotheses can be tested by properly setting all arguments in cm.ARMA(). For instance, the default list of constraint matrices returned by cm.ARMA() allows to test

$Ho: \theta_k = \theta_1^j$

for $k = 2$, in a TS model of order-2 with one response.

Value

A list of constraint matrices with specific restrictions over the AR($p$), MA($q$) or ARMA ($p, q$) coefficients. Each matrix returned is conformable with the VGAM/VGLM framework.

Paragrpah above means that each constraint matrix returned by cm.ARMA() is full-rank with M rows (number of parameters), as required by VGAM. Note that constraint matrices within the VGAM/VGLM framework are M by M identity matrices by default.

Restrictions currently handled by cm.ARMA() are (increasing) arithmetic and geometric progressions.

Warning

Hypotheses above can be tested by properly applying parameter link functions. If the test

$Ho: \theta_k = \theta_1^k,$

arises, then constraint matrices returned by cm.ARMA() are conformable to the use of loglink.

On the other hand, the following hypothesis

$Ho: \theta_k = k * \theta_1,$

properly adapts to the link function identitylink. $k = 2, 3,\ ldots$.

For further details on parameter link functions within VGAM, see CommonVGAMffArguments.

Note

cm.ARMA() can be utilized to compute constraint matrices for many VGLTSM fmaily functions, e.g., ARXff and MAXff in VGAMextra.

More improvements such as restrictions on the drift parameter and white noise standard deviation will be set later.

Author(s)

Victor Miranda and T. W. Yee

References

Yee, T. W. and Hastie, T. J. (2003) Reduced-rank vector generalized linear models. Statistical Modelling, 3, 15–41.

Yee, T. W. (2008) The VGAM Package. R News, 8, 28–39.

See Also

loglink, rhobitlink, CommonVGAMffArguments.

Examples

#############
# Example 1.
#############
# Constraint matrices for a TS family function (AR or MA) 
# with 6 lagged terms.
# Restriction commences at the third position (theta[3]) powered to
# or multiplied by 4. Intercept-only model.
position   <- -3
numberLags <-  6
myfactor   <-  4
cm.ARMA(offset = position, lags.cm = numberLags, factorSeq = myfactor)

# With one covariate
cm.ARMA(Model =  ~ x2, offset = position, 
        lags.cm = numberLags, factorSeq = myfactor)


# Or 2 responses...
cm.ARMA(offset = position, lags.cm = numberLags, 
        factorSeq = myfactor, Resp = 2)


# The following call causes an ERROR.
# cm.ARMA(offset = -1, lags.cm = 6, factorSeq = 2)


##############
# Example 2.
##############

# In this example, the use of constraints via 'cm.ARMA()' is
# included in the 'vglm' call. Here, two AR(2) models are fitted
# in the same call (i.e. two responses), where different constraints
# are set, as follows:
# a) list(ar = c(theta1, theta1^2)) and
# b) list(ar = c(theta2, theta2^2 )).

# 2.0 Generate the data.
set.seed(1001)
nn     <- 100
# A single covariate.
covdata <- data.frame(x2 =  runif(nn)) 

theta1 <- 0.40; theta2 <- 0.55
drift  <- c(0.5, 0.75)
sdAR   <- c(sqrt(2.5), sqrt(2.0))

# Generate AR sequences, TS1 and TS2, considering Gaussian white noise

# Save both in a data.frame object: the data.
tsdata  <- 
  data.frame(covdata, # Not used 
             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(theta2, theta2^2)), 
                             rand.gen = rnorm, 
                             mean = drift[2], sd = sdAR[2]))

# 2.1 Fitting both time series with 'ARXff'... multiple responses case.
fit1 <- vglm(cbind(TS1, TS2) ~ 1, 
             ARXff(order = c(2, 2), type.EIM = "exact"), 
             data = tsdata,  
             trace = TRUE)

Coef(fit1)                
coef(fit1, matrix = TRUE)
summary(fit1)

## Same length for both vectors, i.e. no constraints.
length(Coef(fit1))
length(coef(fit1, matrix = TRUE))




###2.2 Now, fit the same models with suitable constraints via 'cm.ARMA()'
# Most importantly, "loglink" is used as link function to adequately match 
# the relationship between coefficients and constraints. That is:
# theta2 = theta1^2, then log(theta2) = 2 * log(theta1).

fit2 <- vglm(cbind(TS1, TS2) ~ 1, 
             ARXff(order = c(2, 2), type.EIM = "exact", lARcoeff = "loglink"), 
             constraints = cm.ARMA(Model = ~ 1, 
                                   Resp = 2,
                                   lags.cm = c(2, 2),
                                   offset  = -2),
             data = tsdata,  
             trace = TRUE)
Coef(fit2)
coef(fit2, matrix = TRUE)
summary(fit2)

# NOTE, for model 1, Coeff2 = Coeff1^2, then log(Coeff2) = 2 * log(Coeff1)
( mycoef <- coef(fit2, matrix = TRUE)[c(3, 4)] )
2 * mycoef[1] - mycoef[2]    # SHOULD BE ZERO

# Ditto for model 2: 
( mycoef <- coef(fit2, matrix = TRUE)[c(7, 8)] )
2 * mycoef[1] - mycoef[2]    # SHOULD BE ZERO

## Different lengths, due to constraints
length(Coef(fit2))
length(coef(fit2, matrix = TRUE))

---

Density for the multivariate Normal distribution

Description

Density for the multivariate Normal distribution

Usage

dmultinorm(vec.x, vec.mean = c(0, 0),
                 mat.cov = c(1, 1, 0),
                 log = FALSE)

Arguments

vec.x	For the $R$–multivariate Normal, an $R$–vector of quantiles.
vec.mean	The vector of means.
mat.cov	The vector of variances and covariances, arranged in that order. See below for further details.
log	Logical. If TRUE, the logged values are returned.

Details

This implementation of the multivariate (say $R$–dimensional) Normal density handles the variances and covariances, instead of correlation parameters.

For more than one observation, arrange all entries in matrices accordingly.

For each observation, mat.cov is a vector of length $R \times (R + 1) / 2$, where the first $R$ entries are the variances $\sigma^2{i}$, $i = 1, \ldots, R$, and then the covariances arranged as per rows, that is, $cov_{ij}$ $i = 1, \ldots, R, j = i + 1, \ldots, R$.

By default, it returns the density of two independent standard Normal distributions.

Value

The density of the multivariate Normal distribution.

Warning

For observations whose covariance matrix is not positive definite, NaN will be returned.

Author(s)

Victor Miranda

See Also

binormal.

Examples

###
### Two - dimensional Normal density.
###
set.seed(180228)
nn  <- 25
mean1 <- 1; mean2 <- 1.5; mean3 = 2
var1 <- exp(1.5); var2 <- exp(-1.5); var3 <- exp(1); cov12 = 0.75
dmvndata <- rbinorm(nn, mean1 = 1, mean2 = 1.5, var1 = var1, var2 = var2,
                    cov12 = cov12)

## Using dbinorm() from VGAM.
d2norm.data <- dbinorm(x1 = dmvndata[, 1], x2 = dmvndata[, 2],
                        mean1 = mean1, mean2 = mean2, var1 = var1, var2 = var2,
                        cov12 = cov12)
## Using dmultinorm().
d2norm.data2 <- dmultinorm(vec.x = dmvndata, vec.mean = c(mean1, mean2),
                        mat.cov = c(var1, var2, cov12))
summary(d2norm.data)
summary(d2norm.data2)
##
## 3--dimensional Normal.
##
dmvndata <- cbind(dmvndata, rnorm(nn, mean3, sqrt(var3)))

d2norm.data3 <- dmultinorm(dmvndata, vec.mean = c(mean1, mean2, mean3),
                       mat.cov = c(var1, var2, var3, cov12, 0, 0))

hist(d2norm.data3)
summary(d2norm.data3)

---

VGLTSM family function for the Two–dimensional Error–Correction Model (Engle and Granger, 1987) for $I(1)$–variables

Description

Estimates a bidimensional error-correction model of order–($K$, $L$), as proposed by Engle–Granger (Two step–approach; 1987), with bivariate normal errors by maximum likelihood estimation using Fisher scoring.

Usage

ECM.EngleGran(ecm.order = c(1, 1),
                    zero = c("var", "cov"),
                    resids.pattern = c("intercept", "trend",
                                       "neither", "both")[1],
                    lag.res = 1, 
                    lmean = "identitylink",
                    lvar  = "loglink",
                    lcov  = "identitylink",
                    ordtsDyn = 0)

Arguments

ecm.order	Length–2 (positive) integer vector. The order of the ECM model.
zero	Integer or character–string vector. Details at zero.
resids.pattern	Character. How the static linear regression $y_{2, t} \sim y_{1, t}$ must be settle to estimate the residuals $\widehat{z_t}$. The default is a linear model with intercept, and no trend term. See below for details.
lag.res	Numeric, single positive integer. The error term for the long–run equilibrium path is lagged up to order lag.res. See below for further details.
lmean, lvar, lcov	Same as MVNcov.
ordtsDyn	Positive integer. Allows to compare the estimated coefficients with those provided by the package 'tsDyn'. See below for further details.

Details

This is an implementation of the two–step approach as proposed by Engle–Granger [1987] to estimate an order–($K$, $L$) bidimensional error correction model (ECM) with bivariate normal errors.

This ECM class models the dynamic behaviour of two cointegrated $I(1)$-variables, say $y_{1, t}$ and $y_{2, t}$ with, probably, $y_{2, t}$ a function of $y_{1, t}$. Note, the response must be a two–column matrix, where the first entry is the regressor, i.e, $y_{1, t}$ above, and the regressand in the second colum. See Example 2 below.

The general specification of the ECM class described by this family function is

$\Delta y_{1, t} |\Phi_{t - 1} = ~\phi_{0, 1} + \gamma_1 \widehat{z}_{t - k} + \sum_{i = 1}^K \phi_{1, i} \Delta y_{2, t - i} + \sum_{j = 1}^L \phi_{2, j} \Delta y_{1, t - j} + \varepsilon_{1, t},$

$\Delta y_{2, t} |\Phi_{t - 1}= ~\psi_{0, 1} + \gamma_2 \widehat{z}_{t - k} + \sum_{i = 1}^K \psi_{1, i} \Delta y_{1, t - i} + \sum_{j = 1}^L \psi_{2, j} \Delta y_{2, t - j} + \varepsilon_{2, t}.$

Under the binormality assumption on the errors $(\varepsilon_{1, t}, \varepsilon_{2, t})^T$ with covariance matrix $\boldsymbol{\textrm{V}}$, model above can be seen as a VGLM fitting linear models over the conditional means, $\mu_{\Delta y_{1, t} } = E(\Delta y_{1, t} | \Phi_{t - 1})$ and $\mu_{\Delta y_{2, t} } = E(\Delta y_{2, t} | \Phi_{t - 1} )$, producing

$(\Delta y_{1, t} |\Phi_{t - 1} , \Delta y_{2, t} |\Phi_{t - 1} )^T \sim N_{2}(\mu_{\Delta y_{1, t}} , \mu_{\Delta y_{2, t}}, \boldsymbol{\textrm{V}})$

The covariance matrix is assumed to have elements $\sigma_1^2, \sigma_2^2,$ and $\textrm{Cov}_{12}.$

Hence, the parameter vector is

$\boldsymbol{\theta} = (\phi_{0, 1}, \gamma_1, \phi_{1, i}, \phi_{2, j}, \psi_{0, 1}, \gamma_2, \psi_{1, i}, \psi_{2, j}, \sigma_1^2, \sigma_2^2, \textrm{Cov}_{12})^T,$

for $i = 1, \ldots, K$ and $j = 1, \ldots, L$.

The linear predictor is

$\boldsymbol{\eta} = (\mu_{\Delta y_{1, t} }, \mu_{\Delta y_{2, t} }, {\color{blue}\texttt{loglink}}~\sigma_1^2, {\color{blue}\texttt{loglink}}~\sigma_2^2, \textrm{Cov}_{12})^T.$

The estimated cointegrated vector, $\boldsymbol{\widehat{\beta^{\star}}}$ = $(1, -\boldsymbol{\widehat{\beta})}^T$ is obtained by linear regression depending upon resids.pattern, as follows:

1) $y_{2, t} = \beta_0 + \beta_1 y_{1, t} + z_t$, if resids.pattern = "intercept",

2) $y_{2, t} = \beta_1 y_{1, t} + \beta_2 t + z_t$, if resids.pattern = "trend",

3) $y_{2, t} = \beta_1 y_{1, t} + z_t$, if resids.pattern = "neither", or else,

4) $y_{2, t} = \beta_0 + \beta_1 y_{1, t} + \beta_2 t + z_t$, if resids.pattern = "both",

where $\boldsymbol{\widehat{\beta}} = (\widehat{\beta_0}, \widehat{\beta_1}, \widehat{\beta_2})^T,$ and $z_t$ assigns the error term.

Note, the estimated residuals, $\widehat{z_t}$ are (internally) computed from any of the linear models 1) – 4) selected, and then lagged up to order alg.res, and embedded as explanatories in models $\Delta y_{1, t} |\Phi_{t - 1}$ and $\Delta y_{3, t} |\Phi_{t - 1}$ above. By default, $\widehat{z}_{t - 1}$ are considered (as lag.res = 1), although it may be any lag $\widehat{z}_{t - k}$, for $k > 0$. Change this through argument lag.res.

Value

An object of class "vglmff" (see vglmff-class) to be used by VGLM/VGAM modelling functions, e.g., vglm or vgam.

Note

Reduced–Rank VGLMs (RR-VGLMs) can be utilized to aid the increasing number of parameters as $K$ and $L$ grows. See rrvglm.

By default, $\sigma_1^2, \sigma_2^2$ and $\textrm{Cov}_{12}$ are intercept–only. Set argument zero accordingly to change this.

Package tsDyn also has routines to fit ECMs. However, the bivariate–ECM handled (similar to that one above) differs in their parametrization: tsDyn considers the current estimated residual, $\widehat{z}_t$ instead of $\widehat{z}_{t - 1}$ in models $\Delta y_{1, t} |\Phi_{t - 1}$ and $\Delta y_{2, t} |\Phi_{t - 1}$.

See Example 3 below which compares ECMs fitted with VGAMextra and tsDyn.

Author(s)

Victor Miranda

References

Engle, R.F. and Granger C.W.J. (1987) Co-integration and error correction: Representation, estimation and testing. Econometrica, 55(2), 251–276.

Pfaff, B. (2011) Analysis of Integrated and Cointegrated Time Series with R. Seattle, Washington, USA: Springer.

See Also

MVNcov, rrvglm, CommonVGAMffArguments, Links, vglm.

Examples

## Example 1. Comparing the Engle -- Granger procedure carried oud by two procedures.
##            ECM.EngleGran() makes easier the fitting process.
## Here, we will use:
## A) The R code 4.2, in Chapter 4, Pfaff (2011).
##    This code 1) generates artificial data and 2) fits an ECM, following
##    the Engle --Granger procedure. 
## B) The ECM.EngleGran() family function to fit the same model assuming
##    bivariate normal innovations. 
## The downside in the R code 4.2 is the assumption of no--correlation among
## the errors. These are generated indenpendently.
## A)
## STEP 1. Set up the data (R code as in Pfaff (2011)).
nn <- 100
set.seed(123456)
e1 <- rnorm(nn)   # Independent of e2
e2 <- rnorm(nn)
y1 <-  cumsum(e1)
y2 <- 0.6 * y1 + e2
lr.reg <- lm(y2 ~ y1)
error <- residuals(lr.reg)
error.lagged <- error[-c(nn - 1, nn)]
dy1 <- diff(y1)
dy2 <- diff(y2)
diff.dat <- data.frame(embed(cbind(dy1, dy2), 2))
colnames(diff.dat) <- c('dy1', 'dy2', 'dy1.1', 'dy2.1')

##  STEP 2. Fit the ECM model, using lm(), R code as in Pfaff (2011).
ecm.reg <- lm(dy2 ~ error.lagged + dy1.1 + dy2.1, data = diff.dat)

summary(ecm.reg)


## B) Now, using ECM.EngleGran() and VGLMs, the steps at A) can be skipped. 
## Enter the I(1)--variables in the response vector only, putting down the
## the dependent variable from the I(1) set, i.e. y2, in the second column.

coint.data <- data.frame(y1 = y1, y2 = y2)
fit.ECM <- vglm(cbind(y1, y2) ~ 1, ECM.EngleGran, data = coint.data, trace = TRUE)

## Check coefficients ##
coef(fit.ECM, matrix = TRUE)  ## Compare 'Diff2' with summary(ecm.reg)
coef(summary(ecm.reg))

head(depvar(fit.ECM))   # The estimated differences (first order)
vcov(fit.ECM)
constraints(fit.ECM, matrix = TRUE)

## Not run: 
### Example 2.  Here, we compare ECM.EngleGran() from VGAMextra with VECM() from
##              package "tsDyn" when fitting an ECM(1, 1). We will make use of
##              the argument 'ordtsDyn' so that the outcomes can be compared.

library("tsDyn")  # Need to be installed first.
fit.tsDyn1 <- with(coint.data, VECM(cbind(y2, y1), lag = 1, estim = "2OLS")) #  MODEL 1
summary(fit.tsDyn1)

### Fit same model using ECM.EngleGran(). NOTE: Set ordtsDyn = 1 !!          #  MODEL 2
fit.ECM.2 <- vglm(cbind(y1, y2) ~ 1, ECM.EngleGran(ecm.order = c(1, 1),
                  resids.pattern = "neither", ordtsDyn = 1),
                  data = coint.data, trace = TRUE)

coef.ECM.2 <- coef(fit.ECM.2, matrix = TRUE)
fit.tsDyn1$coefficients                      ## From pakage 'tsDyn'.
t(coef.ECM.2[, 1:2][c(2, 1, 4, 3), ][, 2:1]) ## FROM VGAMextra 


### Example 3. An ECM(2, 2), with residuals estimated by OLS, with NO intercept
###            and NO trend term. The data set is 'zeroyld', from package tsDyn.
###            ECM.EngleGran() and with VECM() willbe compared again.
data(zeroyld, package = "tsDyn")

# Fit a VECM with Engle-Granger 2OLS estimator:
vecm.eg <- VECM(zeroyld, lag=2, estim = "2OLS") 
summary(vecm.eg)

# For the same data, fit a VECM with ECM.EngleGran(), from VGAMextra.
# Set ordtsDyn = 1 for compatibility! 
fit.ECM.3 <- vglm(cbind(long.run, short.run) ~ 1, ECM.EngleGran(ecm.order = c(2, 2),
                                  resids.pattern = "neither", ordtsDyn = 1),
                  data = zeroyld, trace = TRUE)
coef.ECM.3 <- coef(fit.ECM.3, matrix = TRUE)

#### Compare results
vecm.eg$coefficients                               # From tsDyn
t(coef.ECM.3[, 1:2][c(2, 1, 5, 3, 6, 4 ),][, 2:1]) # FROM VGAMextra

## End(Not run)

---