Implementation of PCMCI

[bgroenks96]

It would be nice to have a Julia-native implementation of the (relatively) new PCMCI [1] algorithm for causal inference in time series proposed by @jakobrunge and his team.

The algorithm essentially consists of two steps:

(1) PC1 condition selection to identify relevant adjacent links
(2) Momentary conditional independence (MCI) test (eq 6 in the paper) to identify time-shifted causal links.

I highly recommend reading the paper if you haven't already.

There is a reference implementation here, but bear in mind that the code is licensed under GPL 3.0 so it is incompatible with the MIT license for CausalInference.jl. Fortunately, CausalInference.jl already has PC and partial correlation conditional independence tests so really all that needs to be done is adding MCI, which should be doable from the definition given in the paper.

[1] Runge, J., Nowack, P., Kretschmer, M., Flaxman, S., & Sejdinovic, D. (2019). Detecting and quantifying causal associations in large nonlinear time series datasets. Science Advances, 5(11), eaau4996.

[schnirz]

Thanks for the suggestion, I will have a look at the article.

[bgroenks96]

I'm happy to help out with a PR as well.

[schnirz]

Yeah, sounds good. I’ll let you know once I’ve gone through the article and reference implementation.
It’s also worth pointing out that conditional mutual independence tests for conditional independence are also already implemented.

[jakobrunge]

Hi there. Good idea. But it won't be that easy since PCMCI is not just PC-algorithm plus some more tests. It is adapted to time series and you cannot just apply the original PC algorithm to time series. But the pseudo-code is all there in the paper:)

[mschauer]

Perhaps it helps that our implementations are indistinguishable from pseudo-code 😉

function skeleton(n::V, I, par...; kwargs...) where {V}
    g = complete_graph(n)
    S = Dict{edgetype(g),Vector{V}}()
    d = 0 # depth

    while true
        isdone = true
        for e0 in collect(edges(g)) # cannot remove edges while iterating
            for e in (e0, reverse(e0))
                nb = neighbors(g, src(e))
                if length(nb) > d  # i.e. |nb\{dst(e)}| >= d
                    r = searchsorted(nb, dst(e))
                    isempty(r) && continue # the edge does not exist anymore
                    isdone = false
                    for s in combinations_without(nb, d,  first(r)) # do not modify s!                nb0 = neighbors(g, src(e))
                        if I(src(e), dst(e), s, par...; kwargs...)
                            @debug "Removing edge $(e0) given $(s)"
                            rem_edge!(g, e0)
                            if !(e0 in keys(S))
                                S[e0] = s
                            end
                            break
                        end
                    end
                end
            end
        end
        d = d + 1
        if isdone
            return g, S
        end
    end
end