用Excle做数据分析与预测

假期收到导师布置的作业，用时间序列季节指数的方法预测虫情，由于数据量比较少，用python的话有点杀鸡用牛刀了，用Excle简单、方便、又快捷，直接起飞。

时间序列

时间序列（或称动态数列）是指将同一统计指标的数值按其发生的时间先后顺序排列而成的数列。时间序列分析的主要目的是根据已有的历史数据对未来进行预测。经济数据中大多数以时间序列的形式给出。根据观察时间的不同，时间序列中的时间可以是年份、季度、月份或其他任何时间形式。

模型建立

获取数据

选取2017年—2019年某地各月份的昆虫种类及数量（导师瞎给的数据）

数据处理

在时间序列的模型中：

Y=T（长期趋势）*S（季节指数）CI

首先对数据中的昆虫数量Y进行四项平均

做完四项平均后还要对所得到的数据进行中心平均得到TC的值

求出TC以后我们根据公式就能得出带不规则变动的季节指数SI（SI=Y/TC）：

当我们求出SI以后，因为我们是需要得到季节指数S，所以我们需要把里面的不规则变动I去掉，所以我们还需要对SI进行处理：

求出各年同季观察值平均值A
求出历年总季度平均值B
季节指数C=A/B
最后我们就能够得出我们所需要的去除不规则变动I的季节指数S
放入图表中

进行数据分析

点击Excle上方的数据

找到右边的数据分析

选择回归分析

把需要进行分析的数据范围填好，置信度一般选择%95残差项根据自己的需求选择

对未来进行预测

方差分析表，在这里面一共有五个数，分别是df（自由度）、SS（样本数据平方和）、MS（样本数据平均平方和）、F（F统计量的值）、Significance
F（P值）

从图中我们就可以得出我们所计算出来的回归方程，在这个里面的t值和P值也都能够直接看到，最后我们得出回归方程它的表达式为：

Y=-4.43706x+1037.758

由上公式把长期趋势T求出：

当把表格完成以后，我们需要的预测值就等于2019年各个季度的季节指数S*长期趋势T:

分析与预测就做完了

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

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

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

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

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

2.1 计算移动平均值

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

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

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

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

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

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

2.2 计算季节比率

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

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

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

2.3 季节指数调整

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

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

3. 分离季节成分

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

4. 建立预测模型

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

5. 预测2016年度销量

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

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

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

趋势稳定性检测：Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test

null-hypothesis: 时间序列趋势稳定。significance: 0.05. 选择KPSS和不是Dickey Fuller由于KPSS的null-hypothesis是趋势稳定，所以接受条件相对于Dickey Fuller更宽松，差分阶数更少。
$y_t={\beta }^{\prime }{D}_t+{\mu }_t+u_t$
${\mu }_t={\mu }_{t-1}+{ϵ}_t,$
${ϵ}_t\sim WN(0,{\sigma }_{ϵ}^2)$

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

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

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

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

时间序列yi的拟合：
$y_i=\mu +x_i^{\prime }\beta +S_i+e_i$
$S_i={\sum }_{j=1}^{m/2}{f}_{ji}^{\prime }{\gamma }_j$
其中
$f_{ji}^{\prime }=(cos((j/q)\pi i),sin((j/q)\pi i))$

LM statistic:
$$LM=\frac{1}{n^2}trace(({\hat{\Omega }}^f{)}^{-}1\sum _{i=1}^n{\hat{F}}_i{\hat{F}}_i^{\prime })$$

其中
${\hat{F}}_i={\sum }_{t=1}^i{f}_t{\hat{e}}_t$
${\hat{\Omega }}^f={\sum }_{k=-m}^{m}w(\frac{k}{m})\frac{1}{n}{\sum }_i{f}_{i+k}{\hat{e}}_{i+k}{f}_i^{\prime }{\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
pd.DataFrame(seasonalDummy(tsArray, frequency)).head(2)

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=\mu +x_i^{\prime }\beta +S_i+e_i$
$S_i={\sum }_{j=1}^{m/2}{f}_{ji}^{\prime }{\gamma }_j$

其中$f_{ji}^{\prime }=(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：$$\Omega =li{m}_{n\to \infty }\frac{1}{n}E(F_n{F}_n^{\prime })$$
在ei可能有serial correlation的情况下可以用estimate
$$\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(\cdot )$ 是kernel function，来自rob j. Hyndman的forecast包:

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

${\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 = 没有单位根）：
$LMstatistic=\frac{1}{n^2}trace((A^{\prime }{\hat{\Omega }}^{f}A{)}^{-1}{A}^{\prime }{\sum }_{i=1}^n\widehat{F_i}{\widehat{F_i}}^{\prime }A)$

其中${\sum }_{i=1}^n\widehat{F_i}$ 和${\sum }_{i=1}^n{\widehat{F_i}}^{\prime }$由前面步骤中的 $\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=\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)=\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}_{tot}={\sum }_i(y_i-\bar{y}{)}^2$
$S{S}_{res}={\sum }_i(y_i-{\hat{y}}_i{)}^2$
$R^2=1-\frac{S{S}_{res}}{S{S}_{tot}}$

RMSE:
$RMSE=\sqrt{\frac{{\sum }_{t=1}^n({\hat{y}}_t-y_t{)}^2}{n}}$

MAE:
$MAE=\frac{{\sum }_{t=1}^n|{\hat{y}}_t-y_t|}{n}$

MAPE:
$MAPE=\frac{100}{n}{\sum }_{t=1}^n|\frac{{\hat{y}}_t-y_t}{y_t}|$

sMAPE:
$SMAPE=\frac{100}{n}{\sum }_{t=1}^n\frac{|{\hat{y}}_t-y_t|}{(|{\hat{y}}_t|+|y_t|)/2}$

MASE without seasonality:
$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:
$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-\alpha )\ast W_{n-1}+\alpha \ast y_n$
$W_0=y_0$
$\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

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

$ynew=log(y)=log(Trend\ast Seasonality\ast Residual)=log(Trend)+log(Seasonality)+log(Residual)$

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

ARMA模型范式：

$$X_t-{\alpha }_1{X}_{t-1}-\cdots -{\alpha }_{p^{\prime }}{X}_{t-p^{\prime }}={ϵ}_t+{\theta }_1{ϵ}_{t-1}+\cdots +{\theta }_q{ϵ}_{t-q}$$

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

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

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

总结

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