ryu908

### IIPへの当てはめ
x<-read.table("clipboard",header=T)
ty1<-ts(x,start=c(1978,1),freq=12)
plot(ty1)
plot(decompose(ty1))
spectrum(ty1)
spectrum(ty1,method="ar")

par(mfrow=c(2,1))
acf(ty1,lag.max=100)
pacf(ty1,lag.max=100)
par(mfrow=c(1,1))

##### １階差分
dty1<-diff(ty1)
par(mfrow=c(3,1))
plot(dty1)
acf(dty1,lag.max=60)
pacf(dty1,lag.max=60)
par(mfrow=c(1,1))

sarima01<-arima(ty1,order=c(2,1,2),seasonal=list(order=c(1,1,1),period=12))

par(mfrow=c(2,1))
tres<-sarima01$resid
ty1_hat<-ty1-tres
plot(ty1)
lines(ty1_hat,lty=2,col=2)
plot(tres,type="l",col=4)
abline(h=mean(tres))
par(mfrow=c(1,1))


pred1<-predict(sarima01,n.ahead=12)
xt<-c(2005,2009)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)


pred1<-predict(sarima01,n.ahead=24)
xt<-c(1978,2009)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)