R/mat.polyn.R

##' @title pol.inv
##' 
##' @description Calculation of left inverse of matrix polynomial. The
##' leading term
##' is expected to be the (k by k) identity matrix. This is checked
##' and the proper leading unity term is taken into account when the
##' inverse is calculated.
##'
##' phi = matrix polynomial coefficients = I, phi1, phi2, ..., phi(p).
##' 
##' dim(phi) = c(k, k, p+1) where k = dimension of coefficient
##' matrices (k by k), and L  = order of polynomial (length = 1+L ,
##' including the leading unity matrix).
##'
##' @param phi polynomium (an array) to invert
##' @param L  order of inverse polynomium
##'
##' @return  left inverse of phi of order L (L+1 terms including leading
##' unity matrix)
##'
##' @examples
##' set.seed(4711)
##' p2<-check.one(array(rnorm(32),dim=c(4,4,2)))
##' pi2<-pol.inv(p2,L=12)
##' short.form(pi2)
##'
##' @export

pol.inv <- function(phi, L) {
   
  "[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)
    aphi <- check.one(phi)
    k <- dim(aphi)[1]
    # if (dim(phi)[3]==dim(aphi)[3]) {up <- L+1}
    up <- L + 1
    ip <- dim(aphi)[3]
    polinv <- array(0, dim = c(k, k, up))
    polinv[, , 1] <- diag(k)
    if (up > 1) {
        if (ip > 1) {
            for (LL in 2:up) {
                for (J in 2:LL) {
                   if (J <= ip) {
      help <- matrix(polinv[, , LL + 1 - J], nrow=k ) %*%
              matrix(aphi[, , J], nrow = k )
      polinv[, , LL] <- matrix(polinv[, , LL], nrow = k) -
                        matrix(help, nrow = k)
                   }
                }
            }
        }
    }
    return(polinv)
}


##' @title pol.mul
##' 
##' @description Calculation of product of two matrix
##' polynomials (arrays).
##'
##' If one or both leading unity matrices (of eta and theta) are
##' missing, they are (it is)
##' generated (and taken into account).
##'
##' @param eta first matrix polynomial
##' @param theta  second matrix olynomial
##' @param L order of output polynomial (length = L+1)
##'
##' @return  matrix polynomial product af eta and theta
##'
##' @examples
##' set.seed(4711)
##' p1 <- check.one(matrix(rnorm(16), nrow=4))
##' p2 <- check.one(array(rnorm(32),dim=c(4, 4, 2)))
##' p12 <- pol.mul(p1, p2, L=(2+3))
##' short.form(p12)
##'
##' @export

pol.mul <- function(eta, theta, L) {
    
  "[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)
    eta <- check.one(eta)
    theta <- check.one(theta)
    if (L < 1) {
        L <- dim(eta)[3] + dim(theta)[3] - 1
    }
    L <- L + 1
    k <- dim(eta)[1]
    iet <- dim(eta)[3]
    ith <- dim(theta)[3]
    polmul <- array(0, dim = c(k, k, L))
    for (LL in 1:L) {
        polmul[, , LL] <- 0
        for (J in 1:LL) {
            if (J <= iet) {
                if (LL + 1 - J <= ith) {
        help <- matrix(eta[, , J], nrow = k) %*%
            matrix(theta[, , LL + 1 - J], nrow = k)
        polmul[, , LL] <- matrix(polmul[, , LL], nrow = k) +
            matrix(help, nrow = k)
                }
            }
        }
    }
    return(polmul)
}

##' @title rand.shock
##' 
##' @description Calculation of random shock form for arma model
##'
##' @param ar.poly autoregressive matrix part of model
##' @param ma.poly moving average matrix part of model
##' @param L order of return polynomial  (length=L+1 including
##' leading unity matrix)
##'
##' @return random shock form of arma model up to order L (array(k,k,L+1))
##'
##' @examples
##' set.seed(4711)
##' p1 <- check.one(matrix(rnorm(16),nrow=4))
##' p2 <- check.one(array(rnorm(32),dim=c(4, 4, 2)))
##' randshock <- rand.shock(ar.poly=p1, ma.poly=p2, L=6)
##' short.form(randshock)
##'
##' @export

rand.shock <- function(ar.poly, ma.poly, L) {
    
 "[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)    
    ar <- check.one(ar.poly)
    ma <- check.one(ma.poly)
    invar <- pol.inv(ar, (L + 1))
    rand.shock <- pol.mul(invar, ma, L)
    return(rand.shock)
}

##' @title inverse.form
##'
##' @description Calculation of inverse form for arma model
##'
##' @param ar.poly =autoregressive matrix part of model (array(k, k, ar-order)).
##' @param ma.poly =moving average matrix part of model (array(k, k, ma-order)).
##' @param L  =order of return polynomial (length=L+1 including leading
##' unity matrix).
##'
##' @return inverse form for arma model up to order L (array(k, k, L+1)).
##'
##' @examples
##' set.seed(4711)
##' p1  <- check.one(matrix(rnorm(16), nrow=4))
##' p2  <- check.one(array(rnorm(32),dim=c(4, 4, 2)))
##' inverse <- inverse.form(ar.poly=p1, ma.poly=p2, L=6)
##' short.form(inverse)
##'
##' @export

