1 平稳时间序列分析：平稳非白噪声序列

arima.sim()：拟合平稳 AR 序列、MA 序列、平稳 ARMA 序列、ARIMA 序列
filter()：拟合 AR 序列、MA 序列

1.1 AR(p) 模型

1.1.1 平稳性判别

\[ 1. \quad x_t=0.8x_{t-1}+\epsilon_t \]

\[ 2. \quad x_t=-1.1x_{t-1}+\epsilon_t\]

\[ 3. \quad x_t=x_{t-1}-0.5x_{t-2}+\epsilon_t\]

\[ 4. \quad x_t=x_{t-1}+0.5x_{t-2}+\epsilon_t\]

1.1.1.1 图示法

x1<-arima.sim(n=100,list(ar=0.8))
x3<-arima.sim(n=100,list(ar=c(1,-0.5)))
e<-rnorm(100)
x2<-filter(e,filter = -1.1,method = "recursive")
x4<-filter(e,filter = c(1,0.5),method = "recursive")
par(mfrow=c(2,2))
ts.plot(x1)
ts.plot(x2)
ts.plot(x3)
ts.plot(x4)

图示可知，AR 模型 1.3. 平稳，2.4. 非平稳。

而图示判别法只是一种粗糙的直观判别方法，需要有准确的平稳性判别方法。

1.1.1.2 特征根判别法

1.1.1.3 平稳域判别法

1.1.2 自相关图

\[ 1. \quad x_t=0.8x_{t-1}+\epsilon_t \]

\[ 2. \quad x_t=-0.8x_{t-1}+\epsilon_t \]

\[ 3. \quad x_t=x_{t-1}-0.5x_{t-2}+\epsilon_t\]

\[ 4. \quad x_t=-x_{t-1}-0.5x_{t-2}+\epsilon_t\]

x1<-arima.sim(n=1000,list(ar=0.8))
x2<-arima.sim(n=1000,list(ar=-0.8))
x3<-arima.sim(n=1000,list(ar=c(1,-0.5)))
x4<-arima.sim(n=1000,list(ar=c(-1,-0.5)))
par(mfrow=c(2,2))
acf(x1)
acf(x2)
acf(x3)
acf(x4)

平稳 AR 模型的自相关系数均呈现拖尾性（ \(\rho_k\) 始终有非零取值，不会在 k 大于某个常熟之后恒等于 0）和呈指数衰减至 0 附近的性质。

AR(p)模型的自相关函数ACF不能在某一步之后为零（截尾），而是按指数衰减（或成正弦波形式)，称其具有拖尾性。

但由于特征根的不同，其自相关系数衰减的方式也不同，例如：

  1. 按负指数单调收敛到 0

  2. 呈现正负相间衰减到 0

  3. 呈现出类似周期性余弦衰减到 0 （“伪周期”）

1.1.3 偏自相关图

par(mfrow=c(2,2))
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)

AR(p)模型的偏自相关函数PACF在p阶之后应为零，称其具有截尾性。

1.2 MA(q) 模型

1.2.1 自相关图

\[ 1. \quad x_t=\epsilon_t-2\epsilon_{t-1}\]

\[ 2. \quad x_t=\epsilon_t-0.5\epsilon_{t-1}\]

\[ 3. \quad x_t=\epsilon_t-\frac{4}{5}\epsilon_{t-1}+\frac{16}{25}\epsilon_{t-2}\]

\[ 4. x_t=\epsilon_t-\frac{5}{4}\epsilon_{t-1}+\frac{25}{16}\epsilon_{t-2}\]

x1<-arima.sim(n=1000,list(ma=-2))
x2<-arima.sim(n=1000,list(ma=-0.5))
x3<-arima.sim(n=1000,list(ma=c(-4/5,16/25)))
x4<-arima.sim(n=1000,list(ma=c(-5/4,25/16)))
par(mfrow=c(2,2))
acf(x1)
acf(x2)
acf(x3)
acf(x4)

MA(q)模型的自相关函数ACF在q阶之后应为零，称其具有截尾性。

再次观察发现 MA 模型 1.2. 均具有相同的自相关系数。因为一个自相关系数未必对应同一个平稳时间序列模型。

