Create a gist now

Instantly share code, notes, and snippets.

#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
#install r-forge development version of bfast
#install.packages("bfast", repos="")
require(animation) #Yihui Xie surfaces again
#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",
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,
#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) <
#start monitoring from the max or min farthest away from current
mon5y <- bfastmonitor(GSPC.ts,
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment