• 2020-09-16 21:32:36

sarima模型

In this tutorial I will show you how to model a seasonal time series through a SARIMA model.

在本教程中，我将向您展示如何通过SARIMA模型对季节时间序列进行建模。

Here you can download the Jupyter notebook of the code described in this tutorial.

在这里，您可以下载本教程中描述的代码的Jupyter笔记本。

# 入门 (Getting Started)

## 将数据集转换为时间序列 (Convert the dataset into a time series)

In this example we will use the number of tourist arrivals to Italy. Data are extracted from the European Statistics: Annual Data on Tourism Industries. Firstly, we import the dataset related to foreign tourists arrivals in Italy from 2012 to 2019 October and then we convert it into a time series.

在此示例中，我们将使用前往意大利的游客人数。 数据摘自《 欧洲统计：旅游业年度数据》 。 首先，我们导入与2012年至2019年10月在意大利入境的外国游客有关的数据集，然后将其转换为时间序列。

In order to perform the conversion to time series, two steps are needed:

为了执行到时间序列的转换，需要两个步骤：

• the column containing dates must be converted to datetime. This can be done through the function to_datetime(), which converts a string into a datetime.

包含日期的列必须转换为datetime。 这可以通过函数to_datetime()完成，该函数将字符串转换为日期时间。

• set the index of the dataframe to the column containing dates. This can be done through the function set_index() applied to the dataframe.

将数据框的索引设置为包含日期的列。 这可以通过将函数set_index()应用于数据set_index()来完成。

import pandas as pddf = pd.read_csv('../sources/IT_tourists_arrivals.csv')df['date'] = pd.to_datetime(df['date'])df = df[df['date'] > '2012-01-01']df.set_index('date', inplace=True)

We can get some useful statistics related to the time series through the describe() function.

我们可以通过describe()函数获得一些与时间序列有关的有用统计信息。

df.describe()

# 初步分析 (Preliminary analysis)

## 绘制时间序列以检查季节性 (Plot the time series to check the seasonality)

The preliminary analysis involves a visual analysis of the time series, in order to understand its general trend and behaviour. Firstly, we create the time series and we store it in the variable ts.

初步分析包括对时间序列的可视化分析，以便了解其总体趋势和行为。 首先，我们创建时间序列并将其存储在变量ts

ts = df['value']

Then, we plot the ts trend. We use the matplotlib library provided by Python.

然后，我们绘制ts趋势。 我们使用Python提供的matplotlib库。

import matplotlib.pylab as pltplt.plot(ts)plt.ylabel('Total Number of Tourists Arrivals')plt.grid()plt.tight_layout()plt.savefig('plots/IT_tourists_arrivals.png')plt.show()

# 计算模型的参数 (Calculate the parameters for the model)

## 调整模型 (Tune the model)

We build a SARIMA model to represent the time series. SARIMA is an acronym for Seasonal AutoRegressive Integrated Moving Average. It is composed of two models AR and MA. The model is defined by three parameters:

我们建立了一个SARIMA模型来表示时间序列。 SARIMA是“季节性自回归综合移动平均线”的首字母缩写。 它由两个模型AR和MA组成。 该模型由三个参数定义：

• d = degree of first differencing involved

d =涉及的一次微分的程度
• p = order of the AR part

p = AR部分的顺序
• q = order of the moving average part.

q =移动平均线部分的阶数。

The value of p can be determined through the ACF plot, which shows the autocorrelations which measure the relationship between an observation and its previous one. The value of d is the order of integration and can be calculated as the number of transformations needed to make the time series stationary. The value of q can be determined through the PACF plot.

可以通过ACF图确定p的值，该图显示了测量一个观测值与其前一个观测值之间的关系的自相关。 d的值是积分的阶数，可以计算为使时间序列平稳所需的变换次数。 q的值可以通过PACF图确定。

In order determine the value of d, we can perform the Dickey-Fuller test, which is able to verify whether a time series is stationary or not. We can use the adfuller class, contained in the statsmodels library. We define a function, called test_stationarity(), which returns True, if the time series is positive, False otherwise.

为了确定d的值，我们可以执行Dickey-Fuller测试，该测试可以验证时间序列是否稳定。 我们可以使用statsmodels库中包含的adfuller类。 我们定义了一个名为test_stationarity()的函数，如果时间序列为正，则返回True，否则返回False。

from statsmodels.tsa.stattools import adfullerdef test_stationarity(timeseries):    dftest = adfuller(timeseries, autolag='AIC')    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])    for key,value in dftest[4].items():        dfoutput['Critical Value (%s)'%key] = value    critical_value = dftest[4]['5%']    test_statistic = dftest[0]    alpha = 1e-3    pvalue = dftest[1]    if pvalue < alpha and test_statistic < critical_value:  # null hypothesis: x is non stationary        print("X is stationary")        return True    else:        print("X is not stationary")        return False

We transform the time series through the diff() function as many times as the time series becomes stationary.

当时间序列变得平稳时，我们通过diff()函数对时间序列进行了多次转换。

ts_diff = pd.Series(ts)d = 0while test_stationarity(ts_diff) is False:    ts_diff = ts_diff.diff().dropna()    d = d + 1

In order to calculate the value of p and q, we can plot the ACF and PACF graphs, respectively. We can use the plot_acf() and plot_pacf() functions available in the statsmodels library. The value of p corresponds to the maximum value in the ACF graph external to the confidence intervals (shown in light blue). In our case, che correct value of p = 9.

为了计算p和q的值，我们可以分别绘制ACF和PACF图。 我们可以使用statsmodels库中可用的plot_acf()plot_pacf()函数。 p的值对应于置信区间外部的ACF图中的最大值(以浅蓝色显示)。 在我们的情况下，p = 9的正确值。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacfplot_acf(ts_trend, lags =12)plt.savefig('plots/acf.png')plt.show()

Similarly, the value of q corresponds to the maximum value in the PACF graph external to the confidence intervals (shown in light blue). In our case, che correct value of q = 1.

类似地，q的值对应于置信区间之外的PACF图中的最大值(以浅蓝色显示)。 在我们的情况下，q等于1的正确值。

plot_pacf(ts_trend, lags =12)plt.savefig('plots/pacf.png')plt.show()

# 建立SARIMA模型 (Build the SARIMA model)

## 如何训练SARIMA模型 (How to train the SARIMA model)

Now we are ready to build the SARIMA model. We can use the SARIMAX class provided by the statsmodels library. We fit the model and get the prediction through the get_prediction() function. We can retrieve also the confidence intervals through the conf_int() function.

现在，我们准备构建SARIMA模型。 我们可以使用statsmodels库提供的SARIMAX类。 我们拟合模型并通过get_prediction()函数获得预测。 我们还可以通过conf_int()函数检索置信区间。

from statsmodels.tsa.statespace.sarimax import SARIMAXp = 9q = 1model = SARIMAX(ts, order=(p,d,q))model_fit = model.fit(disp=1,solver='powell')fcast = model_fit.get_prediction(start=1, end=len(ts))ts_p = fcast.predicted_meants_ci = fcast.conf_int()

We plot results.

我们绘制结果。

