R/internal.R

In jpburgard/mixtureVAR: mixtureVAR

isomLRinv2<-function (x)
{
  x <- -x
  y = matrix(0, nrow = nrow(x), ncol = ncol(x) + 1)
  D = ncol(x) + 1
  y[, 1] = -sqrt((D - 1)/D) * x[, 1]
  for (i in 2:ncol(y)) {
    for (j in 1:(i - 1)) {
      y[, i] = y[, i] + x[, j]/sqrt((D - j + 1) * (D -
                                                     j))
    }
  }
  if(ncol(y)>2){
  for (i in 2:(ncol(y) - 1)) {
    y[, i] = y[, i] - sqrt((D - i)/(D - i + 1)) * x[, i]
  }
  }
  yexp = exp(y)
  x.back = yexp/apply(yexp, 1, sum)
  if (is.data.frame(x))
    x.back <- data.frame(x.back)
  return(x.back)
}


VARlist <- function(Y, m, EXO = NULL) {
  if(is.null(EXO))  a <- vars::VAR(t(Y), m) else a <- VAR(t(Y), m, exogen = t(EXO))
  t(sapply(coef(a), function(aa) {
    aaa <- aa[, "Estimate"]
    const<-which(rownames(aa)=="const")
    c(aaa[length(aaa)], aaa[-const])
  }))
}

VARlist2 <- function(Y, m, EXO = NULL) {
  if(is.null(EXO))  a <- vars::VAR(t(Y), m) else a <- VAR(t(Y), m, exogen = t(EXO))
  t(sapply(coef(a), function(aa) {
    aaa <- aa[, "Estimate"]
    c(aaa[length(aaa)], aaa[-length(aaa)])
  }))
}

is.formula <- function(form)class(form)=="formula"

thetastepRest <- function(R, r, tau, X, Sigma, Y)
{
  tmp1 <- tau * t(apply(X, 2, tcrossprod))
  A <-
    solve(t(R) %*% (kronecker(Reduce(
      '+', lapply(1:nrow(tmp1), function(xx)
        matrix(tmp1[xx, , drop = FALSE], ncol = sqrt(ncol(
          tmp1
        ))))
    ), solve(Sigma))) %*% R)
  B1 <- solve(Sigma) %*% (Y - r %*% X)
  B2 <-
    Reduce('+', lapply(1:ncol(X), function(xx)
      tau[xx] * as.vector(B1[, xx, drop = FALSE] %*% t(X[, xx, drop = FALSE]))))
  matrix(A %*% (t(R) %*% B2), ncol = nrow(X), nrow = nrow(Y))
}

thetastep <- function(tau, X, Y, pd)
{
  if (pd > 0) {
    X      <- X[,-c(1:pd)]
    Y      <- Y[,-c(1:pd)]
  }
  tmp1   <-
    matrix(apply(tau * t(apply(X, 2, tcrossprod)), 2, sum),
           ncol = nrow(X),
           nrow = nrow(X))
  tmp2a  <- lapply(1:ncol(X), function(xx)
    X[, xx])
  tmp2b  <- lapply(1:ncol(Y), function(yy)
    Y[, yy])
  tmp2c  <-
    tau * t(mapply(function(X, Y)
      X %*% t(Y), tmp2a, tmp2b, SIMPLIFY = TRUE))
  tmp2   <- matrix(tmp2c, ncol = nrow(Y), nrow = nrow(X))
  return(t(solve(tmp1) %*% tmp2))
}

thetastep2 <- function(tau, X, Y, pd)
{
  if (pd > 0) {
    X      <- X[,-c(1:pd)]
    Y      <- Y[,-c(1:pd)]
  }
  tmp1   <- X %*% (tau * t(X))
  tmp2   <- X %*% (tau * t(Y))
  return(t(ginv(tmp1) %*% tmp2))
}

likl <- function(alph, Sig, e, pd) {
  alph * (det(Sig) ^ (-0.5)) * apply(e, 2, function(ee)
    exp(-0.5 * t(ee) %*% ginv(Sig) %*% ee))
}

res <- function(Y, X, theta, pd) {
  if (pd < 1)
    return(Y - theta %*% X)
  else
    return((Y - theta %*% X)[,-c(1:pd)])
}

