时间序列分析:第一章

时间序列的特征

学习笔记
统计方法
时间序列分析
作者

周博霖

发布于

September 6, 2024

重要

第五版年底出版,目前暂停更新

主要参考书籍为Shumway和Stoffer(2017)的Time Series Analysis and Its Applications: With R Examples。

该书第五版正在制作,配套astsa包中的部分数据集会逐渐更新,所以后期也许会换成第五版内容。

该内容为学习笔记,会很碎片化,同时参考了作者在github上给出的代码。

提示

提示框表明这部分完全是个人理解,可能有误。

一切错误和想讨论的内容欢迎联系邮箱(zhoubolin0404@126.com),感谢!

时间序列分析适用情境:

两种不同但不一定互斥的时间序列分析方法:

1 时间序列数据

1.1 强生公司季度每股收益

jj:强生公司的季度每股收益,从1960年第一季度到1980年最后一季度的84个季度(共21年)。

library(astsa)
jj
          Qtr1      Qtr2      Qtr3      Qtr4
1960  0.710000  0.630000  0.850000  0.440000
1961  0.610000  0.690000  0.920000  0.550000
1962  0.720000  0.770000  0.920000  0.600000
1963  0.830000  0.800000  1.000000  0.770000
1964  0.920000  1.000000  1.240000  1.000000
1965  1.160000  1.300000  1.450000  1.250000
1966  1.260000  1.380000  1.860000  1.560000
1967  1.530000  1.590000  1.830000  1.860000
1968  1.530000  2.070000  2.340000  2.250000
1969  2.160000  2.430000  2.700000  2.250000
1970  2.790000  3.420000  3.690000  3.600000
1971  3.600000  4.320000  4.320000  4.050000
1972  4.860000  5.040000  5.040000  4.410000
1973  5.580000  5.850000  6.570000  5.310000
1974  6.030000  6.390000  6.930000  5.850000
1975  6.930000  7.740000  7.830000  6.120000
1976  7.740000  8.910000  8.280000  6.840000
1977  9.540000 10.260000  9.540000  8.729999
1978 11.880000 12.060000 12.150000  8.910000
1979 14.040000 12.960000 14.850000  9.990000
1980 16.200000 14.670000 16.020000 11.610000
tsplot(jj, col = 4, type = "o", xlab = "时间", ylab = "季度每股收益")
图 1: 强生公司季度每股收益

1.2 全球变暖

这里我使用了最新数据,老数据globtemp计划被删除。

gtemp_land:1850年到2023年,每年地球陆地年平均温度偏差(以摄氏度为单位,相对于1991-2020年的平均值)。

gtemp_land
Time Series:
Start = 1850 
End = 2023 
Frequency = 1 
  [1] -0.50 -0.60 -0.50 -0.50 -0.20 -0.50 -0.80 -0.40 -0.40 -0.10 -0.70 -0.20
 [13] -0.50 -0.20 -0.30 -0.60 -0.40 -0.60 -0.30 -0.50 -0.20 -0.10 -0.40 -0.30
 [25] -0.40 -0.50 -0.20 -0.10  0.26 -0.20 -0.50 -0.40 -0.20 -0.40 -0.70 -0.40
 [37] -0.50 -0.30 -0.70  0.00 -0.20 -0.30 -0.50 -0.40  0.10 -0.50 -0.50 -0.70
 [49] -1.10 -0.50  0.00  0.21 -0.40  0.05 -0.60 -0.70 -0.30 -0.30 -0.80 -0.70
 [61] -0.50 -0.80 -0.70 -0.40 -0.10 -0.20 -0.50 -1.00 -0.30 -0.50  0.05 -0.10
 [73]  0.00 -0.30  0.01 -0.20  0.53 -0.60 -0.20 -0.30  0.21 -0.10 -0.10 -0.30
 [85] -0.30  0.13 -0.30 -0.40  0.38 -0.20  0.06 -0.10 -0.10  0.00  0.36  0.00
 [97]  0.00  0.32 -0.30  0.24  0.00 -0.20 -0.30  0.33 -0.30 -0.60 -0.20 -0.30
[109]  0.03  0.49 -0.90  0.08  0.36 -0.40 -0.30 -0.10 -0.10  0.19  0.80 -0.20
[121]  0.00  0.00  0.13  0.61  0.18  0.32 -0.20  0.29  0.37  0.10  0.19  0.85
[133]  0.00  0.64  0.36  0.42  0.67 -0.20  0.51  0.66  1.61  0.45  0.78  0.64
[145]  0.41  0.75  0.44  0.81  0.94  0.36  1.14  0.85  1.58  0.84  1.02  1.36
[157]  1.20  1.27  1.65  0.83  1.54  1.27  0.93  1.10  1.52  1.76  2.50  2.15
[169]  1.63  2.14  1.95  1.58  2.13  2.26
tsplot(gtemp_land, col = 4, type = "o", xlab = "时间", ylab = "全球温度偏差/℃")
图 2: 全球陆地年平均温度偏差(1850-2023)

1.3 语音数据

speech:为“aaa…hhh”这一短语录制的0.1秒(1000采样点)语音样本。

attributes(speech)
$tsp
[1]    1 1020    1

$class
[1] "ts"
tsplot(speech)
图 3: 10000点每秒速度采样的音节aaa···hhh的语音记录,n = 1020

1.4 道琼斯工业平均指数

djia:从2006年4月20日到2016年4月20日,道琼斯工业平均指数的每日数值

djia数据目前只使用其中的收盘价(Close,不确定是昨收还是今收,不过这也不重要)。

假设\(x_t\)是实际值,则\(r_t=\frac{x_t-x_{t-1}}{x_{t-1}}\)是收益率,那么\(1+r_t=\frac{x_t}{x_{t-1}}\),则有\(log(1+r_t)=log(\frac{x_t}{x_{t-1}})=log(x_t)-log(x_{t-1})\approx r_t\)(\(r_t\)很小,泰勒级数展开后成立)。

djiar为2006年4月21日到2016年4月20日的日收益率。

library(xts)
djiar = diff(log(djia$Close))[-1]
plot(djiar, col = 4, main = "DJIA Returns") 
图 4: DJIA日收益率

1.5 厄尔尼诺和鱼群

soi:1950-1987年间453个月的南方涛动指数(SOI)

rec:1950-1987年间453个月的新鱼数量指数(Recruitment)

par(mfrow = c(2, 1)) 
tsplot(soi, col = 4, xlab = "", ylab = "", main = "Southern Oscillation Index")
tsplot(rec, col = 4, xlab = "", ylab = "", main = "Recruitment") 
图 5: 每月SOI和新鱼数量指数

1.6 fMRI成像

