Inconsistency with explicit ARMA simulation

[Datseris]

I'm not sure what I am doing wrong in the source code for ARFIMA.jl, but take the following explicit simulation of an ARMA process and compare it with its ""equivalent" ARFIMA.jl syntax:

using Random
Random.seed!(77163)
η = randn(5000)
s = ones(5000)
for n in 4:5000
    s[n] = 1.625s[n-1] - 0.284s[n-2] - 0.355s[n-3] + η[n] - 0.96η[n-1]
end
s ./= std(s)

using ARFIMA
φ = SVector(1.625, -0.284, -0.355)
θ = SVector(0.96)
Random.seed!(77163)
s2 = arma(5000, 1.0, φ, θ)
s2 ./= std(s2)

using PyPlot
plot(s)
plot(s2)
xlim(0, 1000)

Both simulations have the characteristic oscillatory behavior of the specified ARMA process, but the one coming from ARFIMA.jl has much, much less noise.

Why...?

[ndgnuh]

Hi. I've tried to rewrite the noise and the auto regressive function in ARFIMA, but the result is the same (i.e. less noise in arma). So I tried to generate an AR and an ARMA series to compare.

let N = 1000,
    σ = 1f0,
    φ = SVector(1.625, -0.284, -0.355),
    θ = SVector(0.96)

    Random.seed!(1)
    s1 = arma(N, σ, φ)
    s1 = s1 ./ std(s1)
    Random.seed!(1)
    s2 = arma(N, σ, φ, θ)
    s2 = s2 ./ std(s2)

    p = plot()
    plot!(p, s1; label = "AR(3)")
    plot!(p, s2; label = "ARMA(3, 1)")
end

I suspect something is wrong with MA(q) sampling. Still, I'm no time series expert, just suggesting.

[Datseris]

Thanks a lot for having a look, I just went through the code that generates the moving average noise, but I Don't find the mistake... :(

function generate_noise(rng, N, σ, θ::Nothing)
    noise = randn(rng, N)
    if σ == 1.0
        return noise
    else
        return σ .* noise
    end
end

function generate_noise(rng, N, σ, θ::SVector{Q}) where {Q} # MA
    ε = generate_noise(rng, N+Q, σ, nothing)
    noise = zeros(N)
    θ = -θ # this change is necessary due to the defining equation
    # simply now do the average process
    @inbounds for i in (Q+1):N
        noise[i-Q] = bdp(θ, ε, i) + ε[i]
    end
    return noise
end

[mattiasvillani]

Yes, something is off with MA part. Can you please put a warning on the front page? Lost quite a bit of time here trying to debug my code that used used your package to simulate data. But thanks for doing this, hope that you can get it to work properly soon. Here is my check that AR(1) works, but MA(1) does not give the right autocorrelations:

`using ARFIMA, StatsBase

AR(1)

N = 10000 σ = 1 ϕ = 0.8 x = arfima(N, σ, nothing, SVector(ϕ)) estACF = autocor(x, 1:2) # First and second autocorrelation theoreticalACF = [ϕ,ϕ^2] estACF-theoreticalACF # Estimated Autocorrelation and theoretical are close. OK!

MA(1)

N = 10000 σ = 1 θ = -0.8 x = arfima(N, σ, nothing, nothing, SVector(θ)) estACF = autocor(x, 1:2) # First and second autocorrelation theoreticalACF = [-θ/(1+θ^2),0] estACF-theoreticalACF # Estimated Autocorrelation and theoretical are not close. Not OK!`

[Datseris]

The source code is really small, can you have a look and perhaps you find the mistake? I won't have time to fix this myself at the moment.

The moving average noise is generated here: https://github.com/JuliaDynamics/ARFIMA.jl/blob/master/src/ARFIMA.jl#L120-L129

but sure i'll put a warning

[Datseris]

version 0.3.1 will throw a warning anytime a non-nothing θ is given. Would be nice however if someone does a PR here (I'll review, but I won't fix this myself, as I don't know the problem).

[mattiasvillani]

I have very limited time now, but I will try to find time to look into this. Great with a warning.

[mattiasvillani]

If someone else has a go at this, here is a simulator for ARTFIMA (ARFIMA is the special case with \lambda = 0) using RCall that can be useful for testing:

using RCall R"library(artfima)"

function artsim(n, d, λ, ϕ, θ, μ, σ²) R""" x = artsim(n = $n, d = $d, lambda = $λ, phi = $ϕ, theta = $θ, mean = $μ, sigma2 = $σ²) """ @rget x return x end