Phi.mixvar <- function (x, nstep = 20, ...)
{
  if (!(class(x) == "mixvar")) {
    stop("\nPlease provide an object of class 'mixvar', generated by 'FMVAR()'.\n")
  }

  nstep <- abs(as.integer(nstep))
  K <- lapply(x$Theta,nrow)
  P <- x$P

  As <- mapply(function(theta,p,k){
    theta <- theta[,-1]
    ret <- array(0,dim=c(k,k,nstep+1))
    for(i in 1:p) ret[,,i] <- theta[,(1+(i-1)*k):(i*k)]
    ret
  },x$Theta,P,K,SIMPLIFY = FALSE)

  Phi <- mapply(function(k,as){
    phi <- array(0, dim = c(k, k, nstep + 1))
    phi[, , 1] <- diag(k)
    phi[, , 2] <- phi[, , 1] %*% as[, , 1]
    phi
  },K,As,SIMPLIFY=FALSE)

  Phi <- mapply(function(phi,as,k) {
    for (i in 3:(nstep + 1)) {
      tmp1 <- phi[, , 1] %*% as[, , i - 1]
      tmp2 <- matrix(0, nrow = k, ncol = k)
      idx <- (i - 2):1
      for (j in 1:(i - 2)) {
        tmp2 <- tmp2 + phi[, , j + 1] %*% as[, , idx[j]]
      }
      phi[, , i] <- tmp1 + tmp2
    }
    phi}
    ,Phi,As,K,SIMPLIFY = FALSE)
  return(Phi)
}

.irf.mixvar = function (x, impulse, response, y.names, n.ahead, ortho, cumulative)
{
  #irf  <- Phi.mixvar(x, nstep = n.ahead)# Alte Version, leider nicht orthogonalisiert!!!!!!!
  irf  <- Psi.mixvar(x, nstep = n.ahead)#
  irf  <-lapply(irf,function(tmpirf){
    dimnames(tmpirf) <- list(y.names,y.names,1:(n.ahead+1))
    tmpirf
  })
  idx  <- length(impulse)

  irs <- lapply(irf, function(tmpirf) {
    irs  <- list()
    for (i in 1:idx) {
      irs[[i]] <-
        matrix(t(tmpirf[response, impulse[i], 1:(n.ahead + 1)]), nrow = n.ahead + 1)
      colnames(irs[[i]]) <- response
      if (cumulative) {
        if (length(response) > 1)
          irs[[i]] <- apply(irs[[i]], 2, cumsum)
        if (length(response) == 1) {
          tmp <- matrix(cumsum(irs[[1]]))
          colnames(tmp) <- response
          irs[[1]] <- tmp
        }
      }
    }
    names(irs) <- impulse
    irs
  })
  return(irs)
}

Psi.mixvar <- function (x, nstep = 20, ...)
{
  if (!(class(x) == "mixvar")) {
    stop("\nPlease provide an object of class 'mixvar', generated by 'FMVAR()'.\n")
  }
  nstep   <- abs(as.integer(nstep))
  K       <- lapply(x$Theta,nrow)
  Phi     <- Phi(x, nstep = nstep)
  params  <- mapply(function(theta,k)ncol(theta)-1-k,x$Theta,K)
  DFS     <- mapply(function(param,e)ncol(e)-param,params,x$e,SIMPLIFY = FALSE)
  sigma.u <- mapply(function(e,DF)crossprod(t(e))/(DF),x$e,DFS,SIMPLIFY = FALSE)
  Sig     <- lapply(sigma.u,function(sig)t(chol(sig)))
  Psi     <- mapply(function(phi,sig){
    psi <- array(0,dim=dim(phi))
    for(i in 1:dim(phi)[3]) psi[,,i]<-phi[,,i]%*%sig
    psi
  },Phi,Sig,SIMPLIFY = FALSE)

  return(Psi)
}
# Alternative
# Psi.mixvar = function(x, nstep = 20, ...){
#   nstep = abs(as.integer(nstep))
#   Phi = Phi.mixvar(x, nstep = nstep)
#   params = 0
#   for(i in 1:length(Mx$Theta)){
#     params = params + ncol(x$Theta[[i]])
#   }
#   Sigma = x$Sigma
#   P = lapply(x$Sigma, function(Sigma) t(chol(Sigma)))
#   dim3 = nstep+1
#   lapply(Phi, function(Phi, dim3){
#     for(i in 1:dim3){
#       Psi[[, , i]] = Phi[[, , i]]%*%P
#     }
#     return(Psi)
#   })
# }

ysynthfun <- function(x,theta,tau,pd) {
  if(pd<1)
    return(t(tau * t(theta %*% x[,])))
  else
    return(t(tau * t(theta %*% x[,-c(1:pd)])))
}