Formule ACF et PACF

18

Je veux créer un code pour tracer ACF et PACF à partir de données de séries chronologiques. Tout comme ce graphique généré à partir de minitab (ci-dessous).

Tracé ACF

Tracé PACF

J'ai essayé de rechercher la formule, mais je ne la comprends toujours pas bien. Pourriez-vous me dire la formule et comment l'utiliser, s'il vous plaît? Quelle est la ligne rouge horizontale sur le tracé ACF et PACF ci-dessus? Quelle est la formule?

Je vous remercie,

Surya Dewangga
la source
1
@javlacalle La formule que vous fournissez est-elle correcte? Cela ne fonctionnerait pas si n t = 1 ( y t - ˉ y ) < 0
ρ(k)=1nkt=k+1n(yty¯)(ytky¯)1nt=1n(yty¯)1nkt=k+1n(ytky¯),
droite? Doit-il être comme suit? $$ \ 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(yty¯)<0and/ort=k+1n(yt-k-y¯)<0
conighion
@conighion Vous avez raison, merci. Je ne l'ai pas vu avant. Je l'ai corrigé.
javlacalle

Réponses:

33

Autocorrélations

La corrélation entre deux variables y1,y2 est définie comme:

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

où E est l'opérateur d'espérance, μ1 et μ2 sont les moyennes respectivement pour y1 et y2 et σ1,σ2 sont leurs écarts-types.

Dans le contexte d'une seule variable, c'est-à - dire l' auto- corrélation, y1 est la série d'origine et y2 est une version retardée. Lors de la définition ci - dessus, l' échantillon autocorrélations d'ordre k=0,1,2,...peut être obtenue en calculant l'expression suivante avec la série observée yt , t=1,2,...,n :

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

y¯ est la moyenne de l'échantillon des données.

Autocorrélations partielles

Les autocorrélations partielles mesurent la dépendance linéaire d'une variable après avoir supprimé l'effet des autres variables affectant les deux variables. Par exemple, l'autocorrélation partielle de l'ordre mesure l'effet (dépendance linéaire) de yt2 sur yt après avoir supprimé l'effet de yt1 sur yt et yt2 .

Chaque autocorrélation partielle peut être obtenue sous la forme d'une série de régressions de la forme:

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

y~t est la série d'origine moins la moyenne de l'échantillon, yty¯ . L'estimation de ϕ22 donnera la valeur de l'autocorrélation partielle d'ordre 2. En étendant la régression avec k décalages supplémentaires, l'estimation du dernier terme donnera l'autocorrélation partielle d'ordre k .

Une autre façon de calculer l'échantillon d'autocorrélations partielles consiste à résoudre le système suivant pour chaque ordre k :

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

ρ() sont les autocorrélations de l'échantillon. Cette correspondance entre les autocorrélations d'échantillons et les autocorrélations partielles est connue sous le nom de récursion de Durbin-Levinson . Cette approche est relativement facile à mettre en œuvre à titre d'illustration. Par exemple, dans le logiciel R, on peut obtenir l'autocorrélation partielle d'ordre 5 comme suit:

# 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

Bandes de confiance

Les bandes de confiance peuvent être calculées comme la valeur des autocorrélations de l'échantillon ±z1α/2n , oùz1α/2est le quantile1α/2dans la distribution gaussienne, par exemple 1,96 pour les bandes de confiance à 95%.

Parfois, des bandes de confiance qui augmentent à mesure que l'ordre augmente sont utilisées. Dans ce cas, les bandes peuvent être définies comme ±z1α/21n(1+2i=1kρ(i)2) .

javlacalle
la source
1
(+1) Pourquoi les deux bandes de confiance différentes?
Scortchi - Réintégrer Monica
2
@Scortchi Les bandes constantes sont utilisées lors des tests d'indépendance, tandis que les bandes croissantes sont parfois utilisées lors de l'identification d'un modèle ARIMA.
javlacalle
1
Les deux méthodes de calcul des bandes de confiance sont expliquées un peu plus en détail ici .
Scortchi - Réintégrer Monica
Explication parfaite!
Jan Rothkegel
1
ρ(k)
9

"Je veux créer un code pour tracer ACF et PACF à partir de données de séries chronologiques".

Bien que le PO soit un peu vague, il peut être plus ciblé sur une formulation de codage de style "recette" qu'une formulation de modèle d'algèbre linéaire.


L' ACFtt1tst333

Exemple:

Nous allons concocter une série temporelle avec un motif sinusoïdal cyclique superposé à une ligne de tendance et du bruit, et tracer l'ACF généré par R. J'ai obtenu cet exemple dans une publication en ligne de Christoph Scherber, et je viens d'y ajouter du bruit:

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

entrez la description de l'image ici

Normalement, nous devrions tester la stationnarité des données (ou simplement regarder le graphique ci-dessus), mais nous savons qu'il y a une tendance, alors sautons cette partie et passons directement à l'étape de décroissance:

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

entrez la description de l'image ici

Nous sommes maintenant prêts à reprendre cette série chronologique en générant d'abord l'ACF avec la acf()fonction dans R, puis en comparant les résultats à la boucle de fortune que j'ai composée:

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

entrez la description de l'image ici


tst4tsttst1tst2tst3tst4tsttst1+tst2+tst3+tst4 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.
}

Et enfin, tracer à nouveau côte à côte, les calculs générés par R et manuels:

entrez la description de l'image ici

Que l'idée est correcte, à côté des problèmes de calcul probables, peut être comparé PACFà pacf(st.y, plot = F).


code ici .

Antoni Parellada
la source
1

et les bandes de confiance vous aident à déterminer si un niveau peut être considéré comme du bruit uniquement (car environ 95% des temps seront dans les bandes).

user120580
la source
Bienvenue sur CV, vous voudrez peut-être envisager d'ajouter des informations plus détaillées sur la façon dont OP procéderait spécifiquement. Peut-être aussi ajouter quelques informations sur ce que représente chaque ligne?
Repmat
1

Voici un code python pour calculer 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
la source
Le formatage du code Hmmm était mauvais:
Sada