Testando um efeito causal (com 2 séries temporais)

Alguns dias atrás, voltei com uma frase que encontrei (em um jornal francês), em que alguém estava afirmando que

“… Uma variável antiga explica 85% da mudança em uma nova variável. Então podemos falar sobre causalidade ”

e tentei explicar que era estúpido: se considerarmos a regressão da temperatura no dia t + 1 contra o número de ciclistas no dia t, a R ^ 2 excede 80%… mas é difícil afirmar que o número de ciclistas em um dia específico realmente causa a temperatura no dia seguinte …

No entanto, isso foi frustrante, e eu queria saber se havia uma maneira inteligente de testar a causalidade nesse caso. Uma popular é a causalidade de Granger (posso mencionar um artigo que publicamos há alguns anos atrás, onde usamos esse teste, tendas, tweets e eventos: a interação entre protestos em andamento e mídias sociais). Para explicar esse teste, considere uma série temporal bivariada (como a que temos aqui), boldsymbol {z} _t = (x_t, y_t)e considere algum modelo autoregressivo bivariado
{ displaystyle { begin {bmatrix} x_ {t} \ y_ {t} end {bmatrix}} = { begin {bmatrix} c_ {1} \ c_ {2} end {bmatrix}} + { begin {bmatrix} a_ {1,1} & textcolor {red} {a_ {1,2}} \ textcolor {blue} {a_ {2,1}} e a_ {2,2} end {bmatrix }} { begin {bmatrix} x_ {t-1} \ y_ {t-1} end {bmatrix}} + { begin {bmatrix} u_ {t} \ v_ {t} end {bmatrix} }}Onde boldsymbol { varepsilon} _t = (u_t, v_t) é algum ruído branco bivariado, no sentido de que (i) { displaystyle mathbb {E} ( boldsymbol { varepsilon} _ {t}) = boldsymbol {0}} (o ruído é centralizado) (ii) { displaystyle mathbb {E} ( boldsymbol { varepsilon} _ {t} boldsymbol { varepsilon} _ {t} ^ top) = Omega} , portanto, a matriz de variância é constante, mas possivelmente não diagonal (iii) { displaystyle mathbb {E} ( boldsymbol { varepsilon} _ {t} boldsymbol { varepsilon} _ {t-h} ^ top) = boldsymbol {0}} para todos h neq 0. Observe que podemos usar a expressão simplificada{ displaystyle { boldsymbol {z} _t = boldsymbol {c} + boldsymbol {A} boldsymbol {z} _ {t-1} + boldsymbol { varepsilon} _t}}Agora, o teste Granger é baseado em várias quantidades. Com termos fora da diagonal da matriz Ómega, nós temos o chamado instantâneo causalidade e desde Ómega é simetria, vamos escrever x esquerda. Com termos fora da diagonal da matriz boldsymbol {A}, nós temos o chamado atrasado causalidade, com textcolor {blue} {x rightarrow y} ou textcolor {red} {x leftarrow y} (e possivelmente ambos, se os dois termos forem significativos).

Então, eu queria experimentar o meu problema de duas variáveis.

1
2
3
4

df = read.csv("cyclistsTempHKI.csv")
dfts = cbind(C=ts(df$cyclists,start = c(2014, 1,2), frequency = 365),
             T=ts(df$meanTemp,start = c(2014, 1,2), frequency = 365))
library(vars)

Agora tenho objetos de “série temporal” e podemos ajustar um modelo VAR,

1
2
3
4
5
6
7
8
9
10
11
12
13

var2 = VAR(dfts, p = 1, type = "const")
coefficients(var2)
$C
         Estimate   Std. Error   t value      Pr(>|t|)
C.l1    0.8684009   0.02889424 30.054460 8.080226e-107
T.l1   70.3042012  20.07247411  3.502518  5.102094e-04
const 807.6394001 187.75472482  4.301566  2.110412e-05
 
$T
           Estimate   Std. Error   t value     Pr(>|t|)
C.l1   0.0003865391 6.257596e-05  6.177118 1.540467e-09
T.l1   0.6611135594 4.347074e-02 15.208241 6.086394e-42
const -1.6413074565 4.066184e-01 -4.036481 6.446018e-05

Por um instante, podemos executar uma causalidade, para testar se o número de ciclistas pode causar a temperatura (no dia seguinte)

1
2
3
4
5
6
7