fmri1:由8个位置的平均fMRI 血氧水平依赖(BOLD)信号组成的数据集。刺激持续32秒,然后停止32秒。信号周期为64秒,采样率为每2秒1次,采样时间为256秒(\(n=128\))。

par(mfrow = c(2, 1))  
tsplot(fmri1[, 2:5], col = 1:4, ylab = "BOLD", main = "大脑皮层", spaghetti = TRUE)
tsplot(fmri1[, 6:9], col = 1:4, ylab = "BOLD", main = "丘脑和小脑", spaghetti = TRUE)
图 6: 大脑皮层、丘脑、小脑各部位的fMRI数据
tsplot(fmri1[, 2:9], col = 1:8, lwd = 2, ncol = 2, ylim = c(-.6,.6))
图 7: 分别展示
x = ts(fmri1[, 4:9], start = 0, freq = 32)         
names = c("大脑皮层", "丘脑", "小脑")
u = ts(rep(c(rep(.6, 16), rep(-.6, 16)), 4), start = 0, freq = 32)
par(mfrow = c(3, 1))
for (i in 1:3){ 
  j = 2*i - 1
  tsplot(x[, j:(j+1)], ylab = "BOLD", xlab = "", main = names[i], col = 5:6, ylim = c(-.6, .6), 
         lwd = 2, xaxt = "n", spaghetti = TRUE)
  axis(seq(0, 256, 64), side = 1, at = 0:4)
  lines(u, type = "s", col = gray(.3)) 
}
mtext("seconds", side = 1, line = 1.75, cex = .9)
图 8: 包含刺激信息

1.7 地震和爆炸

EQ5:某地区地震的纵波和横波(两阶段) EXP6:某地区爆炸的纵波和横波(两阶段)

par(mfrow = 2:1)
tsplot(EQ5,  col = 4, main = "地震")
tsplot(EXP6, col = 4, main = "爆炸")
图 9: 地震和爆炸两阶段,1秒40个采样点

2 时间序列统计模型

2.1 白噪声(White Noise)

白噪声:符合\(w_t\sim wn(0,\sigma_w^2)\)的不相关随机变量

独立同分布(independent and identically distributed,iid)的噪声:\(w_t\sim iid(0,\sigma_w^2)\)

高斯白噪声(Gaussian white noise):\(w_t\sim iidN(0,\sigma_w^2)\)

set.seed(1)
w = rnorm(500, 0, 1)                                   # 500 N(0,1) variates
v = stats::filter(w, sides = 2, filter = rep(1/3, 3))  # moving average
par(mfrow = c(2, 1))
tsplot(w, col = 4, main = "white noise")
tsplot(v, col = 4, ylim = c(-3, 3), main = "moving average")
图 10: 高斯白噪声序列(上)和高斯白噪声序列的三点移动平均值(下)

2.2 移动平均值(Moving Averages)和过滤(Filtering)

用移动平均值\(v_t\)替换,从而对\(w_t\)进行平滑处理。图 10下方图片就是三点移动平均值,其中\(v_t=\frac{1}{3}(w_{t-1}+w_t+w_{t+1})\)。从实际效果来讲,实现了较慢的振荡被突出,较快的振荡被掩盖。这种序列一般称为滤波序列(filtered series)。

filter()函数通俗来讲是用来消除噪声的(滤波),生成图 10的代码中,参数filter = rep(1/3, 3)其实就是filter = c(1/3, 1/3, 1/3),分别表示\(w_{t+1}\)\(w_t\)\(w_{t-1}\)的系数。

filter(x,                # ts数据
       filter,           # 一个逆时间顺序的滤波器系数向量
       method = c("convolution", "recursive"), # 卷积和递归,卷积是移动平均,递归是自回归,默认卷积
       sides = 2,        # 只用于卷积,1表示单边滤波(只考虑前面的值),2表示双边滤波(两边都考虑)
       circular = F,     # 只用于卷积,是否末尾和开头相连
       init)             # 只用于递归,初始值,默认0

2.3 自回归(Autoregressions)

可以换种方式替换\(w_t\),即\(x_t=x_{t-1}-0.9x_{t-2}+w_t\)

set.seed(1)
w = rnorm(550,0,1)  # 50 extra to avoid startup problems
x = stats::filter(w, filter = c(1, -.9), method = "recursive")[-(1:50)]
tsplot(x, col = 4, main = "自回归")
图 11: 自回归序列

生成图 11的代码中,init参数默认为0,即\(x_1=w_1\)\(x_2=x_1+w_2=w_1+w_2\),之后\(x_3=x_2-0.9x_1+w_3\)\(x_1\)\(x_2\)显然是不符合要求的,常用方法是去除这一部分,生成图 11的代码中多生成了50个数据,之后去除了。

2.4 带漂移项(Drift)的随机游走(Random Walk)

带漂移项的随机游走模型适合处理图 2所展示的全球变暖相关数据。该模型中

\[ x_t=\delta+x_{t-1}+w_t \tag{1}\]

其中\(t=1,\ 2,\ \cdots\)\(x_0=0\)\(w_t\)为白噪声,常数\(\delta\)为漂移项。当\(\delta=0\)时,式 1为随机游走。

式 1还可以转化为\(x_t=\delta t+\sum_{j=1}^{t}w_j\),可见当\(\delta=0\)时,\(x_t\)就是白噪声累加和。

set.seed(154)
w = rnorm(200)
x = cumsum(w)
wd = w + .2
xd = cumsum(wd)
tsplot(xd, ylim = c(-5, 55), main = "random walk", ylab = '')
lines(x, col = 4) 
clip(0, 200, 0, 50)
abline(h = 0, col = 4, lty = 2)
abline(a = 0, b = .2, lty = 2)
图 12: 随机游走图

图 12\(\sigma_w=1\),黑色实线的漂移项\(\delta=0.2\),黑色虚线是斜率为0.2的直线,蓝色实线为随机游走(无漂移项,即\(\delta=0\)),蓝色虚线表示\(y=0\)的直线。

2.5 噪声信号(Signal in Noise)

很多时间序列模型都是假设由固定周期性变化的信号和随机噪声叠加构成的,图 6所表示的fMRI序列就有很明显的周期性。

现在令:

\[ x_t=2\cos(2\pi\frac{t+15}{50})+w_t \tag{2}\]

其中\(t=1,\ 2,\ \cdots,\ 500\),第一项\(2\cos(2\pi\frac{t+15}{50})\)就是信号。

正弦波形还可以写成:

\[ A\cos(2\pi\omega t+\phi) \]

其中\(A\)是振幅,\(\omega\)是振荡频率,\(\phi\)是相位。式 2\(A=2\)\(\omega=\frac{1}{50}\)(50个时间点为一个周期),\(\phi=\frac{2\pi15}{50}=0.6\pi\)

