*
* @HinichTest(options) x start end
*
* Performs the Hinich bispectrum test for linearity and Gaussianity
* as described in Hinich and Patterson (1989), "Evidence of
* nonlinearity in the trade-by-trade stock market return generating
* process", in Economic Complexity: Chaos, Sunspots, Bubbles and
* Nonlinearity, Barnett, Geweke and Shell, eds. Cambridge
* University Press.
*
* Note that this requires quite a bit of memory if you get above a
* few thousand data points.
*
* Parameters:
*   x            series to test
*   start  end   range to use (by default, the defined range of <<x>>)
*
* Options:
*   FRAME=separation of lattice points [square root of number of data points]
*     The input or computed value is rounded up to an odd number.
*   QUANTILE=quantile of the empirical distribution examined in the linearity
*            test [.8]
*
* Revision Schedule:
*  05/2005 Written by Tom Doan, Estima
*  04/2006 QUANTILE option added
*  02/2018 Add check for negative non-centrality on linearity test.
*
procedure HinichTest x start end
type series x
type integer start end
*
option integer frame
option real    quantile   .8
*
local integer endl startl
local integer n m p j
local integer s1 s2 s3 r1 r2 w1 w2 base half fill
local real scalebi scaleuni uni1 uni2 uni12 ssqr
local rect[complex] bispect
local vector unispect skewsq
local complex bispects
local real xihat lambda0 mean var lstat
*
* Get the sample range, and Fourier transform the data
*
inquire(series=x) startl<<start endl<<end
compute n=endl-startl+1
freq 1 n
dim bispect(n,n) unispect(n)
rtoc startl endl 1
# x
# 1
fft 1
*
* Check the frame size. Round up to an odd number
*
if .not.%defined(frame)
   compute m=fix(sqrt(n))/2*2+1
else
   compute m=frame/2*2+1
*
* Compute the raw bispectrum and raw spectrum
*
compute scalebi=1.0/(n*(2*%pi)**2),scaleuni=1.0/(n*(2*%pi))
do s1=1,n
   do s2=1,n
      compute s3=%clock(s1+s2-1,n)
      compute bispect(s1,s2)=scalebi*%z(s1,1)*%conjg(%z(s2,1))*%z(s3,1)
   end do s2
   compute unispect(s1)=scaleuni*%cabs(%z(s1,1))**2
end do s1
*
* Count the number of test frequencies.
*
compute p=0
compute base=(m+1)/2
compute half=(m-1)/2
do s1=base,n/2,m
   do j=s1+m,n/2-(s1+1),m
      compute p=p+1
   end do j
end do i
*
* Smooth the bispectrum at the lattice points and the spectrum at
* each component and their sum.
*
dim skewsq(p)
*
compute fill=0
do s1=base,n/2,m
   compute uni1=0.0
   do r1=-half,half
      compute uni1=uni1+unispect(%clock(s1+r1,n))
   end do r1
   compute uni1=uni1/m

   do s2=s1+m,n/2-(s1+1),m
      compute uni2=0.0
      do r2=-half,half
         compute uni2=uni2+unispect(%clock(s2+r2,n))
      end do r2
      compute uni2=uni2/m
      compute uni12=0.0
      compute s3=s1+s2-1
      do r2=-half,half
         compute uni12=uni12+unispect(%clock(s3+r2,n))
      end do r2
      compute uni12=uni12/m

      compute bispects=0.0
      do r1=-half,half
         compute w1=%clock(s1+r1,n)
         do r2=-half,half
            compute w2=%clock(s2+r2,n)
            compute bispects=bispects+bispect(w1,w2)
         end do r2
      end do r1
      compute bispects=bispects/m**2
      compute fill=fill+1
      *
      * Compute the skew statistic normalized by the variance of
      * the estimator of the |bispectrum|**2. (There's an extra 2*pi
      * factor relative to the formula in Hinich and Patterson
      * because of the 2*pi scale factors in the estimates of the
      * raw bispectrum and spectrum).
      *
      compute skewsq(fill)=(2.0*%pi)*(2.0*m**2/n)*%cabs(bispects)**2/(uni1*uni2*uni12)
   end do s2
end do s1
*
* Test statistics
*
* For Gaussianity (that is, linearity plus Gaussian joint
* distribution), this is just the sum of the skew statistics (which
* already have the 2 multiplier). This has an asymptotic central
* chisquare distribution with 2p degrees of freedom.
*
* For linearity, the .8 quantile of the empirical distribution of
* the skew statistics is compared with that for independent draws
* from a non-central chisquare with non-centrality parameter
* estimated from the observed values.
*
compute ssqr=%sum(skewsq)
compute xihat=%scalar(%fractiles(skewsq,||quantile||))
compute lambda0=%max(%avg(skewsq)-2,0)
compute mean=%invchisqrnc(quantile,2,lambda0)
compute var =quantile*(1-quantile)/(p*%chisqrncdensity(mean,2,lambda0))
compute lstat=(xihat-mean)/sqrt(var)
report(action=define)
report(row=new,atcol=1,span) "Hinich Bispectral Tests"
report(row=new,atcol=1) "Observations" n
report(row=new,atcol=1) "Frame Size" m
report(row=new,atcol=1) "Lattice Points" p
report(row=new,atcol=1) "Test Quantile" quantile
report(row=new,atcol=1) "Linearity" lstat 1-%cdf(lstat)
report(row=new,atcol=1) "Gaussianity" ssqr %chisqr(ssqr,2*p)
report(action=format,atcol=2,tocol=3,width=10)
report(action=show)
end
