R/SOBIdefl_ascov.R

In BSSasymp: Asymptotic Covariance Matrices of Some BSS Mixing and Unmixing Matrix Estimates

ASCOV_SOBIdefl <- function(psi, taus, Beta=NULL, A=NULL)
{
  p <- dim(psi)[2]
  q <- dim(psi)[1]
  K <- length(taus)  

  if(is.null(A)) A <- diag(p)

  if(q<(3*K)) psi <- rbind(psi,matrix(0,3*K-q,p))

  q <- dim(psi)[1]

  Psi <- matrix(0,q+1,p)
  for(i in 1:p){
    Psi[,i] <- c(1,psi[,i])
    Psi[,i] <- Psi[,i]/sqrt(sum(Psi[,i]^2))
  }
  
  PSI <- array(0,c(p,p,q+1)) 
  for(i in 1:p){
   for(j in 1:(q+1)){
     PSI[i,i,j] <- Psi[j,i]
   }  
  }

  F_tau <- array(0,c(p,p,q+1))
  for(i in 0:q){
   for(j in 1:(q+1-i)){
      F_tau[,,i+1] <- F_tau[,,i+1]+tcrossprod(diag(PSI[,,j]),diag(PSI[,,j+i]))
   }
  }

  if(is.null(Beta)) Beta <- 2*diag(p)+matrix(1,p,p)

  Lambda <- array(0,c(p,p,K))
  for(k in 1:K){
   for(j in 1:p){  
    Lambda[j,j,k] <- F_tau[j,j,taus[k]+1]
   }
  }
  
  Sum_lam <- rep(0,p)
  for(j in 1:p){
    for(k in 1:K){
      Sum_lam[j] <- Sum_lam[j]+Lambda[j,j,k]^2
    }
  }   

  P <- matrix(0,p,p)
  ord <- order(Sum_lam,decreasing=TRUE)
  for(j in 1:p){
    P[j,ord[j]] <- 1 
  }  
  
  for(k in 1:K){
    Lambda[,,k] <- tcrossprod(crossprod(t(P),Lambda[,,k]),P)
  }

  for(i in 1:(q+1)){
    F_tau[,,i] <- tcrossprod(crossprod(t(P),F_tau[,,i]),P)
  }


  ASCOV <- matrix(0,p^2,p^2)
  for(j in 1:p){
   for(i in 1:p){
     if(j<i){
      ASV <- .C("ascov_deflji", as.double(as.vector(F_tau)),as.double(as.vector(Lambda)), as.double(taus),as.integer(c(i-1,j-1,p,q,K)),as.double(as.vector(Beta)),res=double (2), PACKAGE="BSSasymp")$res
      
      ASCOV <- ASCOV+ASV[2]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,j]),tcrossprod(diag(p)[,j],diag(p)[,i]))+ASV[1]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,i]),tcrossprod(diag(p)[,j],diag(p)[,j]))

     }  
     
     if(i==j){
       ASCOV <- ASCOV+0.25*D_lm(F_tau,0,0,Beta)[i,i]*kronecker(tcrossprod(diag(p)[,i], diag(p)[,i]),tcrossprod(diag(p)[,i],diag(p)[,i])) 
     }
     
     if(i<j){
       ASV <- .C("ascov_deflij", as.double(as.vector(F_tau)),as.double(as.vector(Lambda)),as.double(taus),as.integer(c(i-1,j-1,p,q,K)),as.double(as.vector(Beta)), res=double (2),PACKAGE="BSSasymp")$res
       
       ASCOV <- ASCOV+ASV[2]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,j]),tcrossprod(diag(p)[,j],diag(p)[,i]))+ASV[1]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,i]), tcrossprod(diag(p)[,j],diag(p)[,j]))
     }
    
   }
  }  
  
  EMD <- sum(diag(ASCOV)-diag(ASCOV)*as.vector(diag(p)))
  W <- crossprod(t(P),solve(A))
  W <- crossprod(diag(sign(rowMeans(W))),W)
  A <- solve(W)
  COV_A <- crossprod(t(tcrossprod(kronecker(diag(p),A),ASCOV)),kronecker(diag(p),t(A)))
  COV_W <- crossprod(t(tcrossprod(kronecker(t(W),diag(p)),ASCOV)),kronecker(W,diag(p)))
  
  list(W=W, COV_W=COV_W, A=A, COV_A=COV_A, EMD=EMD)
}


