• 季节性ARIMA：时间序列预测
万次阅读
2019-02-17 15:22:36

SARIMAX (seasonal autoregressive integrated moving average with exogenous regressor)是一种常见的时间序列预测方法，可以分为趋势部分和周期性部分；每个部分又可以分为自回归、差分和平滑部分。

null-hypothesis: 时间序列趋势稳定。significance: 0.05. 选择KPSS和不是Dickey Fuller由于KPSS的null-hypothesis是趋势稳定，所以接受条件相对于Dickey Fuller更宽松，差分阶数更少。
y t = β ′ D t + μ t + u t y_t = \beta^{\prime} D_t + \mu_t + u_t
μ t = μ t − 1 + ϵ t , \mu_t = \mu_{t-1} + \epsilon_t,
ϵ t ∼ W N ( 0 , σ ϵ 2 ) \epsilon_t \sim WN(0, \sigma^2_\epsilon)

其中D为确定时间序列趋势/常数。u为随机漫步。epsilon为残差。

LM statistic:
K P S S = ( T − 2 ∑ t = 1 T S ^ t 2 ) / λ ^ 2 KPSS = (T^-2 \sum_{t=1}^{T} \hat S_t^2) / \hat\lambda^2
其中 S ^ t = ∑ j = 1 t u ^ j \hat S_t = \sum_{j=1}^{t} \hat u_j
u ^ \hat{u} 是yt 拟合Dt的残差，$\hat \lambda^2$ 是var(ut)的预估值。问题归结为拉格朗日乘数法证明 σ ϵ 2 = 0 \sigma^2_\epsilon = 0

季节稳定性检测Canova-Hansen方法
null-hypothesis: 时间序列季节稳定性。significance: 0.05

测试频率：给定待测最大频率m以下所有 2pi/m整数倍。

时间序列yi的拟合：
y i = μ + x i ′ β + S i + e i y_i = \mu + x^{\prime}_i \beta + S_i + e_i
S i = ∑ j = 1 m / 2 f j i ′ γ j S_i = \sum_{j=1}^{m/2} f^{\prime}_{ji} \gamma_j
其中
f j i ′ = ( c o s ( ( j / q ) π i ) , s i n ( ( j / q ) π i ) ) f^{\prime}_{ji} = (cos((j/q) \pi i), sin((j/q) \pi i))

LM statistic:
L M = 1 n 2 t r a c e ( ( Ω ^ f ) − 1 ∑ i = 1 n F ^ i F ^ i ′ ) LM = \frac{1}{n^2} trace( (\hat \Omega ^f)^-1 \sum_{i=1}^{n} \hat F_i \hat F^{\prime}_i)

其中
F ^ i = ∑ t = 1 i f t e ^ t \hat F_i = \sum_{t=1}^{i}f_t \hat e_t
Ω ^ f = ∑ k = − m m w ( k m ) 1 n ∑ i f i + k e ^ i + k f i ′ e ^ i \hat \Omega^f = \sum_{k=-m}^{m} w (\frac{k}{m}) \frac{1}{n} \sum_{i} f_{i+k} \hat e_{i+k} f^{\prime}_i \hat e_i

用根据待检测频率m计算f’ji向量：用sample data示例代码如下：

import numpy as np
import pandas as pd

def seasonalDummy(tsArray, frequency):
n = len(tsArray)
m = frequency
#if m == 1: tsArray.reshape([n, m])
tt = np.arange(1, n + 1, 1)
mat = np.zeros([n, 2 * m], dtype = float)
for i in np.arange(0, m):
mat[:, 2 * i] = np.cos(2.0 * np.pi * (i + 1) * tt / m)
mat[:, 2 * i + 1] = np.sin(2.0 * np.pi * (i + 1) *tt / m)
return mat[:, 0:(m-1)]

tsArray = np.array([  4.00195672,   4.99944175,   5.99861146,   7.00000213,
8.00124207,   9.00059539,  10.00029542,  10.99871969,
11.99728933,  12.99965231,  14.00148869,  15.0006378 ,
16.00117112,  17.00159081,  17.99848509,  18.99957668,
20.00066721,  20.99932292,  21.99992471,  23.00099164,
24.00127222,  25.00014385,  26.00014191,  27.00037435,
27.9985619 ,  28.99949718,  29.99844772,  30.99911627,
31.99965086,  33.00211019,  34.00240288,  34.99889972,
36.00240406,  37.0002379 ,  37.99958145,  38.99825111,
39.99932529,  40.9998374 ,  42.00034236,  43.00206289,
43.9994545 ,  45.00141283,  46.00041818,  47.00132581,
48.00216031,  48.99812424,  50.00060522,  51.00049892,
51.99817633,  52.9997362 ])
frequency = 6


Canova-Hansen用sample data示例代码如下：

from numpy.linalg import lstsq as lsq

