时间序列分析（一）单变量时间序列分析

Background

在我们对单变量时间序列进行分析时，为了简化模型，我们往往需要预先作出一些假设。
我们令随机变量序列 \(\{Y_t: t=0,\pm1, \pm2, \pm3,···\}\)为一个随机过程，并将其作为我们所观测的时间序列，整个序列可以看作是 Y 的联合分布。对于一般的时间序列模型，我们可以采用一些相应的统计量来描述整个时间序列（如，Y的联合分布为多元正态的，则一阶和二阶矩就可以确定整个模型），这些模型往往较为简单，但是十分实用。

基本概念

1.1基本统计量

均值、方差、协方差

基础统计量在实际单变量时间序列分析中非常有用，它们不仅可以对时间序列的周期性、平稳性等方面进行分析，甚至还可以辅助确定一些时间序列模型的参数。

均值

对随机过程 \(Y_t\) ，均值定义如下：
$$\mu _t = E(Y _t), t = 0, \pm 1, \pm 2, \pm 3, \cdots $$

自协方差

$$ \gamma _{t,s} = Cov(Y _t, Y _s), t,s = 0, \pm 1, \pm 2, \pm 3, \cdots$$
其中 \(Cov(Y _t, Y _s) = E[(Y _t - \mu _t)(Y _s - \mu _s)] = E(Y _tY _s) - \mu _t\mu _s\)（根据均值定义消去均值乘积）。

自相关函数

$$ \rho _{t,s} = Corr(Y _t,Y _s), t,s = 0, \pm 1, \pm 2, \pm 3, \cdots $$
$$ Corr(Y _t,Y _s) = \frac{Cov(Y _t, Y _s)}{\sqrt{Var(Y _t)Var(Y _s)}} = \frac{\gamma _{t,s}}{\sqrt{\gamma _{t,t}\gamma _{s,s}}}$$

1.2 平稳性

对于一切时间间隔k，如果在时间点 \(t_1 , t_2, t_3, \cdots , t_n\)下均有，\(Y_{t_1}, Y_{t_2}, Y_{t_3}, \cdots , Y_{t_n}\)与\(Y_{t_1 -k}, Y_{t_2 -k}, Y_{t_3 -k}, \cdots , Y_{t_n -k}\) 的联合分布相同，则过程\(\{Y_{t}\}\)称为严平稳的。
而若序列是平稳性的，则：
当 n=1 时，\(Y_{t}\) 与\(Y_{t-k}\) 的分布相同，所以Y序列的每一点具有相同的分布，因此序列均值与方差恒为常数。
进而任意时间点t来说，由于序列的每一点均值方差均相同，所以我们可以化简自相关函数，即
$$ \rho _{k} = Corr(Y _t,Y _{t-k}), \gamma _{k} = Cov(Y _t, Y _{t-k})$$
严平稳性是一个非常强的假设，因此我们引出了弱（二阶矩）平稳性：

1. 均值函数在所有时间上恒为常数。
2. 对于所有时间点t和时间间隔k，有\(\gamma_{t,t-k}=\gamma_{0,k}\)。
   在实际应用中，\(\gamma_{k}\)是一个非常有用的工具，它可以帮助我们判断时间序列的周期性，在某些时间序列模型中，我们使用\(\gamma_{k}\)来帮我们判断模型参数，这些应用即使在非平稳的状态下也有一定借鉴意义。

基本模型

在进行单变量时间序列分析时，我们往往需要对时间序列模型作出一些假设。最常见的假设为线性叠加模型，即模型中各个组成部分是以线性叠加的方式组合到一起的。对于能预测的单变量模型来说，我们往往认为其具有平稳周期性，和非平稳部分。而非平稳部又可以分为趋势部分与局部异常部分。

周期性

在单变量时间序列模型中，我们往往假设模型是具有周期性的，在对用户数据建模时，这一点往往是成立的（比如：分析用户PV、UV等，往往具有明显的7天周期与一年周期，在这不同周期中，用户数据整体波动出现明显相似的变化）。
对于这一部分的建模，有很多方式，如：简单季节性模型根据周期性的经验估计直接建立模型，ARIMA中利用假设过去周期性的时间节点对当前时间点存在残留影响来建立模型，谱分析对时间序列进行傅立叶变换从而拟合周期部分等等。
不论是使用什么方式，在进行周期性估计时，我们往往不考虑序列的整体趋势，也就是说我们认为周期性部分是平稳的（往往只需要弱平稳即可）。

非平稳部分

趋势

对于时间序列模型，我们往往可以将其拆分为两部分，一部分是趋势拟合部分，另一部分是周期性拟合部分。对于趋势部分，一些模型是进行单独拟合的，例如：使用线性模型，或根据经验，使用相应的非线性模型，也有一些模型使用了差分的方式，消除趋势对周期性判断对影响（例如：ARIMA）。
当然也有一些模型不需要对趋势进行单独拟合，模型可以自动学习出整体趋势与周期特征（例如：使用深度学习方法对序列进行拟合）。

干预与异常

干预分析（Intervention Analysis）是指在时间序列中，因为在某一时刻发生了自然产生或人为施加的影响，导致整体序列产生了一定时间内均值或趋势发生了相应的变化，例如：外部广告投放导致一定时间内某应用UV突然增长，一段时间后稳定用户有一个明显的数量提升、商家在某一天投放广告导致商品购买量提升，广告停止投放，几天内购买量下降，但仍高于投放前，直至平稳、双十一大促，大促时购买量激增，但大促后用户消费明显疲软，一段时间后恢复正常等。在实际数据中，不同的干预可能有不同的特征，下文将介绍几种简单常见的干预拟合方式。

产生干预时，我们的干预部分用\(S_t ^{(T)}\)表示，其中T为干预产生的时间点。则我们可以定义脉冲函数即：
$$P _t ^{(T)}=S _t ^{(T)}-S _{t-1} ^{(T)}$$

