基于梯度下降算法自建一种短期有效的自回归模型

前言
一：移动平均模型
二：基于自适应滤波思想的权重优化
三：代码实现
四：实验分析
五：总结与展望

前言

基于时间序列自回归预测模型还是比较多的，简单的有移动平均，灰色预测，AR等等，复杂的有ARIMA，GARCH、LSTM，TCN等等。自回归模型说白了就是“当下的自己”跟“过去的自己”建立回归模型来预测“未来的自己”，它不需要任何其它的自变量，是个易理解与易应用的模型。如果自回归模型想要好的预测效果，那么我们还是希望数据随时间变化是稳定的或缓慢变化的，或者呈周期性季节性变化的短期预测。
移动平均模型就是典型的自回归模型，本文主要讲解移动平均模型的实现与运用，并不断深入最后到基于自适应思想的权重优化移动平均模型。尽管时序有很多现成复杂模型，但实际运用中，移动平均深入优化后也是大有用武之地，尤其对于短期和有规律数据。

虽曰如云，匪我思存，写作不易，走过路过的朋友们，别忘了点赞收藏加关注一下哈，谢谢！

一：移动平均模型

简单移动平均

简单移动平均是用某一时刻点前 n 期的数值之简单平均数来预测该时刻点的数值，简单用数学公式表达为

xt+1=xt+xt−1+...+xt−n+1bx_{t+1}=\frac{x_t+x_{t-1}+...+x_{t-n+1}}{b}xt+1=bxt+xt−1+...+xt−n+1，

其中 xt+1x_{t+1}xt+1 是下一时刻的预测值， xt−ix_{t-i}xt−i 时候前 iii 时刻的具体数值，以后每再预测下一时刻，则就是要从移动平均中去除最早的一个观测值，再加上一个最新的观测值。

加权移动平均

简单移动平均实际上是对每一时刻的数值赋予相同权重，加权平均则就是在权重上变化，其余是跟简单移动平均是一样的，数学公式表达

xt+1=wtxt+wt−1xt−1+...+wt−n+1xt−n+1x_{t+1}=w_tx_t+w_{t-1}x_{t-1}+...+w_{t-n+1}x_{t-n+1}xt+1=wtxt+wt−1xt−1+...+wt−n+1xt−n+1，

其中 wiw_iwi 为第 iii 时刻数值前的权重值，其余参数跟上述相同含义，在实际中，我们一般会把越接近预测点的数据点前的权值设定大一些，反之相反。

指数移动平均

原本指数加权移动平均是不想舍弃过去数据，只是随着时间点的远离，其权值追减递减，而是随着时间间隔的增大成指数衰减，但实际上越往后权重衰减的非常严重，完全可以忽略，所以我们还是正常选定前 nnn 期数据作为建模数据，表达公式如下

xt+1=αxt+α(1−α)xt−1+...+α(1−α)(t−n+1)xt−n+1x_{t+1}=\alpha x_t+\alpha(1-\alpha)x_{t-1}+...+\alpha(1-\alpha)^{(t-n+1)}x_{t-n+1}xt+1=αxt+α(1−α)xt−1+...+α(1−α)(t−n+1)xt−n+1

其中 α\alphaα 为平滑系数为一常数，一般取值介于0.1~0.4之间。

二：基于自适应滤波思想的权重优化

内容简介

自适应滤波器是指根据环境的改变，使用自适应来改变滤波器的参数和结构的滤波器（这里主要涉及控制理论，我们不过多展开，有兴趣可查阅相关资料）。我们主要基于其思想构造权重优化算法如下：

1.设输入参数向量为 Xn=(x1,x2,...xn)X_n=(x_1,x_2,...x_n)Xn=(x1,x2,...xn) ，已知权重向量 Wp=(w1,w2,...wp)W_p=(w_1,w_2,...w_p)Wp=(w1,w2,...wp) ，其中1≤p≤n1\leq p\leq n1≤p≤n ,

2.我们构造误差 Lossp+1=xp+1ˉ−xp+1Loss_{p+1}=\bar{x_{p+1}}-x_{p+1}Lossp+1=xp+1ˉ−xp+1，

其中xp+1ˉ=Xp⋅WpT\bar{x_{p+1}}=X_p\cdot W_p^Txp+1ˉ=Xp⋅WpT ，那么均方误差定义为

