Time Series Analysis Chapter 5

1 非平稳时间序列的随机分析

1.1 差分运算

1.1.1 序列蕴含显著线性趋势，1 阶差分实现趋势平稳

# 这是1964年-1999年我国纱年产量序列，其蕴含近似线性的趋势，对其进行 1 阶差分运算
a<-read.table("file4.csv",sep=",",header = T)
x<-ts(a$output,start = 1964)
x.dif<-diff(x)
plot(x.dif)

时序图显示，1 阶差分运算已从原序列中提取出线性趋势，差分后的序列呈现出非常平稳的随机波动。

1.1.2 序列蕴含曲线趋势，低阶差分提取曲线趋势影响

# 这是1950年-1999年北京市民用车辆拥有量序列
b<-read.table("file15.csv",sep=",",header = T)
x<-ts(b$vehicle,start = 1950)
plot(x)

时序图中我们看到该序列蕴含曲线递增的长期趋势。

# 1 阶差分
x.dif<-diff(x)
# 1 阶差分后序列的时序图
plot(x.dif)

1 阶差分序列时序图表示 1 阶差分提取出原序列的部分长期趋势，但长期趋势的信息提取并不充分，1 阶差分后序列中仍然蕴含长期递增的趋势。

# 2 阶差分
x.dif2<-diff(x,1,2)
# 2 阶差分后序列的时序图
plot(x.dif2)

2 阶差分较为充分地提取出原序列的长期趋势，使得差分后序列不再呈现递增的确定性趋势。

1.1.3 蕴含固定周期的序列

# 这是2001年1月至2014年11月中国月度出口额序列，利用差分运算提取确定性信息
c<-read.table("file16.csv",sep=",",header = T)
x<-ts(c$export,start = c(2001,1),frequency = 12)
plot(x)

时序图显示，该序列具有 一个线性递增的长期趋势 和 一个周期长度为一年的稳定的季节变动。

对原序列先做一阶差分，提取线性递增趋势。

# 1 阶差分，并绘制出差分后序列的时序图
x.dif<-diff(x)
plot(x.dif)

1 阶差分后的序列具有稳定的季节波动和随机波动。

对 1 阶差分后的序列再进行 12 步的周期差分，提取季节波动信息。

#1阶差分加12步差分，并绘制出差分后序列的时序图
x.dif1_12<-diff(diff(x),12)
plot(x.dif1_12)

周期差分可非常好地提取周期信息。此时，差分运算已经较为充分地提取出原序列中蕴含的季节效应和长期确定性趋势等确定性信息。

1.2 ARIMA(p,d,q) 模型

1.2.1 随机游走序列：ARIMA(0,1,0)

\[x_t=x_{t-1}+\epsilon_t,\epsilon_t \sim WN(0,Var(\epsilon_t))\]

# 拟合随机游走序列，其中残差标准差为 10
x<-arima.sim(n=1000,list(order=c(0,1,0)),sd=10)
plot(x)

时序图显示该序列显然非平稳。

1.2.2 ARIMA 模型建模

# 这是1952年-1988年中国农业实际国民收入指数序列
d<-read.table("file17.csv",sep=",",header = T)
x<-ts(d$index,start = 1952)
plot(x)

时序图显示该序列具有线性趋势。

#1阶差分，并绘制出差分后序列的时序图
x.dif<-diff(x)
plot(x.dif)

1 阶差分时序图显示，该序列在均值附近较稳定的波动。

考察该 1 阶差分序列图地平稳性。

# 绘制差分后序列自相关图
acf(x.dif)

1 阶差分序列自相关图显示，除了 lag 1 阶自相关系数显著非 0 外，其他阶数自相关系数均落在 2 倍标准差地范围之内，显示出很强的短期相关性，可认为 1 阶差分后序列平稳。

# 绘制差分后序列偏自相关图
pacf(x.dif)

由自相关图、偏自相关图显示，自相关系数 1 阶截尾，偏自相关系数拖尾。故对原序列拟合 ARIMA(0,1,1) 模型。

# 拟合ARIMA(0,1,1)模型
x.fit<-arima(x,order=c(0,1,1))
x.fit

## 
## Call:
## arima(x = x, order = c(0, 1, 1))
## 
## Coefficients:
##          ma1
##       0.7355
## s.e.  0.1545
## 
## sigma^2 estimated as 61.95:  log likelihood = -125.74,  aic = 255.49

拟合的模型为

\[x_t=x_{t-1}+\epsilon_t+0.7355 \epsilon_{t-1},\epsilon_t \sim WN(0,61.95)\]

# 残差白噪声检验
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 3.3169, df = 6, p-value = 0.7681
## 
## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 6.0284, df = 12, p-value = 0.9146

残差序列白噪声检验说明该模型显著成立，说明此 ARIMA(0,1,1) 模型对该序列拟合成功。

1.2.3 ARIMA 模型预测

# 对1952年-1988年我国农业实际国民收入指数序列作为期 10 年的预测
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

d<-read.table("file17.csv",sep=",",header = T)
x<-ts(d$index,start = 1952)
x.fit<-arima(x,order=c(0,1,1))
x.fore<-forecast(x.fit,h=10)
plot(x.fore)

1.2.4 ARIMA 疏系数模型

模型中有部分参数省缺的模型。R 中用arima() 函数拟合。

• ’transform.parts’参数置顶参数估计是否由系统自动完成。’transform.parts=T’表明系统默认根据’order’参数设置模型阶数自动完成参数估计，’transform.parts=F’表明需要拟合疏系数模型，不能让系统根据模型的最高阶数自动完成参数估计，需要人为干预。

• ’fixed’参数说明对疏系数模型指定疏系数的位置。

# 这是1917年-1975年美国 23 岁妇女每万人生育率序列
e<-read.table("file18.csv",sep=",",header = T)
x<-ts(e$fertility,start = 1917)
plot(x)

时序图可知，序列并非平稳。

# 1 阶差分，并绘制出差分后序列的时序图
x.dif<-diff(x)
plot(x.dif)

1 阶差分后的序列没有明显的趋势和周期特征。

# 绘制差分后序列自相关图
acf(x.dif)

自相关图在 5 阶后均落在 2 倍标准差范围之内，满足短期自相关属性，可认为 1 阶差分后序列平稳。

# 绘制差分后序列偏自相关图
pacf(x.dif)

由自相关图、偏自相关图，均呈现截尾的特性，可尝试拟合多种模型。

由偏自相关系数 1 阶和 4 阶显著非 0 ，4 阶后截尾，构造 1 阶差分后的 AR(4) 模型，且4 个参数中 2、3 参数恒为 0 ，即拟合 ARIMA((1,4),1,0) 模型。

# 拟合疏系数模型ARIMA((1,4),1,0)
x.fit<-arima(x,order = c(4,1,0),transform.pars = F,fixed = c(NA,0,0,NA))
x.fit

## 
## Call:
## arima(x = x, order = c(4, 1, 0), transform.pars = F, fixed = c(NA, 0, 0, NA))
## 
## Coefficients:
##          ar1  ar2  ar3     ar4
##       0.2583    0    0  0.3408
## s.e.  0.1159    0    0  0.1225
## 
## sigma^2 estimated as 118.2:  log likelihood = -221,  aic = 448.01

# 残差白噪声检验
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 4.0909, df = 6, p-value = 0.6644
## 
## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 5.3742, df = 12, p-value = 0.9443

白噪声检验结果显示，该疏系数模型显著成立。

