public

  • Download Gist
bfastmonitor experiments.r
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
#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"))))

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.