L2=(xp+1ˉ−xp+1)2=(xp+1ˉ)2−2xp+1ˉXpWpT+(XpWpT)2=(xp+1ˉ)2−2xp+1ˉXpWpT+Xp(WpTWp)XpTL^2= (\bar{x_{p+1}}-x_{p+1})^2 =(\bar{x_{p+1}})^2-2\bar{x_{p+1}}X_pW_p^T+(X_pW_p^T)^2 =(\bar{x_{p+1}})^2-2\bar{x_{p+1}}X_pW_p^T+X_p(W_p^TW_p)X_p^TL2=(xp+1ˉ−xp+1)2=(xp+1ˉ)2−2xp+1ˉXpWpT+(XpWpT)2=(xp+1ˉ)2−2xp+1ˉXpWpT+Xp(WpTWp)XpT ,

3 为了找出 WpW_pWp 使 L2L^2L2 函数最小，我们一是用 ∂L2∂Wp=0\frac{\partial L^2}{\partial W_p}=0∂Wp∂L2=0 求出 WpW_pWp，但这样我们会涉及大量的方程组运算，以及数据集用不完整的情况，于是我们采用梯度下降算法，

4.根据2中 L2L^2L2 的定义，我们有

∂L2∂Wp=−2xp+1ˉXp+2XpWpTXp=−2(xp+1ˉ−XpWpT)Xp=−2LXp\frac{\partial L ^2}{\partial W_p}= -2\bar{x_{p+1}}X_p+2X_pW_p^TX_p =-2(\bar{x_{p+1}}-X_pW_p^T)X_p =-2LX_p∂Wp∂L2=−2xp+1ˉXp+2XpWpTXp=−2(xp+1ˉ−XpWpT)Xp=−2LXp ,

那么更新 WpW_pWp 的公式为Wp=Wp+2lLXpW_p=W_p+2lLX_pWp=Wp+2lLXp (开始的 WpW_pWp 值来源于初始化的一组已知数据)，至于 WpW_pWp 函数能否收敛，常数 lll的取值至关重要，要满足维纳解的条件，有兴趣的同学可以参考这方面的资料，在算法实现上 lll 取值不能太大，一般太小收敛速率会比较慢，

5.如此上述步骤会更新一组 Wp1W_p^1Wp1 值，再输入一个新的已知数据 xp+2x_{p+2}xp+2 ，把新 Wp1W_p^1Wp1值带入上述2-4中，于是会得到新的误差函数 Lossp+2Loss_{p+2}Lossp+2和新的 Wp2W_p^2Wp2。如此反复遍历完所有数据集 XnX_nXn就会得到一组更新 n−pn-pn−p 次的Wpn−pW_p^{n-p}Wpn−p ，最后一轮下来定义误差函

C=1n−p∑pn−1(xp+1−Xpˉ(Wpn−p)T)2C=\frac{1}{n-p}\sum_{p}^{n-1}{\left( x_{p+1} -\bar{X_p}(W_p^{n-p}\right)^T)^2}C=n−p1∑pn−1(xp+1−Xpˉ(Wpn−p)T)2 ，

这里 Xpˉ\bar{X_p}Xpˉ 为长度为 ppp ，且一次只向前推移1步的样本数据向量，很明显 Xpˉ\bar{X_p}Xpˉ 最后只能推移到数据集的第 n−1n-1n−1 的位置，

