时序预测之三 _ 傅立叶和小波变换

1. 说明

  用傅立叶变换预测时序数据，原理是把时域数据转换到频域，再转换回来．python 的 numpy 和 scipy 里面都有现成的转换工具 fft() 和 ifft()，但使用时会遇到一个问题：比如 25 天的数据转到频域再转回时域，还是 25 天，虽然拟合了数据，但没法直接预测未来，本篇介绍用它实现预测的方法．

2. 傅立叶变换

(1) 相关知识

  之前写过关于傅立叶变换原理的文档，这次就不再重复了，具体请见：https://www.jianshu.com/p/9e786be6dccb

  本篇只从程序的角度看如何使用它．经过 FFT 转换的数据和转换前长度一致，每个数据分为实部和虚部两部分，假设时序时数长度为 N（Ｎ最好是 2 的整数次幂，这样算起来更快），用 fft() 转换后：下标为 0 和 N /2 的两个复数的虚数部分为 0，下标为 i 和 N - i 的两个复数共辄，也就是其虚部数值相同、符号相反。再用 ifft() 从频域转回时域之后，出现了由误差引起的很小的虚部，用 np.real() 取其实部即可．   由于一半是另一半的共轭，因此只需要关心一半数据．fft 转换后下标为 0 的实数表示时域信号中的直流成分（不随时间变化），下标为 i 的复数 a + bj，其中 a 表示余弦成分，b 表示其正弦成分．

(2) 示例功能

  数据是航空乘客数据＂AirPassengers.csv＂，可以从 CSDN 下载，其中包括从 1949-1960 年，每月旅客的数量，程序预测未来几年的旅客数据．

  如图所示，数据为非平稳数据，其趋势向上，且波动加俱，为将其变为平稳数据，先对其做了对数和差分处理．

(3) 示例代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 函数功能：将频域数据转换成时序数据
# bins为频域数据，n设置使用前多少个频域数据，loop设置生成数据的长度

def fft_combine(bins, n, loops=1):
 length = int(len(bins) * loops)
 data = np.zeros(length)
 index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi)
 for k, p in enumerate(bins[:n]):
   if k != 0 : p *= 2 # 除去直流成分之外, 其余的系数都 * 2
   data += np.real(p) * np.cos(k*index) # 余弦成分的系数为实数部分
   data -= np.imag(p) * np.sin(k*index) # 正弦成分的系数为负的虚数部分
 return index, data

if __name__ == '__main__':
 data = pd.read_csv('AirPassengers.csv')
 ts = data['Passengers']

# 平稳化
 ts_log = np.log(ts)
 ts_diff = ts_log.diff(1)
 ts_diff = ts_diff.dropna()
 print(fy[:10]) # 显示前10个频域数据
 fy = np.fft.fft(ts_diff)
 conv1 = np.real(np.fft.ifft(fy)) # 逆变换
 index, conv2 = fft_combine(fy / len(ts_diff), int(len(fy)/2-1), 1.3) # 只关心一半数据
 plt.plot(ts_diff)
 plt.plot(conv1 - 0.5) # 为看清楚，将显示区域下拉0.5
 plt.plot(conv2 - 1)
plt.show()

(4) 运行结果

1
2
3
4
[ 1.34992672+0. j -0.09526905-0.14569535j -0.03664114-0.12007802j
 -0.2670005 +0.24512406j -0.10075074+0.0314084 j -0.26409417+0.04197159j
 0.14411338+0.18703009j 0.07467991+0.05367644j -0.26663142+0.15324939j
 0.03248223+0.14130114j]

(5). 示例分析

  输出的是 fft 转换后的数据，只显示了前十个，形式为复数．复数模 (绝对值) 的两倍为对应频率的余弦波的振幅；复数的辐角表示对应频率的余弦波的相位。第 0 个元素表示直流分量，虚部为 0．在数据中的位置标记了频率大小，值标记了振幅大小．   图中显示的三条曲线分别为原始数据，做了 fft 以及 ifft 逆变换后的数据，以及 fft 后自己实现算法还原并预测了未来的数据，从图中可见，基本拟合了原始曲线，预测曲线看起来也比较合理．   上述方法可实现用傅里叶变换预测时序数据．与 ARMA 算法相比，它没有明显衰减，更适合长时间的预测．   对于随时间变化的波形，比如语音数据，一般使用加窗后做傅立叶变量的方法拟合数据．

3. 小波变换

(1) 相关知识

  有了傅立叶变换，为什么还用小波呢？上面提到，如果波型随时间变化，就需要对波型加窗分段后再处理，而且有时候需要大窗口，有时候需要小窗口，处理起来就更加麻烦．于是引入了更灵活的小波．   傅立叶变换的基是正余弦函数，而小波的基是各种形状小波，也就是说它把整个波形看成是多个位置和宽度不同的小波的叠加．小波有两个变量：尺度 a 和平移量 t，尺度控制小伸的伸缩，平移量控制小波的平移，它不需将数据切分成段，就可以处理时变数据．尤其对突变信号处理得更好．

  下图是几种常见的小波．

  离散小波变换，Discrete Wavelet Transformatio (dwt)，可以说是小波变换中最简单的一种。这里使用 Python 调用 pywt 库实现最简单的功能．   经过变换之后的返回值：cA:Approximation（近似）, cD:Detail（细节），其中近似 cA 是周期性有规律的部分，可以被模拟和预测，而 cD 可看做是噪声。换言之，用此方法可以拆分周期性数据，和其上的扰动数据。

(2) 示例功能

  示例使用的仍然是乘客数据，下面代码是将细节 D 设为 0，然后还原。

(3) 示例代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import pywt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':
 data = pd.read_csv('AirPassengers.csv')
 ts = data['Passengers']
 # 平稳化
 ts_log = np.log(ts)
 ts_diff = ts_log.diff(1)
 ts_diff = ts_diff.dropna()

 cA,cD = pywt.dwt(ts_diff, 'db2')
 cD = np.zeros(len(cD))
 new_data = pywt.idwt(cA, cD, 'db2')

 plt.plot(ts_diff)
 plt.plot(new_data - 0.5)

 plt.show()

(4) 运行结果

(5) 示例分析

  可以看到，用小波拟合的效果也还可以，一般可以使用小波拟合 cA，使用 ARMA 拟合 cD 部分，两种方法配合使用．