BTSM: Bayesian Time Series Models

#' @name .irf.sign.zero
#' @noRd
#' @importFrom MASS Null
.irf.sign.zero <- function(xdat,plag,n.ahead,Amat,Smat,Emat,shockinfo,MaxTries,...){
  bigT     <- nrow(xdat)
  bigK     <- ncol(xdat)
  varNames <- colnames(xdat)

  P0G <- try(t(chol(Smat[,,drop=FALSE])),silent=TRUE)
  if(is(P0G,"try-error")) suppressWarnings(P0G <- t(chol(Smat,pivot=TRUE)))
  colnames(P0G) <- rownames(P0G) <- varNames

  # create dynamic multiplier
  PHIx <- array(0,c(bigK,bigK,plag+n.ahead+1)); dimnames(PHIx)[[1]] <- dimnames(PHIx)[[2]] <- varNames
  PHIx[,,plag+1]  <-  diag(bigK)
  temp <- .gen_compMat(Amat, bigK, plag)
  Cm   <- temp$Cm
  Jm   <- temp$Jm
  Cmat <- diag(bigK*plag)
  for (ihor in (plag+2):(plag+n.ahead+1)){
    Cmat  <- Cmat%*%Cm
    PHIx[,,ihor]  <- t(Jm)%*%Cmat%*%Jm
  }
  PHI  <-  PHIx[,,(plag+1):(plag+n.ahead+1)]
  #-----------------------------------------------------------------------------
  sign.horizon   <- unique(shockinfo$horizon)
  sign.horizon   <- sort(sign.horizon, decreasing=FALSE)
  sign.shockvars <- unique(shockinfo$shock)
  H.restr        <- length(sign.horizon)
  N.restr        <- bigK*H.restr
  S.cube         <- array(NA, c(N.restr, N.restr, bigK)) # sign restrictions
  Z.cube         <- array(NA, c(N.restr, N.restr, bigK)) # zero restrictions
  dimnames(S.cube)[[1]] <- dimnames(Z.cube)[[1]] <-
    dimnames(S.cube)[[2]] <- dimnames(Z.cube)[[2]] <- paste(rep(varNames,H.restr),".",
                                                            rep(sign.horizon,each=bigK),sep="")
  dimnames(S.cube)[[3]] <- dimnames(Z.cube)[[3]] <- varNames
  for(vv in 1:length(varNames)){
    S.temp <- matrix(0, N.restr, N.restr)
    Z.temp <- matrix(0, N.restr, N.restr)
    if(varNames[vv]%in%sign.shockvars){
      idx <- which(shockinfo$shock==varNames[vv])
      sign.restr <- shockinfo$restrictions[idx]
      sign.signs <- shockinfo$sign[idx]
      sign.horiz <- shockinfo$horizon[idx]

      s.point <- which(sign.signs=="<"|sign.signs==">")
      z.point <- seq(1,length(idx))[-s.point]

      if(length(s.point)>0){
        for(ss in 1:length(s.point)){
          grp <- which(sign.horiz[s.point[ss]] == sign.horizon)
          col <- seq(which(sign.restr[s.point[ss]]==varNames),bigK*grp,by=bigK)
          for(ii in 1:length(col)){
            S.temp[col[ii],col[ii]] <- ifelse(sign.signs[s.point[ss]]=="<",-1,1)
          }
        }
      }
      if(length(z.point)>0){
        for(zz in 1:length(z.point)){
          if(sign.signs[z.point[zz]]=="0"){
            grp <- which(sign.horiz[z.point[zz]] == sign.horizon)
            row <- (grp-1)*bigK+which(sign.restr[z.point[zz]]==varNames)
            Z.temp[row,row] <- 1
          }else{ # take row from above
            grp <- which(sign.horiz[z.point[zz]] == sign.horizon)
            col <- (grp-1)*bigK+which(sign.restr[z.point[zz]]==varNames)
            Z.temp[row,col] <- as.numeric(sign.signs[z.point[zz]])
          }
        }
      }
    }
    S.cube[,,vv] <- S.temp
    Z.cube[,,vv] <- Z.temp
  }

  no.zero.restr <- ifelse(base::sum(abs(Z.cube))>0,FALSE,TRUE)
  shock.order   <- rep(NA, bigK)
  search.Znum   <- apply(Z.cube, 3, function(x) base::sum(abs(x)))
  search.Snum   <- apply(S.cube, 3, function(x) base::sum(abs(x)))
  for(mm in 1:bigK){
    if(!no.zero.restr){
      max.Z <- which(search.Znum==max(search.Znum))
      if(length(max.Z)==1){
        shock.order[mm] <- max.Z
      }else{
        shock.order[mm] <- sample(max.Z,1)
      }
    } else {
      shock.order[mm] <- mm
    }
    search.Znum[shock.order[mm]] <- -1
    search.Snum[shock.order[mm]] <- -1
  }
  shock.order <- varNames[shock.order]

  irf.restr         <- matrix(NA, N.restr, bigK)
  for(hh in 1:H.restr){
    # ARRW: Definition 1
    if(sign.horizon[hh]!=Inf) irf.hh<-PHI[,,sign.horizon[hh]]%*%P0G
    # ARRW: Definition 2
    #if(sign.horizon[hh]==Inf) irf.hh <- solve(A0-A0%*%Cm[1:M,]%*%do.call("rbind",rep(list(diag(M)),p)))
    irf.restr[((hh-1)*bigK+1):(bigK*hh),1:bigK] <- irf.hh
  }
  colnames(irf.restr) <- varNames
  rownames(irf.restr) <- paste(rep(varNames,H.restr),".",
                               rep(sign.horizon,each=bigK),sep="")

  Z.cube <- Z.cube[,,shock.order]

  # draw rotation matrix here
  icounter <- 0
  condall <- 0
  max.counter <- MaxTries
  impresp<-Q_bar<-NA
  while(condall == 0 && icounter < max.counter){
    signCheck <- matrix(NA, bigK, 1)
    randMat <- matrix(rnorm(bigK*bigK,0,1),bigK,bigK)
    Q <- matrix(0, bigK, bigK)
    if(no.zero.restr){
      Q <- qr(randMat)
      Q <- qr.Q(Q)
    }else{
      for(mm in 1:bigK){
        Z.temp <- Z.cube[,,mm]
        Z.temp <- Z.temp[rowSums(abs(Z.temp))!=0,,drop=F]
        if(nrow(Z.temp)==0){
          Z.temp <- matrix(0, 1, N.restr)
        }
        if(all(Z.temp==0) && mm>1){
          R <- c()
        }else{
          R <- Z.temp%*%irf.restr
        }
        if(mm > 1){R <- rbind(R, t(Q[,(1:(mm-1)), drop=FALSE]))}

        NU  <- Null(t(R))
        x_j <- randMat[,mm,drop =FALSE]

        q_j <- NU%*%(t(NU)%*%x_j/sqrt(as.numeric(crossprod(t(NU)%*%x_j))))
        Q[,mm] <- q_j
      }
    }
    colnames(Q) <- shock.order; rownames(Q) <- varNames
    Q <- Q[,varNames]
    Q_bar <- Q%*%diag(((diag(Q)>0)-(diag(Q)<0)))

    irf.check <- irf.restr%*%Q_bar
    colnames(irf.check) <- varNames
    rownames(irf.check) <- paste(rep(varNames,H.restr),".",rep(sign.horizon,each=bigK),sep="")

    for(ss in 1:bigK){
      STemp <- as.matrix(diag(S.cube[,,ss]))
      IrfCheckTemp <- irf.check[,ss,drop = FALSE]
      signCheck[ss,] <- t(sign(IrfCheckTemp))%*%STemp==sum(abs(STemp))
    }
    condall <- prod(signCheck)
    icounter <- icounter + 1
  }

  shock <- P0G%*%Q_bar
  # computing impulse responses
  irfa  <- array(0,c(n.ahead,bigK,bigK)); dimnames(irfa)[[2]] <- dimnames(irfa)[[3]] <- varNames
  for (ihor in 1:n.ahead){
    irfa[ihor,,] <- PHI[,,ihor]%*%shock
  }
  # computing structural errors
  eps <- Emat%*%shock

  if(icounter==MaxTries){
    irfa <- eps <- NA
  }
  # end rotation matrix loop ----------------------------------------------------------------------------
  return(list(impl=irfa,rot=Q_bar,eps=eps,icounter=icounter))
}

