library(tseries)? ?? ?? ?#加載算法包tseries
library(forecast)? ?? ???#加載算法包forecast
YEAR <- 2012? ?? ?? ?? ?#設(shè)置初始時(shí)間年份
Frequency<- 12? ?? ?? ? #設(shè)置數(shù)據(jù)類型（月度數(shù)據(jù)為12，季度數(shù)據(jù)為4）
p<-1
q<-1
d<-0
P<-2
Q<-0
D<-1
step<-4
data <-read.table('C:/Users/Administrator/Desktop/cpi.txt',head = T)
#讀入數(shù)據(jù)源文件
Data<-ts(data$CPI, freq = Frequency, start = c(YEAR))
#acf(data)? ?? ?? ?? ?? ? #查看自相關(guān)圖
#pacf(data)? ?? ?? ?? ?? ?#查看偏自相關(guān)圖
model <-arima(Data,order=c(p,q,d),seasonal=list(order=c(P,Q,D),period= Frequency), include.mean= FALSE, method="ML")
#調(diào)用ARIMA模型
#其中參數(shù)：
# ① p、q、d、P、Q、D對(duì)應(yīng)的是非季節(jié)性參數(shù)和季節(jié)性參數(shù)；
# ② include.mean表示是否帶常數(shù)項(xiàng)，TRUE為帶常數(shù)項(xiàng)(默認(rèn)項(xiàng)，為非差分模型)，F(xiàn)ALSE為不帶常數(shù)項(xiàng)。
pred <-forecast(model, h = step, level=c(95))
Fit<- data.frame(pred$fitted)? ?? ?? ?? ???#模型擬合值
Error<-data.frame(model$residuals)? ?? ?#導(dǎo)出殘差
#tsdiag(model)? ?? ?? ?? ?? ?? ?? ?? ? #R中殘差及相關(guān)檢驗(yàn)圖
Fittable<-cbind(data[,1],data$CPI,Fit,Error)? ? #擬合表
names(Fittable)<-c('Time','CPI','CPIFitted','Residual')
#write.table(Fittable,'C:/Users/Administrator/Desktop/Fittable.txt')
#輸出擬合殘差結(jié)果表
Forecast<-data.frame(cbind(predmean,predupper,pred$lower))
#預(yù)測(cè)值、預(yù)測(cè)值95%上限、預(yù)測(cè)值95%下限
names(Forecast)<-c('forecast','upper','lower')
lastDt <- data[nrow(data),1]
if(nchar(lastDt)==6){
lastDt <-as.Date(paste(as.character(lastDt),'01'),format='%Y%m%d')
forecastDt<- seq(lastDt,length.out=step+1,by='1 month')
forecastDt<- format(forecastDt,'%Y%m')
monthkey <-as.numeric(forecastDt)
forecast<- data.frame(monthkey=monthkey[2:length(monthkey)],Forecast)
}
if (nchar(lastDt)==4){
lastDt<-as.numeric(lastDt)
forecastDt<- seq(lastDt,lastDt+step,by=1)
Year<- forecastDt
forecast<- data.frame(Year=Year[2:length(Year)],Forecast)
}
write.table(forecast,'C:/Users/Administrator/Desktop/forecast.txt')