Skip to content

Instantly share code, notes, and snippets.

@dtsmith2001

dtsmith2001/bfast bull and bear.r

Forked from timelyportfolio/bfast bull and bear.r
Created February 3, 2016 14:23
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save dtsmith2001/5d3e044ffee71282d5c3 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save dtsmith2001/5d3e044ffee71282d5c3 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. @timelyportfolio timelyportfolio created this gist Apr 26, 2012.
    47 changes: 47 additions & 0 deletions bfast bull and bear.r
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,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)