<jupyter_start><jupyter_text>______Copyright <NAME>For more information, visit us at www.pieriandata.com  ACF and PACF Autocorrelation Function / Partial Autocorrelation FunctionBefore we can investigate autoregression as a modeling tool, we need to look at covariance and correlation as they relate to lagged (shifted) samples of a time series. Goals * Be able to create ACF and PACF charts * Create these charts for multiple times series, one with seasonality and another without * Be able to calculate Orders PQD terms for ARIMA off these charts (highlight where they cross the x axis) Related Functions:stattools.acovf(x[, unbiased, demean, fft, …])&nbsp;Autocovariance for 1Dstattools.acf(x[, unbiased, nlags, qstat, …])&nbsp;&nbsp;Autocorrelation function for 1d arraysstattools.pacf(x[, nlags, method, alpha])&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Partial autocorrelation estimatedstattools.pacf_yw(x[, nlags, method])&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Partial autocorrelation estimated with non-recursive yule_walkerstattools.pacf_ols(x[, nlags])&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Calculate partial autocorrelations   Related Plot Methods:tsaplots.plot_acf(x)&nbsp;&nbsp;&nbsp;Plot the autocorrelation functiontsaplots.plot_pacf(x)&nbsp;&nbsp;Plot the partial autocorrelation functionFor Further Reading:Wikipedia:&nbsp;&nbsp;AutocovarianceForecasting: Principles and Practice&nbsp;&nbsp;AutocorrelationNIST Statistics Handbook&nbsp;&nbsp;Partial Autocorrelation Plot  Perform standard imports and load datasets<jupyter_code>import pandas as pd
import numpy as np
%matplotlib inline
import statsmodels.api as sm

# Load a non-stationary dataset
df1 = pd.read_csv('../Data/airline_passengers.csv',index_col='Month',parse_dates=True)
df1.index.freq = 'MS'

