计量经济与时间序列_ACF自相关与PACF偏自相关算法解析(Python,TB(交易开拓者))

1   在时间序列中ACF图和PACF图是非常重要的两个概念,如果运用时间序列做建模、交易或者预测的话。这两个概念是必须的。

2   ACF和PACF分别为:自相关函数(系数)和偏自相关函数(系数)。

3   在许多软件中比如Eviews分析软件可以调出某一个序列的ACF图和PACF图,如下:

  

  3.1   有时候这张图是横躺着的,不过这个不重要,反正一侧为小于0的负值范围,一侧为大于0的正值范围,均值(准确的说是坐标y轴为0,有些横着的图,会把x轴和y轴表示出来,值都在x轴上下附近呈现出来)。

  3.2   红色框框部分就是ACF图,青色框框部分就是PACF图,其中对应左边的Autocorrelation就是英文单词自相关的全称;Partial Correlation就是英文单词偏自相关的全称。

  3.3   我们要计算的就是这两列数值。

  3.4   其中紫色箭头标注出来的是指的2倍标准误范围,后面可以用对应的数值是否超过范围来判断截尾、拖尾等信息,进而判断采用哪种模型。

  3.5   这里特别注明一下在默认状态下,这两根线是如何计算的:

    在大样本下(T很大的时候,这里T指的是样本的个数;其实准确的说样本符合均值为0的正太分布)。因此这里的对于ACF或PACF属于一种分位点检验,这个东西在很多baidu可以找到一个正太分布图,然后左右画线会得到99%,90%...的分位点,这里的这两根线就是指的这个。我们这里要做的是双侧对称检验,所以上下两根线分布式0±分位点的值。分位点=2倍×sqrt开方(1/T),这里的T指的是样本的个数,样本个数指的是原始的那个样本的个数,不是ACF或PACF计算完的样本个数。如果某一个值大于2倍标准误,也就是说大于正态分布左/右的95%分位点,于是,在拒绝域,则拒绝0假设(也就是拒绝均值为0的假设)

    例如:样本共10个,ACF或PACF计算完毕后,他们的数量都为9个,其分为点为:2×sqrt(1/10) = 0.6324555320...,双侧检验边沿值为:(0-0.6324555320,0+0.6324555320) = [-0.6324555320,+0.6324555320](这就是两根虚线的值)。    (3.5.1)

       另外我们还通过另外一个公式Sρk ≈ sqrt(1/n * (1+2*ρ12 + ... + 2*ρk-12))来计算acf的双侧检验值。具体在这篇博文中得出。http://www.cnblogs.com/noah0532/p/8453959.html