plt.plot(ts_p,label='prediction')plt.plot(ts,color='red',label='actual')plt.fill_between(ts_ci.index[1:],                ts_ci.iloc[1:, 0],                ts_ci.iloc[1:, 1], color='k', alpha=.2)plt.ylabel('Total Number of Tourists Arrivals')plt.legend()plt.tight_layout()plt.grid()plt.savefig('plots/IT_trend_prediction.png')plt.show()

# 计算一些统计 (Calculate some statistics)

## 检查模型的性能 (Check the performance of the model)

Finally, we can calculate some statistics to evaluate the performance of the model. We calculate the Pearson’s coefficient, through the pearsonr() function provided by the scipy library.

最后，我们可以计算一些统计信息以评估模型的性能。 我们通过scipy库提供的pearsonr()函数来计算Pearson系数。

from scipy import statsstats.pearsonr(ts_trend_p[1:], ts[1:])

We also calculate the R squared metrics.

我们还计算R平方指标。

residuals = ts - ts_trend_pss_res = np.sum(residuals**2)ss_tot = np.sum((ts-np.mean(ts))**2)r_squared = 1 - (ss_res / ss_tot)r_squared

# 学过的知识 (Lesson learned)

In this tutorial I have illustrated how to model a time series through a SARIMA model. Summarising, you should follow the following steps:

在本教程中，我说明了如何通过SARIMA模型对时间序列进行建模。 总结一下，您应该遵循以下步骤：

• convert your data frame into a time series

将您的数据帧转换为时间序列
• calculate the values of p, d and q to tune the SARIMA model

计算p，d和q的值以调整SARIMA模型
• build the SARIMA model with the calculated p, d and q values

使用计算的p，d和q值构建SARIMA模型
• test the performance of the model.

测试模型的性能。

An improvement of the proposed model could involve the splitting of the time series in two parts: training and test sets. Usually, the training set is used to fit the model while the test set is used to calculate the performance of the model.

所提出模型的改进可能涉及将时间序列分为两部分：训练和测试集。 通常，训练集用于拟合模型，而测试集用于计算模型的性能。

# 参考书目 (Bibliography)

sarima模型

更多相关内容
• ARMA/ARIMA/SARIMA模型预测 a基本原理： 这三种模型都是用来预测时序性数据。其中ARIMA和SARIMA是由ARMA模型演变过来的，而ARMA是由AR模型(自回归模型)和MA模型(移动平均模型)组合出来的（AR模型和MA模型会在下文...

# ARMA/ARIMA/SARIMA模型预测

## a基本原理：

这三种模型都是用来预测时序性数据。其中ARIMA和SARIMA是由ARMA模型演变过来的，而ARMA是由AR模型(自回归模型)和MA模型(移动平均模型)组合出来的（AR模型和MA模型会在下文讲解）。

假设有一个时间节点t。AR做的事是利用t之前的随机变量来预测t之后随机变量；MA做的事是利用t之前（包括t）的随机误差项和滞后误差项，来形成一个误差项模型。

AR模型可以理解成ARMA模型的一种特殊的形式，相对而言ARMA功能更强大。

ARMA：针对弱平稳、宽平稳时间序列分析

ARIMA：针对非平稳、非周期性时间序列分析

SARIMA：针对非平稳、周期性时间序列分析

## b模型原理：

### AR模型(自回归模型)

自回归模型（Autoregressive Model）是用自身做回归变量的过程，即利用历史时序数据值的线性组合来预测当前时刻点的线性回归模型，它是时间序列中的一种常见形式。

### MA模型(移动平均模型)

MA模型和AR模型大同小异。移动平均模型(moving average model)使用历史白噪声的线性组合来预测当前时刻点的线性回归模型。与AR最大的不同之处在于，AR模型中历史白噪声的影响是间接影响当前预测值的（通过影响历史时序值）。

### ARMA模型(自回归移动平均模型)

ARMA（Auto-Regressive and Moving Average Model）模型顾名思义就是将AR模型和MA模型结合起来，即历史随机变量和历史白噪声结合起来，其模型可以看成是AR模型+MA模型。

### ARIMA模型(整合自回归移动平均模型)

ARIMA（Autoregressive Integrated Moving Average model）模型就是在ARMA模型的基础上多了一个差分数据处理。该模型的主要功能是将不平稳数据进行d次差分形成一个平稳的时间序列数据，然后采用ARMA模型。模型记为ARIMA(p,d,q),其中p和q就是ARMA(p,q)模型里面的p和q，而d是差分的次数（阶数）。

设当前时间为t，那么先将[t-p,t-1]历史随机变量和[t-q,t-1]历史白噪声进行d次差分形成平稳的时间序列数据，再将数据进行ARMA（p，q），则称该模型为（p,d,q）阶整合自回归移动平均模型，其数学模型和ARMA的模型一致。

### SARIMA模型(季节性整合自回归移动平均模型)

SARIMA(Seasonal Autoregressive Integrated Moving Average)模型在非稳的数据上多了一些对有周期性特征数据的处理。其主要原理是来源于外界对数据周期性的观察，得到季节的长度（s）、季节自回归的阶数（P，取值由PACF判断）、季节移动平均的阶数（Q，取值由ACF判断）季节差分的阶数（D，一般就是0或1），再建立关于季节性的模型，将该模型和ARIMA模型结合，就是SARIMA模型。

## d如何选择p、q值？

1. 利用自身相关函数（ACF）和偏自身相关函数（PACF）画出相应的图。利用图中的趋势 变换情况判断该趋势属于截尾还是拖尾情况，再通过图像判断其是在第几阶拖尾/截尾 的。在自相关函数（ACF）第p阶拖尾，偏自相关函数（PACF）第q阶拖尾。

AR（p）在自相关函数（ACF）第p阶拖尾，偏自相关函数（PACF）第q阶截尾。

MA（q）在自偏自相关函数（PACF）第q阶拖尾，自相关函数（ACF）第阶截尾。

ARMA(p,q)在自相关函数（ACF）第p阶拖尾，偏自相关函数（PACF）第q阶拖尾。

注：模型计算的时候有时还需要考虑一下p+1、p-1、q+1、q-1这样的情况。ACF和PACF 判断只是一种近似判断，所以还需要考虑一下几种情况差不多的概率。

2.第二种方法相对简单，因为通常状态下ｐ和ｑ的阶数不会太大（况且要是ｐ和ｑ的阶数 太大的话ARMA模型会很慢），所以可以暴力枚举所有的ｐ和ｑ的组合，找出ａｉｃ （Akaike information criterion）最小的情况，即为最优的ｐ和ｑ。

## f实例参考：

AR模型(ar_model.AR)

# -*- coding: utf-8 -*-
'''
自回归模型。预测澳大利亚墨尔本市最低日常温度
'''

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.ar_model import AR
from sklearn.metrics import mean_squared_error

path = r'csv/daily-minimum-temperatures-in-me.csv'  # csv路径

def getData(path):  # 获得训练数据和测试数据
X = series.values
train, test = X[:len(X) - 7], X[len(X) - 7:]  # 最后一周为测试数据，其余的为训练数据
return train, test

def getModel(X):  # 模型训练
try:
model = AR(endog=X)
return model.fit()
except ValueError:
print('参数错误')
return None

