*
* @MVJB( options)  u start end
*
* Computes a multivariate version of the Jarque-Bera test for normality.
* Note that there are more sophisticated versions of this (for instance,
* Doornik and Hansen, "An Omnibus Test for Univariate and Multivariate
* Normality"). This just transforms the input residual series to
* uncorrelated components (using an eigenbased factorization if not
* provided by the user) and sums up the univariate JB statistics from
* those.
*
* Parameters:
*   u          VECTOR of SERIES
*   start end  Range for the calculation [by default, the maximum
*              combined range of the u's]
*
* Options:
*   SIGMA=SYMMETRIC covariance matrix of u
*   FACTOR=RECTANGULAR factor of the covariance matrix of u.
*     One of these two must be supplied.
*
* Revision Schedule
*  03/2005 Written by Tom Doan, Estima
*  10/2010 Add title to report.
*
procedure MVJB u start end
type vect[series] u
type integer start end
*
option rect factor
option symm sigma
*
local vect[series] v
local integer i startl endl
local rect eigfac invfact
local integer n
local real jbtotal
*
inquire(reglist) startl<<start endl<<end
# u
*
* Check that we have either sigma or factor
*
if .not.%defined(sigma).and..not.%defined(factor) {
   disp "Need either SIGMA or FACTOR option"
   return
}
*
* Get the inverse factor of the covariance matrix
*
if .not.%defined(factor) {
   eigen(scale) sigma * eigfac
   compute invfact=inv(eigfac)
}
else
   compute invfact=inv(factor)
*
* Transform the input residuals to uncorrelated components
*
compute n=%rows(u)
dim v(n)
do i=1,n
   set v(i) startl endl = 0.0
end do i
*
do t=startl,endl
   compute %pt(v,t,invfact*%xt(u,t))
end do t
*
* Create a table with the individual components and the overall
* statistic
*
report(action=define,hlabels=||"Var","JB","P-Value"||)
compute jbtotal=0.0
do i=1,n
   stats(noprint) v(i) startl endl
   report(row=new,atcol=1) i %jbstat %jbsignif
   compute jbtotal=jbtotal+%jbstat
end do i
report(row=new,atcol=1) "All" jbtotal %chisqr(jbtotal,2*n)
report(action=format,atcol=2,tocol=3,picture="*.###")
report(action=show)
end