阶梯干预

我们可以将阶梯干预用如下阶梯函数表示：
$$
S _t=\begin{cases}
1,\quad x\geq T \\
0,\quad x<T
\end{cases}
$$
对应的脉冲函数\(P_t ^{(T)}\)，当t=T时\(P_t ^{(T)}\)等于1，其他时间为0。
此时我们干预对*时间序列均值的叠加部分\(m_t\)*为：
$$m _t = \omega S _t ^{(T)} \tag{1.1}$$
也就是说在T时刻激增并保持不变，这种干预往往可以看作是一种可加异常，也就是说当前时刻的异常并不影响之后的序列。

逐渐平稳的阶梯干预

阶梯干预的形式非常理想，但实际应用中，全部影响作用其实会持续很长一段时间最终趋于稳定。我们可以对\(m_t\)建立一个AR(1)模型（后续会对其进行介绍）：
$$
m _t=\begin{cases}
\omega\frac{1-\delta ^{t-T}}{1-\delta},\quad x\geq T \\
0,\quad other
\end{cases}
$$
其中0<\(\delta\)<1。此时
均值增量在会在一段时间内激增最终收敛到\(\frac{\delta}{1-\delta}\)，其中达到极限的一半所用时间\(log(0.5)/log(\delta)\)被称为干预的半衰期

逐渐消失的干预

对于一些干预，在干预瞬间，序列均值会激增，但随着时间推移这种干预影响会逐渐消失。令T时刻干预产生，这时当t=T时\(P_t ^{(T)}\)等于1，其他时间为0。我们可以将这种干预混入我们的模型例如，混入AR(1)模型中：
$$m _t=\delta m _{t-1} + \omega P^{(T)} _{t-1}$$

此时T时刻产生的干预成指数级衰减。

判断与检验

受到干预或异常的数据，有很多表示方法，针对不容模型，可能表示形式也不尽相同，而上述方法对于ARIMA模型是适用的，对于其他模型则要因模型而异。面对可能的异常，我们需要进行判断和检验，常用的方法可以使用先猜测异常，然后修正模型，进行假设检验，通过反复修正，最终直至不再发现异常为止。
在进行检验时，常用多重检验，可以使用保守的Bonferroni矫正（如果在同一数据集上同时检验n个独立的假设，那么用于每一假设的统计显著水平，应为仅检验一个假设时的显著水平的1/n。）控制多重检验的总体误差率。

ARIMA（自回归差分滑动平滑方法）

对于时间序列\(\{Y _t\}\)，若我们令\(\{e _t\}\)代表一系列未观测到的白噪声序列数据，若我们将其看作一系列独立同分布的均值为0的随机变量。则对于一个序列而言，我们可以将其转化为过去时间点序列值残留的影响与当前时间点白噪声的叠加。
因此，对于我们可以将上述过程，简化成一个简单的线性过程：
$$Y _t=e _t + \psi _1 e _{t-1} + \psi _2 e _{t-2} + \cdots \tag{2.1}$$
也就是说当前时间点是过去时间点的残留影响加上当前时间点点白噪声。对于这一模型，若表达式右边是一个无穷级数项，则要求：
$$\sum _ {i=1} ^ \infty {\psi _i ^2} < \infty$$
也就是说我们希望整体序列整体是收敛的。显然，我们希望考虑到过去的影响，但又不希望过去对当前的影响是发散的。

##MA
对于公式2.1，若\(\psi \)的数量有限且不为0个，则我们称：
$$Y _t=e _t - \theta _1 e _{t-1} - \cdots - \theta _n e _{t-q} $$

为q阶滑动平均过程，简记MA(q)，其中为方便后续公式表达，我们用\(-\theta\)。
对于MA(q)过程我们有：
$$\gamma _0 = (1 + \theta _1 ^2 + \theta _2 ^2 + \cdots + \theta _n ^2)\delta _e ^2 $$
$$
\rho _k=\begin{cases}
\frac{-\theta _k + \theta _1 \theta _{k+1} + \theta _2 \theta _{k+2} + \theta _{q-k} \theta _{q}}{1+\theta _1 ^2+\theta _2 ^2 + \cdots + \theta _q ^2},\quad k=1,2,\cdot,q \\
0,\quad k>q
\end{cases}
$$

AR

对于时间序列\(\{Y _t\}\)若满足方程：
$$Y _t=\phi _1 Y _{t-1} + \phi _2 Y _{t-2} + \cdots + \phi _n Y _{t-p} + e _t $$
则，我们称其为p阶自回归过程。从方程我们可以看出，自回归过程是用自身序列做回归变量。对于t时刻的\(Y _t\)来说其值等于过去p个时间点与当前时间点噪声\(e _t\)的加权叠加，其中\(e _t\)要求与t时刻前的Y独立。从等式2.2我们可以看出，通过将右侧Y展开，我们可以还原出等式2.1的形式，但参数的表达可能十分复杂，同时该序列项数可能趋于无穷大。
对于MA模型来说，我们较容易保证其平稳，而对于AR模型，平稳显然需要满足条件。令，
$$\phi(x)=1 - \phi _1 x - \phi _2 x^2 - \cdots - \phi _p x ^ p $$
$$0=1-\phi _1 x - \phi _2 x^2 - \cdots - \phi _p x ^ p $$
为AR模型的特征多项式和特征方程。若特征方程存在平稳解，当且仅当AR特征方程每一个根的绝对值都大于1。因此有如下必要不充分的条件：
$$\phi _1 + \phi _2 + \cdots + \phi _p < 1 $$
$$| \phi _p |< 1$$
假定序列平稳且均值为0。则，我们可以在AR模型方程两边同时乘以\(Y _{t-k}\)并求期望，再除以\(\gamma _0\)可以获得以下重要递推关系：