inverse.form <- function(ar.poly, ma.poly, L) {

 "[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)
 
    ar <- check.one(ar.poly)
    ma <- check.one(ma.poly)
    inma <- pol.inv(ma, (L + 1))
    inverse.form <- pol.mul(inma, ar, L)
    return(inverse.form)
}

##' @title check.one
##' 
##' @description Function to check and insert leading unity matrix
##' if NOT present.
##'
##' @param polyn (k, k, ...) matrix polynomium with or without leading
##' unity matrix.
##'
##' @return polyn (array) with a leading unity matrix being
##' inserted if not present.
##'
##' @examples
##' set.seed(4711)
##' X <- array(rnorm(32),dim=c(4, 4, 2))
##' X <- check.one(X)
##' short.form(X)
##'
##' @export

check.one <- function(polyn = NULL) {

    "[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)

    if (is.null(polyn)) {
        stop("No input polynomial to 'check.one'. \n")
    }
    d <- dim(polyn)[1]
    dd <- dim(polyn)
    if (is.null(d)) {
        cat("Empty polynomial in call to check.one \n")
    }
    if (length(dim(polyn)) < 2) {
        cat("Argument not polynomial in call to check.one")
        cat(dd, " ", c(polyn), "\n")
    }
    pol2 <- polyn
    if (length(dim(polyn)) < 3) {
        pol2 <- array(polyn, dim = c(d, d, 1))
    }
    
# cat("pol2",pol2,"\n")
#     print(pol2)
# cat("max(abs(pol2[, , 1]))", max(abs(pol2[, , 1])),"\n")
# cat("diag(1,d)",diag(1,d),"\n")

    Dif <- c(pol2[, , 1])-c(diag(1, d))
        if (max(abs(Dif)) > 1e-16) {
        pol2 <- lead.one(pol2, +1)
    }
    pol2[, , 1] <- diag(1, d)  # set leading polyn = unity exactly.
    return(pol2)
}

##' @title lead.one
##' @description Function to add (or remove) a leading unity matrix to (from)
##' an  array
##' (being an array representation of ar- or ma-polynomial).
##' 
##' @param polyn an input polynomium (an array).
##' @param add indicator for adding or removing unity matrix:
##' +1 = add leading unity matrix,
##' -1 = remove leading matrix.
##' @return changed array (with leading unity matrix inserted or
##' removed).
##'
##' @export

lead.one <- function(polyn = NULL, add = 0) {

    "[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)
    #
    # title: lead.one
    # 
    # description:  Function to add or remove leading identity matrix from
    # matrix polynomial. Is used by function 'check.one'.
    #
    # param: polyn matrix polynomial = array of ar- or ma-coefficients
    # param: add +1 or -1 for addition or removal of leading unity matrix
    #
    # return: polyn (an array) with or without a leading unity matrix
    #
    
    d <- dim(polyn)[1]
    dd <- dim(polyn)
    if (is.null(d)) {
        cat("Empty polynomial in call to lead.one,\n")
    }
    if (length(dim(polyn)) < 2) {
        cat("Argument not polynomial in call to lead.one")
        cat(dd, c(polyn), "\n")
    }
    l <- dim(polyn)[3]
    d <- dim(polyn)[1]
    if (add == 1) {
        out.pol <- array(0, dim = c(d, d, (l + 1)))
        out.pol[, , 1] <- diag(1, d)
        out.pol[, , 2:(l + 1)] <- polyn
    }
    if (add == -1) {
        out.pol <- array(0, dim = c(d, d, (l - 1)))
        out.pol[, , 1:(l - 1)] <- polyn[, , 2:l]
    }
    return(out.pol)
}

##' @title pol.order
##' 
##' @description Function to evaluate (significant)
##' order of matrix polynomium.
##'
##' @param polyn the polynomium the order of which is determined.
##' 
##' @param digits number of significant digits to be considered (values
##' smaller than 10^(-digits) are taken to be 0 (zero)).
##'
##' @return pol.order order of polynomium polyn.
##' (exclusive the leading unity matrix if present.
##' pol.order=0 corresponds to the k by k unity matrix)
##'
##' @examples
##' pol          <- array(1e-8*rnorm(96),dim=c(4,4,6))
##' pol[, , 1:3] <- array(rnorm(48), dim=c(4,4,3))
##' pol.order(polyn=pol, digits=12)
##' pol.order(polyn=pol, digits=4)
##'
##' @export

pol.order <- function(polyn = NULL, digits = 12) {
    
     "[" <- function(x, ...) .Primitive("[")(x, ..., drop = FALSE)
    pol.order <- NULL
    if (is.array(polyn)) {
        if (dim(polyn)[1] == dim(polyn)[2]) {
            polyn <- check.one(polyn)
            d <- dim(polyn)
            for (i in 1:d[3]) {
                S <- round(sum(abs(polyn[, , i])), digits = digits)
                if (S > 0) {
                  pol.order <- (i - 1)
                }
            }
        }
    }
    return(pol.order)
}

Try the marima package in your browser

Any scripts or data that you put into this service are public.

marima documentation built on May 2, 2019, 2:10 p.m.