基于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