# 做5期预测
x.fore<-forecast(x.fit,h=5)
x.fore

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1976      109.18385 95.25187 123.1158 87.87673 130.4910
## 1977      103.43765 81.04480 125.8305 69.19075 137.6846
## 1978      100.93090 71.90947 129.9523 56.54645 145.3154
## 1979       98.64765 64.12267 133.1726 45.84624 151.4491
## 1980       94.88318 53.11326 136.6531 31.00160 158.7648

# 绘制预测图
plot(x.fore)

1.2.5 季节模型：ARIMA 模型可对具有季节效应的序列建模。

1.2.5.1 简单季节模型：季节效应与其他效应之间为加法关系。

R 中arima() 函数中 ’seasonal’指定季节模型阶数、周期。

# 这是1962年-1991年德国工人季度失业率序列
f<-read.table("file19.csv",sep=",",header = T)
x<-ts(f$unemployment_rate,start = c(1962,1),frequency = 4)
plot(x)

时序图显示此序列具有（确定性）趋势和周期，故进行 1 阶差分提取趋势信息，进行 4 阶差分提取周期信息。

# 1 阶 4 步差分，并绘制出差分后序列的时序图
x.dif<-diff(diff(x),4)
plot(x.dif)

差分后的时序图显示，差分后的序列平稳。

# 绘制差分后序列自相关图和偏自相关图
par(mfrow=c(2,1))
acf(x.dif)
pacf(x.dif)

由自相关图、偏自相关图，自相关系数拖尾，偏自相关系数 1、4 阶显著非 0 ，4 阶之后截尾。故拟合加法季节模型 ARIMA((1,4),(1,4),0) 。

# 拟合加法季节模型ARIMA((1,4),(1,4),0)
x.fit<-arima(x,order = c(4,1,0),seasonal = list(order=c(0,1,0),period=4),transform.pars=F,fixed = c(NA,0,0,NA))
x.fit

## 
## Call:
## arima(x = x, order = c(4, 1, 0), seasonal = list(order = c(0, 1, 0), period = 4), 
##     transform.pars = F, fixed = c(NA, 0, 0, NA))
## 
## Coefficients:
##          ar1  ar2  ar3      ar4
##       0.4449    0    0  -0.2720
## s.e.  0.0807    0    0   0.0804
## 
## sigma^2 estimated as 0.09266:  log likelihood = -26.7,  aic = 59.39

这里，得到的拟合模型 存在疑惑，需要解决 。

\[\nabla(\nabla_4 x_t)= 0.4449x_{t-1}-0.2720x_{t-4}+\epsilon_t,\epsilon_t \sim WN(0,0.09266)\]

# 残差白噪声检验
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 2.1452, df = 6, p-value = 0.9059
## 
## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 10.546, df = 12, p-value = 0.5682

残差序列白噪声检验表明该模型拟合良好，对序列相关信息提取充分。

# 做3年期预测
x.fore<-forecast(x.fit,h=12)
x.fore

##         Point Forecast    Lo 80    Hi 80      Lo 95     Hi 95
## 1992 Q1       6.615077 6.224975 7.005178  6.0184679  7.211686
## 1992 Q2       5.856360 5.170875 6.541845  4.8080015  6.904719
## 1992 Q3       5.965918 5.027511 6.904325  4.5307483  7.401088
## 1992 Q4       5.844074 4.687986 7.000161  4.0759905  7.612157
## 1993 Q1       6.146436 4.634425 7.658448  3.8340140  8.458859
## 1993 Q2       5.326343 3.479538 7.173148  2.5018984  8.150787
## 1993 Q3       5.433197 3.294966 7.571428  2.1630545  8.703339
## 1993 Q4       5.343294 2.949603 7.736985  1.6824590  9.004128
## 1994 Q1       5.690528 2.925569 8.455487  1.4618876  9.919168
## 1994 Q2       4.907094 1.762628 8.051559  0.0980490  9.716138
## 1994 Q3       5.030993 1.528164 8.533821 -0.3261208 10.388106
## 1994 Q4       4.939984 1.103738 8.776231 -0.9270484 10.807017

# 绘制预测图
plot(x.fore)

1.2.5.2 乘积季节模型

# 这是1948年-1981年美国女性 20 岁以上月度失业率序列
g<-read.table("file20.csv",sep=",",header = T)
x<-ts(g$unemployment_rate,start = c(1948,1),frequency = 12)
plot(x)

时序图显示该序列既有趋势也有周期特征。

# 作1阶12步差分，并绘制出差分后序列的时序图
x.dif<-diff(diff(x),12)
plot(x.dif)

1 阶 12 步差分后序列呈现平稳特征。

# 绘制差分后序列自相关图和偏自相关图
par(mfrow=c(2,1))
acf(x.dif)
pacf(x.dif)

差分后的序列自相关图、偏自相关图均呈现拖尾属性。

首先拟合加法季节模型 ARIMA(1,(1,12),1) 。

# 拟合ARIMA(1,(1,12),1)模型
x.fit<-arima(x,order = c(1,1,1),seasonal = list(order=c(0,1,0),period=12))
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 11.204, df = 6, p-value = 0.08228
## 
## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 105.78, df = 12, p-value < 2.2e-16

残差序列的白噪声检验显示，该模型拟合效果并不理想，加法季节模型并不适合此序列。

如下考虑乘法季节模型 \(ARIMA(1,1,1)*ARIMA(0,1,1)_{12}\) 。

由于短期相关性与季节效应具有乘积关系，实际上拟合的模型为 ARMA(p,q) 与 ARMA(P,Q) 的乘积。前者表示趋势关系，后者表示季节关系。

首先考虑 1 阶 12 步差分后，序列 12 阶以内的自相关系数和偏自相关系数特征，以确定短期相关模型。 ACF 和 PACF 显示两者均不截尾，故用 ARMA(1,1) 模型拟合。

其次考虑季节自相关特征，考察 lag 12 阶、24 阶等以周期长度为单位的自相关系数和偏自相关系数特征。 ACF 显示 lag 12 阶自相关系数显著非 0 ，但是 lag 24 阶自相关系数落入 2 倍标准差范围之内；PACF 显示 lag 12、24 阶偏自相关系数显著非 0 ，故可认为季节自相关特征为 自相关系数截尾，偏自相关系数拖尾。此时，以 12 步为周期的 \(ARMA(0,1)_{12}\) 模型提取差分后序列的季节自相关信息。

# 拟合ARIMA(1,1,1)*ARIMA(0,1,1)_{12}模型
x.fit<-arima(x,order = c(1,1,1),seasonal = list(order=c(0,1,1),period=12))
x.fit

## 
## Call:
## arima(x = x, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##           ar1     ma1     sma1
##       -0.7290  0.6059  -0.7918
## s.e.   0.1497  0.1728   0.0337
## 
## sigma^2 estimated as 7444:  log likelihood = -2327.14,  aic = 4662.28

# 残差序列白噪声检验
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 4.5564, df = 6, p-value = 0.6018
## 
## 
##  Box-Pierce test
## 
## data:  x.fit$residual
## X-squared = 9.6288, df = 12, p-value = 0.6485

残差序列白噪声检验显示该拟合模型显著成立。

1.3 残差自回归模型：Auto-regressive 模型

Auto-regressive 模型解决差分方法很难对模型进行直观解释的问题，确定性因素分解能够对各种确定性效应解释，但是它对残差信息存在“浪费”。