# Load a stationary dataset
df2 = pd.read_csv('../Data/DailyTotalFemaleBirths.csv',index_col='Date',parse_dates=True)
df2.index.freq = 'D'
# Import the models we'll be using in this section
from statsmodels.tsa.stattools import acovf,acf,pacf,pacf_yw,pacf_ols<jupyter_output><empty_output><jupyter_text>Ignore harmless warningsA quick note before we get started. Many of the models used in this and upcoming sections are likely to raise harmless warnings.For instance, the unbiased partial autocorrelation pacf_yw() performed below may raise a RuntimeWarning: invalid value encountered in sqrt.We don't really need to be concerned with this, and we can avoid it with the following code:<jupyter_code>import warnings
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Autocovariance for 1DIn a deterministic process, like $y=sin(x)$, we always know the value of $y$ for a given value of $x$. However, in a stochastic process there is always some randomness that prevents us from knowing the value of $y$. Instead, we analyze the past (or lagged) behavior of the system to derive a probabilistic estimate for $\hat{y}$.One useful descriptor is covariance. When talking about dependent and independent $x$ and $y$ variables, covariance describes how the variance in $x$ relates to the variance in $y$. Here the size of the covariance isn't really important, as $x$ and $y$ may have very different scales. However, if the covariance is positive it means that $x$ and $y$ are changing in the same direction, and may be related.With a time series, $x$ is a fixed interval. Here we want to look at the variance of $y_t$ against lagged or shifted values of $y_{t+k}$For a stationary time series, the autocovariance function for $\gamma$ (gamma) is given as:${\displaystyle {\gamma}_{XX}(t_{1},t_{2})=\operatorname {Cov} \left[X_{t_{1}},X_{t_{2}}\right]=\operatorname {E} [(X_{t_{1}}-\mu _{t_{1}})(X_{t_{2}}-\mu _{t_{2}})]}$We can calculate a specific $\gamma_k$ with:${\displaystyle \gamma_k = \frac 1 n \sum\limits_{t=1}^{n-k} (y_t - \bar{y})(y_{t+k}-\bar{y})}$ A NOTE ON FORMULA CONVENTIONS: Different texts employ different symbol conventions. For example, in the above autocovariance formula we use $k$ to represent the amount of lag or shift. Some texts use $h$ instead.  Autocovariance Example:Say we have a time series with five observations: {13, 5, 11, 12, 9}.We can quickly see that $n = 5$, the mean $\bar{y} = 10$, and we'll see that the variance $\sigma^2 = 8$.The following calculations give us our covariance values:$\gamma_0 = \frac {(13-10)(13-10)+(5-10)(5-10)+(11-10)(11-10)+(12-10)(12-10)+(9-10)(9-10)} 5 = \frac {40} 5 = 8.0 \\\gamma_1 = \frac {(13-10)(5-10)+(5-10)(11-10)+(11-10)(12-10)+(12-10)(9-10)} 5 = \frac {-20} 5 = -4.0 \\\gamma_2 = \frac {(13-10)(11-10)+(5-10)(12-10)+(11-10)(9-10)} 5 = \frac {-8} 5 = -1.6 \\\gamma_3 = \frac {(13-10)(12-10)+(5-10)(9-10)} 5 = \frac {11} 5 = 2.2 \\\gamma_4 = \frac {(13-10)(9-10)} 5 = \frac {-3} 5 = -0.6$Note that $\gamma_0$ is just the population variance $\sigma^2$Let's see if statsmodels gives us the same results! For this we'll create a fake dataset:<jupyter_code>df = pd.DataFrame({'a':[13, 5, 11, 12, 9]})
arr = acovf(df['a'])
arr<jupyter_output><empty_output><jupyter_text>Unbiased AutocovarianceNote that the number of terms in the calculations above are decreasing.Statsmodels can return an "unbiased" autocovariance where instead of dividing by $n$ we divide by $n-k$.$\gamma_0 = \frac {(13-10)(13-10)+(5-10)(5-10)+(11-10)(11-10)+(12-10)(12-10)+(9-10)(9-10)} {5-0} = \frac {40} 5 = 8.0 \\\gamma_1 = \frac {(13-10)(5-10)+(5-10)(11-10)+(11-10)(12-10)+(12-10)(9-10)} {5-1} = \frac {-20} 4 = -5.0 \\\gamma_2 = \frac {(13-10)(11-10)+(5-10)(12-10)+(11-10)(9-10)} {5-2} = \frac {-8} 3 = -2.67 \\\gamma_3 = \frac {(13-10)(12-10)+(5-10)(9-10)} {5-3} = \frac {11} 2 = 5.5 \\\gamma_4 = \frac {(13-10)(9-10)} {5-4} = \frac {-3} 1 = -3.0$<jupyter_code>arr2 = acovf(df['a'],unbiased=True)
arr2<jupyter_output><empty_output><jupyter_text>Autocorrelation for 1DThe correlation $\rho$ (rho) between two variables $y_1,y_2$ is given as: $\rho = \frac {\operatorname E[(y_1−\mu_1)(y_2−\mu_2)]} {\sigma_{1}\sigma_{2}} = \frac {\operatorname {Cov} (y_1,y_2)} {\sigma_{1}\sigma_{2}}$,where $E$ is the expectation operator, $\mu_{1},\sigma_{1}$ and $\mu_{2},\sigma_{2}$ are the means and standard deviations of $y_1$ and $y_2$.When working with a single variable (i.e. autocorrelation) we would consider $y_1$ to be the original series and $y_2$ a lagged version of it. Note that with autocorrelation we work with $\bar y$, that is, the full population mean, and not the means of the reduced set of lagged factors (see note below).Thus, the formula for $\rho_k$ for a time series at lag $k$ is:${\displaystyle \rho_k = \frac {\sum\limits_{t=1}^{n-k} (y_t - \bar{y})(y_{t+k}-\bar{y})} {\sum\limits_{t=1}^{n} (y_t - \bar{y})^2}}$This can be written in terms of the covariance constant $\gamma_k$ as:${\displaystyle \rho_k = \frac {\gamma_k n} {\gamma_0 n} = \frac {\gamma_k} {\sigma^2}}$For example,$\rho_4 = \frac {\gamma_4} {\sigma^2} = \frac{-0.6} {8} = -0.075$Note that ACF values are bound by -1 and 1. That is, ${\displaystyle -1 \leq \rho_k \leq 1}$<jupyter_code>arr3 = acf(df['a'])
arr3<jupyter_output><empty_output><jupyter_text>Partial AutocorrelationPartial autocorrelations measure the linear dependence of one variable after removing the effect of other variable(s) that affect both variables. That is, the partial autocorrelation at lag $k$ is the autocorrelation between $y_t$ and $y_{t+k}$ that is not accounted for by lags $1$ through $k−1$.A common method employs the non-recursive Yule-Walker Equations:$\phi_0 = 1\\\phi_1 = \rho_1 = -0.50\\\phi_2 = \frac {\rho_2 - {\rho_1}^2} {1-{\rho_1}^2} = \frac {(-0.20) - {(-0.50)}^2} {1-{(-0.50)}^2}= \frac {-0.45} {0.75} = -0.60$ As $k$ increases, we can solve for $\phi_k$ using matrix algebra and the Levinson–Durbin recursion algorithm which maps the sample autocorrelations $\rho$ to a Toeplitz diagonal-constant matrix. The full solution is beyond the scope of this course, but the setup is as follows:$\displaystyle \begin{pmatrix}\rho_0&\rho_1&\cdots &\rho_{k-1}\\\rho_1&\rho_0&\cdots &\rho_{k-2}\\\vdots &\vdots &\ddots &\vdots \\\rho_{k-1}&\rho_{k-2}&\cdots &\rho_0\\\end{pmatrix}\quad \begin{pmatrix}\phi_{k1}\\\phi_{k2}\\\vdots\\\phi_{kk}\end{pmatrix}\mathbf = \begin{pmatrix}\rho_1\\\rho_2\\\vdots\\\rho_k\end{pmatrix}$<jupyter_code>arr4 = pacf_yw(df['a'],nlags=4,method='mle')
arr4<jupyter_output><empty_output><jupyter_text>NOTE: We passed in method='mle' above in order to use biased ACF coefficients. "mle" stands for "maximum likelihood estimation". Alternatively we can pass method='unbiased' (the statsmodels default):<jupyter_code>arr5 = pacf_yw(df['a'],nlags=4,method='unbiased')
arr5<jupyter_output><empty_output><jupyter_text>Partial Autocorrelation with OLSThis provides partial autocorrelations with ordinary least squares (OLS) estimates for each lag instead of Yule-Walker.<jupyter_code>arr6 = pacf_ols(df['a'],nlags=4)
arr6<jupyter_output><empty_output><jupyter_text>PlottingThe arrays returned by .acf(df) and .pacf_yw(df) show the magnitude of the autocorrelation for a given $y$ at time $t$. Before we look at plotting arrays, let's look at the data itself for evidence of autocorrelation.Pandas has a built-in plotting function that plots increasing $y_t$ values on the horizontal axis against lagged versions of the values $y_{t+1}$ on the vertical axis. If a dataset is non-stationary with an upward trend, then neighboring values should trend in the same way. Let's look at the Airline Passengers dataset first.<jupyter_code>from pandas.plotting import lag_plot

