public
Last active

  • Download Gist
bfast bull and bear.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
#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("^GSPC",from="1950-01-01")
#convert to log price
GSPC.monthly <- log(to.monthly(GSPC)[,4])
#get monthly returns for the close price
#not necessary, leave in price form
#GSPC.return <- monthlyReturn(GSPC[,4])
#need ts representation so do some brute force conversion
GSPC.ts <- ts(as.vector(GSPC.monthly["1951-01::"]),start=c(1951,1),frequency=12)
#look at the stl Seasonal-Trend decomposition procedure already in R
GSPC.stl <- stl(GSPC.ts,s.window="periodic")
plot(GSPC.stl,main="STL Decomposition of S&P 500")
 
 
#get the results from bfast
#adjusting h lower will result in more breakpoints
GSPC.bfast <- bfast(GSPC.ts,h=0.2,max.iter=1,season="none")
plot(GSPC.bfast,type="components",ylim=c(3,max(GSPC.monthly)+1),main="S&P 500 with bfast Breakpoints and Components")
plot(GSPC.bfast,type="trend",ylim=c(3,max(GSPC.monthly)+1),main="S&P 500 with bfast Trend Breakpoints")
 
 
#see everything with type="all" but in bfast calculation set seasonal to "none"
#play away with this
#plot(GSPC.bfast,type="all")
 
#do some additional plotting
#[[1]] is used since for speed I only did one iteration
#could plot each iteration if I did multiple
plot(GSPC.bfast$Yt/GSPC.bfast$output[[1]]$Tt-1,
main="bfast Remainder as % of S&P 500 Price",
xlab=NA, ylab="remainder (% of price)",bty="l")
#add vertical line for the breakpoints
abline(v=breakdates(GSPC.bfast$output[[1]]$bp.Vt),col="gray70")
#add horizontal line at 0
abline(h=0,col="black",lty=2)
text(x=breakdates(GSPC.bfast$output[[1]]$bp.Vt),y=par("usr")[3]+.01,
labels=breakdates(GSPC.bfast$output[[1]]$bp.Vt,format.times=TRUE),
srt=90,pos=4,cex=0.75)

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.