最近接触时间序列较多，在借鉴很多人的知识之后，特此总结一下。目前关于时间序列数据分析预测大致有三种主流方法：
1、ARIMA系列方法
2、facebook开源的Prophet模型
3、LSTM时间序列预测

本系列希望在项目和实践的角度，用python实现上述三种方法并做出对比总结。如有不足之处，感谢指出，虚心改正。
所需环境： win10/ubuntu均可，python3.6.x，pandas，numpy，tensorflow2.0，statsmodels，pmdarima…

（1）项目简介
本文项目中所用数据为近一段时间内，间隔30分钟采样的气象数据（包括温度、湿度、风速、风向等数据）。在本文的理解中，arima方法仅支持单变量预测，也就是需要单独取出某列进行该列的预测。若想要多变量输入多变量输出，arima方法要拟合多个模型，再整合输入输出。
（2）数据介绍
数据是我从某气象网站爬取到的气象数据，已经做好数据清洗，提供一组测试数据放在文末，需要自取哈，举例如下：

表中单位：温度(华氏度)，湿度%，风向可转换角度，风速（mph）
读取csv数据，数据清洗：

    def read_temperature():
        # 以首列的日期时间作为时间戳 此处要转换成pandas的DatetimeIndex格式
        data = pd.read_csv('data/test_2020.csv', usecols=[0, 1, 2, 3, 4])
        data['Date Time'] = pd.to_datetime(data['Date Time'])
        data.set_index('Date Time', inplace=True)
        data = data['Temperature']  # DatetimeIndex对象
        # 华氏度转摄氏度
        for i in range(len(data)):
            data[i] = round(5 / 9 * (data[i] - 32))
        # 时间频率重采样
        data = data.resample('30T').mean()  # 更新频率30分钟，缺省值用临近的后一值补充
        data = data.fillna(data.bfill())
		# 显示数据
        data.plot(figsize=(15, 12))
        plt.show()
        return data

当前数据图：

(3) 平稳性检验
arima时间序列数据预测首先要保证数据的平稳性，肉眼可以看到当前数据可能有一定上升趋势。但还是要有数据量化这个平稳性，这里我们使用ADF平稳性检验，基于statsmodels库实现：

    def adf_stability_test(self, data):
        # ADF单位根检验：不存在单位根，表示数据平稳。且不能是白噪声数据(随机数)
        x = data
        res = ts.adfuller(x)
        lb_res = lb_test(x, None, True)[1]

        tag = False
        for i in range(len(lb_res)):
            if lb_res[i] < 0.05:
                continue
            else:
                print('序列为白噪声!')
                tag = True
                break
        if res[0] < res[4]['1%'] and res[0] < res[4]['5%'] and res[0] < res[4]['10%'] and res[1] < 0.05 and not tag:
            print('平稳序列！非白噪声')
            # plt.plot(lb_test(x, None, True)[1])
            # plt.ylabel('p-Value')
            # plt.show()
            return True
        else:
            return False

adfuller方法返回值：

(-6.260374222209651, 4.237742412839035e-08, 5, 906, {'1%': -3.4375883271133243, '5%': -2.8647353885968214, '10%': -2.568471435365895}, 1615.9162778870964)

ADF单位根检验，如果序列平稳，就不存在单位根；否则，就会存在单位根。该检验先假设序列不平稳，存在单位根。如果得到的显著性检验统计量小于三个置信度（10%，5%，1%），则对应有（90%，95，99%）的把握来拒绝原假设 ->序列平稳。
adfuller的返回值意义：

adf：Test statistic，T检验，假设检验值。
p-value：假设检验结果。
usedlag：使用的滞后阶数。
nobs：用于ADF回归和计算临界值用到的观测值数目。
icbest：如果autolag不是None的话，返回最大的信息准则值。
resstore：将结果合并为一个dummy。

具体原理我也不是很懂，只需知道adf值，也就是本文的-6.260374222209651 均小于（10%，5%，1%）三个置信度的值，且p-value < 0.05时能够较好的拒绝原假设，表示序列平稳。
若序列不平稳，需要对序列连续差分，直到序列平稳，对差分后的序列进行预测，再将预测结果逆差分，得到我们需要的值。

    def stability_test(self):
        count = 0
        flag = False
        if not self.adf_stability_test(self.data):
            while not self.adf_stability_test(self.data):
                count += 1
                self.data = self.data.diff(count)			# 差分
                self.data = self.data.fillna(self.data.bfill())	# 缺省值用后一值补充
                flag = True
        print('经过{}次差分，序列平稳'.format(count)) 
        return count, flag

count是让序列变平稳的差分次数，也就是ARIMA(p,d,q)中的d，记录count值便于后面调pdq参数。flag记录是否经过差分，便于后期作逆差分还原数据。
差分操作很容易理解：原数列的后一项减去前一项得到的一组差数列。
相当于序列整体后移，空出第一位是空NAN，再与原序列做差

0 1 2 3 4 5 6
3 6 7 8 9 10 12 
# 一阶差分
0 	 1 2 3 4 5 6
NAN  3 1 1 1 1 2  