causality(var2, cause = "C")
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 38.157, df1 = 1, df2 = 842, p-value = 1.015e-09

Aqui, devemos rejeitar claramente H_0, que é que não há efeito causal. É assim que o estatístico diz que deve haver algum efeito causal entre o número de ciclistas e a temperatura …

Tão claramente, algo está errado aqui. Ou é algum tipo de superpotência que os ciclistas não conhecem. Ou esse teste usado por quarenta anos (Clive Granger chegou a receber um preço Nobel) não está funcionando. Ou perdemos alguma coisa. Na verdade … acho que perdemos alguma coisa aqui. As séries não são estacionárias. Quase podemos vê-lo com

1
2
3
4
5

Phi = matrix(c(coefficients(var2)$C[1:2,1],coefficients(var2)$T[1:2,1]),2,2)
eigen(Phi)
eigen() decomposition
$values
[1] 0.9594810 0.5700335

onde o valor próprio mais alto está muito próximo de um. Mas, na verdade, olhamos aqui para a temperatura …

ou seja, pelo menos, devemos esperar alguma raiz sazonal da unidade aqui. Então, vamos usar duas técnicas. O primeiro é uma diferença clássica de um ano, Delta_ {365} boldsymbol {z} _t = boldsymbol {z} _t- boldsymbol {z} _ {t-365}

1
2
3
4
5
6
7
8
9
10
11
12
13

var2 = VAR(diff(dfts,365), p = 1, type = "const")
coefficients(var2)
$C
          Estimate   Std. Error   t value     Pr(>|t|)
C.l1     0.8376424   0.07259969 11.537823 1.993355e-16
T.l1    42.2638410  28.58783276  1.478386 1.449076e-01
const -507.5514795 219.40240747 -2.313336 2.440042e-02
 
$T
         Estimate   Std. Error   t value     Pr(>|t|)
C.l1  0.000518209 0.0003277295 1.5812096 1.194623e-01
T.l1  0.598425288 0.1290511945 4.6371154 2.162476e-05
const 0.547828079 0.9904263469 0.5531235 5.823804e-01

O teste no modelo VAR adequado produz

1
2
3
4
5
6
7

causality(var2, cause = "C") 
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 2.5002, df1 = 1, df2 = 112, p-value = 0.1167

ou seja, com 11% p, devemos rejeitar a suposição de que o número de ciclistas causa a temperatura (no dia seguinte) e, na verdade, também devemos rejeitar o contrário

1
2
3
4
5
6
7

causality(var2, cause = "T") 
$Granger
 
	Granger causality H0: T do not Granger-cause C
 
data:  VAR object var2
F-Test = 2.1856, df1 = 1, df2 = 112, p-value = 0.1421

No entanto, se olharmos para a causalidade instantânea, essa faz mais sentido

1
2
3
4
5
6

$Instant
 
	H0: No instantaneous causality between: T and C
 
data:  VAR object var2
Chi-squared = 13.081, df = 1, p-value = 0.0002982

A segunda idéia seria usar uma diferença de um dia, Delta_ {1} boldsymbol {z} _t = boldsymbol {z} _t- boldsymbol {z} _ {t-1} e para ajustar um modelo VAR nesse

1
2
3
4

VARselect(diff(dfts,1), lag.max = 4, type="const")
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      3      2      3

mas nesse caso, um modelo VAR (1) – com apenas um atraso – pode não ser suficiente. Talvez seja melhor considerar um VAR (3)

1

var2 = VAR(diff(dfts,1), p = 3, type = "const")

e nessa, mais uma vez, devemos rejeitar o efeito causal do número de ciclistas na temperatura (no dia seguinte)

1
2
3
4
5
6
7

causality(var2, cause = "C")  
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 0.67644, df1 = 3, df2 = 828, p-value = 0.5666

e desta vez, poderia haver um efeito causal (atrasado) da temperatura no número de ciclistas

1
2
3
4
5
6
7
8
9
10
11
12
13
14

causality(var2, cause = "T")  
$Granger
 
	Granger causality H0: T do not Granger-cause C
 
data:  VAR object var2
F-Test = 7.7981, df1 = 3, df2 = 828, p-value = 3.879e-05
 
$Instant
 
	H0: No instantaneous causality between: T and C
 
data:  VAR object var2
Chi-squared = 55.83, df = 1, p-value = 7.905e-14

mas nada instantaneamente … Parece que a causalidade de Granger funciona bem nessa!

Relacionado