ASCOV_SOBIdefl_estN <- function(X, taus, mixed=TRUE, M=100)
{
  p <- dim(X)[2]
  T <- dim(X)[1]
  
  if(length(taus)==1) taus <- 1:taus
  K <- length(taus)  
  
  if(mixed){
    W <- SOBI(X,taus,method="djd")$W
  }else W <- diag(p)
  
  X <- tcrossprod(sweep(X,2,colMeans(X)),W)  

  F_tau <- array(0,c(p,p,M+1))
  for(m in 0:M){
    F_tau[,,m+1] <- tcrossprod(t(X[1:(T-m),]),t(X[(m+1):T,]))/(T-m)
  }

  Beta <- 2*diag(p)+matrix(1,p,p)

  Lambda <- array(0,c(p,p,K))
  for(k in 1:K){
   for(j in 1:p){  
    Lambda[j,j,k] <- F_tau[j,j,taus[k]+1]
   }
  }

  Sum_lam <- rep(0,p)
  for(j in 1:p){
    for(k in 1:K){
      Sum_lam[j] <- Sum_lam[j]+Lambda[j,j,k]^2
    }
  }   

  P <- matrix(0,p,p)
  ord <- order(Sum_lam,decreasing=TRUE)
  for(j in 1:p){
    P[j,ord[j]] <- 1
  }  
  
  for(k in 1:K){
    Lambda[,,k] <- tcrossprod(crossprod(t(P),Lambda[,,k]),P)
  }

  for(i in 1:(M+1)){
    F_tau[,,i] <- tcrossprod(crossprod(t(P),F_tau[,,i]),P)
  }

  W <- crossprod(t(P),W)

  ASCOV <- matrix(0,p^2,p^2)
  for(j in 1:p){
   for(i in 1:p){
     if(j<i){
       ASV <- .C("ascov_deflji", as.double(as.vector(F_tau)),as.double(as.vector(Lambda)),as.double(taus),as.integer(c(i-1,j-1,p,M,K)),as.double(as.vector(Beta)), res=double(2),PACKAGE="BSSasymp")$res
      
      ASCOV <- ASCOV+ASV[2]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,j]),tcrossprod(diag(p)[,j],diag(p)[,i]))+ASV[1]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,i]),tcrossprod(diag(p)[,j],diag(p)[,j]))
     }  
     
     if(i==j){
       ASCOV <- ASCOV+0.25*D_lm(F_tau,0,0,Beta)[i,i]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,i]),tcrossprod(diag(p)[,i],diag(p)[,i]))
     }
     
     if(i<j){
       ASV <- .C("ascov_deflij", as.double(as.vector(F_tau)),as.double(as.vector(Lambda)),as.double(taus),as.integer(c(i-1,j-1,p,M,K)),as.double(as.vector(Beta)), res=double (2),PACKAGE="BSSasymp")$res
       
       ASCOV <- ASCOV+ASV[2]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,j]),tcrossprod(diag(p)[,j],diag(p)[,i]))+ASV[1]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,i]),tcrossprod(diag(p)[,j],diag(p)[,j]))
     }
    
   }
  }  

  A <- solve(W)
  COV_A <- crossprod(t(tcrossprod(kronecker(diag(p),A),ASCOV)),kronecker(diag(p),t(A)))/T  
  COV_W <- crossprod(t(tcrossprod(kronecker(t(W),diag(p)),ASCOV)),kronecker(W,diag(p)))/T
  
  list(W=W, COV_W=COV_W, A=A, COV_A=COV_A)
}



