python通过自相关函数确定周期 python自相关系数

介绍

代码、数据全部免费，都放在我的gitee仓库里面https://gitee.com/yuanzhoulvpi/time_series，想要使用的可以直接到我这个仓库下载。

本文是继第二篇，详细介绍了时间序列的两个重要的相关性系数：acf和pacf，本文将围绕着acf、pacf怎么算、怎么画来具体展开。

这也是你见到的第一篇详细介绍如何用python复现acf、pacf～

本文内容比较多【粗略统计：所有文本接近8000字；代码超过300行】。

具体的jupyter notebook结构如下：

1. 时间序列中最基本的概念；
2. 什么叫lag（滞后）；
3. 时间序列的公式；
4. 统计基础知识；
5. 自相关系数（公式、复现、调包、可视化）；
6. 偏自相关系数（公式、复现、调包、可视化）；

1.时间序列中最基本的概念

很多人看不懂时间序列的公式、很多人看不懂时间序列的数据、代码、数学函数。本文就是将数据、公式、代码、数学函数串联起来，让你不仅仅了解抽象的数学公式，也能将抽象的内容联系到具体的数据和python代码中，并且也不只是单纯的调用python包，而是自己手写算法。

下面的代码是近20年城镇就业人员数据(万人)，数据来自中国统计局。

import numpy as np
import pandas as pd
import datetime
import matplotlib.pyplot as plt

data = pd.read_csv("../datasets/近20年城镇就业人员数据(万人).csv", index_col=['year'])
data.head(4)

# 可视化
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(data.index, data['value'], 'o-')
ax.set_title(u"Data on urban employees in the past 20 years")
ax.set_ylabel("Ten thousand people")
ax.set_xlabel("year")
ax.set_xticks(data.index)
ax.set_xticklabels(data.index, rotation=45)

2.什么叫lag（滞后）

经常在时间序列里面可以看到lag这个东西，大白话就是叫滞后；举个例子来说：

1. lag1代表滞后数据1项，也就是2002年的数据滞后1年得到的就是2001年的数据；2003年的数据滞后1年就是2002年的数据；
2. lag2代表滞后数据2项，也就是2003年的数据滞后2年得到2001年的数据；2004年的数据滞后2年就是2002年的数据；
3. 依次类推，滞后数据就是这么计算的。

下面使用python对上面的【近20年城镇就业人员数据(万人)】数据做个滞后处理。

# 使用pandas的shift可以做lag
data1 = data.copy()
# data1['value_lag1'] = data1['value'].shift(1)
# data1['value_lag2'] =  data1['value'].shift(2)
for i in [1, 2, 3, 4, 5]:
    data1[f"value_lag{i}"] = data1['value'].shift(i)
data1.head(10)

数据解读：

1. 因为2001年之前是没有数据了，所以2001年滞后一项得到的数据是空数据；
2. 因为2002年之前是没有数据，所以2002年滞后两项也是没有数据的。
3. 依次类推。

3.时间序列的公式

经常看到时间序列的公式写的乱七八糟的，这是一堆公式，那是一堆公式，活活把人绕晕。现在为了我们说明需要，我们定义几个简单的公式：

假设我们的时间序列是value = [4,3,2,10,3]。我们把这个时间序列定义为 ；那么这个时间序列里面每一个元素对应的符号依次是：、、等。

如果做了lag1，得到的时间序列为value lag1=[nan, 4,3,2,10]，记这个新的时间序列为。

如果做了lag2，得到的时间序列为value lag2=[nan,nan, 4,3,2]，记这个新的时间序列为。

以此类推，就是一个简单的时间序列公式了。

4.统计基础知识

1. 均值:
2. 方差：
3. 协方差：
4. 相关系数:

5.自相关系数

5.1自相关基础计算

在了解上面的【统计基础知识】后，再理解自相关性基本上就不难了，自相关是变量与自身在不同时间滞后下的相关性。对于延迟阶的时间序列。自协方差为，自相关性系数为

5.2调用包计算acf数值

# 使用statsmodels包的acf函数
from statsmodels.tsa.stattools import acf

5.3自己根据公式写函数

