|
#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) |