Simple autoregressive surrogates (#77)

* initial autoregresive surrogator

* add autoregressive codfe

* add docstring

* runnable code

* fix autoregressive length

* add AR to the docs

* add test for autoregressive

* add changelog / increment version

* fix forgotten < lol

Author: Datseris

CHANGELOG.md

@@ -1,4 +1,6 @@
 *Changelog is kept with respect to version 1.0. This software follows SymVer2.0*
 
+# 1.2
+New surrogate methods: `AutoRegressive`
 # 1.1
 * New surrogate methods: `CycleShuffle`, `ShuffleDimensions`, `CircShift`

Project.toml

@@ -2,7 +2,7 @@ name = "TimeseriesSurrogates"
 uuid = "c804724b-8c18-5caa-8579-6025a0767c70"
 authors = ["Kristian Agasøster Haaga <kahaaga@gmail.com>", "George Datseris"]
 repo = "https://github.com/kahaaga/TimeseriesSurrogates.jl.git"
-version = "1.1.1"
+version = "1.2.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

docs/src/index.md

@@ -32,6 +32,7 @@ TFTS
 AAFT
 TAAFT
 IAAFT
+AutoRegressive
 PseudoPeriodic
 WLS
 ShuffleDimensions

src/TimeseriesSurrogates.jl

@@ -29,6 +29,7 @@ include("methods/truncated_fourier.jl")
 include("methods/wavelet_based.jl")
 include("methods/pseudoperiodic.jl")
 include("methods/multidimensional.jl")
+include("methods/ar.jl")
 
 # Visualization routine for time series + surrogate + periodogram/acf/histogram
 using Requires

src/methods/ar.jl

@@ -0,0 +1,54 @@
+using DSP, Random
+export AutoRegressive
+
+"""
+    AutoRegressive(n, method = LPCLevinson())
+
+Autoregressive surrogates of order-`n`. The autoregressive coefficients `φ` are estimated
+using `DSP.lpc(x, n, method)`, and thus see the documentation of DSP.jl for possible
+`method`s.
+
+While these surrogates are obviously suited to test the null hypothesis whether the data
+are coming from a autoregressive process, the Fourier Transform-based surrogates are
+probably a better option. The current method is more like an explicit way to
+produce surrogates for the same hypothesis by fitting a model.
+It can be used as convient way to estimate
+autoregressive coefficients and automatically generate surrogates based on them.
+
+The coefficients φ of the autoregressive fit can be found by doing
+```julia
+sg = surrogenerator(x, AutoRegressive(n))
+φ = sg.init.φ
+```
+"""
+struct AutoRegressive{M} <: Surrogate
+    n::Int
+    m::M
+end
+AutoRegressive(n::Int) = AutoRegressive(n, LPCLevinson())
+
+function surrogenerator(x, method::AutoRegressive)
+    φ, e = lpc(x, method.n, method.m)
+    init = (d = Normal(0, std(x)), φ = φ)
+    return SurrogateGenerator(method, x, init)
+end
+
+function (sg::SurrogateGenerator{<:AutoRegressive})()
+    N, φ, d = length(sg.x), sg.init.φ, sg.init.d
+    return autoregressive(d, N, φ)
+end
+
+function autoregressive(d::Normal{T}, N, φ) where {T}
+    f = length(φ)
+    z = zeros(T, N+f)
+    @inbounds for i in 1:f; z[i] = rand(d); end
+    @inbounds for i in f+1:length(z)
+        s = zero(T)
+        for j in 1:f
+            s += φ[j]*z[i-j]
+        end
+        s += rand(d)
+        z[i] = s
+    end
+    return view(z, f+1:N+f)
+end

test/runtests.jl

@@ -48,6 +48,14 @@ end
     @test all([s[i] ∈ x for i = 1:N])
 end
 
+@testset "AutoRegressive" begin
+	y = TimeseriesSurrogates.AR1(2000, 0.1, 0.5)
+    sg = surrogenerator(y, AutoRegressive(1))
+	@test 0.4 ≤ abs(sg.init.φ[1]) ≤ 0.6
+	s = sg()
+    @test length(s) == length(y)
+end
+
 @testset "AAFT" begin
     aaft = AAFT()
     s = surrogate(x, aaft)