Durbin algorithm works correctly

Author: kagalenko-m-b

src/signalcorr.jl

@@ -530,7 +530,7 @@ crosscor(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) =
 #
 #######################################
 
-function pacf_regress!(r::RealMatrix, X::RealMatrix, lags::IntegerVector, mk::Integer)
+function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector, mk::Integer)
     lx = size(X, 1)
     tmpX = ones(eltype(X), lx, mk + 1)
     for j = 1 : size(X,2)
@@ -548,8 +548,8 @@ function pacf_regress!(r::RealMatrix, X::RealMatrix, lags::IntegerVector, mk::In
     r
 end
 
-function pacf_yulewalker!(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector, mk::Integer) where T<:RealFP
-    tmp = Vector{T}(undef, mk)
+function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector, mk::Integer)
+    tmp = Vector{eltype(X)}(undef, mk)
     for j = 1 : size(X,2)
         acfs = autocor(X[:,j], 1:mk)
         for i = 1 : length(lags)
@@ -571,7 +571,7 @@ using the Yule-Walker equations.
 
 `r` must be a matrix of size `(length(lags), size(x, 2))`.
 """
-function pacf!(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector; method::Symbol=:regression) where T<:RealFP
+function pacf!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector; method::Symbol=:regression)
     lx = size(X, 1)
     m = length(lags)
     minlag, maxlag = extrema(lags)
@@ -602,11 +602,11 @@ If `x` is a vector, return a vector of the same length as `lags`.
 If `x` is a matrix, return a matrix of size `(length(lags), size(x, 2))`,
 where each column in the result corresponds to a column in `x`.
 """
-function pacf(X::RealMatrix, lags::IntegerVector; method::Symbol=:regression)
+function pacf(X::AbstractMatrix, lags::IntegerVector; method::Symbol=:regression)
     out = Matrix{float(eltype(X))}(undef, length(lags), size(X,2))
     pacf!(out, float(X), lags; method=method)
 end
 
-function pacf(x::RealVector, lags::IntegerVector; method::Symbol=:regression)
+function pacf(x::AbstractVector, lags::IntegerVector; method::Symbol=:regression)
     vec(pacf(reshape(x, length(x), 1), lags, method=method))
 end

src/toeplitzsolvers.jl

@@ -1,28 +1,28 @@
 # Symmetric Toeplitz solver
-function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T<:BlasReal
+function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T
     n = length(r)
     n <= length(y) || throw(DimensionMismatch("Auxiliary vector cannot be shorter than data vector"))
     y[1] = -r[1]
     β = one(T)
     α = -r[1]
     for k = 1:n-1
-        β *= one(T) - α*α
+        β *= one(T) - abs2(α)
         α = -r[k+1]
         for j = 1:k
             α -= r[k-j+1]*y[j]
         end
         α /= β
         for j = 1:div(k,2)
             tmp = y[j]
-            y[j] += α*y[k-j+1]
-            y[k-j+1] += α*tmp
+            y[j] += α*conj(y[k-j+1])
+            y[k-j+1] += α*conj(tmp)
         end
-        if isodd(k) y[div(k,2)+1] *= one(T) + α end
+        if isodd(k) y[div(k,2)+1] += α*conj(y[div(k,2)+1]) end
         y[k+1] = α
     end
     return y
 end
-durbin(r::AbstractVector{T}) where {T<:BlasReal} = durbin!(r, zeros(T, length(r)))
+durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r)))
 
 function levinson!(r::AbstractVector{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T<:BlasReal
     n = length(b)

test/signalcorr.jl

@@ -104,34 +104,17 @@ c22 = crosscor(x2, x2)
 
 
 ## pacf
-x_pacf = [-2.133252557240862
-          .1775816414485478
-          -.6264517920318317
-          -.8809042583216906
-          .09251017186697393
-          -.9271887119115569
-          3.355819743178915
-          -.2834039258495755
-          .5354280026977677
-          .39182285417742585]
-
-rpacfr = [-0.2181581223814188
-          0.1950153168287112
-          0.1443158046061393
-          -0.1997912294497786
-          1.686452299291934
-          -1.446624527965185]
-
-@test pacf(real.(x_pacf), 1:4) ≈ rpacfr[1:4]
-
-rpacfy = [-0.2211730116688735
-          0.1896833143080208
-          0.1118570207337193
-          -0.1750206698354201
-          0.04490978543108597
-          -0.3715806391613901
-          -0.1743836518336153
-          0.07574637899230210
-          0.04970724340234131]
-
-@test pacf(real.(x_pacf), 1:4, method=:yulewalker) ≈ rpacfy[1:4]
+
+pacfr = [-1.598495044296996e-03 - 2.915104118351207e-01im
+         -5.560162016912027e-01 + 2.950837739894279e-01im
+         -2.547001916363494e-02 + 2.326084658014266e-01im
+         -5.427443903358727e-01 + 3.146715147305132e-01im];
+
+@test pacf(x1, 1:4) ≈ pacfr[1:4]
+
+pacfy = [-1.598495044296996e-03 - 2.915104118351207e-01im
+         -5.560162016912027e-01 + 2.950837739894279e-01im
+         -2.547001916363494e-02 + 2.326084658014266e-01im
+         -5.427443903358727e-01 + 3.146715147305132e-01im];
+
+@test pacf(x1, 1:4, method=:yulewalker) ≈ pacfy