# 这是1952年-1988年中国农业实际国民收入的指数序列，对其进行残差自回归模型分析。
# 该序列有显著的线性递增趋势，但是没有季节效应。
d<-read.table("file17.csv",sep=",",header = T)
x<-ts(d$index,start = 1952)
# 首先对趋势效应拟合线性回归模型：T_t=\beta_0+\beta_1 t,t=1,2,3...
t<-c(1:37)
x.fit1<-lm(x~t)
summary(x.fit1)

## 
## Call:
## lm(formula = x ~ t)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -28.71 -20.48 -10.81  26.42  46.17 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  66.1491     8.1197   8.147 1.35e-09 ***
## t             4.5158     0.3726  12.121 4.40e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.2 on 35 degrees of freedom
## Multiple R-squared:  0.8076, Adjusted R-squared:  0.8021 
## F-statistic: 146.9 on 1 and 35 DF,  p-value: 4.404e-14

# 其次尝试对趋势效应拟合关于延迟变量的自回归模型：T_t=\beta_0+\beta_1 x_{t-1},t=1,2,3...
xlag<-x[2:37]
x2<-x[1:36]
x.fit2<-lm(x2~xlag)
summary(x.fit2)

## 
## Call:
## lm(formula = x2 ~ xlag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.764  -5.066  -0.703   5.539  20.424 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.39226    4.00444   1.846   0.0736 .  
## xlag         0.91932    0.02464  37.309   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.936 on 34 degrees of freedom
## Multiple R-squared:  0.9762, Adjusted R-squared:  0.9755 
## F-statistic:  1392 on 1 and 34 DF,  p-value: < 2.2e-16

# 对比两个趋势拟合模型：拟合效果图
fit1<-ts(x.fit1$fitted.value,start = 1952)
fit2<-ts(x.fit2$fitted.value,start = 1952)
plot(x,type = "p",pch=8)
lines(fit1,col=2)
lines(fit2,col=4)

由此，我们可以确定两种确定性趋势拟合模型：

（1）

\[x_t=66.1931+4.5158t+\epsilon_t,\epsilon_t \sim WN(0,24.2^2)\]

（2）

\[x_t=7.3923+0.9193 x_{t-1}+\epsilon_t,\epsilon_t \sim WN(0,7.936^2)\]

1.3.1 残差自相关检验：对模型的拟合效果做检验：Durbin-Watson 检验

# 检验第一个确定性趋势模型
library(lmtest)
dwtest(x.fit1)

## 
##  Durbin-Watson test
## 
## data:  x.fit1
## DW = 0.13784, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

DW 统计量值小于 2 ，且 p 值极小，说明残差序列高度正相关。故为充分提取相关信息，需要对残差序列再次进行信息提取。

# 检验第二个确定性趋势模型
# 修正DW统计量以用于检验延迟因变量
dwtest(x.fit2,order.by=xlag) # order.by 指定延迟因变量

## 
##  Durbin-Watson test
## 
## data:  x.fit2
## DW = 1.8153, p-value = 0.2295
## alternative hypothesis: true autocorrelation is greater than 0

应用Durbin h检验显示，残差序列不存在自相关，故不需要对第二个趋势拟合模型的残差进行再度拟合。

1.3.2 残差自相关模型的拟合

# 对第一个确定性拟合模型的残差序列拟合自相关模型
# 绘制差分后序列自相关图和偏自相关图
par(mfrow=c(2,1))
acf(x.fit1$residual)
pacf(x.fit1$residual)

由残差序列自相关图、偏自相关图，可判断自相关系数拖尾、偏自相关系数 2 阶截尾，故拟合 AR(2) 模型。

# 残差拟合 AR(2) 模型
r.fit<-arima(x.fit1$residual,order=c(2,0,0),include.mean = F)
r.fit

## 
## Call:
## arima(x = x.fit1$residual, order = c(2, 0, 0), include.mean = F)
## 
## Coefficients:
##          ar1      ar2
##       1.4995  -0.6028
## s.e.  0.1274   0.1356
## 
## sigma^2 estimated as 50.6:  log likelihood = -126.59,  aic = 259.17

由上，拟合结果为

\[\epsilon_t=-1.4995 \epsilon_{t-1}+0.6028 \epsilon_{t-2}+e_t,e_t \sim WN(0,50.6)\]

模型显著性检验显示，该拟合模型显著成立。

综合上述分析，对原序列，使用如下残差自回归模型进行拟合：

\[\begin{equation} \begin{cases} x_t = 66.1491+4.5158t+\epsilon_t \\ \epsilon_t =-1.4995\epsilon_{t-1}+0.6028\epsilon_{t-2}+e_t,e_t \sim WN(0,50.6) \end{cases} \end{equation}\]

# 残差自相关模型的显著性检验
for(i in 1:2) print(Box.test(r.fit$residual,lag=6*i))

## 
##  Box-Pierce test
## 
## data:  r.fit$residual
## X-squared = 3.1979, df = 6, p-value = 0.7836
## 
## 
##  Box-Pierce test
## 
## data:  r.fit$residual
## X-squared = 10.661, df = 12, p-value = 0.5582

而该模型相对而言更容易直观解释序列的波动规律：

• 该序列有一个长期的线性递增趋势（斜率 4.5158）。

• 该序列受诸多随机因素的影响，导致随机波动序列具有短期的自相关性。

1.4 异方差性质

忽视异方差的存在会导致残差的方差被严重低估，继而参数的显著性假设容易犯采纳伪错误，从而使得参数的显著性检验失去意义，最终导致模型的拟合精度受到影响。

1.4.1 直观判断

# 这是1963年4月-1971年7月美国短期国库券的月度收益率序列，检验其方差齐性
h<-read.table("file21.csv",sep=",",header = T)
x<-ts(h$yield_rate,start = c(1963,4),frequency = 12)
plot(x)

时序图显示该序列显著非平稳。

# 1阶差分，并绘制出差分后残差序列时序图
x.dif<-diff(x)
plot(x.dif)

1 阶差分后序列时序图显示该序列均值平稳但方差递增的性质。

# 1阶差分后残差平方图
plot(x.dif^2)

1 阶差分后残差平方图显示其具有明显的异方差的特征。

当残差序列异方差时，我们处理的思路有如下方式：

• 方差齐性变换：若已知异方差函数的具体形式

• 拟合条件异方差模型：若未知异方差函数的具体形式。

1.4.2 方差齐性变换

对标准差与水平呈正比的异方差序列，采取对数变换就可以有效实现方差齐性，而这一点常用于金融时间序列。

\[\sigma_t^2=h(\mu_t) \Rightarrow Var[g(x_t)]=\sigma^2 \Rightarrow h(\mu_t)=\mu_t^2, g(\mu_t)=\ln (\mu_t)\]

# 对上述数据使用方差齐性变换方法进行分析
# 作对数变换，并绘制对数变换后时序图
lnx<-log(x) # 对数变换
plot(lnx)

# 1阶差分，并绘制出差分后序列时序图
dif.lnx<-diff(lnx)
plot(dif.lnx)

# 残差白噪声检验
for(i in 1:2) print(Box.test(dif.lnx,lag=6*i))

## 
##  Box-Pierce test
## 
## data:  dif.lnx
## X-squared = 3.4118, df = 6, p-value = 0.7557
## 
## 
##  Box-Pierce test
## 
## data:  dif.lnx
## X-squared = 9.8323, df = 12, p-value = 0.6307

最终的拟合模型为

\[\nabla \ln (x_t)=\epsilon_t\]

其中 \(\epsilon_t\) 为零均值白噪声序列。

1.4.3 条件异方差模型

1.4.3.1 ARCH模型

1.4.3.1.1 集群效应

