【前編】初心者のための時系列間の因果関係検出

function calDist(x)
    exp(-0.5*x^2)
end

function pred(X,Y,σ;cross_validate=true)
    T = length(Y)
    upper = zeros(T)
    lower = zeros(T)
    for t=1:T, i=1:T
        if cross_validate && i == t continue end
        prodw = prod(calDist.((X[i,:] - X[t,:]) ./ σ))
        lower[t] += prodw
        upper[t] += prodw*Y[i]
    end
    upper ./ lower
end

function Cnpmr(Xy,Xxy,Yy;σy=NaN,σxy=NaN)
    Ŷy = pred(Xy,Yy,isnan(σy) ? ones(size(Xy)[2]) : σy)
    Ŷxy = pred(Xxy,Yy,isnan(σxy) ? ones(size(Xxy)[2]) : σxy)
    log(var(Ŷy-Yy)/var(Ŷxy-Yy))
end

function makeDataset(fromX::Array{Float64,2},
                toY::Array{Float64,2},d;Z=Array{Float64}(0,2))
    T = length(toY)-d
    Xy = vcat([toY[i:-1:(i-d+1)]' for i=d:(d+T-1)]...)
    if size(Z)[1] > 0
        Xz = vcat([collect((Z[i:-1:(i-d+1),:]'...))' for i=d:(d+T-1)]...)
        Xy = [Xy Xz]
    end
    Xx = vcat([collect((fromX[i:-1:(i-d+1),:]'...))' for i=d:(d+T-1)]...)
    Xxy = [Xx Xy]
    Yy = toY[(d+1):(d+T)]
    (Xy,Xx,Xxy,Yy[:,:])
end

function nudge(X::Array{Float64,2},col,∇)
    dk = ∇*abs(maximum(X[:,col]) - minimum(X[:,col]))
    nX = repeat(X; inner=[2,1])
    nX[:,col] += repeat([dk; -dk], outer=[size(X)[1],1])
    nX
end

function Q(X,Y,col;σ=NaN,∇=0.05)
    Ŷ = pred(X,Y,isnan(σ) ? ones(size(X)[2]) : σ,cross_validate=false)
    nX = nudge(X,col,∇)
    nY = repeat(Y,inner=[2,1])
    Ŷs = pred(nX,nY,isnan(σ) ? ones(size(X)[2]) : σ,cross_validate=false)
    Qcol = sum(abs(Ŷs - repeat(Ŷ,inner=[2,1]))) /
            (2*length(Y)*∇*abs(maximum(Y)-minimum(Y)))
end

function sensitivity(toY::Array{Float64,2},fromX::Array{Float64,2};d=2)
    (Xy,Xx,Xxy,Yy) = makeDataset(fromX,toY,d)
    xw = size(fromX)[2]
    for col in 1:(d*xw)
        Qcol = Q(Xxy,Yy,col)
        println("Q_x$(col%xw==0?xw:col%xw)(t-$(Int(ceil(col/xw)))): $Qcol")
    end
end

function surrogate(X;lower=33,upper=90)
    circshift(X,Int(floor(size(X)[1]*rand(lower:upper)/100)))
end

function isCausal(
            fromX::Array{Float64,2},
            toY::Array{Float64,2};d=2,α=0.05,Z=Array{Float64}(0,2))
    (Xy,Xx,Xxy,Yy) = makeDataset(fromX,toY,d,Z=Z)
    Cmain = Cnpmr(Xy,Xxy,Yy)
    println("Cnpmr(X->Y): $Cmain")
    if Cmain < 0 return false end
    surnum = Int(floor(1/α - 1))
    for cnt in 1:surnum
        sfromX = surrogate(fromX)
        (sXy,sXx,sXxy,Yy) = makeDataset(sfromX,toY,d)
        Csurr = Cnpmr(sXy,sXxy,Yy)
        println("$cnt: $Csurr")
        if Csurr > Cmain return false end
    end
    true
end

function batchTest(datfn;d=3,trials = 10)
    x1_x2sigcnt = 0
    x2_x1sigcnt = 0
    (X1,X2,Z) = datfn()
    for i = 1:trials
        if isCausal(X1,X2,d=d) x1_x2sigcnt += 1 end
        if isCausal(X2,X1,d=d) x2_x1sigcnt += 1 end
    end
    println("X1 -> X2: $x1_x2sigcnt / $trials significant count(s)")
    println("X2 -> X1: $x2_x1sigcnt / $trials significant count(s)")
end