$$\rho _k=1 - \phi _1 \rho _{k-1} - \phi _2 \rho _{k-2} - \cdots - \phi _p \rho _{k-p} $$
分别将\(k = 1,2,\cdots,p\)带入上述方程，则可以得到方程组（Yule-Walker方程组）通过求解相应方程，我们可以求的任意自相关系数\(\rho _k\)。
对于AR模型而言，若序列平稳，通过对自回归进行逐层展开，我们可以将其还原为式2.1形式。
而对于MA模型，我们可以通过\(Y _{t-k}\) 迭代消去序列中的第k项\(e _{t-k}\),从而将MA模型转化为AR模型。若原时间序列是无穷的，那么我们最终可能转化为一个无穷阶的AR模型。这种MA可逆的转化是有条件的，需要MA特征方程：
$$0=1-\theta _1 x - \theta _2 x^2 - \cdots - \theta _p x ^ p $$
的根的模大于1（类似AR模型平稳性条件）。

ARMA(p,q)模型

若序列心中一部分是自回归，一部分是滑动平均的，那么我们可以得到一个相当普遍的时间序列模型。如下：
$$Y _t=\phi _1 Y _{t-1} + \phi _2 Y _{t-2} + \cdots + \phi _n Y _{t-p} + e _t - \theta _1 e _{t-1} - \theta _2 e _{t-2} - \cdots - \theta _q e _{t-q}$$
对于该模型而言，MA部分我们可以保证解的平稳，因此对于模型整体平稳性，我们只需要保证AR部分过程平稳，即保证模型中AR特征方程根的模大于1。
此时模型可以由下式决定：

$$
\begin{cases}
\psi _0 = 1\\
\psi _1 = -\theta _1 + \phi _1\\
\psi _2 = -\theta _2 + \phi _2 + \phi _1\theta _1\\
\cdots\\
\psi _j = -\theta _j + \phi _p\theta _{j-p} + \phi _{p-1}\theta _{j-p+1} + \cdots + \phi _1\theta _{j-1} \\
\end{cases}
$$

这里我们将模型中AR部分展开得到公式2.1的形式，其中参数\(\psi\)为2.1中系数。
对于ARMA模型，我们通常要求AR和MA分别满足平稳和可逆两个条件（上文提到，各自特征方程根模大于一）

非平稳模型

上述描述的AR与MA模型都存在平稳的前提，而实际上，我们的时间序列往往具有趋势。这时我们可以不再使用的时间序列\(\{Y _t\}\)本身来进行预测，而使用查分的方式用变化率作为序列来进行预测。如：\( \Delta Y _t = Y _t - Y _{t-1} \)为序列 \(\{Y _t\}\)的一阶差分。

如果一个时间序列 \(\{Y _t\}\)的d次差分 \(W _t = \Delta ^d Y _t\)是一个平稳的ARMA(p,q)过程，则称\(\{Y _t\}\)为自回归滑动平均求和模型ARIMA(p,d,q)。通常取d=1或最多2。
考虑ARIMA(p,1,q)，令\(W _t = \Delta ^d Y _t\)，我们有：
$$W _t=\phi _1 W _{t-1} + \phi _2 W _{t-2} + \cdots + \phi _n W _{t-p} + e _t - \theta _1 e _{t-1} - \theta _2 e _{t-2} - \cdots - \theta _q e _{t-q}$$
通过将W用Y替换并合并，可得：

$$Y _t=(1 + \phi _1)W _{t-1} + (\phi _2 -\phi _1)W _{t-2} + \cdots + (\phi _{t-p} -\phi _{t-p-1}) W _{t-p} - \phi _p W _{t-p-1} \\+ e _t - \theta _1 e _{t-1} - \theta _2 e _{t-2} - \cdots - \theta _q e _{t-q}$$
我们称上述形式为模型的差分方程形式。我们分析其AR特征方程：
$$1-(1+\phi _1)x-(\phi _2 - \phi _1)x ^2-\cdots-(\phi _{t-p} -\phi _{t-p-1})x ^p + \phi _p x ^{p+1}=(1-\phi _1 -\phi _2x ^2-/cdots-\phi _px ^p)(1-x)$$
此时 x=1是一个根，因此这个过程非平稳，但其余的根是平稳过程\( \Delta ^d Y _t\)的特征多项式的根。这也说明了为什么当原序列非平稳，我们使用一阶差分来得到平稳序列的原因。

差分是实现平稳较好的方式，出了差分外，我们还可以根据序列特征对序列进行其他的变换。例如：如果发现序列散度随着时间增加，可以尝试使用对数方法对序列进行变换等。

季节性ARIMA

对于普通的ARIMA模型，我们考虑的是当前相近的时间节点对现在的影响，然而实际上数据的周期可能会更长。例如：对于用户数据而言，往往存在以周和年为周期的pattern。对于这种较长跨度的周期性影响，ARIMA模型不能进行很好的表达，这时我们可以引入季节性ARIMA模型。
简单的季节性ARIMA模型我们可以通过经验知识，选出相应周期，考虑序列影响时，只考虑倍数于周期的相应时间节点。例如以s（s大于P和Q）为周期，则我们的AR(P)和MA(Q)以及差分I(D)模型如下：

$$
Y _t = Y _t - \Phi _1Y _{t-s} - \Phi _2Y _{t-2s} - \cdots - \Phi _PY _{t-Ps}
$$

$$
Y _t = e _t - \Theta _1e _{t-s} - \Theta _2e _{t-2s} - \cdots - \Theta _Qe _{t-Qs}
$$

