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