Untitled

# 加载库
library(tseries)

Warning: package 'tseries' was built under R version 4.4.3

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

library(forecast)

Warning: package 'forecast' was built under R version 4.4.3

# 1. 数据加载
data <- c(63.2, 67.9, 55.8, 49.5, 50.2, 55.4, 49.9, 45.3, 48.1, 61.7, 55.2, 53.1,
          49.5, 59.9, 30.6, 30.4, 33.8, 42.1, 35.8, 28.4, 32.9, 44.1, 45.5, 36.6,
          39.5, 49.8, 48.8, 29, 37.3, 34.2, 47.6, 37.3, 39.2, 47.6, 43.9, 49,
          51.2, 60.8, 67, 48.9, 65.4, 65.4, 67.6, 62.5, 55.1, 49.6, 57.3, 47.3,
          45.5, 44.5, 48, 47.9, 49.1, 48.8, 59.4, 51.6, 51.4, 60.9, 60.9, 56.8,
          58.6, 62.1, 64, 60.3, 64.6, 71, 79.4, 59.9, 83.4, 75.4, 80.2, 55.9,
          58.5, 65.2, 69.5, 59.1, 21.5, 62.5, 170, -47.4, 62.2, 60, 33.1, 35.3,
          43.4, 42.7, 58.4, 34.4)
ts_data <- ts(data, start = c(1971, 3), frequency = 4)

# 2. 平稳性检验
adf_test <- adf.test(ts_data)
print(adf_test)

    Augmented Dickey-Fuller Test

data:  ts_data
Dickey-Fuller = -1.9213, Lag order = 4, p-value = 0.6084
alternative hypothesis: stationary

#p=0.6084大于显著水平，无法拒绝原假设，所以序列是非平稳的。

3. 纯随机性检验

Box.test(ts_data, lag = 20, type = "Ljung-Box")

    Box-Ljung test

data:  ts_data
X-squared = 25.997, df = 20, p-value = 0.1659

#p = 0.1659：大于0.05，无法拒绝原假设，说明序列没有显著的自相关。

# 4. 拟合ARIMA模型
fit <- auto.arima(ts_data)
summary(fit)

Series: ts_data 
ARIMA(2,1,3) 

Coefficients:
         ar1      ar2      ma1     ma2      ma3
      0.1726  -0.5524  -1.4150  1.2054  -0.4778
s.e.  0.2276   0.1414   0.2418  0.3131   0.1894

sigma^2 = 336.8:  log likelihood = -375.15
AIC=762.31   AICc=763.36   BIC=777.1

Training set error measures:
                     ME     RMSE      MAE       MPE     MAPE      MASE
Training set -0.9105263 17.71563 10.78488 -4.048162 21.70858 0.7901698
                    ACF1
Training set -0.01058567

# 5. 绘制拟合图
plot(ts_data, main = "澳大利亚季度人口变动（1971-1993）", ylab = "变动（千人）")
lines(fitted(fit), col = "red")
legend("topright", legend = c("实际值", "拟合值"), col = c("black", "red"), lty = 1)

# 6. 预测未来5年
forecast_values <- forecast(fit, h = 20)
plot(forecast_values, main = "未来5年预测", xlab = "年份", ylab = "变动（千人）")