#' @name .irf.chol
#' @noRd
.irf.chol <- function(xdat,plag,n.ahead,Amat,Smat,Emat,type,...){
  bigT      <- nrow(xdat)
  bigK      <- ncol(xdat)
  varNames  <- colnames(xdat)
  cons      <- FALSE
  if(nrow(Amat)%%plag!=0) cons <- TRUE

  if(type=="long-run"){
    Asum <- matrix(0,bigK,bigK)
    for(pp in 1:plag){
      Asum <- Asum + Amat[((pp-1)*bigK+1):(pp*bigK),]
    }
    A1inv <- solve(diag(bigK)-t(Asum))
    C1 <- try(t(chol(A1inv%*%Smat%*%t(A1inv))),silent=TRUE)
    if(is(C1,"try-error")) suppressWarnings(C1 <- t(chol(A1inv%*%Smat%*%t(A1inv),pivot=TRUE)))
    shock <- solve(A1inv)%*%C1
  }else if(type=="short-run"){
    shock <- try(t(chol(Smat)),silent=TRUE)
    if(is(shock,"try-error")) suppressWarnings(shock <- t(chol(Smat,pivot=TRUE)))
  }
  colnames(shock) <- rownames(shock) <- varNames

  diagonal <- diag(diag(shock))
  shock    <- solve(diagonal)%*%shock

  # create dynamic multiplier
  PHIx <- array(0,c(bigK,bigK,plag+n.ahead+1)); dimnames(PHIx)[[1]] <- dimnames(PHIx)[[2]] <- varNames
  PHIx[,,plag+1]  <-  diag(bigK)
  temp <- .gen_compMat(Amat, bigK, plag)
  Cm   <- temp$Cm
  Jm   <- temp$Jm
  Cmat <- diag(bigK*plag)
  for (ihor in (plag+2):(plag+n.ahead+1)){
    Cmat  <- Cmat%*%Cm
    PHIx[,,ihor]  <- t(Jm)%*%Cmat%*%Jm
  }
  PHI  <-  PHIx[,,(plag+1):(plag+n.ahead+1)]

  # computing impulse response function
  irfa  <- array(0,c(n.ahead,bigK,bigK)); dimnames(irfa)[[2]] <- dimnames(irfa)[[3]] <- varNames
  for (ihor in 1:n.ahead){
    irfa[ihor,,] <- PHI[,,ihor]%*%shock
  }
  # computing structural errors
  eps <- Emat%*%shock

  out <- list(impl=irfa,rot=shock,eps=eps)
  return(out)
}