def cal_acf(x, nlags):
    """
    按照公式自己写个acf函数
    只是用来帮助理解，不建议用在别的环境中
    """
    x = np.array(x)
    mean_x = np.mean(x)
    length_x = x.shape[0]
    c_0 = np.mean((x-mean_x) **2)
    c_k = np.sum((x[:(length_x-nlags)] - mean_x) * (x[nlags:] - mean_x)) / length_x
    r_k = c_k / c_0
    return r_k

5.4两种结果对比

# 结果汇总
pd.DataFrame({'index':np.arange(11),
             'value_by_myself':[cal_acf(x=data['value'], nlags=i) for i in range(11)],
             'value_by_statsmodels':acf(data.value,nlags=10)})

结果如下：

5.5自相关系数画图

画图有两种方法：

1. 直接使用statsmodels的plot_acf;
2. 使用acf函数来计算得到数值，得到数值后再进行可视化，这种一般便于自己定制，也方便自己做更多的细节调整之类的。

# 自己画一个
# 自己计算
acf_value, acf_interval, _, _ = acf(data.value,nlags=14,qstat=True,alpha=0.05, fft=False)

xlabel = np.arange(start=0, stop=acf_value.shape[0], dtype='float')

fig, ax = plt.subplots(nrows=2, figsize=(10,6), sharex=True, dpi=100)
ax[0].hlines(y=0, xmin=np.min(xlabel)-2, xmax=np.max(xlabel)+2, colors='red', linewidth=0.5)
ax[0].scatter(x=xlabel, y=acf_value, c='red')
ax[0].vlines(x = xlabel, ymin=0, ymax=acf_value, colors='red', linewidth=0.5)
xlabel[1] -= 0.5
xlabel[-1] += 0.5
ax[0].fill_between(x=xlabel[1:], y1=acf_interval[1:,0] - acf_value[1:], y2=acf_interval[1:, 1]-acf_value[1:], alpha=0.25, linewidth=0, color='red')
ax[0].set_title("acf plot by myself")


# 使用别人写好的函数

from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data['value'], ax=ax[1])
ax[1].set_title("acf plot by statsmodels")
ax[1].set_xlim(-1, np.max(xlabel)+1)

6.偏自相关系数

偏自相关系数可以被看成一种条件相关，两个变量之间有相关性是基于一系列的其他的条件变量。想象一个这样的问题：

我们有一个Y变量，这个Y变量和、相关，还有一个变量。那我要知道和是否相关应该怎么计算？

肯定是应该要剔除掉、对、对的影响。那么计算和的相关性就是为：

那么对于一个时间序列来说， 和之间的偏自相关系数：就是基于 、、、…… 对 和的影响，计算 和之间相关系数。

6.1偏自相关计算公式：

1. 当h=1的时候，我们考虑的就是 和之间的偏自相关系数,这个时候中间是没有别的条件序列的，所以偏自相关系数就是等于自相关系数。
2. 当h=2的时候，我们考虑的就是 和之间的偏自相关系数，这个时候中间的条件序列就是，所以偏自相关系数就是为:

.

1. 当h=3的时候，我们考虑的就是 和之间的偏自相关系数，这个时候中间的条件序列就是、，所以偏自相关系数就是为:

，，，，.

1. 依次类推。

6.2偏自相关计算：

1. 使用statsmodels的pacf函数，这个函数可以选择不同的计算方法，有使用ols方法的，也有使用Yule-Walker方法的。
2. 除了方法1调用包，我们还可以使用公式方法，分别使用ols方法、使用Yule-Walker方法来计算pacf的系数。实际上算出来的结果和statsmodels的pacf函数算出来的结果是一模一样的。

在下面的代码块中，我们使用statsmodels包里面的pacf、pacf_ols、pacf_yw函数计算偏自相关系数，另外我们自己通过ols理论和Yule-Walker理论计算偏自相关系数。

# 使用statsmodels包的pacf函数
from statsmodels.tsa.stattools import pacf,pacf_ols,pacf_yw


# 自己实现用ols求解pacf的函数
from numpy.dual import lstsq
from statsmodels.tools import add_constant
from statsmodels.tsa.tsatools import lagmat

