*
* @FLUX( options ) start end
* # list of score series
*
* Computes the Nyblom fluctuations test on a set of series of scores. Computes
* both individual and joint test statistics.
*
* Parameters:
*   start   end    Range over which to compute test (by default the common range
*                  of the score series.
*
* Options:
*   [PRINT]/NOPRINT
*     Controls whether or not the stability test results are printed.
*   SMPL=SMPL Series
*
* Variables Defined:
*   %HJOINT  real with the overall test
*   %HSTATS  VECTOR with the tests on the individual components
*
* Reference:
*   Nyblom(1989), "Testing for Constancy of Parameters Over Time", JASA, vol 84, pp
*   223-230.

* Revision Schedule:
*   03/2006 Written by Tom Doan, Estima
*
procedure flux start end
type integer   start end
*
option switch print 1
option series smpl

local series     smpll
local integer    i startl endl time n
local vector     s x
local symm       vinv
local equation   fluxeq
*
declare vector %HSTATS
declare real   %HJOINT
*
equation fluxeq *
inquire(equation=fluxeq,valid=smpll) startl<<start endl<<end
if %defined(smpl)
   set smpll startl endl = smpl
cmom(smpl=smpll,equation=fluxeq) startl endl
compute vinv=inv(%nobs*%cmom)
compute n=%ncmom
*
dim %hstats(n) s(n)
compute %hjoint=0.0,%hstats=%const(0.0),s=%const(0.0)
do time=startl,endl
   if .not.smpll(time)
      next
   compute x=%eqnxvector(fluxeq,time)
   compute s=s+x
   ewise %hstats(i)=%hstats(i)+s(i)**2
   compute %hjoint=%hjoint+%qform(vinv,s)
end do time
ewise %hstats(i)=%hstats(i)/(%nobs*%cmom(i,i))

if .not.print
   return

report(action=define,hlabels=||"Test","Statistic","P-Value"||)
report(atrow=1,atcol=1) "Joint" %hjoint %nyblomtest(%hjoint,n)
do i=1,n
   report(atrow=2+i,atcol=1) i %hstats(i) %nyblomtest(%hstats(i),1)
end do i
report(action=format,atcol=2,tocol=2,width=10)
report(action=format,atcol=3,tocol=3,picture="#.##")
report(action=show,window="Stability Test")
end flux