此时需要给 MA 模型施加可逆性条件。符合的模型为 2.3.

数学推导可得，MA 模型的可逆概念与 AR 模型的平稳概念完全对偶。

1.2.2 偏自相关图

par(mfrow=c(2,2))
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)

MA(q)模型的偏自相关函数PACF不能在某一步之后为零（截尾），而是按指数衰减（或成正弦波形式)，称其具有拖尾性。

1.3 ARMA(p,q) 模型

拟合 ARMA(1,1) 模型：

\[x_t=0.5x_{t-1}+\epsilon_t-0.8\epsilon_{t-1}\]

x<-arima.sim(n=1000,list(ar=0.5,ma=-0.8))
acf(x)

pacf(x)

ARMA模型自相关函数和偏自相关函数均拖尾。

1.4 平稳时间序列建模

1.4.1 综合应用 1.

这是1950年-2008年我国邮路及农村投递路线每年新增里程数序列

a<-read.table("file8.csv",sep=",",header = T)
x<-ts(a$kilometer,start = 1950)
plot(x)

序列没有显著的非平稳特征。

# 1. 白噪声检验
for(i in 1:2) print(Box.test(x,type = "Ljung-Box",lag=6*i))
## 
##  Box-Ljung test
## 
## data:  x
## X-squared = 37.754, df = 6, p-value = 1.255e-06
## 
## 
##  Box-Ljung test
## 
## data:  x
## X-squared = 44.62, df = 12, p-value = 1.197e-05

白噪声检验显示，序列值之间蕴含着相关关系（p值很小，拒绝原假设），故此为非白噪声序列。

# 2. 绘制自相关图
acf(x)

自相关图同时显示，除lag 1-3 阶在 2 倍标准差范围之外，其余的阶数基本都在 2 倍标准差之内，即该序列具有短期相关性，进一步确定序列平稳。

同时，自相关系数衰减到 0 时，可看到有明显的正弦波动轨迹，说明自相关系数并非突然衰减到 0 ，而是一个连续渐变的过程，这是自相关系数典型的拖尾特征。

# 3. 绘制偏自相关图
pacf(x)

偏自相关系数衰减到 0 的过程中，除lag 1-2 阶偏自相关系数在 2 倍标准差范围外，其余阶偏自相关系数均在 2 倍标准差范围之内，这是偏自相关系数 2 阶截尾的典型特征。

由上分析，可初步判定模型为 AR(2) 模型。

1.4.2 综合应用 2.

这是美国科罗拉多州某加油站连续57天的 Overshort 序列。

overshort<-read.table("file9.csv",sep=",",header = T)
overshort<-ts(overshort$overshort)
plot(overshort)

时间序列图显示，该序列没有显著的非平稳特征。

# 1. 白噪声检验
for(i in 1:2) print(Box.test(overshort,type = "Ljung-Box",lag=6*i))
## 
##  Box-Ljung test
## 
## data:  overshort
## X-squared = 20.239, df = 6, p-value = 0.002511
## 
## 
##  Box-Ljung test
## 
## data:  overshort
## X-squared = 31.366, df = 12, p-value = 0.001732
# 2. 自相关图
acf(overshort)

自相关图同时显示，除lag 1 阶在 2 倍标准差范围之外，其余的阶数基本都在 2 倍标准差之内，即该序列具有短期相关性，进一步确定序列平稳。

同时，该自相关系数 1 阶截尾。

# 3. 偏自相关图
pacf(overshort)

偏自相关图显示非截尾的性质。

由上分析，拟合该模型为 MA(1) 模型。

1.4.3 综合应用 3.

这是1880年-1985年全球气表平均温度改变值差分序列

b<-read.table("file10.csv",sep=",",header = T)
dif_x<-ts(diff(b$change_temp),start = 1880)
plot(dif_x)

# 1. 白噪声检验
for(i in 1:2) print(Box.test(dif_x,type = "Ljung-Box",lag=6*i))
## 
##  Box-Ljung test
## 
## data:  dif_x
## X-squared = 17.436, df = 6, p-value = 0.007807
## 
## 
##  Box-Ljung test
## 
## data:  dif_x
## X-squared = 26.266, df = 12, p-value = 0.009841
# 2. 自相关图
acf(dif_x)