n = len(tsArray)
frec = np.ones(int((frequency + 1) / 2), dtype = int)
ltrunc = int(np.round(frequency * np.power(n / 100.0, 0.25)))


y i = μ + x i ′ β + S i + e i y_i = \mu + x^{\prime}_i \beta + S_i + e_i
S i = ∑ j = 1 m / 2 f j i ′ γ j S_i = \sum_{j=1}^{m/2} f^{\prime}_{ji} \gamma_j

其中 f j i ′ = ( c o s ( ( j / q ) π i ) , s i n ( ( j / q ) π i ) ) f^{\prime}_{ji} = (cos((j/q) \pi i), sin((j/q) \pi i))

# create dummy column f'ji
r1 = seasonalDummy(tsArray, frequency)
#create intercept column for regression
r1wInterceptCol = np.column_stack([np.ones(r1.shape[0], dtype = float), r1])

# residual ei:
result = lsq(a = r1wInterceptCol, b = tsArray)
residual = tsArray - np.matmul(r1wInterceptCol, result[0])


long-run covariance matrix： Ω = l i m n → ∞ 1 n E ( F n F n ′ ) \Omega = lim_{n \to \infty}\frac{1}{n}E(F_n F_n^{\prime})
在ei可能有serial correlation的情况下可以用estimate
Ω ^ = ∑ k = − m m w ( k m ) 1 n ∑ i F i + k F i ′ \hat{\Omega} = \sum_{k=-m}^{m}w(\frac{k}{m})\frac{1}{n}\sum_{i}F_{i+k}F_i^{\prime}

fhat = np.zeros([n, frequency - 1], dtype = float)
fhataux = np.zeros([n, frequency - 1], dtype = float)

for i in np.arange(0, frequency - 1):
fhataux[:, i] = r1[:, i] * residual

for i in np.arange(0, n):
for j in np.arange(0, frequency - 1):
mySum = sum(fhataux[0:(i + 1), j])
fhat[i, j] = mySum


w ( ⋅ ) w(\cdot) 是kernel function，来自rob j. Hyndman的forecast包:

wnw = np.ones(ltrunc, dtype = float) - np.arange(1, ltrunc + 1, 1) / (ltrunc + 1)


Ω ^ f \hat{\Omega}^f 的计算：

Ne = fhataux.shape[0]
omnw = np.zeros([fhataux.shape[1], fhataux.shape[1]], dtype = float)
for k in np.arange(0, ltrunc):
omnw = omnw + np.matmul(fhataux.T[:, (k+1):Ne], fhataux[0:(Ne-(k+1)), :]) * float(wnw[k])
cross = np.matmul(fhataux.T, fhataux)
omnwplusTranspose = omnw + omnw.T
omfhat = (cross + omnwplusTranspose) / float(Ne)


Generalized Hannan’s model：通过设置矩阵A选择需要测试的频率的子集。在程序中用的A = eye:

sq = np.arange(0, frequency - 1, 2)
frecob = np.zeros(frequency - 1, dtype = int)
for i in np.arange(0, len(frec)):
if (frec[i] == 1) & (i + 1 == int(frequency / 2.0)):
frecob[sq[i]] = 1
if (frec[i] == 1) & (i + 1 < int(frequency / 2.0)):
frecob[sq[i]] = 1
frecob[sq[i] + 1] = 1

a = frecob.tolist().count(1)  # find nr of 1's in frecob
A = np.zeros([frequency - 1, a], dtype = float)
j = 0
for i in np.arange(0, frequency - 1):
if frecob[i] == 1:
A[i, j] = 1
j = j + 1