cs = 2 * cos(2 * pi * (1:500) / 50 + .6 * pi)
w = rnorm(500, 0, 1)
par(mfrow = c(3, 1))
tsplot(cs, ylab = "", main = expression(x[t] == 2 * cos(2 * pi * t / 50 + .6 * pi)))
tsplot(cs + w, ylab = "", main = expression(x[t] == 2 * cos(2 * pi * t / 50 + .6 * pi) + N(0, 1)))
tsplot(cs + 5 * w, ylab = "", main = expression(x[t] == 2 * cos(2 * pi * t / 50 + .6 * pi) + N(0, 25)))
图 13: 余弦波和可加性高斯白噪声污染的余弦图

图 13中间的图和下面的图分别加上了\(\sigma_w=1\)\(\sigma_w=5\)的白噪声。信号被遮掩的程度取决于信号幅度和\(\sigma_w\)的大小。信号幅度和\(\sigma_w\)之比称为信噪比(signal-to-noise ratio,SNR)。SNR越大,检测信号就越容易,反之越难。例如图 13中,中间的图就比较容易看出信号,下面的图就比较困难。通常我们关心的是信号。

3 相关性测量

相关性是时间序列分析的基本特征,所以最常用的描述性度量是协方差和相关函数。边际分布函数\(F_t(x)=P\{x_t\leq x\}\)、相应的边际密度函数\(f_t(x)=\frac{\partial F_t(x)}{\partial x}\)和平均函数都是有用的边际描述指标。

3.1 平均函数

\[ \mu_{xt}=\operatorname{E}(x_t)=\int_{-\infty}^{\infty}xf_t(x)\operatorname{d}x \tag{3}\]

3.1.1 移动平均序列的平均函数

平均函数和白噪声相同,以图 10里的情况为例。

白噪声序列\(w_t\)\(\mu_{wt}=\operatorname{E}(w_t)=0\),进行平滑后有\(\mu_{vt}=\operatorname{E}(v_t)=\frac{1}{3}[\operatorname{E}(w_{t-1})+\operatorname{E}(w_t)+\operatorname{E}(w_{t-1})]=0\),平均函数不变。

3.1.2 带漂移项的随机游走的平均函数

平均函数为一条直线,以式 1为例。

已知\(\operatorname{E}(w_t)=0\)\(\delta\)为常数,有:

\[ \mu_{xt}=\operatorname{E}(x_t)=\delta t+\sum_{j=1}^{t}\operatorname{E}(w_j)=\delta t \]

这是一条斜率为\(\delta\)的直线,图 12中的虚线就是其平均函数。

3.1.3 信号加噪声的平均函数

平均函数为余弦波,以式 2为例。

\[ \mu_{xt}=\operatorname{E}(x_t)=\operatorname{E}[2cos(2\pi\frac{t+15}{50})+w_t]=2cos(2\pi\frac{t+15}{50})+\operatorname{E}(w_t)=2cos(2\pi\frac{t+15}{50}) \]

3.2 自协方差函数(Autocovariance Function)

对于任意\(s\)\(t\)自协方差函数被定义为二阶矩函数:

\[ \gamma_x(s,t)=\operatorname{cov}(x_s,x_t)=\operatorname{E}[(x_s-\mu_s)(x_t-\mu_t)] \tag{4}\]

自协方差衡量在不同时间观察到的同一时间序列的两点间的线性依赖性。平滑的时间序列,即使\(t\)\(s\)相隔很远,自协方差函数也会比较大;而波动大的时间序列,自协方差函数几乎为0。

从经典统计数据角度来看,如果\(\gamma_x(s,t)=0\)\(x_s\)\(x_t\)非线性相关(可能有其他依赖结构)。但如果\(x_s\)\(x_t\)是二元正态分布,\(\gamma_x(s,t)=0\)说明相互独立。对于\(s=t\),自协方差就是方差,即:

\[ \gamma_x(t,t)=\operatorname{E}[(x_t-\mu_t)^2]=\operatorname{var}(x_t) \]

3.2.1 白噪声的自协方差

白噪声序列\(w_t\)\(\operatorname{E}(w_t)=0\),且

\[ \gamma_w(s,t)=\operatorname{cov}(w_s,w_t)=\begin{cases} \sigma_w^2 & s=t\\ 0 & s\neq t \end{cases} \]

3.2.2 线性组合的协方差

我们经常要计算滤波序列之间的自协方差,可以用到下面的性质:

如果有随机变量

\[ U=\sum_{j=1}^ma_jX_j\ 和\ V=\sum_{k=1}^rb_kY_k \]

分别是随机变量\(\{X_j\}\)\(\{Y_k\}\)的线性组合,则它们的协方差为:

\[ \operatorname{cov}(U,V)=\sum_{j=1}^m\sum_{k=1}^ra_jb_k\operatorname{cov}(X_j,Y_k) \]

\(\operatorname{var}(U)=\operatorname{cov}(U,U)\)

3.2.3 移动平均值的自协方差

小节 2.2中的例子为例。变换后的滤波序列有:

\[ \gamma_v(s,t)=\operatorname{cov}(v_s,v_t)=\operatorname{cov}\{\frac{1}{3}(w_{s-1}+w_s+w_{s+1}),\frac{1}{3}(w_{t-1}+w_t+w_{t+1})\} \]

\(s=t\)时:

\[ \begin{aligned} \gamma_v(t,t)&=\frac{1}{9}\operatorname{cov}\{(w_{t-1}+w_t+w_{t+1}),(w_{t-1}+w_t+w_{t+1})\}\\ &=\frac{1}{9}[\operatorname{cov}(w_{t-1},w_{t-1})+\operatorname{cov}(w_{t},w_{t})+\operatorname{cov}(w_{t+1},w_{t+1})]\\ &=\frac{3}{9}\sigma_w^2 \end{aligned} \]

\(s=t+1\)时:

\[ \begin{aligned} \gamma_v(t+1,t)&=\frac{1}{9}\operatorname{cov}\{(w_t+w_{t+1}+w_{t+2}),(w_{t-1}+w_t+w_{t+1})\}\\ &=\frac{1}{9}[\operatorname{cov}(w_t,w_t)+\operatorname{cov}(w_{t+1},w_{t+1})+\operatorname{cov}(w_{t+2},w_{t-1})]\\ &=\frac{1}{9}[\operatorname{cov}(w_t,w_t)+\operatorname{cov}(w_{t+1},w_{t+1})+0]\\ &=\frac{2}{9}\sigma_w^2 \end{aligned} \]

通过类似计算,我们得到\(\gamma_v(t-1,t)=\gamma_v(t+1,t)=\frac{2}{9}\sigma_w^2\)\(\gamma_v(t-2,t)=\gamma_v(t+2,t)=\frac{1}{9}\sigma_w^2\),当\(|t-s|>2\)时为0。即:

