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"))