LM statistic的计算 （Nyblom(1989), Hansen(1990)）：当LM statistic值超过对应自由度的critical value时，拒绝 （null hypothesis = 没有单位根）：
L M s t a t i s t i c = 1 n 2 t r a c e ( ( A ′ Ω ^ f A ) − 1 A ′ ∑ i = 1 n F i ^ F i ^ ′ A ) LM statistic = \frac{1}{n^2} trace((A&#x27; \hat{\Omega}^f A)^{-1} A&#x27; \sum_{i=1}^{n}\hat{F_i} \hat{F_i}^{\prime}A)

其中 ∑ i = 1 n F i ^ \sum_{i=1}^{n}\hat{F_i} ∑ i = 1 n F i ^ ′ \sum_{i=1}^{n} \hat{F_i}^{\prime} 由前面步骤中的 f ^ \hat{f} 给出:

from numpy.linalg import svd

aTomfhat = np.matmul(A.T, omfhat)
tmp = np.matmul(aTomfhat, A)

machineDoubleEps = 2.220446e-16

problems = min(svd(tmp)[1]) < machineDoubleEps # svd[1] are the singular values
if problems:
stL = 0.0
else:
solved = np.linalg.solve(tmp, np.eye(tmp.shape[1], dtype = float))
step1 = np.matmul(solved, A.T)
step2 = np.matmul(step1, fhat.T)
step3 = np.matmul(step2, fhat)
step4 = np.matmul(step3, A)
stL = (1.0 / np.power(n, 2.0)) * sum(np.diag(step4))


在Canova-Hansen test接受alternative hypothesis的情况下对时间序列进行lag = 待测频率的差分（季节差分）。

季节性检测
在进行季节性时间序列稳定性检测之前，首先判断a.时间序列是否有季节性，和b.时间序列在什么频率上有季节性。结果会作为时间序列稳定性检测的参数输入。

季节性检测根据离散傅里叶变换和自相关函数的“与”关系得出结论（只有两个method都返回真值，才会判定时间序列有季节性）。

离散傅里叶变换吧时间序列从时域变为频域。变换后频域的新序列为：

X k = ∑ n = 0 N − 1 x n ⋅ [ c o s ( 2 π k n / N ) − i ⋅ s i n ( 2 π k n / N ) ] X_k = \sum_{n=0}^{N-1}x_n \cdot [cos(2 \pi kn/N) - i \cdot sin(2 \pi kn/N)]

在待检测频率上如果能量为最大值，则返回真值。

自相关函数检测最大lag = 待检频率各阶上的correlation coefficient。时间点t和s之间的自相关性R的定义为：

R ( s , t ) = E [ ( X t − μ t ) ( X s − μ s ) σ t σ s R(s, t) = \frac{E[(X_t - \mu_t)(X_s - \mu_s)}{\sigma_t \sigma_s}

如果在待检频率上的相关系数超过双边confidence interval在0.05的临界值

clim = norm.ppf((1 + ci) / 2) / np.sqrt(n)


则method返回真值。

DFT调用了numpy.fft.fft方法。ACF调用了statsmodels.tsa.stattools.acf方法。

测量函数：根据不同情况采用不同模型测量方法。算法使用了Rob J. Hyndman的 MASE (mean absolute scaled error)。与其他测量方法的优劣对比。

定义：
R-squared:
S S t o t = ∑ i ( y i − y ˉ ) 2 SS_{tot} = \sum_{i}(y_i - \bar{y})^2
S S r e s = ∑ i ( y i − y ^ i ) 2 SS_{res} = \sum_{i}(y_i - \hat{y}_i)^2
R 2 = 1 − S S r e s S S t o t R^2 = 1 - \frac{SS_{res}}{SS_{tot}}

RMSE:
R M S E = ∑ t = 1 n ( y ^ t − y t ) 2 n RMSE = \sqrt{\frac{\sum_{t=1}^{n} (\hat{y}_t - y_t)^2}{n}}

MAE:
M A E = ∑ t = 1 n ∣ y ^ t − y t ∣ n MAE = \frac{\sum_{t=1}^{n} |\hat{y}_t - y_t|}{n}

MAPE:
M A P E = 100 n ∑ t = 1 n ∣ y ^ t − y t y t ∣ MAPE = \frac{100}{n} \sum_{t=1}^{n} | \frac{\hat{y}_t - y_t}{y_t} |

sMAPE:
S M A P E = 100 n ∑ t = 1 n ∣ y ^ t − y t ∣ ( ∣ y ^ t ∣ + ∣ y t ∣ ) / 2 SMAPE = \frac{100}{n} \sum_{t=1}^{n} \frac{|\hat{y}_t - y_t|}{(|\hat{y}_t| + |y_t|)/2}

MASE without seasonality:
M A S E = ∑ t = 1 n ∣ y ^ t − y t ∣ n n − 1 ∑ t = 2 n ∣ y t − y t − 1 ∣ MASE = \frac{\sum_{t = 1}^{n} |\hat{y}_t - y_t|}{\frac{n}{n-1} \sum_{t=2}^{n} |y_t - y_{t-1}|}

MASE with seasonality:
M A S E = ∑ t = 1 n ∣ y ^ t − y t ∣ n n − m ∑ t = m + 1 n ∣ y t − y t − m ∣ MASE = \frac{\sum_{t = 1}^{n} |\hat{y}_t - y_t|}{\frac{n}{n-m} \sum_{t=m+1}^{n} |y_t - y_{t-m}|}

趋势平滑：在SARIMA模型中引入时间序列的趋势作为exogenous regressor(X)，有几种算法可以选择：

Lowess (locally weighted scatterplot smoothing): 基于KNN的非参数拟合方法。代码调用了 statsmodels.api.nonparametric.lowess

RANSAC
Random sample consensus，一种robust regression方法，可以探测异常值并使拟合对于异常值的敏感度降低。代码调用sklearn.linear_model.RANSACRegressor

Weighted moving average:

指数平滑递归表达：
W n = ( 1 − α ) ∗ W n − 1 + α ∗ y n W_n = (1-\alpha) * W_{n-1} + \alpha * y_n
W 0 = y 0 W_0 = y_0
α = 2 / ( s p a n + 1 ) \alpha = 2 / (span + 1)

调用了 pandas.ewma

关于脉冲响应

如果有另外的（多维）exogenous regressor Xi 影响预测模型，比如类似离散脉冲波形的机会点数据：假设各个脉冲regressor之间是独立的，并不受时间序列本身的影响：可以用多元线性回归首先发现Xi 和时间序列趋势yhat的关系。

考虑到历史上时间序列对于脉冲输入的响应不同，在算法中会测试三种不同脉冲响应与时间序列的相关性，并挑选相关性最强的脉冲响应作为Xi向量。这三种分别是a. 原始脉冲，b. 经指数平滑的脉冲，c. 经指数平滑并累加的脉冲（cumulative)。

示例：

%matplotlib inline
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings('ignore')

def myEwma(x, histPeriod = 6, fcstPeriod = 18):
from pandas import ewma
xfit = ewma(x, span = histPeriod)
xpred = np.zeros(fcstPeriod)
tmp = np.zeros(histPeriod + 1)
tmp[:histPeriod] = x[-histPeriod:].copy()
tmp[histPeriod] = xfit[-1]
for i in np.arange(0, fcstPeriod):
xpred[i] = ewma(tmp, span = histPeriod)[-1]
tmp = shiftLeft(tmp)
tmp[-1] = xpred[i]
return xfit, xpred

def shiftLeft(_ar, step = 1):
if step == 0:
return _ar
ar = _ar.copy()
ar[0:-step] = ar[step:]
ar[-step:] = 0
return ar

original = np.array([3000,0,0,0,0,0,0,0,0,0,1000,0,0,0,0,0,0,0,0,1000,0,0,0,0,0,0,0,0,0,0])
pd.DataFrame(original).plot(title = 'original')
movingAverage = myEwma(original)[0]
pd.DataFrame(movingAverage).plot(title = 'moving average')
cumulative = movingAverage.cumsum()
pd.DataFrame(cumulative).plot(title = 'cumulative')


检测相加性：时间序列拆分成趋势，季节性和残差的方式有相加和相乘两种。

y = Trend + Seasonality + Residual, 或者

y = Trend * Seasonality * Residual

如果是相乘的情况，残差的分布是不稳定的。所以如果时间序列没有通过相加性检测，会对时间序列做对数处理，变为相加：

y n e w = l o g ( y ) = l o g ( T r e n d ∗ S e a s o n a l i t y ∗ R e s i d u a l ) = l o g ( T r e n d ) + l o g ( S e a s o n a l i t y ) + l o g ( R e s i d u a l ) ynew = log(y) = log(Trend * Seasonality * Residual) = log(Trend) + log(Seasonality) + log(Residual)

自动搜索SARIMA参数空间 Auto ARIMA：搜索差分阶数，检测季节性的存在，并搜索能给出Akaike Information Criterion最小值的ARMA模型维度 ARMA(p, q, P, Q)

ARMA模型范式：

X t − α 1 X t − 1 − ⋯ − α p ′ X t − p ′ = ϵ t + θ 1 ϵ t − 1 + ⋯ + θ q ϵ t − q X_t - \alpha_1 X_{t-1} - \dots - \alpha_{p^{\prime}} X_{t-p^{\prime}} = \epsilon_t + \theta_1 \epsilon_{t-1} + \dots + \theta_q \epsilon_{t-q}

自动搜索参数空间返回的结果是p和q的个数。作为输入给SARIMAX核心算法

搜索参数空间的初始值来自Rob J. Hyndman的paper section3.2.

SARIMAX 带外部变量的季节性自回归差分平滑算法：预测是基于SARIMAX的state-space representation。 核心算法调用了statsmodels.api.tsa.statespace.SARIMAX.

总结

以上是SARIMAX各个功能模块的详细数学方法和编程实现。由于一些语言（包括python）的时间序列预测开源算法包里缺少比如稳定性检测、季节性差分和自动搜索参数空间的功能，需要自己依照数学公式实现。

更多相关内容
• 季节性时间序列建模与预测 ，孟玲清，王晓雨，所谓所谓时间序列，就是各种社会、经济、自然现象的数量指标按照时间次序排列起来的统计数据，其有多种构成因素，每种因素对系统
• 测试序列稳定：看以看到整体的序列并没有到达稳定要求，要将时间序列转为平稳序列，有如下几种方法：DeflationbyCPILogarithmic（取对数）FirstDifference（一阶差分）SeasonalDifference（季节差分）...
• 用SARIMA预测周游客数，时间序列的周期应该为多少
• 希望我整理的内容对路过的你有所帮助，点赞或评论，都是相互的鼓励~ 【问题】根据下图中某啤酒生产企业2010-2015年各季度的销售量数据，预测2016年各季度产量 ...可以认定啤酒销售量序列是一个含有季节性成分和趋...

希望我整理的内容对路过的你有所帮助，点赞或评论，都是相互的鼓励~

【问题】根据下图中某啤酒生产企业2010-2015年各季度的销售量数据，预测2016年各季度产量

1. 绘制时间序列图，观察啤酒销售量的构成要素

从上图可以明显看出，啤酒销售量具有明显季节成分，而且后面年份销量比前面年份高，因此其中含有趋势成分，但其周期性难以判断。可以认定啤酒销售量序列是一个含有季节性成分和趋势成分的时间序列。

2. 确定季节成分，计算季节指数

2.1 计算移动平均值

-- 对于季节数据，从2010年1季度开始，每4个季度计算4项移动平均，如：

 年份/季度 4项移动平均计算 4项移动平均值 4项移动平均    对应季度 2010/1, 2010/2,2010/3, 2010/4 (25.0+32.0+37.0+26.0) / 4 30.00 2.50 2010/2,2010/3, 2010/4, 2011/1 (32.0+37.0+26.0+30.0) / 4 31.25 3.50 2010/3, 2010/4, 2011/1, 2011/2 (37.0+26.0+30.0+38.0) / 4 32.75 4.50 … … … …

这里出现的问题是，计算出的4项移动平均，没有对应着具体的某个季度，而是在季度之间！

为了解决这个问题，需要进行中心化处理。

-- 对计算结果进行中心化处理，也就是再进行一次二项移动平均，得出中心化移动平均值CMA。

这样处理之后，移动平均值便对应具体季度。思路如下（给我自己做的图点赞❤）：

按照上述思路，计算出的中心化移动平均值CMA情况如下：

 年份 时间代码 销售量 4项移动平均 中心化移动平均值 CMA 2010/1 1 25.0 1.5 2 2 32.0 2.5 30.000 3 3 37.0 30.625 3.5 31.250 4 4 26.0 32.000 4.5 32.750 2011/1 5 30.0 33.375 5.5 34.000 2 6 38.0 34.500 6.5 35.000 3 7 42.0 34.875 7.5 34.750 4 8 30.0 34.875 8.5 35.000 2012/1 9 29.0 36.000 9.5 37.000 2 10 39.0 37.625 10.5 38.250 3 11 50.0 38.375 11.5 38.500 4 12 35.0 38.500 12.5 38.500 2013/1 13 30.0 38.625 13.5 38.750 2 14 39.0 39 14.5 39.250 3 15 51.0 39.125 15.5 39.000 4 16 37.0 39.375 16.5 39.750 2014/1 17 29.0 40.250 17.5 40.750 2 18 42.0 40.875 18.5 41.000 3 19 55.0 41.250 19.5 41.500 4 20 38.0 41.625 20.5 41.750 2015/1 21 31.0 41.625 21.5 41.500 2 22 43.0 41.875 22.5 42.250 3 23 54.0 23.5 4 24 41.0

2.2 计算季节比率

销售量 同 中心化移动平均值CMA 的比值 = 季节比率

在乘法模型中，季节指数是以其平均数等于100%为条件而构成的，它反映了某一季度的数值占全年平均数值的大小。

这里，我们计算出的四个季节比率的平均数为0.9963，不等于1，需进行调整。

2.3 季节指数调整

将每个季节比率的平均值除以四个季节比率的总平均值，得到季节指数

从季节指数变动图可以看出，啤酒销售量的旺季是3季度，淡季是1季度。

3. 分离季节成分

将实际销售量分别除以相应的季节指数，将季节成分从时间序列中分离出去，得到分离季节成分的序列。

4. 建立预测模型

剔除季节成分后，可以观察到啤酒销量有明显的线性增长趋势。用一元线性模型进行回归分析，得到分离季节因素后的序列对应的线性趋势方程为：$\widehat{Y_{t}} = 30.6067 + 0.5592 * t$

5. 预测2016年度销量

根据趋势方程，带入t=25，可以求得2016年1季度销售量（不含季节因素），再乘以对应的季度指数，就可以求得最终的销售量预测值。

将实际销售量和最终预测值进行做图比对，可以看出，预测效果非常好。

展开全文
• ARMA是一种平稳时间序列模型，即均值和协方差不随时间的平移而改变。 ARMA有三种类型 AR序列 MA序列 ARMA序列 但是由于ARMA只能处理平稳序列，而现实中的问题往往有趋势或周期等。为了得到平稳序列，我们对...

ARMA是一种平稳时间序列模型，即均值和协方差不随时间的平移而改变。
ARMA有三种类型

AR序列

MA序列

ARMA序列

但是由于ARMA只能处理平稳序列，而现实中的问题往往有趋势性或周期性等。为了得到平稳序列，我们对数据进行差分运算，使得新序列成为平稳序列，就能够进行ARMA分析，因此ARIMA模型，是在ARMA的基础上多了差分运算，使得其能够处理的序列范围增加了。
ARIMA序列

例题1：

clc,clear
a = [17.0 16.6 16.3 16.1 17.1 16.9 16.8 17.4 17.1 17.0
16.7 17.4 17.2 17.4 17.4 17.0 17.3 17.2 17.4 16.8
17.1 17.4 17.4 17.5 17.4 17.6 17.4 17.3 17.0 17.8
17.5 18.1 17.5 17.4 17.4 17.1 17.6 17.7 17.4 17.8
17.6 17.5 16.5 17.8 17.3 17.3 17.1 17.4 16.9 17.3
17.6 16.9 16.7 16.8 16.8 17.2 16.8 17.6 17.2 16.6
17.1 16.9 16.6 18.0 17.2 17.3 17.0 16.9 17.3 16.8
17.0 16.6 16.3 16.1 17.1 16.9 16.8 17.4 17.1 0
17.3 17.4 17.7 16.8 16.9 17.0 16.9 17.0 16.6 16.7
16.8 16.7 16.4 16.5 16.4 16.6 16.5 16.7 16.4 16.4
16.2 16.4 16.3 16.4 17.0 16.9 17.1 17.1 16.7 16.9
16.5 17.2 16.4 17.0 17.0 16.7 16.2 16.6 16.9 16.5
16.6 16.6 17.0 17.1 17.1 16.7 16.8 16.3 16.6 16.8
16.9 17.1 16.8 17.0 17.2 17.3 17.2 17.3 17.2 17.2
17.5 16.9 16.9 16.9 17.0 16.5 16.7 16.8 16.7 16.7
16.6 16.5 17.0 16.7 16.7 16.9 17.4 17.1 17.0 16.8
17.2 17.2 17.4 17.2 16.9 16.8 17.0 17.4 17.2 17.2
17.1 17.1 17.1 17.4 17.2 16.9 16.9 17.0 16.7 16.9
17.3 17.8 17.8 17.6 17.5 17.0 16.9 17.1 17.2 17.4
17.5 17.9 17.0 17.0 17.0 17.2 17.3 17.4 17.4 17.0
18.0 18.2 17.6 17.8 17.7 17.2 17.4 0 0 0]';

a = nonzeros(a);
r11 = autocorr(a);
r12 = parcorr(a);
da = diff(a);
r21 = autocorr(da)
r22 = parcorr(da)
n = length(da);
k = 0;
for i = 0:3
for j = 0:3
if i == 0 & j == 0
continue
elseif i == 0
ToEstMd = arima('MALags', 1:j, 'Constant', 0);
elseif j == 0
ToEstMd = arima('ARLags', 1:i, 'Constant', 0);
else
ToEstMd = arima('ARLags', 1:i, 'MALags', 1:j, 'Constant', 0);
end
k = k + 1;
R(k) = i;
M(k) = j;
[EstMd, EstParamCov, logL, info] = estimate(ToEstMd, da);
numParams = sum(any(EstParamCov));
[aic(k), bic(k)] = aicbic(logL, numParams, n);
end
end
fprintf('R, M, AIC, BIC的对应值如下\n ');
check = [R', M', aic', bic']
r = input('输入阶数R=');
m = input('输入阶数M=');
ToEstMd = arima('ARLags', 1:r, 'MALags', 1:m, 'Constant', 0);
[EstMd, EstParamCov, logL, info] = estimate(ToEstMd, da);
dx_Forecast = forecast(EstMd, 10, 'Y0', da)
x_Frecast = a(end) + cumsum(dx_Forecast)



我代码敲出来结果有一点偏差，可能数据打错了几个吧，太多了，不想看了。

例题2

clc,clear
x = [9.40 8.81 8.65 10.01 11.07 11.54 12.73 12.43 11.64 11.39 11.1 10.85
10.71 10.24 8.48 9.88 10.31 10.53 9.55 6.51 7.75 7.8 5.96 5.21
6.39 6.38 6.51 7.14 7.26 8.49 9.39 9.71 9.65 9.26 8.84 8.29
7.21 6.93 7.21 7.82 8.57 9.59 8.77 8.61 8.94 8.4 8.35 7.95
7.66 7.68 7.85 8.53 9.38 10.09 10.59 10.83 10.49 9.21 8.66 8.39
8.27 8.14 8.71 10.43 11.47 11.73 11.61 11.93 11.55 11.35 11.11 10.49
10.16 9.96 10.47 11.70 10.1 10.37 12.47 11.91 10.83 10.64 10.29 10.34];
x = x';
x = x(:);
s = 12; % 周期
n = 12; % 预测数
m1 = length(x); % 原始数据数
for i = s+1:m1
y(i-s) = x(i) - x(i-s); % 进行周期差分变换
end
w = diff(y); % 消除趋势性差分运算
m2 = length(w); % 计算最终差分后数据个数
k = 0; % 模型个数

for i = 0:3
for j = 0:3
if i == 0 & j == 0
continue
elseif i == 0
ToEstMd = arima('MALags', 1:j, 'Constant', 0);
elseif j == 0
ToEstMd = arima('ARLags', 1:i, 'Constant', 0);
else
ToEstMd = arima('ARLags', 1:i, 'MALags', 1:j, 'Constant', 0);
end
k = k + 1;
R(k) = i;
M(k) = j;
[EstMd, EstParamCov, logL, info] = estimate(ToEstMd, w'); % 模型拟合
numParams = sum(any(EstParamCov)); % 计算拟合参数个数
[aic(k), bic(k)] = aicbic(logL, numParams, m2);
end
end
fprintf('R, M, AIC, BIC的对应值如下\n ');
check = [R', M', aic', bic']
r = input('输入阶数R=');
m = input('输入阶数M=');
ToEstMd = arima('ARLags', 1:r, 'MALags', 1:m, 'Constant', 0); % 指定模型结构
[EstMd, EstParamCov, logL, info] = estimate(ToEstMd, w'); % 拟合模型
w_Forecast = forecast(EstMd, n, 'Y0', w') % 计算12步预测值，注意已知数据为列向量
yhat = y(end) + cumsum(w_Forecast) % 求一阶差分还原值
for j = 1:n
x(m1+j) = yhat(j) + x(m1 + j - s);% 求x的预测值
end
xhat = x(m1 + 1: end) % 截取n个预报值


分析：

这题我和他的结果是对的上的。

使用ARIMA模型简单来说就是先差分使得序列成为平稳序列，然后就是遍历AR和MA的参数的事了，看aic和bic越小的，就选哪个参数就好了。

本文参考的是司守奎，孙兆亮主编的数学建模算法与应用（第二版）

展开全文
• 即我们先对数据的趋势作预测预测完了之后再在原来数据的基础上加上这个趋势值就能够较好地预测出需要的值了。 一阶差分指数平滑法 例题： clear, clc yt = [24 26 27 30 32 33 36 40 41 44]'

差分指数平滑法

差分就是后一个数据减去前一个数据得到的差值，为什么说如果有直线趋势便改变了趋势呢？因为例如y = a * x + b，那么每个差分的值就相当于a了，当然因为误差存在，肯定是围绕a上下波动的，因此便转换为一种围绕a变化的趋势而不是直线趋势了，便改变了数据的趋势。即我们先对数据的趋势作预测，预测完了之后再在原来数据的基础上加上这个趋势值就能够较好地预测出需要的值了。

一阶差分指数平滑法

例题：

clear, clc
yt = [24 26 27 30 32 33 36 40 41 44]';
n = length(yt);
alpha = 0.4;
dyt = diff(yt); % 向前求一阶差分，即后面的数-前面数
dyt = [0 ; dyt]; % 第一个数没得求，用0补
dyhat(2) = dyt(2); % 指数平滑的初始值
for i = 2:n
dyhat(i + 1) = alpha * dyt(i) + (1 - alpha) * dyhat(i);
end
for i = 1:n
yhat(i + 1) = dyhat(i + 1) + yt(i);
end
yt
dyt
yhat' % 预测值


具有季节性特点的时间序列预测
季节可以指自然季节，也可以指销售季节等等。例如羽绒服肯定是冬季销量夏季的。下面介绍一种季节系数法，计算步骤如下：

例题：

clc,clear
format long g
data = [137920 186742 274561 175433
142814 198423 265419 183521
131002 193987 247556 169847
157436 200144 283002 194319
149827 214301 276333 185204];
[m, n] = size(data);
a_mean = mean(mean(data)); % 计算所有数据的算术平均值
aj_mean = mean(data); % 计算同季节的算术平均值
bj = aj_mean/a_mean % 计算季节系数
w = 1:m;
yhat = w * sum(data,2)/sum(w) % 预测下一年的年加权平均值，这里是求行和
yjmean = yhat/n % 计算预测年份的季节平均值
yjhat = yjmean * bj % 预测年份的季节预测值
format


本文参考的是司守奎，孙兆亮主编的数学建模算法与应用（第二版）

展开全文
• 利用SPSS 和Matlab 进行时间序列预测1.移动平均和滑动平均计算例1：下表给出了某地区1990～2004年粮食产量数据(表1)。试分别用Matlab 和SPSS 软件，对该地区的粮食产量进行移动平均和和滑动平均计算。表1 某地区1990...
• 本篇文章将总结时间序列预测方法，并将所有方法分类介绍并提供相应的python代码示例，以下是本文将要介绍的方法列表： 1、使用平滑技术进行时间序列预测 指数平滑 Holt-Winters 法 2、单变量时间序列预测 自...
• 这段是实现预测的关键代码，需要大家根据个人情况编写剩余代码 参数定义 ： x为原始序列数据 y为进行周期差分变换后的x w为差分运算之后的y 话不多说，直接上代码，后面有分段的代码讲解 figure(1); plot...
• 序列分解1、非季节性时间序列分解 移动平均MA(Moving Average)①SAM（Simple Moving Average） 简单移动平均，将时间序列上前n个数值做简单的算术平均。 SMAn=(x1+x2+…xn)/n②WMA(Weighted Moving Average) ...
• 时间序列预测，非季节性ARIMA及季节性SARIMA

万次阅读 多人点赞 2019-03-24 21:55:00
我们将首先介绍和讨论自相关，平稳性和季节性的概念，并继续应用最常用的时间序列预测方法之一，称为ARIMA。 介绍 时间序列提供了预测未来价值的机会。基于以前的价值观，可以使用时间序列来预测经济，天气和...
• 季节性自回归综合移动平均（SARIMA）或季节性ARIMA是ARIMA的一个扩展，它明确支持具有季节性分量的单变量时间序列数据，它增加了三个新的超参数来指定序列季节性成分的自回归（AR）、差分（I）和移动平均（MA），...
• 在真实世界季节性时间序列数据集的实验中，我们的模型始终优于传统的线性模型和 GRU 递归神经网络。 【本文剩余部分内容的组织】 本文的其余部分安排如下。第 2 节概述了相关背景，包括具有代表性的自回归方法和...
• 时间序列预测方法最全总结！

万次阅读 多人点赞 2021-03-12 00:15:38
时间序列预测就是利用过去一段时间的数据来预测未来一段时间内的信息，包括连续型预测（数值预测，范围估计）与离散型预测（事件预测）等，具有非常高的商业价值。需要明确一点的是，与回归分析预测模型...
• 这本来是我回答的一个问题：有什么好的模型可以做高精度的时间序列预测呢？ - BINGO Hong的回答 - 知乎 https://www.zhihu.com/question/21229371/answer/533770345但觉得在那个答案下一直更新好麻烦，干脆就移到...
• 该工具箱提供了深入的文档系统和联机帮助，并且其中包含许多演示，这些演示将指导您完成时间序列建模的过程。 Matlab的ECOTOOL工具箱已发布在PLOS ONE中，您可以在其中找到一些示例以及工具箱概述（）。
• 本文介绍了布朗单一参数线性指数平滑法、霍特双参数指数平滑法、布朗三参数指数平滑法及温特线性和季节性指数平滑法四种时间序列平滑法在产品产量预测中的应用，并对这四种方法的适用范围进行了总结。
• 在matlab中实现ARIMA时间序列预测。函数形式如下： function [result] = ARIMA_algorithm(data, Periodicity, ACF_P, PACF_Q, n) 其中data为预测所用的数据，为一维列向量；Periodicity为数据的周期；ACF_P和PACF_Q...
• 这是一个提供时间序列预测功能的 Java 开源库。 它是加法模型的实现。 该库由 Workday 的 Syman 团队发布，用于支持某些 Workday 产品中的基本时间序列预测功能。 如何使用 为了使用这个库，您需要提供输入时间序列...
• 在我们的研究中，我们考虑了对印度旁遮普省降雨数据进行统计分析的季节性和周期性时间序列模型。 在本研究论文中，我们应用季节性自回归综合移动平均和周期自回归模型来分析旁遮普省的降雨数据。 为了评估模型识别...
• 基于增量线性递减趋势和多正弦季节性加权组合的时间序列预测，劳稳超，张雨浓，本文基于趋势、季节性以及其他组成成分的加权组合模型，提出了一种新型的成分加权组合（weighted-combination-of-components，WCC）...
• 1.季节性时间序列 包含：长期趋势Trend，季节趋势Seasonal，周期循环Circle，随机项Random 这里分解为相加模型X=T+S+C+R   在对时间序列进行分解之前，应该对序列进行检验：（下次写） 2.decompose()函数 将...
• 时间序列季节性分解 可将一个序列分解成一个季节性成分、一个组合趋势和循环成分、一个“误差”成分。 时间序列分两种模型：乘法模型和加法模型，对上面3个分量做乘法或加法。 如果季节性波动的幅度不随着序列水平...
• SARIMA(p,d,q)(P,D,Q,s) 季节性自回归移动平均模型，结构参数有七个 AR(p) 自回归模型，即用自己回归自己。基本假设是，当前序列值取决于序列的历史值。p 表示用多少个历史值来回归出预测值。 要确定初始 p，需要...
• 1.时间序列预测简介 时间序列是在定期时间间隔内记录度量的序列。 根据频率，时间序列可以是每年（例如：年度预算），每季度（例如：支出），每周（例如：销售数量），每天（例如天气），每小时（例如：股票价格）...
• 概述： 本文演示了 11 种不同的经典时间序列预测方法，它们分别是1）自回归（AR） 2）移动平均线3) 自回归移动平均线4) 自回归综合移动平均线 (ARIMA) 5) 季节性自回归综合移动平均 (SARIMA) 6) 带外生回归量的季节...

...