6.对5中的 CCC 设定条件，如果 CCC 满足 C<α(αC<\alpha( \alphaC<α(α为人为阈值），则最终退出程序，或者设定一个最大循环次数，满足循环结束也可退出程序，否则继续把循环过一数据集得到的Wpn−pW_p^{n-p}Wpn−p带再入数据集 XnX_nXn ，重新开始计算 LossLossLoss 和 CCC ，继续重复2-5步骤的计算即可，

7.如果对深度学习熟悉的朋友，其实这里的每次给一个新的已知数据相当于深度学习中的batchsize=1，iteration= n−pn-pn−p即计算剩下 n−pn-pn−p 个新数据还需要iteration次，那么最后的epoch相当于我们彻底完整更新了多少次WpW_pWp。

三：代码实现

先初始化一些参数如下

###初始化数据
n = 5##阶数这里设为5
d = 0.1
init = 0
E = np.inf##开始误差默认无穷大
E_ls = []##误差列表
count_ls = []##计数列表
undate_w_ls = []##权重更新列表
winit_ls = np.random.uniform(0,1,n)##初始权重为标准均匀分布随机数
steps = n##滑动步数

再建立数据集用于算法建模使用。我们建立三种数据类型，一是有周期性的三角函数，二是有平稳增长的对数函数，三是具有趋势的实际销量数据

def bulid_experiment_data(choose,yv=0,xv=0):# yv = 0if choose == 'sin':yv = np.sin(xv) + 1.5elif choose == 'log':yv = 2 + np.log(xv) + 1else:xv = [i for i in range(len(yv))]return xv,yv# choose = 'sin'
# l = 1.0e-3###学习率
# min_error = 1.0e-8#设置误差标准
# x_sin = np.arange(-2 * np.pi, 2 * np.pi, d)
# sin_or_log = bulid_experiment_data(choose,xv=x_sin)# choose = 'log'
# l = 1.0e-3
# min_error = 5.0e-6#设置误差标准
# x_log = np.arange(1, 10, d)
# sin_or_log = bulid_experiment_data(choose,xv=x_log)choose = 'else'
l = 3.0e-9
min_error = 20000#设置误差标准
yv_all = np.array([3023,3039,3056,3138,3188,3224,3226,3029,2859,2870,2910,3012,3142,3252,3342,3365,3339,3345,3421,3443,3428,3554,3615,3646,3614,3574,3635,3738,3707,3827,4039,4210,4493,4560,4637,4755,4817])
yv = yv_all[:-n]
sin_or_log = bulid_experiment_data(choose,yv=yv)

接着开始我们的自适应算法主逻辑编写

start = time.time()
def auto_filter_measure(v_ls, w_ls, steps, init, E_ls,l,n):if steps == len(v_ls):for i in range(len(v_ls) - n):res = np.dot(v_ls[i:i + n], w_ls)E_ls.append(res)S = mean_squared_error(E_ls,v_ls[n:])##误差平方和return S, w_lselse:v_temp = np.dot(v_ls[init:steps], w_ls)##内积运算# 计算预测误差error_temp = v_ls[steps] - v_tempw_newls = w_ls + 2 * l * error_temp * v_ls[init:steps]###更新权重steps = steps + 1init = init + 1return auto_filter_measure(v_ls, w_newls, steps, init, E_ls,l,n)##递归调用函数

然后不停迭代上述算法直到损失函数达到阈值，并记录损失函数图和相关结果

# ###执行主逻辑
R_afm = auto_filter_measure(sin_or_log[1], winit_ls, steps, init, E_ls,l,n)
error_ls = []while abs(E) > min_error:E_ls = []#每次都要初始化误差列表if type(R_afm) == tuple:##第一次res的输出r = auto_filter_measure(sin_or_log[1], R_afm[1], steps, init, E_ls,l,n)else:r = auto_filter_measure(sin_or_log[1], R_afm, steps, init, E_ls,l,n)E = r[0]count_ls.append(E)R_afm = r[1]####R_afm在这更新undate_w_ls.append(R_afm)error_ls.append(E)print('均方误差值：', E)stop = time.time()
plt.figure(figsize=(15,8))
plt.plot(error_ls,label='%s函数下，学习率为%s且要求最小绝对精确误差为%s'%(choose,l,min_error))
plt.ylabel('MSE',fontsize=15)
plt.xlabel('迭代次数',fontsize=15)
plt.tick_params(labelsize=15)
plt.legend(fontsize=18)
plt.show()print('\033[1;31m选择的是%s函数，默认设定自回归阶数为：%s阶\033[0m'%(choose,n))
print('\033[1;32m原始数据长度：%s，最小条件误差为%s\033[0m'%(len(sin_or_log[1]),min_error))
print('\033[1;33m最终迭代次数：%s，才达到最佳权值\033[0m'%(len(count_ls)))
print('\033[1;34m最佳权值组合：%s\033[0m'%(np.round(undate_w_ls[-1],3)))
print('\033[1;35m计算用时：%.3f秒！\033[0m'%(stop-start))

最后我们再把另外三种移动平均算法加入，进行对比预测分析并做出可视化图形

#################样本内预测
def predict_inner(data,step,w):pyinner = []for i in range(len(data)-step):res = np.dot(data[i:i+step],w)pyinner.append(res)return pyinnerpreinner_afm = predict_inner(sin_or_log[1],n,undate_w_ls[-1])##简单移动平均的权重值
w_sma= [1/n for i in range(n)]
w_sma = np.array(w_sma)##加权移动平均的权重值
w_wma = np.arange(0.008,n * 0.008,0.0085)
w_wma = list(w_wma)
w_wma.append(1 - sum(w_wma))
w_wma = np.array(w_wma)###指数移动平均的权重值
a = np.arange(0.5,1,0.05)
rd = np.random.choice(a)
w_ema = [rd]for i in range(1,n):w_ema.append(rd * np.power((1 - rd),i))preinner_sma = predict_inner(sin_or_log[1],n,w_sma)
preinner_wma = predict_inner(sin_or_log[1],n,w_wma)
preinner_ema = predict_inner(sin_or_log[1],n,w_ema)####################样本外预测
w_afmfin = undate_w_ls[-1]##取最后更新好的权重向量
data_afmfin = list(sin_or_log[1][-n:])##自适应方法下取原始样本的最后n个数据
data_smafin = list(sin_or_log[1][-n:])##简单平均方法
data_wmafin = list(sin_or_log[1][-n:])##加权平均方法
data_emafin = list(sin_or_log[1][-n:])##指数平均方法def predict_outer(data,w):pyouter = np.dot(data,w）return pyouterdef update_predict_outer(data, pred, pls, w, count=0):if count == pred:print('\033[1;38m未来%s期预测结果：\033[0m' % (len(pls)))print('\033[1;38m%s\033[0m'% (np.round(pls,5)))return plselse:rp = predict_outer(data,w)pls.append(rp)##预测列表追加一个预测值data.append(rp)##原始样本追加一个预测值data_update = data[-n:]##再取倒数len(w_fin)个数据作为自回归模型的自变量count = count + 1##计数加1直到等于我们设定的预测长度为止，if逻辑改遍历循环也行return update_predict_outer(data_update, pred, pls, w, count)try:prefd = int(input('\033[1;38m请输入未来预测期数（正整数）：\033[0m'))由于篇幅原因，作图代码这里就不再上传，如有需要请别客气联系作者~~~~~~

四：实验分析

我们知道，回归预测类算法对于短期预测是非常有效的，长期的话有太多的扰动影响。此次实验数据会带有周期和趋势的特点，如果数据集有异常，是可以先进行异常点分析，如果数据没有任何规律而言，那么自回归模型并不是非常合适的选择。

sinsinsin 周期函数实验

图1

说明：图1展示的是学习率和误差精度阈值以及阶数ppp的确定下， sinsinsin 函数在训练时期的 MSEMSEMSE 随迭代次数的变化值，我们知道学习率，阈值的取值不同，迭代过程也会有变化的，但趋势是一样的，我们就不再过多实验。基于经验而言，阶数 ppp 的取值相对于数据样本，最好不要过大。

图2

说明：图2就是我们四种算法的具体实验结果，很明显自适应时序算法对周期函数的预测明显有效，其他三种算法（尤其是加权移动和指数移动，并多次调节权重下）在样本内预测还有些周期性，但样本外都根本不适合做中长期预测。

图3

说明：图3就是我们这次实验的样本内外预测的 MSE 分布情况。一样很明显自适应算法的是比其他三种算法有效的。

logloglog 对数趋势函数实验

图4

图5

图6

说明：logloglog函数数据的测试方法跟 sinsinsin 很类似，不过最终结果也是符合我们预期的，自适应算法在中短期预测方面无论是样本内还是样本外，明显是优于其他三种。但是我们还是强调一点，自回归预测方法对长期预测并不友好，一是从算法角度考虑，由于自身过去的数据结构原因，过长预测会出现数据值平稳不振荡，或者就是另一个极端数据太过振荡，二是从实际考虑，长期影响因素太多。

基于实际销量数据的短期预测

图7

图8

图9

说明：由于我们的原始数据比较大，对应的 MSE 的值也会很大，但不影响算法之间的评估。在图7中，我们设定的 MSE 不能过小，一是会大量耗时，越往后训练， MSE 下降速率越慢，二是防止样本内训练过于太好而不能泛化模型，所以要多次调参来选择。从图8和图9明显观察到自适应方法在短期预测中的又一次“强大”，样本外的预测非常接近原始值了。

五：总结与展望

从本篇内容看，自己搭建的自适应自回归模型对有趋势和周期的数据，在中短期预测上表现不俗，在理解上较于已存在的高级模型也容易的多，所以说善于利用模型也是很重要的部分。
自适应自回归模型在一些方面也是可以改进的。如果不考虑长期预测外界因素的扰动，阶数的取定是个值得思考的问题，当然这里可以参考ARIMA模型的Pacf图或者信息准则等方法，有兴趣的朋友可以试一试，相互交流。