def cal_my_pacf_ols(x, nlags=5):
    """
    自己实现pacf，原理使用的就是ols(最小二乘法)
    :param x:
    :param nlags:
    :return:
    """
    pacf = np.empty(nlags+1) * 0
    pacf[0] = 1.0

    xlags, x0 = lagmat(x, nlags, original="sep")
    xlags = add_constant(xlags)

    for k in range(1, nlags+1):
        params = lstsq(xlags[k:, :(k+1)], x0[k:], rcond=None)[0]
        pacf[k] = params[-1]

    return pacf




def cal_my_yule_walker(x, nlags=5):
    """
    自己实现yule_walker理论
    :param x:
    :param nlags:
    :return:
    """
    x = np.array(x, dtype=np.float64)
    x -= x.mean()
    n = x.shape[0]

    r = np.zeros(shape=nlags+1, dtype=np.float64)
    r[0] = (x ** 2).sum()/n

    for k in range(1, nlags+1):
        r[k] = (x[0:-k] * x[k:]).sum() / (n-k*1)

    from scipy.linalg import toeplitz
    R = toeplitz(c=r[:-1])
    result = np.linalg.solve(R, r[1:])
    return result

def cal_my_pacf_yw(x, nlags=5):
    """
    自己通过yule_walker方法求出pacf的值
    :param x:
    :param nlags:
    :return:
    """
    pacf = np.empty(nlags+1) * 0
    pacf[0] = 1.0
    for k in range(1, nlags+1):
        pacf[k] = cal_my_yule_walker(x,nlags=k)[-1]

    return pacf

cal_my_pacf_yw(x=data.values)

# 结果比较

pd.DataFrame({'pacf':pacf(data.value,nlags=5, method='ols'),
              'pacf_ols':pacf_ols(data.value,nlags=5),
              'pacf_yw':pacf_yw(data.value,nlags=5),
              'pacf_ols_bymyself':cal_my_pacf_ols(x=data.values),
              'pacf_yw_bymyself':cal_my_pacf_yw(x=data.values)})

效果对比：

6.3偏自相关系数画图

画图有两种方法：

1. 直接使用statsmodels的plot_pacf;
2. 使用pacf函数来计算数据，得到数据后再进行可视化，这种一般便于自己定制，也方便自己做更多的细节调整之类的。

#%%
sunspotarea = pd.read_csv("../datasets/sunspotarea.csv",parse_dates=['date'])


# 自己画一个
# 自己计算
pacf_value, pacf_interval = pacf(sunspotarea.value,nlags=15,alpha=0.05)

xlabel = np.arange(start=0, stop=pacf_value.shape[0], dtype='float')

fig, ax = plt.subplots(nrows=2, figsize=(10,6), sharex=True, dpi=100)
ax[0].hlines(y=0, xmin=np.min(xlabel)-2, xmax=np.max(xlabel)+2, colors='red', linewidth=0.5)
ax[0].scatter(x=xlabel, y=pacf_value, c='red')
ax[0].vlines(x = xlabel, ymin=0, ymax=pacf_value, colors='red', linewidth=0.5)
xlabel[1] -= 0.5
xlabel[-1] += 0.5
ax[0].fill_between(x=xlabel[1:], y1=pacf_interval[1:,0] - pacf_value[1:], y2=pacf_interval[1:, 1]-pacf_value[1:], alpha=0.25, linewidth=0, color='red')
ax[0].set_title("pacf plot by myself")


# 使用别人写好的函数


from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(sunspotarea.value, ax=ax[1], lags=15)
ax[1].set_title("pacf plot by statsmodels")
ax[1].set_xlim(-1, np.max(xlabel)+1)

写到后面

1. 本文主要强调的是如何手写代码复现acf和pacf计算以及可视化。个人建议多多看看Yule-Walker方法计算pacf的代码。
2. 后续会出一个教程，介绍Yule-Walker方程的意义、复现、代码等。
3. 这期之所以会隔这么多天，是在计算pacf的时候，卡住了。后来看了很多资料，终于搞懂了原理和复现出代码。