The RATS Software Forum

Code: Select all

*
* @DurbinLevinson( options)  series  start  end
*
* Uses the Durbin-Levinson recursion to estimate the coefficients in a sequence of
* autoregressive representations for a stationary series.
*
* Parameters:
*   series     series to analyze
*   start end  range of <<series>> to use. By default, the defined range of <<series>>
*
* Options:
*   M=number of autoregressive coefficients
*   METHOD=[YULE]/BURG
*     Sets the method used for computing the covariances. If you use METHOD=BURG, you'll
*     get the sequence of Burg AR estimates; if METHOD=YULE (default), they'll be the
*     Yule-Walker AR estimates.
*   COVARIANCES=input series of covariances. If you don't use this, the covariances
*     are computed using the input <<X>> series. If you do use this, the <<X>> series
*     is ignored.
*   PHI=(output) m x m matrix of autoregressive representations. Row i has the AR(i)
*   V=(output) m vector of estimated residual variances. Row i has the residual
*     variance for the AR(i)
*
* Defines:
*   %BETA=AR coefficients in the final AR
*   %SIGMASQ=Estimated residual variance for the final AR
*
procedure DurbinLevinson series start end
type series  series
type integer start end
*
option integer m      13
option choice  method 1   yule burg
option series  covar
option rect    *phi
option vect    *v
*
local rect     phil
local vect     vl
local real     s
local series   covarx
local integer  j n
*
local integer startl endl
*
* If the <<covar>> option is used, just copy the covariances out of it.
* Otherwise, compute the covariances from the input series.
*
if %defined(covar)
   set covarx 1 m+1 = covar
else {
   if .not.%defined(series) {
      disp "Syntax: @DurbinLevinson(options)  SERIES start end"
      return
   }
   inquire(series=series) startl<<start endl<<end
   corr(covariances,number=m+1,method=method,noprint) series startl endl covarx
}
dim phil(m,m) vl(m+1)
compute phil=%zeros(m,m)
compute vl(1)=covarx(1)
*
do n=1,m
   compute s=covarx(n+1)
   do j=1,n-1
      compute s=s-phil(n-1,j)*covarx(n+1-j)
   end do j
   compute phil(n,n)=s/vl(n)
   do j=1,n-1
      compute phil(n,j)=phil(n-1,j)-phil(n,n)*phil(n-1,n-j)
   end do j
   compute vl(n+1)=vl(n)*(1-phil(n,n)**2)
end do n
if %defined(phi)
   compute phi=phil
if %defined(v)
   compute v=%xsubvec(vl,2,m+1)
*
compute %beta   =%xsubmat(phil,m,m,1,m)
compute %sigmasq=vl(m+1)
*
end