$$
\Delta _s Y _t = Y _t + Y _{t-s}
$$
可以看到AR、MA、以及差分，我们都只考虑不同周期s上相同位置的时间节点（只考虑最近P和Q个周期内的时刻）。注意这里我们使用参数：\(P,Q,\Phi,\Theta\)均为大写以表示季节模型。
简单季节性ARIMA加入引入一个周期s，整体模型可以看作是对时间序列间隔s时刻，从而抽出s个子序列，并使用普通ARIMA模型进行拟合。然而这种方法只考虑了一个固定周期的影响，实际上真实序列相近时刻可能也会对当前时刻造成影响（基本的ARIMA模型的假设），因此我们需要同时考虑季节性影响和临近影响。因此引入了人乘法季节ARIMA模型。
乘法季节ARIMA模型结合了季节和非季节ARIMA，利用类似卷积的思路，将二者结合。例如：我们考虑结合\(MA(1)\) 与 \(MA(1) _s\)两个模型，得到如下特征多项式：
$$
(1 -\theta x)(1-\Theta x ^{12}) = 1 - \theta x - \Theta x ^{12} - \theta\Theta x^{13}
$$
通过特征多项式我们可以写出相应的时间序列表达式：
$$
Y _t = e _t - \theta e _{t-1} - \Theta e _{t-12} - \theta\Theta e _{t-13}
$$

因此我们定义周期为s的乘法季节\(ARMA(p,q)\times(P,Q) _s\)模型，通过分别求AR特征多项式\(\phi\Phi(x)\)和MA特征多项式\(\theta(x)\Theta(x)\)得到相应AR与MA模型序列通过相加合并两个子模型。
对于差分模型而言我们可以定义混合差分操作：
$$
W _t = \Delta ^d \Delta _s ^D Y _t
$$
因此我们可以得到新的季节ARIMA模型：\(ARMA(p,d,q)\times(P,D,Q) _s\)

模型识别

k阶滞后自相关函数（ACF）

回顾第一节，我们有k阶滞后的自相关函数:
$$
\rho _{k} = Corr(Y _t,Y _{t-k}) = \frac{Cov(Y _t, Y _{t-k})}{\sqrt{Var(Y _t)Var(Y _{t-k})}} = \frac{\gamma _{k}}{\gamma _{0}}
$$

$$
\gamma _{k} = Cov(Y _t, Y _{t-k}) = \frac{\sum _{t=k+1} ^{n} {(Y _t - \overline{Y})(Y _{t-k} - \overline{Y})} }{\sum _{t=1} ^n {(Y _t - \overline{Y}) ^2}} k=1,2,\cdots
$$
对于MA(q)模型而言当k>q时，\(\rho _k = 0\)，这里r为\(\rho \)的估计量。因此对于MA模型而言我们可以使用K阶之后自相关函数确定起q的取值。

k阶滞后偏自相关函数（PACF）

对于AR模型而言，由于模型往往可以展开成无穷阶MA模型，因此k阶滞后自相关函数不适用。从另一方面考虑，若我们希望计算\(Y _{t-k} \)对\(Y _t \)的影响程度从而判断模型是否含有第k阶，直接计算\(Y _{t-k} \)与\(Y _t \)的相关性是不合理的，这是因为对于\(Y _t \)而言其自身还包含了
\(Y _{t-1},Y _{t-2}, \cdots ,Y _{t-k+1} \)项。因此我们考虑消除\(Y _t \)中的
\(Y _{t-1},Y _{t-2}, \cdots ,Y _{t-k+1} \)影响后，再计算消除影响后的\(Y _t \)与\(Y _{t-k} \)的相关系数，我们称这个系数为k阶滞后偏自相关系数。
由于AR模型为回归方法，为了达到消除影响的目的，因此我们可以考虑先对序列进行拟合，然后为了计算得到样本的k阶滞后偏自相关函数。对\(Y _t \)我们使用\(Y _{t-1},Y _{t-2}, \cdots ,Y _{t-k+1} \)项进行预测：

$$
Y _t=
$$
对\(\beta \)我们可以使用均方误差最小化，利用最小二乘法估计。
对于\(Y _{t-k} \)，我们可以利用t-k之后的k-1个项来进行向前预测，由序列平稳可得（？？？）：
$$
Y _k=\beta _1 Y _{t-K+1} + \beta _2 Y _{t-K+2} + \cdots + \beta _{k-1} Y _{t-1}
$$
因此我们定义K阶偏子相关系数为预测误差之间的相关系数：
$$
\phi _{kk}=Corr(Y _t - \beta _1 Y _{t-1} - \beta _2 Y _{t-2} - \cdots - \beta _{k-1} Y _{t-k+1}, \\
Y _{t-k} - \beta _1 Y _{t-K+1} - \beta _2 Y _{t-K+2} - \cdots - \beta _{k-1} Y _{t-1} )
$$

对于若序列服从AR(p)模型，显然当预测的模型为p项时，对模型的预测将最准确，此时k阶之后偏子相关系数有：
$$
\phi _{kk} = 0, k>p
$$
而对于MA模型来说，其偏自相关函数从不为0，但随着滞后的增加，会指数级快速衰减到0。这一点与AR(1)过程的子相关函数类似。
对于K阶偏自相关系数，我们可以用最小二乘法分别求不同阶的预测模型，然后再求其相关系数。同时我们也可以使用Yule-Walker方程来进行求解。
具有自相关函数\(\rho _k\)的任意平稳过程，对于k阶滞后，我们有：
$$
\rho _j = \phi _{k1} \rho _{j-1} +\phi _{k2} \rho _{j-2} +\phi _{k3} \rho _{j-3} + \cdots+\phi _{kk} \rho _{j-k} , j=1,2,\cdots,k
$$
我们可以先算出各阶\(\rho\)然后带入方程组从而求解出\(\phi _{k1},\phi _{k2},\cdots,\phi _{kk} \)，这里我们只需要最后一项。对于不同阶k，反复求解方程组从而得到\(\phi _{kk}\)，我们有\(\phi _{kk}\)递归方程：
$$
\phi _{kk} = \frac{\phi _{k} - \sum _{j=1} ^{k-1}{\phi _{k-1,j}\rho _{k-j}}}{1 - \sum _{j=1} ^{k-1}{\phi _{k-1,j} \rho _{j}}}
$$
对于大于p时样本偏自相关函数近似服从均值为0，方差为1/n的正态分布. 因此，当 k>p 时，我们可以用\(\pm 2/\sqrt{n}\)作为\(\hat \phi _{kk}\)的临界极值来检验AR(p)模型的正确0假设。

