confirming a theorem regarding the asymptotic expectation and variance blah blah blah

confirming.r

library(matrixcalc)

# linear algebra utilities
# quadratic form x' A x
qform <- function(A,x) { t(x) %*% (A %*% x) }
# quadratic form x A x'
qoform <- function(A,x) { qform(A,t(x)) }
# outer gram: x x'
ogram <- function(x) { x %*% t(x) }
# A kron A
AkronA <- function(A) { kronecker(A,A,FUN='*') }
# matrix trace
matrace <- function(A) { matrixcalc::matrix.trace(A) }

# duplication, elimination, commutation, symmetrizer
# commutation matrix;
# is p^2 x p^2
Comm <- function(p) {
  Ko <- diag(p^2)
  dummy <- diag(p)
  newidx <- (row(dummy) - 1) * ncol(dummy) + col(dummy)
  Ko[newidx,,drop=FALSE]
}
# Symmetrizing matrix, N
# is p^2 x p^2
Symm <- function(p) { 0.5 * (Comm(p) + diag(p^2)) }
# Duplication & Elimination matrices
Dupp <- function(p) { matrixcalc::duplication.matrix(n=p) }
Elim <- function(p) { matrixcalc::elimination.matrix(n=p) }
# vector function
fvec <- function(x) {
  dim(x) <- c(length(x),1)
  x
}

# compute Theta from mu, Sigma
make_Theta <- function(mu,Sigma) {
  stopifnot(nrow(Sigma) == ncol(Sigma),
            nrow(Sigma) == length(mu))
  mu_twid <- c(1,mu)
  Sg_twid <- cbind(0,rbind(0,Sigma))
  Theta   <- Sg_twid + ogram(mu_twid)
  Theta_i <- solve(Theta)
  zeta_sq <- Theta_i[1,1] - 1
  list(pp1=nrow(Theta),
       p=nrow(Theta)-1,
       mu_twid=mu_twid,
       Sg_twid=Sg_twid,
       Theta=Theta,
       Theta_i=Theta_i,
       zeta_sq=zeta_sq,
       zeta=sqrt(zeta_sq))
}

# construct four parts of Omega_0 from the Theorem;
Omega_bits <- function(mu,Sigma,kurtf=1) {
  tvals <- make_Theta(mu,Sigma)

  Nmat <- Symm(tvals$pp1)
  P1 <- (kurtf-1) * ogram(fvec(tvals$Sg_twid))
  P2 <- (kurtf-1) * 2 * Nmat %*% AkronA(tvals$Sg_twid)
  P3 <- 2 * Nmat  %*% AkronA(tvals$Theta)
  P4 <- 2 * Nmat  %*% AkronA(ogram(tvals$mu_twid))
  list(P1=P1,P2=P2,P3=P3,P4=P4)
}
Omega_0 <- function(mu,Sigma,kurtf=1) {
  obits <- Omega_bits(mu=mu,Sigma=Sigma,kurtf=kurtf)
  obits$P1 + obits$P2 + obits$P3 - obits$P4
}
# construct matrices F and H and vector h 
FHh_values <- function(mu,Sigma,R=1) {
  tvals <- make_Theta(mu,Sigma)
  zeta_sq <- tvals$zeta_sq
  zeta <- sqrt(zeta_sq)
  mp <- t(t(-tvals$Theta_i[2:tvals$pp1,1]))

  Fmat <- (1 / R^2) * ((ogram(mu) / zeta) - (zeta * Sigma))
  H1 <- - cbind(cbind((1 / (2 *zeta_sq)) * mp, 
                      (R / zeta) * diag(tvals$p)),
                matrix(0,nrow=tvals$p,ncol=tvals$p * tvals$pp1))
  H2 <- - AkronA(tvals$Theta_i)
  Hmat <- H1 %*% H2 %*% Dupp(tvals$pp1)

  hvec <- - AkronA(tvals$Theta_i[1,,drop=FALSE]) %*% Dupp(tvals$pp1)
  list(Fmat=Fmat, Hmat=Hmat, hvec=hvec)
}

# compute the expected bias and variance
# in two ways, one directly from the identity of
# H, F, h and Omega_0
# the other from the theorem
# 
# then compute relative errors of the theorem's approximation.
testit <- function(mu,Sigma,kurtf=1,R=1) {
  p <- length(mu)
  tvals <- make_Theta(mu,Sigma)

  OmegMat <- Omega_0(mu,Sigma,kurtf=kurtf)
  Omega <- qoform(A=OmegMat,Elim(p+1))
  FHh <- FHh_values(mu,Sigma,R=R)

  # now test corollary 'true' values
  HtFH <- qform(A=FHh$Fmat,x=FHh$Hmat)
  Vari <- as.numeric(qoform(Omega,FHh$hvec) / (4 * tvals$zeta_sq))
  Eb1 <- as.numeric((1/2) * matrace((HtFH) %*% Omega))
  Eb2 <- Vari / (2 * tvals$zeta)
  Ebias <- Eb1 + Eb2
  
  # Theorem says:
  T_Vari  <- (((3 * kurtf - 1) / 4) * tvals$zeta_sq + 1) 
  T_Eb1 <- ((kurtf * tvals$zeta_sq + 1) * 
            (1 - tvals$p)) / (2 * tvals$zeta)
  T_Ebias <- ((kurtf * tvals$zeta_sq + 1) * 
              (1 - tvals$p) + T_Vari) / (2 * tvals$zeta)

  ERR_Ebias <- T_Ebias - Ebias
  ERR_Vari <- T_Vari - Vari
  max(c(max(abs(ERR_Ebias)) / (abs(Ebias)),
        max(abs(ERR_Vari)) / (abs(Vari))))
}

# create a population
myp <- 5
set.seed(1234)
mu <- matrix(rnorm(myp,1),ncol=1)
ZZ <- matrix(rnorm(100,myp),ncol=myp)
Sigma <- t(ZZ) %*% ZZ

# test it:
testit(mu,Sigma,kurtf=2,R=1)