ASCOV_SOBIdefl_est <- function(X, taus, arp=NULL, maq=NULL, mixed=TRUE, M=100, ...)
{
  p <- dim(X)[2]
  T <- dim(X)[1]
  K <- length(taus)  
 
  if(mixed){
    W <- SOBI(X,taus,method="djd")$W
  }else W <- diag(p)
    
  X <- tcrossprod(sweep(X,2,colMeans(X)),W)  

  if(is.null(arp)) arp <- rep(3,p) 
  if(is.null(maq)) maq <- rep(4,p)

  Psi <- matrix(0,M+1,p)
  for(i in 1:p){
    arma <- arima(X[,i],c(arp[i],0,maq[i]),...)
    psi0 <- ARMAtoMA(ar=arma$coef[min(1,arp[i]):arp[i]],
              ma=arma$coef[(arp[i]+1):(arp[i]+maq[i])],M)
    Psi[,i] <- c(1,psi0)
    Psi[,i] <- Psi[,i]/sqrt(sum(Psi[,i]^2))
  }

  Beta <- matrix(0,p,p)
  for(i in 1:p){
    Ex4 <- mean(X[,i]^4)
    SumPsi4 <- sum(Psi[,i]^4)
    Beta[i,i] <- (Ex4-3)/SumPsi4+3
  }  
  
  for(i in 1:(p-1)){
    for(j in (i+1):p){
      Exi2xj2 <- mean(X[,i]^2*X[,j]^2)
      SumPsii2j2 <- sum(Psi[,i]^2*Psi[,j]^2)
      Beta[i,j] <- (Exi2xj2-1)/SumPsii2j2+1    
      Beta[j,i] <- Beta[i,j]
    }  
  }

  PSI <- array(0,c(p,p,M+1)) 
  for(i in 1:p){
   for(j in 1:(M+1)){
     PSI[i,i,j] <- Psi[j,i]
   }  
  }

  F_tau <- array(0,c(p,p,M+1))
  for(i in 0:M){
   for(j in 1:(M+1-i)){
      F_tau[,,i+1] <- F_tau[,,i+1]+tcrossprod(diag(PSI[,,j]),diag(PSI[,,j+i]))
   }
  }

  
  for(m in 0:M){
     diag(F_tau[,,m+1]) <- diag(tcrossprod(t(X[1:(T-m),]),t(X[(m+1):T,]))/(T-m))
  }
  
  Lambda <- array(0,c(p,p,K))
  for(k in 1:K){
   for(j in 1:p){  
    Lambda[j,j,k] <- F_tau[j,j,taus[k]+1]
   }
  }
  
  Sum_lam <- rep(0,p)
  for(j in 1:p){
    for(k in 1:K){
      Sum_lam[j] <- Sum_lam[j]+Lambda[j,j,k]^2
    }
  }   

  P <- matrix(0,p,p)
  ord <- order(Sum_lam,decreasing=TRUE)
  for(j in 1:p){
    P[j,ord[j]] <- 1
  }  
  
  for(k in 1:K){
    Lambda[,,k] <- tcrossprod(crossprod(t(P),Lambda[,,k]),P)
  }

  for(i in 1:(M+1)){
    F_tau[,,i] <- tcrossprod(crossprod(t(P),F_tau[,,i]),P)
  }

  W <- crossprod(t(P),W)

  ASCOV <- matrix(0,p^2,p^2)
  for(j in 1:p){
   for(i in 1:p){
     if(j<i){
       ASV <- .C("ascov_deflji", as.double(as.vector(F_tau)),as.double(as.vector(Lambda)),as.double(taus),as.integer(c(i-1,j-1,p,M,K)),as.double(as.vector(Beta)), res=double (2),PACKAGE="BSSasymp")$res
      
      ASCOV <- ASCOV+ASV[2]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,j]),tcrossprod(diag(p)[,j],diag(p)[,i]))+ASV[1]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,i]),tcrossprod(diag(p)[,j],diag(p)[,j]))

     }  
     
     if(i==j){
       ASCOV <- ASCOV+0.25*D_lm(F_tau,0,0,Beta)[i,i]*kronecker(tcrossprod(diag(p)[,i], diag(p)[,i]),tcrossprod(diag(p)[,i],diag(p)[,i])) 
     }
     
     if(i<j){
       ASV <- .C("ascov_deflij", as.double(as.vector(F_tau)),as.double(as.vector(Lambda)),as.double(taus),as.integer(c(i-1,j-1,p,M,K)),as.double(as.vector(Beta)), res=double (2),PACKAGE="BSSasymp")$res
       
       ASCOV <- ASCOV+ASV[2]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,j]),tcrossprod(diag(p)[,j],diag(p)[,i]))+ASV[1]*kronecker(tcrossprod(diag(p)[,i],diag(p)[,i]), tcrossprod(diag(p)[,j],diag(p)[,j]))
     }
    
   }
  }  
  
  A <- solve(W)
  COV_A <- crossprod(t(tcrossprod(kronecker(diag(p),A),ASCOV)),kronecker(diag(p),t(A)))/T  
  COV_W <- crossprod(t(tcrossprod(kronecker(t(W),diag(p)),ASCOV)),kronecker(W,diag(p)))/T
  
  list(W=W, COV_W=COV_W, A=A, COV_A=COV_A)
}