# 3. 偏自相关图
pacf(dif_x)

由上，我们可以初步拟合该序列为 ARMA(1,1) 模型。

1.4.4 综合应用 4.

auto.arima() 函数可自动拟合 ARMA(p,q) 模型，基于 信息量最小原则 自动识别模型阶数，并给出模型的参数估计值。

# 对上面的综合应用 1. 重新定阶
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
auto.arima(x)
## Series: x 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     mean
##       0.7185  -0.5294  11.0226
## s.e.  0.1083   0.1067   3.0906
## 
## sigma^2 estimated as 384.7:  log likelihood=-258.23
## AIC=524.46   AICc=525.2   BIC=532.77
# 对上面的综合应用 2. 重新定阶
auto.arima(overshort)
## Series: overshort 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##           ma1     mean
##       -0.8477  -4.7945
## s.e.   0.1206   1.0252
## 
## sigma^2 estimated as 2093:  log likelihood=-298.42
## AIC=602.84   AICc=603.29   BIC=608.97
# 对上面的综合应用 3. 重新定阶
auto.arima(dif_x)
## Series: dif_x 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    mean
##       0.3926  -0.8867  0.0053
## s.e.  0.1180   0.0604  0.0024
## 
## sigma^2 estimated as 0.01586:  log likelihood=69.66
## AIC=-131.32   AICc=-130.92   BIC=-120.71

1.5 参数估计：arima()

# 对上面的综合应用 1. 进行参数估计
a<-read.table("file8.csv",sep=",",header = T)
x<-ts(a$kilometer,start=1950)
x.fit<-arima(x,order = c(2,0,0),method = "ML") # 极大似然估计
x.fit
## 
## Call:
## arima(x = x, order = c(2, 0, 0), method = "ML")
## 
## Coefficients:
##          ar1      ar2  intercept
##       0.7185  -0.5294    11.0223
## s.e.  0.1083   0.1067     3.0906
## 
## sigma^2 estimated as 365.2:  log likelihood = -258.23,  aic = 524.46

由上，我们可以确定模型：

\[x_t=11.0223+0.7185x_{t-1}-0.5294x_{t-2}+\epsilon_t\]

其中 \(\epsilon_t \sim WN(0,365.2)\) ，AIC 值为 524.46 。

# 对上面的综合应用 2. 进行参数估计
overshort<-read.table("file9.csv",sep=",",header = T)
overshort<-ts(overshort$overshort)
overshort.fit<-arima(overshort,order = c(0,0,1))
overshort.fit
## 
## Call:
## arima(x = overshort, order = c(0, 0, 1))
## 
## Coefficients:
##           ma1  intercept
##       -0.8477    -4.7945
## s.e.   0.1206     1.0252
## 
## sigma^2 estimated as 2020:  log likelihood = -298.42,  aic = 602.84

由上，我们可以确定模型：

\[x_t=-4.7945+\epsilon_t-0.8477\epsilon_{t-1}\]

其中 \(\epsilon_t \sim WN(0,2020)\) ，AIC值为 602.84 。

# 对上面的综合应用 3. 进行参数估计
b<-read.table("file10.csv",sep=",",header = T)
dif_x<-ts(diff(b$change_temp),start = 1880)
dif_x.fit<-arima(dif_x,order = c(1,0,1))
dif_x.fit
## 
## Call:
## arima(x = dif_x, order = c(1, 0, 1))
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.3926  -0.8867     0.0053
## s.e.  0.1180   0.0604     0.0024
## 
## sigma^2 estimated as 0.01541:  log likelihood = 69.66,  aic = -131.32

由上，我们可以确定模型：

\[x_t=0.0053+0.3926x_{t-1}+\epsilon_t-0.8867\epsilon_{t-1}\]

其中 \(\epsilon_t \sim WN(0,0.01541)\) ，AIC 值为 -131.32 。（理论上说 AIC 负的是好的，越负越拟合的好，但是拟合的太好了，模型说服力就不够了，总之 AIC 最小原则是看实际值，不是绝对值）

