R语言实战——基于KNN聚类的时间序列分析预测

---
title: "Summary of Reading"
author: "ChenWei"
date: "2020/6/24"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

    这一次我阅读的是Francisco Martinez等著的《Time Series Forecasting with KNN in R: the tsfknn Package》，这篇文章描述了如何使用$k$最近邻回归的tsfknn软件包进行时间序列预测。

一、背景

    传统上，使用统计方法(例如ARIMA模型或指数平滑)进行时间序列预测。但是，最近几十年出现了使用计算智能技术来预测时间序列。尽管人工神经网络是用于时间序列预测的最杰出的机器学习技术，但其他方法如高斯过程或KNN等也已应用。本文重点讲讲用KNN方法做时间序列预测。

    KNN算法又称$k$近邻分类(k-nearest neighbor classification)算法。它是根据不同特征值之间的距离来进行分类的一种简单的机器学习方法，它是一种简单但是懒惰的算法。他的训练数据都是有标签的数据，即训练的数据都有自己的类别。KNN算法主要应用领域是对未知事物进行分类，即判断未知事物属于哪一类，判断思想是，给出一个新的样本，KNN根据如欧氏距离等距离度量找到它最相似的$k$个样本，称为最近的邻居，以这$k$个样本中数目最多的类别作为新样本的归属类别。

二、原理分析

    先讲一讲用KNN做预测的原理。假定第$i$个训练样本的$n$个特征向量为 $(f_1^i,f_2^i,\cdots ,f_n^i)$，它们对应的$m$维特征属性目标为$(t_1^i,t_2^i,\cdots ,t_n^i)$ 。当给定一个特征向量为$(q_1,q_2,\cdots ,q_n)$ 的新样本时，我们用其特征寻找$k$ 个最相似的训练样本，近邻判定的标准为特征向量之间的欧式距离，新样本与第$i$个训练样本之间的欧氏距离的计算公式如下：

$$\sqrt{\sum _{x=1}^n{\left(f_x^i-q_x\right)}^2}$$

    假定找到的$k$个最近邻训练样本的目标向量分别为$(t^1,t^2,\cdots ,t^k)$ ，以它们的平均值作为需要预测的新样本的目标向量：

$$\sum _{i=1}^k\frac{t^i}{k}$$

    而KNN在时间序列中的应用方法为，以时间序列数据为目标值，目标值的滞后值为特征值，构造一个类似于自回归的模型。以12为滞后位数，得到的特征值与目标值如下：

Features	Target
x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12	x13
x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13	x14
x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14	x15
…	…
x120,x121,x122,x123,x124,x125,x126,x127,x128,x129,x130,x131	x132

三、单值预测实例分析

    接下来讲讲如何用在R中用KNN做时间序列预测。在tsfknn包当中只需要使用 knn_forecasting() 函数即可做时间序列预测，非常方便。该函数的调用方法为

$$pred=knn\_forecasting(time\_series,h=1,lags=1:12,k=2)$$

其中time_series是时间序列数据，h为预测的数据，lags代表用以预测的数据个数(即滞后数)，k为KNN中用到的最近邻个数；返回的pred包含了KNN模型的信息以及预测值，我们可以用它来做预测或者可视化模型数据。

    实例中使用的数据是包自带的nottem数据集，部分输入的特征数据及预测数据如下：

library(tsfknn)  # 载入包
timeS <- window(nottem,end=c(1930,12))  # 准备数据
pred <- knn_forecasting(timeS, h=1, lags=1:12, k=2) # 预测
head(knn_examples(pred))  # 展示前面的部分数据

    将这些数据可视化，得：

plot(pred,lwd=1.5,type='b',pch=20)  # 可视化数据

    如果需要查看预测数据，我们可以借助knnForecast对象pred实现：

pred$prediction

    想要将预测数据可视化也很方便，只需要调用ggplot2绘图包里面的 auplot(pred) 即可。此时我们可以在图中清楚地看到原始数据、特征数据、标签数据、实例数据以及预测数据，不同类型的数据有不同的颜色及点的类型。使用 knn_prediction() 函数做预测，得到预测的值为37.4，对应下图中的红点：

# 高级可视化
library(ggplot2)
autoplot(pred, highlight="neighbors", faceting=FALSE)

    此外，我们还可以用tsfknn这个包看到它是如何使用KNN做预测的，也就是说，我们可以使用 nearest_neighbors(pred) 函数清楚地看到被预测的数据对应的输入数据以及其对应的最近邻。其中的原理也很简单，我们使用1到12个特征数据做预测，而这些特征数据与所预测的数据之间是存在一定关联性的。

nearest_neighbors(pred)   # 最近邻

四、多值预测实例分析

    论文中除了预测未来的一个数据之外，还可以预测多个数据，这种情况也是现实情况中最常见的。文中讲到了两种方法——MIMO策略与递归策略。调用的方法也很简单，只需要在knn_forecasting函数中添加一个 ”msas”参数，然后将参数取值为 “MIMO” 或者 “recursive” 即可。