\[ \gamma_v(s,t)=\begin{cases} \frac{3}{9}\sigma_w^2&s=t\\ \frac{2}{9}\sigma_w^2&|s-t|=1\\ \frac{1}{9}\sigma_w^2&|s-t|=2\\ 0&|s-t|>2 \end{cases} \]

平滑后序列的协方差函数随两时间点间的间隔增加而减小,其协方差只取决于时间间隔或滞后,不取决于时间点在时间序列中的绝对位置。

3.2.4 随机游走的自协方差

随机游走模型\(x_t=\sum_{j=1}^tw_j\),有

\[ \gamma_x(s,t)=\operatorname{cov}(x_s,x_t)=\operatorname{cov}(\sum_{j=1}^sw_j,\sum_{k=1}^tw_k)=\operatorname{min}\{s,t\}\sigma_w^2 \]

随机游走的自协方差函数取决于特定的时间值\(s\)\(t\),不取决于时间间隔或滞后。且\(\operatorname{var}(x_t)=\gamma_x(t,t)=t\sigma_w^2\)

3.3 自相关函数(Autocorrelation Function,ACF)

\[ \rho(s,t)=\frac{\gamma(s,t)}{\sqrt{\gamma(s,s)\gamma(t,t)}} \]

ACF用来衡量序列中不同时间点之间的依赖性,常用来检测周期性、趋势和模式。

3.3.1 交叉协方差函数(Cross-covariance Function)

\[ \gamma_{xy}(s,t)=\operatorname{cov}(x_s,y_t)=\operatorname{E}[(x_s-\mu_{xs})(y_t-\mu_{yt})] \]

交叉协方差函数用于衡量两个时间序列在不同滞后下的协方差。

3.3.2 交叉相关函数(cross-correlation function,CCF)

\[ \rho_{xy}(s,t)=\frac{\gamma_{xy}(s,t)}{\sqrt{\gamma_x(s,s)\gamma_y(t,t)}} \]

CCF用于衡量两个时间序列之间在不同滞后下的线性关系。它类似于自相关函数,但用于评估两个不同时间序列之间的相关性,而不仅仅是同一个序列的内部相关性。

4 平稳时间序列(Stationary Time Series)

4.1 严格平稳(Strictly Stationary)

如果时间序列每个值的集合\(\{x_{t_1},\ x_{t_2},\ \cdots,\ x_{t_k}\}\)的概率等于时间移动后值的集合\(\{x_{t_1+h},\ x_{t_2+h},\ \cdots,\ x_{t_k+h}\}\)的概率,即

\[ \Pr\{x_{t_1}\le c_1,\cdots,x_{t_k}\le c_k\}=\Pr\{x_{t_1+h}\le c_1,\cdots,x_{t_k+h}\le c_k\} \tag{5}\]

其中所有\(k=1,\ 2,\ \cdots\),所有时间点\(t_1,\ t_2,\ \cdots,\ t_k\),所有数值\(c_1,\ c_2,\ \cdots,\ c_k\)和所有时间移动\(h=0,\ \pm1,\ \pm2,\ \cdots\)

\(k=1\)时,式 5变为对任意\(s\)\(t\)

\[ \Pr\{x_s\le c\}=\Pr\{x_t\le c\} \tag{6}\]

可见严格平稳的时间序列对于所有\(s\)\(t\)\(\mu_s=\mu_t\),因此\(\mu_t\)必定为常数。

\(k=2\)时,式 5又变为对任意\(s\)\(t\)以及位移\(h\)

\[ \Pr\{x_s\le c_1,x_t\le c_2\}=\Pr\{x_{s+h}\le c_1,x_{t+h}\le c_2\} \tag{7}\]

根据式 6式 7可推知,对任意\(s\)\(t\)以及位移\(h\)

\[ \gamma(s,t)=\gamma(s+h,t+h) \]

严格平稳的时间序列自协方差函数取决于\(s\)\(t\)之间的时间差异,而不是具体时间点。

4.2 弱平稳(Weakly Stationary)

弱平稳时间序列\(x_t\)是一个方差有限的随机过程,满足以下条件:

  1. 式 3中定义的均值函数\(\mu_t\)是常数,不依赖于时间\(t\)

  2. 式 4中定义的自协方差函数\(\gamma(s,t)\)仅取决于\(s\)\(t\)的差值\(|s-t|\)而不是具体的\(s\)或者\(t\)的取值。

提示

条件1和2可推导出方差恒定,同时方差恒定和自协方差函数仅取决于时差可以推导出均值恒定,但均值恒定和方差恒定不能推导出自协方差函数仅取决于时差。

严格平稳要求序列的所有阶矩(如均值、方差、偏度、峰度等)都不随时间变化(如抛硬币,理想情况下随时随刻抛正反面分布情况不变),弱平稳的要求则宽松一些,只要求前两阶矩(均值、方差)恒定。

阶矩用来衡量数据的分布情况。

一阶矩\(\operatorname{E}[X]\)衡量均值(\(\mu\));二阶矩\(\operatorname{E}[(X-\mu)^2]\)衡量方差(\(\sigma^2\));三阶矩\(\operatorname{E}[(X-\mu)^3]\)衡量偏度(\(\gamma_1\));四阶矩\(\operatorname{E}[(X-\mu)^4]\)衡量峰度(\(\gamma_2\));五阶矩表达式为\(\operatorname{E}[(X-\mu)^5]\),衡量更高的维度。

我们很难从单个数据集中评估出严格平稳,更多的是使用弱平稳。

注释

为了方便,之后使用平稳表示弱平稳,严平稳表示严格平稳。

方差有限的严平稳也是弱平稳。

注释

平稳时间序列的均值\(\operatorname{E}(x_t)=\mu_t\)与时间\(t\)无关,方便起见用\(\mu\)表示\(\mu_t\)

平稳时间序列\(x_t\)的自协方差函数\(\gamma(s,t)\)仅取决于时差\(|s-t|\),令\(s=t+h\)\(h\)表示时移或滞后,则

\[ \gamma(t+h,t)=\operatorname{cov}(x_{t+h},x_t)=\operatorname{cov}(x_h,x_0)=\gamma(h,0) \]

方便起见简写为\(\gamma(h)\)

平稳时间序列的自协方差函数为:

\[ \gamma(h)=\operatorname{cov}(x_{t+h},x_t)=\operatorname{E}[(x_{t+h}-\mu)(x_t-\mu)] \]

平稳时间序列的自相关函数(ACF)为:

\[ \rho(h)=\frac{\gamma(t+h,t)}{\sqrt{\gamma(t+h,t+h)\gamma(t,t)}}=\frac{\gamma(h)}{\gamma(0)} \tag{8}\]

4.2.1 白噪声的平稳性

