R/ECMEngleGran.R

In VGAMextra: Additions and Extensions of the 'VGAM' Package

##########################################################################
# These functions are
# Copyright (C) 2014-2020 V. Miranda & T. Yee
# Auckland University of Technology & University of Auckland
# All rights reserved.
#
# Links renamed on Jan-2019 conforming with VGAM_1.1-0


ECM.EngleGran.control <- function(save.weights = TRUE,
                                  summary.HDEtest = FALSE,...) { 
  #criterion <- "loglikelihood"
  list(save.weights = save.weights,
       summary.HDEtest = summary.HDEtest,...)
}


ECM.EngleGran <- function(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) {
  
  lmean <- as.list(substitute(lmean))
  emean <- link2list(lmean)
  lmean <- attr(emean, "function.name")
  
  lvar <- as.list(substitute(lvar))
  evar <- link2list(lvar)
  lvar <- attr(evar, "function.name")
  
  loffd <- lcov; rm(lcov)
  loffd <- as.list(substitute(loffd))
  eoffd <- link2list(loffd)
  loffd <- attr(eoffd, "function.name")
  
  V.class  <- c("ECM", "VAR")[1]
  trinorm  <- FALSE
  noRemove <- FALSE
  
  if (!Is.Numeric(ordtsDyn, length.arg =  1) || (ordtsDyn < 0))
    stop("Bad input for 'ordtsDyn'.")
  
  resids.pattern <- match.arg(resids.pattern,  c("intercept", "trend",
                                                 "neither", "both"))[1]
  r.pat <- resids.pattern; rm(resids.pattern)
  V.class <- match.arg(V.class, c("ECM", "VAR"))[1]
  
  if (!Is.Numeric(lag.res, length.arg = 1, 
                  isInteger = TRUE) || lag.res < 1)
    stop("Bad enter for argument 'lag.res'.")
  
  if (!is.logical(noRemove))
    stop("Wrong input for argument 'noRemove'.")
  
  if ((length(my.ord <- ecm.order ) != 2 ) || any(ecm.order <= 0))
    stop("Wrong input for argument 'ecm.order'.")
  
  temp1 <- switch(V.class,
                  ECM = "Error Correction Model (Engle -- Granger)\n",
                  VAR = "Vector Autoregressive Model\n")
  alo <-
    new("vglmff",
      blurb = c(temp1,
                "Links:   ",
                namesof("Cond.Means", lmean, earg = emean), ", ",
                namesof("variance", lvar, earg = evar), ", ",
                namesof("ErrorCovariances", loffd, earg = eoffd)),
      
      
    constraints = eval(substitute(expression({
      M <- M1 <- extra$n.means + ( extra$n.means *
                                     (extra$n.means + 1) )/2
      constraints <- 
        cm.zero.VGAM(constraints, x = x, zero = .zero , M = M ,
                     predictors.names = parameters.names, 
                     M1 = M1)
    }), list( .zero = zero, .trinorm = trinorm ))),
      
      
      
      
    first = eval(substitute(expression({
      if (!length(ncol(y)))
        stop("\n Refer to uninormal() from VGAM to fit the",
             " univariate Normal.")
      
      if (ncol(y) != 2)
        stop("Currently, only bivariate responses handled.")
      
      class.in <- .V.class; res.p <- .r.pat
      NOS.f <- NCOL(y); nn <- NROW(y)
      int.ord <- .my.ord
      T.trend <- 1:nn
      
      if (class.in %in% c("ECM")) {
        my.form <- switch(res.p,
                          "neither" = y[, 2] ~ y[, 1] - 1,
                          "intercept" = y[, 2] ~ y[, 1],
                          "trend" = y[, 2] ~ y[, 1] + T.trend - 1,
                          "both"  = y[, 2] ~ y[, 1] + T.trend )
        int.regre  <- lm(my.form)
        int.resids <- extra$coint.res <- cbind(residuals(int.regre))
        extra$coint.vec <- coef(int.regre)
        names(extra$coint.vec) <- switch(res.p,
                        "neither" = "betaY1",
                        "intercept" = c("(Intercept)", "betaY1"),
                        "trend" = c("betaY1", "T.trend"),
                        "both"  =  c("(Intercept)", "betaY1", "T.trend"))
        rm(T.trend, int.regre)
        
        rem.res <- .lag.res - 1
        temp1   <- if ( !int.ord[1] ) NULL else 
          embed(diff(y[, 1]), 1 + max(int.ord))
        temp2   <- if ( !int.ord[2] ) NULL else
          embed(diff(y[, 2]), 1 + max(int.ord))
        x.temp  <- embed(x, 2)[-c(1:max(int.ord)), 1:NCOL(x), 
                               drop = FALSE]
        
        if ( .ordtsDyn ) {
          temp3 <-
            int.resids[c((.ordtsDyn + 1):(nn - .lag.res - 1 + .ordtsDyn))]
        } else {
          # As Pfaff(2011)
          temp3 <- int.resids[c(1:(nn - .lag.res - 1))]
        }
        
        if (max(int.ord) - .lag.res > 0) {
          temp3 <- temp3[-c(1:(max(int.ord) - .lag.res))]
        } else {
          if (max(int.ord) - .lag.res < 0) {
            temp1 <- temp1[-c(1:abs(max(int.ord) - .lag.res)), ,
                           drop = FALSE]
            temp2 <- temp2[-c(1:abs(max(int.ord) - .lag.res)), , 
                           drop = FALSE]
            x.temp <- x.temp[-c(1:abs(max(int.ord) - .lag.res)), , 
                             drop = FALSE]
          } else {
            
          }
        }
        
        x.matrix <- cbind(x.temp, temp3,
                          if (!int.ord[1]) NULL else 
                            temp1[, 2:(1 + int.ord[1]), drop = FALSE],
                          if (!int.ord[2]) NULL else 
                            temp2[,  2:(1 + int.ord[2]), drop = FALSE])
        
        colnames(x.matrix) <- 
          c(colnames(x), paste("ErrorsLag", .lag.res , sep = ""),
            #paste("dify", 1:NOS.f, sep = ""), 
            if ( !int.ord[1] ) NULL else
              paste("diffy1Lag", 1:int.ord[1], sep = ""),
            if ( !int.ord[2] ) NULL else
              paste("diffy2Lag", 1:int.ord[2], sep = ""))
        
        list.names <- vector("list", NCOL(x) + 1 + sum(int.ord))
        names(list.names) <- colnames(x.matrix)
        for (ii in 1:( NCOL(x) + 1 + sum(int.ord)))
          list.names[[ii]] <- ii
        attr(x.matrix, "assign") <- list.names
        y <- cbind(temp1[, 1], temp2[, 1])
        colnames(y) <- paste("diffy", 1:NOS.f, sep = "")
        
        x <- x.matrix
        w <- w[c(1:NROW(y))]
        
      } else {
        stop("If you see this, then there is a bug. Please, report it.")
        temp1 <-  if ( !int.ord[1] ) NULL else 
          embed(y[, 1], 1 + max(int.ord))
        temp2 <-   if ( !int.ord[2] ) NULL else
          embed(y[, 2], 1 + max(int.ord))
        
        x.matrix <- cbind(x[-c(1:max(int.ord)), ],
                          #int.resids[-c(1:(max(int.ord)))],
                          if (!int.ord[1]) NULL else
                            temp1[, 2:(1 + int.ord[1])],
                          if (!int.ord[2]) NULL else
                            temp2[, 2:(1 + int.ord[2])])
        
        colnames(x.matrix) <- c(colnames(x),
                                paste("LagY1", 1:int.ord[1], sep = ""),
                                paste("LagY2", 1:int.ord[2], sep = ""))
        
        list.names <- vector("list", NCOL(x) + sum(int.ord))
        names(list.names) <- colnames(x.matrix)
        for (ii in 1:( NCOL(x) + sum(int.ord)))
          list.names[[ii]] <- ii
        attr(x.matrix, "assign") <- list.names
        y <- cbind(y[-c(1:max(int.ord)), 1], y[-c(1:max(int.ord)), 2])
        colnames(y) <- paste("y", 1:NOS.f, sep = "")
        
        x <- x.matrix
        w <- w[c(1:NROW(y))] 
      }
      
    }), list( .my.ord = my.ord , .lag.res = lag.res ,
              .V.class = V.class , .r.pat = r.pat ,
              .ordtsDyn = ordtsDyn ))),
      
      
      
      
      
      infos = eval(substitute(function(...){
        list(M1 = NULL, # Can't set the expression for M1, Q1. No variables
             Q1 = NULL, # are allowed to be read in this slot,
             expected = TRUE,
             multipleResponses = FALSE,
             lmean = .lmean ,
             lvar  = .lvar ,
             loffd = .loffd ,
             #parameters.names = parameters.names,
             zero  = .zero )
        
      }, list( .zero = zero , .lmean = lmean , .lvar = lvar ,
               .loffd = loffd , .trinorm = trinorm ))),
      
      
      
      
      
   initialize = eval(substitute(expression({
        
        Q1 <- ncol(y)
        M <- M1 <- Q1 + Q1 * (Q1 + 1) / 2
        
        check <- w.y.check(w = w, y =y, 
                           Is.positive.y = FALSE, 
                           ncol.w.max = 1, 
                           ncol.y.max = Q1, 
                           ncol.y.min = Q1,
                           out.wy = TRUE, 
                           colsyperw = Q1, 
                           maximize = TRUE)
        
        w <- check$w
        y <- check$y
        NOS.b <- extra$n.means  <- extra$NOS.b <- Q1
        nn <- nrow(y)
        extra$M <- M
        
        der.list <- vector("list", extra$M)
        der.MUs <- vector("list", extra$M)
        der.Vars <- vector("list", extra$M)
        
        dMUs <- array(0, dim = c(extra$NOS.b, 1, extra$NOS.b))
        for (ii in 1:extra$NOS.b) 
          der.list[[ii]] <- dMUs[, , ii] <- 
                     diag(extra$NOS.b)[, ii, drop = FALSE]
        
        dVs <- array(0, dim = c(extra$NOS.b, extra$NOS.b, extra$NOS.b))
        
        for(ii in 1:extra$NOS.b) {
          dVs[, , ii][ii, ii] <- 1
          der.list[[ii + extra$NOS.b]] <- dVs[, ,  ii]
        }
        
        dTs <- array(0, dim = c(extra$NOS.b, extra$NOS.b,
                                M1 - 2 * extra$NOS.b))
        
        hlcount <- 1
        for (ii in 1:(extra$NOS.b - 1)) {
          for (jj in (ii + 1):extra$NOS.b) {
            dTs[, , hlcount][ii, jj] <- 1
            dTs[, , hlcount][jj, ii] <- 1
            hlcount <- hlcount + 1
          }
        }; rm(hlcount)
        
        for (ii in 1:dim(dTs)[3])
          der.list[[ii + 2 * extra$NOS.b]] <- dTs[, , ii]
        
        for (ii in 1:extra$NOS.b)
          der.MUs[[ii]] <- cbind(dMUs[, , ii])
        for (ii in (1 + extra$NOS.b):M)
          der.MUs[[ii]] <- matrix(0, extra$NOS.b, 1)
        
        for (ii in 1:extra$NOS.b)
          der.Vars[[ii]] <- matrix(0, extra$NOS.b, extra$NOS.b)
        
        for (ii in (1 + extra$NOS.b):(2 * extra$NOS.b))
          der.Vars[[ii]] <- der.list[[ii]]
          
        for (ii in (2 * extra$NOS.b + 1):M)
          der.Vars[[ii]] <- der.list[[ii]]
        
        extra$derVs  <- dVs
        extra$derTs  <- dTs
        extra$derMUs <- dMUs
        rm(dMUs)
        extra$der.MUs <- der.MUs
        extra$der.Vars <- der.Vars
        
        parameters.names <-  c( if ( .V.class %in% c("ECM") )
                   paste("Diff", 1:NOS.b, sep = "") else 
                     paste("Mean", 1:NOS.b, sep = ""),
                              paste("var", 1:NOS.b, sep = ""))
        pre.name <- character(NOS.b * (NOS.b + 1) / 2 - NOS.b)
        
        for (ii in 1:(NOS.b - 1)) 
          pre.name[ii] <- paste(paste("cov", ii, sep = ""), 
                                (ii + 1):NOS.b, sep = "")
        parameters.names <- c(parameters.names, pre.name)
        
        predictors.names <- 
          c(namesof(c( if ( .V.class %in% c("ECM") )
                       paste("Diff", 1:NOS.b, sep = "") else 
                            paste("Mean", 1:NOS.b, sep = "")),
                    link = .lmean , earg = .emean , tag = FALSE),
            namesof(paste("var", 1:NOS.b, sep = ""),
                    link = .lvar , earg = .evar , tag = FALSE))
        
        pre.name <- character(0)
        for (ii in 1:(NOS.b - 1)) 
          pre.name <- c(pre.name,
                        paste(paste("cov", ii, sep = ""), 
                              (ii + 1):NOS.b, sep = ""))
        
        predictors.names <- c(predictors.names,
                              namesof(pre.name, .loffd , .eoffd ))
        extra$colnames.y <- colnames(y)
        
        ini.mean <- matrix(apply(y, 2, function(x) weighted.mean(x, w = w)),
                           nn, NOS.b , byrow = TRUE)
        ini.var <- cov(y)
        
        real.var <- ini.var[row(ini.var) > col(ini.var)]
        real.var <- matrix(real.var, nn, M - 2 * Q1, byrow = TRUE)
        ini.var <- matrix(diag(ini.var), nn, NOS.b, byrow = TRUE)
        
        etastart <- cbind(theta2eta(ini.mean, .lmean , earg = .emean ),
                      theta2eta(ini.var, .lvar , earg = .evar ),
                      theta2eta(real.var, .loffd , earg = .eoffd ))
        etastart
        
      }), list( .trinorm = trinorm , 
                .lmean = lmean , .lvar = lvar , .loffd = loffd ,
                .emean = emean , .evar = evar , .eoffd = eoffd ,
                .V.class = V.class ))),
      
      
      
      
      linkinv = eval(substitute(function(eta, extra = NULL) {
        eta2theta(eta[, 1:extra$NOS.b, drop = FALSE],
                              .lmean , earg = .emean )
        
      }, list( .trinorm = trinorm ,
                .lmean =lmean , .emean = emean ))),
      
      
      
      
      last = eval(substitute(expression({
        
        int.pat <- .r.pat
        misc$link <- c(rep( .lmean , extra$NOS.b),
                       rep( .lvar , extra$NOS.b),
                       rep( .loffd , M - 2*extra$NOS.b))
        names(misc$link) <- parameters.names
        
        misc$earg <- vector("list", length = M)
        for (ii in 1:extra$NOS.b) {
          misc$earg[[ii]] <- .emean
          misc$earg[[extra$NOS.b + ii]] <- .evar
        }
          
        for (ii in (2 * extra$NOS.b):(M - extra$NOS.b))
          misc$earg[[ii]] <- .eoffd
        names(misc$earg) <- parameters.names

        misc$expected <- TRUE
        misc$multipleResponses <- FALSE
        misc$coint.vector <- extra$coint.vec
        misc$coint.res <- extra$coint.res
        
        my.sample <- matrix(nrow(y), 1, 1)
        rownames(my.sample) <- "Final sample size:"
        colnames(my.sample) <- c("")
        c.vec <- c(1, -extra$coint.vec)
        names(c.vec) <- c("betaY2", names(extra$coint.vec))
        
        int.type  <- ifelse((int.pat == "intercept") || 
                             (int.pat == "intercept"), "level", "trend") 
        test.coint <- KPSS.test(extra$coint.res, type.H0 = int.type,
                                show.output = FALSE)$pvalue
        if (test.coint < 5e-2) {
          cat("\nWarning:",
              "The residuals appear to possess a unit root.")
        } else {
          cat("\nThe residuals appear to be free of unit roots.")
        }
          
        if ( .V.class %in% c("ECM")) {
          cat("\nEstimated co-integrating vector:\n")
          print(c.vec)
          print(my.sample)
        }
        rm(c.vec, my.sample)
        
        if ( .ordtsDyn )
          cat("\n Argument 'ordtsDyn' is on. Results should be similar\n",
              "to those from VECM(), from package 'tsDyn'.")
          
      }), list( .trinorm = trinorm , 
                .lmean = lmean , .lvar = lvar , .loffd = loffd ,
                .emean = emean , .evar = evar , .eoffd = eoffd ,
                .V.class = V.class , .r.pat = r.pat ,
                .ordtsDyn = ordtsDyn ))),
      
      
      
      loglikelihood = eval(substitute(
        function(mu, y , w,
                 residuals = FALSE, eta,
                 extra = NULL,
                 summation = TRUE) {
          
          nn <- nrow(y)
          
          vec.mean <- eta2theta(eta[, 1:extra$NOS.b],
                                .lmean , earg = .emean )
          var.only <- eta2theta(eta[, (1 + extra$NOS.b):(2 * extra$NOS.b)],
                                .lvar  , earg = .evar )
          offds <- eta2theta(eta[, -c(1:(2 * extra$NOS))],
                             .loffd , earg = .eoffd )
          
          if (residuals) {
            stop("loglikelihood residuals not implemented yet")
          } else {
            ll.elts <- c(w) * dmultinorm(vec.x = y, vec.mean = vec.mean,
                                mat.cov = cbind(var.only, offds),
                                log = TRUE)
            if (summation) {
              sum(ll.elts)
            } else {
              ll.elts
            }
          }
        }, list( .trinorm = trinorm , 
                 .lmean = lmean , .lvar = lvar , .loffd = loffd ,
                 .emean = emean , .evar = evar , .eoffd = eoffd  ))),
      
      
      
      vfamily = c("ECM.EngleGran"),
      
      
      
      deriv = eval(substitute(expression({
        nn <- nrow(eta)
        counter <- cbind(1:nn)
        M1 <- extra$n.means + ( extra$n.means *
                                  (extra$n.means + 1) )/2
        
        dVs  <- extra$derVs 
        dTs <- extra$derTs
        dMUs <- extra$derMUs
        der.MUs  <- extra$der.MUs
        der.Vars <- extra$der.Vars
        
        vec.mean <- eta2theta(eta[, 1:extra$NOS.b, drop = FALSE],
                              .lmean , earg = .emean )
        var.only <- eta2theta(eta[, (1 + extra$NOS.b):(2 * extra$NOS.b),
                                  drop = FALSE], .lvar  , earg = .evar )
        offds <- eta2theta(eta[, -c(1:(2 * extra$NOS)), drop = FALSE],
                           .loffd , earg = .eoffd )
        
        bigMat.der <- matrix(NA_real_, nn, M)
        nos.int <- extra$NOS.b
        
        mat11 <- cbind(var.only, offds)
        temp1 <-  apply(mat11, 1 , function(x) {
          int.mat <- matrix(0, nos.int, nos.int)
          diag(int.mat) <- x[1:nos.int]
          int.mat[row(int.mat) > col(int.mat)] <- x[-c(1:nos.int)]
          int.mat <- t(int.mat)
          int.mat[row(int.mat) > col(int.mat)] <- x[-c(1:nos.int)]
          int.mat
        }); rm(nos.int)
        
        attr(temp1, "dim") <- c(extra$NOS.b, extra$NOS.b, nn)
        
        temp1 <- apply(temp1, 3, function(x) solve(x))
        attr(temp1, "dim") <- c(extra$NOS.b, extra$NOS.b, nn)
        
        # From @initialize
        dVs <- extra$derVs
        dTs <- extra$derTs
        
        for (ii in 1:(extra$NOS.b)) {
          temp2 <- apply(temp1, 3, function(x) {
            (-0.5) *  sum(diag(( x %*% dVs[, , ii] )))
          })
          
          temp3 <- apply(temp1, 3, function(x) {
            x %*% dVs[, , ii] %*% x 
          })
          
          dim(temp3) <- c(extra$NOS.b, extra$NOS.b, nn)
          
          temp3 <- apply(counter, 1, function(x) {
            (0.5) * (y - mu)[x, ] %*% temp3[, , x] %*% cbind((y - mu)[x, ])
          })
          
          bigMat.der[, extra$NOS.b + ii] <- temp2 + temp3
          
        }
        
        for (ii in 1:(M1 - 2 * extra$NOS.b)) {
          
          temp2 <- apply(temp1, 3, function(x) {
            (-0.5) *  sum(diag(( x %*% dTs[, , ii] )))
          })
          
          temp3 <- apply(temp1, 3, function(x) {
            x %*% dTs[, , ii] %*% x 
          })
          
          dim(temp3) <- c(extra$NOS.b, extra$NOS.b, nn)
          
          temp3 <- apply(counter, 1, function(x) {
            (0.5) * (y - mu)[x, ] %*% temp3[, , x] %*% cbind((y - mu)[x, ])
          })
          
          bigMat.der[, 2 * extra$NOS.b + ii] <- temp2 + temp3
          
          
        }
        
        temp3 <- apply(counter, 1, function(x) {
          temp1[, , x] %*% cbind((y - mu))[x, ]
        })
        
        bigMat.der[, 1:(extra$NOS.b)] <- t(temp3)
       
        ### Der theta / d.eta
        dmeans.deta <- dtheta.deta(vec.mean, .lmean , .emean )
        dvars.deta  <- dtheta.deta(var.only, .lvar, .evar )
        doffds.deta <- dtheta.deta(offds, .loffd , .eoffd )
        
        c(w) * bigMat.der * cbind(dmeans.deta,
                                  dvars.deta,
                                  doffds.deta)
        
      }), list( .trinorm = trinorm , 
                .lmean = lmean , .lvar = lvar , .loffd = loffd ,
                .emean = emean , .evar = evar , .eoffd = eoffd  ) )),
      
      
      
      
      weight = eval(substitute(expression({
        
        wz <- matrix(0.0, nn, dimm(extra$M))
        
        
        for (ii in 1:(extra$NOS.b)) 
          wz[, ii] <- temp1[ii, ii, ] * dmeans.deta[, ii]^2
        
        for (ii in 1:extra$NOS.b)
          wz[, (ii + extra$NOS.b)] <- (0.5) * temp1[ii, ii, ]^2 *
                                                dvars.deta[, ii]^2
          
        for (ii in 1:(M1 - 2 * extra$NOS.b))
             wz[, (2 * extra$NOS.b + ii)] <- apply(counter, 1, function(x) {
               (0.5) * sum(diag(temp1[, , x] %*% dTs[, , ii] %*% 
                                temp1[, , x] %*% dTs[, , ii])) 
             }) * doffds.deta[, ii]^2
        
        pos.thetas <- combVGAMextra(x = 1:extra$M)[-c(1:extra$M), ]
        Der.dETAS <- cbind(dmeans.deta, dvars.deta, doffds.deta)
        
        for (ii in 1:nrow(pos.thetas)) {
          x.row <- pos.thetas[ii, ]
          wz[, extra$M + ii] <- apply(counter, 1, function(x){
            ((0.5) * sum(diag(temp1[, , 1] %*% der.Vars[[x.row[1]]] %*% 
                       temp1[, , 1] %*% der.Vars[[x.row[2]]])) +
                            t(der.MUs[[x.row[1]]]) %*%
                         temp1[, , 1] %*% der.MUs[[x.row[2]]]) 
          })
          wz[, extra$M + ii] <- wz[, extra$M + ii] * 
                        Der.dETAS[, x.row[1]] * Der.dETAS[, x.row[2]] 
        }
       
        c(w) * wz
        
      }), list( .trinorm = trinorm , 
                .lmean = lmean , .lvar = lvar , .loffd = loffd ,
                .emean = emean , .evar = evar , .eoffd = eoffd  ))) )

  alo   
}