*
* @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]
*
* Variables Defined:
*   %CDSTAT    test statistic
*   %SIGNIF    approximate significance level
*
* 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
*  07/2012 Switch to use REPORT. Add TITLE and PRINT options
*  05/2015 Fail-safe against log's of 0 and 1.
*
proc ADTest x start end
type series x
type integer start end
*
option real    mean
option real    variance
option string  title
option switch  print  	1
*
local integer  startl endl i
local series   ox fx fxc
local real     meanl stddev z
local real     z z0 k kpr k2pr v w
local real     target weight
local real     ad ks gap
local report   areport
local string   ltitle
*
inquire(series=x) startl<<start endl<<end
*
* ox is a copy of x which we sort
*
clear ox
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. Because we need to take
* logs, we use the left tails, which are computed more accurately than
* the right tails.
*
set fx  startl endl  = z=(ox-meanl)/stddev,%if(z<0,%cdf(z),1-%cdf(-z))
set fxc startl endl  = z=(ox-meanl)/stddev,%if(z<0,1-%cdf(z),%cdf(-z))
*
* 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.0)-1)/%nobs*(log(fx(i))+log(fxc(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 %cdstat=ad
compute %signif=1-%cdf(w)+%density(w)*(1.0/v-1.0/w)
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Anderson-Darling Test for Normality"
*
report(use=areport,action=define,title=ltitle)
report(use=areport,atrow=1,atcol=1,span) ltitle
report(use=areport,atrow=2,atcol=1) "Series" %l(x)
report(use=areport,atrow=3,atcol=1) "Observations" %nobs
report(use=areport,atrow=4,atcol=1) "AD Statistic" ad
report(use=areport,atrow=5,atcol=1) "Signif Level" %signif
report(use=areport,atrow=4,action=format,picture="*.####")
if print
   report(use=areport,action=show)
end
