Fórmula ACF e PACF

18

Quero criar um código para plotagem do ACF e PACF a partir de dados de séries temporais. Assim como esse gráfico gerado no minitab (abaixo).

Plotagem ACF

Plotagem PACF

Eu tentei pesquisar a fórmula, mas ainda não a entendo bem. Você se importaria de me dizer a fórmula e como usá-la, por favor? Qual é a linha vermelha horizontal no gráfico ACF e PACF acima? Qual é a fórmula?

Obrigado,

Surya Dewangga
fonte
11
@javlacalle A fórmula que você forneceu está correta? Não funcionaria se certo? Deve ser o seguinte? $$ \ rho (k) = \ frac {\ frac {1} {nk} \ sum_ {t = k + 1} ^ n (y_t - \ bar {y}) (y_ {tk} - \ bar {y} )} {\ sqrt {\ frac {1} {n} \ sum_ {t = 1} ^ n (y_t - \ bar {y}) ^ 2} \ sqrt {\ frac {1} {nk} \ sum_ {t = k + 1} ^ n (y_ {tk} - \ bar {y}) ^ 2}} \ ,,n t = 1 ( y t - ˉ y ) < 0
ρ(k)=1 1n-kt=k+1 1n(yt-y¯)(yt-k-y¯)1 1nt=1 1n(yt-y¯)1 1n-kt=k+1 1n(yt-k-y¯),
t=1n(yty¯)<0and/ort=k+1n(ytky¯)<0 0
conighion
@Conighion Você está certo, obrigado. Eu não vi isso antes. Eu consertei isso.
Javlacalle 13/10/19

Respostas:

33

Autocorrelações

A correlação entre duas variáveis y1,y2 é definida como:

ρ=E[(y1μ1)(y2μ2)]σ1σ2=Cov(y1,y2)σ1σ2,

onde E é o operador de expectativa, μ1 e μ2 são as médias respectivamente para y1 e e são seus desvios padrão.y2σ1,σ2

No contexto de uma única variável, ou seja, correlação automática , é a série original e é uma versão atrasada dela. Após a definição acima, autocorrelações amostra de ordem k = 0 , 1 , 2 , . . . pode ser obtido utilizando a seguinte expressão, com a série observada y t , t = 1 , 2 , . . . , n :y1y2k=0,1,2,...ytt=1,2,...,n

ρ(k)=1nkt=k+1n(yty¯)(ytky¯)1nt=1n(yty¯)21nkt=k+1n(ytky¯)2,

onde y¯ é a média da amostra dos dados.

Autocorrelações parciais

Autocorrelações parciais medem a dependência linear de uma variável após remover o efeito de outras variáveis ​​que afetam as duas variáveis. Por exemplo, a autocorrelação parcial de ordem mede o efeito (dependência linear) de yt2 em yt após remover o efeito de yt1 em yt e yt2 .

Cada autocorrelação parcial pode ser obtida como uma série de regressões da forma:

y~t=ϕ21y~t1+ϕ22y~t2+et,

onde y~t é a série original menos a média da amostra, yty¯ . A estimativa de ϕ22 fornecerá o valor da autocorrelação parcial da ordem 2. Estendendo a regressão com k atrasos adicionais, a estimativa do último termo fornecerá a autocorrelação parcial da ordem k .

Uma maneira alternativa de calcular as autocorrelações parciais da amostra é resolver o seguinte sistema para cada ordem k :

(ρ(0)ρ(1)ρ(k1)ρ(1)ρ(0)ρ(k2)ρ(k1)ρ(k2)ρ(0))(ϕk1ϕk2ϕkk)=(ρ(1)ρ(2)ρ(k)),

onde ρ() são as autocorrelações da amostra. Esse mapeamento entre as autocorrelações da amostra e as autocorrelações parciais é conhecido como recursão de Durbin-Levinson . Essa abordagem é relativamente fácil de implementar para ilustração. Por exemplo, no software R, podemos obter a autocorrelação parcial da ordem 5 da seguinte maneira:

# sample data
x <- diff(AirPassengers)
# autocorrelations
sacf <- acf(x, lag.max = 10, plot = FALSE)$acf[,,1]
# solve the system of equations
res1 <- solve(toeplitz(sacf[1:5]), sacf[2:6])
res1
# [1]  0.29992688 -0.18784728 -0.08468517 -0.22463189  0.01008379
# benchmark result
res2 <- pacf(x, lag.max = 5, plot = FALSE)$acf[,,1]
res2
# [1]  0.30285526 -0.21344644 -0.16044680 -0.22163003  0.01008379
all.equal(res1[5], res2[5])
# [1] TRUE

Faixas de confiança

As faixas de confiança podem ser calculadas como o valor das autocorrelações da amostra ±z1α/2n , ondez1α/2é o quantil1α/2na distribuição gaussiana, por exemplo, 1,96 para faixas de confiança de 95%.