4    ACF算法的实现(python

  4.1   首先:给定一组数据,这里的是线性数据,不用太多10个元素,一组的时间序列。

1 TimeSeries = [13, 8, 15, 4, 4, 12, 11, 7, 14, 12]  # 列表形式
2 print(TimeSeries)
3 # 显示结果:[13, 8, 15, 4, 4, 12, 11, 7, 14, 12]

  4.2   算法表达式。

  

  4.3   算法解析:

  计算ρ1:

  给定一阶滞后序列:其中滞后用L来表示,具体表示方法见链接:http://www.cnblogs.com/noah0532/p/8449991.html

1 LZt = TimeSeries[:-1:]
2 print(LZt)
3 # 显示结果:[13, 8, 15, 4, 4, 12, 11, 7, 14]

  对应的这两个序列分别为:

  [13, 8, 15, 4, 4, 12, 11, 7, 14, 12]
  [13, 8, 15, 4, 4, 12, 11, 7, 14]

  我们在把原始序列与一阶滞后序列进行对齐:

1 Zt = TimeSeries[1::]
2 print(Zt)
3 # 显示结果:[8, 15, 4, 4, 12, 11, 7, 14, 12]

  这样两个序列变为:

  [ 8, 15, 4, 4, 12, 11, 7, 14, 12]
  [13, 8, 15, 4, 4, 12, 11, 7, 14]  

  对这两个序列进行计算就会得出ρ1的结果(也就是ACF的第一个值)

  4.4   其实我们可以得到这样一个结果,就是这样不断的比较下去,计算每一个滞后,原始序列会少一个值来对应滞后的序列,然后分别来进行计算。我们把这个过程化简一下,代码如下(以如上10个序列为例,最后会得出9对)

 1 TimeSeries = [13, 8, 15, 4, 4, 12, 11, 7, 14, 12]  # 列表形式
 2 Zt = []
 3 LZt = []
 4 i = 1
 5 while i < len(TimeSeries):
 6     L = TimeSeries[i::]
 7     LL = TimeSeries[:-i:]
 8     Zt.append(L)
 9     LZt.append(LL)
10     i += 1
11 print(TimeSeries)
12 print(Zt)
13 print(LZt)
14 # 显示结果:
15 # [13, 8, 15, 4, 4, 12, 11, 7, 14, 12]
16 # [[8, 15, 4, 4, 12, 11, 7, 14, 12], [15, 4, 4, 12, 11, 7, 14, 12], [4, 4, 12, 11, 7, 14, 12], [4, 12, 11, 7, 14, 12], [12, 11, 7, 14, 12], [11, 7, 14, 12], [7, 14, 12], [14, 12], [12]]
17 # [[13, 8, 15, 4, 4, 12, 11, 7, 14], [13, 8, 15, 4, 4, 12, 11, 7], [13, 8, 15, 4, 4, 12, 11], [13, 8, 15, 4, 4, 12], [13, 8, 15, 4, 4], [13, 8, 15, 4], [13, 8, 15], [13, 8], [13]]

  4.5   通过这种遍历加减方式,如果原始列表为10个元素的话会得到9对ACF带计算序列,以此类推,这样做的目的就是为了把序列进行对齐计算。把这9对列出来看一下:

  [8, 15, 4, 4, 12, 11, 7, 14, 12]
  [13, 8, 15, 4, 4, 12, 11, 7, 14]

  [15, 4, 4, 12, 11, 7, 14, 12]
  [13, 8, 15, 4, 4, 12, 11, 7]

  [4, 4, 12, 11, 7, 14, 12]
  [13, 8, 15, 4, 4, 12, 11]

  [4, 12, 11, 7, 14, 12]
  [13, 8, 15, 4, 4, 12]

  [12, 11, 7, 14, 12]
  [13, 8, 15, 4, 4]

  [11, 7, 14, 12]
  [13, 8, 15, 4]

  [7, 14, 12]
  [13, 8, 15]

  [14, 12]
  [13, 8]

  [12]
  [13]

   4.6   根据算法公式,其分母为原始序列每一个元素距离平均值的平方。因此我们也可以把分母计算出来。

 1 TimeSeries = [13, 8, 15, 4, 4, 12, 11, 7, 14, 12]  # 列表形式
 2 Zt = []
 3 LZt = []
 4 sum = 0
 5 i = 1
 6 while i < len(TimeSeries):
 7     L = TimeSeries[i::]
 8     LL = TimeSeries[:-i:]
 9     sum = sum + TimeSeries[i - 1]
10     Zt.append(L)
11     LZt.append(LL)
12     i += 1
13 sum = sum + TimeSeries[-1]
14 avg = sum / len(TimeSeries)
15 print(TimeSeries)
16 print(Zt)
17 print(LZt)
18 print(avg)
19 # 显示结果:
20 # 10.0

  4.7   我们具备了这些所有的元素后开始计算这些所有ACF的ρ值,最终我们还是用列表来表示

 1 k = 0
 2 result_Deno = 0
 3 # 首先计算分母=分母都为通用
 4 while k < len(TimeSeries):
 5     result_Deno = result_Deno + pow((TimeSeries[k] - avg), 2)
 6     k += 1
 7 print(result_Deno)
 8 # 显示结果:144.0
 9 
10 # 然后计算分子
11 p = 0
12 q = 0
13 14 acf = []
15 while p < len(Zt):
16     # print(Zt[p])
17     # print(LZt[p])
18     q = 0
19     result_Mole = 0
20     while q < len(Zt[p]):
21         result_Mole = result_Mole + (Zt[p][q] - avg) * (LZt[p][q] - avg)
22         q += 1
23     acf.append(round(result_Mole / result_Deno, 3))  # 保留小数点后三位
24     p += 1
25 print(acf)
26 # 显示结果:
27 # [-0.188, -0.201, 0.181, -0.132, -0.326, 0.118, -0.049, 0.056, 0.042]

  4.7   结果验证:与用Eviews软件计算的结果一致,算法结束:

[-0.188, -0.201, 0.181, -0.132, -0.326, 0.118, -0.049, 0.056, 0.042]

5   PACF算法的实现(python

    5.1   PACF算法源于ACF算法,必须先算出ACF的值来,然后再计算PACF,PACF这里用的是一个递推形式的公式,计算每一个φ值。在书中有些φ值用φkk来表示,有时候也用一个大Pk来表示,其实道理都一样。

  5.2   递推式如下:

  5.3   这个递推式,紫色框内为每期值,下面的属于过度值,用于递推累计。

  5.4   这个递推式算法看起来无从下手比较麻烦。我们先分解一下:

  5.5   j和k的下标解析

    如果k = 1,那么 j = 1,2

    如果k = 2,那么 j = 1,2

    如果k = 3,那么 j = 1,2,3

    如果k = 4,那么 j = 1,2,3,4

    如果k = 5,那么 j = 1,2,3,4,5

    如果k = 6,那么 j = 1,2,3,4,5,6

    如果k = 7,那么 j = 1,2,3,4,5,6,7

    如果k = 8,那么 j = 1,2,3,4,5,6,7,8

    如果k = 9,那么 j = 1,2,3,4,5,6,7,8,9

    ... ....

    从下标可以看出这个递推关系,如果处于k的某种状态(数值),那么 j 将从当前状态递减到1。这也就是 j = 1,...,k的解析。

  5.6   从这个下标解析可以看出,这两部分的递推式,第一部分(紫色框框中)的是计算PACF当前值,而第二部分是计算每阶段的PACF值的缺省值(所需值)的计算。因此每期的PACF值的计算,需要递推第二个公式中的值,然后再去计算下一期PACF值。这个逻辑是这样

  5.7   从5.5的这个规律可以看出,j 是 k递减数,也就是当前值如果为k = 3的话,我们需要有3个第二个公式的值;如果k = 7的话,我们需要有7个第二个公式的值;这个逻辑是没有问题的,我们以一个多的为例子

    k = 8 通式如下:

    φ9,j = φ8j - φ998,9-j

    j = 1

    φ91 = φ81 - φ9988

    j = 2

    φ92 = φ82 - φ9987

    j = 3

    φ93 = φ83 - φ9986

    j = 4

    φ94 = φ84 - φ9985

    j = 5

    φ95 = φ85 - φ9984

    j = 6

    φ96 = φ86 - φ9983

    j = 7

    φ97 = φ87 - φ9982

    j = 8

    φ98 = φ88 - φ9981

   5.8   正如上面所示,我们会得到8个表达式。其每一期的表达式的值我们需要用列表来进行记录;这里的规律是φ99来自于第一个部分的公式,其φ8x和后面的φ8x是来自第二部分公式计算的上一期的值。这样这个第二部分的表达式我们就可以计算出来了。

   5.9   我们再来递推第一部分的那个公式。还是以上面的为例子。我们举一个小点儿的数。

    k = 2   通式如下

    φ33 = φ3 - sigma(φ2j3-j) / 1 - sigma(φ2jj)

    j = 1

    φ33 = φ3 - sigma(φ212) / 1 - sigma(φ211)

    j = 2

    φ33 = φ3 - sigma(φ221) / 1 - sigma(φ222)

    我们把sigma中的两部分进行合并得:

    φ33 = φ3 - (φ212 +  φ221) / 1 - (φ211 φ222)

   5.10   我们发现一个规律。如果要计算33的时候,其括号内的值φ的顺序不变,其ρ值为反向求积然后再求和。

   5.11   我们现在试着来求一下这段代码的写法,以上面计算出来的ACF值,接着来计算PACF值:

  ACF = [-0.188, -0.201, 0.181, -0.132, -0.326, 0.118, -0.049, 0.056, 0.042]

  ACF[0] = ρ1 = -0.188

  ACF[1] = ρ2 = -0.201

  ACF[2] = ρ4 = 0.181

  ACF[3] = ρ5 = -0.132

  ACF[4] = ρ6 = -0.326

  ACF[5] = ρ7 = 0.118

  ACF[6] = ρ8 = -0.049

  ACF[7] = ρ9 = 0.056

  ACF[8] = ρ9 = 0.042

  5.12   (具体过程稍显复杂,所以下面用这个颜色写,这样更清楚一些):

    (1):递推算法,在编写的时候可能会对初学者烧脑一些。重要的是对于每一种算法,我们首先要找到他们的规律,再不会算,可以找张纸先手工验算一下。

    (2):acf值与pacf值的长度都是一样,且第一个值都相等。φkk如果是两个整数,比如φ11,φ33,这就是我们要求的pacf对应每期的值。

    (3):正如上面所说的,都是整数的话是所求的值,那么不是整数的话,类似于φ41,φ42,φ43...等等这些,就是第二个递推式所需要计算的值。这是值是为了求下一组数值所用。

    (4):因此我们得出这样一个规律,除去第一个值(pacf和acf第一个值都相等)。后面每期值分两组来进行计算,然后一组循环结束,以此类推,知道pacf长度减-1都完成为止。

    (5):那么一组是什么样子的?比如第一组是:φ22,φ21;第二组:φ33,φ32,φ31;第三组:φ44,φ43,φ42,φ41...以此类推。一组等于一个大循环。

    (6):每组值分第一部分整值计算和第二部分非整值计算,每一组所需的φ值(看递推公式),是来自于上一组的值。因此我们计算完一组后,把所需要的值(代码中是用个过渡变量bridge来记录),送给下一组计算用,计算完毕删除,再赋新值,再给下一组。循环往复,直到结束。

    (7):那每组需要的值是什么样子,以计算到φ44,这个为例子,所需要的bridge值为[φ31,φ32,φ33]。用完了后,bridge更新为[φ41,φ42,φ43,φ44],供下一组φ55使用,其他每组都是这个原理。在这里引入了一个桥变量

    (8):由于精度可能略有差别,这样我们把之前计算的acf值先不要保留小数点后三位,计算完毕后,我们统一保留小数点后三位。

    (9):最后再把上边界和下边界的值算出来。这个公式上面有。

  5.13   具体代码如下:

  1 TimeSeries = [13, 8, 15, 4, 4, 12, 11, 7, 14, 12]  # 列表形式
  2 Zt = []
  3 LZt = []
  4 total = 0
  5 i = 1
  6 while i < len(TimeSeries):
  7     L = TimeSeries[i::]
  8     LL = TimeSeries[:-i:]
  9     total = total + TimeSeries[i-1]
 10     Zt.append(L)
 11     LZt.append(LL)
 12     i += 1
 13 total = total + TimeSeries[-1]
 14 avg = total / len(TimeSeries)
 15 
 16 k = 0
 17 result_Deno = 0
 18 # 首先计算分母=分母都为通用
 19 while k < len(TimeSeries):
 20     result_Deno = result_Deno + pow((TimeSeries[k] - avg), 2)
 21     k += 1
 22 # print(result_Deno)
 23 # 显示结果:144.0
 24 
 25 # 然后计算分子
 26 p = 0
 27 q = 0
 28 acf = []
 29 while p < len(Zt):
 30     # print(Zt[p])
 31     # print(LZt[p])
 32     q = 0
 33     result_Mole = 0
 34     while q < len(Zt[p]):
 35         result_Mole = result_Mole + (Zt[p][q] - avg) * (LZt[p][q] - avg)
 36         q += 1
 37     acf.append(result_Mole / result_Deno)  # 我们把前面计算的acf值先不要保留小数点后三位。
 38     # acf.append(result_Mole / result_Deno)
 39     p += 1
 40 # print(acf)
 41 
 42 # 初始化pacf
 43 pacf = []
 44 bridge = []
 45 bridge.append(acf[0])
 46 pacf.append(acf[0])   # 第一个φ11等于ρ1,再初始化pacf的第一个值。
 47 # 这里采用的是bridge这个循环变量。计算的逻辑是每次重新赋值一遍共下一轮计算所用。
 48 
 49 
 50 T = 0  # 初始化大循环的值为0,9个数值一共循环7次。原因是下面。
 51 while T < len(acf)-1:   # pacf的初始值已经计算过一个了,acf的长度和pacf的长度是一样的,因此,len(acf)此时要减去一个1
 52     # 每进行一轮循环,这些变量都重新初始化一遍
 53     Mole = 0 # 计算紫框中Σ累加的分子值
 54     Deno = 0 # 计算紫框中Σ累计的分母值
 55     cross = 0 # 每计算完一遍第一部分的值,暂时用cross这个变量先保存起来。
 56     UnderCross = [] # 每计算完第二步的值,用UnderCross的值先保存起来,因为有时候第二步的值不是一个,所以用列表形式。
 57 
 58     # 计算第一部分
 59     t = 0
 60     while t < len(bridge):    # 计算φ22,φ33,φ44,φ55,φ66,φ22,φ77,φ88...这些值
 61         Mole = Mole + (bridge[t] * acf[len(bridge)-1-t])
 62         Deno = Deno + (bridge[t] * acf[t])
 63         t += 1
 64     cross = (acf[t] - Mole) / (1 - Deno)  # acf[t] 正确,每次迭代完成,都是分子的第一个值。
 65     # print(cross)
 66     # print(t)
 67 
 68     # 计算第二部分
 69     p = 0
 70     while p < len(bridge):   # 计算像φ21,φ32,φ31...这样的值
 71         UnderCross.append(bridge[p] - cross * bridge[len(bridge)-1-p])
 72         p += 1
 73 
 74     bridge = []  # bridge使用完毕,初始化,再次赋值供下一次计算
 75     bridge.extend(UnderCross)   # 把计算完毕的Under值赋值到bridge
 76     bridge.append(cross)   # 把计算完毕的cross值存放到最后的位置
 77     pacf.append(cross)   # 把pacf值累计进行添加并记录
 78     # print(bridge)
 79     # print(UnderCross)
 80 
 81     T += 1
 82 # print(pacf)
 83 # 整个的一段循环完毕。
 84 
 85 
 86 # 最后我们重新遍历一遍,然后把acf和pacf都保留小数点后三位。
 87 AcfValue = []
 88 PacfValue = []
 89 lag = 0
 90 for sht in acf:
 91     lag = round(sht, 3)
 92     AcfValue.append(lag)
 93     lag = 0
 94 print("acf值(保留小数点后三位)为:", AcfValue)
 95 
 96 for sht in pacf:
 97     lag = round(sht, 3)
 98     PacfValue.append(lag)
 99     lag = 0
100 print("pacf值(保留小数点后三位)为:", PacfValue)
101 
102 
103 # 最后把bounds值算出来。
104 import math
105 bound = []
106 bound.append(-2*math.sqrt(1/len(TimeSeries)))
107 bound.append(2*math.sqrt(1/len(TimeSeries)))
108 print("上下边界值分别为:", bound)
109 
110 # 最终显示结果:
111 # acf值(保留小数点后三位)为: [-0.188, -0.201, 0.181, -0.132, -0.326, 0.118, -0.049, 0.056, 0.042]
112 # pacf值(保留小数点后三位)为: [-0.188, -0.245, 0.097, -0.134, -0.361, -0.126, -0.213, 0.036, -0.119]
113 # 上下边界值分别为: [-0.6324555320336759, 0.6324555320336759]

  5.14   最终结果验证:

acf值(保留小数点后三位)为: [-0.188, -0.201, 0.181, -0.132, -0.326, 0.118, -0.049, 0.056, 0.042]
pacf值(保留小数点后三位)为: [-0.188, -0.245, 0.097, -0.134, -0.361, -0.126, -0.213, 0.036, -0.119]
上下边界值分别为: [-0.6324555320336759, 0.6324555320336759]这里是计算的pacf的bound值。在看图的时候直接通用pacf的bound值。但是实际过程中这两个标准差是不太一样。具体见这篇博文http://www.cnblogs.com/noah0532/p/8453959.html

  对照上面给出的Eviews软件自动统计的值

  结果完全一致  

 6   TB(交易开拓者代码)

(持续编辑中。。。。。。。。。。)

原文地址:https://www.cnblogs.com/noah0532/p/8451375.html