diff --git a/src/signalcorr.jl b/src/signalcorr.jl
index 06c83ba1d..99dd0385f 100644
--- a/src/signalcorr.jl
+++ b/src/signalcorr.jl
@@ -14,7 +14,7 @@
 default_laglen(lx::Int) = min(lx-1, round(Int,10*log10(lx)))
 check_lags(lx::Int, lags::AbstractVector) = (maximum(lags) < lx || error("lags must be less than the sample length."))
 
-function demean_col!(z::AbstractVector{<:Real}, x::AbstractMatrix{<:Real}, j::Int, demean::Bool)
+function demean_col!(z::AbstractVector, x::AbstractMatrix, j::Int, demean::Bool)
     T = eltype(z)
     m = size(x, 1)
     @assert m == length(z)
@@ -44,7 +44,7 @@ end
 default_autolags(lx::Int) = 0 : default_laglen(lx)
 
 _autodot(x::AbstractVector{<:Union{Float32, Float64}}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx)
-_autodot(x::AbstractVector{<:Real}, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx))
+_autodot(x::AbstractVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx))
 
 
 ## autocov
@@ -61,7 +61,7 @@ where each column in the result will correspond to a column in `x`.
 
 The output is not normalized. See [`autocor!`](@ref) for a method with normalization.
 """
-function autocov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocov!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     m = length(lags)
     length(r) == m || throw(DimensionMismatch())
@@ -75,7 +75,7 @@ function autocov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::Ab
     return r
 end
 
-function autocov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocov!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     ns = size(x, 2)
     m = length(lags)
@@ -110,17 +110,17 @@ When left unspecified, the lags used are the integers from 0 to
 
 The output is not normalized. See [`autocor`](@ref) for a function with normalization.
 """
-function autocov(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocov(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(eltype(x))}(undef, length(lags))
     autocov!(out, x, lags; demean=demean)
 end
 
-function autocov(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocov(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2))
     autocov!(out, x, lags; demean=demean)
 end
 
-autocov(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
+autocov(x::AbstractVecOrMat; demean::Bool=true) =
     autocov(x, default_autolags(size(x,1)); demean=demean)
 
 ## autocor
@@ -139,7 +139,7 @@ where each column in the result will correspond to a column in `x`.
 The output is normalized by the variance of `x`, i.e. so that the lag 0
 autocorrelation is 1. See [`autocov!`](@ref) for the unnormalized form.
 """
-function autocor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocor!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     m = length(lags)
     length(r) == m || throw(DimensionMismatch())
@@ -154,7 +154,7 @@ function autocor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::Ab
     return r
 end
 
-function autocor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocor!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     ns = size(x, 2)
     m = length(lags)
@@ -191,17 +191,17 @@ When left unspecified, the lags used are the integers from 0 to
 The output is normalized by the variance of `x`, i.e. so that the lag 0
 autocorrelation is 1. See [`autocov`](@ref) for the unnormalized form.
 """
-function autocor(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocor(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(eltype(x))}(undef, length(lags))
     autocor!(out, x, lags; demean=demean)
 end
 
-function autocor(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function autocor(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2))
     autocor!(out, x, lags; demean=demean)
 end
 
-autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
+autocor(x::AbstractVecOrMat; demean::Bool=true) =
     autocor(x, default_autolags(size(x,1)); demean=demean)
 
 
@@ -213,14 +213,14 @@ autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
 
 default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l)
 
-function _crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:Union{Float32, Float64}}
+function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int)
     if l >= 0
         dot(x, 1:(lx-l), y, (1+l):lx)
     else
         dot(x, (1-l):lx, y, 1:(lx+l))
     end
 end
-function _crossdot(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lx::Int, l::Int)
+function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int)
     if l >= 0
         dot(view(x, 1:(lx-l)), view(y, (1+l):lx))
     else
@@ -233,7 +233,7 @@ end
 """
     crosscov!(r, x, y, lags; demean=true)
 
-Compute the cross covariance function (CCF) between real-valued vectors or matrices
+Compute the cross covariance function (CCF) between vectors or matrices
 `x` and `y` at `lags` and store the result in `r`. `demean` specifies whether the
 respective means of `x` and `y` should be subtracted from them before computing their
 CCF.
@@ -246,7 +246,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`.
 
 The output is not normalized. See [`crosscor!`](@ref) for a function with normalization.
 """
