【时间序列模型】基本性质判断及相关计算---附代码

---

前言

本文主要简单整理几种线性时间序列模型的表达式，时间序列过程性质的判断（如判断稳定性等）和模型参数的选择，并提供一些手动计算部分模型ACF和PACF的方法。文章最后附有简单的R相关代码。还在学习中，如有错误，欢迎指正。

---

一、常用时间序列模型

Prerequisite
时间序列： { x t } \left\{ x_t \right\} {xt​}
白噪声(white noise)：$w_t$ $$w_t\sim wn(0,{\sigma }_w^2)$$ 后移算子(Backshift operator): $B$ $$\begin{gathered} x_{t-1}=B{x}_t \\ {x}_{t-p}=B^p{x}_t \end{gathered}$$ 差分算子(Difference operator): $\nabla$ $$\nabla x_t=x_t-x_{t-1}=(1-B)x_t$$ d阶差分(Differenc of order d)$${\nabla }^d=(1-B{)}^d$$

1.自回归模型 —AR($p$)

$$x_t={\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+...+{\varphi }_p{x}_{t-p}+w_t，w_t\sim wn(0,{\sigma }_w^2)$$ $$\varphi (B)=1-{\varphi }_{1}B-{\varphi }_2{B}^2-...-{\varphi }_p{B}^p$$ $$\varphi (B)x_t=w_t$$

1.1 参数p的确定

①若已知时间序列过程表达式 $x_t$，则根据表达式 $x_t$中p的数值确定。
②若有一组时间序列数据，则画出其PACF图，PACF在滞后数为p的地方截尾，即PACF在 lag=p 处之后在2个标准差内。

1.2 平稳性判断（stationary）

将 $\varphi (B)$写成 $\varphi (z)$,则有$$\varphi (z)=1-{\varphi }_{1}z-{\varphi }_2{z}^2-...-{\varphi }_p{z}^p$$ 若具有平稳性，则 $\varphi (z)$的根不在单位圆上
令 $\varphi (z)=0$，并求解。 $|z|\ne 1⟺$ 具有平稳性(stationary)

1.3 因果性判断（causal）

若具有因果性，则 $\varphi (z)$的根在单位圆外
令 $\varphi (z)=0$，并求解。 $|z|>1⟺$ 具有因果性(causal)

---

2. 移动平均模型 —MA(q)

$$x_t=w_t+{\theta }_1{w}_{t-1}+{\theta }_2{w}_{t-2}+...+{\theta }_q{w}_{t-q}，w_t\sim wn(0,{\sigma }_w^2)$$ $$\theta (B)=1+{\theta }_{1}B+{\theta }_2{B}^2+...+{\theta }_q{B}^q$$ $$x_t=\theta (B)w_t$$

2.1 参数q的确定

①若已知时间序列过程表达式 $x_t$，则根据表达式 $x_t$中q的数值确定。
②若有一组时间序列数据，则画出其ACF图，ACF在滞后数为p的地方截尾，即ACF在 lag=q 处之后在2个标准差内。

2.2 平稳性判断（stationary）

对于移动平均模型，总是平稳的。

2.3 可逆性判断（invertible）

将 $\theta (B)$写成 $\theta (z)$,则有$$\theta (z)=1+{\theta }_{1}z+{\theta }_2{z}^2+...+{\theta }_q{z}^q$$ 若具有可逆性性，则 $\theta (z)$的根在单位圆外
令 $\theta (z)=0$，并求解。 $|z|>1⟺$ 具有可逆性(invertible)

---

3. 自回归移动平均模型 —ARMA(p,q)

$$x_t-{\varphi }_1{x}_{t-1}-{\varphi }_2{x}_{t-2}-...-{\varphi }_p{x}_{t-p}=w_t+{\theta }_1{w}_{t-1}+{\theta }_2{w}_{t-2}+...+{\theta }_q{w}_{t-q}，w_t\sim wn(0,{\sigma }_w^2)$$ $$\varphi (B)x_t=\theta (B)w_t$$

3.1 p,q

①根据 $\varphi (B)$和 $\theta (B)$的阶数确定ARMA(p,q)
避免参数冗余（redundancy）
确定 $\varphi (B)$和 $\theta (B)$的阶数时，两个表达式不能有相同的因子，即两个特征方程不能有相同的根。$\varphi (B)$和 $\theta (B)$有相同因子则模型表达式存在参数冗余。
②根据ACF图和PACF图酌情选择

3.2 稳定性，因果性，可逆性

$$\varphi (z)=0⟹\{\begin{array}{c}|z|\ne 1⟺具有平稳性(stationary)\\ |z|>1⟺具有因果性(causal)\end{array}$$ $$\theta (z)=0⟹|z|>1⟺具有可逆性(invertible)$$

3.3 转换系数 $\psi ,\pi$

①对于具有因果性的模型，可计算 $\psi$ $$\begin{gathered} x_t=\sum _{j=0}^{\infty }{\psi }_j{w}_{t-j},\psi (B)={\psi }_0+{\psi }_{1}B+{\psi }_2{B}^2-... \\ \{\begin{array}{c}\varphi (B)x_t=\theta (B)w_t\\ x_t=\psi (B)w_t\end{array}⟹\varphi (B)\psi (B)=\theta (B) \end{gathered}$$ $$(1-{\varphi }_{1}B-{\varphi }_2{B}^2-...)({\psi }_0+{\psi }_{1}B+{\psi }_2{B}^2+...)=(1+{\theta }_{1}B+{\theta }_2{B}^2+...)$$ 根据不同阶B前面的系数，可以写出 $\varphi ,\psi$和 $\theta$ 的等式，进而求出系数 $\psi$。
②对于具有可逆性的模型，可计算 $\pi$ $$\begin{gathered} w_t=\sum _{j=0}^{\infty }{\pi }_j{x}_{t-j},\pi (B)={\pi }_0+{\pi }_{1}B+{\pi }_2{B}^2-... \\ \{\begin{array}{c}\varphi (B)x_t=\theta (B)w_t\\ w_t=\pi (B)x_t\end{array}⟹\varphi (B)=\pi (B)\theta (B) \end{gathered}$$ $$(1-{\varphi }_{1}B-{\varphi }_2{B}^2-...)=(1+{\theta }_{1}B+{\theta }_2{B}^2+...)({\pi }_0+{\pi }_{1}B+{\pi }_2{B}^2+...)$$ 根据不同阶B前面的系数，可以写出 $\varphi ,\pi$和 $\theta$ 的等式，进而求出系数 $\pi$

---

4. ARIMA(p,d,q)

时间序列不平稳，进行d阶差分后为平稳序列ARMA
$$\begin{gathered} {\nabla }^d{x}_t=(1-B{)}^d{x}_t \\ \varphi (B)(1-B{)}^d{x}_t=\theta (B)w_t \end{gathered}$$ 一阶差分------消除线性趋势
二阶差分------消除二次趋势
差分应适度，过少差分模型仍不平稳，过度差分会引入额外依赖。

---

5. 季节性模型

5.1 Pure seasonal ARMA

周期为s,记作$ARMA(P,Q{)}_s$ $$\begin{gathered} {\Phi }_P(B^s)x_t={\Theta }_Q(B^s)w_t \\ {\Phi }_P(B^s)=1-{\Phi }_1{B}^s-{\Phi }_2{B}^{2s}-...-{\Phi }_P{B}^{Ps} \\ {\Theta }_Q(B^s)=1+{\Theta }_1{B}^s+{\Theta }_2{B}^{2s}+...+{\Theta }_Q{B}^{Qs} \end{gathered}$$ 因果性----$\Phi (z^s)$ 根在单位圆外
可逆性----$\Theta (z^s)$ 根在单位圆外

5.2 SARIMA

记作$ARIMA(p,d,q)\times (P,D,Q{)}_s$ $$\begin{gathered} {\nabla }_s^D=(1-B^s{)}^D \\ {\Phi }_P(B^s)\varphi (B){\nabla }_s^D{\nabla }^d{x}_t={\Theta }_Q(B^s)\theta (B)w_t \end{gathered}$$ 因果性----$\Phi (z^s),\varphi (z)$ 根在单位圆外
可逆性----$\Theta (z^s),\theta (z)$ 根在单位圆外

---

---

二、相关系数计算

自协方差函数Autocovariance function of {$x_t$} 记作$\gamma (h)$
根据 $E(x_t)=0$
$\gamma (h)=\gamma (t+h,t)=Cov(x_{t+h},x_t)=E[(x_{t+h}-E(x_{t+h})(x_t-E(x_t)]=E(x_{t+h}{x}_t)$
对于白噪声$w_t$ $${\gamma }_w(h)=\{\begin{array}{c}{\sigma }_w^2,h=0\\ 0,otherwise\end{array}$$

1. 自相关ACF

自相关函数Autoccorelation function of {$x_t$} 记作$\rho (h)$
$$\rho (h)=\frac{\gamma (h)}{\gamma (0)}$$ 显然，有 $$\rho (0)=\frac{\gamma (0)}{\gamma (0)}=1$$

1.1 AR(1）

表达式：$$\begin{gathered} x_t=\varphi x_{t-1}+w_t,w_t\sim wn(0,{\sigma }_w^2) \\ {x}_t=\sum _{j=0}^{\infty }{\varphi }^j{w}_{t-j} \end{gathered}$$ 自协方差函数$\gamma (h)$及ACF$\rho (h)$ $$\begin{gathered} \gamma (h)=\frac{{\sigma }_w^2{\varphi }^h}{1-{\varphi }^2} \\ \rho (h)={\varphi }^h \end{gathered}$$

1.2 AR(2)

表达式：$$x_t={\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+w_t,w_t\sim wn(0,{\sigma }_w^2)$$ 对两边同乘$x_{t-h}$并取期望 $$\begin{gathered} E(x_t{x}_{t-h})={\varphi }_{1}E(x_{t-1}{x}_{t-h})+{\varphi }_{2}E(x_{t-2}{x}_{t-h}) \\ \gamma (h)={\varphi }_1\gamma (h-1)+{\varphi }_2\gamma (h-2)，h=1,2... \\ \rho (h)={\varphi }_1\rho (h-1)+{\varphi }_2\rho (h-2)，h=1,2... \end{gathered}$$ 可以推出：$$\begin{gathered} \rho (1)={\varphi }_1\rho (0)+{\varphi }_2\rho (1) \\ \rho (1)=\frac{{\varphi }_1\rho (0)}{1-{\varphi }_2}=\frac{{\varphi }_1}{1-{\varphi }_2} \end{gathered}$$ 由齐次差分方程，得通解 ：$$\begin{gathered} \varphi (z)=0⟹z_1,z_2 \\ \rho (h)=c_1(z_1{)}^{-h}+c_2(z_2{)}^{-h} \end{gathered}$$根据得到的$\rho (0),\rho (1)$ 代入通解，求得 $c_1,c_2$，最终得到$\rho (h)$

1.3 MA(1)

表达式：$$\begin{gathered} x_t=w_t+\theta w_{t-1},w_t\sim wn(0,{\sigma }_w^2) \\ {w}_t=\sum _{j=0}^{\infty }(-\theta {)}^j{x}_{t-j} \end{gathered}$$ 自协方差函数$\gamma (h)$：$$\gamma (h)=\{\begin{array}{c}(1+{\theta }^2){\sigma }_w^2，h=0\\ \theta {\sigma }_w^2，h=1\\ 0，h>1\end{array}$$ 计算ACF $\rho (h)$ $$\rho (h)=\{\begin{array}{c}\frac{\theta }{1+{\theta }^2}，h=1\\ 0，h>1\end{array}$$

1.4 MA(q)

自协方差函数：$$\gamma (h)=\{\begin{array}{c}{\sigma }_w^2{\sum }_{j=0}^{q-h}{\theta }_j{\theta }_{j+h}，0\le h\le q\\ 0，h>q\end{array}$$ 计算ACF $\rho (h)$ $$\rho (h)=\{\begin{array}{c}\frac{{\sum }_{j=0}^{q-h}{\theta }_j{\theta }_{j+h}}{1+{\theta }_1^2+{\theta }_2^2+...+{\theta }_q^2}，1\le h\le q\\ \\ 0，h>q\end{array}$$

1.5 ARMA(p,q)

与求AR($p$)的ACF思路相同，将ARMA(p,q)转为AR($p$)的形式
已知系数$\psi ，x_t={\sum }_{j=0}^{\infty }{\psi }_j{w}_{t-j}$
先由p阶齐次方程得到自相关函数通解：$$\gamma (h)-{\varphi }_1\gamma (h-1)-...-{\varphi }_p\gamma (h-p)=0，h\ge max(p,q+1)$$ 再根据初始条件 $$\gamma (h)-\sum _{j=1}^p{\varphi }_j\gamma (h-j)={\sigma }_w^2\sum _{j=h}^q{\theta }_j{\psi }_{j-h}，0\le h<max(p,q+1)$$ 将得到的自相关函数等式分别除以$\gamma (0)$，最终得出$\rho (h)$

---

2. 偏自相关PACF

平稳过程的偏自相关函数 Partial Autoccorelation function of {$x_t$} 记作${\varphi }_{hh}$
$\widehat{x_{t+h}}$和$\widehat{x_t}$分别是{$x_{t+1},x_{t+2},...,x_{t+h-1}$} 对$x_{t+h}$和$x_t$的线性回归
${\varphi }_{hh}$是$x_{t+h}$和$x_t$在移去{$x_{t+1},x_{t+2},...,x_{t+h-1}$}的影响后的相关系数
$$\begin{gathered} {\varphi }_{11}=corr(x_{t+1},x_t)=\rho (1) \\ {\varphi }_{hh}=corr(x_{t+h}-\widehat{x_{t+h}},x_t-\widehat{x_t})，h\ge 2 \end{gathered}$$

2.1 AR(1)

$$\begin{gathered} {\varphi }_{11}=1 \\ {\varphi }_{hh}=0,h>1 \end{gathered}$$

2.2 AR(2)

$$\begin{gathered} {\varphi }_{11}=1 \\ {\varphi }_{22}={\varphi }_2 \\ {\varphi }_{hh}=0,h>2 \end{gathered}$$

2.3 MA(1)

$${\varphi }_{hh}=-\frac{(-\theta {)}^h(1-{\theta }^2)}{1-{\theta }^{2(h+1)}}，h\ge 1$$

2.5 Durbin-Levinson算法

$${\varphi }_{hh}=\frac{\rho (h)-{\sum }_{k=1}^{h-1}{\varphi }_{h-1,k}\rho (h-k)}{1-{\sum }_{k=1}^{h-1}{\varphi }_{h-1,k}\rho (k)}$$ 其中，当$h\ge 2$时，${\varphi }_{hk}={\varphi }_{h-1,k}-{\varphi }_{hh}{\varphi }_{h-1,h-k}$
对于AR($p$), ${\varphi }_{pp}={\varphi }_p$，且 $h>p$ 时，${\varphi }_{hh}=0$

2.4 Yule-Walker Equation

$\left[\begin{array}{c}\rho (1)\\ \rho (2)\\ \rho (3)\\ .\\ .\\ .\\ \rho (h)\end{array}\right]$ = $\left[\begin{array}{cccccc}1& \rho (1)& \rho (2)& .& .& \rho (h-1)\\ \rho (1)& 1& \rho (1)& .& .& \rho (h-2)\\ \rho (2)& \rho (1)& 1& .& .& \rho (h-3)\\ .& .& .& & & .\\ .& .& .& & & .\\ .& .& .& & & .\\ \rho (h-1)& \rho (h-2)& \rho (h-3)& .& .& 1\end{array}\right]$ $\left[\begin{array}{c}{\varphi }_{h,1}\\ {\varphi }_{h,2}\\ {\varphi }_{h,3}\\ .\\ .\\ .\\ {\varphi }_{h,h}\end{array}\right]$
循环计算，每一组方程组的 ${\varphi }_{h,h}$即为 ${\varphi }_{hh}$

---

3. 小结

	AR($p$)	MA(q)	ARMA(p,q)
ACF	拖尾	在lag=q之后截尾	拖尾
PACF	在lag=p之后截尾	拖尾	拖尾

SARMA	AR($P$)${}_s$	MA(Q)${}_s$	ARMA(P,Q)${}_s$
ACF	在lag=ks 处拖尾 (k=1,2,…)	在lag=Qs处截尾	在lag=ks处拖尾
PACF	在lag=Ps处截尾	在lag=ks处拖尾(k=1,2,…)	在lag=ks处拖尾

在非季节性滞后点（h$\ne$ks）处，值为0

---

---

三. R 代码

## E.g. AR(2) phi1=0.2,phi2=0.8
#生成100个模拟数
ar2 = arima.sim(list(order=c(2,0,0), ar=c(-.2,-.8)), n = 100)
#计算该模型的acf
ACF_AR2 = ARMAacf(ar=c(-.2,-.8), ma=0, 50)
#画出acf图
plot(ACF_AR2, type="h", xlab="lag")
abline(h=0)
title("AR(2)")
#计算该模型的pacf
PACF_AR2 = ARMAacf(ar=c(-.2,-.8), ma=0, 10,pacf=TRUE)
#画出pacf图
plot(PACF_AR2, type="h", xlab="lag")
title("AR(2)")
abline(h=0)

#使用 R包中的时间序列数据
library(astsa)
plot.ts(globtemp)  #画出时间序列，观察趋势
acf(globtemp) #画出acf
pacf(globtemp) #画出pacf
#进行一阶差分
globtemp_d <- diff(globtemp)
plot.ts(globtemp_d)
#画acf和pacf
par(mfrow = c(2,1))
acf(globtemp_d)
pacf(globtemp_d)
#拟合ARIMA(0,1,2)
sarima(globtemp, 0, 1, 2)
#预测之后10个序列值
sarima.for(globtemp, 10, 0, 1, 2) #输出预测值，SE和预测图

#若为季节性时间序列模型
#使用 UnempRate 数据(月度数据)，选择模型为(3,1,2)×(0,1,1),s=12
sarima(UnempRate,3,1,2,0,1,1,12) #拟合模型
sarima.for(UnempRate,12, 3,1,2,0,1,1,12)  #预测下一年(12个月)

---

---

参考书目

Time Series Analysis and Its Applications with R examples
Time Series Analysis With Applications in R