def prediction(model_fit, train, test):  # 预测
return model_fit.predict(start=len(train), end=len(train) + len(test) - 1, dynamic=False)  # 要输入预测开始和结束的位置

if __name__ == '__main__':
train, test = getData(path)
model_fit = getModel(train[:, 1])
print('常数项+模型参数:{}'.format(model_fit.params))

predict = prediction(model_fit, train, test)
for i in range(len(predict)):
print('预测值:{}，期望值:{}'.format(predict[i], test[:, 1][i]))

MSE = mean_squared_error(test[:, 1], predict)  # 输出均方误差
plt.plot(test[:, 1], color='blue', label='expect')
plt.plot(predict, color='red', label='predict')
plt.legend()
plt.show()

运行结果如下

数据在信息已删除

### ARMA模型

'''
自回归滑动平均模型
'''
from statsmodels.tsa.arima_model import ARMA
from itertools import product

def myARMA(data):
p = range(0, 9)
q = range(0, 9)
parameters = list(product(p, q))  # 生成(p,q)从(0,0)到(9,9)的枚举
best_aic = float('inf')
result = None
for param in parameters:
try:
model = ARMA(endog=data, order=(param[0], param[1])).fit()
except ValueError:
print("参数错误：", param)
continue
aic = model.aic
if aic < best_aic:  # 选取最优的aic
best_aic = model.aic
result = (model, param)
return result

## h参考文献：

数据分析——时间序列分析模型(AR,MA,ARMA,ARIMA)

数据分析——时间序列分析模型(AR,MA,ARMA,ARIMA)_人的抱怨源自，自我无能的愤怒-CSDN博客_时间序列分析模型

知乎 移动平均模型 8.4 移动平均模型 - 知乎

知乎 自回归模型 自回归模型 - 知乎

自回归模型（AR Model）自回归模型（AR Model）_shigangzwy的博客-CSDN博客_自回归模型

为什么常常用AR模型代替ARMA模型

为什么常常用AR模型代替ARMA模型_deco1515的专栏-CSDN博客

时间序列预测基础教程系列(7)_如何用自回归模型(AR)预测时间序列预测(Python)

时间序列预测基础教程系列(7)_如何用自回归模型(AR)预测时间序列预测(Python)_日拱一卒-CSDN博客_python自回归模型

statsmodels AR官方文档

https://www.statsmodels.org/stable/generated/statsmodels.tsa.ar_model.AR.html?highlight=ar#statsmodels.tsa.ar_model.AR

statsmodels ARMA官方文档

https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima_model.ARMA.html?highlight=arma#statsmodels.tsa.arima_model.ARMA

知乎 AR、MA及ARMA模型AR、MA及ARMA模型 - 知乎

如何确定ARIMA模型中参数p、d、q 如何确定ARIMA模型中参数p、d、q | 码农网

基于sarima模型的分析及预测 基于sarima模型的分析及预测 - 豆丁网

python 时间序列预测 —— SARIMA python 时间序列预测 —— SARIMA_颹蕭蕭-CSDN博客_python时间序列预测

展开全文
• 文章目录1、移动平均moving average方法weighted average方法2、指数平滑单指数平滑 exponential_smoothing双指数平滑三指数平滑 Triple exponential smoothing3、平稳性以及时间序列建模SARIMA模型4、时间序列的...

# 1、移动平均

## moving average方法

moving average方法不适合我们进行长期的预测（为了预测下一个值，我们需要实际观察的上一个值）。但是 移动平均数 还有另一个应用场景，即对原始的时间序列数据进行平滑处理，以找到数据的变化趋势
pandas 提供了一个实现接口 DataFrame.rolling(window).mean() ，滑动窗口 window 的值越大，意味着变化趋势将会越平滑，对于那些噪音点很多的数据集（尤其是金融数据），使用 pandas 的这个接口，有助于探测到数据中存在的共性（common pattern）
由于数据采集gap为10min，接下来我们调整window大小，看看分别有什么效果。
window_size = 6 -> 1h:

window_size = 24 -> 4h:

window_size = 144 -> 24h:

这里我们针对window_size = 6加上置信区间（scale=1.96）

我们可以基于 moving average 创建一个简单的异常检测系统（即如果数据点在置信区间之外，则认为是异常值
这里我们进行数据造假+扩大置信区间（scale=2.96）：

rtt_anomaly = df.rtt.copy()
rtt_anomaly.iloc[-20] = rtt_anomaly.iloc[-20] * 0.5 # say we have 80% drop of ads


## weighted average方法

上面提到了用 移动平均值对原始数据做平滑处理，接下来要说的是加权平均值，它是对上面 移动平均值 的简单改良。
也就是说，前面 k 个观测数据的值，不再是直接求和再取平均值，而是计算它们的加权和（权重和为1）。通常来说，时间距离越近的观测点，权重越大。数学表达式为：

def weighted_average(series, weights):
"""
Calculate weighter average on series
"""
result = 0.0
weights.reverse()
for n in range(len(weights)):
result += series.iloc[-n-1] * weights[n]
return float(result)

weighted_average(ads, [0.6, 0.3, 0.1])


# 2、指数平滑

## 单指数平滑 exponential_smoothing

如果不用上面提到的， /每次对当前序列中的前k个数/ 的加权平均值作为模型的预测值，而是直接对 /目前所有的已观测数据/ 进行加权处理，并且每一个数据点的权重，呈指数形式递减。
这个就是指数平滑的策略，具体怎么做的呢？一个简单的数学式为：

式子中的 α \alphaα 表示平滑因子，它定义我们“遗忘”当前真实观测值的速度有多快。α \alphaα 越小，表示当前真实观测值的影响力越小，而前一个模型预测值的影响力越大，最终得到的时间序列将会越平滑。
下面是选择不同的权重得到的曲线：

1. 单指数平滑的特点： 能够追踪数据变化。预测过程中，添加了最新的样本数据之后，新数据逐渐取代老数据的地位，最终老数据被淘汰。
2. 单指数平滑的局限性： 第一，预测值不能反映趋势变动、季节波动等有规律的变动；第二，这个方法多适用于短期预测，而不适合中长期的预测；第三，由于预测值是历史数据的均值，因此与实际序列相比，有滞后的现象。
3. 单指数平滑的系数： EViews提供两种确定指数平滑系数的方法：自动给定和人工确定。一般来说，如果序列变化比较平缓，平滑系数值应该比较小，比如小于0.l；如果序列变化比较剧烈，平滑系数值可以取得大一些，如0.3～0.5。若平滑系数值大于0.5才能跟上序列的变化，表明序列有很强的趋势，不能采用一次指数平滑进行预测。

## 双指数平滑

单指数平滑在产生新的序列的时候，考虑了前面的 K 条历史数据，但是仅仅考虑其静态值，即没有考虑时间序列当前的变化趋势。
如果当前的时间序列数据处于上升趋势，那么当我们对明天的数据进行预测的时候，就不应该仅仅是对历史数据进行”平均“，还应考虑到当前数据变化的上升趋势。同时考虑历史平均和变化趋势，这个就是我们的双指数平滑法
通过 序列分解法 (series decomposition)，我们可以得到两个分量，一个叫 intercept (also, level) ℓ \ellℓ ，另一个叫 trend (also, slope，斜率) b bb. 我们根据前面学习的方法，知道了如何预测 intercept （截距，即序列数据的期望值），我们可以将同样的指数平滑法应用到 trend (趋势)上。时间序列未来变化的方向取决于先前加权的变化。