非平稳性

在实际使用ARIMA模型时，我们往往仅使用1～2阶差分，过度的差分，会导致差分后的序列产生不必要的相关性，从而使模型变得复杂，最终导致过拟合。同时过度差分还可能导致差分后序列模型不可逆。
在对差分阶数进行判断时，ACF仍然是一个非常好用的工具。如果ACF的图像呈现明显近似线性的衰减，说明非平稳性影响了序列的自相关性。此时序列需要进行差分。更具体的方法还有Dickey-Fuller单位根检验（检验AR方程是否存在单位根）等。

ARIMA模型的判别

对于ARMA模型来说，由于ACF和PACF都只能针对AR或MA其中一个模型进行判别，若二者混合，不能很好判别。因此我们引入了一些扩展的方法来识别。

EACF（扩展自相关法）：

如果AR模型是已知的，则可以从ARMA模型中去掉自回归部分从而得到一个纯MA过程。

AIC（Akaike’s Information Criterion,赤池信息准则）：

AIC要求模型最小化：
$$
AIC = -2\log {(D(p,q _\theta))+2k}
$$
如果模型包含常数项或者截距项，则k=p+q+1，否则k=p+q。
AIC是估计模型与真实模型等平均KL偏离度的估计量。另\(p(y _1,y _2,\cdots,y _n)\)是\(Y _1,Y _2,\cdots,Y _n\)的真实pdf，\(q _\theta(y _1,y _2,\cdots,y _n)\) 是参数为 \(\theta\)的模型相应的pdf。则q和p的KL偏离由下面公式定义：
$$
D(p,q _\theta) = \int _{-\infty} ^{\infty} \int _{-\infty} ^{\infty} \cdots\int _{-\infty} ^{\infty}{p(y _1,y _2,\cdots,y _n) \log{[\frac{p(y _1,y _2,\cdots,y _n)}{q _\theta(y _1,y _2,\cdots,y _n)}]}}dy _1dy _2\cdots dy _n
$$
AIC估计了\(E(D(p,q _\hat \theta))\)，这里\(\hat \theta\)是参数\(\theta\)的极大似然估计。但是AIC是有偏估计量，若参数数量和数据容量的比率较大，则偏差会很大。因此引入\(AIC _c\)方法对AIC进行修正：
$$
AIC _c = AIC + \frac{2(k+1)(k+2)}{n-k-2}
$$
其中n是样本容量，k是出去噪声方差后的总参数量，当k/n>0.1时，\(AIC _c\)表现优于其他模型选择准则（包括AIC、BIC）

BIC（Bayesian Information Criterion,贝叶斯信息准则）：

$$
BIC = -2\log {(D(p,q _\theta))+k\log{(n)}}
$$

参数估计

当确定相应的ARIMA阶数后，我们既可以进行参数估计。常用的参数估计有最小二乘法与极大似然函数。

最小二乘法

对于AR模型而言，其满足标准最小二乘法形式，因此用最小二乘法估计十分简单。但对MA模型，由于其中存在e的噪声项，因此无法直接求的，所以我们需要对其进行转换：

$$
e _t = Y _t + \theta _1 e _{t-1}+ \theta _2 e _{t-2} + \cdots + \theta _{q} e _{t-q}
$$
其中\(e _0 = e _{-1}=\cdots=e _{-q} = 0\)，因此：

$$
e _1 = Y _1
$$

$$
e _2 = Y _2 + \theta _1 e _{t-1}
$$

$$
e _3 = Y _3 + \theta _1 e _{2}+ \theta _2 e _{1}
$$

$$
\cdots
$$
此时，我们最小化目标函数：

$$
S _c (\theta) = \sum{(e _t)^2}
$$

通过递推的计算\(e _1 , e _{2} , \cdots , e _{q}\)，带入方程组，利用多元数值算法，可以就\(\theta _1 ,\theta _{2} , \cdots , \theta _{q}\)联合的求取平方和最小值。
对于ARMA(p,q)模型，通过对模型公式进行整理，我们也可以使用同样方法：
$$
e _t = Y _t - \phi _1 Y _{t-1} - \phi _2 Y _{t-2} - \cdots - \phi _p Y _{t-p} + \theta _1 e _{t-1}+ \theta _2 e _{t-2} + \cdots + \theta _{q} e _{t-q}
$$
注意对于大样本而言，初始值\(e _{p} , e _{p-1},\cdots,e _{p-q+1}\)对于可逆模型影响甚。因此这里\(e _{p} = e _{p-1}=\cdots=e _{p-q+1} = 0\),序列从p+1开始递推。
####自回归逼近法
自回归逼近法，先将数据看做纯\(AR(\hat p)\)模型，利用AIC定阶确定\(\hat p\)，利用最小二乘法估计出相应的参数。用得到的\(AR(\hat p)\)模型对序列进行一步预测得到预测变量\(\hat Y _t\)，从而进一步估计出\(\hat e _t = Y _t - \hat Y _t\)。将得到的\(\hat e _t \)与真实数据\( Y _t\) 带入ARMA模型，于是就可以使用最小二乘法直接对AR与MA的参数进行估计。

极大似然估计