我们知道arima或者是seasonal arima最重要的是（p,d,q）和（P,D,Q,S）这四个参数的确定。
而相对准确的方法是计算acf,pacf图，通过观察自相关图和偏相关图，确定拖尾还是截尾等什么鬼的…我到现在还似懂非懂。总之了解下来是要肉眼观察参数，不能自动化调参…
但我们的数据是实时变化的，参数可能变动很大。自动化训练参数很有必要：
(4) 模型参数训练
目前了解两种自动确定参数方法：
1）网格搜索的方法：
感觉就是暴力搜索，用一些模型评价标准确定最优的参数，计算量很大复杂度高。不在乎效率也值得一试，此处用AIC准则评价模型，找出AIC评分最小的模型参数

    def temp_param_optimization(self, data):
        paramBest = []
        warnings.filterwarnings("ignore")
        p = q = range(0, 3)		# 限制pq范围
        pdq = [(x[0], self.d, x[1]) for x in list(itertools.product(p, q))]	# 此处的d我们已经得到
        seasonal_pdq = [(x[0], self.d, x[1], s) for x in list(itertools.product(p, q))]		# 此处的s是季节性参数，可根据数据指定，变化不大
        
        for param in pdq:
            for param_seasonal in seasonal_pdq:
                try:
                    mod = sm.tsa.statespace.SARIMAX(data,
                                                    order=param,
                                                    seasonal_order=param_seasonal,
                                                    enforce_stationarity=False,
                                                    enforce_invertibility=False)

                    results = mod.fit()
                    paramBest.append([param, param_seasonal, results.aic])
                    print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
                except:
                    continue
        # print('paramBest:', paramBest)
        minAic = sys.maxsize
        for i in np.arange(len(paramBest)):
            if paramBest[i][2] < minAic:
                minAic = paramBest[i][2]
        # print("minAic:", minAic)
        for j in np.arange(len(paramBest)):
            if paramBest[j][2] == minAic:
                return paramBest[j][0], paramBest[j][1]

2）pmdarima
了解过时间序列发现，很多用R语言作时间序列预测的。因为R提供自动训练参数的库auto_arima。也是在项目后期才发现python也有类似的开源库pmdarima，貌似就是模仿R而来的。
直接pip install pmdarima -i https://pypi.tuna.tsinghua.edu.cn/simple
，pip安装要记得换源

    def auto_parameters(self, data, s_num):
        kpss_diff = arima.ndiffs(data, alpha=0.05, test='kpss', max_d=s_num)
        adf_diff = arima.ndiffs(data, alpha=0.05, test='adf', max_d=s_num)
        d = max(kpss_diff, adf_diff)
        D = arima.nsdiffs(data, s_num)
        
        stepwise_model = auto_arima(data, start_p=1, start_q=1,
                                    max_p=9, max_q=9, max_d=3, m=s_num,
                                    seasonal=True, d=d, D=D, trace=True,
                                    error_action='ignore',
                                    suppress_warnings=True,
                                    stepwise=True)
        print("AIC: ", stepwise_model.aic())
        print(stepwise_model.order)		# (p,d,q)
        print(stepwise_model.seasonal_order)	# (P,D,Q,S)
        print(stepwise_model.summary())		# 详细模型
        return stepwise_model.order, stepwise_model.seasonal_order

十分方便的库，arima.ndiffs()，arima.nsdiffs()可以方便的求出最合适的d,和D参数。代码中s_num是自己指定的季节性参数。

auto_arima(data, start_p=1, start_q=1,
		     max_p=9, max_q=9, max_d=3, m=s_num,
		     seasonal=True, d=d, D=D, trace=True,
		     error_action='ignore',
		     suppress_warnings=True,
		     stepwise=True)

start_p，start_q设置初始p,q参数，max_p，max_q最大值。stepwise=True貌似会直接选择最优参数。参数还有很多，请查阅官方api：https://alkaline-ml.com/pmdarima/index.html
对比网格搜索的方法，本文觉得两种方法确定的参数几乎相同，但auto_arima训练参数特别快。同样的数据auto_arima要快一半以上。
（5）模型预测

    def model_prediction(self, data, param, s_param, n_steps, flag):
        mod = sm.tsa.statespace.SARIMAX(data,
                                        order=param,
                                        seasonal_order=s_param,
                                        enforce_stationarity=False,
                                        enforce_invertibility=False)
        results = mod.fit()

        pred_uc = results.get_forecast(steps=n_steps)		# n_steps可指定预测的步数(多少时间间隔)
        pred_ci = pred_uc.conf_int()

        if flag:  # 还原差分
            pred_res = pd.Series([data[0]], index=[data.index[0]]).append(pred_uc.predicted_mean).cumsum()
            print("预测结果(℃):  ", pred_res)
            return pred_res
        else:
            print("预测结果(℃):  ", pred_uc.predicted_mean)
            return pred_uc.predicted_mean

传入参数param，s_param分别对应（pdq）(PDQS)
测试数据预测结果：

---

4.21更新
感觉对于数据的季节性参数设置可能有点迷。比如文中提到的网格搜索的参数s和auto_arima的参数m都是要自己设置的，当然要根据所用数据的频率间隔来设置季节性。我的理解是数据规律的周期性呈现方式。比如，我用的是每30分钟的观测数据，一个月的数据。就可以设置为日周期性30*48=1440，一天共48个间隔，m或者s参数设置为48可能比较合适，但还是要针对数据具体分析。
如果各位同学可以科学上网的话，https://robjhyndman.com/hyndsight/seasonal-periods/ 这位博主关于季节性参数写的比较好！
针对一般情况的设置：

链接：https://pan.baidu.com/s/1QiXhYYCOti2PiNcTTrraBg
提取码：w2di
但有些数据单位需要转换哈，像温度是华氏度