#' @name .irf.proxy
#' @importFrom stats lm
#' @noRd
.irf.proxy <- function(xdat,plag,n.ahead,Amat,Smat,Emat,proxy,shockinfo,...){
  bigT      <- nrow(xdat)
  bigK      <- ncol(xdat)
  varNames  <- colnames(xdat)
  shockvars <- shockinfo$shock
  instrvars <- shockinfo$instr
  scale     <- shockinfo$scale

  # create dynamic multiplier
  PHIx <- array(0,c(bigK,bigK,plag+n.ahead+1)); dimnames(PHIx)[[1]] <- dimnames(PHIx)[[2]] <- varNames
  PHIx[,,plag+1]  <-  diag(bigK)
  temp <- .gen_compMat(Amat, bigK, plag)
  Cm   <- temp$Cm
  Jm   <- temp$Jm
  Cmat <- diag(bigK*plag)
  for (ihor in (plag+2):(plag+n.ahead+1)){
    Cmat  <- Cmat%*%Cm
    PHIx[,,ihor]  <- t(Jm)%*%Cmat%*%Jm
  }
  PHI  <-  PHIx[,,(plag+1):(plag+n.ahead+1)]

  randMat <- matrix(rnorm(bigK*bigK,0,1),bigK,bigK)
  Q <- qr(randMat)
  Q <- qr.Q(Q)

  # identification of shocks
  option <- 1

  F_stats <- matrix(NA, 4, bigK, dimnames=list(c("F_test","rob-F_test","F_test_lag","rob-F_test_lag"), varNames))

  # option 1
  if(option == 1){
    for(ss in 1:length(shockvars)){
      sd <- ifelse(scale[ss]=="sd", TRUE, FALSE)
      proxyVar <- proxy[,instrvars[ss],drop=FALSE]
      iP  <- which(varNames==shockvars[ss])
      niP <- (1:bigK)[-iP]
      res <- cbind(Emat[,iP],Emat[,niP])
      # adjust length in case of NAs
      if(any(is.na(proxyVar))){
        idx <- which(is.na(proxyVar))
        proxyVar <- proxyVar[-idx,,drop=FALSE]
        res <- res[-idx,,drop=FALSE]
      }
      Xdum   <- kronecker(diag(bigK), cbind(1,proxyVar))
      betaIV <- solve(crossprod(Xdum))%*%crossprod(Xdum,matrix(res,ncol=1))
      betaIV <- t(matrix(betaIV,nrow=nrow(betaIV)/bigK,bigK))

      beta_11 <- betaIV[1,2]
      beta_21 <- betaIV[-1,2,drop=FALSE]
      B21B11  <- beta_21/beta_11

      # fitted.err <- lm(res[,iP] ~ proxyVar)$fitted
      # b21ib11    <- t(lm(res[,niP] ~ fitted.err)$coef)

      SigmaU <- Smat[c(iP,niP),c(iP,niP)]
      Zeta   <- B21B11%*%SigmaU[1,1]%*%t(B21B11) - SigmaU[2:bigK,1]%*%t(B21B11) + B21B11%*%t(SigmaU[1,2:bigK]) + SigmaU[2:bigK,2:bigK]
      B12B12 <- t(SigmaU[2:bigK,1]-B21B11%*%SigmaU[1,1])%*%solve(Zeta)%*%(SigmaU[2:bigK,1]-B21B11%*%SigmaU[1,1])
      B11B11 <- SigmaU[1,1] - B12B12
      if(B11B11<0)
        return(list(impl=NA,rot=NA,eps=NA))
      B11    <- sqrt(B11B11)
      shock  <- c(B11, B21B11*c(B11))
      if(!sd)
        shock  <- shock/shock[1]

      Q[c(iP,niP),iP] <- shock

      # # F stat (regression on instruments of relevant innovations)
      # tempX  <- cbind(1, proxyVar)
      # tempU  <- res[,iP] - tempX%*%t(betaIV[iP,,drop=FALSE])
      # tempY  <- tempX%*%t(betaIV[iP,,drop=FALSE]) - matrix(mean(res[,iP]),nrow(res),1)
      # k <- length(betaIV[iP,])-1
      # F_stat <- (t(tempY)%*%tempY/k)%*%solve(t(tempU)%*%tempU/(nrow(tempU)-k-1))

      # F-test
      reg1 <- lm(res[,iP]~proxyVar)
      reg0 <- lm(res[,iP]~1)
      F_stats[1,iP] = anova(reg0, reg1)$F[2]
      F_stats[2,iP] = waldtest(reg0, reg1, vcov=vcovHC(reg1, type="HC3"))$F[2]

      # lags as controls
      proxyVarlag <- cbind(proxyVar, .mlag(proxyVar, plag))
      proxyVarlag <- proxyVarlag[(plag+1):nrow(proxyVarlag),,drop=FALSE]
      resuse      <- res[(plag+1):nrow(res),iP,drop=FALSE]

      reg1  <- lm(resuse~proxyVarlag)
      reg0 <- lm(resuse~1)
      F_stats[3,iP] = anova(reg0, reg1)$F[2]
      F_stats[4,iP] = waldtest(reg0, reg1, vcov=vcovHC(reg1, type="HC3"))$F[2]
    }
  }
  ## option 2
  if(option == 2){
    m   <- ncol(proxy)
    iP  <- which(varNames%in%shockvars)
    niP <- (1:bigK)[-iP]
    res <- cbind(Emat[,iP],Emat[,niP])
    # adjust length in case of NAs
    if(any(is.na(proxy))){
      idx <- which(is.na(proxy))
      proxy <- proxy[-idx,,drop=FALSE]
      res <- res[-idx,,drop=FALSE]
    }
    Xdum   <- kronecker(diag(bigK), cbind(1,proxy))
    betaIV <- solve(crossprod(Xdum))%*%crossprod(Xdum,matrix(res,ncol=1))
    betaIV <- t(matrix(betaIV,nrow=nrow(betaIV)/bigK,bigK))

    beta_11 <- betaIV[1:m,2:(m+1)]
    beta_21 <- betaIV[(m+1):nrow(betaIV),2:(m+1),drop=FALSE]
    B21B11  <- beta_21%*%solve(beta_11)

    # fitted.err <- lm(res[,iP] ~ proxyVar)$fitted
    # b21ib11    <- t(lm(res[,niP] ~ fitted.err)$coef)

    SigmaU <- Smat[c(iP,niP),c(iP,niP)]
    Zeta   <- B21B11%*%SigmaU[1:m,1:m]%*%t(B21B11) - SigmaU[(m+1):bigK,1:m]%*%t(B21B11) + B21B11%*%t(SigmaU[(m+1):bigK,1:m]) + SigmaU[(m+1):bigK,(m+1):bigK]
    B12B12 <- t(SigmaU[(m+1):bigK,1:m]-B21B11%*%SigmaU[1:m,1:m])%*%solve(Zeta)%*%(SigmaU[(m+1):bigK,1:m]-B21B11%*%SigmaU[1:m,1:m])
    B11B11 <- SigmaU[1:m,1:m] - B12B12

    if(m == 1){
      B11    <- sqrt(B11B11)
      shock  <- c(B11, B21B11*c(B11))
      shock  <- shock/shock[1]

      Q[c(iP,niP),iP] <- shock
    }else{
      B22B22 <- SigmaU[(m+1):bigK,(m+1):bigK] + B21B11%*%(B12B12 - SigmaU[1:m,1:m])%*%t(B21B11)
      B12B22 <- (B12B12%*%t(B21B11) + t(SigmaU[(m+1):bigK,1:m] - B21B11%*%SigmaU[1:m,1:m]))%*%solve(B22B22)
      B11S1  <- diag(m) - B12B22%*%B21B11
      B21S1  <- B21B11%*%solve(B11S1)
      S1S1   <- B11S1%*%B11B11%*%t(B11S1)
      S1     <- try(t(chol(S1S1)),silent=TRUE)
      if(is(S1,"try-error"))
        return(list(impl=NA,rot=NA,eps=NA))
      shock  <- rbind(solve(B11S1), B21S1)%*%S1

      Q[c(iP,niP),iP] <- shock
    }

    # # F stat (regression on instruments of relevant innovations)
    # tempX  <- cbind(1, proxy)
    # tempU  <- res[,1:m] - tempX%*%t(betaIV[1:m,])
    # tempY  <- tempX%*%t(betaIV[1:m,]) - matrix(mean(res[,1:m]),nrow(res),m); k <- length(betaIV[1:m,])-1
    # F_stat <- (t(tempY)%*%tempY/k)%*%solve(t(tempU)%*%tempU/(nrow(tempU)-k-1))
  }

  # fitted.err <- lm(Eest[,1] ~ proxy)$fitted
  # b21ib11    <- t(lm(Eest[,-1] ~ fitted.err-1)$coef)
  # Sig11      <- matrix(Smat[1,1], 1, 1)
  # Sig21      <- matrix(Smat[2:bigK,1],bigK-1,1)
  # Sig12      <- matrix(Smat[1,2:bigK],1,bigK-1)
  # Sig22      <- matrix(Smat[2:bigK,2:bigK],bigK-1,bigK-1)
  # ZZp        <- b21ib11%*%Sig11%*%t(b21ib11) - Sig21%*%t(b21ib11)+b21ib11%*%t(Sig21)+Sig22
  # b12b12p    <- t(Sig21-b21ib11%*%Sig11)%*%solve(ZZp)%*%(Sig21-b21ib11%*%Sig11)
  # b11b11p    <- Sig11 - b12b12p
  # if(b11b11p<0){
  #   return(list(impl=NA,rot=NA,eps=NA))
  # }
  # b11        <- sqrt(b11b11p)
  # shock      <- c(b11, b21ib11*c(b11))
  # shock      <- shock/shock[1] # normalize to unit shock
  #
  # randMat <- matrix(rnorm(bigK*bigK,0,1),bigK,bigK)
  # Q <- qr(randMat)
  # Q <- qr.Q(Q)
  # Q[,1] <- shock

  # computing impulse response function
  irfa  <- array(0,c(n.ahead,bigK,bigK)); dimnames(irfa)[[2]] <- dimnames(irfa)[[3]] <- varNames
  for (ihor in 1:n.ahead){
    irfa[ihor,,] <- PHI[,,ihor]%*%Q
  }
  # computing structural errors
  eps <- Emat%*%Q

  return(list(impl=irfa,rot=Q,eps=eps,F_stats=F_stats))
}