4.1 MIMO策略

     MIMO(Mutiple Input Mutiple Output) 策略在KNN中广泛使用。在该策略下，输入变量与输出变量的长度应当保持一致，并且它们在时间上应该是相邻的。

    在本节中我们使用过去12个月因车祸丧生的司机数预测未来12个月的因同原因丧生的司机数，数据集的名称为USAccDeaths。原始数据与预测数据的可视化结果如下图所示：

# MIMO策略
timeS <- window(UKDriverDeaths, end=c(1979,12))
pred <- knn_forecasting(timeS, h=12, lags=1:12, k=2, msas = "MIMO")  # 预测
autoplot(pred, highlight = "neighbors", faceting=FALSE)

4.2 递归策略

    递归策略常常用于ARIMA模型或者指数平滑预测模型。与MIMO策略不同的是，递归策略采用单步预测的方法。值得注意的是，当过去时间的特征数据不可用于新的实例的预测时，常常使用对应的过去的预测数据作为替代。

    在本节中，我们使用USAccDeaths中的数据作为例，预测未来两个月的数据。同样的，仍然使用过去的12个数据作为特征数据做预测，预测的结果如下图的第一幅图所示。

pred <- knn_forecasting(USAccDeaths,h =2,lags=1:12,k=2,msas = "recursive")
autoplot(pred,highlight = "neighbors")

上图中，第二幅图向我们展示了已经预测过的数据也会被当做特征值用于预测未来的数据。

五、预测精度分析

    既然KNN预测是一种预测方法，那么自然会相较真实值有所差距。而tsfknn包可以很方便地得到预测误差，只需要调用**ro=rolling_origin(pred,h)**即可。该函数将数据集分为训练集和测试集，并且在时间序列中测试集包含了最新的观测值。

    当不采用滚动规则的时候，使用最新的六个值预测六个值。使用 rolling_origin 函数可以得到测试集数据(test_sets)、预测数据(predictions)以及误差(errors)等。同时也可获得其他的预测精度判别指标，如RMSE、MAE及MAPE等。在不采用滚动规则时得到的预测结果如下表：

	h=1	h=2	h=3	h=4	h=5	h=6
Test_sets	1461	1354	1333	1492	1781	1915
Predictions	1513.5	1363.5	1351.5	1567	1587.5	2392
Errors	-52.5	-9.5	-18.5	-75	193.5	-477
RMSE=213.6137, MAE=137.6667, MAPE=7.7472

    当采用滚动规则的时候，得到的预测结果与不采用时一致，只是分为了$h=1,2,\dots 6$六种情况。使用测试集预测得到的两个观察值如下图中红点所示。

# 滚动计算误差
pred <- knn_forecasting(ldeaths, h=12, lags=1:12, k=2) # 拟合模型
ro <- rolling_origin(pred, h=6, rolling=FALSE)   # 不带滚动的
print(ro$test_sets)   # 测试集，h=6因此测试集有6个数
print(ro$predictions)    # 测试集数据对应的预测数据
print(ro$errors)     # 测试数据与预测数据之间的误差
print(ro$global_accu)   # 主要为RMSE MAE MAPE
plot(ro)

ro <- rolling_origin(pred, h=6)   # 滚动预测
print(ro$test_sets)   # 测试集
print(ro$predictions)  # 预测值
print(ro$errors)   # 预测误差
print(ro$global_accu)  
print(ro$h_accu)  # 不同h对应的accu，精确度
plot(ro,h=6)

六、KNN参数设置

    这一部分主要讲解一下包中与模型选择有关的几个特征参数，这几个特征参数主要为：

• 1、寻找最近邻的距离函数；

• 2、整合目标值的联合函数；

• 3、最近邻的个数，即参数$k$；

• 4、自回归的外生变量，即作为输入的特征的个数；

6.1 距离和联合函数

    包中默认使用的距离函数是欧氏距离，当然也可以使用其他的距离函数。

    包中的联合函数默认为平均值函数，但是也可以用其他的联合函数。不久前该包已经可以从平均值函数(mean)、中值函数(median)、加权函数(weighted combination)。如果要修改整合函数，可以在knn_forecasting函数中添加cb参数即可。

    接下来讲讲怎么使用加权函数作为联合函数，总体思想是距离新样本越近的点的权重越大。不妨设$d_i$和$t_i$依次为第$i$个最近邻与新样本之间的距离和第$i$个最近邻的标签值。接着，定义第$i$个最近邻的权重为$w_i=\frac{1}{d_i^2}$，则新样本的预测标签的形式如下所示：

$$target=\frac{\sum _{i=1}^k{w}_i{t}_i}{\sum _{i=1}^k{w}_i}$$

    特别的，当新样本与训练样本的距离为0时，使用上式计算会出现问题。在这种情况下，新样本被判为与其距离为0的训练样本对应的标签类别。

6.2 最近邻数量$k$

    对于特定的KNN模型，如果$k$很小，模型很容易受噪声的影响。相反的，如果$k$很大，距离时间$t$很遥远的样本仍然会影响到它的值。

    一般来说有几种比较好的选择参数$k$的方法。