消除确定性非平稳因素的影响后，残差序列的波动在某些时段持续偏大，某些时段波动持续偏小，呈现出集群效应。

# 这是1926年-1991年标准普尔500股票价值甲醛月度收益率序列，考察其集群效应特征
k<-read.table("file22.csv",sep=",",header = T)
x<-ts(k$returns,start = c(1926,1),frequency = 12)
plot(x)

时序图显示该序列并没有明显的非平稳特征，序列围绕在 0 值波动，但在一些特殊时段（1930年、1940年、1975年、1990年），序列的波动很大，这就是集群效应的特征。

# 绘制序列平方图
plot(x^2)

集群效应意味着在整个序列观察期，序列的方差基本是齐性的，但在某段时间方差显著异于期望方差。这种情况尤其对于资产持有者而言，他们并不关心资产收益率在所有时段的综合表现，只关心在他们持有资产的这段时间，资产收益率是否有较大波动。

1.4.3.1.2 ARCH 模型结构

假设历史数据已知，残差平方序列的自相关系数

\[\rho_k=\frac{Cov(\epsilon_t^2,\epsilon_{t-k}^2)}{Var(\epsilon_t)}\]

当其自相关系数存在某个不为 0 时，意味着残差平方序列蕴含某种相关信息，构建 ARCH(q) q 阶自回归条件异方差模型以获取异方差波动特征。

\[h_t=E(\epsilon_t^2)=\omega+\sum_{j=1}^q \lambda_j \epsilon_{t-j}^2\]

当其自相关系数恒为 0 时，历史数据对未来异方差的估计毫无作用，至今仍无有效方法对其进行分析。

1.4.3.1.3 ARCH 检验：LM 拉格朗日乘子检验、Portmanteau Q 检验

# 这是1926年-1991年标准普尔500股票价值加权月度收益率序列，对其进行 ARCH 检验，并拟合该序列的波动特征
library(tseries)
library(FinTS)

## 
## Attaching package: 'FinTS'

## The following object is masked from 'package:forecast':
## 
##     Acf

for(i in 1:5) print(ArchTest(x,lag=i)) # LM检验

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  x
## Chi-squared = 55.521, df = 1, p-value = 9.244e-14
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  x
## Chi-squared = 68.292, df = 2, p-value = 1.481e-15
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  x
## Chi-squared = 87.129, df = 3, p-value < 2.2e-16
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  x
## Chi-squared = 87.282, df = 4, p-value < 2.2e-16
## 
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  x
## Chi-squared = 87.492, df = 5, p-value < 2.2e-16

for(i in 1:2) print(Box.test(x^2,lag=i)) # Portmanteau Q检验

## 
##  Box-Pierce test
## 
## data:  x^2
## X-squared = 55.545, df = 1, p-value = 9.137e-14
## 
## 
##  Box-Pierce test
## 
## data:  x^2
## X-squared = 85.09, df = 2, p-value < 2.2e-16

两种检验均显示该序列显著方差非齐性，且残差平方序列具有显著自相关关系，用 ARCH 模型提取残差平方序列中蕴含的相关信息。

同时两种检验均显示 1 阶至 5 阶 ARCH 模型均显著成立，即残差平方序列具有长期相关。

拟合 ARCH(5)、ARCH(4) 模型后，其参数检验结果显示有参数不显著，拟合 ARCH(3) 模型尝试。

# 拟合ARCH(3)模型
x.fit<-garch(x,order=c(0,3))

## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     2.904816e-03     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
##      4     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1 -1.904e+03
##      1    5 -1.909e+03  2.72e-03  1.20e-02  1.0e-02  2.3e+07  1.0e-03  1.37e+05
##      2    6 -1.912e+03  1.68e-03  3.42e-03  6.4e-03  2.0e+00  1.0e-03  4.09e+01
##      3    7 -1.913e+03  6.03e-04  5.11e-04  7.3e-03  2.0e+00  1.0e-03  4.33e+01
##      4   11 -1.930e+03  8.56e-03  1.17e-02  2.9e-01  2.0e+00  5.6e-02  4.25e+01
##      5   15 -1.930e+03  3.11e-05  2.99e-04  7.0e-04  7.7e+00  1.3e-04  2.58e-01
##      6   16 -1.930e+03  5.76e-05  5.35e-05  5.3e-04  2.0e+00  1.3e-04  2.42e-01
##      7   21 -1.930e+03  1.81e-04  2.69e-04  7.6e-02  2.0e+00  1.6e-02  2.40e-01
##      8   23 -1.930e+03  7.75e-06  1.49e-04  4.7e-02  2.0e+00  1.2e-02  2.13e-02
##      9   24 -1.931e+03  2.01e-04  2.22e-04  2.1e-02  2.0e+00  6.0e-03  3.53e-02
##     10   25 -1.931e+03  8.77e-05  1.22e-04  2.2e-02  2.0e+00  6.0e-03  1.15e-01
##     11   26 -1.931e+03  1.78e-04  2.33e-04  2.4e-02  2.0e+00  6.0e-03  2.20e-02
##     12   29 -1.935e+03  1.75e-03  3.45e-03  2.1e-01  1.8e+00  5.9e-02  4.92e-02
##     13   30 -1.942e+03  3.85e-03  5.60e-03  1.8e-01  1.2e+00  5.9e-02  1.54e-02
##     14   31 -1.945e+03  1.61e-03  1.22e-03  1.2e-01  1.9e-02  5.9e-02  1.22e-03
##     15   32 -1.946e+03  6.13e-04  4.70e-04  9.4e-02  0.0e+00  5.1e-02  4.70e-04
##     16   33 -1.947e+03  1.02e-04  8.61e-05  4.3e-02  0.0e+00  2.7e-02  8.61e-05
##     17   34 -1.947e+03  5.54e-06  4.97e-06  1.2e-02  0.0e+00  7.4e-03  4.97e-06
##     18   35 -1.947e+03  3.35e-07  2.31e-07  2.5e-03  0.0e+00  1.4e-03  2.31e-07
##     19   36 -1.947e+03  2.25e-07  1.75e-07  2.2e-03  0.0e+00  1.4e-03  1.75e-07
##     20   37 -1.947e+03  6.16e-08  4.86e-08  1.0e-03  0.0e+00  6.4e-04  4.86e-08
##     21   38 -1.947e+03  1.16e-08  9.94e-09  4.5e-04  0.0e+00  2.5e-04  9.94e-09
##     22   39 -1.947e+03  8.74e-10  7.68e-10  8.2e-05  0.0e+00  6.2e-05  7.68e-10
##     23   40 -1.947e+03  3.90e-11  3.57e-11  2.0e-05  0.0e+00  1.3e-05  3.57e-11
## 
##  ***** RELATIVE FUNCTION CONVERGENCE *****
## 
##  FUNCTION    -1.946628e+03   RELDX        2.006e-05
##  FUNC. EVALS      40         GRAD. EVALS      24
##  PRELDF       3.567e-11      NPRELDF      3.567e-11
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    1.436507e-03     1.000e+00    -1.670e-01
##      2    7.953834e-02     1.000e+00     8.530e-04
##      3    2.231052e-01     1.000e+00    -7.066e-04
##      4    2.732265e-01     1.000e+00    -5.770e-04

summary(x.fit)