-function crosscov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov!(r::AbstractVector, x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     m = length(lags)
     (length(y) == lx && length(r) == m) || throw(DimensionMismatch())
@@ -262,7 +262,7 @@ function crosscov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::Abst
     return r
 end
 
-function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov!(r::AbstractMatrix, x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     ns = size(x, 2)
     m = length(lags)
@@ -282,7 +282,7 @@ function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::Abst
     return r
 end
 
-function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     ns = size(y, 2)
     m = length(lags)
@@ -302,7 +302,7 @@ function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::Abst
     return r
 end
 
-function crosscov!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov!(r::AbstractArray{<:Any,3}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     nx = size(x, 2)
     ny = size(y, 2)
@@ -343,7 +343,7 @@ end
 """
     crosscov(x, y, [lags]; demean=true)
 
-Compute the cross covariance function (CCF) between real-valued vectors or
+Compute the cross covariance function (CCF) between vectors or
 matrices `x` and `y`, optionally specifying the `lags`. `demean` specifies
 whether the respective means of `x` and `y` should be subtracted from them
 before computing their CCF.
@@ -356,27 +356,27 @@ When left unspecified, the lags used are the integers from
 
 The output is not normalized. See [`crosscor`](@ref) for a function with normalization.
 """
-function crosscov(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov(x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-function crosscov(x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov(x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-function crosscov(x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov(x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-function crosscov(x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscov(x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-crosscov(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=true) =
+crosscov(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) =
     crosscov(x, y, default_crosslags(size(x,1)); demean=demean)
 
 
@@ -384,7 +384,7 @@ crosscov(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=
 """
     crosscor!(r, x, y, lags; demean=true)
 
-Compute the cross correlation between real-valued vectors or matrices `x` and `y` at
+Compute the cross correlation between vectors or matrices `x` and `y` at
 `lags` and store the result in `r`. `demean` specifies whether the respective means of
 `x` and `y` should be subtracted from them before computing their cross correlation.
 
@@ -397,7 +397,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`.
 The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov!`](@ref) for the
 unnormalized form.
 """
-function crosscor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor!(r::AbstractVector, x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     m = length(lags)
     (length(y) == lx && length(r) == m) || throw(DimensionMismatch())
@@ -414,7 +414,7 @@ function crosscor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::Abst
     return r
 end
 
-function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor!(r::AbstractMatrix, x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     ns = size(x, 2)
     m = length(lags)
@@ -436,7 +436,7 @@ function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::Abst
     return r
 end
 
-function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     ns = size(y, 2)
     m = length(lags)
@@ -458,7 +458,7 @@ function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::Abst
     return r
 end
 
-function crosscor!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor!(r::AbstractArray{<:Any, 3}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     nx = size(x, 2)
     ny = size(y, 2)
@@ -504,7 +504,7 @@ end
 """
     crosscor(x, y, [lags]; demean=true)
 
-Compute the cross correlation between real-valued vectors or matrices `x` and `y`,
+Compute the cross correlation between vectors or matrices `x` and `y`,
 optionally specifying the `lags`. `demean` specifies whether the respective means of
 `x` and `y` should be subtracted from them before computing their cross correlation.
 
@@ -517,27 +517,27 @@ When left unspecified, the lags used are the integers from
 The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov`](@ref) for the
 unnormalized form.
 """
-function crosscor(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor(x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-function crosscor(x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor(x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-function crosscor(x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor(x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-function crosscor(x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
+function crosscor(x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-crosscor(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=true) =
+crosscor(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) =
     crosscor(x, y, default_crosslags(size(x,1)); demean=demean)
 
 
@@ -549,7 +549,7 @@ crosscor(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=
 #
 #######################################
 
-function pacf_regress!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}, mk::Integer)
+function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer)
     lx = size(X, 1)
     tmpX = ones(eltype(X), lx, mk + 1)
     for j = 1 : size(X,2)
@@ -561,19 +561,22 @@ function pacf_regress!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}, lag
         for i = 1 : length(lags)
             l = lags[i]
             sX = view(tmpX, 1+l:lx, 1:l+1)
-            r[i,j] = l == 0 ? 1 : (cholesky!(sX'sX, Val(false)) \ (sX'view(X, 1+l:lx, j)))[end]
+            r[i,j] = l == 0 ? 1 : (cholesky!(sX'sX) \ (sX'view(X, 1+l:lx, j)))[end]
         end
     end
     r
 end
 
-function pacf_yulewalker!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}, mk::Integer) where T<:Union{Float32, Float64}
-    tmp = Vector{T}(undef, mk)
+function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer)
+    T = eltype(X)
+    p = Vector{T}(undef, mk)
+    y = Vector{T}(undef, mk)
     for j = 1 : size(X,2)
         acfs = autocor(X[:,j], 1:mk)
+        durbin!(acfs, y, p)
         for i = 1 : length(lags)
             l = lags[i]
-            r[i,j] = l == 0 ? 1 : l == 1 ? acfs[i] : -durbin!(view(acfs, 1:l), tmp)[l]
+            r[i,j] = l == 0 ? 1 : -p[l]
         end
     end
 end
@@ -590,7 +593,7 @@ using the Yule-Walker equations.
 
 `r` must be a matrix of size `(length(lags), size(x, 2))`.
 """
-function pacf!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) where T<:Union{Float32, Float64}
+function pacf!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
     lx = size(X, 1)
     m = length(lags)
     minlag, maxlag = extrema(lags)
@@ -611,7 +614,7 @@ end
 """
     pacf(X, lags; method=:regression)
 
-Compute the partial autocorrelation function (PACF) of a real-valued vector
+Compute the partial autocorrelation function (PACF) of a vector
 or matrix `X` at `lags`. `method` designates the estimation method. Recognized
 values are `:regression`, which computes the partial autocorrelations via successive
 regression models, and `:yulewalker`, which computes the partial autocorrelations
@@ -621,11 +624,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::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
+function pacf(X::AbstractMatrix, lags::AbstractVector{<:Integer}; 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::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
+function pacf(x::AbstractVector, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
     vec(pacf(reshape(x, length(x), 1), lags, method=method))
 end
diff --git a/src/toeplitzsolvers.jl b/src/toeplitzsolvers.jl
index 19146bf7a..16a905a6c 100644
--- a/src/toeplitzsolvers.jl
+++ b/src/toeplitzsolvers.jl
@@ -1,12 +1,34 @@
-# Symmetric Toeplitz solver
-function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T<:BlasReal
+"""
+    durbin!(r::AbstractVector, y::AbstractVector, p::AbstractVector) -> Vector
+
+Solve Yule-Walker equations using Durbin algorithm.
+
+Solution of N×N system is obtained iteratively, by solving 1×1, 2×2, ...
+in succession. For use in computing partial autocorrelation matrix,
+return the last elements of the successive solutions in vector `p`.
+
+Reference: Golub, G. H., and C. F. Van Loan. "Matrix computations 4th edition the johns hopkins university press." Baltimore, MD (2013), section 4.7
+
+# Arguments
+- `r::AbstractVector`: first column of a Herimitian positive definite
+    Toeplitz matrix, excluding the diagonal element (equal to one).
+- `y::AbstractVector`: work vector for solution, should have the same length
+    as `r`
+- `p::AbstractVector`: last elements of the successive solutions, should
+     have the same length  as `r`
+
+# Returns
+- `y::AbstractVector`: return the solution vector (same as the second argument)
+"""
+function durbin!(r::AbstractVector{T}, y::AbstractVector{T}, p::AbstractVector{T}) where {T}
     n = length(r)
-    n <= length(y) || throw(DimensionMismatch("Auxiliary vector cannot be shorter than data vector"))
+    n <= length(p) || n <= length(y) || throw(DimensionMismatch("Auxiliary vectors cannot be shorter than data vector"))
     y[1] = -r[1]
+    p[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]
@@ -14,15 +36,18 @@ function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T<:BlasReal
         α /= β
         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] += α*conj(y[div(k,2)+1]) 
         end
-        if isodd(k) y[div(k,2)+1] *= one(T) + α end
         y[k+1] = α
+        p[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)), zeros(T, length(r)))
 
 function levinson!(r::AbstractVector{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T<:BlasReal
     n = length(b)
diff --git a/test/runtests.jl b/test/runtests.jl
index 94986d240..1caa2269a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,7 @@ tests = ["ambiguous",
          "hist",
          "rankcorr",
          "signalcorr",
+         "signalcorr_real",
          "misc",
          "pairwise",
          "reliability",
diff --git a/test/signalcorr.jl b/test/signalcorr.jl
index bfbe90fed..22da15f3e 100644
--- a/test/signalcorr.jl
+++ b/test/signalcorr.jl
@@ -9,88 +9,89 @@ using Test
 
 # random data for testing
 
-x = [-2.133252557240862    -.7445937365828654;
-       .1775816414485478   -.5834801838041446;
-      -.6264517920318317   -.68444205333293;
-      -.8809042583216906    .9071671734302398;
-       .09251017186697393 -1.0404476733379926;
-      -.9271887119115569   -.620728578941385;
-      3.355819743178915    -.8325051361909978;
-      -.2834039258495755   -.22394811874731657;
-       .5354280026977677    .7481337671592626;
-       .39182285417742585   .3085762550821047]
+x = [0.5316732188954801 + 0.015032026010878084im -0.8051436916816284 - 0.730238071058262im
+     -0.26236085703683865 + 0.1561582684776427im -0.7733975759329091 + 0.7517417038318899im
+     -0.12750317530769195 + 0.8359780109572358im -0.1426497126379726 - 0.1976368614625515im
+     -0.3779012097765515 - 0.15588720329528907im  0.4625694946681199 - 0.6687893051682382im
+     -0.14198553599197566 - 0.911613885474317im  -0.38003834312962825 - 0.16244828932805647im
+     0.6809910683979451 - 1.1951424604960583im    0.36127363787490463 + 1.3578919531660747im
+     -0.7024909626802256 - 0.2915841698461933im   0.8432078641246675 - 1.373657703019187im
+     -0.5684618031428246 + 0.3729428437838556im  -0.18853209977543475 - 0.9698921354142658im
+     1.2992421672413557 - 0.17943637636017656im  -1.7801079183259816 - 1.3063359194617465im
+     -0.26481424468298875 - 1.1636864449113618im -0.2605299632484854 - 0.7356290016171836im]
 
 x1 = view(x, :, 1)
 x2 = view(x, :, 2)
-realx = convert(AbstractMatrix{Real}, x)
-realx1 = convert(AbstractVector{Real}, x1)
-realx2 = convert(AbstractVector{Real}, x2)
-
+cmplx_x = convert(AbstractMatrix{Complex}, x)
+cmplx_x1 = convert(AbstractVector{Complex}, x1)
+cmplx_x2 = convert(AbstractVector{Complex}, x2)
 # autocov & autocorr
 
 @test autocov([1:5;])          ≈ [2.0, 0.8, -0.2, -0.8, -0.8]
 @test autocor([1, 2, 3, 4, 5]) ≈ [1.0, 0.4, -0.1, -0.4, -0.4]
-
-racovx1 =  [1.839214242630635709475,
-           -0.406784553146903871124,
-            0.421772254824993531042,
-            0.035874943792884653182,
-           -0.255679775928512320604,
-            0.231154400105831353551,
-           -0.787016960267425180753,
-            0.039909287349160660341,
-           -0.110149697877911914579,
-           -0.088687020167434751916]
-
-@test autocov(x1) ≈ racovx1
-@test autocov(realx1) ≈ racovx1
+@test_throws MethodError autocov([1 missing 2 3 4 5])
+@test_throws MethodError autocor([1 missing 2 3 4 5])
+
+acovx1 =  [0.755284179631112 + 0.0im,
+           -0.005333112170365584 - 0.18633291805458002im,
+           -0.28638133506011604 + 0.1397225175604478im,
+           -0.05476759020476188 + 0.12991227531087618im,
+           -0.00499902418309206 + 0.023463421186162473im,
+           0.11854989361702409 - 0.0028772373657594066im,
+           -0.005123812003979093 - 0.08041878599400383im,
+           -0.14090692583679154 + 0.035230042290069444im,
+           0.039899180711674635 + 0.004918236900383659im,
+           -0.03857936468514856 - 0.04063999068714769im]
+
+@test autocov(x1) ≈ acovx1
+@test autocov(cmplx_x1) ≈ acovx1
 @test autocov(x)  ≈ [autocov(x1) autocov(x2)]
-@test autocov(realx)  ≈ [autocov(realx1) autocov(realx2)]
-
-racorx1 = [0.999999999999999888978,
-          -0.221173011668873431557,
-           0.229321981664153962122,
-           0.019505581764945757045,
-          -0.139015765538446717242,
-           0.125681062460244019618,
-          -0.427909344123907742219,
-           0.021699096507690283225,
-          -0.059889541590524189574,
-          -0.048220059475281865091]
-
-@test autocor(x1) ≈ racorx1
-@test autocor(realx1) ≈ racorx1
+@test autocov(cmplx_x)  ≈ [autocov(cmplx_x1) autocov(cmplx_x2)]
+
+acorx1 = [1.0 + 0.0im,
+          -0.007061066965510023 - 0.24670570770539213im,
+          -0.3791703080554228 + 0.18499330626611243im,
+          -0.07251256107537019 + 0.17200449686941222im,
+          -0.006618732813301651 + 0.031065686027770673im,
+          0.1569606471499575 - 0.0038094765432061303im,
+          -0.00678395250709688 - 0.106474871528861im,
+          -0.18656146869859214 + 0.04664475073114351im,
+          0.05282671316002114 + 0.0065117700502952056im,
+          -0.051079270194684986 - 0.0538075492419246im]
+
+@test autocor(x1) ≈ acorx1
+@test autocor(cmplx_x1) ≈ acorx1
 @test autocor(x)  ≈ [autocor(x1) autocor(x2)]
-@test autocor(realx)  ≈ [autocor(realx1) autocor(realx2)]
+@test autocor(cmplx_x)  ≈ [autocor(cmplx_x1) autocor(cmplx_x2)]
 
 
 # crosscov & crosscor
 
-rcov0 = [0.320000000000000006661,
-        -0.319999999999999951150,
-         0.080000000000000029421,
-        -0.479999999999999982236,
-         0.000000000000000000000,
-         0.479999999999999982236,
-        -0.080000000000000029421,
-         0.319999999999999951150,
-        -0.320000000000000006661]
+rcov0 = [0.96 + 0.32im,
+         -0.96 - 0.32im,
+         0.24 + 0.08im,
+         -1.44 - 0.48im,
+         0.0 + 0.0im,
+         1.44 + 0.48im,
+         -0.24 - 0.08im,
+         0.96 + 0.32im,
+         -0.96 - 0.32im]
 
-@test crosscov([1, 2, 3, 4, 5], [1, -1, 1, -1, 1]) ≈ rcov0
-@test crosscov([1:5;], [1:5;]) ≈ [-0.8, -0.8, -0.2, 0.8, 2.0, 0.8, -0.2, -0.8, -0.8]
+@test crosscov([1+9im, 2+7im, 3+5im, 4+3im, 5+1im], [1-im, -1+im, 1-im, -1+im, 1-im]) ≈ rcov0
+@test crosscov([1:5 ;], [1:5;]) ≈ [-0.8, -0.8, -0.2, 0.8, 2.0, 0.8, -0.2, -0.8, -0.8]
 
 c11 = crosscov(x1, x1)
 c12 = crosscov(x1, x2)
 c21 = crosscov(x2, x1)
 c22 = crosscov(x2, x2)
-@test crosscov(realx1, realx2) ≈ c12
+@test crosscov(cmplx_x1, cmplx_x2) ≈ c12
 
 @test crosscov(x,  x1) ≈ [c11 c21]
-@test crosscov(realx, realx1) ≈ [c11 c21]
+@test crosscov(cmplx_x, cmplx_x1) ≈ [c11 c21]
 @test crosscov(x1, x)  ≈ [c11 c12]
-@test crosscov(realx1, realx)  ≈ [c11 c12]
+@test crosscov(cmplx_x1, cmplx_x)  ≈ [c11 c12]
 @test crosscov(x,  x)  ≈ cat([c11 c21], [c12 c22], dims=3)
-@test crosscov(realx,  realx)  ≈ cat([c11 c21], [c12 c22], dims=3)
+@test crosscov(cmplx_x,  cmplx_x)  ≈ cat([c11 c21], [c12 c22], dims=3)
 
 # issue #805: avoid converting one input to the other's eltype
 @test crosscov([34566.5345, 3466.4566], Float16[1, 10]) ≈
@@ -114,14 +115,14 @@ c11 = crosscor(x1, x1)
 c12 = crosscor(x1, x2)
 c21 = crosscor(x2, x1)
 c22 = crosscor(x2, x2)
-@test crosscor(realx1, realx2) ≈ c12
+@test crosscor(cmplx_x1, cmplx_x2) ≈ c12
 
 @test crosscor(x,  x1) ≈ [c11 c21]
-@test crosscor(realx, realx1) ≈ [c11 c21]
+@test crosscor(cmplx_x, cmplx_x1) ≈ [c11 c21]
 @test crosscor(x1, x)  ≈ [c11 c12]
-@test crosscor(realx1, realx)  ≈ [c11 c12]
+@test crosscor(cmplx_x1, cmplx_x)  ≈ [c11 c12]
 @test crosscor(x,  x)  ≈ cat([c11 c21], [c12 c22], dims=3)
-@test crosscor(realx, realx)  ≈ cat([c11 c21], [c12 c22], dims=3)
+@test crosscor(cmplx_x, cmplx_x)  ≈ cat([c11 c21], [c12 c22], dims=3)
 
 # issue #805: avoid converting one input to the other's eltype
 @test crosscor([34566.5345, 3466.4566], Float16[1, 10]) ≈
@@ -129,18 +130,36 @@ c22 = crosscor(x2, x2)
     crosscor([34566.5345, 3466.4566], Float16[1, 10])
 
 
-## pacf
-
-rpacfr =  [-0.218158122381419,
-            0.195015316828711,
-            0.144315804606139,
-           -0.199791229449779]
-
-@test pacf(x[:,1], 1:4) ≈ rpacfr
-
-rpacfy = [-0.221173011668873,
-           0.189683314308021,
-           0.111857020733719,
-          -0.175020669835420]
-
-@test pacf(x[:,1], 1:4, method=:yulewalker) ≈ rpacfy
+## pacf least squares
+pacf_ls = [-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) ≈ pacf_ls[1:4]
+
+## pacf Yule-Walker
+
+function yulewalker_qr(v::AbstractVector)
+    A = toeplitz(v)
+    b = v[2:end]
+    x = -A\b
+end
+function toeplitz(v::AbstractVector{T}) where T
+    N = length(v)
+    A = zeros(T, N - 1, N - 1)
+    for n in 1:N-1
+        A[n, n+1:end] = conj(v[2:N-n])
+        A[n, 1:n] = reverse(v[1:n])
+    end
+    return A
+end
+# durbin solver 
+acf = autocor(x1)
+p_qr = [yulewalker_qr(acf[1:n])[n-1] for n in 2:length(acf)]
+y = similar(acf[2:end])
+pd = similar(acf[2:end])
+StatsBase.durbin!(acf[2:end], y, pd)
+@test p_qr ≈ pd
+
+@test pacf(x1, 1:4, method=:yulewalker) ≈ -p_qr[1:4]
diff --git a/test/signalcorr_real.jl b/test/signalcorr_real.jl
new file mode 100644
index 000000000..bce1c83a3
--- /dev/null
+++ b/test/signalcorr_real.jl
@@ -0,0 +1,136 @@
+# Test corr.jl
+#
+# Many test cases are migrated from test/01.jl in the old version
+# The reference results are generated from R.
+#
+
+using StatsBase
+using Test
+
+# random data for testing
+
+x = [-2.133252557240862    -.7445937365828654;
+       .1775816414485478   -.5834801838041446;
+      -.6264517920318317   -.68444205333293;
+      -.8809042583216906    .9071671734302398;
+       .09251017186697393 -1.0404476733379926;
+      -.9271887119115569   -.620728578941385;
+      3.355819743178915    -.8325051361909978;
+      -.2834039258495755   -.22394811874731657;
+       .5354280026977677    .7481337671592626;
+       .39182285417742585   .3085762550821047]
+
+x1 = view(x, :, 1)
+x2 = view(x, :, 2)
+realx = convert(AbstractMatrix{Real}, x)
+realx1 = convert(AbstractVector{Real}, x1)
+realx2 = convert(AbstractVector{Real}, x2)
+
+# autocov & autocorr
+
+@test autocov([1:5;])          ≈ [2.0, 0.8, -0.2, -0.8, -0.8]
+@test autocor([1, 2, 3, 4, 5]) ≈ [1.0, 0.4, -0.1, -0.4, -0.4]
+
+racovx1 =  [1.839214242630635709475,
+           -0.406784553146903871124,
+            0.421772254824993531042,
+            0.035874943792884653182,
+           -0.255679775928512320604,
+            0.231154400105831353551,
+           -0.787016960267425180753,
+            0.039909287349160660341,
+           -0.110149697877911914579,
+           -0.088687020167434751916]
+
+@test autocov(x1) ≈ racovx1
+@test autocov(realx1) ≈ racovx1
+@test autocov(x)  ≈ [autocov(x1) autocov(x2)]
+@test autocov(realx)  ≈ [autocov(realx1) autocov(realx2)]
+
+racorx1 = [0.999999999999999888978,
+          -0.221173011668873431557,
+           0.229321981664153962122,
+           0.019505581764945757045,
+          -0.139015765538446717242,
+           0.125681062460244019618,
+          -0.427909344123907742219,
+           0.021699096507690283225,
+          -0.059889541590524189574,
+          -0.048220059475281865091]
+
+@test autocor(x1) ≈ racorx1
+@test autocor(realx1) ≈ racorx1
+@test autocor(x)  ≈ [autocor(x1) autocor(x2)]
+@test autocor(realx)  ≈ [autocor(realx1) autocor(realx2)]
+
+
+# crosscov & crosscor
+
+rcov0 = [0.320000000000000006661,
+        -0.319999999999999951150,
+         0.080000000000000029421,
+        -0.479999999999999982236,
+         0.000000000000000000000,
+         0.479999999999999982236,
+        -0.080000000000000029421,
+         0.319999999999999951150,
+        -0.320000000000000006661]
+
+@test crosscov([1, 2, 3, 4, 5], [1, -1, 1, -1, 1]) ≈ rcov0
+@test crosscov([1:5;], [1:5;]) ≈ [-0.8, -0.8, -0.2, 0.8, 2.0, 0.8, -0.2, -0.8, -0.8]
+
+c11 = crosscov(x1, x1)
+c12 = crosscov(x1, x2)
+c21 = crosscov(x2, x1)
+c22 = crosscov(x2, x2)
+@test crosscov(realx1, realx2) ≈ c12
+
+@test crosscov(x,  x1) ≈ [c11 c21]
+@test crosscov(realx, realx1) ≈ [c11 c21]
+@test crosscov(x1, x)  ≈ [c11 c12]
+@test crosscov(realx1, realx)  ≈ [c11 c12]
+@test crosscov(x,  x)  ≈ cat([c11 c21], [c12 c22], dims=3)
+@test crosscov(realx,  realx)  ≈ cat([c11 c21], [c12 c22], dims=3)
+
+rcor0 = [0.230940107675850,
+        -0.230940107675850,
+         0.057735026918963,
+        -0.346410161513775,
+         0.000000000000000,
+         0.346410161513775,
+        -0.057735026918963,
+         0.230940107675850,
+        -0.230940107675850]
+
+@test crosscor([1, 2, 3, 4, 5], [1, -1, 1, -1, 1]) ≈ rcor0
+@test crosscor([1:5;], [1:5;]) ≈ [-0.4, -0.4, -0.1, 0.4, 1.0, 0.4, -0.1, -0.4, -0.4]
+
+c11 = crosscor(x1, x1)
+c12 = crosscor(x1, x2)
+c21 = crosscor(x2, x1)
+c22 = crosscor(x2, x2)
+@test crosscor(realx1, realx2) ≈ c12
+
+@test crosscor(x,  x1) ≈ [c11 c21]
+@test crosscor(realx, realx1) ≈ [c11 c21]
+@test crosscor(x1, x)  ≈ [c11 c12]
+@test crosscor(realx1, realx)  ≈ [c11 c12]
+@test crosscor(x,  x)  ≈ cat([c11 c21], [c12 c22], dims=3)
+@test crosscor(realx, realx)  ≈ cat([c11 c21], [c12 c22], dims=3)
+
+
+## pacf
+
+rpacfr =  [-0.218158122381419,
+            0.195015316828711,
+            0.144315804606139,
+           -0.199791229449779]
+
+@test pacf(x[:,1], 1:4) ≈ rpacfr
+
+rpacfy = [-0.221173011668873,
+           0.189683314308021,
+           0.111857020733719,
+          -0.175020669835420]
+
+@test pacf(x[:,1], 1:4, method=:yulewalker) ≈ rpacfy