在不同平滑因子的组合下的时序图:

调整两个参数 α \alphaα 和 β \betaβ 。前者决定时间序列数据自身变化趋势的平滑程度，后者决定趋势的平滑程度

## 三指数平滑 Triple exponential smoothing

三指数平滑，也叫 Holt-Winters 平滑，与前两种平滑方法相比，我们这次多考虑了一个因素，seasonality （季节性）。这其实也意味着，如果我们的时间序列数据，不存在季节性变化，就不适合用三指数平滑
模型中的 季节性 分量，用来解释 截距趋势 的重复变化，并且由季节长度来描述，也就是变化重复的周期来描述。
对于一个周期内的每一个观测点，都有一个单独的组成部分。

我们根据常识可知，这个数据集中，存在一个明显的季节性变化，变化周期为24小时，因此我们设置 slen = 24*6 = 144 :
Holt-Winters 模型以及其他指数平滑模型中，对平滑参数的大小有一个限制，每个参数都在0到1之间。因此我们必须选择支持模型参数约束的最优化算法，在这里，我们使用 Truncated Newton conjugate gradient (截断牛顿共轭梯度法)

# 3、平稳性以及时间序列建模

在我们开始建模之前，我们需要提到时间序列的一个重要特性，如平稳性 (stationarity)。
我们称一个时间序列是平稳的，是指它不会随时间而改变其统计特性，即平均值和方差不会随时间而改变。
因为现在大多数的时间序列模型，或多或少都是基于未来序列与目前已观测到的序列数据有着相同的统计特性(均值、方差等) 的假设。也就是说，如果原始序列（已观测序列）是不平稳的，那么我们现有模型的预测结果，就可能会出错。
糟糕的是，我们在教科书之外所接触到的时间序列数据，大多都是不平稳的，不过还好，我们有办法把它改变成平稳分布。
如果我们可以用一阶差分从非平稳序列中得到一个平稳序列，我们称这个非平稳序列为一阶积分。我们可以使用不同的方法来对抗非平稳性，如 d阶差分、趋势和季节性消除、平滑处理，也可以使用像box cox或对数这样的转换。

## SARIMA模型

1、绘制时间序列图、ACF 图和 PACF

Autocorrelation Function (ACF)自相关函数，指任意时间 t（t=1,2,3…n)的 序列值X^t^ 与其自身的滞后（这里取滞后一阶，即lag=1）值X^t-1^之间的线性关系。
这个图看不太懂，先试试其他的方法

# 4、时间序列的（非）线性模型

## 时间序列的滞后值

将时间序列来回移动 n 步，我们能得到一个特征，其中时序的当前值和其t-n时刻的值对齐。如果我们移动1时差，并训练模型预测未来，那么模型将能够提前预测1步。增加时差，比如，增加到6，可以让模型提前预测6步，不过它需要在观测到数据的6步之后才能利用

里面有大量不必要的特征。

# 5、一些疑惑以及技术选型

1、时间序列的（非）线性模型中将时间序列来回移动获取的特征的意义还不太懂，感觉是几阶差分。但是从拟合结果来看感觉还可以。还需要更细致的研究一下。但是无论是线性回归还是随机森林得到的pred，时间序列开始的部分匹配率都很差，得研究一下为什么。
2、SARIMA模型的ACF 图和 PACF 图看不懂，并且在工作中，构建模型的原则是快、好、省。 这也就意味着有些模型并不适合用于生产环境。因为它们需要过长的数据准备时间，或者需要频繁地重新训练新数据，或者很难调整参数（SARIMA 模型就包含了着三个缺点）。
3、双指数平滑、三指数平滑调参比较复杂，不过也有一定实现意义
4、moving average简单，不过容易把一些周期性的峰值判断成异常数据，这是因为它没有在我们的数据中捕捉到周期中出现的季节性变化
5、划分数据集和训练集的代码还得再看看。
6、当数据库中数据被填满的时候（也就是7天的数据量：7 x 24 x 6 = 1,008），每插入一条新数据，最旧的一条数据就被会抛弃掉。所以序列整体就是会平移。这样的话使用旧模型就会有问题了，可以通过编程手段，拟合出一周的模型：周一0点->周日24点的模型。然后每次获取db数据的时候进行剪裁，换算出当前数据在本周的位置。
7、考虑对哪些维度的数据进行模型拟合。

# 6、文章参考

【Python】时间序列分析完整过程：https://blog.csdn.net/jh1137921986/article/details/90257764
时间序列与时间序列分析：https://www.cnblogs.com/tianqizhi/p/9277376.html
如何根据自相关（ACF）图和偏自相关（PACF）图选择ARIMA模型的p、q值：https://blog.csdn.net/weixin_41013322/article/details/108801516

# 7、代码附录

前两个py文件放在/var/www/html/NewTest/InternShipProject/my_pylib目录下作为库文件，所以细节代码就在这里面了。

## 传统模型法

my_time_series_algorithm.py

import numpy as np                               # 向量和矩阵运算
import pandas as pd                              # 表格与数据处理
import matplotlib.pyplot as plt                  # 绘图
import seaborn as sns                            # 更多绘图功能
sns.set()

from dateutil.relativedelta import relativedelta # 日期数据处理
from scipy.optimize import minimize              # 优化函数

import statsmodels.formula.api as smf            # 数理统计
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs

from itertools import product                    # 一些有用的函数
from tqdm import tqdm_notebook

import warnings                                  # 勿扰模式
warnings.filterwarnings('ignore')

from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error
from sklearn.model_selection import TimeSeriesSplit # you have everything done for you

def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

# 加权平均
def weighted_average(series, weights):
"""
Calculate weighter average on series
"""
result = 0.0
weights.reverse()
for n in range(len(weights)):
result += series.iloc[-n-1] * weights[n]
return float(result)

def move_average(series, window):
return series.rolling(window=window).mean()

# 滑动平均
def plotMovingAverage(series, window, plot_intervals=False, scale=2.96, plot_anomalies=False, pic_name="plotMovingAverage"):

"""
series - dataframe with timeseries
window - rolling window size
plot_intervals - show confidence intervals
plot_anomalies - show anomalies
"""
rolling_mean = move_average(series, window)

plt.figure(figsize=(15,5))
plt.title("Moving average\n window size = {}".format(window))
plt.plot(rolling_mean, "g", label="Rolling mean trend")
plt.savefig('/var/www/html/NewTest/pics/%s.png' % "rolling_mean")

# Plot confidence intervals for smoothed values
if plot_intervals:
mae = mean_absolute_error(series[window:], rolling_mean[window:])
deviation = np.std(series[window:] - rolling_mean[window:])
lower_bond = rolling_mean - (mae + scale * deviation)
upper_bond = rolling_mean + (mae + scale * deviation)
plt.plot(upper_bond, "r--", label="Upper Bond / Lower Bond")
plt.plot(lower_bond, "r--")

