*
* @BDStest( options)   series  start  end
* performs the Brock, Dechert, Scheinkman test for i.i.d.
*
* Parameters:
*   series        series to analyze
*   start  end    range to use. By default, the defined range of series
*
* Options:
*   P  =multiple of standard error (or sample range) to use for the closeness threshold
*   DIM=embedding dimension [depends upon number of data points]
*   RANGE/[NORANGE]
*     If RANGE, the threshold is the multiple P of the observed range (max-min) of the
*     series. In that case, P should be less than 1.
*   [PRINT]/NOPRINT
*
* Variables Defined:
*   %CDSTAT    The BDS test statistic (asymptotically N(0,1))
*   %SIGNIF    The asymptotic p-value
*
* Reference:
*   Brock, Dechert, & Scheinkman(1993), "A test of independence based on the
*    correlation dimension" (working paper versions), included in "Nonlinear
*    Dynamics,Chaos,and Instability", MIT Press,1993,Chap.2.
*
* Revision Schedule:
*   Orignally written by Hideaki ISHIGAMI, Nippon Sports Science Univ.
*   08/2004 Updated to version 6. General code cleanup by Tom Doan, Estima.
*   01/2007 Lower minimum number of observations to 100.
*   05/2008 Add RANGE and PRINT options.
*
****************************************************
Procedure BDSTest series start end
type series series
type integer start end

option real    p      1.0
option integer dim
option switch  range  0
option switch  print  1

local integer j r s nobs startl endl m nobs
local real std eps cmt c1t K n var
local real dif indic indics sumh sumhh sumz

inquire(series=series) startl<<start endl<<end


********** Sets eps and M ********************

compute std=eps=n=0.
compute m=0

statistics(noprint,fract) series startl endl
*
* eps is the threshold value for the indicator function. If <<range>>, it is set
* as a fraction of the series range. If not (which is the  default, it is set as a
* multiple of the series standard deviation. It's recommended that the multiple be
* between .5 and 2.0
*
if range
   compute eps=p*(%maximum-%minimum)
else {
   compute std=sqrt(%variance)
   compute eps=p*std
   }

compute nobs=endl-startl+1

if nobs<=100 {
   dis "Sample of" nobs "too small to do BDS test!"
   return
}
if %defined(dim)
   compute m=dim
else
if nobs<=600
   compute m=2 ; * T/m should be greater than 2
else
if nobs<=800
   compute m=3
else
if nobs<=1000
   compute m=4
else
   compute m=5

compute n=nobs-m+1


************ Computation of Correlation Integral ****************

*
* Computation of C_mt(eps)
*

compute indics=0.
do r=startl,endl-m
   do s=r+1,endl-m+1
      do j=0,m-1
         compute dif =abs(series(r+j)-series(s+j))
         if dif>eps
            break
      end do j
      compute indics=indics+(dif<=eps)
   end do s
end do r
*
compute cmt=indics/(n*(n-1)/2.)

*
* Computation of C_1t
*

compute indic=0.
do r=startl,endl-1
   do s=r+1,endl
      compute indic=indic+(abs(series(r)-series(s))<eps)
   end do s
end do r
*
compute c1t=indic/(nobs*(nobs-1)/2.)


*********** Computation of Variance of D_mt*********************

*
*  Computation of K
*

compute sumhh=0.0
do r=startl,endl
   compute sumh=0.0
   do s=startl,endl
      compute sumh  = sumh+(abs(series(r)-series(s))<eps)
   end do s
   compute sumhh=sumhh+(sumh**2)
end do r
*
compute K=sumhh/(nobs**3.)


*
*  Computation of Variance
*

compute sumz=0.
do j=1,m-1
   compute sumz=sumz+(k**(m-j))*(c1t**(2*j))
end do j

compute var=4*((k**m)+(2*sumz)+((m-1)**2)*(c1t**(2*m)) $
              -(m**2)*k*(c1t**(2*m-2)))


************** Computation of BDS statistic ********************

compute %cdstat=sqrt(nobs)*(cmt-c1t**m)/sqrt(var)
cdf(noprint) normal %cdstat

if .not.print
   return

dis
dis
dis
dis @6 "************************************"
dis @12 "BDS Test for" %label(series)
dis @8 "************************************"
dis @11 "     P =" p
dis @11 "   eps =" eps
dis @11 "     m =" @+6 # m
dis @11 "  C_mT =" cmt
dis @11 "   C_T =" c1t
dis @11 "     K =" k
dis @11 "sigma2 =" var
dis @11 "  BDS  =" %cdstat
dis @11 "signif =" %signif
dis @6 "************************************"
dis
dis
dis
End BDSTest