## 
## Call:
## garch(x = x, order = c(0, 3))
## 
## Model:
## GARCH(0,3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2420 -0.3985  0.1671  0.7501  4.4193 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 1.437e-03   6.903e-05   20.809  < 2e-16 ***
## a1 7.954e-02   2.821e-02    2.820   0.0048 ** 
## a2 2.231e-01   2.587e-02    8.624  < 2e-16 ***
## a3 2.732e-01   4.384e-02    6.232 4.61e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 585.27, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 8.8831e-06, df = 1, p-value = 0.9976

ARCH(3) 模型参数均显著，最后拟合的是 ARCH(3) 模型。模型的形式：

\[E(\epsilon_t^2|\epsilon_{t-1},\epsilon_{t-2},...)=0.001437+0.07954\epsilon_{t-1}^2+0.2231\epsilon_{t-2}^2+0.2732\epsilon_{t-3}^2\]

# 绘制条件异方差模型拟合的95%置信区间
x.pred<-predict(x.fit)
plot(x.pred)

考虑到原序列的波动特征，显然 ARCH 模型更好拟合原序列的集群效应波动特征。

#条件异方差置信区间和方差齐性置信区间比较图示
plot(x)
lines(x.pred[,1],col=2)
lines(x.pred[,2],col=2)
abline(h=1.96*sd(x),col=4,lty=2)
abline(h=-1.96*sd(x),col=4,lty=2)

由上图，ARCH 模型拟合的置信区间比无条件方差两条平行线给出的 95% 置信区间更加符合原序列真实波动情况，即 ARCH 模型对序列波动预测更加准确。

1.4.3.2 GARCH 模型

ARCH 模型适用于异方差函数短期自相关过程，而实践中残差序列异方差函数有些是具有长期自相关性，此时采用 GARCH(p,q) 广义自回归条件异方差模型:

\[h_t=\omega+\sum_{i=1}^p \eta_i h_{t-i}+\sum_{j=1}^q \lambda_j \epsilon_{t-j}^2\]

实质上是在 ARCH 模型基础上增加考虑了异方差函数的 p 阶自相关性而形成的。

# 这是1979年12月31日-1991年12月31日外币兑美元的日兑换率序列
w<-read.table("file23.csv",sep=",",header = T)
x<-ts(w$exchange_rates,start=c(1979,12,31),frequency = 365)
plot(x)

序列时序图显示该序列非平稳，有明显的趋势特征。

# 差分序列时序图
plot(diff(x))

差分后原序列的趋势消除，但有明显的集群效应，故需要分析该序列时同时提取水平相关信息和波动相关信息。

# 考察差分序列的自相关、偏自相关
par(mfrow=c(2,1))
acf(diff(x))
pacf(diff(x))

由上，对于水平信息的提取，拟合 ARIMA(0,1,1) 模型。

#水平相关信息提取，拟合ARIMA(0,1,1)模型
x.fit<-arima(x,order = c(0,1,1))
x.fit

## 
## Call:
## arima(x = x, order = c(0, 1, 1))
## 
## Coefficients:
##          ma1
##       0.0357
## s.e.  0.0143
## 
## sigma^2 estimated as 0.0002007:  log likelihood = 13545.61,  aic = -27087.22

#残差白噪声检验
for (i in 1:6) print(Box.test(x.fit$residual,type = "Ljung-Box",lag=i))

## 
##  Box-Ljung test
## 
## data:  x.fit$residual
## X-squared = 0.0005354, df = 1, p-value = 0.9815
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual
## X-squared = 0.55102, df = 2, p-value = 0.7592
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual
## X-squared = 2.6528, df = 3, p-value = 0.4483
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual
## X-squared = 3.3062, df = 4, p-value = 0.5079
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual
## X-squared = 6.8276, df = 5, p-value = 0.2338
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual
## X-squared = 6.8306, df = 6, p-value = 0.3368

白噪声检验显示该 ARIMA(0,1,1) 模型显著成立，利用该拟合模型还可以预测序列未来的水平。

#水平预测，并绘制预测图
library(forecast)
x.fore<-forecast(x.fit,h=365)
plot(x.fore)

波动信息的提取需要首先考察 ARIMA(0,1,1) 模型残差平方序列的异方差特征。

#条件异方差检验(Portmanteau Q检验)
for (i in 1:6) print(Box.test(x.fit$residual^2,type = "Ljung-Box",lag=i))

## 
##  Box-Ljung test
## 
## data:  x.fit$residual^2
## X-squared = 82.803, df = 1, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual^2
## X-squared = 237.9, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual^2
## X-squared = 343.33, df = 3, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual^2
## X-squared = 490.84, df = 4, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual^2
## X-squared = 602.1, df = 5, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  x.fit$residual^2
## X-squared = 841.96, df = 6, p-value < 2.2e-16

检验显示残差序列显著方差非齐性，且具有长期相关性，构造 GARCH(1,1) 模型。

#拟合GARCH(1,1)模型
r.fit<-garch(x.fit$residual,order=c(1,1))

## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     1.806109e-04     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1 -1.800e+04
##      1    7 -1.800e+04  1.28e-04  3.05e-04  1.0e-04  5.5e+10  1.0e-05  8.36e+06
##      2    8 -1.800e+04  5.34e-06  5.87e-06  9.8e-05  2.0e+00  1.0e-05  2.84e+01
##      3   16 -1.807e+04  3.63e-03  6.18e-03  5.2e-01  2.0e+00  1.1e-01  2.83e+01
##      4   18 -1.819e+04  6.73e-03  9.78e-03  7.7e-01  2.0e+00  4.4e-01  5.00e+00
##      5   22 -1.819e+04  1.22e-04  1.77e-02  1.9e-02  2.0e+00  2.0e-02  5.70e-01
##      6   26 -1.831e+04  6.22e-03  2.96e-03  1.3e-01  1.7e+00  1.6e-01  1.42e-02
##      7   32 -1.833e+04  1.49e-03  1.58e-03  2.4e-06  4.5e+01  3.2e-06  2.29e-01
##      8   33 -1.837e+04  2.05e-03  2.29e-03  4.8e-06  4.0e+00  6.3e-06  2.84e+00
##      9   34 -1.837e+04  1.06e-04  4.42e-04  3.6e-06  2.0e+00  6.3e-06  4.46e-01
##     10   35 -1.838e+04  1.63e-04  2.09e-04  4.3e-06  2.0e+00  6.3e-06  9.74e-02
##     11   36 -1.838e+04  5.21e-06  4.66e-06  4.6e-06  2.0e+00  6.3e-06  1.36e-01
##     12   37 -1.838e+04  1.48e-07  1.63e-07  4.6e-06  2.0e+00  6.3e-06  1.42e-01
##     13   44 -1.838e+04  1.48e-04  2.91e-04  1.9e-02  2.0e+00  2.6e-02  1.42e-01
##     14   46 -1.842e+04  2.33e-03  1.15e-03  4.6e-02  0.0e+00  8.6e-02  1.15e-03
##     15   48 -1.844e+04  9.88e-04  9.84e-04  1.8e-02  1.8e+00  3.4e-02  1.46e-02
##     16   50 -1.848e+04  1.89e-03  1.90e-03  3.4e-02  1.5e+00  6.9e-02  3.78e-02
##     17   51 -1.849e+04  8.65e-04  3.36e-03  5.8e-02  1.6e+00  1.4e-01  2.18e-02
##     18   53 -1.850e+04  2.20e-04  5.04e-03  4.1e-03  1.5e+00  9.4e-03  5.09e-03
##     19   56 -1.851e+04  9.04e-04  7.21e-04  5.0e-03  1.5e+00  1.1e-02  1.25e-03
##     20   57 -1.852e+04  5.17e-04  8.81e-04  5.9e-03  1.5e+00  1.1e-02  1.70e-03
##     21   58 -1.852e+04  1.15e-04  1.24e-04  4.7e-03  0.0e+00  1.0e-02  1.24e-04
##     22   59 -1.853e+04  3.66e-05  3.82e-05  4.4e-03  3.9e-01  1.0e-02  4.13e-05
##     23   61 -1.853e+04  1.72e-06  2.65e-06  2.1e-04  1.5e+00  4.8e-04  1.16e-05
##     24   71 -1.853e+04  2.00e-07  4.83e-07  1.2e-08  3.5e+00  2.3e-08  1.82e-06
##     25   81 -1.853e+04  4.30e-07  5.19e-07  3.9e-04  6.4e-01  8.9e-04  6.97e-07
##     26   82 -1.853e+04  5.83e-08  9.89e-08  2.6e-04  0.0e+00  5.8e-04  9.89e-08
##     27   83 -1.853e+04  8.04e-10  3.34e-11  3.5e-06  0.0e+00  7.1e-06  3.34e-11
##     28   84 -1.853e+04 -9.81e-12  1.56e-14  4.2e-08  0.0e+00  7.7e-08  1.56e-14
## 
##  ***** RELATIVE FUNCTION CONVERGENCE *****
## 
##  FUNCTION    -1.852505e+04   RELDX        4.199e-08
##  FUNC. EVALS      84         GRAD. EVALS      28
##  PRELDF       1.558e-14      NPRELDF      1.558e-14
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    2.133129e-06     1.000e+00    -1.886e+00
##      2    7.623127e-02     1.000e+00    -8.062e-03
##      3    9.143603e-01     1.000e+00    -6.761e-03

