对未来5期CPI数据进行预测,RPI、PPI为解释变量,用R语言实现,很shallow,还请路过的大神指教
实现VAR模型对数据有较多要求,主要分为以下几步
1)平稳性检验,平稳最棒了,不平稳需要检查是否为同阶单整,不是同阶单整可以考虑对数差分化处理,使数据保持平稳,但同时也要注意差分以后的数据是否还具有意义。
2)对于平稳的数据可以做格兰杰因果检验,对于不平稳但同阶单整的数据可以做协整检验。格兰杰因并不等于格兰杰果。
3)确定模型的滞后阶数,不同的信息准则构造不同,但都是值越小越好,一般常用的有AIC、SC但选择时也要综合考虑,尽量选择使模型比较简介的吧(滞后阶数较少的)。
4)拟合模型
5)稳定性检验
6)脉冲响应分析
7)方差分解
8)预测

library(vars)
library(urca)
library(MASS)
library(sandwich)
library(tseries)
library(lmtest)
library(readxl)
CPI_DATA <- read_excel("D:/CPI DATA.xlsx")
data.ts<-ts(CPI_DATA[,2:4],start = c(2001,1),end = c(2019,5),frequency = 12)
plot.ts(data.ts,lty=1,type="o",main="2001.01-2019.05 CPI RPI PPI 走势图")
#legend("topright",legend = colnames(data.ts),col=c("green","red","blue"),pch = 17,cex = 0.8)#没必要加图例
#对数化处理
ldata.ts<-log(data.ts)
plot.ts(ldata.ts,type="o",main="对数时间序列图",)
#平稳性检验ur.df,adf.test都可用  ur.df函数中,none无,drift含常数项,trend含常数项和趋势项
#summary(ur.df(ldata.ts[,1],type="trend",selectlags = "AIC"))#接受原假设,非平稳
#summary(ur.df(ldata.ts[,2],type="none",selectlags = "AIC"))#接受原假设,非平稳
#summary(ur.df(ldata.ts[,3],type="none",selectlags = "AIC"))#拒绝原假设,平稳
adf.test(ldata.ts[,1])#不平稳
adf.test(ldata.ts[,2])#平稳
adf.test(ldata.ts[,3])#平稳
#平稳化处理
cpi.data<-diff(ldata.ts[,1])
plot.ts(cpi.data,lty=1,type="o",main="CPI变化率")
adf.test(cpi.data)对数差分化后CPI平稳
ndata.ts<-data.frame(cpi.data,ldata.ts[2:221,2],ldata.ts[2:221,3])#建立新数据
colnames(ndata.ts)<-c("CPI","RPI","PPI")
ndata.ts<-ts(ndata.ts,start = c(2001,2),end = c(2019,5),frequency = 12)
plot.ts(ndata.ts,type="o",main="时间序列图")
#jh.test<-ca.jo(ndata.ts,type = "trace",ecdet = "none",K=3,spec="transitory")#协整检验,因为数据是平稳序列,所以直接做格兰杰因果检验就OK
#summary(jh.test)
#格兰杰因果检验
grangertest(ndata.ts[,1]~ndata.ts[,2],order=3)
grangertest(ndata.ts[,1]~ndata.ts[,3],order=3)
grangertest(ndata.ts[,2]~ndata.ts[,1],order=3)
grangertest(ndata.ts[,3]~ndata.ts[,1],order=3)
#模型滞后阶数的确定
VARselect(ndata.ts,lag.max=10,type = "const")#选择最优滞后阶数
 #模型拟合