小节 2.1小节 2.3所展示的白噪声序列均值为\(\mu_{wt}=0\),自协方差函数为:

\[ \gamma_w(h)=\operatorname{cov}(w_{t+h},w_t)=\begin{cases} \sigma_w^2&h=0\\ 0&h\neq0\\ \end{cases} \]

白噪声为弱平稳时间序列,如果白噪声符合正态分布或高斯分布,则还是严格平稳的。因此自相关函数为\(\rho_w(0)=1\),且当\(h\neq0\)时,\(\rho(h)=0\)

提示

书中似乎错误引用式 5,应该是根据式 8计算,自相关函数为: \[ \rho_w(h)=\begin{cases} \frac{\gamma_w(0)}{\gamma_w(0)}=\frac{\sigma_w^2}{\sigma_w^2}=1&h=0\\ \frac{\gamma_w(h)}{\gamma_w(0)}=\frac{0}{\sigma_w^2}=0&h\neq0 \end{cases} \]

4.2.2 移动平均序列的平稳性

小节 2.2中的三点移动平均过程是平稳的。根据小节 3.1.1小节 3.2.3中的内容可知其均值为\(\mu_{vt}=0\),自协方差函数为:

\[ \gamma_v(h)=\begin{cases} \frac{3}{9}\sigma_w^2&h=0\\ \frac{2}{9}\sigma_w^2&h\pm1\\ \frac{3}{9}\sigma_w^2&h\pm2\\ \frac{3}{9}\sigma_w^2&|h|>2\\ \end{cases} \]

与时间\(t\)无关。自相关函数为:

\[ \rho_v(h)=\begin{cases} 1&h=0\\ \frac{2}{3}&h\pm1\\ \frac{1}{3}&h\pm2\\ 0&|h|>2\\ \end{cases} \]

其关于滞后值0对称。

4.2.3 随机游走是非平稳的

小节 3.1.2小节 3.2.4可知,随机游走的平均函数为\(\mu_{xt}=\delta t\),自协方差函数为\(\gamma_x(s,t)=\operatorname{min}\{s,t\}\sigma_w^2\),都和时间有关,故随机游走不是平稳的。

4.2.4 趋势平稳性(Trend Stationarity)

趋势平稳性指模型围绕线性趋势的平稳行为,即均值函数不独立于时间,但自协方差函数独立于时间。

\(x_t=\alpha+\beta t+y_t\),其中\(y_t\)是平稳的,则均值函数为\(\mu_{x,t}=\operatorname{E}(x_t)=\alpha+\beta t+\mu_y\),不独立于时间;自协方差函数为\(\gamma_x(h)=\operatorname{cov}(x_{t+h},x_t)=\operatorname{E}[(x_{t+h}-\mu_{x,t+h})(x_t-\mu_{x,t})]=\operatorname{E}[(y_{t+h}-\mu_y)(y_t-\mu_y)]=\gamma_y(h)\),独立于时间。

例子如图 14所示。

图 14: 鸡肉价格

这种数据有以下几个特殊属性:

  1. \(\gamma(h)\)非负;
  2. \(\gamma(0)=\operatorname{E}[(x_t-\mu)^2]\)
  3. 自协方差函数关于原点对称,\(\gamma(h)=\gamma(-h)\)

4.3 联合平稳(Jointly Stationary)

两个时间序列\(x_t\)\(y_t\)均是平稳的,且它们的交叉协方差函数为:

\[ \gamma_{xy}(h)=\operatorname{cov}(x_{t+h},y_t)=\operatorname{E}[(x_{t+h}-\mu_x)(y_t-\mu_y)] \]

我们称\(x_t\)\(y_t\)是联合平稳的,且交叉协方差函数仅和滞后值\(h\)有关。

4.3.1 交叉相关函数(Cross-Correlation Function,CCF)

联合平稳时间序列\(x_t\)\(y_t\)的交叉相关函数为:

\[ \rho_{xy}(h)=\frac{\gamma_{xy}(h)}{\sqrt{\gamma_x(0)\gamma_y(0)}} \tag{9}\]

其交叉相关函数通常关于0不对称,即\(\rho_{xy}(h)\neq \rho_{xy}(-h)\),即\(\operatorname{cov}(x_2,y_1)\)\(\operatorname{cov}(x_1,y_2)\)不必相同。

\(\rho_{xy}(h)=\rho_{yx}(-h)\)

提示

这里给出推导过程。

\(\rho_{xy}(h)\neq \rho_{xy}(-h)\)是因为:

\[ \rho_{xy}(h)=\frac{\gamma_{xy}(h)}{\sqrt{\gamma_x(0)\gamma_y(0)}} \]

\[ \rho_{xy}(-h)=\frac{\gamma_{xy}(-h)}{\sqrt{\gamma_x(0)\gamma_y(0)}} \]

中,\(\gamma_{xy}(h)\neq \gamma_{xy}(-h)\),即:

\[ \operatorname{cov}(x_{t+h},y_t)\neq\operatorname{cov}(x_t,y_{t+h}) \]

\(\rho_{xy}(h)=\rho_{yx}(-h)\)是因为:

\[ \rho_{yx}(-h)=\frac{\gamma_{yx}(-h)}{\sqrt{\gamma_y(0)\gamma_x(0)}} \]

根据协方差的性质:

\[ \operatorname{cov}(x_t,y_{t+h})=\operatorname{E}[(x_t-\mu_x)(y_{t+h}-\mu_y)]=\operatorname{E}[(y_{t+h}-\mu_y)(x_t-\mu_x)]=\operatorname{cov}(y_{t+h},x_t) \]

即:

\[ \operatorname{cov}(A,B)=\operatorname{cov}(B,A) \]

可得:

\[ \gamma_{yx}(-h)=\operatorname{cov}(y_{t-h},x_t)=\operatorname{cov}(x_t,y_{t-h})=\gamma_{xy}(h) \]

4.3.2 联合平稳性(Joint Stationarity)

两个序列\(x_t\)\(y_t\),分别是白噪声过程两个连续值的和与差,即\(x_t=w_t+w_{t-1}\)\(y_t=w_t-w_{t-1}\),其中\(w_t\)是均值为0,方差为\(\sigma_w^2\)的独立随机变量。易证\(\gamma_x(0)=\gamma_y(0)=2\sigma_w^2\),以及\(\gamma_x(1)=\gamma_y(-1)=\sigma_w^2\)\(\gamma_y(1)=\gamma_x(-1)=-\sigma_w^2\)

同时有:

\[ \gamma_{xy}(1)=\operatorname{cov}(x_{t+1},y_t)=\operatorname{cov}(w_{t+1}+w_t,w_t-w_{t-1})=\sigma_w^2 \tag{10}\]