summary(r.fit)

## 
## Call:
## garch(x = x.fit$residual, order = c(1, 1))
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.83074 -0.58407  0.02616  0.58758  4.54060 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 2.133e-06   3.014e-07    7.077 1.48e-12 ***
## a1 7.623e-02   5.456e-03   13.972  < 2e-16 ***
## b1 9.144e-01   6.015e-03  152.009  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 319.23, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.28019, df = 1, p-value = 0.5966

# 绘制波动 95% 置信区间
r.pred<-predict(r.fit)
plot(r.pred)

综合水平模型与波动模型，最终拟合的完整拟合模型为：

\[\begin{equation} \begin{cases} & \nabla x_t=\epsilon_t+0.0357\epsilon_{t-1}+v_t,v_t \sim WN(0,0.0002) \\ & v_t=\sqrt{h_t}e_t,e_t \overset{i.i.d}{\sim} N(0,\sigma^2) \\ & h_t=0.9144h_{t-1}+0.07623v_{t-1}^2 \end{cases} \end{equation}\]

1.5 Homework

1.5.1 1.

a1=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/5th_homework/5.1.txt')
a1=ts(a1)
ts.plot(a1)

a1.diff_1=diff(a1)
plot(a1.diff_1)

由于一阶差分提取出线性趋势，差分后序列表现出非常平稳的随机波动，故我们选取ARIMA(0,1,0)模型进行拟合，即随机游走模型。

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

估计下一天收盘价为：

a1.fit=arima(a1,order=c(0,1,0))
a1.fit

## 
## Call:
## arima(x = a1, order = c(0, 1, 0))
## 
## 
## sigma^2 estimated as 13.08:  log likelihood = -286.69,  aic = 575.39

library(zoo)
library(forecast)
a1.fore=forecast(a1.fit,h=1)
a1.fore

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 108            289 284.3642 293.6358 281.9102 296.0898

故下一次收盘价预测为289。

1.5.2 2.

a2=read.table('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/5th_homework/5.2.txt',header=T)
a2_1=c(a2[,1],a2[,3],a2[,5])
a2_2=c(a2[,2],a2[,4],a2[,6])
a2=data.frame(年=a2_1,铁路货运量=a2_2)
a2=ts(a2[,2],start = 1949)
a2.fit=auto.arima(a2)
a2.fit

## Series: a2 
## ARIMA(0,2,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.4427  -0.3315
## s.e.   0.1190   0.1212
## 
## sigma^2 estimated as 55215220:  log likelihood=-598.66
## AIC=1203.32   AICc=1203.77   BIC=1209.5

根据auto.arima函数，采用ARIMA(0,2,2)进行拟合。

a2.fore=forecast(a2.fit,h=5)
a2.fore

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2009       344539.9 335017.1 354062.7 329976.0 359103.8
## 2010       359930.8 342306.4 377555.1 332976.7 386884.9
## 2011       375321.7 350848.2 399795.2 337892.7 412750.7
## 2012       390712.6 359649.2 421776.1 343205.1 438220.1
## 2013       406103.5 368449.6 443757.5 348516.7 463690.3

1.5.3 3.

a3=read.table('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/5th_homework/5.3.txt',header=T)
a3_1=c(a3[,1],a3[,3],a3[,5])
a3_2=c(a3[,2],a3[,4],a3[,6])
a3=data.frame(月度=a3_1,死亡人数=a3_2)
a3=ts(a3[,2],start = 1973,frequency = 12)
a3.fit=auto.arima(a3)
a3.fit

## Series: a3 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4264  -0.5584
## s.e.   0.1226   0.1787
## 
## sigma^2 estimated as 102999:  log likelihood=-425.53
## AIC=857.06   AICc=857.5   BIC=863.3

根据auto.arima函数，采用\(ARIMA(0,1,1)*ARIMA(0,1,1)_{12}\)模型进行拟合。

1.5.4 4.

a4=read.table('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/5th_homework/5.4.txt',header=T)
a4_1=c(a4[,1],a4[,3],a4[,5],a4[,7])
a4_2=c(a4[,2],a4[,4],a4[,6],a4[,8])
a4=data.frame(年=a4_1,出生率=a4_2)
a4=ts(a4[,2],start = 1750)
a4.fit=auto.arima(a4)
a4.fit

## Series: a4 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.4616  6.7627
## s.e.  0.0885  0.9550
## 
## sigma^2 estimated as 27.43:  log likelihood=-306.58
## AIC=619.15   AICc=619.4   BIC=626.97

根据auto.arima函数，采用AR(1)模型进行拟合。

plot(1750:1849,a4.fit$residuals)

根据残差图，我们认为本序列存在(综合型)异方差，由于不知道异方差函数具体形式，故拟合条件异方差模型。

或者：

for (i in 1:6){
  print(Box.test(a4.fit$residuals^2,type = 'Ljung-Box',lag=i))
}

## 
##  Box-Ljung test
## 
## data:  a4.fit$residuals^2
## X-squared = 32.488, df = 1, p-value = 1.199e-08
## 
## 
##  Box-Ljung test
## 
## data:  a4.fit$residuals^2
## X-squared = 32.773, df = 2, p-value = 7.647e-08
## 
## 
##  Box-Ljung test
## 
## data:  a4.fit$residuals^2
## X-squared = 33.257, df = 3, p-value = 2.843e-07
## 
## 
##  Box-Ljung test
## 
## data:  a4.fit$residuals^2
## X-squared = 33.877, df = 4, p-value = 7.896e-07
## 
## 
##  Box-Ljung test
## 
## data:  a4.fit$residuals^2
## X-squared = 34.467, df = 5, p-value = 1.922e-06
## 
## 
##  Box-Ljung test
## 
## data:  a4.fit$residuals^2
## X-squared = 34.642, df = 6, p-value = 5.055e-06

