Skip to content

Instantly share code, notes, and snippets.

@timelyportfolio
Created January 12, 2012 21:22
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save timelyportfolio/1603209 to your computer and use it in GitHub Desktop.
stocks under extreme bonds
require(quantmod)
require(PerformanceAnalytics)
require(RQuantLib)
require(lattice)
getSymbols("AAA",src="FRED") # load Moody's AAA from Fed Fred
#Fed monthly series of yields is the monthly average of daily yields
#set index to yyyy-mm-dd format rather than to.monthly mmm yyy for better merging later
index(AAA)<-as.Date(index(AAA))
AAApricereturn<-AAA
AAApricereturn[1,1]<-0
colnames(AAApricereturn)<-"PriceReturn-monthly avg AAA"
#use quantlib to price the AAA and BAA bonds from monthly yields
#AAA and BAA series are 20-30 year bonds so will advance date by 25 years
for (i in 1:(NROW(AAA)-1)) {
AAApricereturn[i+1,1]<-FixedRateBondPriceByYield(yield=AAA[i+1,1]/100,issueDate=Sys.Date(),
maturityDate= advance("UnitedStates/GovernmentBond", Sys.Date(), 25, 3),
rates=AAA[i,1]/100,period=2)[1]/100-1
}
#total return will be the price return + yield/12 for one month
AAAtotalreturn<-AAApricereturn+lag(AAA,k=1)/12/100
colnames(AAAtotalreturn)<-"TotalReturn-monthly avg AAA"
AAAtotalreturn[1,1] <- 0
AAAcumul <- as.xts(apply(AAAtotalreturn+1,MARGIN=2,cumprod))
index(AAAcumul) <- as.Date(format(index(AAAcumul),"%Y-%m-31"))
#get DJ Industrial from FRED
getSymbols("DJIA",src="FRED")
roc <- na.omit(merge(annualReturn(DJIA),annualReturn(AAAcumul)))
colnames(roc) <- c("DJ Industrial","Moodys AAA")
par(mfrow=c(1,2))
plot(coredata(roc[,1])~coredata(roc[,2]),type="p",pch=20,
xlab=paste(colnames(roc)[2]," Annual Return (%)", sep=""),
ylab=paste(colnames(roc)[1]," Annual Return (%)", sep=""),
main="DJ Industrial and Moody's AAA
Annual Returns (%)",
cex.axis=0.85,
cex.lab=0.85)
abline(v=0.15,col="grey70")
abline(h=0,col="black")
text(x=roc[which(roc[,2] > 0.15),2],y=roc[which(roc[,2] > 0.15),1],
labels = format(index(roc[which(roc[,2] > 0.15)]),"%Y"),
pos=1, offset=-0.02,
cex = 0.7)
rect(xleft = 0.15,
xright = par()$usr[2],
ybottom = par()$usr[3],
ytop = par()$usr[4],
col = rgb(100,179,179,alpha=100,maxColorValue=255),
border = "gray70",
density = NA)
par(mar = c(3,4,5,2))
barplot(roc[which(roc[,2] > 0.15),],beside=TRUE,
names.arg = format(index(roc[which(roc[,2] > 0.15),]),"%Y"),
cex.names = 0.85, cex.axis = 0.85,
ylim=c(min(roc[which(roc[,2] > 0.15),]-.10),max(roc[which(roc[,2] > 0.15),]) + 0.05),
col = c("darkolivegreen4","slateblue4"),
main = "when Moody's AAA Annual Returns > 15%", cex.main = 0.85,
ylab = "Annual Return (%)", cex.lab = 0.85,
legend.text=TRUE, args.legend = list(x="top",bty="n",cex=0.85,horiz=TRUE))
abline(h=0,col="black")
colors = ifelse(roc[,1]<0,"red","green")
dotplot(factor(round(coredata(roc[,2]),1))~coredata(roc[,1]),
# type=c("p","a"),fun=median,
main="DJ Industrial Annual Returns
Grouped By Moody's AAA Annual Returns",
col=colors, xlab = "DJ Industrial Annual Return (%)",
ylab = "Moody's AAA Annual Return (rounded %)")
#trying to make negative half of graph red but this does not work
#grid.rect(x = unit(0.083,"npc"), y = unit(0.105,"npc"),
# just = c("left", "bottom"),
# width = unit(0.445,"npc"),height = unit(0.77,"npc"),
# draw=TRUE, gp=gpar(fill=rgb(1,0,0,0.3)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment