Rでout-of-sampleテスト（ローリング）

http://agecon.ucdavis.edu/graduate-program/class-websites/class/cid_426/inflation.1.do
をRに翻訳したもの

rm(list=ls(all=TRUE))

library(foreign)
library(dynlm)
library(sandwich)

phdata <- read.dta("inflation.dta")

phdata <- ts(phdata,freq=12,start=c(1959,1))

ln.cpiaucns <- log(phdata[,"cpiaucns"])
inf1 <- 1200*diff(ln.cpiaucns)
inf12 <- 100*diff(ln.cpiaucns,lag=12)

yvar <- inf12-lag(inf1,-12)

ln.indpro <- log(phdata[,"indpro"])
ip=1200*diff(ln.indpro)

ln.m2ns <- log(phdata[,"m2ns"])
m2 <- diff(ln.m2ns)

ln.ppicrm <- log(phdata[,"ppicrm"])
comm <- diff(ln.ppicrm)

ys <- phdata[,"gs10"] - phdata[,"tb3ms"]

unrate <- phdata[,"unrate"]

plot(ts.union(yvar,ip,m2,comm,ys,unrate))


# some basic regressions
reg1 <- dynlm(yvar~diff(lag(inf1,-12))+lag(unrate,-12),start=c(1960,1),end=c(1969,12))
summary(reg1)

reg2 <- dynlm(yvar~diff(lag(inf1,-12)),start=c(1960,1),end=c(1969,12))
summary(reg2)


# generate forecasts from first model
yhat <- NULL
e_beg <- c(1960,1)
e_end <- c(1969,12)
for (i in 1:(length(inf12)-12*10+1)) {
    reg <- dynlm(yvar~diff(lag(inf1,-12))+lag(unrate,-12),start=e_beg,end=e_end)
    yhat[i] <- tail(predict(reg),1)
    e_beg <- e_beg + c(0,1)
    e_end <- e_end + c(0,1)
}
yhat <- ts(yhat,freq=12,start=c(1970,12))
f_unemp <- yhat+lag(inf1,-12)


# generate forecasts from second model
yhat <- NULL
e_beg <- c(1960,1)
e_end <- c(1969,12)
for (i in 1:(length(inf12)-12*10+1)) {
    reg <- dynlm(yvar~diff(lag(inf1,-12)),start=e_beg,end=e_end)
    yhat[i] <- tail(predict(reg),1)
    e_beg <- e_beg + c(0,1)
    e_end <- e_end + c(0,1)
}
yhat <- ts(yhat,freq=12,start=c(1970,12))
f_linf <- yhat+lag(inf1,-12)
rm(yhat)


# forecast evaluation
e_unemp <- inf12 - f_unemp
e_linf <- inf12 - f_linf
e2 <- e_unemp^2 - e_linf^2

res_all <- dynlm(e2~1,start=c(1970,12),end=c(2011,8))
summary(res_all)
NeweyWest(res_all,lag=12)

res_7083 <- dynlm(e2~1,start=c(1970,12),end=c(1983,12))
summary(res_7083)
NeweyWest(res_7083,lag=12)

res_8498 <- dynlm(e2~1,start=c(1984,1),end=c(1998,12))
summary(res_8498)
NeweyWest(res_8498,lag=12)

res_9911 <- dynlm(e2~1,start=c(1999,1),end=c(2011,8))
summary(res_9911)
NeweyWest(res_9911,lag=12)