****************************************************
*  Brock, Dechert, & Scheinkman Test for i.i.d.
*    Brock,W.A., D.A. Hsieh, and B.Lebaron
*      _Nonlinear Dynamics,Chaos,and Instability_, MIT Press,1993,Chap.2
*
*  Syntax:
*
*  @bdstest( options)   series  start  end
*
*  Parameters:
*   series        series to analyze
*   start  end    range to use. By default, the defined range of series
*
*  Options:
*   P=multiple of standard error to use for the closeness threshold
*   DIM=embedding dimension [depends upon number of data points]
*
*  by Hideaki ISHIGAMI
*  Nippon Sports Science Univ.
*
* Revision Schedule:
*  08/04 Updated to version 6. General code cleanup by
*        Tom Doan, Estima.
****************************************************


Procedure BDSTest series start end
type series series
type integer start end

option real    p   1.0
option integer dim

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) series startl endl
compute std=sqrt(%variance)
compute eps=p*std                ;* eps is threshold value for indicator function,
                                  * it should be .5*std < eps < 2*std

compute nobs=endl-startl+1

if nobs<=400 {
   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


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


