Time Series Analysis Chapter 2
TextBook:《时间序列分析——基于R》王燕
0.1 时序图、自相关图
0.1.1 时序图
yield<-c(15.2,16.9,15.3,14.9,15.7,15.1,16.7)
yield<-ts(yield,start = 1884)
plot(yield,main='1884-1890年英格兰和威尔士地区小麦平均亩产量时序图',xlab = "年份",ylab="亩产量") # main:a main title for the plot
接下来研究 ‘plot’ 参数的影响。
plot(yield,type="p") # 散点图
plot(yield,type = "o") # 点线图
plot(yield,type = "o",pch=17) # 点线图，其中点用17代替，详见参数pch的解释
plot(yield,lty=2) # 连线参数设置为2，参见参数lty的解释
plot(yield,lwd=2) # 线宽度设置，2表示为默认宽度的2倍
plot(yield,xlim = c(1886,1890)) # 指定坐标轴范围，输出横轴范围为1886到1890
plot(yield,ylim=c(15,16)) # 指定坐标轴范围，输出纵轴范围为15到16
plot(yield)
abline(v=1887,lty=2) # 添加参照线为一条垂直线
plot(yield)
abline(v=c(1885,1889),lty=2) #添加参照线为多条垂直参照线
plot(yield)
abline(h=c(15.5,16.5),lty=2) # 添加参照先为多条水平线
0.1.2 自相关图
acf(yield) # 绘制序列自相关图
0.2 平稳性检验
0.2.1 时序图检验
# 1. 明显递增的时间序列
# 这是1964年-1999年的中国纱年产量序列
sha<-read.table("file4.csv",sep=",",header = T)
output<-ts(sha$output,start=1964)
plot(output)
很明显，此时间序列具有明显递增的趋势，故一定不是平稳时间序列。
# 2. 呈现规则周期性、呈现明显逐年递增趋势
# 这是1962年1月至1975年12月平均每头奶牛月产奶量序列
a<-read.table("file5.csv",sep=",",header = T)
milk<-ts(a$milk,start = c(1962,1),frequency = 12)
plot(milk)
很明显，此时间序列具有规则的周期性、明显递增的趋势，故一定不是平稳时间序列。
# 3. 时序图似乎可认为是平稳序列
# 这是1949年-1998年北京市最高气温序列
b<-read.table("file6.csv",sep=",",header = T)
temp<-ts(b$temp,start = 1949)
plot(temp)
0.2.2 自相关图检验
平稳序列的自相关系数 \(\hat{\rho_k}\) 会很快地衰减为 0，而非平稳序列的自相关系数 \(\hat{\rho_k}\) 衰减向 0 的速度会非常慢。
# 1.
acf(output,lag=25)
序列的自相关系数递减到0的速度非常慢，自相关图上明显呈现三角对称性，这是具有单调趋势的非平稳序列的一种典型自相关图的趋势。
# 2.
acf(milk)
序列的自相关系数一直都大于 0 ，这是具有单调趋势序列的典型特征，
同时自相关图呈现明显的正弦波动规律，这是具有周期变化的非平稳序列的典型特征。
# 3.
acf(temp)
序列的自相关系数一直都很小，始终控制在2倍标准差的范围之内，可认为是随机性非常强的平稳时间序列具有的自相关图特征。
0.3 纯随机性检验
序列平稳的前提下，并非所有的平稳序列都值得建模，只有序列值之间具有密切的相关关系、历史数据对未来的发展具有一定影响的序列才值得建模，用于预测未来的发展。
若序列之间没有任何相关性，即意味着该序列为“没有记忆”的序列，过去的行为对未来没有任何影响，此时序列称为纯随机性序列，又称白噪声序列（white noise），记作 \(X_t \sim WN(\mu,\sigma^2)\) 。
统计角度而言，纯随机序列并没有任何分析的价值。
#随机生成1000个服从标准正态分布的白噪声序列观察值，并绘制时间序列图
white_noise<-rnorm(1000)
white_noise<-ts(white_noise)
plot(white_noise)
0.3.1 白噪声检验
0.3.1.1 自相关图检验
acf(white_noise)
纯随机序列的序列值之间应当没有任何相关关系。
0.3.1.2 Bartlett 定理用于检验白噪声：Box.test()
假设检验
\[H_0:\rho_0=\rho_1=...=\rho_m=0 \quad v.s.\quad H_1: \exists k \leq m \quad s.t. \quad \rho_k \neq 0,\quad \forall m \geq 1\]
Q 统计量：大样本场合下检验效果很好，小样本场合下不太精确。 \[Q=n\sum_{k=1}^m \hat{\rho_k}^2 \dot{\sim} \chi^2(m)\]
LB 统计量：修正 Q 统计量，各种场合下的统计量其实指的是 LB 统计量。 \[LB=n(n+2)\sum_{k=1}^m (\frac{\hat{\rho_k}^2}{n-k}) \dot{\sim} \chi^2(m)\]
Box.test(white_noise,lag = 6) # 默认情况下为 Q 统计量##
## Box-Pierce test
##
## data: white_noise
## X-squared = 5.8224, df = 6, p-value = 0.4434Box.test(white_noise,lag = 12)##
## Box-Pierce test
##
## data: white_noise
## X-squared = 14.325, df = 12, p-value = 0.2805# 为一次得到多个白噪声检验的结果，编写一个循环命令
for(i in 1:2) print(Box.test(white_noise,lag=6*i)) # 输出6、12阶延迟的结果##
## Box-Pierce test
##
## data: white_noise
## X-squared = 5.8224, df = 6, p-value = 0.4434
##
##
## Box-Pierce test
##
## data: white_noise
## X-squared = 14.325, df = 12, p-value = 0.2805由于p值较大，故不能拒绝纯随机性的假设，即我们可以认为该序列并没有任何统计规律可循，故停止对该序列的统计分析。
由于平稳序列通常具有短期相关性，若一个平稳序列短期lag序列值之间都不存在显著的相关关系，通常长期lag之间就更不会存在显著的相关关系了。
# 检验 3. 的白噪声序列
for(i in 1:2)print(Box.test(temp,lag=6*i))##
## Box-Pierce test
##
## data: temp
## X-squared = 5.0166, df = 6, p-value = 0.5417
##
##
## Box-Pierce test
##
## data: temp
## X-squared = 5.8565, df = 12, p-value = 0.9231同上分析，我们认为北京市最高气温的变动为纯随机性变动，说明我们很难根据历史信息预测未来年份的最高气温。至此对该序列的分析结束。
0.4 综合应用
这是1950年-1998年北京市城乡居民定期储蓄所占比例序列，下面对其进行平稳性与纯随机性检验
c<-read.table("file7.csv",sep=",",header = T)
prop<-ts(c$prop,start = 1950)
plot(prop)
由时间序列图可知，北京市城乡居民的定期储蓄始终占储蓄存款余额的 80% 左右，波动比较平稳。
acf(prop)
样本自相关图中，延迟 7 阶后，自相关系数均落入 2 倍标准差之内，且自相关系数向 0 衰减的速度非常快，延迟 8 阶后自相关系数在 0 值附近波动。我们可以认为这是一个典型的短期相关的样本自相关图。
由前，我们认为该序列为平稳序列。
下面进行纯随机性检验。
for(i in 1:2) print(Box.test(prop,lag=6*i))##
## Box-Pierce test
##
## data: prop
## X-squared = 68.724, df = 6, p-value = 7.467e-13
##
##
## Box-Pierce test
##
## data: prop
## X-squared = 74.74, df = 12, p-value = 4.115e-11纯随机性检验显示，p值较小，故拒绝原假设，即我们认为该序列并非为白噪声序列。
此时我们应当对该序列进行统计分析。
0.5 Homeworks：
0.5.1 1.
a1 <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
b1 <-ts(a1)
plot(b1)
由时序图可知，此序列有明显的递增趋势，故一定不是平稳序列。
acf(b1)
上述自相关图表明，随着延迟数（lag）k的增加，平稳序列的自相关系数 \(\hat{\rho_k}\) 会很快的衰减到0，然而这个序列并不是，且在相关图上显示出明显的倒三角对称性，故序列具有单调趋势，且非平稳。
0.5.2 2.
a2=c(t(read.table('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/1st_homework/test 2.2.txt',header=FALSE)))
b2 <-ts(a2,frequency=12,start=c(1975,1))
plot(b2)
由时间序列图可知，此时间序列数据具有逐年递增趋势，且呈现出规则的周期性，显然不为平稳的时间序列。
acf(b2,lag=24)
如上为样本自相关图，由图可知，该序列具有单调趋势，且由于图中呈现出明显的正弦波动规律，故此序列为具有周期变化规律的非平稳序列。
0.5.3 3.
a3=c(t(read.table('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/1st_homework/test 2.3.txt',header=FALSE)))
b3 <-ts(a3,frequency=12,start=c(1945,1))
c3=acf(b3,lag=24)
c3##
## Autocorrelations of series 'b3', by lag
##
## 0.0000 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333
## 1.000 0.013 0.042 -0.043 -0.179 -0.251 -0.094 -0.068 -0.072 0.014 0.109
## 0.9167 1.0000 1.0833 1.1667 1.2500 1.3333 1.4167 1.5000 1.5833 1.6667 1.7500
## 0.217 0.316 -0.025 0.075 -0.141 -0.204 -0.245 0.066 -0.139 -0.034 0.206
## 1.8333 1.9167 2.0000
## -0.010 0.080 0.118如上为k=1,2,…,24的样本自相关系数值，与样本自相关图，由图可知，该序列随着延迟数（lag）k的增加，平稳序列的自相关系数 \(\hat{\rho_k}\) 会很快的衰减到0,且落在横带内部，故可以判断为平稳序列。
for(i in 1:12) print(Box.test(b3,lag=i))##
## Box-Pierce test
##
## data: b3
## X-squared = 0.011742, df = 1, p-value = 0.9137
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 0.13635, df = 2, p-value = 0.9341
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 0.2709, df = 3, p-value = 0.9654
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 2.5699, df = 4, p-value = 0.6322
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 7.1168, df = 5, p-value = 0.2121
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 7.7505, df = 6, p-value = 0.257
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 8.0812, df = 7, p-value = 0.3255
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 8.4542, df = 8, p-value = 0.3904
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 8.4681, df = 9, p-value = 0.4877
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 9.3306, df = 10, p-value = 0.501
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 12.73, df = 11, p-value = 0.3113
##
##
## Box-Pierce test
##
## data: b3
## X-squared = 19.914, df = 12, p-value = 0.06873在滞后1-12期中，由于p-value>\(\alpha\)=0.05，故不能拒绝该序列纯随机的假设，我们可以认为此序列为white noise，具有纯随机性。
0.5.4 4.
a4<-c(0.02,0.05,0.1,-0.02,0.05,0.01,0.12,-0.06,0.08,-0.05,0.02,-0.05)
barplot(a4)
计算LB统计量，
\[LB=n(n+2)\sum_{k=1}^{m} (\frac{\hat{\rho_k^2}}{n-k})\] 其中n=100,m=12
LBStatistics = function(n,m,series){
sums=0
for(i in 1:m){
sums=sums+series[i]**2/(n-i)
}
LB=sums*n*(n+2)
return(LB)
}
LBStatistics(n=100,m=12,a4)## [1] 4.989531LB统计量对应的分位点为0.9634，p-value=0.0363<\(\alpha\)=0.05，
故拒绝原假设，即序列不能视为纯随机序列。
0.5.5 5.
a5=read.table('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/1st_homework/test 2.5.txt',header=FALSE)
a5_1=c(t(t(a5[2:nrow(a5),2:ncol(a5)])))
a5_1<-as.numeric(a5_1)
#用read.csv()命令读取后为数据框格式，两次使用t()命令后转化为矩阵格式，但其中的数字成为字符格式，电脑上显示均带有""号。
b5=ts(a5_1,frequency=12,start=2000)
plot(b5)
acf(b5)
由时间序列图可知，此时间序列数据呈现出规则的周期性。
由样本自相关图可知，此序列为具有周期变化规律的非平稳序列。
for(i in 1:12) print(Box.test(b5,lag=i))##
## Box-Pierce test
##
## data: b5
## X-squared = 26.249, df = 1, p-value = 3.001e-07
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 36.276, df = 2, p-value = 1.327e-08
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 36.592, df = 3, p-value = 5.612e-08
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 41.509, df = 4, p-value = 2.108e-08
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 59.889, df = 5, p-value = 1.281e-11
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 84.84, df = 6, p-value = 3.331e-16
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 103.85, df = 7, p-value < 2.2e-16
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 109.39, df = 8, p-value < 2.2e-16
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 109.39, df = 9, p-value < 2.2e-16
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 114.74, df = 10, p-value < 2.2e-16
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 130.87, df = 11, p-value < 2.2e-16
##
##
## Box-Pierce test
##
## data: b5
## X-squared = 156.51, df = 12, p-value < 2.2e-16在滞后1-12期中，由于p-value<\(\alpha\)=0.05，故拒绝该序列纯随机的假设，我们可以认为此序列不具有纯随机性，可以进行分析。
0.5.6 6.
a6=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/1st_homework/test 2.6.txt')
b6=ts(a6,frequency=12,start=1969)
b6## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1969 10 15 10 10 12 10 7 7 10 14 8 17
## 1970 14 18 3 9 11 10 6 12 14 10 25 29
## 1971 33 33 12 19 16 19 19 12 34 15 36 29
## 1972 26 21 17 19 13 20 24 12 6 14 6 12
## 1973 9 11 17 12 8 14 14 12 5 8 10 3
## 1974 16 8 8 7 12 6 10 8 10 5acf(b6)
由样本自相关图可知，此序列不具有平稳性。
下面用ADF检验来检验序列的平稳性。
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooadf.test(b6)##
## Augmented Dickey-Fuller Test
##
## data: b6
## Dickey-Fuller = -2.168, Lag order = 4, p-value = 0.5069
## alternative hypothesis: stationary由于p-value>\(\alpha\)=0.05，故不拒绝原假设（原假设认为时间序列是非平稳的），即可认为序列不是平稳的。
for(i in 1:12) print(Box.test(b6,lag=i))##
## Box-Pierce test
##
## data: b6
## X-squared = 17.921, df = 1, p-value = 2.303e-05
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 38.22, df = 2, p-value = 5.019e-09
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 47.99, df = 3, p-value = 2.14e-10
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 53.906, df = 4, p-value = 5.505e-11
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 58.56, df = 5, p-value = 2.41e-11
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 60.084, df = 6, p-value = 4.327e-11
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 65.174, df = 7, p-value = 1.388e-11
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 67.6, df = 8, p-value = 1.474e-11
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 69.808, df = 9, p-value = 1.66e-11
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 74.482, df = 10, p-value = 6e-12
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 77.482, df = 11, p-value = 4.519e-12
##
##
## Box-Pierce test
##
## data: b6
## X-squared = 81.066, df = 12, p-value = 2.583e-12在滞后1-12期中，由于p-value<\(\alpha\)=0.05，故拒绝该序列纯随机的假设，我们可以认为此序列不具有纯随机性。
a6_2=diff(a6)
b6_2=ts(a6_2,frequency=12)
b6_2## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1 5 -5 0 2 -2 -3 0 3 4 -6 9 -3
## 2 4 -15 6 2 -1 -4 6 2 -4 15 4 4
## 3 0 -21 7 -3 3 0 -7 22 -19 21 -7 -3
## 4 -5 -4 2 -6 7 4 -12 -6 8 -8 6 -3
## 5 2 6 -5 -4 6 0 -2 -7 3 2 -7 13
## 6 -8 0 -1 5 -6 4 -2 2 -5acf(b6_2)
adf.test(b6_2)## Warning in adf.test(b6_2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: b6_2
## Dickey-Fuller = -4.3625, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary检验可知，序列是平稳的。
for(i in 1:12) print(Box.test(b6_2,lag=i))##
## Box-Pierce test
##
## data: b6_2
## X-squared = 19.342, df = 1, p-value = 1.093e-05
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 21.965, df = 2, p-value = 1.699e-05
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 22.412, df = 3, p-value = 5.355e-05
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 22.654, df = 4, p-value = 0.0001485
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 23.242, df = 5, p-value = 0.0003035
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 27.755, df = 6, p-value = 0.0001045
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 30.96, df = 7, p-value = 6.324e-05
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 31.345, df = 8, p-value = 0.000122
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 31.682, df = 9, p-value = 0.0002259
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 32.38, df = 10, p-value = 0.000346
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 32.539, df = 11, p-value = 0.0006247
##
##
## Box-Pierce test
##
## data: b6_2
## X-squared = 33.288, df = 12, p-value = 0.0008727在滞后1-12期中，由于p-value<\(\alpha\)=0.05，故拒绝该序列纯随机的假设，我们可以认为此序列不具有纯随机性。