类似的:\(\gamma_{xy}(0)=0\)\(\gamma_{xy}(-1)=-\sigma_w^2\)。根据式 9,可得:

\[ \rho_{xy}(h)=\begin{cases} 0&h=0\\ 1/2&h=1\\ -1/2&h=-1\\ 0&|h|\ge2 \end{cases} \]

自协方差和交叉协方差函数仅取决于滞后值\(h\),因此序列是联合平稳的。

提示

给出补充的推导过程帮助记忆和理解。

注意\(\gamma_{xy}(h)\)等于\(\operatorname{cov}(x_{t+h},y_t)\),而不是等于\(\operatorname{cov}(x_t,y_{t+h})\),在自协方差里面混淆了不影响结果,但在交叉协方差里面是有影响的。

  1. 方差

因为是独立随机变量,所以\(\operatorname{cov}(w_t,w_{t-1})=0\),所以:

\[ \gamma_x(0)=\operatorname{var}(x_t)=\operatorname{var}(w_t+w_{t-1})=\operatorname{var}(w_t)+\operatorname{var}(w_{t-1})+2\operatorname{cov}(w_t,w_{t-1})=2\sigma_w^2 \]

同理:

\[ \gamma_y(0)=\operatorname{var}(y_t)=\operatorname{var}(w_t-w_{t-1})=\operatorname{var}(w_t)+\operatorname{var}(w_{t-1})-2\operatorname{cov}(w_t,w_{t-1})=2\sigma_w^2 \]

  1. 滞后自协方差

\[ \begin{aligned} \gamma_x(1)=\operatorname{cov}(x_{t+1},x_t) &=\operatorname{cov}(w_{t+1}+w_t,w_t+w_{t-1})\\ &=\operatorname{cov}(w_{t+1},w_t)+\operatorname{cov}(w_{t+1},w_{t-1})+\operatorname{cov}(w_t,w_t)+\operatorname{cov}(w_t,w_{t-1})\\ &=0+0+\sigma_w^2+0\\ &=\sigma_w^2 \end{aligned} \]

类似的:

\[ \begin{aligned} \gamma_y(1)=\operatorname{cov}(y_{t+1},y_t) &=\operatorname{cov}(w_{t+1}-w_t,w_t-w_{t-1})\\ &=\operatorname{cov}(w_{t+1},w_t)-\operatorname{cov}(w_{t+1},w_{t-1})-\operatorname{cov}(w_t,w_t)+\operatorname{cov}(w_t,w_{t-1})\\ &=0-0-\sigma_w^2+0\\ &=-\sigma_w^2 \end{aligned} \]

  1. 自协方差

\[ \gamma_x(h)=\begin{cases} 2\sigma_w^2&h=0\\ \sigma_w^2&|h|=1\\ 0&|h|\ge2\\ \end{cases} \]

\[ \gamma_y(h)=\begin{cases} 2\sigma_w^2&h=0\\ -\sigma_w^2&|h|=1\\ 0&|h|\ge2\\ \end{cases} \]

交叉相关函数是在式 10的基础上,简单的公式代入,不做补充推导。

4.3.3 使用交叉相关进行预测

有模型:

\[ y_t=Ax_{t-\ell}+w_t \]

前导(lead):当\(\ell>0\)时,序列\(x_t\)是序列\(y_t\)的前导;

滞后(lag):当\(\ell<0\)时,序列\(x_t\)是序列\(y_t\)的滞后。

在用\(x_t\)预测\(y_t\)时,判断清楚前导和滞后关系很重要。

假设噪声\(w_t\)\(x_t\)不相关,则交叉协方差函数如下:

\[ \begin{aligned} \gamma_{yx}(h) &=\operatorname{cov}(y_{t+h},x_t)=\operatorname{cov}(Ax_{t+h-\ell}+w_{t+h},x_t)\\ &=\operatorname{cov}(Ax_{t+h-\ell},x_t)=A\gamma_x(h-\ell) \end{aligned} \]

根据柯西不等式(Cauchy-Schwarz),当\(h=\ell\)时,即\(\gamma_x(0)\)时,\(\gamma_{yx}(h)\)取最大绝对值。如果\(h>0\)时取最大绝对值,说明\(x_t\)\(y_t\)的先导;如果\(h<0\)时取最大绝对值,说明\(x_t\)滞后于\(y_t\)

提示

因为是相关,不涉及因果,所以先导和滞后都是相对的。

考虑到标准差一般为固定值,协方差最大的时候,就是\(x_t\)\(y_t\)相关性最强的时候。

set.seed(2)
x = rnorm(100)
y = stats::lag(x, -5) + rnorm(100) # h = 5
astsa::ccf2(y, x, ylab = 'CCovF', type = 'covariance')
text( 9, 1.1, 'x leads')
text(-8, 1.1, 'y leads')
图 15: 交叉协方差函数

4.4 线性(随机)过程(Linear Process)

\(x_t\)是白噪声\(w_t\)的线性组合,具体如下:

\[ x_t=\mu+\sum_{j=-\infty}^{\infty}\psi_jw_{t-j},\ \sum_{j=-\infty}^{\infty}|\psi_j|<\infty \tag{11}\]

\(h\ge0\)时,其自协方差函数为:

\[ \gamma_x(h)=\sigma_w^2\sum_{j=-\infty}^{\infty}\psi_{j+h}\psi_j \tag{12}\]

提示
  1. 模型解释

式 11表示一种广义的线性随机过程。

\[ x_t=\mu+\sum_{j=-\infty}^{\infty}\psi_jw_{t-j},\ \sum_{j=-\infty}^{\infty}|\psi_j|<\infty \]

其中:

  • \(x_t\):当前时间点\(t\)的观测值;

  • \(\mu\):序列的均值;

  • \(\psi_j\):权重系数,描述了过去(\(j>0\))、现在(\(j=0\))和未来(\(j < 0\))的随机扰动对当前观测值的影响;

  • \(w_t\):随机扰动项。

  1. 自协方差函数推导

\[ \gamma_x(h)=\operatorname{cov}(x_{t+h},x_t) \]

将给定的\(x_t\)代入:

\[ \begin{aligned} \gamma_x(h) &=\operatorname{cov}(\mu+\sum_{j=-\infty}^{\infty}\psi_jw_{t+h-j},\mu+\sum_{k=-\infty}^{\infty}\psi_kw_{t-k})\\ &=\operatorname{cov}(\sum_{j=-\infty}^{\infty}\psi_jw_{t+h-j},\sum_{k=-\infty}^{\infty}\psi_kw_{t-k})\\ &=\sum_{j=-\infty}^{\infty}\sum_{k=-\infty}^{\infty}\psi_j\psi_k\operatorname{cov}(w_{t+h-j},w_{t-k}) \end{aligned} \]

