SXtask
ZHL_2023311075
2025-11-13
R Markdown
1.加载所需包
library(readxl)## Warning: package 'readxl' was built under R version 4.5.1library(aTSA)##
## Attaching package: 'aTSA'## The following object is masked from 'package:graphics':
##
## identifylibrary(forecast)## Warning: package 'forecast' was built under R version 4.5.1## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'forecast'## The following object is masked from 'package:aTSA':
##
## forecastex4.1–> 1900-1998年全球7级以上地震发生次数序列
#加载数据集A1_7
A1_7 <- read_excel("C:/Users/LunaLin/Desktop/时间序列数据/时间序列分析——基于R（第2版）案例数据/时间序列分析——基于R（第2版）案例数据/A1_7.xlsx")
number<-ts(A1_7$number,start=c(1900,1),frequency=12)
plot(number)
从时序图显示，该序列无明显的趋势或周期特征，即无显著的非平稳特征。
#绘制该序列自相关图和偏自相关图
par(mfrow=c(1,2))
acf(number)
pacf(number)
从本例的自相关图可看出，自相关系数是以一种有规律的方式，按指数函数轨迹衰减的，这说明自相关系数衰减到零下不是一个突然截尾的过程，这是自相关系数拖尾的典型特征；偏自相关系数图中，除了1阶偏自相关系数都在2倍标准差范围内，这表明偏自相关系数1阶结尾。综上，可用AR(1)模型来拟合该序列
#拟合AR(1)模型——例4-1续（1）
fit1<-arima(number,order=c(1,0,0),method="ML")
fit1##
## Call:
## arima(x = number, order = c(1, 0, 0), method = "ML")
##
## Coefficients:
## ar1 intercept
## 0.5432 19.8911
## s.e. 0.0840 1.3179
##
## sigma^2 estimated as 36.7: log likelihood = -318.98, aic = 643.97得到AR(1)模型口径：x_t–19.8911=e_t/(1-0.5432B) ，Var(e_t)=36.7
#拟合模型AR(1)显著性检验——例4-1续（2）
ts.diag(fit1)
考察残差序列的白噪声检验结果（Residual Diagnostics Plots 中左下图），可发现各阶延迟下白噪声检验统计量的P值都显著大于0.05(红色虚线为参考线)。我们可以认为这个拟合模型的残差序列属于白噪声序列，也就是拟合模型显著成立。
#参数显著性检验（近似）——例4-1续（3）
fit1##
## Call:
## arima(x = number, order = c(1, 0, 0), method = "ML")
##
## Coefficients:
## ar1 intercept
## 0.5432 19.8911
## s.e. 0.0840 1.3179
##
## sigma^2 estimated as 36.7: log likelihood = -318.98, aic = 643.97因为参数1的估计值大于它的两倍标准差（0.5432》2*0.084），所以参数1显著非零；同理，因为19.8911>2*1.3179，所以参数μ也显著非零。也就是拟合的AR(1)模型两参数都显著非零。
#构造t统计量，调用pt函数求P值
t=abs(fit1$coef)/sqrt(diag(fit1$var.coef))
pt(t,length(number)-length(fit1$coef),lower.tail=F)## ar1 intercept
## 2.002272e-09 1.687781e-27因为两参数t统计量的P值都远小于0.05，所以可以判断我们为例4.1序列拟合的AR(1)模型两参数都显著非零
#调用forecast函数，做10期预测——例4-1 续（4）
fore1<-forecast::forecast(fit1,h=10)
fore1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Apr 1908 17.77725 10.013811 25.54069 5.904094 29.65041
## May 1908 18.74274 9.907694 27.57778 5.230705 32.25477
## Jun 1908 19.26723 10.139955 28.39451 5.308265 33.22620
## Jul 1908 19.55217 10.340414 28.76392 5.464007 33.64032
## Aug 1908 19.70695 10.470420 28.94349 5.580895 33.83301
## Sep 1908 19.79104 10.547207 29.03487 5.653817 33.92826
## Oct 1908 19.83672 10.590734 29.08271 5.696204 33.97724
## Nov 1908 19.86154 10.614915 29.10816 5.720048 34.00303
## Dec 1908 19.87502 10.628208 29.12183 5.733243 34.01679
## Jan 1909 19.88234 10.635476 29.12921 5.740482 34.02420#输出预测图
plot(fore1,lty=2)
lines(fore1$fitted,col=4)
图中，虚线为观测值，实线为拟合值。深色阴影部分为置信水平为80%的预测值置信区间，浅色阴影部分为置信水平95%的预测值置信区间。
ex4.2 美国科罗拉多州某个加油站连续57天的盈亏序列 #读入数据文件，绘时序图
A1_8 <- read_excel("C:/Users/LunaLin/Desktop/时间序列数据/时间序列分析——基于R（第2版）案例数据/时间序列分析——基于R（第2版）案例数据/A1_8.xlsx")
print(A1_8)## # A tibble: 57 × 2
## day overshort
## <dbl> <dbl>
## 1 1 78
## 2 2 -58
## 3 3 53
## 4 4 -63
## 5 5 13
## 6 6 -6
## 7 7 -16
## 8 8 -14
## 9 9 3
## 10 10 -74
## # ℹ 47 more rowsovershort<-ts(A1_8$overshort)
plot(overshort)
时序图显示该序列没有明显的趋势或周期特征，说明该序列没有显著的非平稳特征，进一步进行ADF检验，判断该序列的平稳性#ADF检验及输出结果
adf.test(overshort)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -13.04 0.01
## [2,] 1 -7.32 0.01
## [3,] 2 -7.01 0.01
## [4,] 3 -6.04 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -13.10 0.01
## [2,] 1 -7.46 0.01
## [3,] 2 -7.38 0.01
## [4,] 3 -6.61 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -12.99 0.01
## [2,] 1 -7.39 0.01
## [3,] 2 -7.34 0.01
## [4,] 3 -6.60 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01 ADF检验结果显示，该序列的所有ADF检验统计量的P值都小于显著性水平，所以可以认为该序列为平稳序列。再对平稳序列进行纯随机性检验。#纯随机性检验及输出结果
for (k in 1:2)
print(Box.test(overshort,lag=6*k,type="Ljung-Box"))##
## 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随机性检验结果显示，6阶和12阶延迟的LB统计量的P值均小于0.05，所以可以判断该序列为平稳非白噪声序列，可使用ARIMA模型来拟合。
#绘制自相关图和偏自相关图
par(mfrow=c(1,2))
acf(overshort)
pacf(overshort)
自相关图显示除了延迟1阶的自相关系数再2倍标准差范围之外，其他阶数的自相关系数都在2倍标准差范围内波动，且自相关系数衰减没有显著的规律性；偏自相关系数显示有规律的衰减，即拖尾特征。综合该序列自相关系数1阶截尾和偏自相关系数拖尾特征，可将该序列的拟合模型定阶为MA(1)
#拟合MA(1)模型——例4-2续（1）
fit2<-arima(overshort,order=c(0,0,1),method="CSS")
fit2##
## Call:
## arima(x = overshort, order = c(0, 0, 1), method = "CSS")
##
## Coefficients:
## ma1 intercept
## -0.8208 -4.4095
## s.e. 0.0996 1.1655
##
## sigma^2 estimated as 2105: part log likelihood = -298.96可得到该模型口径为： x_t=-4.4095+e_t-0.8202e_(t-1)，Var(e_t)=2105
模型MA(1)显著性检验——例4-2续（2）
ts.diag(fit2)
模型的显著性检验结果显示，残差序列可以视为白噪声序列，所以拟合模型显著成立。
#参数显著性检验（近似方法）
fit2##
## Call:
## arima(x = overshort, order = c(0, 0, 1), method = "CSS")
##
## Coefficients:
## ma1 intercept
## -0.8208 -4.4095
## s.e. 0.0996 1.1655
##
## sigma^2 estimated as 2105: part log likelihood = -298.96#参数显著性检验（精确方法）
t=abs(fit2$coef)/sqrt(diag(fit2$var.coef))
pt(t,length(overshort)-length(fit2$coef),lower.tail = F)## ma1 intercept
## 1.775795e-11 1.919653e-04#输出t统计量p值参数的显著性检验结果显示，无论使用近似方法还是精确方法，都能得出两参数显著非零的结论，因其P值皆远小于0.05。
ex5.1 尝试提取1964-1999年中国沙年产量序列中的确定性信息
A1_11 <- read_excel("C:/Users/LunaLin/Desktop/时间序列数据/时间序列分析——基于R（第2版）案例数据/时间序列分析——基于R（第2版）案例数据/A1_11.xlsx")
sha<-ts(A1_11$sha,start=1964)
plot(sha,type="o")
从时序图发现，该序列呈明显线性递增趋势，需对该序列进行1阶差分提取线性趋势信息#一阶差分后序列
dif_sha<-diff(sha)
plot(dif_sha,type="o")
1阶差分运算非常成功地从原序列中提取出线性趋势，差分后序列呈现出非常平稳的波动特征.ex5.3 利用差分运算提取1962年1月至1975年12月平均每头奶牛地月产量序列中的确定性信息
#加载文件，绘制时序图
A1_13<-read_excel("C:/Users/LunaLin/Desktop/时间序列数据/时间序列分析——基于R（第2版）案例数据/时间序列分析——基于R（第2版）案例数据/A1_13.xlsx")
print(A1_13)## # A tibble: 168 × 2
## time milk
## <dttm> <dbl>
## 1 1962-01-01 00:00:00 589
## 2 1962-02-01 00:00:00 561
## 3 1962-03-01 00:00:00 640
## 4 1962-04-01 00:00:00 656
## 5 1962-05-01 00:00:00 727
## 6 1962-06-01 00:00:00 697
## 7 1962-07-01 00:00:00 640
## 8 1962-08-01 00:00:00 599
## 9 1962-09-01 00:00:00 568
## 10 1962-10-01 00:00:00 577
## # ℹ 158 more rowsmilk<-ts(A1_13$milk)
plot(milk)
时序图显示该序列具有一个线性递增的长期趋势和一个周期长度为一年的稳定的季节性波动。所以对原序列先做1阶差分，提取线性递增趋势。
#1阶差分后时序图
dif_milk<-diff(milk)
plot(dif_milk)
1阶差分后序列的时序图显示1阶差分能提取原序列中蕴含的线性递增趋势，但还残留季节变动和随机波动，故对1阶差分后序列再进行12步的周期差分，提取季节波动性息。
#1阶12步差分后时序图
dif_milk_12<-diff(dif_milk,12)
plot(dif_milk_12)
上图显示，周期差分可以非常好地提取周期信息。至此，1阶12步差分运算比较充分地提取了原序列中蕴含的长期趋势和季节效应等确定性信息。
ex5.6 对1898-1970年美国国民生产总值平减指数（GNP)序列建模 #读取数据，绘制时序图
A1_14<-read_excel("C:/Users/LunaLin/Desktop/时间序列数据/时间序列分析——基于R（第2版）案例数据/时间序列分析——基于R（第2版）案例数据/A1_14.xlsx")
GNP<-ts(A1_14$GNP,start=1889)
plot(GNP)
时序图显示该序列有显著线性递增趋势，这是典型的非平稳序列特征。#一阶差分图
dif_GNP<-diff(GNP)
plot(dif_GNP)
对序列进行1阶差分后可见，差分序列基本围绕在0值附近波动，无明显的趋势特征。为进一步确认差分后序列的平稳性，还需对差分后序列做ADF检验。 #差分后序列平稳性检验
adf.test(dif_GNP)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -4.38 0.0100
## [2,] 1 -3.25 0.0100
## [3,] 2 -2.61 0.0101
## [4,] 3 -2.26 0.0245
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -5.14 0.0100
## [2,] 1 -4.03 0.0100
## [3,] 2 -3.44 0.0143
## [4,] 3 -3.11 0.0326
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -5.73 0.0100
## [2,] 1 -4.62 0.0100
## [3,] 2 -4.03 0.0127
## [4,] 3 -3.78 0.0239
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01检验结果显示，该序列所以ADF检验统计量的P值均小于0.05，所以可确认1阶差分后序列平稳。再对其进行纯随机性检验 #差分后序列纯随机性检验
for (k in 1:2)
print(Box.test(dif_GNP,lag=6*k,type="Ljung-Box"))##
## Box-Ljung test
##
## data: dif_GNP
## X-squared = 25.334, df = 6, p-value = 0.000296
##
##
## Box-Ljung test
##
## data: dif_GNP
## X-squared = 28.09, df = 12, p-value = 0.005367检验结果显示，各阶延迟下LB统计量的P值小于0.05，这说明差分后序列不是白噪声序列。综上可知，1阶差分后序列为平稳非白噪声序列。
#差分后序列自相关图和偏自相关图
par(mfrow=c(1,2))
acf(dif_GNP)
pacf(dif_GNP)
1阶差分后序列的自相关图显示拖尾特征，偏自相关图显示1阶结尾特征，所以可考虑用AR(1)模型拟合1阶差分后序列，考虑到前面已经进行了1阶差分运算，所有最终我们为美国国民生产总值平减指数序列(GNP)拟合ARIMA(1,1,0)模型。
#提取训练集数据,选择1960年前的数据。基于该测试集，分别拟合无漂移项ARIMA(1,1,0)模型和有漂移项ARIMA(1,1,0)。
library(stats)
x<-window(GNP,stat=1889,end=1960)
fit_x<-arima(x,order=c(1,1,0))
fit_x##
## Call:
## arima(x = x, order = c(1, 1, 0))
##
## Coefficients:
## ar1
## 0.4761
## s.e. 0.1031
##
## sigma^2 estimated as 6.929: log likelihood = -169.59, aic = 343.18#无漂移项拟合模型ARIMA(1,1,0)显著性检验
ts.diag(fit_x)
#调用forecast包中的Arima函数，对序列拟合ARIMA（1，1，0）
fit_x2<-Arima(x,order=c(1,1,0),include.drift = T)
fit_x2## Series: x
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.3983 1.0812
## s.e. 0.1078 0.5006
##
## sigma^2 = 6.752: log likelihood = -167.61
## AIC=341.23 AICc=341.59 BIC=348.02#有漂移项拟合模型显著性检验
tsdiag(fit_x2)
模型的残差检验显示，两个模型都显著成立，它们最大的差异会体现在预测上，有无漂移项对预测趋势会有显著影响。对GNP做为期10年的预测——例5.6续 #基于无漂移项模型的预测
fore_x1<-forecast::forecast(fit_x,h=10)
fore_x1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1961 104.1094 100.73602 107.4828 98.95025 109.2686
## 1962 104.4948 98.48016 110.5094 95.29621 113.6934
## 1963 104.6783 96.36127 112.9953 91.95850 117.3981
## 1964 104.7657 94.44647 115.0848 88.98383 120.5475
## 1965 104.8073 92.72643 116.8881 86.33123 123.2833
## 1966 104.8271 91.17236 118.4818 83.94399 125.7101
## 1967 104.8365 89.75480 119.9182 81.77104 127.9019
## 1968 104.8410 88.44880 121.2332 79.77130 129.9107
## 1969 104.8431 87.23447 122.4518 77.91301 131.7732
## 1970 104.8441 86.09634 123.5919 76.17185 133.5164将序列最后10期的真实值作为测试集，与预测集进行比较，得到预测误差序列。 #基于无漂移项拟合模型的预测误差
test<-window(GNP,start=1961)
error1<-test-fore_x1$mean
file1<-data.frame(fore_x1$mean,test,error1)
file1## fore_x1.mean test error1
## 1 104.1094 104.6 0.4905844
## 2 104.4948 105.8 1.3052000
## 3 104.6783 107.2 2.5217082
## 4 104.7657 108.8 4.0343428
## 5 104.8073 110.9 6.0927458
## 6 104.8271 113.9 9.0729403
## 7 104.8365 117.6 12.7635104
## 8 104.8410 122.3 17.4590206
## 9 104.8431 128.2 23.3568828
## 10 104.8441 135.3 30.4558650结果显示，无漂移项模型的预测值几乎为常数，预测误差很大。
#基于无漂移项拟合模型的预测效果图
plot(fore_x1)
lines(fore_x1$fitted,col=2,lty=2)
拟合模型的预测效果图中，实线为序列观测值，虚线为模型拟合值，阴影部分实线为预测值，深色阴影为序列80%置信区间，浅色阴影为序列95%置信区间。
#基于有漂移项拟合模型的预测
fore_x2<-forecast::forecast(fit_x2,h=10)
fore_x2## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1961 104.6277 101.29768 107.9577 99.53489 109.7205
## 1962 105.8070 100.08240 111.5317 97.05196 114.5621
## 1963 106.9273 99.20370 114.6509 95.11507 118.7396
## 1964 108.0241 98.60261 117.4455 93.61519 122.4329
## 1965 109.1114 98.21258 120.0103 92.44307 125.7798
## 1966 110.1951 97.98141 122.4087 91.51588 128.8743
## 1967 111.2772 97.87167 124.6828 90.77519 131.7793
## 1968 112.3588 97.85696 126.8606 90.18016 134.5374
## 1969 113.4401 97.91850 128.9617 89.70186 137.1784
## 1970 114.5213 98.04266 131.0000 89.31938 139.7233#基于有漂移项拟合模型的预测误差
error2<-test-fore_x2$mean
file2<-data.frame(fore_x2$mean,test,error2)
file2## fore_x2.mean test error2
## 1 104.6277 104.6 -0.02767376
## 2 105.8070 105.8 -0.00703442
## 3 106.9273 107.2 0.27268423
## 4 108.0241 108.8 0.77593664
## 5 109.1114 110.9 1.78856353
## 6 110.1951 113.9 3.70492467
## 7 111.2772 117.6 6.32277331
## 8 112.3588 122.3 9.94121449
## 9 113.4401 128.2 14.75989170
## 10 114.5213 135.3 20.77866293有漂移项模型的预测值呈线性递增，短期内预测误差很小，但预测期数变大，预测误差也变大。但总体而言，此例中有漂移项模型的预测效果显著优于无漂移项模型的预测效果。
#基于有漂移项拟合模型的预测图效果
plot(fore_x2)
lines(fore_x2$fitted,col=2,lty=2)
预测效果图显示：这两个拟合模型对测试集数据的拟合效果都不错，但是它们的预测效果有显著差别。
无漂移项模型的预测均值几乎为常数，基本改变了原序列的发展趋势，预测误差很大。有漂移项模型的预测效果较好，预测值基本延续了原序列的发展趋势。