Formuła ACF i PACF

18

Chcę utworzyć kod do rysowania ACF i PACF na podstawie danych szeregów czasowych. Podobnie jak ten wygenerowany wykres z minitabu (poniżej).

Rysowanie ACF

Rysowanie PACF

Próbowałem przeszukać formułę, ale nadal nie rozumiem jej dobrze. Czy mógłbyś mi powiedzieć formułę i jak jej użyć, proszę? Jaka jest pozioma czerwona linia na wykresie ACF i PACF powyżej? Jaka jest formuła?

Dziękuję Ci,

Surya Dewangga
źródło
1
@javlacalle Czy podana formuła jest poprawna? Nie działałoby, jeśli n t = 1 ( y t - ˉ y ) < 0
ρ(k)=1nkt=k+1n(yty¯)(ytky¯)1nt=1n(yty¯)1nkt=k+1n(ytky¯),
prawda? Czy powinno to wyglądać następująco? $$ \ 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}} \ ,,
t=1n(yt-y¯)<0i / lubt=k+1n(yt-k-y¯)<0
zjazd
@conighion Masz rację, dzięki. Nie widziałem tego wcześniej. Naprawiłem to.
javlacalle

Odpowiedzi:

33

Autokorelacje

Korelacja między dwiema zmiennymi y1,y2) jest zdefiniowana jako:

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

gdzie E jest operatorem oczekiwanym, μ1 i μ2) są średnimi odpowiednio dla y1 i y2) oraz σ1,σ2) są ich odchyleniami standardowymi.

W kontekście pojedynczej zmiennej, czyli auto -correlation, y1 jest oryginalna seria i y2) jest opóźniony wersja niego. Zgodnie z powyższą definicją przykładowe autokorelacje rzędu k=0,1,2),...można otrzymać przez obliczanie następującego wyrażenia z obserwowanej seriiyt ,t=1,2),...,n :

ρ(k)=1n-kt=k+1n(yt-y¯)(yt-k-y¯)1nt=1n(yt-y¯)2)1n-kt=k+1n(yt-k-y¯)2),

gdzie y¯ jest średnią próbkową danych.

Częściowe autokorelacje

Częściowe autokorelacje mierzą zależność liniową jednej zmiennej po usunięciu efektu innych zmiennych wpływających na obie zmienne. Na przykład częściowa autokorelacja rzędu mierzy wpływ (zależność liniowa) yt-2) na yt po usunięciu efektu yt-1 zarówno na yt i yt-2) .

Każdą częściową autokorelację można uzyskać jako serię regresji postaci:

y~t=ϕ21y~t-1+ϕ22y~t-2)+mit,

gdzie y~t jest serią oryginalną minus średnia próbki, yt-y¯ . Oszacowanie ϕ22 da wartość częściowej autokorelacji rzędu 2. Po rozszerzeniu regresji o k dodatkowych opóźnień, oszacowanie ostatniego terminu da częściową autokorelację rzędu k .

Alternatywnym sposobem obliczenia przykładowej częściowej autokorelacji jest rozwiązanie następującego układu dla każdego rzędu k :

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

ρ()

# 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

Przedziały ufności

±z1-α/2)nz1-α/2)1-α/2)

±z1-α/2)1n(1+2)ja=1kρ(ja)2))

javlacalle
źródło
1
(+1) Dlaczego dwa różne przedziały ufności?
Scortchi - Przywróć Monikę
2
@Scortchi Stałe pasma są używane podczas testowania niezależności, podczas gdy rosnące pasma są czasami używane podczas identyfikacji modelu ARIMA.
javlacalle
1
Obie metody obliczania pasma ufności są wyjaśnione w nieco bardziej szczegółowo tutaj .
Scortchi - Przywróć Monikę
Idealne wyjaśnienie!
Jan Rothkegel
1
@javlacalle, robi wyrażenie dla ρ(k)przegapić kwadraty w mianowniku?
Christoph Hanck
10

„Chcę utworzyć kod do rysowania ACF i PACF na podstawie danych szeregów czasowych”.

Chociaż OP jest nieco niejasny, może być bardziej ukierunkowany na formułę kodującą w stylu „receptury” niż na formułę modelu algebry liniowej.


ACF jest dość prosta: mamy szereg czasowy, a w zasadzie zrobić wiele „kopie” (jak w „kopiuj i wklej”) o tym, rozumiejąc, że każda kopia ma być przesunięty o jeden wpis z wcześniejszej kopii, ponieważ początkowe dane zawierająt punkty danych, podczas gdy poprzednia długość szeregu czasowego (która wyklucza ostatni punkt danych) jest tylko t-1. Możemy wykonać praktycznie tyle kopii, ile jest rzędów. Każda kopia jest skorelowana z oryginałem, pamiętając, że potrzebujemy identycznych długości, i w tym celu będziemy musieli przycinać tylny koniec początkowej serii danych, aby były porównywalne. Na przykład, aby skorelować dane początkowetst-3) musimy pozbyć się ostatniego 3) punkty danych pierwotnego szeregu czasowego (pierwszy 3) chronologicznie).

Przykład:

Opracujemy szereg czasowy z cyklicznym wzorem sinusoidalnym nałożonym na linię trendu i hałas oraz wykreślymy ACF wygenerowany przez R. Ten przykład wziąłem z postu online autorstwa Christopha Scherbera i właśnie dodałem do niego szum:

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

wprowadź opis zdjęcia tutaj

Zwykle musielibyśmy przetestować dane pod kątem stacjonarności (lub po prostu spojrzeć na wykres powyżej), ale wiemy, że jest w tym trend, więc pomińmy tę część i przejdźmy bezpośrednio do etapu usuwania trendów:

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

wprowadź opis zdjęcia tutaj

Teraz jesteśmy gotowi przejąć te szeregi czasowe, najpierw generując ACF z acf()funkcją w R, a następnie porównując wyniki z prowizoryczną pętlą, którą zestawiłem:

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).

wprowadź opis zdjęcia tutaj


DOBRZE. To się udało. Do PACF . O wiele trudniejsze do zhakowania ... Chodzi tutaj o to, aby ponownie sklonować początkowe ts kilka razy, a następnie wybrać wiele punktów czasowych. Jednak zamiast po prostu korelować z początkowymi szeregami czasowymi, zebraliśmy wszystkie opóźnienia pomiędzy nimi i przeprowadziliśmy analizę regresji, aby wariancję wyjaśnioną przez poprzednie punkty czasowe można było wykluczyć (kontrolować). Na przykład, jeśli skupiamy się na PACF kończącym się na czastst-4, trzymamy tst, tst-1, tst-2) i tst-3), jak również tst-4i cofamy się tsttst-1+tst-2)+tst-3)+tst-4 poprzez pochodzenie i zachowując jedynie współczynnik dlatst-4:

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.
}

I na koniec rysowanie ponownie obok siebie, generowane R i obliczenia ręczne:

wprowadź opis zdjęcia tutaj

To, że pomysł jest poprawny, oprócz prawdopodobnych problemów obliczeniowych, można zauważyć w porównaniu PACFdo pacf(st.y, plot = F).


kod tutaj .

Antoni Parellada
źródło
1

Cóż, w praktyce znaleźliśmy błąd (szum), który jest reprezentowany przez mit przedziały ufności pomagają ustalić, czy poziom można uznać za jedynie hałas (ponieważ około 95% razy przypada na przedziały).

użytkownik120580
źródło
Witamy w CV, możesz rozważyć dodanie bardziej szczegółowych informacji na temat tego, w jaki sposób OP może to zrobić. Może także dodać informacje o tym, co reprezentuje każda linia?
Repmat,
1

Oto kod python do obliczenia 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
źródło
Hmmm Formatowanie kodu było złe:
Sada,