Untitled

1绘制该序列时序图：

# 加载必要的包
library(forecast)

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

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

library(tseries)

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

# 读取数据并处理缺失值
data <- read.csv("zy4.csv", header = FALSE)
values <- as.numeric(as.vector(t(as.matrix(data))))  # 转换为数值向量
ts_data <- ts(na.omit(values), start = 1, frequency = 1)  # 创建时间序列对象

# 绘制时序图
plot(ts_data, 
     main = "年度人口死亡率时序图", 
     xlab = "年份", 
     ylab = "死亡率",
     col = "blue", lwd = 2)

2判断该序列的平稳性与纯随机性： 结论： 若ADF p < 0.05，序列平稳；若KPSS p < 0.05，序列非平稳。 若Ljung-Box p < 0.05，序列非纯随机。

Warning in kpss.test(ts_data): p-value greater than printed p-value

ADF检验p值: 0.08757116 

KPSS检验p值: 0.1 

Ljung-Box检验p值: 0.01584384 

3考察该序列的自相关系数和偏自相关系数的性质： 分析： ACF/PACF截尾或拖尾可帮助选择ARIMA模型的参数。

# 绘制ACF和PACF图
par(mfrow = c(1, 2))
acf(ts_data, main = "自相关系数 (ACF)", lag.max = 20)
pacf(ts_data, main = "偏自相关系数 (PACF)", lag.max = 20)

par(mfrow = c(1, 1))

4. 选择AIC最小的模型作为最优模型

# 自动选择ARIMA模型
auto_model <- auto.arima(ts_data)
summary(auto_model)

Series: ts_data 
ARIMA(0,0,1) with non-zero mean 

Coefficients:
         ma1    mean
      0.4495  4.9444
s.e.  0.1203  0.2020

sigma^2 = 0.9269:  log likelihood = -61.23
AIC=128.47   AICc=129.05   BIC=133.89

Training set error measures:
                      ME      RMSE      MAE       MPE     MAPE      MASE
Training set 0.007253807 0.9411272 0.728869 -4.898082 17.12819 0.7452129
                   ACF1
Training set 0.01208157

# 手动尝试其他模型（示例）
model_ar1 <- Arima(ts_data, order = c(1, 0, 0))  # AR(1)
model_ma1 <- Arima(ts_data, order = c(0, 0, 1))  # MA(1)
model_arma <- Arima(ts_data, order = c(1, 0, 1)) # ARMA(1,1)

# 比较AIC值
cat("模型AIC对比:\n")

模型AIC对比:

cat("Auto ARIMA:", AIC(auto_model), "\n")

Auto ARIMA: 128.4692 

cat("AR(1):", AIC(model_ar1), "\n")

AR(1): 131.223 

cat("MA(1):", AIC(model_ma1), "\n")

MA(1): 128.4692 

cat("ARMA(1,1):", AIC(model_arma), "\n")

ARMA(1,1): 130.4451 

# 使用最优模型预测
forecast_result <- forecast(auto_model, h = 5)
plot(forecast_result, 
     main = "未来5年人口死亡率预测",
     xlab = "年份",
     ylab = "死亡率",
     col = "darkred", lwd = 2)
lines(fitted(auto_model), col = "green")  # 添加拟合值

# 输出预测数值
print(forecast_result)

   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
46       4.842503 3.608670 6.076336 2.955518 6.729488
47       4.944402 3.591665 6.297138 2.875569 7.013234
48       4.944402 3.591665 6.297138 2.875569 7.013234
49       4.944402 3.591665 6.297138 2.875569 7.013234
50       4.944402 3.591665 6.297138 2.875569 7.013234