julia - 如何在 Julia 中获得滚动窗口回归

Question

假设我有股票的价格，并且我想在给定窗口大小的情况下以滚动方式找到回归线的斜率。我怎样才能在 Julia 中完成它？我希望它非常快，因此不想使用 for 循环。

Answer 1

一般来说，您不应该担心 for Julia 中的循环，因为它们没有 R 或 Python for 循环的开销。因此，您只需要担心渐进复杂性，而不用担心解释器开销引入的潜在大常数因子。

尽管如此，与使用朴素的 O(n²) 切片和回归方法相比，使用卷积可以更有效地(渐近地)完成此操作。 DSP.jl包提供卷积功能。以下是一个没有截距的示例(它计算滚动贝塔值)；通过修改公式应该可以支持拦截。

using DSP

# Create some example x (signal) and y (stock prices)
# such that strength of signal goes up over time
const x = randn(100)
const y = (1:100) .* x .+ 100 .* randn(100)

# Create the rolling window
const window = Window.rect(20)

# Compute linear least squares estimate (X^T X)^-1 X^T Y
const xᵗx = conv(x .* x, window)[length(window):end-length(window)+1]
const xᵗy = conv(x .* y, window)[length(window):end-length(window)+1]
const lls = xᵗy ./ xᵗx  # desired beta

# Check result against naïve for loop
const βref = [dot(x[i:i+19], y[i:i+19]) / dot(x[i:i+19], x[i:i+19]) for i = 1:81]
@assert isapprox(βref, lls)

编辑添加:支持拦截，即 X = [x 1] ，所以X^T X = [dot(x, x) sum(x); sum(x) w]哪里w是窗口大小，利用二维矩阵的逆公式可以得到(X^T X)^-1 = [w -sum(x); -sum(x) dot(x, x)]/(w * dot(x, x) - sum(x)^2) 。因此，[β, α] = [w dot(x, y) - sum(x) * sum(y), dot(x, x) * sum(y) - sum(x) * dot(x, y)] / (w * dot(x, x) - sum(x)^2) 。这可以翻译为以下卷积代码:

# Compute linear least squares estimate with intercept
const w = length(window)
const xᵗx = conv(x .* x, window)[w:end-w+1]
const xᵗy = conv(x .* y, window)[w:end-w+1]
const 𝟙ᵗx = conv(x, window)[w:end-w+1]
const 𝟙ᵗy = conv(y, window)[w:end-w+1]
const denom = w .* xᵗx - 𝟙ᵗx .^ 2
const α = (xᵗx .* 𝟙ᵗy .- 𝟙ᵗx .* xᵗy) ./ denom
const β = (w .* xᵗy .- 𝟙ᵗx .* 𝟙ᵗy) ./ denom

# Check vs. naive solution
const ref = vcat([([x[i:i+19] ones(20)] \ y[i:i+19])' for i = 1:81]...)
@assert isapprox([β α], ref)

请注意，对于具有不同窗口形状的加权最小二乘法，需要进行一些小的修改才能解开 length(window)和sum(window)它们在上面的代码中可以互换使用。