*
* @MVQStat( options )  start end
* # list of series
*
* Computes the Hosking variant on the multivariate Q statistic
*
* Parameters:
*   start  end    Range over which to compute statistic. By default, the
*                 common range of the listed series.
*
* Options:
*   LAGS=# of lags [12]
*   DFC=degrees of freedom correction [0]
*
*   PRINT/NOPRINT
*   TITLE=title of report ["Multivariate Q Statistic"]
*
* Variables Defined:
*   %CDSTAT    test statistic
*   %SIGNIF    significance level treating statistic as chi-squared
*   %NDFQ      degrees of freedom of the chi-squared
*   %NOBS      number of observations
*
* Reference:
*   Hosking(1981), "Equivalent Forms of the Multivariate Portmanteau
*     Statistic", JRSS-B, vol 43, no. 2, pp 261-262.
*
* Revision Schedule:
*   05/2005 Written by Tom Doan, Estima
*   03/2018 Revise to create REPORT. Add PRINT, TITLE options
*
procedure MVQStat start end
type integer start end
*
option integer lags  12
option integer dfc   0
option switch  print 1
option string  title
*
local vect[integer]  list
local integer        startl endl
local integer        i j l n
local vect[rect]     mcorr
local series         ccorr
local symm           ginv
local real           qstat
local report         qreport
local string         ltitle
*
* Get the list of variables
*
enter(varying,entries=n) list
*
* Prepare the matrices for the correlations, one nxn matrix
* of each lag.
*
dim mcorr(lags+1)
ewise mcorr(i)=%zeros(n,n)
*
* Set the sample range
*
inquire(reglist) startl<<start endl<<end
# list
*
* Use CROSS to compute the pairwise autocorrelations, and stuff
* them into the proper slots in the MCORR matrices.
*
do i=1,n
   do j=1,n
      cross(noprint) list(i) list(j) startl endl 0 lags ccorr
      do l=0,lags
         compute mcorr(l+1)(i,j)=ccorr(l+1)
      end do l
   end do j
end do i
*
* Compute the test statistic
*
compute ginv=inv(mcorr(1))
compute qstat=0.0
do l=1,lags
   compute qstat=qstat+%dot(%mqform(ginv,mcorr(l+1)),ginv)/(%nobs-l)
end do li
compute qstat=qstat*%nobs^2
compute %cdstat=qstat
compute %ndfq  =lags*n*n-dfc
compute %signif=%chisqr(qstat,%ndfq)
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Multivariate Q Test"
*
report(use=qreport,action=define,title=ltitle)
report(use=qreport,atrow=1,atcol=1,span) ltitle
report(use=qreport,atrow=2,atcol=1,span) "Test Run Over "+%datelabel(startl)+" to "+%datelabel(endl)
report(use=qreport,atrow=3,atcol=1) "Lags Tested"  lags
report(use=qreport,atrow=4,atcol=1) "Degrees of Freedom" %ndfq
if dfc>0
   report(use=qreport,row=new,atcol=1) "D of F Correction" dfc
report(use=qreport,row=new,atcol=1) "Q Statistic" qstat
report(use=qreport,row=new,atcol=1) "Signif Level" %signif
report(use=qreport,action=format,atcol=2,width=8,align=decimal)
if print
   report(use=qreport,action=show)
end
