Current section
Files
Jump to
Current section
Files
lib/ljung_box.ex
defmodule LjungBox do
@moduledoc """
The Ljung-Box test for autocorrelation in time series residuals.
The Ljung-Box Q-statistic tests the null hypothesis that a series of
residuals are independently distributed (i.e., exhibit no autocorrelation).
It is commonly used to validate ARIMA model residuals.
## Formula
The Q-statistic is defined as:
Q = n(n+2) * Σ_{k=1}^{h} ρ̂_k² / (n - k)
where `n` is the number of observations, `h` is the number of lags tested,
and `ρ̂_k` is the sample autocorrelation at lag `k`.
Under the null hypothesis of no autocorrelation, `Q` follows a chi-squared
distribution with `h` degrees of freedom (or `h - p - q` when applied to
ARMA(p, q) residuals).
## Interpretation
* **p-value < 0.05**: Reject the null hypothesis; there is significant
autocorrelation at the tested lags.
* **p-value ≥ 0.05**: Fail to reject the null hypothesis; the series is
consistent with white noise at the tested lags.
See `test/2` for usage examples.
"""
@doc """
Runs the Ljung-Box test on `series` up to `lags` lags.
Returns a map with:
* `:statistic` – the Q value
* `:p_value` – probability of observing a Q this large under H₀
* `:lags` – number of lags tested
* `:n` – number of observations
A **low p-value** (typically < 0.05) indicates significant autocorrelation.
## Examples
iex> result = LjungBox.test([1, 2, 1, 2, 1, 2, 1, 2, 1, 2], 2)
iex> result.lags
2
iex> result.n
10
"""
@spec test(list(number()), pos_integer()) :: %{
statistic: float(),
p_value: float(),
lags: pos_integer(),
n: pos_integer()
}
def test(series, lags \\ 10) when is_list(series) and is_integer(lags) and lags > 0 do
n = length(series)
if lags >= n do
raise ArgumentError,
"lags (#{lags}) must be less than the length of the series (#{n})"
end
q = statistic(series, lags)
p = 1.0 - chi_squared_cdf(q, lags)
%{statistic: q, p_value: p, lags: lags, n: n}
end
@doc """
Computes the Ljung-Box Q-statistic for `series` at up to `lags` lags.
## Examples
iex> q = LjungBox.statistic([1, 2, 3, 4], 2)
iex> Float.round(q, 2)
1.58
"""
@spec statistic(list(number()), pos_integer()) :: float()
def statistic(series, lags) when is_list(series) and is_integer(lags) and lags > 0 do
n = length(series)
rhos = autocorrelations(series, lags)
sum =
rhos
|> Enum.with_index(1)
|> Enum.reduce(0.0, fn {rho, k}, acc -> acc + rho * rho / (n - k) end)
n * (n + 2) * sum
end
@doc """
Computes sample autocorrelations of `series` at lags 1 through `max_lag`.
Returns a list of `max_lag` floats.
## Examples
iex> length(LjungBox.autocorrelations([1, 2, 3, 4, 5], 3))
3
"""
@spec autocorrelations(list(number()), pos_integer()) :: list(float())
def autocorrelations(series, max_lag)
when is_list(series) and is_integer(max_lag) and max_lag > 0 do
for lag <- 1..max_lag, do: autocorrelation(series, lag)
end
@doc """
Computes the sample autocorrelation of `series` at the given `lag`.
Returns `0.0` when the series has zero variance (constant series).
## Examples
iex> LjungBox.autocorrelation([1, 2, 3, 4], 1)
0.25
iex> LjungBox.autocorrelation([1, 2, 3, 4], 0)
1.0
iex> LjungBox.autocorrelation([5, 5, 5, 5], 1)
0.0
"""
@spec autocorrelation(list(number()), non_neg_integer()) :: float()
def autocorrelation(series, lag)
when is_list(series) and is_integer(lag) and lag >= 0 do
n = length(series)
mean = Enum.sum(series) / n
variance =
Enum.reduce(series, 0.0, fn x, acc -> acc + (x - mean) * (x - mean) end) / n
if variance == 0.0 do
0.0
else
covariance =
series
|> Enum.zip(Enum.drop(series, lag))
|> Enum.reduce(0.0, fn {x, y}, acc -> acc + (x - mean) * (y - mean) end)
|> Kernel./(n)
covariance / variance
end
end
# ---------------------------------------------------------------------------
# Private helpers: chi-squared CDF via regularized incomplete gamma function
# ---------------------------------------------------------------------------
# P(X ≤ x | df = k) = P(k/2, x/2) where P is the regularized lower
# incomplete gamma function.
defp chi_squared_cdf(x, _k) when x <= 0, do: 0.0
defp chi_squared_cdf(x, k) do
regularized_gamma(k / 2.0, x / 2.0)
end
# Regularized lower incomplete gamma P(a, x).
# Uses a series expansion when x < a + 1, and a continued fraction otherwise.
defp regularized_gamma(_a, x) when x <= 0, do: 0.0
defp regularized_gamma(a, x) do
if x < a + 1.0 do
gamma_series(a, x)
else
1.0 - gamma_continued_fraction(a, x)
end
end
# Series expansion for P(a, x).
defp gamma_series(a, x) do
initial = 1.0 / a
{sum, _} =
Enum.reduce_while(1..300, {initial, initial}, fn n, {sum, term} ->
new_term = term * x / (a + n)
new_sum = sum + new_term
if abs(new_term) < abs(new_sum) * 1.0e-12 do
{:halt, {new_sum, new_term}}
else
{:cont, {new_sum, new_term}}
end
end)
sum * :math.exp(-x + a * :math.log(x) - log_gamma(a))
end
# Continued fraction expansion for Q(a, x) = 1 - P(a, x) via Lentz's method.
defp gamma_continued_fraction(a, x) do
fpmin = 1.0e-300
b0 = x + 1.0 - a
c0 = 1.0 / fpmin
d0 = if abs(b0) < fpmin, do: fpmin, else: 1.0 / b0
h0 = d0
{h, _, _} =
Enum.reduce_while(1..300, {h0, c0, d0}, fn i, {h, c, d} ->
an = -i * (i - a)
b = x + 1.0 - a + 2.0 * i
d = an * d + b
d = if abs(d) < fpmin, do: fpmin, else: d
c = b + an / c
c = if abs(c) < fpmin, do: fpmin, else: c
d = 1.0 / d
delta = d * c
new_h = h * delta
if abs(delta - 1.0) < 1.0e-12 do
{:halt, {new_h, c, d}}
else
{:cont, {new_h, c, d}}
end
end)
h * :math.exp(-x + a * :math.log(x) - log_gamma(a))
end
# Log-gamma function via the Lanczos approximation (g = 7, 9 terms).
# Coefficients from Numerical Recipes, 3rd ed.
@lanczos_g 7
@lanczos_c [
0.99999999999980993,
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7
]
defp log_gamma(z) when z < 0.5 do
:math.log(:math.pi() / :math.sin(:math.pi() * z)) - log_gamma(1.0 - z)
end
defp log_gamma(z) do
z = z - 1.0
x =
@lanczos_c
|> Enum.with_index()
|> Enum.reduce(0.0, fn
{c, 0}, _acc -> c
{c, i}, acc -> acc + c / (z + i)
end)
t = z + @lanczos_g + 0.5
0.5 * :math.log(2 * :math.pi()) + (z + 0.5) * :math.log(t) - t + :math.log(x)
end
end