# Having the intervals, find abnormal values
if plot_anomalies:
anomalies = series[window:].copy()
print(anomalies)
for i in series[window:].index:
if (series[window:][i] >= lower_bond[i]) and (series[window:][i] <= upper_bond[i]):
anomalies[i] = None
print(anomalies)
plt.plot(anomalies, "ro", markersize=10)

plt.plot(series[window:], label="Actual values")
plt.legend(loc="upper left")
plt.savefig('/var/www/html/NewTest/pics/%s.png' % pic_name)
plt.grid(True)

# 单指数平滑
def exponential_smoothing(series, alpha):
"""
series - dataset with timestamps
alpha - float [0.0, 1.0], smoothing parameter
"""
result = [series[0]] # first value is same as series
for n in range(1, len(series)):
result.append(alpha * series[n] + (1 - alpha) * result[n-1])
return result

# 绘制单指数平滑曲线
def plotExponentialSmoothing(series, alphas, pic_name="plotExponentialSmoothing"):
"""
Plots exponential smoothing with different alphas

series - dataset with timestamps
alphas - list of floats, smoothing parameters

"""
with plt.style.context('seaborn-white'):
plt.figure(figsize=(15, 7))
for alpha in alphas:
plt.plot(exponential_smoothing(series, alpha), label="Alpha {}".format(alpha))
plt.plot(series.values, "c", label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Exponential Smoothing")
plt.savefig('/var/www/html/NewTest/pics/%s.png' % pic_name)
plt.grid(True);

# 双指数平滑
def double_exponential_smoothing(series, alpha, beta):
"""
series - dataset with timeseries
alpha - float [0.0, 1.0], smoothing parameter for level
beta - float [0.0, 1.0], smoothing parameter for trend
"""
# first value is same as series
result = [series[0]]
for n in range(1, len(series)+1):
if n == 1:
level, trend = series[0], series[1] - series[0]
if n >= len(series): # forecasting
value = result[-1]
else:
value = series[n]
last_level, level = level, alpha*value + (1-alpha)*(level+trend)
trend = beta*(level-last_level) + (1-beta)*trend
result.append(level+trend)
return result

# 绘制双指数平滑曲线
def plotDoubleExponentialSmoothing(series, alphas, betas, pic_name="plotExponentialSmoothing"):
"""
Plots double exponential smoothing with different alphas and betas

series - dataset with timestamps
alphas - list of floats, smoothing parameters for level
betas - list of floats, smoothing parameters for trend
"""

with plt.style.context('seaborn-white'):
plt.figure(figsize=(20, 8))
for alpha in alphas:
for beta in betas:
plt.plot(double_exponential_smoothing(series, alpha, beta), label="Alpha {}, beta {}".format(alpha, beta))
plt.plot(series.values, label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Double Exponential Smoothing")
plt.savefig('/var/www/html/NewTest/pics/%s.png' % pic_name)
plt.grid(True)

# Holt-Winters model
class HoltWinters:
"""
Holt-Winters model with the anomalies detection using Brutlag method

# series - initial time series
# slen - length of a season
# alpha, beta, gamma - Holt-Winters model coefficients
# n_preds - predictions horizon
# scaling_factor - sets the width of the confidence interval by Brutlag (usually takes values from 2 to 3)

"""

def __init__(self, series, slen, alpha, beta, gamma, n_preds, scaling_factor=1.96):
self.series = series
self.slen = slen
self.alpha = alpha
self.beta = beta
self.gamma = gamma
self.n_preds = n_preds
self.scaling_factor = scaling_factor

def initial_trend(self):
sum = 0.0
for i in range(self.slen):
sum += float(self.series[i+self.slen] - self.series[i]) / self.slen
return sum / self.slen

def initial_seasonal_components(self):
seasonals = {}
season_averages = []
n_seasons = int(len(self.series)/self.slen)
# let's calculate season averages
for j in range(n_seasons):
season_averages.append(sum(self.series[self.slen*j:self.slen*j+self.slen])/float(self.slen))
# let's calculate initial values
for i in range(self.slen):
sum_of_vals_over_avg = 0.0
for j in range(n_seasons):
sum_of_vals_over_avg += self.series[self.slen*j+i]-season_averages[j]
seasonals[i] = sum_of_vals_over_avg/n_seasons
return seasonals

def triple_exponential_smoothing(self):
self.result = []
self.Smooth = []
self.Season = []
self.Trend = []
self.PredictedDeviation = []
self.UpperBond = []
self.LowerBond = []

seasonals = self.initial_seasonal_components()

for i in range(len(self.series)+self.n_preds):
if i == 0: # components initialization
smooth = self.series[0]
trend = self.initial_trend()
self.result.append(self.series[0])
self.Smooth.append(smooth)
self.Trend.append(trend)
self.Season.append(seasonals[i%self.slen])

self.PredictedDeviation.append(0)

self.UpperBond.append(self.result[0] +
self.scaling_factor *
self.PredictedDeviation[0])

self.LowerBond.append(self.result[0] -
self.scaling_factor *
self.PredictedDeviation[0])
continue

if i >= len(self.series): # predicting
m = i - len(self.series) + 1
self.result.append((smooth + m*trend) + seasonals[i%self.slen])

# when predicting we increase uncertainty on each step
self.PredictedDeviation.append(self.PredictedDeviation[-1]*1.01)

else:
val = self.series[i]
last_smooth, smooth = smooth, self.alpha*(val-seasonals[i%self.slen]) + (1-self.alpha)*(smooth+trend)
trend = self.beta * (smooth-last_smooth) + (1-self.beta)*trend
seasonals[i%self.slen] = self.gamma*(val-smooth) + (1-self.gamma)*seasonals[i%self.slen]
self.result.append(smooth+trend+seasonals[i%self.slen])

# Deviation is calculated according to Brutlag algorithm.
self.PredictedDeviation.append(self.gamma * np.abs(self.series[i] - self.result[i])
+ (1-self.gamma)*self.PredictedDeviation[-1])

self.UpperBond.append(self.result[-1] +
self.scaling_factor *
self.PredictedDeviation[-1])

self.LowerBond.append(self.result[-1] -
self.scaling_factor *
self.PredictedDeviation[-1])

self.Smooth.append(smooth)
self.Trend.append(trend)
self.Season.append(seasonals[i%self.slen])

# 用于参数搜索
def timeseriesCVscore(params, series, loss_function=mean_squared_error, slen=144):
"""
Returns error on CV

params - vector of parameters for optimization
series - dataset with timeseries
slen - season length for Holt-Winters model
"""
# errors array
errors = []

values = series.values
alpha, beta, gamma = params

# set the number of folds for cross-validation
tscv = TimeSeriesSplit(n_splits=3)

# iterating over folds, train model on each, forecast and calculate error
for train, test in tscv.split(values):

model = HoltWinters(series=values[train], slen=slen,
alpha=alpha, beta=beta, gamma=gamma, n_preds=len(test))
model.triple_exponential_smoothing()

predictions = model.result[-len(test):]
actual = values[test]
error = loss_function(predictions, actual)
errors.append(error)

return np.mean(np.array(errors))

# 绘制三指数平滑曲线
def plotHoltWinters(series, model, plot_intervals=False, plot_anomalies=False):
"""
series - dataset with timeseries
plot_intervals - show confidence intervals
plot_anomalies - show anomalies
"""

plt.figure(figsize=(20, 10))
plt.plot(model.result, label = "Model")
plt.plot(series.values, label = "Actual")
error = mean_absolute_percentage_error(series.values, model.result[:len(series)])
plt.title("Mean Absolute Percentage Error: {0:.2f}%".format(error))

if plot_anomalies:
anomalies = np.array([np.NaN]*len(series))
anomalies[series.values<model.LowerBond[:len(series)]] = \
series.values[series.values<model.LowerBond[:len(series)]]
anomalies[series.values>model.UpperBond[:len(series)]] = \
series.values[series.values>model.UpperBond[:len(series)]]
plt.plot(anomalies, "o", markersize=10, label = "Anomalies")

if plot_intervals:
plt.plot(model.UpperBond, "r--", alpha=0.5, label = "Up/Low confidence")
plt.plot(model.LowerBond, "r--", alpha=0.5)
plt.fill_between(x=range(0,len(model.result)), y1=model.UpperBond,
y2=model.LowerBond, alpha=0.2, color = "grey")

plt.vlines(len(series), ymin=min(model.LowerBond), ymax=max(model.UpperBond), linestyles='dashed')
plt.axvspan(len(series)-20, len(model.result), alpha=0.3, color='lightgrey')
plt.grid(True)
plt.axis('tight')
plt.legend(loc="best", fontsize=13);
plt.savefig('/var/www/html/NewTest/pics/%s.png' % "plotHoltWinters")

# 绘制时间序列图、ACF 图和 PACF 图
def tsplot(y, lags=None, figsize=(12, 7), style='bmh'):
"""
Plot time series, its ACF and PACF, calculate Dickey–Fuller test

y - timeseries
lags - how many lags to include in ACF, PACF calculation
"""
if not isinstance(y, pd.Series):
y = pd.Series(y)

with plt.style.context(style):
fig = plt.figure(figsize=figsize)
layout = (2, 2)
ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
acf_ax = plt.subplot2grid(layout, (1, 0))
pacf_ax = plt.subplot2grid(layout, (1, 1))

y.plot(ax=ts_ax)
ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p={0:.5f}'.format(p_value))
smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
plt.tight_layout()
plt.savefig('/var/www/html/NewTest/pics/%s.png' % "tsplot")



## 机器学习的方法（线性回归、XGB）

my_machine_learning_algorithm.py

import numpy as np                               # 向量和矩阵运算
import pandas as pd                              # 表格与数据处理
import matplotlib.pyplot as plt                  # 绘图
import seaborn as sns                            # 更多绘图功能
sns.set()

from dateutil.relativedelta import relativedelta # 日期数据处理
from scipy.optimize import minimize              # 优化函数

import statsmodels.formula.api as smf            # 数理统计
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs

from itertools import product                    # 一些有用的函数
from tqdm import tqdm_notebook

import warnings                                  # 勿扰模式
warnings.filterwarnings('ignore')

from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error
from sklearn.model_selection import TimeSeriesSplit # you have everything done for you

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

from xgboost import XGBRegressor

def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

def timeseries_train_test_split(X, y, test_size):
"""
Perform train-test split with respect to time series structure
"""

# get the index after which test set starts
test_index = int(len(X)*(1-test_size))

X_train = X.iloc[:test_index]
y_train = y.iloc[:test_index]
X_test = X.iloc[test_index:]
y_test = y.iloc[test_index:]

return X_train, X_test, y_train, y_test

def plotModelResults(model, X_train, X_test, y_test, y_train, tscv, plot_intervals=False, plot_anomalies=False, scale=1.96, pic_name="plotModelResults"):
"""
Plots modelled vs fact values, prediction intervals and anomalies

"""

prediction = model.predict(X_test)

plt.figure(figsize=(15, 7))
plt.plot(prediction, "g", label="prediction", linewidth=2.0)
plt.plot(y_test.values, label="actual", linewidth=2.0)

if plot_intervals:
cv = cross_val_score(model, X_train, y_train,
cv=tscv,
scoring="neg_mean_squared_error")
#mae = cv.mean() * (-1)
deviation = np.sqrt(cv.std())

lower = prediction - (scale * deviation)
upper = prediction + (scale * deviation)

plt.plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)
plt.plot(upper, "r--", alpha=0.5)

if plot_anomalies:
anomalies = np.array([np.NaN]*len(y_test))
anomalies[y_test<lower] = y_test[y_test<lower]
anomalies[y_test>upper] = y_test[y_test>upper]
plt.plot(anomalies, "o", markersize=10, label = "Anomalies")
# print(prediction.flatten())
# print(y_test)
error = mean_absolute_percentage_error(prediction.flatten(), y_test)
plt.title("Mean absolute percentage error {0:.2f}%".format(error))
plt.legend(loc="best")
plt.tight_layout()
plt.grid(True);
plt.savefig('/var/www/html/NewTest/pics/%s.png' % pic_name)

def plotCoefficients(model, x_train, pic_name="plotCoefficients"):
"""
Plots sorted coefficient values of the model
"""

coefs = pd.DataFrame(model.coef_.flatten(), x_train.columns)
coefs.columns = ["coef"]
coefs["abs"] = coefs.coef.apply(np.abs)
coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)

