arimax的r语言代码可直接运行

arimax是arima带外部变量的形式。一般的arima模型是单变量模型，而arimax是多变量模型。

#下载包
install.packages("Metrics")
install.packages("forecast")
install.packages("readxl")
#install.packages("openxlsx")
install.packages("broom")
install.packages("stargazer")

# 加载必要的库
library(stats)
library(tseries)
library(readxl)
library(forecast)
library(Metrics)
#library(openxlsx)
library(broom)
library(stargazer)
# 导入数据
excel_data <- read_excel("C:/Users/Administrator/Desktop/price/数据.xlsx")
excel_data


time_series <- excel_data$'lnfp'
exogenous <- excel_data$'om'

# 将数据分为训练集和测试集，训练集样本内回归
train_size<-3871
train_data <- time_series[1:train_size]
test_data <- time_series[(train_size+1):length(time_series)]
train_exog <- exogenous[1:train_size]
test_exog <- exogenous[(train_size+1):length(exogenous)]

# 初始化一个向量来存储预测结果
forecast_results <- numeric(length(test_data))

#自动确定p和dq参数并样本内回归
fit<-auto.arima(train_data,max.p = 10,max.q = 10,start.p = 1,
                start.q = 1,xreg=train_exog)
#fit<-auto.arima(train_data,max.p = 10,max.q = 10,start.p = 1,
#                start.q = 1,xreg=train_exog,ic = c("aicc"))
#fit<-auto.arima(train_data,max.p = 10,max.q = 10,start.p = 1,
#                start.q = 1,xreg=train_exog,ic = c("bic"))
summary(fit)
#help("auto.arima")
#提取并汇报回归系数
# 提取系数
# 从模型的方差协方差矩阵中提取标准误差
stderr <- sqrt(diag(vcov(fit)))

# 计算t值
t_values <- fit$coef / stderr

# 计算p值
p_values <- 2 * pt(-abs(t_values), df=length(residuals(fit)) - length(fit$coef))

# 创建一个数据框，将系数、t值和p值放入其中
results_df <- data.frame(
  Term = names(fit$coef),
  Estimate = fit$coef,
  `Std. Error` = stderr,
  `t value` = t_values,
  `Pr(>|t|)` = p_values
)

# 使用stargazer函数生成表格
stargazer(results_df, type="text", summary=FALSE)

#模型检验
qqnorm(fit$residuals)
qqline(fit$residuals)
Box.test(fit$residuals,type="Ljung-Box",lag=10)

#help(stargazer)

# 滚动预测
for (i in 1:length(test_data)) {
  # 使用ARIMAX模型进行拟合
  fit <- arima(train_data, order = c(4,1,0), xreg = train_exog)
  #order中是前面确定的pdq参数值
  # 进行一步预测
  next_pred <- predict(fit, n.ahead = 1, newxreg = test_exog[i])$pred
  
  # 将预测值存储起来
  forecast_results[i] <- next_pred
  
  # 更新训练集，加入新的观测值
  train_data <- c(train_data, test_data[i])
  train_exog <- c(train_exog, test_exog[i])
}

# 比较预测结果和实际值
comparison <- data.frame(
  Actual = test_data,
  Forecast = forecast_results
)

# 打印比较结果
print(comparison)

# 使用Metrics包中的函数计算MAE, MSE, RMSE
mae_value <- mae(test_data, forecast_results)
mse_value <- mse(test_data, forecast_results)
rmse_value <- rmse(test_data, forecast_results)

# 打印误差结果
print(paste("MSE:", mse_value))
print(paste("RMSE:", rmse_value))
print(paste("MAE:", mae_value))

#输出结果
my_data<- cbind(test_data,forecast_results)
write.csv(my_data, file="C:/Users/Administrator/Desktop/oilprice1/my_data.csv")