🍉CSDN小墨&晓末:https://blog.csdn.net/jd1813346972
个人介绍: 研一|统计学|干货分享
擅长Python、Matlab、R等主流编程软件
累计十余项国家级比赛奖项,参与研究经费10w、40w级横向
1 目的
有1971年9月至1993年6月澳大利亚季度常住人口变动情况构成的时间序列 x t x_t xt,判断该时序的平稳性和纯随机性,并拟合适当的ARMA模型预测其未来五年情况。
2 平稳性检验
2.1 原始时序图
运行程序:
data<-scan("F:\\时间序列分析\\数据.txt")#读取数据
data1<-ts(data,start = c(1971,3),frequency =4) #转化为时序数据
plot(data1) #绘制时序图
运行结果:
由图1可以看出,根据时序图显示,该时序数据没有明显的趋势和周期特征,进一步进行ADF检验。
2.2 ADF检验
运行程序:
library(aTSA) #加载包
adf.test(data1) #ADF检验
运行结果:
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -2.719 0.010
## [2,] 1 -1.531 0.128
## [3,] 2 -0.928 0.345
## [4,] 3 -0.698 0.428
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -10.12 0.010
## [2,] 1 -6.41 0.010
## [3,] 2 -3.56 0.010
## [4,] 3 -2.32 0.207
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -10.48 0.0100
## [2,] 1 -6.88 0.0100
## [3,] 2 -3.92 0.0172
## [4,] 3 -2.57 0.3362
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
根据ADF检验结果显示该序列平稳。接下来检验序列的纯随机性。
3 纯随机检验
运行程序:
for (k in 1:2){ #纯随机检验
print(Box.test(data1,lag=2*k,type="Ljung-Box"))
}
运行结果:
##
## Box-Ljung test
##
## data: data1
## X-squared = 1.211, df = 2, p-value = 0.5458
##
##
## Box-Ljung test
##
## data: data1
## X-squared = 16.325, df = 4, p-value = 0.002612
通过纯随机检验的LB检验结果显示该序列为非白噪声序列,所以1971年9月至1993年6月澳大利亚常住人口变动情况时序数据为平稳非白噪声序列。
4 ARMA模型拟合
4.1 自相关图
运行程序:
acf(data1) #自相关图
运行结果:
4.2 偏自相关图
运行程序:
pacf(data1) #偏自相关图
运行结果:
4.3 ARMA模型阶数判定
根据自相关图和偏自相关图显示结果,自相关图和偏自相关图均显示拖尾特点,故考虑建立ARMA模型,尝试建立ARMA(1,1)-ARMA(4,4),通过AIC准则选择最优模型。
运行程序:
#拟合ARMA(1,1)-ARMA(4,4)
k<-matrix(rep(0,16),nrow = 4,ncol = 4) #创建0矩阵
for(i in 1:4){
for (j in 1:4) {
k[i,j]<-arima(data1,order = c(i,0,j),method="ML")$aic
}
}
which(k==min(k),arr.ind = T) #提取aic最小值位置
运行结果:
## row col
## [1,] 3 2
运行程序:
min(k) #最小aic
运行结果:
## [1] 768.0453
根据结果显示,最小AIC模型为ARMA(3,2),其AIC值为768.0453。即选址最优模型ARMA(3,2)。
4.4 模拟拟合
运行程序:
model<-arima(data1,order = c(3,0,2),method="ML") #最优模型
4.5 模型预测
运行程序:
library("forecast")
fore1<-forecast::forecast(model,h=20) #模型预测
fore1 #预测结果及0.05和0.5显著性水平下的置信区间
运行结果:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1993 Q3 55.84804 33.51665 78.17942 21.69513 90.00094
## 1993 Q4 54.40362 31.20745 77.59979 18.92814 89.87910
## 1994 Q1 47.24138 23.85678 70.62597 11.47772 83.00503
## 1994 Q2 49.48305 23.78938 75.17673 10.18797 88.77814
## 1994 Q3 52.73426 26.70517 78.76335 12.92620 92.54232
## 1994 Q4 52.00451 25.97331 78.03571 12.19323 91.81580
## 1995 Q1 51.19748 25.07768 77.31729 11.25070 91.14427
## 1995 Q2 51.87825 25.63253 78.12396 11.73889 92.01761
## 1995 Q3 52.35380 26.08425 78.62335 12.17799 92.52960
## 1995 Q4 52.21059 25.93444 78.48674 12.02468 92.39649
## 1996 Q1 52.17543 25.88658 78.46428 11.97010 92.38076
## 1996 Q2 52.34066 26.04270 78.63862 12.12140 92.55992
## 1996 Q3 52.41657 26.11588 78.71727 12.19313 92.64001
## 1996 Q4 52.40455 26.10234 78.70675 12.17880 92.63029
## 1997 Q1 52.42316 26.11950 78.72682 12.19519 92.65113
## 1997 Q2 52.46207 26.15758 78.76655 12.23283 92.69131
## 1997 Q3 52.47818 26.17334 78.78302 12.24840 92.70796
## 1997 Q4 52.48233 26.17727 78.78740 12.25220 92.71246
## 1998 Q1 52.49184 26.18661 78.79707 12.26147 92.72222
## 1998 Q2 52.50167 26.19635 78.80699 12.27116 92.73218
运行程序:
plot(fore1,lty=2) #绘制预测图
lines(fore1$fitted,col=4) #添加拟合值
运行结果:
通过图4可以看到该序列的拟合图及未来五年的预测图,拟合效果较好。