var.sc<-VAR(ndata.ts,p=2,lag.max = 10,ic="SC")
var.hq<-VAR(ndata.ts,p=3,lag.max = 10,ic="HQ")
var.fpe<-VAR(ndata.ts,p=4,lag.max = 10,ic="FPE")
var.aic<-VAR(ndata.ts,p=10,lag.max = 10,ic="AIC")
plot(var.sc)
plot(var.hq)
plot(var.fpe)
plot(var.aic)
summary(var.sc)
#causality(var,cause = c('RPI','PPI')#拟合模型前和拟合模型后做格兰杰因果检验是有差别的
#稳健性检验
sta.sc<-stability(var.sc,type = "OLS-CUSUM",h=0.15)
sta.hq<-stability(var.hq,type = "OLS-CUSUM",h=0.15)
sta.fpe<-stability(var.fpe,type = "OLS-CUSUM",h=0.15)
sta.aic<-stability(var.aic,type = "OLS-CUSUM",h=0.15)
plot(sta.sc)
plot(sta.hq)
plot(sta.fpe)
plot(sta.aic)
#脉冲响应分析
var.irf<-irf(var.sc,n.ahead = 30)
plot(var.irf)
#方差分解
var.fevd<-fevd(var.sc,n.ahead = 10)
plot(var.fevd)
var.fevd$CPI
var.fevd$RPI
var.fevd$PPI
#预测
var.predict<-predict(var.sc,n.ahead = 5,ci=0.95)
plot(var.predict)
fanchart(var.predict)

神经网络算法代码

library(neuralnet)
library(nnet)
library(AMORE)
library(lattice)
library(grid)
library(DMwR)
library(readxl)
CPI_DATA <- read_excel("D:/CPI DATA.xlsx")
#将输入数据归一化
#maxs<-apply(CPI_DATA,2,max)
#mins<-apply(CPI_DATA,2,min)
#bp.data<-as.data.frame(scale(CPI_DATA,center = mins,scale = maxs-mins))
#传统方法
cmax<-c(max(CPI_DATA[,2]),max(CPI_DATA[,3]),max(CPI_DATA[,4]))
cmin<-c(min(CPI_DATA[,2]),min(CPI_DATA[,3]),min(CPI_DATA[,4]))
bp.data<-data.frame((CPI_DATA[,2]-cmin[1])/(cmax[1]-cmin[1]),(CPI_DATA[,3]-cmin[2])/(cmax[2]-cmin[2]),(CPI_DATA[,4]-cmin[3])/(cmax[3]-cmin[3]))
colnames(bp.data)<-c("CPI","RPI","PPI")
#分组
train<-bp.data[1:216,]
test<-bp.data[217:221,]
#没找到可视化功能,所以放弃了
#net<-newff(n.neurons = c(2,7,5,1),learning.rate.global = 1e-3,momentum.global = 0.4,error.criterium = 'LMS',Stao = NA,hidden.layer = "sigmoid",output.layer = 'purelin',method="ADAPTgdwm")
#result<-train(net,train[,2:3],train[,1],error.criterium="LMS",report=TRUE,show.step=1000,n.shows=5)
#y<-sim(result$net,test[,2:3])
#org<-y*(cmax[1]-cmin[1])+cmin[1]#逆归一化,有争议不好用
#以下都是寻找最优隐藏层
net1<-neuralnet(CPI~RPI+PPI,data = train,hidden = c(5,3),threshold = 0.05,learningrate = 0.1,algorithm = "rprop+",err.fct = "sse",act.fct = 'tanh',linear.output = T)
plot(net1)
net.result1<-compute(net1,test[,2:3])
org1<-net.result1$net.result*(cmax[1]-cmin[1])+cmin[1]
org1
net2<-neuralnet(CPI~RPI+PPI,data = train,hidden = c(7,5),threshold = 0.05,learningrate = 0.1,algorithm = "rprop+",err.fct = "sse",act.fct = 'tanh',linear.output = T)
plot(net2)
org2<-net.result2$net.result*(cmax[1]-cmin[1])+cmin[1]
org2
net3<-neuralnet(CPI~RPI+PPI,data = train,hidden = c(7,4),threshold = 0.05,learningrate = 0.1,algorithm = "rprop+",err.fct = "sse",act.fct = 'tanh',linear.output = T)
plot(net3)
net.result2<-compute(net2,test[,2:3])
org3<-net.result2$net.result*(cmax[1]-cmin[1])+cmin[1]
org3