*
* @ADTest( options )  x start end
* Anderson-Darling test for normality
*
* Parameters:
*   x           series to analyze
*   start end   range to analyze [range of x]
*
* Options:
*   MEAN=mean of assumed distribution [sample mean]
*   VARIANCE=variance of assumed distribution [sample variance]
*
* Reference: Stephens, M. A. (1974). "EDF Statistics for Goodness of Fit and Some Comparisons",
* Journal of the American Statistical Association, Vol. 69, pp. 730-737
*
* Revision Schedule
*  07/2005 Written by Tom Doan, Estima
*
proc ADTest x start end
type series x
type integer start end
*
option real mean
option real variance
*
local integer startl endl i
local series ox fx
local real meanl stddev
local real z z0 k kpr k2pr v w
local real target weight
local real ad ks gap
*
inquire(series=x) startl>>start endl>>end
*
* ox is a copy of x which we sort
*
set ox startl endl = x
order ox startl endl
stats(noprint) ox startl endl
compute meanl=%if(%defined(mean),mean,%mean)
compute stddev=sqrt(%if(%defined(variance),variance,%variance))
*
* Get the CDF's of the standardized values
*
set fx startl endl = %cdf((ox-meanl)/stddev)
*
* Compute the Anderson-Darling and Kolmogorov-Smirnov statistics
* based upon a comparison of the hypothesized and empirical CDF's
*
compute ad=0.0,ks=0.0
do i=startl,endl
   compute ad =ad+(2.0*(i-startl+1)-1)/%nobs*(log(fx(i))+log(1-fx(endl+startl-i)))
   compute gap=abs(fx(i)-(i-startl+1.0)/%nobs)
   compute ks =%max(ks,gap)
end do i
compute ad=-%nobs-ad
*
*  Evaluating the saddle point approximation for the critical value of the
*  Anderson-Darling test. Use Newton's method to find the saddlepoint, which
*  is the point where the derivative of the cumulant generating function is
*  equal to the observed test statistic. The CGF of the AD test is an infinite
*  sum, so we truncate it at 500 terms.
*
compute z=0.0
compute target=ad
loop
   compute k=kpr=k2pr=0.0
   do i=1,500
      compute weight=1.0/(i*(i+1.0))
      compute k   =k    -.5  *  log(1-2*z*weight)
      compute kpr =kpr  +   weight/(1-2*z*weight)
      compute k2pr=k2pr +2*(weight/(1-2*z*weight))**2
   end do i
   *
   *  If step was too big, reduce it and try again
   *
   if .not.%valid(k) {
      compute z=z+.5*(z0-z)
      next
   }
   if kpr==target
      break
   compute z0=z
   compute z=z-(kpr-target)/k2pr
   if abs(z0-z)<1.e-8
      break
end loop
*

compute w=sqrt(2*(z*target-k))*%sign(z)
compute v=z*sqrt(k2pr)
compute %signif=1-%cdf(w)+%density(w)*(1.0/v-1.0/w)
*
disp "Anderson-Darling Test for Normality"
disp "Series" %l(x) "Observations" %nobs
disp "AD Statistic" ad
disp "Approximate P-value" #.##### %signif
compute %cdstat=ad
end