plt.figure(figsize=(15, 7))
coefs.coef.plot(kind='bar')
plt.grid(True, axis='y')
plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed');
plt.savefig('/var/www/html/NewTest/pics/%s.png' % pic_name)



## demo

### sarima_demo

import pymysql                                   # 数据库
# 导入mylib
import sys
sys.path.append('/var/www/html/NewTest/InternShipProject/my_pylib')
from my_time_series_algorithm import *

# read_sql(sql, con, index_col=None, coerce_float=True, params=None, parse_dates=None, columns=None, chunksize=None)
conn = pymysql.connect(
host = 'localhost',
user = 'root',
passwd = '123456',
db = 'prod_179',
port=3306,
charset = 'utf8'
)
df = pd.read_sql('select * from 北京市',conn)

tsplot(df.rtt, lags=60)



### exp_smooth_demo

import pymysql                                   # 数据库
# 导入mylib
import sys
sys.path.append('/var/www/html/NewTest/InternShipProject/my_pylib')
from my_time_series_algorithm import *

# read_sql(sql, con, index_col=None, coerce_float=True, params=None, parse_dates=None, columns=None, chunksize=None)
conn = pymysql.connect(
host = 'localhost',
user = 'root',
passwd = '123456',
db = 'prod_179',
port=3306,
charset = 'utf8'
)
df = pd.read_sql('select * from all_nums where id between 864 and 1008',conn)

# 双指数
# 绘制双指数平滑曲线
def plotDoubleExponentialSmoothingXXX(series, alpha, beta, pic_name="plotExponentialSmoothing"):
"""
Plots double exponential smoothing with different alphas and betas

series - dataset with timestamps
alphas - list of floats, smoothing parameters for level
betas - list of floats, smoothing parameters for trend
"""

