初探时间序列建模ARMA

knitr::opts_chunk$set(echo = TRUE)

数据清洗

包含趋势图、异常值处理、差分的方法

library(readxl)
OIL <- read_excel("C:\\Users\\213yi\\Desktop\\商务时间序列\\大作业\\清洗数据\\RCLC1d.xls", sheet = 2,skip = 2)
View(OIL)

##########时间趋势图##########
plot(OIL$Date, OIL$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)`,
     xlab = "时间", ylab = "WTI原油期货价格", type = "l", xaxt = "n")
tick_indices <- seq(1, length(OIL$Date), by = 1000)
x_ticks <- OIL$Date[tick_indices]
axis(1, at = x_ticks, labels = format(x_ticks, "%Y/%m/%d"))

##########一阶差分##########
plot(OIL$Date[-1], diff(OIL$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)`),
     type='l',ylab = '一阶差分',xlab='时间',  xaxt = "n")

tick_indices_diff <- seq(1, length(OIL$Date[-1]), by = 1000)
x_ticks_diff <- OIL$Date[-1][tick_indices]
axis(1, at = x_ticks_diff, labels = format(x_ticks_diff, "%Y/%m/%d"))

## 去掉异常值
OIL[OIL$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)`< 0,]

# 去掉2020-4-20的异常值
dat <- OIL[OIL$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)` >= 0, ]

# 时间趋势图(去负值)
par(mfrow = c(3, 1), mar = c(5, 4, 2, 1))
plot(dat$Date, dat$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)`,
     xlab = "时间", ylab = "价格", type = "l", xaxt = "n")
tick_indices <- seq(1, length(dat$Date), by = 1000)
x_ticks <- dat$Date[tick_indices]
axis(1, at = x_ticks, labels = format(x_ticks, "%Y/%m/%d"))

# 一阶差分（去负值）
plot(dat$Date[-1], diff(dat$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)`),
     type='l',ylab = '一阶差分',xlab='时间',  xaxt = "n")
tick_indices_diff <- seq(1, length(dat$Date[-1]), by = 1000)
x_ticks_diff <- dat$Date[-1][tick_indices]
axis(1, at = x_ticks_diff, labels = format(x_ticks_diff, "%Y/%m/%d"))

# 一阶对数差分（去负值）
plot(dat$Date[-1], diff(log(dat$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)`)),
     type='l',ylab = '一阶对数差分',xlab='时间',  xaxt = "n")
axis(1, at = x_ticks_diff, labels = format(x_ticks_diff, "%Y/%m/%d"))


price_WTI <- dat$`Cushing, OK Crude Oil Future Contract 1 (Dollars per Barrel)`
##plot(price_WTI,type='l',ylab = 'WTI原油期货价格',xlab='时间',xaxt = "n")
##axis(1,seq(1,length(price_WTI),length.out = 9),las = 1,label = seq(1983, 2024, by = 5),cex.axis = 1.3)


return <- diff(log(price_WTI))
#write.csv(return,file='w.csv')

定阶方法

ARMA定阶法
BIC定阶法（信息准则法）
自动定阶法auto.arima，便捷但是不一定最优

####ARMA定阶####
AICs <- matrix(0, max_lag, max_lag)
BICs <- matrix(0, max_lag, max_lag)
Models <- matrix(0, max_lag, max_lag)
for(p in 0:max_lag){
  for(q in 0:max_lag){
    fit.arma <- stats::arima(return, order = c(p,0,q),method = "ML")
    BICs[p,q] <- BIC(fit.arma)
    AICs[p,q] <- AIC(fit.arma)
    a=Box.test(fit.arma$residuals,lag=10,type = "Ljung-Box",fitdf=p+q)
    Models[p,q]<-a$p.value
  }
}
a <- which(BICs==min(BICs),arr.ind = TRUE)        ## ARMA(3,2)
b <- which(AICs==min(AICs), arr.ind = TRUE)       ## ARMA(10,4)
c <-  which(Models>0.05, arr.ind = TRUE)[1,]      ## ARMA(4,1)
eacf(return)


####自动定阶####
auto.arima(return, include.mean = TRUE)           ## ARMA(3,3)


### 使用BIC确定模型
par(mfrow = c(2, 5))
## AR模型：AR(4)
fit.ar <- arima(return, order = c(4, 0, 0))
summary(fit.ar)
BIC(fit.ar)
p_ar=fit.ar$coef/sqrt(diag(fit.ar$var.coef))
acf(fit.ar$residuals, xlab = "阶数", ylab = "自相关函数", main = "AR(4)")
pacf(fit.ar$residuals,  xlab = "阶数", ylab = "偏自相关函数", main = "AR(4)")
Box.test(fit.ar$residuals,lag = 10, type=c("Box-Pierce"))
Box.test(fit.ar$residuals,lag = 10, type=c("Ljung"))