1.6 模型检验

1.6.1 模型的显著性检验：残差序列的白噪声序列

一个模型是否显著有效，主要是看它提取的信息是否充分。一个好的拟合模型应当能够提取观察值序列中几乎所有的样本相关信息，此时残差应当不含任何信息，为白噪声序列。

\[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\]

# 对上面的综合应用 1. 检验拟合模型的显著性
a<-read.table("file8.csv",sep=",",header = T)
x<-ts(a$kilometer,start=1950)
x.fit<-arima(x,order = c(2,0,0),method = "ML")
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))
## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 2.0949, df = 6, p-value = 0.9108
## 
## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 2.8341, df = 12, p-value = 0.9966

由于各阶延迟下p值较大，故不拒绝原假设，可认为该拟合模型残差序列属于白噪声序列，拟合模型显著有效。

# 对上面的综合应用 2. 检验拟合模型的显著性
overshort<-read.table("file9.csv",sep=",",header = T)
overshort<-ts(overshort$overshort)
overshort.fit<-arima(overshort,order = c(0,0,1))
for(i in 1:2) print(Box.test(overshort.fit$residual,lag=6*i))
## 
##  Box-Pierce test
## 
## data:  overshort.fit$residual
## X-squared = 2.984, df = 6, p-value = 0.8108
## 
## 
##  Box-Pierce test
## 
## data:  overshort.fit$residual
## X-squared = 8.4545, df = 12, p-value = 0.7487

由于各阶延迟下p值较大，故不拒绝原假设，可认为该拟合模型残差序列属于白噪声序列，拟合模型显著有效。

# 对上面的综合应用 3. 检验拟合模型的显著性
b<-read.table("file10.csv",sep=",",header = T)
dif_x<-ts(diff(b$change_temp),start = 1880)
dif_x.fit<-arima(dif_x,order = c(1,0,1),method = "CSS")
for(i in 1:2) print(Box.test(dif_x.fit$residual,lag=6*i))
## 
##  Box-Pierce test
## 
## data:  dif_x.fit$residual
## X-squared = 4.593, df = 6, p-value = 0.597
## 
## 
##  Box-Pierce test
## 
## data:  dif_x.fit$residual
## X-squared = 9.1007, df = 12, p-value = 0.6943

由于各阶延迟下p值较大，故不拒绝原假设，可认为该拟合模型残差序列属于白噪声序列，拟合模型显著有效。

1.6.2 参数的显著性检验：检验模型参数是否显著为0，目的是为模型最精简

R 中并不提供参数的显著性检验的结果，默认输出参数显著非 0 ，需要自行计算p值、t统计量。

# 对上面的综合应用 1. 检验参数的显著性
a<-read.table("file8.csv",sep=",",header = T)
x<-ts(a$kilometer,start=1950)
x.fit<-arima(x,order = c(2,0,0),method = "ML")
#ar1系数显著性检验
t1<-0.7185/0.1083
pt(t1,df=56,lower.tail = F) # lower.tail=F 表明计算 Pr(X > x)。对于参数显著性检验，若参数估计值为正，选择此。
## [1] 6.94276e-09
#ar2系数显著性检验
t2<-0.5294/0.1067
pt(t2,df=56,lower.tail = F)
## [1] 3.43633e-06
#常数数显著性检验
t0=11.0223/3.0906
pt(t0,df=56,lower.tail = F)
## [1] 0.0003748601

检验结果显示，三个系数均显著非 0 。

1.7 模型优化

对于某一个时间序列，可能有两个不同的拟合模型均显著有效，此时需要判断对比模型优劣。

常用的准则为 AIC 信息准则、 BIC 信息准则。

# 问题的引入
# 这是等时间间隔、连续读取某化学反应的70个过程数据
x<-read.table(file = "file11.csv",sep=",",header = T)
x<-ts(x$yield)
plot(x)

# 1. 白噪声检验
for(i in 1:2) print(Box.test(x,lag=6*i))
## 
##  Box-Pierce test
## 
## data:  x
## X-squared = 20.209, df = 6, p-value = 0.002542
## 
## 
##  Box-Pierce test
## 
## data:  x
## X-squared = 21.622, df = 12, p-value = 0.04198
# 2. 自相关图
acf(x)