with plt.style.context('seaborn-white'):
plt.figure(figsize=(20, 8))
res_line = double_exponential_smoothing(series, alpha, beta)
print(res_line[-1])
plt.plot(res_line, label="Alpha {}, beta {}".format(alpha, beta))
plt.plot(series.values, label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Double Exponential Smoothing")
plt.savefig('/var/www/html/NewTest/pics/%s.png' % pic_name)
plt.grid(True)

# plotDoubleExponentialSmoothingXXX(df.user_nums,0.9, 0.49, pic_name="双指数平滑效果")
plotDoubleExponentialSmoothing(df.user_nums, [0.78], [0.125, 0.1], pic_name="双指数平滑效果")
conn.close()
# # 单指数
# plotExponentialSmoothing(df.rtt, [0.3, 0.05])



### linear_reg_demo

import pymysql                                   # 数据库
# 导入mylib
import sys
sys.path.append('/var/www/html/NewTest/InternShipProject/my_pylib')
from my_machine_learning_algorithm import *

# read_sql(sql, con, index_col=None, coerce_float=True, params=None, parse_dates=None, columns=None, chunksize=None)
conn = pymysql.connect(
host = 'localhost',
user = 'root',
passwd = '123456',
db = 'prod_179',
port=3306,
charset = 'utf8'
)
df = pd.read_sql('select * from 北京市',conn)

# Creating a copy of the initial datagrame to make various transformations
data = pd.DataFrame(df.rtt.copy())
data.columns = ["y"]

# Adding the lag of the target variable from 6 steps back up to 24
for i in range(6, 25):
data["lag_{}".format(i)] = data.y.shift(i)

# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)

y = data.dropna().y
x = data.dropna().drop(['y'], axis=1)

# reserve 30% of data for testing
x_train, x_test, y_train, y_test = timeseries_train_test_split(x, y, test_size=0.3)

# # machine learning in two lines
lr = LinearRegression()
lr.fit(x_train, y_train.values.reshape(-1,1))
# # print(x_train.shape)
# # print(y_train.shape)
# # (527, 19)
# # (527,)
plotModelResults(lr, x_train, x_test, y_test, y_train, tscv, plot_intervals=True)
plotCoefficients(lr, x_train)



### xgb_demo

import pymysql                                   # 数据库
# 导入mylib
import sys
sys.path.append('/var/www/html/NewTest/InternShipProject/my_pylib')
from my_machine_learning_algorithm import *

# read_sql(sql, con, index_col=None, coerce_float=True, params=None, parse_dates=None, columns=None, chunksize=None)
conn = pymysql.connect(
host = 'localhost',
user = 'root',
passwd = '123456',
db = 'prod_179',
port=3306,
charset = 'utf8'
)
df = pd.read_sql('select * from all_nums',conn)

# Creating a copy of the initial datagrame to make various transformations
data = pd.DataFrame(df.flow_up.copy())
data.columns = ["y"]

# Adding the lag of the target variable from 6 steps back up to 24
# for i in range(6, 25):
#     data["lag_{}".format(i)] = data.y.shift(i)

# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)

y = data.dropna().y
x = data.dropna().drop(['y'], axis=1)

# reserve 30% of data for testing
x_train, x_test, y_train, y_test = timeseries_train_test_split(x, y, test_size=0.3)

X_train_scaled = scaler.fit_transform(x_train)
X_test_scaled = scaler.transform(x_test)

xgb = XGBRegressor()
xgb.fit(X_train_scaled, y_train)

plotModelResults(xgb,
X_train=X_train_scaled,
X_test=X_test_scaled,
y_test=y_test,
y_train=y_train,
tscv=tscv,
plot_intervals=True, plot_anomalies=True)


展开全文
• 第五章、时间序列计量经济学模型5.1时间序列模型的序列相关性5.1.1序列相关概念对一时间序列模型 关于随机误差项有一基本假定是它们之间不相关，即 当这一假定不再满足时，如果其他基本假定都满足，则 ，此时 如果...

## 5.1时间序列模型的序列相关性

### 5.1.1序列相关概念

对一时间序列模型

关于随机误差项有一基本假定是它们之间不相关，即

当这一假定不再满足时，如果其他基本假定都满足，则

，此时

如果相关的解释变量只相差一期

，称为
一阶序列相关一阶自相关。一般我们 设这一相关是线性的，即
，其中
自相关系数
是满足OLS假定的误差项。

如果是

，称为
s阶序列相关s阶自相关

### （一）经济现象固有惯性

• 当期投资与前几年的投资有关；
• 当期消费受原有消费水平制约；
• 企业当期产量与前几期产量密切相关。

因为被解释变量

自相关，引起随机误差项序列相关。

### （二）随机误差项本身自相关

一些影响因素，如战争、自然灾害，会持续若干时期。

例如在消费模型中，消费习惯处在误差项中，并且会持续一段时间，引起误差项的序列相关。

### （1）遗漏重要变量：

误差项中包含该重要变量，会随这一重要变量的变化而变化，进而产生序列相关。

这是遗漏变量的又一后果

### （2）函数形式不当：

如这是模型中带有平方项，误设模型中没有，那么平方项就在误差中，误差就会随平方项变动，进而产生序列相关。

又如真实模型中带有滞后项，误设模型中没有。类似地也会造成序列相关。

### （四）数据处理不当

如按月数据转化成按季度数据，需要进行平滑处理：一二三月、二三四月、三四五月取平均...这样其中相邻的两个数据包含相同月份，产生序列相关。

经验：因为误差项中包含的因素往往在时间上连续，所以时间序列数据往往存在序列相关

### （一）参数估计量非有效

原因：有效性证明用到

大样本具有一致性，但不具有渐近有效性。

### （二）变量显著性检验失效

原因：t统计量正确

OLS估计量的方差正确估计
误差项同方差+相互独立

在一阶序列相关

下，OLS估计量方差为

如果按基本假定做，只有第一项，即出现偏误。

### （三）预测失效

原因：预测的区间中包含OLS估计量的方差。

### 5.1.4检验序列相关

思路：通过分析残差序列的相关性来判断随机误差项是否序列相关。

做出

的图像

典型负相关：蛛网模型。

### （二）回归检验法

思路：残差关于其滞后项做各种回归，看是否出现显著。

优点：能确定几阶序列相关；适用于任何类型的序列相关检验。

### （三）D.W.检验

假定条件：

• 解释变量非随机（不能内生）
• 误差项一阶自相关（D.W.检验只能检验是否存在一阶自相关）
• Y的滞后不作解释变量
• 含有截距项

检验模型

，假设为
，构造检验统计量为

比较大时，
，结合
得到

检验步骤：

• 计算DW检验统计量的值（不用手算）
• 查表得临界值
注意k是包含常数项的解释变量个数！
• 比较判断

根据DW分子

可以这样
记忆
• 正相关是正数减正数平方，DW值偏小；
• 负相关是正数减负数，DW值偏大。

落入不能确定区域：改变样本或改变模型函数形式。

一般经验：不存在一阶自相关，也就不存在高阶自相关。

### （四）拉格朗日乘数检验--GB检验

优点：能检验高阶自相关；能处理滞后应变量做解释变量的模型。

怀疑有

阶自相关

要检验的模型为

检验的假设为

，取原模型的残差

做辅助回归

检验统计量为

其中

辅助回归的样本容量（因为滞后p期，所以少了p个数据），
辅助回归的可决系数

比较判断：如果LM统计量大于临界值，那么拒绝原假设，表明可能存在直到p阶自相关。

说明：实际可以从1阶开始逐步增加阶数，并通过辅助回归中变量的显著性来判断自相关阶数

### （一）广义最小二乘GLS

如果序列相关和异方差同时存在，那么协方差矩阵为

其中

是对称正定矩阵，所以
可逆。

原模型

，两边左乘
的逆

变换后的模型

，计算该模型误差项的协方差矩阵

此时已经消除异方差和序列相关！利用OLS估计量的形式可得GLS估计量为

形式上就是在OLS估计量中间插了两个

的逆

需要对模型进行设定才能得到

，如假设一阶自相关

此时可以证明

于是协方差矩阵为

后面用代数方法将

分解成
并得到
的逆即可，最终结果为

也就是把误差项

变换成了满足基本假定的

### （二）广义差分法

思想：通过配凑把序列相关的随机误差项

变为

可以用于多元模型、多阶相关。以一元模型，一阶相关为例：

对于一阶相关的一元模型

原模型滞后一期

得到

如何得到

？利用

注意：广义差分是广义最小二乘损失部分的观测值！

如果我们把损失的第一次观测值设为

，那么广义差分就和广义最小二乘完全等同了，这一变化称为
普莱斯-温斯特变换

### （三）随机误差项相关系数的估计

问题：DW只能告诉我们一个

，对于多元如何估计其他相关系数？

### 科克伦-奥科特( Cochrane-Orcutt)迭代法

思路：用残差代替误差做回归，用回归系数估计相关系数。

步骤：以一元模型一阶相关为例

回归得到

计算残差

辅助回归

得到

利用估计值建立差分模型

将差分模型的OLS估计量

带回原模型（计算残差，辅助回归），重复迭代

停止条件：相邻两个

的估计值相差很小或模型已不存在序列相关。

消除一阶序列相关的Eviews命令：LS Y C X AR(1)

注意：如果相关系数是被估计出来的，那么不是广义最小二乘，而是可行的广义最小二乘FGLS。FGLS估计量不再是无偏的，但具有一致性，并且用科克伦-奥科特法得到的FGLS估计量还有渐近有效性。

### （四）序列相关稳健标准误法

思路：仍用OLS，但改变OLS估计量的方差。（消除后果）

如对于一元一阶相关，斜率项用正确的方差估计

### 5.1.6虚假序列相关

概念：由于模型设定偏误引起的随机误差项序列相关。

解决：在模型上处理，而不是用前面序列相关的修正方法。

展开全文
• 看过掌柜前几篇关于时序文章： 利用ARIMA模型对时间序列进行分析的经典案例（详细代码） “利用ARIMA模型对时间序列进行分析的经典...Python是如何使用SARIMA模型进行时序分析的？（SARIMA模型的步骤有哪些？） SA
• 背景：SARIMA，简单说就是AR+MA+差分+季节性因素+趋势。所以参数在statsmodels.tsa.statespace.sarimax.SARIMAX里边，用3个指标涵盖核心参数，order(p,d,q)、seasonal_order(P,D,Q,s)和trend. Seasonal ...
• 问题1：计量经济学是分析啥的？包含些什么内容？计量经济学的主要用途或目的主要有两个方面：1.理论检验。这是计量经济学用途最为主要的和可靠的方面。这也是计量经济学本身的一个主要内容。2.预测应用。...
• 指数平滑模型是基于对数据趋势和季节性的描述，而ARIMA模型则是为了描述数据的自相关性。 在讨论ARIMA模型之前，我们先来讨论平稳性的概念和时间序列的差分技术。 平稳性 平稳时间序列数据的性质不依赖于时间，...
• 使用ARIMA模型，您可以使用序列过去的值预测时间序列。在本文中，我们从头开始构建了一个最佳ARIMA模型，并将其扩展到Seasonal ARIMA（SARIMA）和SARIMAX模型
• 基于SARIMA模型的分析及预测.docx
• 基于SARIMA模型的分析及预测.pdf
• 基于SARIMA模型分析与外汇管理知识分析成本.doc
• 基于SARIMA模型分析与外汇管理知识分析成本.docx
• SARIMA模型在新疆布鲁氏菌病发病预测中的应用.pdf
• 日本地震间接经济损失的SARIMA模型估计，梁淇俊，王斌会，过往关于地震间接经济损失的研究大多停留在社会经济关联型损失上，导致对间接经济损失的评估不够全面。本文通过典型相关分析以及
• 基于SARIMA模型分析与外汇管理知识分析成本(共39页).doc
• SARIMA(p,d,q)(P,D,Q,s) 季节性自回归移动平均模型，结构参数有七个 AR(p) 自回归模型，即用自己回归自己。基本假设是，当前序列值取决于序列的历史值。p 表示用多少个历史值来回归出预测值。 要确定初始 p，需要...
• 昨天做完ARIMA的模型，由于对于其中seasonal部分的不理解，外加我的用电负荷模型应该满足一个季节特征，因此得到建议去了解一下SARIMA模型 先看相关教程： 《R语言时间序列分析》 memo： （1）将数据从文本...
• SARIMA时间序列模型预测城市房价数据 数据清洗 文件中含有大量城市的房价数据，考虑到此次为学习性质的练习，为了节省数据处理的繁琐步骤。我截取了北京的2010-2021房价数据作为样例，并将价格的数据格式改为数值，...
• 中国快递包裹总量：基于SARIMA模型的预测 shangfr 2015年10月28日 目前，我国快递业务总量已成为全球第一。国家邮政局发布的数据显示，2015年4月底，快递业务量完成15亿件，同比增长50.9%。目前，快递业务量增速已...
• 去年用过EViews来做SARIMA模型，今天又要用到，结果操作全忘光了，好不容易摸索出来，故在此小记下。 1. File->New->Workfile，出来窗口，Frequency，Start date和End date都是对等会要处理的时间序列说的，...
• 然而在实际还会遇到一些其他的模型ARIMA、SARIMA。而模型参数也会也会变的复杂，从pq到pdq到pdqPDQS。本文主要通过R语言生成符合三种模型的随机数，进行建模预测。 在写模型之前先给大家介绍一个时间序列总结的不错...
• SARIMA 本系列前面我们介绍过MA(q)Moving Average/移动平均模型、AR/ Autoregressive/自回归、ARMA(p,q)Autoregressive Moving Average/自回归移动平均、ARIMA(p,q,d),Autoregressive Integrated Moving Average/差...
• 但是，使用具有365天周期的输入数据来计算SARIMA模型需要大量资源，并且在大多数计算机中可能会失败。 What is FB Prophet? 什么是FB先知？ Prophet is an open Python/R library for time series forecasting ...
• ARMA、ARIMA、AR、MA均是时间序列的重要方法。此例程中包含了以上所有的实现过程，java实现的，且含有main函数供自行调试，已亲测可用！
• 开只选取了120即10的数据来进行分析，但是到最后发现模型有很多的波动之后，去问老师，老师说这是数据太少导致波动太大造成的，所以建议我们再多训练一些数据。就之后进行模型的定阶而言，至少需要三四百的数据。...

...