statistics of drawdowns.r

require(fGarch)
require(ttrTests)
require(quantmod)
require(PerformanceAnalytics)

getSymbols("^GSPC",from="1900-01-01")
#get monthly close prices from daily
price.monthly<-to.monthly(GSPC)[,4]
#get dates in yyyy-mm-01 format
index(price.monthly) <- as.Date(index(price.monthly))
#get arithmetic one month rate of change
roc.monthly <-ROC(price.monthly,n=1,type="discrete")
#change NA from first month to 0
roc.monthly[1,] <- 0
#get vector of prices for use with ttrTests
price.vector<-as.vector(price.monthly)

#fit historical to skewed student t 
sstd.fit <- sstdFit(as.vector(roc.monthly))
#fit historical to student t
std.fit <- stdFit(as.vector(roc.monthly))

#this takes a while; choose smaller for playing
bsamples=10000
#prepopulate drawdowns matrix with NA
#for quicker filling later
drawdowns <- matrix(nrow=bsamples,ncol=3)
returns <- drawdowns
#will generate samples for bootstrap and stationary bootstrap
#to evaluate the maxDrawdown differences between the distributions
colnames(drawdowns) <- c("student.t","skew.student.t","ttrTests.stat.boot")
colnames(returns) <- colnames(drawdowns)

for (i in 1:bsamples) {
  std.sample <- as.xts(rstd(n=NROW(roc.monthly),mean=std.fit$par[1],sd=std.fit$par[2],nu=std.fit$par[3]),
                       order.by=index(roc.monthly))
  sstd.sample <- as.xts(rsstd(n=NROW(roc.monthly),mean=sstd.fit$estimate[1],sd=sstd.fit$estimate[2],nu=sstd.fit$estimate[3],xi=sstd.fit$estimate[4]),
                        order.by=index(roc.monthly))
  boot.sample <- ROC(
                  as.xts(generateSample(price.vector,"stationaryBootstrap"),
                        order.by=index(price.monthly)),
                  n=1,type="discrete")
  drawdowns[i,1] <- maxDrawdown(std.sample)
  returns[i,1] <- Return.cumulative(std.sample)
  drawdowns[i,2] <- maxDrawdown(sstd.sample)
  returns[i,2] <- Return.cumulative(sstd.sample)
  drawdowns[i,3] <- maxDrawdown(boot.sample)
  returns[i,3] <- Return.cumulative(boot.sample)
}

#do the density plot for drawdown
d1 <- density(drawdowns[,1])
d2 <- density(drawdowns[,2])
d3 <- density(drawdowns[,3])
plot( d1, col=2, lwd=3, main="Density Plot of Drawdown by Distribution Method")
lines( d2, col=4, lwd=3)
lines( d3, col=3, lwd=3)
abline(v=maxDrawdown(ROC(price.monthly,type="discrete",n=1)),col="grey70")
#label experienced drawdown for the historical price series
text(x=maxDrawdown(ROC(price.monthly,type="discrete",n=1)), y=3.7, pos=3,
     labels="SP500",srt=90,col="grey70", cex=0.75)
legend("topright",legend=colnames(drawdowns),col=c(2,4,3),lwd=3,bty="n")

#do the density plot for cumulative return
d1 <- density(returns[,1])
d2 <- density(returns[,2])
d3 <- density(returns[,3])
plot( d1, col=2, lwd=3, main="Density Plot of Cumulative Return by Distribution Method")
lines( d2, col=4, lwd=3)
lines( d3, col=3, lwd=3)
abline(v=Return.cumulative(ROC(price.monthly,type="discrete",n=1)),col="grey70")
#label experienced drawdown for the historical price series
text(x=Return.cumulative(ROC(price.monthly,type="discrete",n=1)), y=3.7, pos=3,
     labels="SP500",srt=90,col="grey70", cex=0.75)
legend("topright",legend=colnames(returns),col=c(2,4,3),lwd=3,bty="n")

#do the density plot for annualized cumulative return
d1 <- density((returns[,1]+1)^(1/(NROW(price.monthly)/12))-1)
d2 <- density((returns[,2]+1)^(1/(NROW(price.monthly)/12))-1)
d3 <- density((returns[,3]+1)^(1/(NROW(price.monthly)/12))-1)
plot( d1, col=2, lwd=3, main="Density Plot of Cumulative Ann. Return by Distribution Method")
lines( d2, col=4, lwd=3)
lines( d3, col=3, lwd=3)
abline(v=Return.cumulative(ROC(price.monthly,type="discrete",n=1)),col="grey70")
#label experienced drawdown for the historical price series
text(x=Return.cumulative(ROC(price.monthly,type="discrete",n=1)), y=3.7, pos=3,
     labels="SP500",srt=90,col="grey70", cex=0.75)
legend("topright",legend=colnames(returns),col=c(2,4,3),lwd=3,bty="n")

#plot Annualized Return vs Drawdown
plot((((returns[,1]+1)^(1/(NROW(price.monthly)/12))-1)~drawdowns[,1]),col=2,
     main="Annualized Return and Max Drawdown",ylab="Annualized Return",xlab="Max Drawdown")
points(((returns[,2]+1)^(1/(NROW(price.monthly)/12))-1)~drawdowns[,2],col=4)
points(((returns[,3]+1)^(1/(NROW(price.monthly)/12))-1)~drawdowns[,3],col=3)
legend("topright",legend=colnames(returns),col=c(2,4,3),lwd=3,bty="n")