#' @name .irf.girf
#' @noRd
.irf.girf <- function(xdat,plag,n.ahead,Amat,Smat,Emat,...){
  bigT      <- nrow(xdat)
  bigK      <- ncol(xdat)
  varNames  <- colnames(xdat)

  # create dynamic multiplier
  PHIx <- array(0,c(bigK,bigK,plag+n.ahead+1)); dimnames(PHIx)[[1]] <- dimnames(PHIx)[[2]] <- varNames
  PHIx[,,plag+1]  <-  diag(bigK)
  temp <- .gen_compMat(Amat, bigK, plag)
  Cm   <- temp$Cm
  Jm   <- temp$Jm
  Cmat <- diag(bigK*plag)
  for (ihor in (plag+2):(plag+n.ahead+1)){
    Cmat  <- Cmat%*%Cm
    PHIx[,,ihor]  <- t(Jm)%*%Cmat%*%Jm
  }
  PHI  <-  PHIx[,,(plag+1):(plag+n.ahead+1)]

  # computing impulse response function
  irfa  <- array(0,c(n.ahead,bigK,bigK)); dimnames(irfa)[[2]] <- dimnames(irfa)[[3]] <- varNames
  for (ihor in 1:n.ahead){
    irfa[ihor,,] <- PHI[,,ihor]%*%Smat
  }
  # computing structural errors
  eps <- Emat%*%Smat

  return(list(impl=irfa,rot=Smat,eps=eps))
}

