1 目的
有1971年9月至1993年6月澳大利亚季度常住人口变动情况构成的时间序列x t x_txt，判断该时序的平稳性和纯随机性，并拟合适当的ARMA模型预测其未来五年情况。
2 平稳性检验
2.1 原始时序图
运行程序：
data<-scan("F:\\时间序列分析\\数据.txt")#读取数据 data1<-ts(data,start = c(1971,3),frequency =4) #转化为时序数据 plot(data1) #绘制时序图
运行结果：
图1 1971年9月-1933年6月澳大利亚季度常住人口时序图
由图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) #自相关图
运行结果：
图2 自相关图
4.2 偏自相关图
运行程序：
pacf(data1) #偏自相关图
运行结果：
图3 偏自相关图
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 预测图
通过图4可以看到该序列的拟合图及未来五年的预测图，拟合效果较好。
