# 4.1 自回归模型（Autoregressive) 自回归模型（AR）基于的基础是可以用历史预测未来，也就是时刻t的数据能用之前时刻的函数来表示。 自回归模型是很多人在拟合时间序列数据时最先尝试的模型，尤其是在对该数据没有额外的了解时。它的基本公式可以用以下公式表示， $$ y\_t = b\_0 + b\_1 \* y\_{t -1} + e\_t $$ 这个最简单的模型称为AR(1)，其中括号中的1表示时间间隔为1，也就是当前时刻的数据值的计算只考虑到前一个时刻的值。从公式可以看出，它和常规的线性回归十分类似，$$b\_0,b\_1$$分别代表截距项和回归系数，$$ e\_t$$表示t时刻的错误项，其中错误项均值为0，具有固定方差。这个公式表示的是用t-1时刻的时间序列值来拟合t时刻的时间序列值。 AR模型可以扩展到p个近邻时间值，此时称为AR(p)。 $$ y\_t = \phi\_0 + \phi\_1 \* y\_{t -1} + \phi\_2 \* y\_{t -2} + ... + \phi\_p \* y\_{t - p} + e\_t $$ 其中$$\phi$$表示一系列回归系数。 **python实战演练** 在实战中我们应当注意使用AR模型的两个前提假定： * 过去发生的数据能够用来预测未来的数据（也就是数据之间是非独立的） * 时间序列数据时平稳的 如果数据集是不平稳的，需要使用一些操作将数据去除趋势项和季节项使其变得平稳。 *冷知识：平稳性非为强平稳性和弱平稳性，其中强平稳性要求数据的分布不随时间变化，而弱平稳性仅仅要求数据的一阶距和二阶矩（均值和方差）不随时间变化。* ```python import matplotlib.pyplot as plt import pandas as pd from statsmodels.tsa.ar_model import AutoReg from statsmodels.tsa.seasonal import seasonal_decompose %matplotlib inline ``` ```python # 导入销售数据，绘制原始图像 sales_data = pd.read_csv('data/retail_sales.csv') sales_data['date']=pd.to_datetime(sales_data['date']) sales_data.set_index('date', inplace=True) sales_data.plot() ``` ![](https://3993849477-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-Mhv1ams1lvf_fMUn_is%2F-MiFWhiLoM9kclvUqRe1%2F-MiFWorUzgdts1WFS-Br%2F4_1.png?alt=media\&token=0d782f6d-f8f9-4ccb-8e8d-88ed746a4422) 我们提到使用AR模型的两个前提假设，相关性和平稳性。 因此先来检验相关性。 检验相关性有两种方法 * pandas提供了autocorrelation\_plot方法用于检验整体自相关性 * 使用pandas的lag\_plot方法检查单个时间间隔的相关性 ```python pd.plotting.autocorrelation_plot(sales_data['sales']) # 从自相关性图中可以看到在lag<=12时，acf值在临界点以外，表现出极强的相关性。 ``` ![](https://3993849477-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-Mhv1ams1lvf_fMUn_is%2F-MiFWhiLoM9kclvUqRe1%2F-MiFWxkhEZHEmXV8wVvr%2F4_2.png?alt=media\&token=d4797157-0720-4d85-819f-3324731e71b5) ```python # 将上图中相关性最强的lag=12单独画散点图，观察相关性 pd.plotting.lag_plot(sales_data['sales'],lag=12) ``` ![](https://3993849477-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-Mhv1ams1lvf_fMUn_is%2F-MiFWhiLoM9kclvUqRe1%2F-MiFX10_lCea2UIHFQh4%2F4_3.png?alt=media\&token=7d3b649d-e526-4d35-8bfd-7691badec5ee) 第二步检查平稳性，一个最快捷的方法之一是使用statsmodels中的seasonal\_decompose方法进行趋势项和季节项的分解。 ```python # 数据集存在明显趋势项（图二单调递增）和季节项（图三周期性变化） decomposed = seasonal_decompose(sales_data['sales'], model='additive') x = decomposed.plot() ``` ![](https://3993849477-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-Mhv1ams1lvf_fMUn_is%2F-MiFWhiLoM9kclvUqRe1%2F-MiFX7xRNfCeyVEUdTNP%2F4_4.png?alt=media\&token=3e22e981-1459-4905-9c16-abcfcb287454) 幸好statsmodel包的AutoReg方法增加了对趋势和季节项特征的处理，直接使用该方法即可 ```python # 划分训练集合测试集 X = sales_data['sales'] train_data = X[1:len(X)-12] test_data = X[len(X)-12:] #以最后十二个点作为待预测值 # 训练AR模型 model = AutoReg(train_data,lags=15,missing='drop',trend='t',seasonal=True) model_fitted = model.fit() ``` ```python # 查看模型结果 model_fitted.summary() ``` | AutoReg Model Results | | | | | --------------------- | ----------------- | ------------------- | -------- | | Dep. Variable: | sales | No. Observations: | 59 | | Model: | Seas. AutoReg(15) | Log Likelihood | -422.998 | | Method: | Conditional MLE | S.D. of innovations | 3621.518 | | Date: | Tue, 22 Jun 2021 | AIC | 17.707 | | Time: | 22:47:53 | BIC | 18.883 | | Sample: | 02-01-2011 | HQIC | 18.144 | | | - 09-01-2014 | | | | | coef | std err | z | P>\|z\| | \[0.025 | 0.975] | | ----------- | --------- | -------- | ------ | ------- | --------- | -------- | | trend | 839.9681 | 362.378 | 2.318 | 0.020 | 129.721 | 1550.216 | | seasonal.0 | 2.213e+05 | 8.34e+04 | 2.653 | 0.008 | 5.78e+04 | 3.85e+05 | | seasonal.1 | 2.808e+05 | 8.08e+04 | 3.476 | 0.001 | 1.22e+05 | 4.39e+05 | | seasonal.2 | 1.924e+05 | 8.7e+04 | 2.212 | 0.027 | 2.19e+04 | 3.63e+05 | | seasonal.3 | 1.689e+05 | 8.69e+04 | 1.945 | 0.052 | -1331.228 | 3.39e+05 | | seasonal.4 | 2.263e+05 | 8.29e+04 | 2.728 | 0.006 | 6.37e+04 | 3.89e+05 | | seasonal.5 | 2.335e+05 | 7.98e+04 | 2.926 | 0.003 | 7.71e+04 | 3.9e+05 | | seasonal.6 | 2.09e+05 | 8.26e+04 | 2.530 | 0.011 | 4.71e+04 | 3.71e+05 | | seasonal.7 | 2.318e+05 | 8.32e+04 | 2.785 | 0.005 | 6.87e+04 | 3.95e+05 | | seasonal.8 | 2.33e+05 | 8.41e+04 | 2.771 | 0.006 | 6.82e+04 | 3.98e+05 | | seasonal.9 | 2.224e+05 | 8.4e+04 | 2.648 | 0.008 | 5.78e+04 | 3.87e+05 | | seasonal.10 | 1.901e+05 | 8.46e+04 | 2.246 | 0.025 | 2.42e+04 | 3.56e+05 | | seasonal.11 | 2.123e+05 | 8.59e+04 | 2.473 | 0.013 | 4.4e+04 | 3.81e+05 | | sales.L1 | 0.2604 | 0.155 | 1.685 | 0.092 | -0.042 | 0.563 | | sales.L2 | 0.1237 | 0.158 | 0.785 | 0.433 | -0.185 | 0.433 | | sales.L3 | 0.0379 | 0.150 | 0.252 | 0.801 | -0.256 | 0.332 | | sales.L4 | -0.2515 | 0.149 | -1.691 | 0.091 | -0.543 | 0.040 | | sales.L5 | 0.2431 | 0.163 | 1.496 | 0.135 | -0.075 | 0.562 | | sales.L6 | -0.1179 | 0.163 | -0.722 | 0.470 | -0.438 | 0.202 | | sales.L7 | -0.1311 | 0.164 | -0.799 | 0.424 | -0.453 | 0.190 | | sales.L8 | -0.0212 | 0.196 | -0.108 | 0.914 | -0.406 | 0.363 | | sales.L9 | 0.2763 | 0.201 | 1.374 | 0.169 | -0.118 | 0.670 | | sales.L10 | 0.0443 | 0.200 | 0.222 | 0.825 | -0.347 | 0.436 | | sales.L11 | 0.1980 | 0.203 | 0.975 | 0.330 | -0.200 | 0.596 | | sales.L12 | 0.1034 | 0.202 | 0.512 | 0.609 | -0.292 | 0.499 | | sales.L13 | -0.4173 | 0.206 | -2.029 | 0.042 | -0.820 | -0.014 | | sales.L14 | 0.0320 | 0.192 | 0.166 | 0.868 | -0.345 | 0.409 | | sales.L15 | 0.0085 | 0.206 | 0.041 | 0.967 | -0.396 | 0.413 | | Roots | | | | | | ----- | ------- | --------- | ------- | --------- | | | Real | Imaginary | Modulus | Frequency | | AR.1 | -0.9178 | -0.4910j | 1.0409 | -0.4218 | | AR.2 | -0.9178 | +0.4910j | 1.0409 | 0.4218 | | AR.3 | -1.1752 | -0.0000j | 1.1752 | -0.5000 | | AR.4 | -0.6070 | -0.8247j | 1.0241 | -0.3510 | | AR.5 | -0.6070 | +0.8247j | 1.0241 | 0.3510 | | AR.6 | -0.0849 | -1.1181j | 1.1213 | -0.2621 | | AR.7 | -0.0849 | +1.1181j | 1.1213 | 0.2621 | | AR.8 | 0.3804 | -0.9597j | 1.0323 | -0.1899 | | AR.9 | 0.3804 | +0.9597j | 1.0323 | 0.1899 | | AR.10 | 0.8226 | -0.5979j | 1.0170 | -0.1000 | | AR.11 | 0.8226 | +0.5979j | 1.0170 | 0.1000 | | AR.12 | 1.1469 | -0.1677j | 1.1591 | -0.0231 | | AR.13 | 1.1469 | +0.1677j | 1.1591 | 0.0231 | | AR.14 | 5.1284 | -0.0000j | 5.1284 | -0.0000 | | AR.15 | -9.1838 | -0.0000j | 9.1838 | -0.5000 | ```python # 预测测试集 predictions = model_fitted.predict( start=len(train_data), end=len(train_data) + len(test_data)-1, dynamic=False) # dynamic参数表示是否用预测值动态预测下一个时刻的值 # 比较真实值和预测值 compare_df = pd.concat( [sales_data['sales'].tail(12), predictions], axis=1).rename( columns={'sales': 'actual', 0:'predicted'}) compare_df.plot() ``` ![](https://3993849477-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-Mhv1ams1lvf_fMUn_is%2F-MiFWhiLoM9kclvUqRe1%2F-MiFXDxEivc4_A9HRFun%2F4_5.png?alt=media\&token=e210afca-5735-4a52-8acd-e22afa364741)