*
* @MVJB( options )  start end
* # list of series
*
* 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:
*   start end  Range for the calculation [by default, the maximum
*              combined range of the listed series]
*
* Options:
*   SIGMA=SYMMETRIC covariance matrix of u
*   FACTOR=RECTANGULAR factor of the covariance matrix of u.
*   HMATRICES=SERIES[SYMM] of time-varying covariance matrices.
*     One of these three must be supplied.
*
*   [PRINT]/NOPRINT
*   TITLE="title of report" ["Multivariate JB Test"]
*
* Variables Defines:
*   %CDSTAT    Joint test statistic
*   %SIGNIF    Significance of joint statistic
*   %NDFTEST   Degrees of freedom of test
*
* Revision Schedule
*  03/2005 Written by Tom Doan, Estima
*  04/2012 Added HMATRICES option
*  03/2018 Add PRINT and TITLE options. Switch to LOCAL REPORT.
*  03/2018 Switch syntax to more standard supplementary card.
*
procedure MVJB start end
type integer start end
*
option rect          factor
option symm          sigma
option series[symm]  hmatrices
option switch        print    1
option string        title
*
local vect[series]   v
local integer        i startl endl
local rect           eigfac invfact
local integer        n
local real           jbtotal
local report         jreport
local string         ltitle
local equation       xeq
*
equation xeq *
compute n=%eqnsize(xeq)
*
inquire(equation=xeq) startl<<start endl<<end
*
* Check that we have sigma, factor or hmatrices
*
if .not.(%defined(sigma).or.%defined(factor).or.%defined(hmatrices)) {
   disp "###@MVJB. Need SIGMA, FACTOR or HMATRICES option"
   return
}
*
* Get the inverse factor of the covariance matrix
*
if %defined(factor)
   compute invfact=inv(factor)
else
if %defined(sigma) {
   eigen(scale) sigma * eigfac
   compute invfact=inv(eigfac)
}
*
* Transform the input residuals to uncorrelated components
*
dim v(n)
do i=1,n
   set v(i) startl endl = 0.0
end do i
*
do t=startl,endl
   if %defined(hmatrices) {
      eigen(scale) hmatrices(t) * eigfac
      compute invfact=inv(eigfac)
   }
   compute %pt(v,t,invfact*%eqnxvector(xeq,t))
end do t
*
* Create a table with the individual components and the overall
* statistic
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Multivariate JB Test"
*
report(use=jreport,action=define,title=ltitle,$
    hlabels=||"Var","JB","P-Value"||)
compute jbtotal=0.0
do i=1,n
   stats(noprint) v(i) startl endl
   report(use=jreport,row=new,atcol=1) i %jbstat %jbsignif
   compute jbtotal=jbtotal+%jbstat
end do i
compute %cdstat=jbtotal,%ndftest=2*n,%signif=%chisqr(%cdstat,%ndftest)
report(use=jreport,row=new,atcol=1) "All" %cdstat %signif
report(use=jreport,action=format,atcol=2,tocol=3,picture="*.###")
if print
   report(use=jreport,action=show)
end