# 3. 偏自相关图
pacf(x)

# 4. 拟合 MA(2) 模型
x.fit1<-arima(x,order = c(0,0,2))
x.fit1
## 
## Call:
## arima(x = x, order = c(0, 0, 2))
## 
## Coefficients:
##           ma1     ma2  intercept
##       -0.3194  0.3019    51.1695
## s.e.   0.1160  0.1233     1.2516
## 
## sigma^2 estimated as 114.4:  log likelihood = -265.35,  aic = 538.71
# 4.1 MA(2) 模型显著性检验
for(i in 1:2) print(Box.test(x.fit1$residual,lag=6*i))
## 
##  Box-Pierce test
## 
## data:  x.fit1$residual
## X-squared = 2.1105, df = 6, p-value = 0.9093
## 
## 
##  Box-Pierce test
## 
## data:  x.fit1$residual
## X-squared = 3.9217, df = 12, p-value = 0.9848
# 5. 拟合 AR(1) 模型
x.fit2<-arima(x,order = c(1,0,0))
x.fit2
## 
## Call:
## arima(x = x, order = c(1, 0, 0))
## 
## Coefficients:
##           ar1  intercept
##       -0.4191    51.2658
## s.e.   0.1129     0.9137
## 
## sigma^2 estimated as 116.6:  log likelihood = -265.98,  aic = 537.96
# 5.1 AR(1) 模型显著性检验
for(i in 1:2) print(Box.test(x.fit2$residual,lag=6*i))
## 
##  Box-Pierce test
## 
## data:  x.fit2$residual
## X-squared = 4.1678, df = 6, p-value = 0.654
## 
## 
##  Box-Pierce test
## 
## data:  x.fit2$residual
## X-squared = 6.1411, df = 12, p-value = 0.9088

对比可得。后一模型的 AIC 值相对较低，故应当选择 AR(1) 模型。

注意到并没有 BIC 值，我们可以用 ‘BIC(x.fit1)’ 与 ‘BIC(x.fit2)’ 计算两个模型的 BIC 值。

1.8 时间序列预测

# 对上面的综合应用 1. 预测2009年-2013年我国邮路及农村的投递路线每年新增里程数
a<-read.table("file8.csv",sep=",",header = T)
x<-ts(a$kilometer,start=1950)
x.fit<-arima(x,order = c(2,0,0))
x.fore<-forecast(x.fit,h=5)
x.fore
##      Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 2009       9.465515 -15.02434 33.95537 -27.98848 46.91951
## 2010       6.215194 -23.94089 36.37128 -39.90456 52.33495
## 2011       8.392736 -21.76507 38.55054 -37.72964 54.51512
## 2012      11.678108 -19.95476 43.31098 -36.70019 60.05640
## 2013      12.885903 -19.44653 45.21833 -36.56228 62.33409
# 输出预测图
plot(x.fore)

# 个性化输出预测图
L1<-x.fore$fitted-1.96*sqrt(x.fit$sigma2)
U1<-x.fore$fitted+1.96*sqrt(x.fit$sigma2)
L2<-ts(x.fore$lower[,2],start = 2009)
U2<-ts(x.fore$upper[,2],start = 2009)
c1<-min(x,L1,L2)
c2<-max(x,L2,U2)
plot(x,type = "p",pch=8,xlim = c(1950,2013),ylim = c(c1,c2))
lines(x.fore$fitted,col=2,lwd=2)
lines(x.fore$mean,col=2,lwd=2)
lines(L1,col=4,lty=2)
lines(L2,col=4,lty=2)
lines(U1,col=4,lty=2)
lines(U2,col=4,lty=2)

2 Homework

2.1 1.

a1=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/3rd_homework/3.17.txt')
b1=ts(a1,frequency=1)
ts.plot(b1)

acf(b1,lag=20)

通过acf图，我们可认为序列是平稳的。

