Testing for normality in non-normal distributions with zero skewness and zero excess kurtosis

Question

[This question was formerly called "On Non-normal distributions with zero skewness and zero excess kurtosis" and relabeled to better reflect its focus.]

I am trying to write a little simulation using @Glen_b's very nice answer in this post: Non-normal distributions with zero skewness and zero excess kurtosis?

In particular, I want to show that (what is known in econometrics as) the Jarque-Bera test (which was actually considered earlier in for example Bowman and Shenton (1975). "Omnibus contours for departures from normality based on √b1 and b2". Biometrika 62 (2): 243–250, see @Glen_b's comment) lacks power against non-normal but symmetric distributions without excess kurtosis. The figure shows the distribution of the p-values of the simulation.

For his uniform (left panel) and Poisson (right panel) example I do get a distribution of the p-values that leads to no power beyond size/a conservative test (depending on whether you call JB a test of normality or a test of the two moments), but in the gamma example (middle panel) there does even seem to be some power.

In neither case do I get a uniform distribution of the p-values though, although I (believe to) simulate data under (what I think is) the null - symmetry and no excess kurtosis.

Thoughts on how/why that happens?

CODE:

library(tseries)
library(MASS)

n <- 1e5
a <- sqrt(5+sqrt(24))
b <- (sqrt(13)+1)/2
lambda <- .5

reps <- 1000
JBpval <- matrix(rep(NA,3*reps),ncol=3)

for (i in 1:reps) {
#(a)
u1a <- runif(n/2,-1,1)  
u2a <- runif(n/2,-a,a)

#(b)
u1b <- rgamma(n/2,shape = b, scale = 1)  
u2b <- -rgamma(n/2,shape = b, scale = 1)

#(c)
u1c <- sqrt(rpois(n/2,lambda = lambda))  
u2c <- -sqrt(rpois(n/2,lambda = lambda))

ua <- c(u1a,u2a)
ub <- c(u1b,u2b)
uc <- c(u1c,u2c)

JBpval[i,1] <- jarque.bera.test(ua)$p.value
JBpval[i,2] <- jarque.bera.test(ub)$p.value
JBpval[i,3] <- jarque.bera.test(uc)$p.value
}

par(mfrow=c(1,3))
truehist(JBpval[,1])  
truehist(JBpval[,2])  
truehist(JBpval[,3])  

Answer 1

Let my try to answer my own question:

If $X\sim\mathcal{N}\left(0,\sigma^2\right)$ (zero mean is w.l.o.g.). Define the centered moment $\mu_k=E\left(X-E X\right)^k$ and its empirical counterpart $m_k=\frac{1}{n}\sum_{i=1}^{n}\bigl(X_i-\overline{X}\bigr)^k$.

Analyze the two (method of moments) estimators $$ b_1=\frac{m_3}{m_2^{3/2}}\qquad\text{and}\qquad b_2=\frac{m_4}{m_2^2}, $$ for the skewness $\beta_1(X)=E(X-E(X))^3/\left[E(X-E(X))^{2}\right]^{3/2}$ and for the kurtosis $\beta_2(X)=E(X-E(X))^4/\left[E(X-E(X))^2\right]^2$. A CLT will ensure that $$ \sqrt{n}\left(\begin{pmatrix}m_2\\m_3\\m_4\end{pmatrix}-\begin{pmatrix}\mu_2\\ \mu_3\\ \mu_4\end{pmatrix}\right) $$ converges in distribution to a centered multivariate normal with covariance matrix $\Sigma$.

Now, we want to show that, when $X$ is normal, $$ \sqrt{n}\left(\begin{pmatrix}b_1\\b_2\end{pmatrix}-\begin{pmatrix}\beta_1\\ \beta_2\end{pmatrix}\right)\to_d\mathcal{N}\left(\begin{pmatrix}0\\0\end{pmatrix},\begin{pmatrix}6 & 0\\0 & 24\end{pmatrix}\right). $$ The follow result is helpful (see, e.g., Rao C.R., Linear Statistical Inference and its applications (1972), Section 6h2): $$\text{Acov}(\sqrt{n}m_j,\sqrt{n}m_k)=\mu_{j+k}-\mu_j\mu_k+jk\mu_2\mu_{j-1}\mu_{k-1}-j\mu_{j-1}\mu_{k+1}-k\mu_{k-1}\mu_{j+1},$$ For the normal distribution, $$ \mu_{2k}=\sigma^{2k}\frac{(2k)!}{k!2^k},\qquad\mu_{2k+1}=0,\qquad k=0,1,\ldots. $$ Thus, $\mu_2=\sigma^2$, $\mu_3=0$, $\mu_4=3\sigma^4$, $\mu_5=0$, $\mu_6=15\sigma^6$ and so forth.

The limiting distribution of $(b_1,b_2)$ then is an application of the multivariate delta method.

The $(2\times3)$ matrix of derivatives of $b_1$ and $b_2$ w.r.t. $m_2$, $m_3$, $m_4$, evaluated at the population moments is $$ J=\begin{pmatrix}0&\sigma^{-3}&0\\-6\sigma^{-2}&0&\sigma^{-4}\end{pmatrix} $$ Also, $$ \Sigma=\begin{pmatrix}2\sigma^4&*&15\sigma^6-\sigma^23\sigma^4\\*&6\sigma^6&*\\15\sigma^6-\sigma^23\sigma^4&*&105\sigma^8-(3\sigma^4)^2\end{pmatrix}, $$ where $*$ omits terms not needed (because they will be multiplied with zeros in $J$ when evaluating the $(2\times 2)$ variance matrix of interest $J\Sigma J'$).

The variance-covariance matrix of the delta method then gives $$ J\Sigma J'=\begin{pmatrix}6 & 0\\0 & 24\end{pmatrix}, $$ as desired.

Thus, the limiting null distribution of the Jarque-Bera test $$ JB=n\left(\frac{b_1^2}{6}+\frac{(b_2-3)^2}{24}\right) $$ follows directly because, under the null of normality, $\sqrt{n}(b_1/\sqrt{6})\to_d N(0,1)$ and $\sqrt{n}((b_2-3)/\sqrt{24})\to_d N(0,1)$ such that $n(b_1^2/6)\to_d \chi^2_1$ and $n(b_2-3)^2/24\to_d\chi^2_1$. By asymptotic independence, $$ JB\to_d\chi^2_2 $$ But, and that was the point of the original post, these arguments inherently required normality for the limiting distribution to obtain.

Hence, even if distributions share the skewness and kurtosis of the normal, there is no reason to believe that the JB test will be $\chi^2_2$, and thus, no reason to believe that the $p$-values will be uniform. So, my intuition was false.

(As also mentioned by @whuber in the comments, what would be required would be a distribution that shares the first eight moments with the normal, as $\text{Acov}(\sqrt{n}m_4,\sqrt{n}m_4)$ contains $\mu_{8}$. I am not aware of such a distribution though - examples would be appreciated, if existent!)

Answer 2

It's not so much that the test would be without any power, but you might find it interesting to see what happens to the relative power as $n$ increases when you compare the "Jarque-Bera" with (say) Shapiro-Wilk. Actually, here you go:

The light grey line marks the type I error rate; below that (with a false null) we have a biased test (biased tests are common with goodness of fit).

[Note to self: I should make a beta example that does the same trick as the gamma one. That comparison should make for a more spectacular comparison.]

R code, if anyone cares to have it --

library(tseries)
f <- function(x) mean(x<.05)

reps <- 5000
b <- (sqrt(13)+1)/2

nseq <- c(16,25,64,100,225,400,625)
nres <- matrix(rep(NA,length(nseq)*2),nc=2)

for (i in seq_along(nseq)){
 res <- replicate(reps,{
     g1 <- rgamma(nseq[i],b)*sign(runif(nseq[i])-.5)
     c(jarque.bera.test(g1)$p.value,shapiro.test(g1)$p.value)
   })
 nres[i,] <- apply(res,1,f)
}

and code to do the plot:

plot(qnorm(nres[,1])~sqrt(nseq),ylim=c(qnorm(.02),qnorm(.997)),
         type="o",col=4,ylab="zscore of power",xlab="sqrt(n)")
points(qnorm(nres[,2])~sqrt(nseq),type="o",col="orange3")
abline(h=qnorm(0.05),col=8)