bfastmonitor experiments.r

#analyze breakpoints in realtime with the R package bfast
#please read the paper 
#Jan Verbesselt, Achim Zeileis, Martin Herold (2011). Near Real-Time Disturbance
#Detection in Terrestrial Ecosystems Using Satellite Image Time Series: Drought
#Detection in Somalia. Working Paper 2011-18. Working Papers in Economics and
#Statistics, Research Platform Empirical and Experimental Economics, Universitaet
#Innsbruck. URL http://EconPapers.RePEc.org/RePEc:inn:wpaper:2011-18
#http://scholar.google.com/scholar?cluster=9016488513865299942&hl=en&as_sdt=0,1

#install r-forge development version of bfast
#install.packages("bfast", repos="http://R-Forge.R-project.org")

require(bfast)
require(quantmod)
require(animation) #Yihui Xie surfaces again

getSymbols("^GSPC",from="1950-01-01")
#convert to log price
#for the sake of this example, let's assume we are in January 2009
#and we will see how this progresses with a for loop to go through
#the painful months of late 2008 and early 2009
#to see when our real time monitoring detects a break

#get a vector of the dates
evaldates <- c("2008-09","2008-10","2008-11","2008-12",
               "2009-01","2009-02","2009-03","2009-04",
               "2009-05","2009-06")

saveGIF(
for(i in 1:length(evaldates)) {

#notice this removes the foresight bias by only going to the current month
GSPC.monthly <- log(to.monthly(GSPC)[paste("1990::",evaldates[i],sep=""),4])

#need ts representation so do some brute force conversion
GSPC.ts <- ts(as.vector(GSPC.monthly),start=as.numeric(c(format(index(GSPC.monthly)[1],"%Y"),format(index(GSPC.monthly)[1],"%m"))),frequency=12)

#since we know the peak was October of 2010 let's use that
#as our start point for real time monitoring
mon5y <- bfastmonitor(GSPC.ts,
                      start=c(2007,10))
plot(mon5y)
mtext(evaldates[i],col="green",cex=1.5)
}
)


#really does not work, but always nice to have more examples


#let's get the minimum and maximum for our first experiment
GSPC.max <- GSPC.monthly[which(GSPC.monthly==max(last(GSPC.monthly,"10 years")))]
GSPC.min <- GSPC.monthly[which(GSPC.monthly==min(last(GSPC.monthly,"10 years")))]
#for evaluation let's pick whatever point is farthest from most recent close
GSPC.eval <- GSPC.monthly[GSPC.monthly == ifelse(as.numeric(last(GSPC.monthly))-as.numeric(GSPC.min) < 
  abs(as.numeric(last(GSPC.monthly)-as.numeric(GSPC.max))),GSPC.max,GSPC.min)]


#start monitoring from the max or min farthest away from current
mon5y <- bfastmonitor(GSPC.ts,
                      start=as.numeric(c(format(index(GSPC.eval),"%Y"),
                                         format(index(GSPC.eval),"%m"))))