扫码加入数据分析学习群

library(tseries)         #加载算法包tseries
library(forecast)        #加载算法包forecast
YEAR <- 2012            #设置初始时间年份
Frequency<- 12          #设置数据类型（月度数据为12，季度数据为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)
#读入数据源文件
Data<-ts(data$CPI, freq = Frequency, start = c(YEAR))
#acf(data)                #查看自相关图
#pacf(data)               #查看偏自相关图
model <-arima(Data,order=c(p,q,d),seasonal=list(order=c(P,Q,D),period= Frequency), include.mean= FALSE, method="ML")
#调用ARIMA模型
#其中参数：
# ① p、q、d、P、Q、D对应的是非季节性参数和季节性参数；
# ② include.mean表示是否带常数项，TRUE为带常数项(默认项，为非差分模型)，FALSE为不带常数项。
pred <-forecast(model, h = step, level=c(95))
Fit<- data.frame(pred$fitted)              #模型拟合值
Error<-data.frame(model$residuals)      #导出残差
#tsdiag(model)                         #R中残差及相关检验图
Fittable<-cbind(data[,1],data$CPI,Fit,Error)    #拟合表
names(Fittable)<-c('Time','CPI','CPIFitted','Residual')
#write.table(Fittable,'C:/Users/Administrator/Desktop/Fittable.txt')
#输出拟合残差结果表
Forecast<-data.frame(cbind(predmean,predupper,pred$lower))
#预测值、预测值95%上限、预测值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')