对于最小二乘法，我们利用了分布的一阶矩和二阶矩，而对于最大似然法，我们使用了数据的全部信息。
模型定阶完成后，对于ARIMA(p,d,q)实际需要估计的参数有\(\phi,\theta,\mu,\sigma _e\)，对于序列中的\(e _t\)服从于\( Normal(0,\sigma ^2)\)，且彼此独立。因此我们可以得到\(e _2,e _3,\cdots,e _n\)联合概率密度：
$$
(2\pi\sigma _e ^2) ^{-(n-1)/2}exp(-e _t ^2 /2\sigma _e ^2\sum _{t=2} ^{n} e _t ^2)
$$
ARIMA的极大似然估计法通过逐步递推预测公式导出，极大似然函数，具体过程较为复杂，可以参考PKU李东风老师的备课笔记CH21[^PKUTimeSerious]

模型预测

在获得模型之后，既可以用模型对未来数据进行预测，对于预测结果，我们同时还可以估计预测的精度。
我们令：\(\hat Y _t(l)\)为序列前t项已知状态下对t之后第l项进行预测的预测值。

AR模型的预测

则对于AR模型而言，有：
$$
\hat Y _t(l) =\sum _{i=1} ^{p} {\phi _i Y _{t+l-i}}\\
Y _{t+l-i} = \hat Y _t(l-i) \qquad when:\quad l-i>0
$$
可以采用一步预测的方式，从预测\(\hat Y _t(1)\)开始逐步递推到l项。
可以看到与AR模型相比，预测缺少了\(e _{t+l}\) 项。对于一步预测\(\hat Y _t(1)\)而言,预测误差\(e _t(1) = Y _{t+l} - \hat Y _t(l) = e _{t+1}\)。
对于
$$
\hat e _t(l) =\sum _{i=1} ^{\infty} {\psi _i e _{t+l-i}}\\
e _{t+l-i} = 0 \qquad when:\quad l-i<=0
$$
其中，\(\psi\)是由AR模型进行MA展开得到的相应参数。
因此，\(E(e _t(l)) = 0\)，相应预测估计是无偏估计。对于误差方差：
$$
Var(e _t(l)) = \delta _e ^ 2 (1 + \psi _{1} ^2 + \psi _{2} ^2 + \cdots + \psi _{l-1} ^2 )
$$
由于我们假设序列平稳，因此这里：
$$
Var(e _t(l)) \approx Var(Y _t) = \gamma _0 , 当l较大时
$$
证明略，感兴趣的同学，可以尝试使用AR(1)模型进行验证。
同时，对于arima模型，若序列非平稳则方差随预测时间窗口l的增大而无限增大。

MA模型的预测

对于MA模型而言，由于预测模型与e相关，因此需要先对e进行求解。
我们使用模型估计时用到的递推估计方法，求得\(E(e _{t+l-1}|Y _1,Y _2, \cdots, Y _t)\)。因此，可以得到：

$$
\hat Y _t(l) =\sum _{i=1} ^{p} {\phi _i E(e _{t+l-i} | Y _1,Y _2, \cdots, Y _t)}\\
E(e _{t+l-i} | Y _1,Y _2, \cdots, Y _t) = \begin{cases}
0, \qquad when\quad l-i>0\\
e _{t+l-i},\qquad when \quad l-i<=0
\end{cases}
$$

其中对于l-i>0的情况，由于相应的Y值是我们预测得到，因此我们无法通过递推关系求得相应残差的期望，另一方面，我们希望我们的预测尽可能准确，因此另残差期望为零也是合理的。
这里值得注意的一点，我们使用递推估计的方法逼近残差，需要满足模型可逆的前提。
同时，在对\(E(e _{t+l-1}|Y _1,Y _2, \cdots, Y _t)\)进行估计的时候，因为我们将e序列转化成了无穷序列的Y只和，因此，我们可以考虑，当逼近项数j>t-q时，对应Y的参数项\(\pi _j\)可以忽略。

ARIMA模型的预测

对arima模型的预测可分为AR部分与MA部分，值得注意的是，若MA的截断j值选取t-q时，当l>q时，ARIMA模型中MA部分将不起作用，因此此时模型只包含AR部分。
对ARMA模型，以\(\{e _t\}\)为基进行展开，得到相应参数\(\{\psi _t\}\)，则此时预测对应残差为：
$$
e _t(l) = Y _{t+1} - \hat Y(l) \\
= e _{t+l} + \psi _{1} e _{t+l-1} + \psi _{2} e _{t+l-2} + \cdots + \psi _{l-1} e _{t+1}
$$
也就是说对于真实的\(Y _{t+l}\)，对l周期内的误差对l点的预测具有累计效应。当序列平稳时，l趋于无穷，序列误差对方差趋于\(\gamma _0 \)。
这里由于我们误差服从正态分布，因此通过计算置信区间，利用误差的方差，我们可以对预测结果的范围进行估计。

模型诊断

为判断学习出的ARIMA模型是否有效，我们需要对得到的模型进行诊断。

残差分析

对于ARIMA模型，我们假设我们的数据满足模型:
$$Y _t=\pi _1 Y _{t-1} + \pi _2 Y _{t-2} + \cdots + e _{t-q}$$
这里假设模型平稳可逆，\(\pi\)是模型差分和MA部分展开并与AR部分合并后得到的合并参数。
在模型建立好后，对每一个\(\pi\)，我们有一个估计量\(\hat \pi\)，因此使用过去序列对现在时间点t进行预测后，我们可以得到残差（实际值-预测值）：
$$\hat e _{t-q} = Y _t - \hat \pi _1 Y _{t-1} - \hat \pi _2 Y _{t-2} + \cdots $$
如果我们的ARIMA模型预测无误，则\(\hat e _{t} = e _{t}\)，也就是说我们的残差应：1. 服从正态分布 2. 每一时刻的残差与前序时刻不相关。通过判断残差这两点特性，我们可以分析出ARIMA模型对数据的拟合程度。
在实际判别中，我们主要判断方法有：残差图，QQ图，残差自相关性判定，统计量检验等方法。