Às vezes, são usadas faixas de confiança que aumentam à medida que a ordem aumenta. Nesse caso, as bandas podem ser definidas como ±z1α/21n(1+2i=1kρ(i)2).

javlacalle
fonte
11
(+1) Por que as duas faixas de confiança diferentes?
Scortchi - Restabelece Monica
2
@ Scortchi As bandas constantes são usadas no teste de independência, enquanto as bandas crescentes às vezes são usadas na identificação de um modelo ARIMA.
Javlacalle
11
Os dois métodos para calcular faixas de confiança são explicados um pouco mais detalhadamente aqui .
Scortchi - Restabelece Monica
Explicação perfeita!
Jan Rothkegel 13/11/19
11
ρ(k)
9

"Quero criar um código para plotagem do ACF e PACF a partir de dados de séries temporais".

Embora o OP seja um pouco vago, ele pode ser mais direcionado a uma formulação de codificação no estilo "receita" do que a uma formulação de modelo de álgebra linear.


tt1tst333

Exemplo:

Vamos inventar uma série temporal com um padrão senoidal cíclico sobreposto a uma linha de tendência e ruído, e traçar o ACF gerado por R. Eu peguei este exemplo em um post online de Christoph Scherber e apenas adicionei o ruído a ele:

x=seq(pi, 10 * pi, 0.1)
y = 0.1 * x + sin(x) + rnorm(x)
y = ts(y, start=1800)

enter image description here

Normalmente, teríamos que testar os dados quanto à estacionariedade (ou apenas olhar para o gráfico acima), mas sabemos que há uma tendência, então vamos pular esta parte e ir diretamente para a etapa de descendência:

model=lm(y ~ I(1801:2083))
st.y = y - predict(model)

enter image description here

Agora, estamos prontos para resolver essa série cronológica, primeiro gerando o ACF com a acf()função em R e comparando os resultados com o loop improvisado que montei:

ACF = 0                  # Starting an empty vector to capture the auto-correlations.
ACF[1] = cor(st.y, st.y) # The first entry in the ACF is the correlation with itself (1).
for(i in 1:30){          # Took 30 points to parallel the output of `acf()`
  lag = st.y[-c(1:i)]    # Introducing lags in the stationary ts.
  clipped.y = st.y[1:length(lag)]    # Compensating by reducing length of ts.
  ACF[i + 1] = cor(clipped.y, lag)   # Storing each correlation.
}
acf(st.y)                            # Plotting the built-in function (left)
plot(ACF, type="h", main="ACF Manual calculation"); abline(h = 0) # and my results (right).

enter image description here


tst4tsttst1tst2tst3, as well as tst4, and we regress tsttst1+tst2+tst3+tst4 through the origin and keeping only the coefficient for tst4:

PACF = 0          # Starting up an empty storage vector.
for(j in 2:25){   # Picked up 25 lag points to parallel R `pacf()` output.
  cols = j        
  rows = length(st.y) - j + 1 # To end up with equal length vectors we clip.

  lag = matrix(0, rows, j)    # The storage matrix for different groups of lagged vectors.

for(i in 1:cols){
  lag[ ,i] = st.y[i : (i + rows - 1)]  #Clipping progressively to get lagged ts's.
}
  lag = as.data.frame(lag)
  fit = lm(lag$V1 ~ . - 1, data = lag) # Running an OLS for every group.
  PACF[j] = coef(fit)[j - 1]           # Getting the slope for the last lagged ts.
}

And finally plotting again side-by-side, R-generated and manual calculations:

enter image description here

That the idea is correct, beside probable computational issues, can be seen comparing PACF to pacf(st.y, plot = F).


code here.

Antoni Parellada
fonte
1

Well, in the practise we found error (noise) which is represented by et the confidence bands help you to figure out if a level can be considerate as only noise (because about the 95% times will be into the bands).

user120580
fonte
Welcome to CV, you might want to consider adding some more detailed information on how OP would go about do this specifically. Maybe also add some information on what each line represents?
Repmat
1

Here is a python code to compute ACF:

def shift(x,b):
    if ( b <= 0 ):
        return x
    d = np.array(x);
    d1 = d
    d1[b:] = d[:-b]
    d1[0:b] = 0
    return d1

# One way of doing it using bare bones
# - you divide by first to normalize - because corr(x,x) = 1
x = np.arange(0,10)
xo = x - x.mean()

cors = [ np.correlate(xo,shift(xo,i))[0]  for i in range(len(x1)) ]
print (cors/cors[0] )

#-- Here is another way - you divide by first to normalize
cors = np.correlate(xo,xo,'full')[n-1:]
cors/cors[0]
Sada
fonte
Hmmm Code formatting was bad:
Sada