stocks under extreme bonds

stocks with extreme bonds.r

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