• 第一种，也是最直白最快的方法，将$k$设置为训练样本个数的平方根；

• 第二种方法，选择能使测试集达到最优的$k$，即能够使得如MAPE等衡量预测准确度的统计量达到最小；

• 第三种方法，同时训练多个函数不同参数$k$的模型，然后取所有模型的预测值的平均值作为最终的预测值。以下为该种方法的一个示例。

pred <- knn_forecasting(ldeaths,h=12,lags=1:12,k=c(2,4)) # 多个k参数的预测模型
pred$prediction  # 预测值

七、与其他时间序列包的比较

    在R语言中有很多包可以做时间序列预测，比如GMDH、NTS等，在这其中forecast包算是做时间序列预测的最好的一个包。forecast包使用了stats包和其他许多包中所能找到的所有的方法，同时，里面还包含绘制季节性时间序列图像的函数以及做Box-Cox变换的函数。

    在论文中，作者从是否能够使用任意的自回归滞后时间、是否只需要时间序列数据和预测水平、是否有绘制模型相关信息图像的函数、是否应用了滚动预测的方法、是否生成了预测间隔、是否使用外生变量做预测等方面对GMDH、tsDyn如下表所示。

	GMDH	tsDyn	nnfor	forecast	tsfknn
任意滞后期数	no	no	yes	no	yes
默认参数	yes	no	yes	yes	yes
绘图函数	no	yes	yes	yes	yes
滚动估计	no	no	yes	yes	yes
区间预测	yes	no	no	yes	no
外生变量	no	no	yes	yes	no

    论文同时还从预测精度与运行时间两个方面比较了以上几个包中的预测方法。使用的数据集仍然为USAccDeaths数据集，当中包含111个月的工业数据，需要预测未来18个月的数据。此外，该数据集还有季节性和非季节性数据。

记NN3时间序列的预测数据为$F$，真实值为$X$。使用平均绝对误差百分比(MAPE)作为精度的衡量指标，其计算公式如下

$$\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}=\frac{100}{18}\sum _{t=1}^{18}\left|\frac{X_t-F_t}{X_t}\right|$$

    在具体比较中，因为GMDH最多只能够预测5个数据，因此不再考虑。对于nnfor包而言，分别使用其中的elm和mlp函数做极限学习机与多层感知器时间序列预测。对于forecast包而言，则分别使用当中专门求解ARIMA模型的auto.arima和求解指数平滑模型的ets函数。以上函数说需要的参数都选择默认设置的。比较结果如下表所示

	elm	mlp	auto.arima	ets	tsfknn
MAPE	0.07744391	0.01523495	0.06068908	0.06748664	0.08383561
Time	2.79s	4.75s	5.31s	0.97s	0.02s

很明显，不同的包中不同的时间序列预测方法的预测结果显著不同，从时间角度来看tsfknn的性能在以上几种方法中属于最好的，但是从精度上来看存在一定的缺陷。

相关代码与结果(注：运行时间可能与上表不完全一致，上表取的某一次运行得到的时间)：

library(forecast)
library(nnfor)
library(greybox)   # 计算MAPE的包
t.start <- proc.time()   # 开始时间
ts.elm <- elm(ldeaths)   # 极限学习机
print(cat('极限学习机的运行时间为:',proc.time()-t.start,'\t'))  # 运行时间
print(cat('极限学习机的MAPE值为：',MAPE(ldeaths[13:length(ldeaths)],ts.elm$fitted),'\t'))  # MAPE值

t.start <- proc.time()   
ts.mlp <- mlp(ldeaths)   # 多层感知机
print(cat('多层感知器的运行时间为:',proc.time()-t.start,'\t'))  # 运行时间
print(cat('多层感知器的MAPE值为：',MAPE(ldeaths[13:length(ldeaths)],ts.mlp$fitted),'\t'))  # MAPE值

t.start <- proc.time()   
ts.arima <- auto.arima(ldeaths)   # arima时间序列
print(cat('arima模型的运行时间为:',proc.time()-t.start,'\t'))  # 运行时间
print(cat('arima模型的MAPE值为：',MAPE(ldeaths,ts.arima$fitted),'\t'))  # MAPE值

t.start <- proc.time()   
ts.ets <- ets(ldeaths)   # arima时间序列
print(cat('指数平滑模型的运行时间为:',proc.time()-t.start,'\t'))  # 运行时间
print(cat('指数平滑模型的MAPE值为：',MAPE(ldeaths,ts.ets$fitted),'\t'))  # MAPE值

library(tsfknn)  
t.start <- proc.time()  
ts.tsfknn <- knn_forecasting(ldeaths[1:60], h=12, lags=1:12)   # knn聚类
print(cat('k最近邻回归模型的运行时间为:',proc.time()-t.start,'\t'))  # 运行时间
print(cat('k最近邻回归模型的MAPE值为：',MAPE(ldeaths[61:72],ts.tsfknn$prediction),'\t'))  # MAPE值