for(i in 1:12) print(Box.test(b1,lag=i))
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 5.9049, df = 1, p-value = 0.0151
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 11.417, df = 2, p-value = 0.003318
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 11.505, df = 3, p-value = 0.009288
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 12.516, df = 4, p-value = 0.0139
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 12.521, df = 5, p-value = 0.02831
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 12.525, df = 6, p-value = 0.05122
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 12.839, df = 7, p-value = 0.07613
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 13.112, df = 8, p-value = 0.108
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 13.268, df = 9, p-value = 0.1509
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 13.282, df = 10, p-value = 0.2083
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 13.485, df = 11, p-value = 0.2628
## 
## 
##  Box-Pierce test
## 
## data:  b1
## X-squared = 13.849, df = 12, p-value = 0.3105

滞后6期以前，由于p-value<\(\alpha\)=0.05，故拒绝该序列纯随机的假设，我们可以认为此序列不具有纯随机性。

综上，我们认为序列平稳且非白噪声。

pacf(b1,lag=20)

根据自相关图拖尾、偏自相关图1阶截尾的特性，我们选择\(AR(1)\)模型进行模拟。

library(zoo)
library(forecast)
b1_fit=arima(b1,order=c(1,0,0),method='ML')
b1_fore=forecast(b1_fit,h=5)
b1_fore
##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 64       90.49459 61.93121 119.0580 46.81067 134.1785
## 65       84.05415 53.97399 114.1343 38.05052 130.0578
## 66       81.92759 51.68667 112.1685 35.67809 128.1771
## 67       81.22543 50.96703 111.4838 34.94920 127.5017
## 68       80.99358 50.73328 111.2539 34.71445 127.2727

以上即为5年预测的降雪量。

2.2 2.

a2=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/3rd_homework/3.18.txt')
b2=ts(a2,frequency=1)
ts.plot(b2)

acf(b2,lag=20)

通过acf图，我们可认为序列是平稳的。

for(i in 1:12) print(Box.test(b2,lag=i))
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 9.7736, df = 1, p-value = 0.00177
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 15.017, df = 2, p-value = 0.0005484
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 18.819, df = 3, p-value = 0.0002981
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 22.152, df = 4, p-value = 0.0001869
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 24.8, df = 5, p-value = 0.0001523
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 27.983, df = 6, p-value = 9.466e-05
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 29.543, df = 7, p-value = 0.0001151
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 29.642, df = 8, p-value = 0.0002445
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 31.521, df = 9, p-value = 0.0002409
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 33.492, df = 10, p-value = 0.0002251
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 34.609, df = 11, p-value = 0.0002874
## 
## 
##  Box-Pierce test
## 
## data:  b2
## X-squared = 35.375, df = 12, p-value = 0.0004081

由于p-value<\(\alpha\)=0.05，故拒绝该序列纯随机的假设，我们可以认为此序列不具有纯随机性。
综上，我们认为序列平稳且非白噪声。

pacf(b2,lag=20)

根据自相关图拖尾、偏自相关图1阶截尾的特性，我们选择\(AR(1)\)模型进行模拟。

b2_fit=arima(b2,order=c(1,0,0),method='ML')
b2_fore=forecast(b2_fit,h=5)
b2_fore
##    Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 75      0.7058602 0.3556749 1.056046 0.1702980 1.241422
## 76      0.7963705 0.4232095 1.169532 0.2256700 1.367071
## 77      0.8296908 0.4535239 1.205858 0.2543932 1.404988
## 78      0.8419572 0.4653848 1.218530 0.2660394 1.417875
## 79      0.8464730 0.4698456 1.223100 0.2704712 1.422475

以上即为5年预测的谷物产量。

2.3 3.

a3=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/3rd_homework/3.19.txt')
b3=ts(a3,frequency=1)
ts.plot(b3)

acf(b3,lag=20)

通过acf图，我们可认为序列是平稳的。

