Created
May 29, 2020 06:44
-
-
Save dill/20fbb893c1e1194eb462843156115106 to your computer and use it in GitHub Desktop.
Forecasting with B-splines
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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 with extended penalty\nred=regular bs\nblue=tprs") | |
# 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) | |
# other models... | |
# bs with penalty over the data range only | |
b0 <- gam(diam~s(year,bs="bs",m=c(3,1),k=25),method="ML", data=skirtseries) | |
# 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 thin plate? | |
b1 <- gam(diam~s(year,bs="tp",k=25),method="ML", data=skirtseries) | |
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) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Resulting plot:
(Probably not a brilliant example, but gives some indication of how to do it.)