lag_plot(df1['Thousands of Passengers']);<jupyter_output><empty_output><jupyter_text>Visually this shows evidence of a very strong autocorrelation; as $y_t$ values increase, nearby (lagged) values also increase.Now let's look at the stationary Daily Total Female Births dataset:<jupyter_code>lag_plot(df2['Births']);<jupyter_output><empty_output><jupyter_text>As expected, there is little evidence of autocorrelation here.  ACF PlotsPlotting the magnitude of the autocorrelations over the first few (20-40) lags can say a lot about a time series.For example, consider the stationary Daily Total Female Births dataset:<jupyter_code>from statsmodels.graphics.tsaplots import plot_acf,plot_pacf
# Let's look first at the ACF array. By default acf() returns 40 lags
acf(df2['Births'])
# Now let's plot the autocorrelation at different lags
title = 'Autocorrelation: Daily Female Births'
lags = 40
plot_acf(df2,title=title,lags=lags);<jupyter_output><empty_output><jupyter_text>This is a typical ACF plot for stationary data, with lags on the horizontal axis and correlations on the vertical axis. The first value $y_0$ is always 1. A sharp dropoff indicates that there is no AR component in the ARIMA model.Next we'll look at non-stationary data with the Airline Passengers dataset:<jupyter_code>acf(df1['Thousands of Passengers'])
title = 'Autocorrelation: Airline Passengers'
lags = 40
plot_acf(df1,title=title,lags=lags);<jupyter_output><empty_output><jupyter_text>This plot indicates non-stationary data, as there are a large number of lags before ACF values drop off.  PACF PlotsPartial autocorrelations work best with stationary data. Let's look first at Daily Total Female Births:<jupyter_code>title='Partial Autocorrelation: Daily Female Births'
lags=40
plot_pacf(df2,title=title,lags=lags);<jupyter_output><empty_output><jupyter_text>To make the Airline Passengers data stationary, we'll first apply differencing:<jupyter_code>from statsmodels.tsa.statespace.tools import diff

df1['d1'] = diff(df1['Thousands of Passengers'],k_diff=1)
df1['d1'].plot(figsize=(12,5));
title='PACF: Airline Passengers First Difference'
lags=40
plot_pacf(df1['d1'].dropna(),title=title,lags=np.arange(lags));  # be sure to add .dropna() here!<jupyter_output><empty_output><jupyter_text>Plot ResizingIn case you want to display the full autocorrelation plot, it helps to increase the figure size using matplotlib:<jupyter_code>import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12,5))

plot_acf(df1,ax=ax);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13729051432345413, 'fraction_numerical': 0.03880128787253364, 'mean_word_length': 5.685430463576159, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '249926', 'license': 'MIT', 'n_tokens_mistral': 4387, 'n_tokens_neox': 3983, 'n_words': 1310, 'path': 'tsa/jose/TSA_COURSE_NOTEBOOKS/06-General-Forecasting-Models/01-ACF-and-PACF.ipynb', 'repo_name': 'juspreet51/ml_templates', 'size': '12119'}