for(i in 1:12) print(Box.test(b3,lag=i))
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 17.879, df = 1, p-value = 2.354e-05
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 20.795, df = 2, p-value = 3.05e-05
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 21.255, df = 3, p-value = 9.32e-05
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 23.588, df = 4, p-value = 9.659e-05
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 25.392, df = 5, p-value = 0.000117
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 30.403, df = 6, p-value = 3.294e-05
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 31.093, df = 7, p-value = 5.977e-05
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 31.133, df = 8, p-value = 0.000133
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 31.909, df = 9, p-value = 0.0002065
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 32.34, df = 10, p-value = 0.0003514
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 32.34, df = 11, p-value = 0.0006725
## 
## 
##  Box-Pierce test
## 
## data:  b3
## X-squared = 33.119, df = 12, p-value = 0.0009273

由于p-value<\(\alpha\)=0.05，故拒绝该序列纯随机的假设，我们可以认为此序列不具有纯随机性。

综上，我们认为序列平稳且非白噪声。

pacf(b3,lag=20)

acf图1阶截尾，pacf显示出非截尾的性质，故我们选择\(MA(1)\)模型拟合。

b3_fit=arima(b3,order=c(0,0,1),method='ML')
b3_fore=forecast(b3_fit,h=1)
b3_fore
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 202       85.68041 82.25054 89.11027 80.43488 90.92593

上述显示了序列下一刻95%的置信区间。

2.4 4.

a4=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/3rd_homework/3.20.txt')
b4=ts(a4,frequency=1)
ts.plot(b4)

acf(b4,lag=20)

通过acf图，我们可认为序列是平稳的。

for(i in 1:12) print(Box.test(b4,lag=i))
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 0.79671, df = 1, p-value = 0.3721
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 1.1663, df = 2, p-value = 0.5581
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 11.705, df = 3, p-value = 0.008464
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 15.397, df = 4, p-value = 0.003945
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 16.378, df = 5, p-value = 0.005844
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 16.805, df = 6, p-value = 0.01003
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 16.89, df = 7, p-value = 0.01812
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 21.301, df = 8, p-value = 0.006389
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 21.305, df = 9, p-value = 0.01136
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 21.417, df = 10, p-value = 0.01837
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 21.432, df = 11, p-value = 0.02916
## 
## 
##  Box-Pierce test
## 
## data:  b4
## X-squared = 22.453, df = 12, p-value = 0.03274

由于p-value<\(\alpha\)=0.05，故拒绝该序列纯随机的假设，我们可以认为此序列不具有纯随机性。

综上，我们认为序列平稳且非白噪声。

pacf(b4,lag=20)

易知acf与pacf均显示出不截尾的性质，故我们选择\(ARMA(1,3)\)拟合模型。

b4_fit=arima(b4,order=c(1,0,3),method='ML')
b4_fore=forecast(b4_fit,h=5)
b4_fore
##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 89       63.65994 41.04997 86.26991 29.08097 98.23890
## 90       57.74120 34.29004 81.19237 21.87575 93.60666
## 91       48.18084 24.47474 71.88695 11.92549 84.43620
## 92       50.92139 25.31851 76.52427 11.76516 90.07762
## 93       51.93917 26.08560 77.79274 12.39954 91.47880
plot(b4_fore)

L1<-b4_fore$fitted-1.96*sqrt(b4_fit$sigma2)
U1<-b4_fore$fitted+1.96*sqrt(b4_fit$sigma2)
L2<-ts(b4_fore$lower[,2],start=89)
U2<-ts(b4_fore$upper[,2],start=89)
c1<-min(b4,L1,L2)
c2<-max(b4,L2,U2)
plot(b4,type="p",pch=8,xlim=c(1,88),ylim=c(c1,c2))
lines(b4_fore$fitted,col=2,lwd=2)
lines(b4_fore$mean,col=2,lwd=2)
lines(L1,col=4,lty=2)
lines(U1,col=4,lty=2)
lines(L1,col=4,lty=2)
lines(L2,col=4,lty=2)
lines(U2,col=4,lty=2)

星号为观察值序列，实线为拟合值，虚线为置信水平95%的置信水平线。

实际上，根据auto.arima()函数的计算，这些模型都不是最优模型，例如：

  • 1 中应当建模为 ARIMA(0,1,1)
  • 2 中应当建模为 ARIMA(0,1,1) with drift
  • 3 中应当建模为 ARIMA(0,0,2) with non-zero mean
  • 4 中应当建模为 ARIMA(2,1,3)