残差图

由于残差应近似服从正态分布，因此我们绘制出残差的时间序列散点图（即残差图，X轴为时间，Y轴为残差值），通过比较散点图的随机性，与正态分布的特点来判断其正态性。这里我们可以关注：

1. 在相邻时刻残差是否产生了明显的相关特点.
2. 不同时间段内，残差的变化幅度是否相近.
3. 残差值超过\( 1\sigma, 2\sigma, 3\sigma \)的数量.

如图：
$$此处应有图+解释$$

残差正态性与QQ图（分位数-分位数图）

使用残差图我们可以很好的看出残差在不同时间段是否有特殊的特征，但其对正态分布的实际分布情况并不能完全的刻画，因此我们引入了QQ图（分位数-分位数图）。
QQ图横坐标是理论分位数，纵坐标是实际分位数，我们通过算出每一个时刻t，残差在实际残差序列中的分位数位置作为Y轴，其在模型中残差的理论上的高斯分布应处于的分位数位置作为X轴，从而绘制出散点图。如果残差满足正态分布，则其应满足：

1. 其散点应该分布在y=x这一条直线上
2. 散点越靠近原点(0,0)点越密集，越远离原点越稀疏。
3. 稀疏变化程度应该与正态分布的分布一致

如图：
$$此处应有图+解释$$

残差自相关性

除了残差的正态性，我们还可以从残差相关性入手。通过计算残差序列不同阶的相关性，我们可以画出残差的ACF图。若无相关性，则各阶滞后均应该处于一个较低的水平。通过对相关性显著性分析，我们可以定量的判断模型是否有异常。
$$ACF显著性判断$$
如图：
$$此处应有图+解释$$

通过观察ACF图，我们可以找到是否有单个的滞后相关性具有异常。而另一方面，我们还可以对各个滞后进行综合考虑，可能单个相关系数没有异常，但是整体来看，具有异常。因此我们引入Ljung-Box检验，其统计量为：
$$Q=n(\hat r _1 ^2 + \hat r _2 ^2 + \cdots +\hat r _k ^2 )$$
若我们估计的模型正确，则Q应该近似服从自由度为k-p-q的卡方分布，若假设检验不能通过，则可以判断模型拟合有误。5t t

过拟合和参数冗余

数据缺失http://rpubs.com/englianhu/handle-missing-value

谱表示与谱分析

ARIMA模型主要从滑动窗口和序列自相关性两方面对序列的周期性进行拟合。而考虑到对周期性问题进行拟合，我们很容易想到信号处理中的谱表示，即通过使用不同频率的正弦余弦函数，拟合出我们希望得到的周期变化。同时，我们还引入了傅立叶变换，将序列从空域变换到频域，通过分析频域图像，从而对周期性，相关性等进行更好的分析，这是谱分析的基本思想。

基本概念

基本公式与概念

一下为一些简单常用的公式：
$$
R\ cos(2\pi ft + \Phi)=A\ cos(2\pi ft )+B\ sin(2\pi ft)
$$
$$
R = \sqrt{(A^2 + B^2)},\ \Phi = artan(-B/A)
$$
$$
B = R * cos(\Phi),\ A = -R * sin(\Phi)
$$

对于谱表示来说，我们可以将通项公式写成如下形式：
$$
Y _t = A _0 + \sum _{j=1} ^m {[A _j cos(2 \pi f _j t) + B _jsin(2 \pi f _j t)]}
\tag{谱1.1}
$$

这是一种利用正弦余弦对，对时间序列的拟合，对于其中的各项A与B我们可以用最小二乘法进行回归，对于其中频率f 如果我们取特殊形式，则回归将会变得特别简单，例如：
1）n=2k+1：
则频率为依次为：1/n,2/n,…,k/n (k/n=1/2-1/(2n))。这种频率我们称之为傅立叶频率，在这些频率下，正弦和余弦的预测变量是正交的。
因此，最小二乘估计可得：

$$
\hat A _0 = \overline Y
$$

$$
\hat A _j = \frac{2}{n} \sum _{t=1} ^n {Y _t cos(2 \pi t j/n)},\ \hat B _j = \frac{2}{n} \sum _{t=1} ^n {Y _t sin(2 \pi t j/n)}
$$

2）n=2k：
则：
$$
\hat A _k = \frac{1}{n} \sum _{t=1} ^n (-1) ^t Y _t,\ \hat B _j = 0
$$

对于任意长度的序列，无论其是确定性的还是随机的，都可以用通项公式进行完美拟合。
同时，对于长序列，通常来说可以使用快速傅立叶变换进行快速求解。

周期图

通过上述通项公式，我们可以对时间序列进行拟合。但实际拟合中，我们需要确定实际需要拟合哪一些周期。这时周期图是一个较好的工具。
周期图定义为:
$$
I(j/n) = \frac{n}{2} (\hat A _j ^2 + \hat B _j ^2 )
$$
对于样本量为偶数时，即 \(n=2k)\)时，在极端频率\(f=k/b=1/2)\)时，有：

$$
I(1/2) = n(\hat A_k)^2
$$

周期图的高峰显示了序列整体行为中不同频率上余弦-正弦对的相对强度。我们可以选取周期图中峰值较高的周期进行建模。
实际上对于一般的时间序列而言，周期图是时间序列从空域向频域变换后频域的图像。
对于较长序列，我们可以使用快速傅立叶变换的方式（FFT）进行快速有效的数值计算。

谱表分析

谱表示

