use of bfast in R to look at structural breakpoints in Retail / SPX relative strength

code.R

#analyze breakpoints with the R package bfast
#please read the paper 
#Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010)
#Detecting Trend and Seasonal Changes in Satellite Image Time Series.
#Remote Sensing of Environment, 114(1), 106–115.
#http://dx.doi.org/10.1016/j.rse.2009.08.014

require(bfast)
require(quantmod)

getSymbols("XRT",from="2000-01-01")
getSymbols("SPY",from="2000-01-01")
#get relative strenght as weekly log price
xrt_rs <- log(to.weekly(XRT[,6] / SPY[,6])[,4])

#need ts representation so do some brute force conversion
xrt_rs.ts <- ts(
  as.vector(xrt_rs),
  start=c(
    as.numeric(format(index(xrt_rs[1,]),"%Y")), #year
    as.numeric(format(index(xrt_rs[1,]),"%m"))  #month
  ),
  frequency=52
)
#look at the stl Seasonal-Trend decomposition procedure already in R
xrt_rs.stl <- stl(xrt_rs.ts,s.window="periodic")
plot(xrt_rs.stl,main="STL Decomposition of S&P 500")


#get the results from bfast
#adjusting h lower will result in more breakpoints
xrt_rs.bfast <- bfast(xrt_rs.ts,h=0.15,max.iter=1,season="none")
plot(
  xrt_rs.bfast,
  type="components",
  ylim=c(3,max(xrt_rs)+1),
  main="Retail Rel Strength with bfast Breakpoints and Components"
)
plot(
  xrt_rs.bfast,
  type="trend",
  ylim=c(3,max(xrt_rs)+1),
  main="Retail Rel Strength with bfast Trend Breakpoints"
)

#do some additional plotting
#[[1]] is used since for speed I only did one iteration
#could plot each iteration if I did multiple
plot(xrt_rs.bfast$Yt/xrt_rs.bfast$output[[1]]$Tt-1,
     main="bfast Remainder as % of Retail Rel Strength",
     xlab=NA, ylab="remainder (% of rel strength)",bty="l")
#add vertical line for the breakpoints
abline(v=breakdates(xrt_rs.bfast$output[[1]]$bp.Vt),col="gray70")
#add horizontal line at 0
abline(h=0,col="black",lty=2)
text(x=breakdates(xrt_rs.bfast$output[[1]]$bp.Vt),y=par("usr")[3]+.01,
     labels=breakdates(xrt_rs.bfast$output[[1]]$bp.Vt,format.times=TRUE),
     srt=90,pos=4,cex=0.75)