Addition of Yule-Walker Analysis Method

[josemanuel22]

This pull request proposes the addition of a method to perform Yule-Walker analysis in Julia, analogous to the yule_walker method available in Python's statsmodels.regression.linear_model. The aim is to extend Julia's statistical analysis capabilities, specifically for autoregressive processes, by implementing a method that estimates the parameters of an AR model using the Yule-Walker equations.

https://stackoverflow.com/questions/77107822/how-to-do-a-yule-walker-analysis-in-julia/77117406#77117406

Details:

1. The new method will calculate the autoregressive coefficients and the variance of the error term.
2. The implementation will focus on efficiency and accuracy, leveraging Julia's strengths in handling mathematical operations.
3. Unit tests will be added to ensure the method's reliability and to maintain the integrity of the package.
   I look forward to feedback and suggestions for improvements.

[andreasnoack]

Hello @josemanuel22. Thanks for the PR. However, please see

StatsBase.jl/src/signalcorr.jl

Lines 570 to 608 in 77e63df

	function pacf_yulewalker!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}, mk::Integer) where T<:Union{Float32, Float64}
	tmp = Vector{T}(undef, mk)
	for j = 1 : size(X,2)
	acfs = autocor(X[:,j], 1:mk)
	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]
	end
	end
	end
	"""
	pacf!(r, X, lags; method=:regression)
	Compute the partial autocorrelation function (PACF) of a matrix `X` at `lags` and
	store the result in `r`. `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
	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}
	lx = size(X, 1)
	m = length(lags)
	minlag, maxlag = extrema(lags)
	(0 <= minlag && 2maxlag < lx) || error("Invalid lag value.")
	size(r) == (m, size(X,2)) || throw(DimensionMismatch())
	if method == :regression
	pacf_regress!(r, X, lags, maxlag)
	elseif method == :yulewalker
	pacf_yulewalker!(r, X, lags, maxlag)
	else
	error("Invalid method: $method")
	end
	return r
	end

[josemanuel22]

@andreasnoack Thank you very much for your comment. I had not seen it, but it seems that this function is not documented, nor is it exported. It is for internal use, maybe it could be documented and exported. What do you think?

[andreasnoack]

Sorry for being a bit brief here. I made the comment a bit in a hurry. The point was mostly that you can compute the Yule-Walker estimates efficiently by reusing existing functions, i.e.

Yule-Walker

julia> -StatsBase.durbin(autocor(x, 1:4))
4-element Vector{Float64}:
  1.2600631466569998
 -0.4558455606978694
 -0.17590708455408563
  0.011391047373061776

Least-squares

julia> [ones(length(x) - k) [x[k - j + i] for i in 1:length(x)-k, j in 1:k]]\x[k+1:end]
5-element Vector{Float64}:
 28.267780791960924
  1.2864783021078277
 -0.48686426630545226
 -0.16910686502522418
  0.017892893292120397

Thinking a bit more about it, I think a function for estimating coefficients might also be better located in a dedicated time series package. The acf and pacf functions are a bit more like summary statistics but if they were added today, I guess they might also go in to a dedicated time series package instead of this one.

[josemanuel22]

Perfect, thank you for the clarification. I am not aware of any time series package in Julia that has this functionality, as it seems they resort to StatsBase.jl. I'm not sure if it would make sense to open a PR there. With this, I believe we can close this PR.