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

[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])  • I don't see anything in your code that would correspond to simulating data according to the null hypothesis of a Normal distribution. Could you show us where that would be? – whuber Aug 5, 2015 at 15:09 • In the link, @Glen_b shows that the distributions I am sampling from (f.ex. (a), a 50-50 mixture of$U[-1,1]$and$U[-a,a]$) are symmetric distributions with excess kurtosis 0, s.th. I believe the data to come from a distribution satisfying the actual null hypothesis of the JB test, viz. symmetry and zero excess kurtosis. Hence, I expect the p-values of the JB test to be$U[0,1]$. Aug 5, 2015 at 15:12 • The p-value for the JB test is computed assuming the underlying distribution is Normal. No p-value can be computed for the nonparametric assumption of zero skewness and zero excess kurtosis, because under that assumption the statistic has no definite distribution. – whuber Aug 5, 2015 at 15:40 • My understanding is that the$\chi^2$limiting distribution is a consequence of a CLT, CMT and delta method applied to sample skewness and kurtosis. As such, the distribution of the data in the population should, provided the CLT works, not matter. I am struggling to find a reference for my claim, though. Aug 5, 2015 at 15:56 • That argument might have merit. But notice that the J-B statistic involves up to the eighth moment of the sample (it includes the square of the sample kurtosis). Such moments cannot be expected to stabilize until sample sizes are enormous.$10^5$is likely far too small. This is one reason the J-B test performs so poorly in the first place. – whuber Aug 5, 2015 at 16:23 ## 2 Answers 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!) • Even when the Jarque Bera test is applied to the normal, the asymptotic distribution doesn't really "kick in" until sample sizes are very large. My own simulations of the bivariate distribution showed its quite clearly not bivariate normal until$n$is well over 500, and I'd hesitate to use the resulting univariate statistic until$n\$ is at least 300. Bowman and Shenton's 1975 paper ( "Omnibus contours for departures from normality based on √b1 and b2". Biometrika 62 (2): 243–250. ) -- five years before Jarque and Bera -- discuss the problem with the (widely misnamed) statistic and ... ctd Aug 7, 2015 at 21:18
• ctd ... describe small sample improvements. However, Fig 1 in their paper suggests either their approximation goes the wrong way (the contours of the approximation are a worse fit than the asymptotic ones) -- or perhaps an error was made in drawing the figure. In any case, the dot diagram itself (easily simulated) shows the characteristic "smile" shape of a sample (skewness,kurtosis) plot seen for many of the common symmetric distributions. The asymptotic test itself was widely known before Bowman and Shenton's paper, so it's odd that it's named after authors of a paper published much later. Aug 7, 2015 at 21:32

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 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)