Formuła ACF i PACF

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

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,

Autokorelacje

Korelacja między dwiema zmiennymi $y_1, y_2$ jest zdefiniowana jako:

$$\rho = \frac{\hbox{E}\left[(y_1-\mu_1)(y_2-\mu_2)\right]}{\sigma_1 \sigma_2} = \frac{\hbox{Cov}(y_1, y_2)}{\sigma_1 \sigma_2} \,,$$

gdzie E jest operatorem oczekiwanym, $\mu_1$ i $\mu_2$ są średnimi odpowiednio dla $y_1$ i $y_2$ oraz $\sigma_1, \sigma_2$ są ich odchyleniami standardowymi.

W kontekście pojedynczej zmiennej, czyli auto -correlation, $y_1$ jest oryginalna seria i $y_2$ 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 serii$y_t$ ,$t=1,2,...,n$ :

$$\rho(k) = \frac{\frac{1}{n-k}\sum_{t=k+1}^n (y_t - \bar{y})(y_{t-k} - \bar{y})}{ \sqrt{\frac{1}{n}\sum_{t=1}^n (y_t - \bar{y})^2}\sqrt{\frac{1}{n-k}\sum_{t=k+1}^n (y_{t-k} - \bar{y})^2}} \,,$$

gdzie $\bar{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) $y_{t-2}$ na $y_t$ po usunięciu efektu $y_{t-1}$ zarówno na $y_t$ i $y_{t-2}$ .

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

$$\tilde{y}_t = \phi_{21} \tilde{y}_{t-1} + \phi_{22} \tilde{y}_{t-2} + e_t \,,$$

gdzie $\tilde{y}_t$ jest serią oryginalną minus średnia próbki, $y_t - \bar{y}$ . Oszacowanie $\phi_{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$ :

$$\begin{eqnarray} \left(\begin{array}{cccc} \rho(0) & \rho(1) & \cdots & \rho(k-1) \\ \rho(1) & \rho(0) & \cdots & \rho(k-2) \\ \vdots & \vdots & \vdots & \vdots \\ \rho(k-1) & \rho(k-2) & \cdots & \rho(0) \\ \end{array}\right) \left(\begin{array}{c} \phi_{k1} \\ \phi_{k2} \\ \vdots \\ \phi_{kk} \\ \end{array}\right) = \left(\begin{array}{c} \rho(1) \\ \rho(2) \\ \vdots \\ \rho(k) \\ \end{array}\right) \,, \end{eqnarray}$$

$\rho(\cdot)$

# 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

$\pm \frac{z_{1-\alpha/2}}{\sqrt{n}}$$z_{1-\alpha/2}$$1-\alpha/2$

$\pm z_{1-\alpha/2}\sqrt{\frac{1}{n}\left(1 + 2\sum_{i=1}^k \rho(i)^2\right)}$

„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ątkowe$ts_{t-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)

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)

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

---

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 czas$ts_{t-4}$, trzymamy $ts_t$, $ts_{t-1}$, $ts_{t-2}$ i $ts_{t-3}$, jak również $ts_{t-4}$i cofamy się $ts_t \sim ts_{t-1} + ts_{t-2} + ts_{t-3}+ts_{t-4}$ poprzez pochodzenie i zachowując jedynie współczynnik dla$ts_{t-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:

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 .

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

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]