##The case study of linear time series analysis
#汽油和原油的回归加入到时间序列中，探究汽油的时间序列和对汽油价格增长率的预测。
#汽油和原油的回归加入到时间序列中，探究汽油的时间序列和对汽油价格增长率的预测。
```{r}
da<-read.table("w-petroprice.txt",header = T)
da1<-read.table("w-gasoline.txt")
pgs<-log(da1[,1])
pus<-log(da$US)
tdx<-c(1:717)/52+1997
par(mfcol=c(2,1))
plot(tdx,pgs,xlab = "year",ylab = "ln(price)",type = "l")
title(main = "(a)gasoline")
plot(tdx,pus,xlab = "year",ylab = "ln(price)",type = "l")
title(main = "(b)petroprice")
#gas纯时间序列
dpg<-diff(pgs)
acf(dpg,lags=20)
pacf(dpg,lags=20)#acf拖尾，pacf截尾，使用AR(p)模型
m1<-ar(dpg,method = "mle")
m1#通过系数判断模型的平稳性和阶数
t.test(dpg)#t检验，检验acf值，是否系数为0显著。
m1<-arima(dpg,order = c(5,0,0),include.mean = F)
m1
#调参
m1<-arima(dpg,order= c(5,0,0),include.mean = F,fixed =c(NA,NA,NA,0,NA))
m1
tsdiag(m1,gof=20)#拟合优度检验
#建立原油和汽油之间的回归方程
dus<-diff(pus)
m3<-lm(dpg~-1+dus)#表示探究回归系数：t期，因变量dps,自变量dus之间的系数关系
summary(m3)
#探究回归中的残差符合的模型
acf(m3$residuals,lag=20)
pacf(m3$residuals,lag=20)#acf拖尾，pacf截尾
m4<-ar(m3$residuals,method = "mle")
m4$order
m4<-arima(m3$residuals,order=c(6,0,0),include.mean = F,xreg = dus)###加入回归自变量影响的残差模型。
m4
#调参,因为第6阶不显著，所以尝试arima，5阶，看AIC的值变化情况
m4<-arima(m3$residuals,order=c(5,0,0),include.mean = F,xreg = dus)
m4
m4<-arima(m3$residuals,order=c(5,0,0),include.mean = F,xreg = dus,fixed = c(NA,NA,NA,0,NA,NA))
m4
tsdiag(m4,gof=20)
#从AIC的减小可以看出，回归模型给出了更精确的样本外预测。
#回测
source("backtest.R")
pm1<-backtest(m1,dpg,316,1,fixed=c(NA,NA,NA,0,NA),inc.mean=F)#用316个数据点去预测317-716的数据，用观察值-预测值
pm4<-backtest(m4,dpg,316,1,xreg=dus,inc.mean=F,fixed=c(NA,NA,NA,0,NA,NA))
tdx<-tdx[2:717]
pm4fit<-dpg[317:716]-pm4$error
pm1fit<-dpg[317:716]-pm1$error
plot(tdx[317:716],dpg[317:716],xlab = "year",ylab = "growth",type = "l")
points(tdx[317:716],pm1fit,pch="*")
plot(tdx[317:716],dpg[317:716],xlab = "year",ylab = "growth",type = "l")#实际的图形
points(tdx[317:716],pm4fit,pch="*")#预测的图形
#Use lagged growrh rate of crude oil price
m6<-lm(dpg[2:716]~-1+dus[1:715])
summary(m6)
par(mfcol=c(2,1))
acf(m6$residuals,lag=20)
pacf(m6$residuals,lag=20)
m7<-ar(m6$residuals,method = "mle")
m7$order
m7<-arima(dpg[2:716],order = c(9,0,0),xreg = dus[1:715],include.mean = F)
m7
m7<-arima(dpg[2:716],order = c(9,0,0),xreg = dus[1:715],include.mean = F,fixed = c(NA,NA,NA,0,NA,0,0,0,NA,NA))
m7
tsdiag(m7,gof=20)
pm7<-backtest(m7,dpg[2:716],315,1,xreg = dus[1:715],inc.mean = F,fixed = c(NA,NA,NA,0,NA,0,0,0,NA,NA))
m6<-lm(dpg[2:716]~-1+dus[1:715])
summary(m6)
par(mfcol=c(2,1))
acf(m6$residuals,lag=20)
pacf(m6$residuals,lag=20)
m7<-ar(m6$residuals,method = "mle")
m7$order
m7<-arima(dpg[2:716],order = c(9,0,0),xreg = dus[1:715],include.mean = F)
m7
m7<-arima(dpg[2:716],order = c(9,0,0),xreg = dus[1:715],include.mean = F,fixed = c(NA,NA,NA,0,NA,0,0,0,NA,NA))
m7
tsdiag(m7,gof=20)
pm7<-backtest(m7,dpg[2:716],315,1,xreg = dus[1:715],inc.mean = F,fixed = c(NA,NA,NA,0,NA,0,0,0,NA,NA))