\(k=j-h\)时,\(\operatorname{cov}(w_{t+h-j},w_{t-k})=\sigma_w^2\)

\(k\neq j-h\)时,\(\operatorname{cov}(w_{t+h-j},w_{t-k})=0\)

因此,非零项只发生在\(k=j-h\)处,表达式可简化为:

\[ \gamma_x(h)=\sigma_w^2\sum_{j=-\infty}^{\infty}\psi_j\psi_{j-h} \]

这和式 12是一致的。

提示

正态时间序列(高斯时间序列)之后也会涉及,目前来看在理解上没那么重要,如有需求会在后续补上。

5 相关系数的估计

上面论述的内容只适用于假设的模型,但大部分研究都是通过样本数据来分析整体的,因此我们要在平稳性假设的基础上进行分析。

如果时间序列是平稳的,均值函数是一个常数,可以用样本均值来估计:

\[ \bar x=\frac{1}{n}\sum_{t=1}^nx_t \]

方差为:

\[ \begin{aligned} \operatorname{var}(\bar{x}) &=\operatorname{var}(\frac{1}{n}\sum_{t=1}^nx_t)=\frac{1}{n^2}\operatorname{cov}(\sum_{t=1}^{n}x_t,\sum_{s=1}^{n}x_s)\\ &\begin{aligned}= &\frac{1}{n^2}(n\gamma_x(0)+(n-1)\gamma_x(1)+(n-2)\gamma_x(2)+\cdots+\gamma_x(n-1)\\ &+(n-1)\gamma_x(-1)+(n-2)\gamma_x(-2)+\cdots+\gamma_x(1-n)) \end{aligned}\\ &=\frac{1}{n}\sum_{h=-n}^{n}(1-\frac{|h|}{n})\gamma_x(h) \end{aligned} \tag{13}\]

提示

式 13补充部分推导过程。

\[ \frac{1}{n^2}\operatorname{cov}(\sum_{t=1}^{n}x_t,\sum_{s=1}^{n}x_s)=\frac{1}{n^2}\sum_{t=1}^{n}\sum_{s=1}^{n}\operatorname{cov}(x_t,x_s)=\frac{1}{n^2}\sum_{t=1}^{n}\sum_{s=1}^{n}\gamma_x(t-s) \]

\(n\)\(t-s=0\)\(n-1\)\(t-s=1\)\(t-s=-1\)\(n-2\)\(t-s=2\)\(t-s=-2\);……

展开后,因为\(\gamma_x(n)=\gamma_x(-n)\),令\(h=t-s\),合并得到:

\[ \frac{1}{n}\sum_{h=-n}^{n}(1-\frac{|h|}{n})\gamma_x(h) \]

5.1 样本自协方差函数(Sample Autocovariance Function)

\[ \hat\gamma(h)=n^{-1}\sum_{t=1}^{n-h}(x_{t+h}-\bar x)(x_t-\bar x) \]

其中\(h=0,1,\cdots,n-1\),且\(\hat\gamma(-h)=\hat\gamma(h)\)

5.2 样本自相关函数(Sample Autocorrelation Function)

\[ \hat\rho(h)=\frac{\hat\gamma(h)}{\hat\gamma(0)} \]

使用小节 1.5中的SOI序列数据,展示不同滞后下得到的不同估计相关系数,以及其散点图。图片左上角为估计相关系数。

r = round(acf1(soi, 6, plot = FALSE), 2) # sample acf values
par(mfrow = c(1, 2))
tsplot(stats::lag(soi, -1), soi, col = 4, type = 'p', xlab = 'lag(soi,-1)')
 legend("topleft", legend = r[1], bg = "white", adj = .45, cex = 0.85)
tsplot(stats::lag(soi, -6), soi, col = 4, type = 'p', xlab = 'lag(soi,-6)')
 legend("topleft", legend = r[6], bg = "white", adj = .25, cex = 0.8)
图 16: SOI序列与其滞后1个月(左)和滞后6个月(右)的相关散点图

如果在大样本下,样本ACF近似服从正态分布,该正态分布均值为0,标准差为\(\frac{1}{\sqrt{n}}\)

在此基础上,我们还是可以根据\(\hat\rho(h)\)是否在95%置信区间外来检验显著性。

通过模拟数据,我们来验证一下是不是大样本数据估计出来的ACF更接近理论值。

用抛硬币的方法构建一组数据\(x_t\),硬币正面令\(x_t=1\),反面令\(x_t=-1\)。然后构建\(y_t=5+x_t-0.7x_{t-1}\)

set.seed(666)
x1 = sample(c(-1,1), 11, replace=TRUE)  # simulated sequence of coin tosses
x2 = sample(c(-1,1), 1001, replace=TRUE)
y1 = 5 + stats::filter(x1, sides=1, filter=c(1,-.7))[-1]
y2 = 5 + stats::filter(x2, sides=1, filter=c(1,-.7))[-1]
c(mean(y1), mean(y2))  # the sample means
[1] 4.8600 5.0018
acf(y1, lag.max=4, plot=FALSE) 

Autocorrelations of series 'y1', by lag

     0      1      2      3      4 
 1.000 -0.226 -0.630  0.221  0.419 
acf(y2, lag.max=4, plot=FALSE)

Autocorrelations of series 'y2', by lag

     0      1      2      3      4 
 1.000 -0.474 -0.003  0.001 -0.008 

计算一下理论值。

\(x_t\)均值为0,方差为1,则\(y_t\)有:

\[ \rho_x(1)=\frac{-0.7}{1+0.7^2}\approx-0.47 \]

可见样本更大的更接近理论值。

提示

理论值的推导。

\[ \begin{aligned} \gamma_y(h)&=\operatorname{cov}(y_{t+h},y_t)=\operatorname{cov}(5+x_{t+h}-0.7x_{t+h-1},5+x_t-0.7x_{t-1})\\ &=\operatorname{cov}(x_{t+h}-0.7x_{t+h-1},x_t-0.7x_{t-1})\\ &=\operatorname{cov}(x_{t+h},x_t)-0.7\operatorname{cov}(x_{t+h},x_{t-1})-0.7\operatorname{cov}(x_{t+h-1},x_t)+0.49\operatorname{cov}(x_{t+h-1},x_{t-1}) \end{aligned} \]

\(x_t\)是均值为0,方差为1的独立同分布,所以当\(h=0\)时,\(\operatorname{cov}(x_t,x_t)=\operatorname{var}(x_t)=1\);当\(h\neq0\)时,\(\operatorname{cov}(x_t,x_{t+h})=0\)

所以:

\[ \gamma_y(h)=\begin{cases} 1-0-0+0.49=1.49&h=0\\ 0-0-0.7+0=0.7&h=1\\ 0-0-0+0=0&h\ge2 \end{cases} \]

所以:

