Doing a Jarque-Bera Test In R

Skewness of \(x\) is measured as \[S = \frac{\left( E[X - \mu]^{3} \right)^{2}}{\left(E[X - \mu]^{2} \right)^{3}}\]

Because the normal distribution is symmetric, the skewness (deviation from symmetry) should be zero.

Kurtosis of \(x\) is measured as \[\kappa = \frac{E[X - \mu]^{4}}{\left( E[X - \mu]^{2} \right)^{2}}\]

and \(\kappa = 3\) for a normal distribution. The Jarque-Bera statistic is \[jb = T\left[ \frac{S}{6} + \frac{(\kappa - 3)^{2}}{24} \right]\]

where \(T\) is the sample size. Under the null hypothesis of normality, \(jb \sim \chi^{2}(2)\).

The Jarque-Bera test is available in R through the package tseries. Example:

library(tseries)
set.seed(100)
x <- rnorm(5000)
jarque.bera.test(x)

The output is

    Jarque Bera Test

data:  x
X-squared = 0.046, df = 2, p-value = 0.9773

We do not reject the null hypothesis of normality for this series. That is a good thing, otherwise we would want to check if R’s random number generating functions are working properly.

Here is the implementation, with some comments that I’ve added myself:

n <- length(x)         ## Number of observations
m1 <- sum(x)/n         ## Mean
m2 <- sum((x-m1)^2)/n  ## Used in denominator of both
m3 <- sum((x-m1)^3)/n  ## For numerator of S
m4 <- sum((x-m1)^4)/n  ## For numerator of K
b1 <- (m3/m2^(3/2))^2  ## S
b2 <- (m4/m2^2)        ## K
STATISTIC <- n*b1/6+n*(b2-3)^2/24

The test statistic (what I called \(jb\) above) is reported as x.squared (not sure why that name was chosen), the degrees of freedom parameter is always 2, and the p-value is calculated as 1 - pchisq(STATISTIC,df = 2).

Last updated: 2015-10-07