对于谱分析的通项公式（谱1.1）来说，我们先考虑\(A_0=0\)的情况，又或者说，我们可以移除样本均值，令\(Y _t\)为样本到样本均值的偏差。此时，我们仅考虑\(0<f<1/2\)的情况（ 对于f在(0,1/2)和(1/2,1)区间，图像实际上是重叠的 ）。
我们令：
$$
a(f)=\sum _{\{j|f _j<f\}} A _j, b(f)=\sum _{\{j|f _j<f\}} B _j
$$
从而构造出两个阶梯函数。从而我们可以对通项公式进行改写,得到谱表示：
$$
Y _t = \int _0 ^{1/2} cos(2\pi ft) \rm{d}a(f) + \int _0 ^{1/2} sin(2\pi ft)
\rm{d}b(f) \tag{谱1.2}
$$
事实上，任何零均值平稳的过程都能表示成上述方程的形式。即，平稳过程可以用连续频带上无穷多个余弦正弦对的线形组合表示。

周期图与谱密度

正如上文提到的，要分析时间序列，周期图是一个很好的工具。实际上周期图是对谱表示的时间序列过程进行傅立叶变换的频域图像。这里我们引入谱密度（功率谱）概念，一个信号的功率谱密度就是该信号自相关函数的傅里叶变换。谱密度是过程经过傅立叶变换后在频域的密度函数。

平稳随机过程的功率谱是一个确定函数。通过计算我们可以得到其公式为：

$$
S(f)=\gamma_0 + 2\sum _{k=1} ^\infty \gamma _k cos(2\pi fk)
$$

其中，\(\gamma _k\)为函数的k阶滞后样本协方差。根据傅立叶理论（QA：Planchere定理？是否需要过程满足先决条件，例如函数属于L2空间?），必存在一个相反关系，即：
$$
\gamma _k = \int _{-1/2} ^{1/2} S(f) cos(2\pi kf) \rm{d}f
$$
从数学上来讲，S(f)是序列 \(\{\gamma _k\}\)的离散时间傅立叶变换，而 \(\{\gamma _k\}\)是谱密度的傅立叶逆变换。
同时还有：
$$
F(f) = \int _{0} ^{f} S(x) \rm{d}x ,0<= f <=1/2
$$
可以看出，谱密度下方区域面积的两倍是过程方差对应频率区间上构成该过程的余弦正弦对的那部分。
对于周期图I(f)，我们有：
$$
\hat{S}(f)=1/2I(f),\hat{S}(1/2)=I(1/2),
$$

###附：
傅立叶变换公式：
$$
F(\omega) = \int _{-\infty} ^{\infty} f(t)e^{-i\omega t} \rm{d}t
$$

综合方法

Prophet

prophet是facebook开源的python预测库，该库的api设计与sklearn很像，也是分为fit方法和predict方法。prophet将时间序列分解为三个不同的成分：趋势成分，周期性成分与特殊点（如节假日），并采用加法模型，将三个成分相结合：
$$
y(t)=g(t)+s(t)+h(t)+\epsilon_t
$$
其中g(t)为趋势函数，代表序列中非周期性的改变，s(t)为周期性成分，h(t)为节假日等特殊点。\(\epsilon _t\) 为随机成分，我们一般令其为正态分布。
我们知道对于ARIMA模型而言我们可以利用差分方法来拟合趋势，而利用这种加性模型，对于趋势我们可以更加自由的进行拟合。

趋势拟合

prophet主要使用了两种趋势拟合模型，增长模型与线性模型

增长模型

增长模型常被用于生态系统模型中（被广为熟知的增长s曲线）如下：

$$
g(t)=\frac{C}{1+exp(-k(t-m))}
$$

然而在实际使用过程中，这种完美的曲线是不存在的，以用户数据为例：当新功能发布，用户增长率（k）与曲线容量（C）可能发生变化，而实际情况中，这种变化是经常发生的。因此，我们需要对模型进行适应。
对于曲线容量（C）而言，我们可以使用一个随时间变化的动态容量C(t)。
而对变化率（k）来说，我们可以假定，变化率是在某些时刻发生突然变化的（这也符合我们上述例子中的实际情况）。我们将这些变化的点用序列\(\boldsymbol{S}=\{s _j\}\)表示，\(\boldsymbol{\delta}=\{\delta _j\}\)表示第\(s _j\)时间点上的变化率增量，令k为基础增长率，令向量\(\boldsymbol{a}(t)=\{0,1\} ^{|S|}\)，有

$$
a _j (t)=\begin{cases}
1,\quad t \geq s _j \\
0,\quad other
\end{cases}
$$

则对增长率，我们有\(k + a _j(t)^T\delta\)，而当增长率发生变化的时候，我们的偏移参数m也需要发生变化，举发生一次变化的图像为例：

f曲线是变化前曲线，在s点处，增长率发生变化，变化后曲线为f‘。则我们有：

$$
f(t)=\frac{C}{1+exp(-k(t-m))}
$$

$$
f(t)=\frac{C}{1+exp(-(k+\delta)(t-(m + \gamma)))}
$$
其中\(\delta,\gamma\)分别为k和m的增量，又由于曲线在t=s点相交，我们有f(s)=f’(s)。通过整理等式我们可得：

$$
\gamma=\frac{\delta(s-m)}{k+\delta}
$$

推广到多点的情况，则我们的得到：

$$
\gamma _j=\frac{(s _j -m -\sum _{l<j}\gamma _l)\delta _j}{k+\sum _{l \leq j}\gamma _l}
$$

令\(\boldsymbol{\gamma}=\{\gamma _j\}\)，则，我们得到最终模型：

$$
g(t)=\frac{C(t)}{1+exp(-(k+\boldsymbol{a}(t)^T\delta)(t-(m+\boldsymbol{a}(t)^t\boldsymbol{\gamma})))}
$$