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

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) 示例代码

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 # 除去直流成分之外, 其余的系数都 * 2data += np.real(p) * np.cos(k*index) # 余弦成分的系数为实数部分data -= np.imag(p) * np.sin(k*index) # 正弦成分的系数为负的虚数部分return index, dataif __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.5plt.plot(conv2 - 1)
plt.show()

(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.04197159j0.14411338+0.18703009j 0.07467991+0.05367644j -0.26663142+0.15324939j0.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) 示例代码

import pywt
import pandas as pd
import numpy as np
import matplotlib.pyplot as pltif __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部分，两种方法配合使用．