library(tseries)
a4.fit_1=garch(a4.fit$residuals,order=c(1,1))

## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     2.443331e+01     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1  2.049e+02
##      1    3  2.005e+02  2.16e-02  4.07e-02  2.0e-03  8.4e+02  1.0e-01  1.70e+01
##      2    5  2.002e+02  1.42e-03  1.39e-03  1.5e-04  1.4e+01  1.0e-02  1.53e+00
##      3    7  1.997e+02  2.73e-03  2.73e-03  3.1e-04  5.3e+00  2.0e-02  9.13e-01
##      4    9  1.996e+02  5.32e-04  5.32e-04  6.5e-05  1.0e+02  4.0e-03  7.96e-01
##      5   11  1.994e+02  1.05e-03  1.05e-03  1.3e-04  1.5e+01  8.0e-03  6.25e-01
##      6   13  1.993e+02  2.08e-04  2.08e-04  2.7e-05  3.0e+02  1.6e-03  5.87e-01
##      7   15  1.993e+02  4.15e-05  4.15e-05  5.4e-06  1.5e+03  3.2e-04  5.32e-01
##      8   17  1.993e+02  8.31e-06  8.31e-06  1.1e-06  7.6e+03  6.4e-05  5.22e-01
##      9   19  1.993e+02  1.66e-05  1.66e-05  2.1e-06  9.6e+02  1.3e-04  5.20e-01
##     10   21  1.993e+02  3.32e-06  3.32e-06  4.3e-07  1.9e+04  2.6e-05  5.20e-01
##     11   23  1.993e+02  6.64e-06  6.64e-06  8.6e-07  2.4e+03  5.1e-05  5.19e-01
##     12   25  1.993e+02  1.33e-06  1.33e-06  1.7e-07  4.8e+04  1.0e-05  5.19e-01
##     13   27  1.993e+02  2.66e-06  2.66e-06  3.4e-07  6.0e+03  2.0e-05  5.18e-01
##     14   29  1.993e+02  5.31e-07  5.31e-07  6.9e-08  1.2e+05  4.1e-06  5.18e-01
##     15   31  1.993e+02  1.06e-07  1.06e-07  1.4e-08  6.0e+05  8.2e-07  5.18e-01
##     16   33  1.993e+02  2.13e-07  2.13e-07  2.7e-08  7.5e+04  1.6e-06  5.18e-01
##     17   35  1.993e+02  4.25e-07  4.25e-07  5.5e-08  3.7e+04  3.3e-06  5.18e-01
##     18   37  1.993e+02  8.50e-08  8.50e-08  1.1e-08  7.5e+05  6.6e-07  5.18e-01
##     19   39  1.993e+02  1.70e-08  1.70e-08  2.2e-09  3.7e+06  1.3e-07  5.18e-01
##     20   41  1.993e+02  3.40e-08  3.40e-08  4.4e-09  4.7e+05  2.6e-07  5.18e-01
##     21   43  1.993e+02  6.80e-09  6.80e-09  8.8e-10  9.3e+06  5.2e-08  5.18e-01
##     22   45  1.993e+02  1.36e-08  1.36e-08  1.8e-09  1.2e+06  1.0e-07  5.18e-01
##     23   47  1.993e+02  2.72e-09  2.72e-09  3.5e-10  2.3e+07  2.1e-08  5.18e-01
##     24   49  1.993e+02  5.44e-09  5.44e-09  7.0e-10  2.9e+06  4.2e-08  5.18e-01
##     25   51  1.993e+02  1.09e-08  1.09e-08  1.4e-09  1.5e+06  8.4e-08  5.18e-01
##     26   54  1.993e+02  2.18e-10  2.18e-10  2.8e-11  2.9e+08  1.7e-09  5.18e-01
##     27   57  1.993e+02  1.74e-09  1.74e-09  2.3e-10  9.1e+06  1.3e-08  5.18e-01
##     28   60  1.993e+02  3.48e-11  3.48e-11  4.5e-12  1.7e+00  2.7e-10 -9.19e-02
##     29   62  1.993e+02  6.96e-12  6.96e-12  9.0e-13  1.7e+00  5.4e-11 -9.19e-02
##     30   64  1.993e+02  1.39e-11  1.39e-11  1.8e-12  1.7e+00  1.1e-10 -9.19e-02
##     31   66  1.993e+02  2.79e-12  2.79e-12  3.6e-13  1.7e+00  2.1e-11 -9.19e-02
##     32   68  1.993e+02  5.57e-13  5.57e-13  7.2e-14  1.7e+00  4.3e-12 -9.19e-02
##     33   71  1.993e+02  4.46e-12  4.46e-12  5.8e-13  1.7e+00  3.4e-11 -9.19e-02
##     34   74  1.993e+02  8.86e-14  8.91e-14  1.2e-14  1.7e+00  6.9e-13 -9.19e-02
##     35   76  1.993e+02  1.79e-13  1.78e-13  2.3e-14  1.7e+00  1.4e-12 -9.19e-02
##     36   78  1.993e+02 -5.02e+07  3.57e-14  4.6e-15  1.7e+00  2.7e-13 -9.19e-02
## 
##  ***** FALSE CONVERGENCE *****
## 
##  FUNCTION     1.993107e+02   RELDX        4.609e-15
##  FUNC. EVALS      78         GRAD. EVALS      36
##  PRELDF       3.565e-14      NPRELDF     -9.193e-02
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    2.443162e+01     1.000e+00     6.993e-01
##      2    1.768209e-01     1.000e+00    -1.480e+01
##      3    7.550488e-15     1.000e+00     2.118e+01

summary(a4.fit_1)

## 
## Call:
## garch(x = a4.fit$residuals, order = c(1, 1))
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3763 -0.4697  0.1184  0.4924  1.7087 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)   
## a0 2.443e+01   1.622e+01    1.506  0.13201   
## a1 1.768e-01   5.803e-02    3.047  0.00231 **
## b1 7.550e-15   5.959e-01    0.000  1.00000   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 59.059, df = 2, p-value = 1.498e-13
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 23.6, df = 1, p-value = 1.186e-06

综合来看，我们得到的模型为：

\[\begin{cases} & x_t=0.4616*x_{t-1}+\epsilon_t+v_t \quad v_t\sim N(6.7627,27.43)\\ & v_t=\sqrt{h_t}*e_t\\ & h_t=7.55*10^{-15}h_{t-1}+0.1768*v_{t-1}^2 \end{cases}\]

1.5.5 5.

a5=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/5th_homework/5.5.txt')
a5=ts(a5,start=1867)
ts.plot(a5)

adf.test(a5)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  a5
## Dickey-Fuller = -1.8653, Lag order = 4, p-value = 0.6302
## alternative hypothesis: stationary

由于p值为0.6302>\(\alpha\)，故不拒绝原假设，没有理由认为序列平稳。

a5.fit=auto.arima(a5)
a5.fit

## Series: a5 
## ARIMA(3,1,0) 
## 
## Coefficients:
##          ar1      ar2      ar3
##       0.4641  -0.2375  -0.2796
## s.e.  0.1268   0.1409   0.1267
## 
## sigma^2 estimated as 4991:  log likelihood=-401.89
## AIC=811.79   AICc=812.39   BIC=820.84

根据auto.arima函数，采用ARIMA(3,1,0)模型进行拟合。

