library(openxlsx)
library(forecast)

data1<-read.xlsx("http://kanggc.iptime.org/time/R/djtinx.xlsx")
number<-data1$number
djtinx<-data1$djtinx
djtinx.ts<-ts(djtinx, start=c(2013,1), frequency=12)

par(mfrow=c(3,1))
ts.plot(djtinx.ts, ylab="djtinx")
acf(djtinx.ts)
pacf(djtinx.ts)

(ddjtinx<-djtinx.ts-lag(djtinx.ts,k=-1))
ts.plot(ddjtinx, ylab="ddjtinx")
acf(ddjtinx)
pacf(ddjtinx)

arma11=arima(ddjtinx, order=c(1,0,1), include.mean=F)
arma11
arma12=arima(ddjtinx, order=c(1,0,2), include.mean=F)
arma21=arima(ddjtinx, order=c(2,0,1), include.mean=F)
AIC(arma11,arma12,arma21, k = 2)
BIC(arma11,arma12,arma21)

res<-residuals(arma11)
ts.plot(res)
acf(res)
pacf(res)

Box.test(arma11$resid,type="Ljung-Box")

par(mfrow=c(1,1))
plot(ddjtinx,lty=1,main="ddjtinx: raw data vs fitted values",ylab="djtinx",xlab="date")
lines(fitted(arma11),lty=2,lwd=1,col="red")
predict(arma11, n.ahead = 6)
(f3<-forecast(arma11, h=6))
plot(f3)

arma11_n=arima(ddjtinx[1:59], order=c(1,0,1), include.mean=F)
res1<-residuals(arma11_n)
res1
(mean(res1))
(mean(abs(res1)))
f11<-forecast(arma11_n, h = 5)
accuracy(f11)
accuracy(f11,ddjtinx[60:64])