#' @name .impulsdtrf
#' @noRd
.impulsdtrf <- function(B,smat,nstep)
  ### By:             As emerges from rfvar, neqn x nvar x lags array of rf VAR coefficients.
  ### smat:           nshock x nvar matrix of initial shock vectors.  To produce "orthogonalized
  ###                 impulse responses" it should have the property that crossprod(t(smat))=sigma,
  ###                 where sigma is the Var(u(t)) matrix and u(t) is the rf residual vector.  One
  ###                 way to get such a smat is to set smat=t(chol(sigma)).  To get the smat
  ###                 corresponding to a different ordering, use
  ###                 smat = t(chol(P %*% Sigma %*% t(P)) %*% P), where P is a permutation matrix.
  ###                 To get impulse responses for a structural VAR in the form A(L)y=eps, with
  ###                 Var(eps)=I, use B(L)=-A_0^(-1)A_+(L) (where A_+ is the coefficients on strictly
  ###                 positive powers of L in A), smat=A_0^(-1).
  ###                 In general, though, it is not required that smat be invertible.
### response:       nvar x nshocks x nstep array of impulse responses.
###
### Code written by Christopher Sims,mat based on 6/03 matlab code.  This version 3/27/04.
### Added dimension labeling, 8/02/04.
{

  neq <- dim(B)[1]
  nvar <- dim(B)[2]
  lags <- dim(B)[3]
  dimnB <- dimnames(B)
  if(dim(smat)[2] != dim(B)[2]) stop("B and smat conflict on # of variables")
  response <- array(0,dim=c(neq,nvar,nstep+lags-1));
  response[ , , lags] <- smat
  response <- aperm(response, c(1,3,2))
  irhs <- 1:(lags*nvar)
  ilhs <- lags * nvar + (1:nvar)
  response <- matrix(response, ncol=neq)
  B <- B[, , seq(from=lags, to=1, by=-1)]  #reverse time index to allow matrix mult instead of loop
  B <- matrix(B,nrow=nvar)
  for (it in 1:(nstep-1)) {
    response[ilhs, ] <- B %*% response[irhs, ]
    irhs <- irhs + nvar
    ilhs <- ilhs + nvar
  }
  dim(response) <- c(nvar, nstep + lags - 1, nvar)
  #drop the zero initial conditions; array in usual format
  if(lags>1){
    response<-response[,-(1:(lags-1)),]
  }
  response <- aperm(response, c(1, 3, 2))
  dimnames(response) <- list(dimnB[[1]], dimnames(smat)[[2]], NULL)
  ## dimnames(response)[2] <- dimnames(smat)[1]
  ## dimnames(response)[1] <- dimnames(B)[2]
  return(response)
}