a5.fore=forecast(a5.fit,h=7)
a5.fore

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1939       1871.380 1780.845 1961.914 1732.919 2009.841
## 1940       1880.349 1719.833 2040.866 1634.860 2125.839
## 1941       1819.559 1612.662 2026.456 1503.138 2135.980
## 1942       1766.741 1539.363 1994.119 1418.996 2114.486
## 1943       1754.162 1517.336 1990.989 1391.967 2116.358
## 1944       1777.871 1532.971 2022.771 1403.329 2152.413
## 1945       1806.630 1549.799 2063.461 1413.841 2199.420

1.5.6 6.

a6=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/5th_homework/5.6.txt')
a6=ts(a6,start = 1969,frequency = 12)
ts.plot(a6)

a6.diff=diff(a6)
plot(a6.diff)

plot(a6.diff^2)

ln_a6=log(a6)
ln_a6.diff=diff(ln_a6)
plot(ln_a6.diff)

原序列方差非齐，差分序列方差非齐，对数变换后，差分序列方差齐性。

a6.fit=auto.arima(a6)
a6.fit

## Series: a6 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1     ma1
##       -0.8112  0.9247
## s.e.   0.0626  0.0390
## 
## sigma^2 estimated as 0.266:  log likelihood=-231.45
## AIC=468.89   AICc=468.97   BIC=480.07

采用ARIMA(1,1,1)模型拟合，由于方差非齐，故使用条件异方差模型拟合。

a6.fit_1=garch(a6.fit$residuals,order=c(1,1))

## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     2.377185e-01     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1 -4.974e+01
##      1    4 -4.993e+01  3.68e-03  4.70e-03  1.3e-02  2.3e+03  1.0e-02  5.52e+00
##      2    7 -5.054e+01  1.21e-02  1.42e-02  7.1e-02  1.1e+01  4.0e-02  8.26e+00
##      3    9 -5.062e+01  1.61e-03  1.68e-03  6.1e-03  8.6e+00  4.0e-03  7.26e+00
##      4   11 -5.076e+01  2.89e-03  2.80e-03  1.4e-02  3.0e+00  8.0e-03  1.53e+01
##      5   14 -5.077e+01  6.32e-05  6.28e-05  3.3e-04  1.2e+02  1.6e-04  8.81e+01
##      6   16 -5.077e+01  1.27e-04  1.26e-04  6.6e-04  1.5e+02  3.2e-04  1.36e+03
##      7   18 -5.078e+01  2.55e-05  2.55e-05  1.3e-04  8.1e+04  6.4e-05  9.98e+04
##      8   20 -5.078e+01  5.11e-05  5.10e-05  2.6e-04  2.9e+06  1.3e-04  9.85e+06
##      9   22 -5.078e+01  1.03e-04  1.02e-04  5.3e-04  9.2e+02  2.6e-04  9.88e+08
##     10   27 -5.078e+01  2.06e-08  2.06e-08  1.1e-07  4.9e+05  5.1e-08  9.08e+10
##     11   29 -5.078e+01  4.11e-08  4.11e-08  2.1e-07  6.1e+05  1.0e-07  1.49e+12
##     12   31 -5.078e+01  8.23e-08  8.22e-08  4.2e-07  3.0e+06  2.0e-07  1.08e+14
##     13   33 -5.078e+01  1.65e-08  1.64e-08  8.4e-08  3.8e+08  4.1e-08  1.04e+16
##     14   35 -5.078e+01  3.29e-09  3.29e-09  1.7e-08  1.0e+10  8.2e-09  1.04e+18
##     15   37 -5.078e+01  6.58e-10  6.58e-10  3.4e-09  2.0e+11  1.6e-09  1.04e+20
##     16   39 -5.078e+01  1.32e-09  1.32e-09  6.7e-09  8.6e+10  3.3e-09  1.04e+22
##     17   41 -5.078e+01  2.63e-10  2.63e-10  1.3e-09  2.0e+11  6.6e-10  1.04e+24
##     18   43 -5.078e+01  5.26e-10  5.26e-10  2.7e-09  3.2e+10  1.3e-09  1.04e+26
##     19   45 -5.078e+01  1.05e-10  1.05e-10  5.4e-10  3.3e+10  2.6e-10  1.03e+28
##     20   47 -5.078e+01  2.11e-11  2.10e-11  1.1e-10  1.2e+01  5.2e-11 -3.83e-01
##     21   49 -5.078e+01  4.21e-11  4.21e-11  2.2e-10  2.1e+01  1.0e-10 -3.92e-01
##     22   51 -5.078e+01  8.42e-12  8.42e-12  4.3e-11  2.4e+00  2.1e-11 -2.54e-01
##     23   53 -5.078e+01  1.69e-12  1.68e-12  8.6e-12  2.8e+00  4.2e-12 -2.88e-01
##     24   55 -5.078e+01  3.37e-12  3.37e-12  1.7e-11  4.5e+00  8.4e-12 -3.44e-01
##     25   57 -5.078e+01  6.73e-13  6.73e-13  3.4e-12  2.0e+00  1.7e-12 -3.53e-01
##     26   59 -5.078e+01  1.35e-12  1.35e-12  6.9e-12  2.0e+00  3.4e-12 -3.27e-01
##     27   62 -5.078e+01  2.80e-14  2.69e-14  1.4e-13  2.0e+00  6.7e-14 -3.96e-01
##     28   64 -5.078e+01  5.43e-14  5.39e-14  2.8e-13  2.0e+00  1.3e-13 -3.56e-01
##     29   66 -5.078e+01  9.65e-15  1.08e-14  5.5e-14  2.0e+00  2.7e-14 -3.06e-01
##     30   68 -5.078e+01  3.08e-15  2.16e-15  1.1e-14  2.2e+00  5.4e-15 -2.30e-01
##     31   70 -5.078e+01  3.78e-15  4.31e-15  2.2e-14  3.7e+00  1.1e-14 -3.27e-01
##     32   72 -5.078e+01  9.09e-15  8.62e-15  4.4e-14  2.0e+00  2.1e-14 -3.39e-01
##     33   74 -5.078e+01  8.40e-16  1.72e-15  8.8e-15  2.0e+00  4.3e-15 -4.00e-01
##     34   75 -5.078e+01 -1.97e+08  3.45e-15  1.8e-14  2.0e+00  8.6e-15 -3.93e-01
## 
##  ***** FALSE CONVERGENCE *****
## 
##  FUNCTION    -5.078284e+01   RELDX        1.764e-14
##  FUNC. EVALS      75         GRAD. EVALS      34
##  PRELDF       3.448e-15      NPRELDF     -3.927e-01
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    2.426948e-01     1.000e+00    -1.507e+00
##      2    1.500637e-15     1.000e+00     2.032e+01
##      3    7.942475e-02     1.000e+00    -4.411e-01

summary(a6.fit_1)

## 
## Call:
## garch(x = a6.fit$residuals, order = c(1, 1))
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.73025 -0.37266 -0.01044  0.30144  4.59196 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)
## a0 2.427e-01   6.548e+00    0.037    0.970
## a1 1.501e-15   3.922e-02    0.000    1.000
## b1 7.942e-02   2.484e+01    0.003    0.997
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 302.6, df = 2, p-value < 2.2e-16
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.4351, df = 1, p-value = 0.5095

使用GARCH(1,1)模型，最终的模型如下：

\[\begin{cases} & \Delta x_t=-0.8112*\Delta x_{t-1}+\epsilon_t+0.9247*\epsilon_{t-1}+v_t \quad v_t\sim N(6.7627,27.43)\\ & v_t=\sqrt{h_t}*e_t\\ & h_t=0.07942*h_{t-1}+1.501*10^{-15}*v_{t-1}^2 \end{cases}\]