Forecasting with B-splines

extr2.R

# extrapolating into the future with B-splines
#  based on code in ?smooth.construct.bs.smooth.spec

library(mgcv)

# annual diameter of women’s skirts at the hem 1866 to 1911
# Hipel and McLeod, 1994
skirts <- scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat",skip=5)
skirtseries <- data.frame(year=1866:1911, diam=skirts)

# prediction grid
pd <- data.frame(year=seq(1864, 1960, by=1))

# plot data
plot(skirtseries$year, skirtseries$diam, xlim=c(1864, 1950),
     xlab="Year", ylab="Skirt diameter",
     main="black= bs k=25, m=3 \nred= bs k= 46, n =3\nblue= k=25, m = 2")

# bs with extended penalty
# setup knots to be ~
# c(start of series minus a bit, start of series, "today" (end), "the future")
b <- gam(diam~s(year,bs="bs",m=c(3,1),k=25),
         knots=list(year=c(1864,1865,1912,1960)),
         data=skirtseries,
         method="ML")

# predict and get approx intervals
fv <- predict(b,pd,se=TRUE)
ul <- fv$fit + fv$se.fit*2
ll <- fv$fit - fv$se.fit*2
# plot
lines(pd$year,fv$fit)
lines(pd$year,ul,lty=2)
lines(pd$year,ll,lty=2)


# How sensitive is this to the number of knots? 

# bs with penalty over the data range only
b0 <- gam(diam~s(year,bs="bs",m=c(3,1),k=46),
         knots=list(year=c(1864,1865,1912,1960)),
         data=skirtseries,
         method="ML")


# predict and get approx intervals
fv <- predict(b0,pd,se=TRUE)
ul <- fv$fit + fv$se.fit*2
ll <- fv$fit - fv$se.fit*2
# plot
lines(pd$year, fv$fit,col=2)
lines(pd$year, ul,lty=2,col=2)
lines(pd$year,ll,lty=2,col=2)



# What about the order of the smooth?
b1 <- gam(diam~s(year,bs="bs",m=c(1,1),k=25),
         knots=list(year=c(1864,1865,1912,1960)),
         data=skirtseries,
         method="ML")


fv <- predict(b1,pd,se=TRUE)
ul <- fv$fit + fv$se.fit*2
ll <- fv$fit - fv$se.fit*2

lines(pd$year,fv$fit,col=4)
lines(pd$year,ul,lty=2,col=4)
lines(pd$year,ll,lty=2,col=4)

the output: