自回归模型
p阶自回归模型AR(p)是一种在描述时间序列方面特别有效的随机时间序列模型。在这个模型中,时间序列的现在值$x_t$是用该序列过去数值的线性组合加上一个白噪声扰动项$\epsilon_t$来表示。
$$x_t=\phi_0+\phi_1x_{t-1}+\phi_2x_{t-2}+\cdots+\phi_px_1+\epsilon_t\tag{1}$$
其中包含三个限制条件:模型的最高阶数为p,即$\phi_p \ne 0$;随机干扰序列$\epsilon_t$为零均值的白噪声序列,即$\epsilon_t \sim WN(0,\sigma_\epsilon^2)$;当期的随机干扰与过去的序列值无关,即$Ex_tx_s=0,s<t$。
当$\phi_0=0$时,其便为中心化的AR(p)模型,对于非中心化的AR(p)模型也可以通过使序列减去$\phi_0$转化为中心化AR(p)模型。
算子形式:
$$\phi(B)x_t=\epsilon_t$$
其中
$$\phi(B)=1-\phi_1B-\phi_2B^2-\cdots-\phi_pB^p$$
移动平均模型
另一种描述观察时间序列的重要模型叫做q阶移动平均模型MA(q)。在这类模型中,时间序列的现在值$x_t$是用白噪声扰动项的线性组合来表示的,即
$$x_t = \theta_0+\epsilon_t-\theta_1\epsilon_{t-1}-\theta_2\epsilon_{t-2}-\cdots-\theta_q\epsilon_{t-q}\tag{2}$$
其中包含两个限制条件:模型最高阶数为q,即$\theta_q \ne 0$;随机干扰序列$\epsilon_t$为零均值的白噪声序列,即$\epsilon_t \sim WN(0,\sigma_\epsilon^2)$。
中心化后的算子形式
$$x_t=\theta(B)\epsilon_t$$
其中
$$\theta(B)=1-\theta_1B-\theta_2B^2-\cdots-\theta_qB^q$$
自回归移动平均模型
时间序列的现在值$x_t$不仅与其以前时刻的自身值有关,而且还与其过去时刻进入系统的扰动存在一定依存关系,形成了ARMA(p,q)模型,即:
$$x_t = \mu+\phi_1x_{t-1}+\cdots+\phi_px_{t-p}+\epsilon_t-\theta_1\epsilon_{t-1}-\cdots-\theta_q\epsilon_{t-q}\tag{3}$$
模型的限制条件与AR(p)模型、MA(q)模型相同。
中心化后的算子形式:
$$\phi(B)x_t=\theta(B)\epsilon_t$$
$\tag{3}$式中含有$p+q+1$个未知参数$\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q,\sigma_\epsilon^2$,这些未知参数需要利用观测数据进行估计。
AR(p)模型可以视为ARMA(p,q)模型当q=0时的形式;MA(q)模型是ARMA(p,q)模型当p=0时的形式。
求和自回归移动平均模型
若$\{x_t,t \in T\}$在$d$阶差分后平稳,则满足$\tag{4}$式的模型称为自回归求和移动平均模型,简记为$ARIMA(p,d,q)$模型,其中心化后的算子形式为:
\begin{cases}
\phi(B)\nabla^dx_t=\theta(B)\epsilon_t \\
E(\epsilon_t)=0,Var(\epsilon_t)=\sigma_{\epsilon}^2,E(\epsilon_t]epsilon_s)=0,s\ne t\\
Ex_sx_t=0,\forall s<t\\
\end{cases}
其中$|B|\le 1$且$\phi(B)$与$\theta(B)$互质,$\phi_p\theta_q\ne 0$,$d$为差分阶数,$p$为自回归阶数,$q$为移动平均阶数,$\{\epsilon_t\}$为零均值白噪声序列。
由(4)式可知,$ARIMA$模型实质上就是差分运算与$ARMA$模型的组合。因此任何非平稳的时间序列若经过适当阶数的差分后平稳,就可以对差分后的序列进行$ARMA$模型拟合。目前$ARMA$模型的分析方法非常成熟,因此对差分后的平稳序列的分析也将非常简单可靠了。
模型选择
平稳非白噪声序列阶数识别
如果一个通过预处理的序列是平稳非白噪声序列,则可以对该序列进行建模。模型识别是时间序列建模的第一个阶段,是根据样本自相关系数和偏相关系数的性质来选择阶数适当的模型,也称为模型定阶。$ARMA(p,q)$模型定阶的基本原则如下表所示:
模型 | 自相关图 | 偏自相关图 |
AR(p) | 拖尾 | p阶截尾 |
MA(q) | q阶截尾 | 拖尾 |
ARMA(p,q) | 拖尾 | 拖尾 |
模型不适合 | 截尾 | 截尾 |
在实践中,由于样本的随机性,样本自相关系数和偏相关系数不会呈现理论截尾的完美情况,本应截尾处还会呈现出零值附近作小值震荡的情形,对于是否截尾及其结束判断通常依靠主观经验,不过可以利用倍标准差范围辅助判断。
取显著性水平$\alpha=0.05$,如果样本自相关系数和偏相关系数在最初的阶明显大于2倍标准差,而后几乎95%的系数都落在2倍标准差的范围内,且非零系数衰减为小值波动的过程非常突然,通常可视为k阶截尾;如果有超过5%样本自相关系数和偏相关系数落入2倍标准差的范围外,或者非零系数衰减为小值波动的过程比较缓慢或连续,通常视为拖尾。
差分阶数识别
对于非平稳的时间序列要通过差分变得平稳,还要识别差分阶数d和季节差分因子s,对于差分阶数的识别一般通过时序图和样本自相关系数来判别,见下表:
模型参数估计
选择好了拟合模型后,下一步就是利用序列值来确定该模型的口径,即估计模型中未知参数的值。采用不同的估计方法会带来不同的误差,其大小直接影响模型的精度。因此对于一个非中心化的$ARMA(p,q)$模型:
$$\phi(B)x_t = \mu + \theta(B)\epsilon_t,\epsilon_t\sim WN(0,\sigma_\epsilon^2)$$
共含有$p+q+2$个参数:$\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q,\mu,\sigma_\epsilon^2$。其中参数$\mu$是序列均值,样本均值是它的无偏估计,即:
$$\hat\mu=\overline x=\frac 1n\sum_{i=1}^nx_i$$
因此可以对模型进行中心化后估计剩余的$p+q+1$的参数,常用的估计方法有矩估计、最小二乘估计、极大似然估计、最小方差估计、最大熵估计等。对于$ARIMA$模型只需要将其差分后序列按照$ARMA$模型定阶估计参数后再整理即可。
下仅以最小二乘估计举例说明:
在$ARMA(p,q)$模型场合,记
$$\beta=(\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q)’$$
$$F_t(\beta)=\phi_1x_{t-1}+\cdots+\phi_px_{t-p}-\theta_1\epsilon_{t-1}-\cdots-\theta_q\epsilon_{t-q}$$
残差平方和为:
$$
Q(\beta)=\sum_{t=1}^n\epsilon_t^2=\sum_{t=1}^n(x_t-F_t(\beta))^2\\
=\sum_{t=1}^n(x_t-\phi_1x_{t-1}-\cdots-\phi_px_{t-p}+\theta_1\epsilon_{t-1}+\cdots+\theta_q\epsilon_{t-q})^2
$$
由于残差项$\epsilon_{t-i}$无法事先给出,因此可以利用逆函数将其化为$x_{t-i}$,$x_{t-i-1}$的线性组合形式,但是由于逆函数本身就是$\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q$的函数,因此最终得到的$Q(\beta)$对参数来说必然是非线性的,因此不能根据使$Q(\beta)$达到极小来求$\beta$,有必要采用求解非线性最小二乘问题。通常采取线性化的方法,将非线性函数展开后仅取线性部分来近似,然后使残差平方和最小。
最小二乘估计法是常用的最佳估计法,由于充分利用了序列值的信息,因而精度很高。
模型检验及优化
在对时间序列建模之后,还必须进行模型检验来判断拟合模型是否适当,需要检验两个内容:模型的平稳性和可逆性、模型适用性。
模型的平稳性和可逆性检验
模型的适用性检验主要包括两部分:模型的显著性检验和参数的显著性检验。
模型的显著性检验是检验模型的有效性。一个模型是否显著有效,主要看它对序列值提取的信息是否充分。一个好的拟合模型应该将序列值中的样本相关信息提取完全,即残差序列应该是纯粹由随机干扰产生的,即它应当是白噪声序列。反之,如果残差序列不是白噪声序列,那就意味着残差序列中还有相关信息可供提取,说明拟合模型不够有效,通常需要选择其他模型重新拟合。因此,模型的显著性检验也就是残差序列的白噪声检验。
参数的显著性检验就是要检验每一个未知参数是否显著非零,目的是使模型最精简。如果某个参数不显著,也就是说该参数所对应的自变量对因变量的影响不明显,本着节约的原则,应将该自变量从拟合模型中删除,最终模型将由一系列参数显著非零的自变量表示。
设$\beta=(\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q)’$为未知参数向量,$\beta_j$为其中第$j$个分量,$\hat\beta_j$为第$j$个分量的最小二乘估计,$Q(\beta)$为残差平方和,$X$为序列值向量,$(X’X)^{-1}=(\alpha_{ij})_{N\star N}$,检验假设:
$$H_0:\beta_j=0\leftrightarrow H_1:\beta_j\ne 0,\forall 1 \le j \le p+q$$
构造$t$检验统计两:
$$T=\sqrt{n-m}\frac{\hat\beta_j-\beta_j}{\sqrt{a_{jj}Q(\beta)}}\sim t(n-m)$$
当该检验统计量的$P$值小于$\frac{\alpha}{2}$或大于$1-\frac{\alpha}{2}$时,拒绝原假设,认为该参数显著;否则,认为该参数不显著。此时,应剔除不显著参数所对应的自变量重新拟合模型,构造出结构更精炼的拟合模型。
模型优化
能够通过检验的模型往往不止一个,最终选择哪个拟合模型可以遵循一定的准则。一个拟合模型的好坏可以从两方面考察:一方面是似然函数值,通常似然函数值越大说明拟合的效果越好;另一方面是模型中未知参数的个数,自然在精度没有显著变化时未知参数的个数越少越好。但是这二者是矛盾的,未知参数个数越多,说明模型中自变量越多,模型变化越灵活,拟合的精度越高,似然函数值越大。但未知参数越多,参数估计的难度就越大,估计的精度越差。因此,一个相对最优的拟合模型应该是拟合精度和未知参数个数的综合最优配置。常用的准则有AIC、BIC和SBC准则,可称为最小信息量检验。
AIC准则的思想是使拟合精度和参数个数的加权函数达到最小的模型被认为相对最优,AIC函数为:
AIC=-2ln(模型的极大似然函数值)+2(模型中未知参数的个数)
但是在样本容量趋于无穷大时,根据AIC准则选择的模型不收敛于真实模型,为了弥补AIC准则的不足,相继提出了BIC准则和SBC准则,其中SBC准则是很据Bayes理论得出的,其定义为:
SBC=-2ln(模型的极大似然函数值)+ln(n)(模型中未知参数的个数)
因为不可能比较所有模型的AIC和SBC函数值,因此只能在尽可能全面的范围内考察有限个模型的AIC和SBC函数值,认为其中AIC和SBC函数值达到最小的模型相对最优,可作为最终的拟合模型。
线性最小方差预测
在对平稳非白噪声序列做了模型识别、参数估计、模型检验及优化这一系列工作后,最终的目的是利用选定的拟合模型对随机序列的未来发展进行预测。时间序列法建立的模型必须满足平稳性条件和可逆性条件,不满足这两个条件的模型不能用来作为预测模型。目前对平稳序列最常用的预测方法是线性最小方差预测,设$\hat x_{t+l}$表示在t时刻向前预测$l$期得到预测值,$l$是预测长度,预测的准则是使得预测方差达到最小,即
$$Var_{\hat x(l)}[x_{t+l}-\hat x_{t+l}]=min\{Var_{\hat x(l)}[x_{t+l}-\hat x_{t+l}]\}$$
经推导可得$ARMA$模型的最小方差预测计算公式为
$$\hat x_{t+l} = \hat {\phi_1}[x_{t+l-1}]+\cdots+ \hat {\phi_p}[x_{t+l-p}]-\hat{\theta_1}[\epsilon_{t+l-1}]-\cdots-\hat{\theta_q}[\epsilon_{t+l-q}]$$
其中$[x_t]$,$[\epsilon_t]$分别为$x_t$,$\epsilon_t$条件期望的简写。在预测计算中,条件期望的取值如下:
$$
[x_{t+j}]=\begin{cases}
x_{t+j},j\le 0\\
\hat {x_{t+j}},j\gt 0
\end{cases}
$$
$$
[\epsilon_{t+j}]=\begin{cases}
0,j\gt 0\\
x_{t+j}-\hat x_{t+j},j\le 0
\end{cases}
$$
两式表明:现在或过去的序列观察值的条件期望就是本身,未来值的条件期望就是其预测值;现在或过去残差的条件期望是此残差的估计值,残差未来值的条件期望等于零。
时序模型参数的估计过程是复杂的非线性回归过程,从而导致时序模型阶数的确定与参数的估计非常复杂。实际应用中,将有关理论、方法与时序模型的参数估计相结合,发展了大量的参数估计方法模型参数估计的经典方法有矩法、非线性最小二乘法、极大似然估计法等。这些估算方法互有优缺点,如矩估计法计算简单,但精度较低,后两种估计法虽然精度较高,但计算过程繁琐,限制条件也很多。
基于ARIMA实现平安银行收盘价预测实战
加载库
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
import numpy as np import pandas as pd import tushare as ts import datetime import matplotlib.pylab as plt from matplotlib.pylab import style from arch.unitroot import ADF import statsmodels.api as sm import statsmodels.formula.api as smf import statsmodels.tsa.api as smt from statsmodels.tsa.stattools import adfuller from statsmodels.stats.diagnostic import acorr_ljungbox from statsmodels.graphics.api import qqplot import pylab as mpl mpl.rcParams['axes.unicode_minus'] = False # 避免符号‘-’为乱码 mpl.rcParams['font.sans-serif'] = ['SimHei']# 设置字体为黑体 |
获取数据
1 2 3 |
pro = ts.pro_api('在此键入你的tokey') df = pro.daily(ts_code='000001.SZ') df |
此时得到的数据如下所示:
ts_code | trade_date | open | high | low | close | pre_close | change | pct_chg | vol | amount | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 000001.SZ | 20220120 | 16.47 | 17.46 | 16.42 | 17.33 | 16.50 | 0.83 | 5.0303 | 3031193.77 | 5.195806e+06 |
1 | 000001.SZ | 20220119 | 16.54 | 16.69 | 16.36 | 16.50 | 16.52 | -0.02 | -0.1211 | 988391.85 | 1.630974e+06 |
2 | 000001.SZ | 20220118 | 16.25 | 16.54 | 16.16 | 16.52 | 16.22 | 0.30 | 1.8496 | 1152974.31 | 1.891095e+06 |
3 | 000001.SZ | 20220117 | 16.35 | 16.48 | 16.12 | 16.22 | 16.33 | -0.11 | -0.6736 | 1143689.43 | 1.855634e+06 |
4 | 000001.SZ | 20220114 | 17.00 | 17.02 | 16.30 | 16.33 | 16.98 | -0.65 | -3.8280 | 2103094.91 | 3.459951e+06 |
… | … | … | … | … | … | … | … | … | … | … | … |
4995 | 000001.SZ | 20001023 | 17.44 | 17.69 | 17.29 | 17.41 | 17.05 | 0.36 | 2.1100 | 101335.00 | 1.771737e+05 |
4996 | 000001.SZ | 20001020 | 16.80 | 17.18 | 16.77 | 17.05 | 16.79 | 0.26 | 1.5500 | 22719.00 | 3.846377e+04 |
4997 | 000001.SZ | 20001019 | 17.00 | 17.06 | 16.70 | 16.79 | 16.99 | -0.20 | -1.1800 | 16769.00 | 2.826112e+04 |
4998 | 000001.SZ | 20001018 | 17.00 | 17.25 | 16.90 | 16.99 | 16.93 | 0.06 | 0.3500 | 18161.00 | 3.097700e+04 |
4999 | 000001.SZ | 20001017 | 16.80 | 17.08 | 16.71 | 16.93 | 16.78 | 0.15 | 0.8900 | 11648.00 | 1.962019e+04 |
5000 rows × 11 columns
数据预处理
1 2 3 4 5 6 7 |
df['date'] = pd.to_datetime(df['trade_date'], format='%Y%m%d') df = df.iloc[179::-1,:] df["diff1"] = df["close"].diff(1).dropna() df["diff2"] = df["diff1"].diff(1).dropna() data = df.loc[:,["close",'diff1','diff2']] data.index = df['date'] data |
此时数据如下所示:
close | diff1 | diff2 | |
---|---|---|---|
date | |||
2021-04-28 | 23.35 | NaN | NaN |
2021-04-29 | 23.59 | 0.24 | NaN |
2021-04-30 | 23.29 | -0.30 | -0.54 |
2021-05-06 | 23.50 | 0.21 | 0.51 |
2021-05-07 | 24.05 | 0.55 | 0.34 |
… | … | … | … |
2022-01-14 | 16.33 | -0.65 | -0.63 |
2022-01-17 | 16.22 | -0.11 | 0.54 |
2022-01-18 | 16.52 | 0.30 | 0.41 |
2022-01-19 | 16.50 | -0.02 | -0.32 |
2022-01-20 | 17.33 | 0.83 | 0.85 |
180 rows × 3 columns
平稳性检验
时序图检验
1 |
data.plot(subplots=True, figsize=(18, 12),title="差分图") |
得到:
时序图检验全靠肉眼的判断和判断人的经验,不同的人看到同样的图形,很可能会给出不同的判断。因此我们需要一个更有说服力、更加客观的统计方法来帮助我们检验时间序列的平稳性,即单位根检验。
单位根检验
1 2 |
print("单位根检验:\n") print(ADF(data.diff1.dropna())) |
1 2 3 4 5 6 7 8 9 10 11 12 13 |
单位根检验: Augmented Dickey-Fuller Results ===================================== Test Statistic -6.274 P-value 0.000 Lags 4 ------------------------------------- Trend: Constant Critical Values: -3.47 (1%), -2.88 (5%), -2.58 (10%) Null Hypothesis: The process contains a unit root. Alternative Hypothesis: The process is weakly stationary. |
单位根检验:对其一阶差分进行单位根检验,得到:1%、%5、%10不同程度拒绝原假设的统计值和ADF Test result的比较,本数据中,P-value 为0,ADF Test result同时小于1%、5%、10%即说明很好地拒绝该假设,本数据中,ADF结果为-6.274,拒绝原假设,即一阶差分后数据是平稳的。
白噪声检验
判断序列是否为非白噪声序列
1 2 |
from statsmodels.stats.diagnostic import acorr_ljungbox acorr_ljungbox(data.diff1.dropna(), lags = [i for i in range(1,12)],boxpierce=True) |
acorr_ljungbox函数返回值:
lbvalue:测试的统计量
pvalue:基于卡方分布的p统计量
bpvalue:((optionsal), float or array) – 基于 Box-Pierce 的检验的p统计量
bppvalue:((optional), float or array) – 基于卡方分布下的Box-Pierce检验的p统计量
本次实验所得结果:
1 2 3 4 5 6 7 8 9 10 11 12 |
(array([ 0.60092286, 1.23084601, 1.97640996, 3.69327488, 9.65987062, 9.67409189, 9.67458357, 14.13452123, 18.69586261, 18.95205152, 18.96251431]), array([0.43822613, 0.54041225, 0.57731731, 0.44910481, 0.08546666, 0.13906476, 0.20777663, 0.07832585, 0.02790675, 0.04087594, 0.06176972]), array([ 0.59096281, 1.20696501, 1.93193327, 3.59188554, 9.32772896, 9.34132167, 9.3417889 , 13.55532116, 17.83945395, 18.07865796, 18.08836927]), array([0.44204705, 0.54690372, 0.58665369, 0.46404516, 0.09668478, 0.15527269, 0.22903941, 0.09411775, 0.03708251, 0.0536512 , 0.07955265])) |
通过P<α,拒绝原假设,故差分后的序列是平稳的非白噪声序列,可以进行下一步建模
模型定阶
现在我们已经得到一个平稳的时间序列,接来下就是选择合适的ARIMA模型,即ARIMA模型中合适的p,q。
第一步我们要先检查平稳时间序列的自相关图和偏自相关图。通过sm.graphics.tsa.plot_acf和sm.graphics.tsa.plot_pacf得到图形。
1 2 |
sm.graphics.tsa.plot_acf(data.iloc[1:,1]) sm.graphics.tsa.plot_pacf(data.iloc[1:,1]) |
从一阶差分序列的自相关图和偏自相关图可以发现:
- 自相关图拖尾或一阶截尾
- 偏自相关图一阶截尾,
所以我们可以建立ARIMA(1,1,0)、ARIMA(1,1,1)、ARIMA(0,1,1)模型。
模型优化
其中L是在该模型下的最大似然,n是数据数量,k是模型的变量个数。
1 2 3 4 5 6 |
arma_mod20 = sm.tsa.ARIMA(data["close"],(1,1,0)).fit() arma_mod30 = sm.tsa.ARIMA(data["close"],(0,1,1)).fit() arma_mod40 = sm.tsa.ARIMA(data["close"],(1,1,1)).fit() values = [[arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic],[arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic],[arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic]] df = pd.DataFrame(values,index=["AR(1,1,0)","MA(0,1,1)","ARMA(1,1,1)"],columns=["AIC","BIC","hqic"]) df |
结果如下:
AIC | BIC | hqic | |
---|---|---|---|
AR(1,1,0) | 217.207469 | 226.769626 | 221.084848 |
MA(0,1,1) | 217.115875 | 226.678032 | 220.993254 |
ARMA(1,1,1) | 218.735952 | 231.485496 | 223.905791 |
MA(0,1,1)的AIC为最小,因此我们选择此模型。
参数估计
1 2 3 4 |
from statsmodels.tsa.arima_model import ARIMA model = ARIMA(data["close"], order=(0,1,1)) result = model.fit() print(result.summary()) |
得到:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
ARIMA Model Results ============================================================================== Dep. Variable: D.close No. Observations: 179 Model: ARIMA(0, 1, 1) Log Likelihood -105.558 Method: css-mle S.D. of innovations 0.436 Date: Fri, 21 Jan 2022 AIC 217.116 Time: 15:54:52 BIC 226.678 Sample: 1 HQIC 220.993 ================================================================================= coef std err z P>|z| [0.025 0.975] --------------------------------------------------------------------------------- const -0.0332 0.035 -0.953 0.340 -0.101 0.035 ma.L1.D.close 0.0677 0.081 0.833 0.405 -0.092 0.227 Roots ============================================================================= Real Imaginary Modulus Frequency ----------------------------------------------------------------------------- MA.1 -14.7711 +0.0000j 14.7711 0.5000 ----------------------------------------------------------------------------- |
模型的显著性检验
1 2 3 4 |
resid = result.resid#残差 fig = plt.figure(figsize=(12,8)) ax = fig.add_subplot(111) fig = qqplot(resid, line='q', ax=ax, fit=True) |
qq图显示,我们看到红色的KDE线与N(0,1)平行,这是残留物正太分布的良好指标,说明残差序列是白噪声序列,模型的信息的提取充分
模型预测
1 2 3 4 |
data['pred'] = [None]*180 for i in range(100,180): pred = result.predict(i,dynamic=True, typ='levels') data.iloc[i,-1]=pred[i] |
参考文献
[1]魏来. ARIMA理论及其在网络环境下绿色农产品定价的应用[D].电子科技大学,2003.
[2]赵国顺. 基于时间序列分析的股票价格趋势预测研究[D].厦门大学,2009.
[3]汤岩. 时间序列分析的研究与应用[D].东北农业大学,2007.
[4]杨宇塬,张梅.基于ARIMA模型的股票价格实证分析[J].科技资讯,2021,19(29):121-123+127.DOI:10.16661/j.cnki.1672- 3791.2110-5042-7272.
[5]https://blog.csdn.net/sunnyxidian/article/details/92946542
[6]https://blog.csdn.net/weixin_41988628/article/details/83149849
[7]https://blog.csdn.net/qq_45176548/article/details/111504846#comments_14302892
暂无评论