基于WIN10的64位系统演示

一、写在前面

这一期，我们使用R进行SARIMA模型的构建。

同样，这里使用这个数据：

《PLoS One》2015年一篇题目为《Comparison of Two Hybrid Models for Forecasting the Incidence of Hemorrhagic Fever with Renal Syndrome in Jiangsu Province, China》文章的公开数据做演示。数据为江苏省2004年1月至2012年12月肾综合症出血热月发病率。运用2004年1月至2011年12月的数据预测2012年12个月的发病率数据。

R语言，这个大家都很熟悉吧。说实话，我不熟悉。

这里用的版本是：R-4.3.1，搭配RStudio食用。

二、R建立SARIMA实战

（1）导入数据

打开Rstudio，如下图操作，导入数据：

（2）单位根（ADF）检验

library(tseries)
# 假设你的时间序列数据存储在变量data$incidence中
test_result <- adf.test(data$incidence, alternative = "stationary")
# 打印测试结果
print(test_result)

结果显示是平稳的，6666。那就平稳吧：

（3）设置为时间序列格式

# 将时间列转换为日期格式
time_series <- ts(data$incidence, frequency=12, start=c(2004, 1))

（4）差分

# 进行一次差分
first_difference <- diff(time_series, differences = 1)
plot(first_difference, main="一次差分后的时间序列")
# 进行季节性差分
seasonal_difference <- diff(time_series, lag = 12)
plot(seasonal_difference, main="季节性差分后的时间序列")

如图：

看起来，一次差分也得了。

（5）自相关图和偏相关图

绘制自相关图 (ACF)
acf(time_series, main="自相关图 (ACF)")
# 绘制偏自相关图 (PACF)
pacf(time_series, main="偏自相关图 (PACF)")

如图，有点丑：

（6）建模

（6.1）数据拆分

# 划分训练集和验证集
train_series <- window(time_series, start=c(2004,1), end=c(2011,12))
validation_series <- window(time_series, start=c(2012,1), end=c(2012,12))

（6.2）搭建SARIMA

看代码，自行体会：

案例：SARIMA(0,1,1)(0,1,1)12：

sarima_model <- Arima(train_series, order=c(0,1,1), seasonal=list(order=c(0,1,1), period=12))

解读：这里，order=c(0,1,1)定义了非季节部分的阶数，而seasonal=list(order=c(0,1,1), period=12)定义了季节部分的阶数和季节周期（在这种情况下为12）。

（6.3）看模型参数

# 显示模型摘要
summary(sarima_model)

如图：

解读如下：

参数没有统计学差异，这模型不得。大家自己试了，我继续。

（6）预测

（6.1）拟合数据

# 获取拟合数据
fitted_values <- fitted(sarima_model)
# 保存拟合数据到CSV文件
write.csv(data.frame(time=as.character(time(fitted_values)), fitted_values=fitted_values), file="fitted_values.csv", row.names=FALSE)
# 打印消息
print("拟合数据已保存到fitted_values.csv文件中")

知道存在哪里不，看R的工作路径：

# 获取当前工作目录
current_working_directory <- getwd()
# 打印当前工作目录
print(current_working_directory)

（6.1）预测数据

# 预测未来12个月的数据
future_forecast <- forecast(sarima_model, h=12)
# 将预测结果与时间戳合并为数据框
forecast_data <- data.frame(time = time(future_forecast$mean),forecast_values = as.numeric(future_forecast$mean)
)
# 保存预测结果到CSV文件
write.csv(forecast_data, file="forecast_values.csv", row.names=FALSE)
# 打印消息
print("预测数据已保存到forecast_values.csv文件中")

收工！！

四、数据

链接：https://pan.baidu.com/s/1qOpPi9pfzKR8TVmpOZaZcg?pwd=tc2z

提取码：tc2z