\[ \rho_y(1)=\frac{\gamma_y(1)}{\gamma_y(0)}=\frac{-0.7}{1+0.7^2}\approx-0.47 \]

再展示小节 1.3的ACF图。

acf1(speech, 250)
图 17: 语音数据的ACF

可见每隔100多点的间隔,峰值重复出现,说明这可能是一系列重复的短信号。

5.3 样本交叉协方差函数

\[ \hat\gamma_{xy}(h)=n^{-1}\sum_{t=1}^{n-h}(x_{t+h}-\bar x)(y_t-\bar y) \]

5.4 样本交叉相关函数

\[ \hat\rho_{xy}(h)=\frac{\hat\gamma_{xy}(h)}{\sqrt{\hat\gamma_x(0)\hat\gamma_y(0)}} \]

如果至少有一个过程是独立的白噪声,则\(\hat\rho_{xy}(h)\)的大样本分布是正态的,均值为0,标准差\(\sigma_{\hat\rho_(xy)}=\frac{1}{\sqrt{n}}\)

还是用小节 1.5中的数据举例。

par(mfrow=c(3,1))
acf1(soi, main="SOI")
acf1(rec, main="新鱼数量")
ccf2(soi, rec, main="SOI vs 新鱼数量")
图 18: SOI序列(上图)和新鱼数量序列(中图)的ACF,以及它们的CCF(下图)

图 18上图和中图表明12个月的间隔表现出正相关,6个月的间隔表现出负相关。下图则表现出两组数据的一些偏移,在\(h=-6\)处取得峰值,说明时间\(t-6\)测量的SOI和时间\(t\)的新鱼数量相关。可以说SOI领先新鱼数量6个月,且SOI的增加可能导致新鱼数量减少。虚线表示\(\pm2/\sqrt{453}\),但现在没有意义,因为两个序列都不是白噪声。

5.5 预白化(Prewhitening)和交叉相关分析

白化是指将一个时间序列转换为白噪声序列的过程,也就是去除其相关性,生成一个没有自相关性的序列,便于之后进行统计分析。

set.seed(1492)
num = 120
t = 1:num
X = ts(2 * cos(2 * pi * t / 12) + rnorm(num), freq = 12)
Y = ts(2 * cos(2 * pi * (t + 5) / 12) + rnorm(num), freq = 12)
Yw = resid(lm(Y ~ cos(2 * pi * t / 12) + sin(2 * pi * t / 12), na.action = NULL))
par(mfrow = c(3, 2))
tsplot(X)
tsplot(Y)
acf1(X, 48, ylab = 'ACF(X)')
acf1(Y, 48, ylab = 'ACF(Y)')
ccf2(X, Y, 24)
ccf2(X, Yw, 24, ylim = c(-.6, .6))
图 19: 白化处理展示

5.6 向量值(Vector-Valued)和多维序列(Multidimensional Series)

当考虑多个变量的时间序列之间的关系时,经常用向量的方法表示,即\(x_t=(x_{t1},x_{t2},\cdots,x_{tp})'\),这表示\(p\)个一元时间序列,是个\(p\times1\)的列向量。

平稳情况下均值为:

\[ \mu=\operatorname{E}(x_t)=(\mu_{t1},\mu_{t2},\cdots,\mu_{tp})' \]

\(p\times p\)协方差矩阵为:

\[ \mathbf{\Gamma}(h)=\operatorname{E}[(x_{x+h}-\mu)(x_t-\mu)']=\begin{bmatrix} \gamma_{11}(h)&\gamma_{12}(h)&\cdots&\gamma_{1p}(h)\\ \gamma_{21}(h)&\gamma_{22}(h)&\cdots&\gamma_{2p}(h)\\ \vdots&\vdots&\ddots&\vdots\\ \gamma_{p1}(h)&\gamma_{p2}(h)&\cdots&\gamma_{pp}(h)\\ \end{bmatrix} \]

因为\(\gamma_{ij}(h)=\gamma_{ji}(-h)\),所有\(\mathbf{\Gamma}(h)=\mathbf{\Gamma}'(-h)\)

样本自协方差矩阵为:

\[ \hat{\mathbf{\Gamma}}(h)=n^{-1}\sum_{t=1}^{n-h}(x_{t+h}-\bar x)(x_t-\bar x)' \]

有时候\(x_t\)本身就是个多维数据(如坐标)。

图 20展示的就是一块\(64\times36\)分割的土地上温度。\(x_s=x_{s1,s2}\)表示\(s_1\)\(s_2\)列的温度。

par(mar = rep(1, 4))
persp(
    1:64,          # x轴范围
    1:36,          # y轴范围
    soiltemp,      # 数据集
    phi = 30,      # X-Y平面角度
    theta = 30,    # X-Z平面角度
    scale = FALSE, # 数据标准化,Z轴的自动缩放
    expand = 4,    # 图的扩展比例
    ticktype = "detailed",  # 刻度线类型
    xlab = "rows",
    ylab = "cols",
    zlab = "temperature"
)
图 20: 坐标上的温度序列

图中可以看出,从第\(40\)行开始,温度沿着行坐标开始有平稳波动。对36列进行平均,绘制行均值图 21

tsplot(rowMeans(soiltemp), xlab = "row", ylab = "Average Temperature")
图 21: 行均值

可发现前面40行的噪声被平滑掉了,整体是一个稳定的温度信号。

可将一个平稳多维随机过程\(x_s\)的自协方差函数定义为多维滞后向量\(h=(h_1,h_2,\cdots,h_r)'\)的一个函数,即:

\[ \gamma(h)=\operatorname{E}[(x_{s+h}-\mu)(x_s-\mu)] \]

迁移到坐标上就是:

\[ \gamma(h_1,h_2)=\operatorname{E}[(x_{s_1+h_1,s_2+h_2}-\mu)(x_{s_1,s_2}-\mu)] \]

表示行方向滞后\(h_1\),列方向滞后\(h_2\)的自协方差函数。

多维样本自协方差函数为:

\[ \hat\gamma(h)=(S_1S_2\cdots S_r)^{-1}\sum_{s_1}\sum_{s_2}\cdots\sum_{s_r}(x_{s+h}-\bar x)(x_s-\bar x) \]

其中:

\[ \bar x=(S_1S_2\cdots S_r)^{-1}\sum_{s_1}\sum_{s_2}\cdots\sum_{s_r}x_{s_1,s_2,\cdots,s_r} \]

样本自相关函数为:

\[ \hat\rho(h)=\frac{\hat\gamma(h)}{\hat\gamma(0)} \]

提示

后面示例在之后也会出现,暂时放下。

参考文献

Shumway, R. H., & Stoffer, D. S. (2017). Time series analysis and its applications: With r examples (4th ed.). Springer Cham. https://doi.org/10.1007/978-3-319-52452-8