Friday, April 29, 2011

Another Use of LSPM in Tactical Portfolio Allocation

After the slightly unconventional use of LSPM presented in Slightly Different Use of Ralph Vince’s Leverage Space Trading Model, I thought I should follow up with something that more closely resembles my interpretation of Ralph Vince’s book.

LSPM seems to work well for portfolio allocation problems.  In this tactical allocation system, I will use optimal f derived in R with the LSPM package to build a portfolio with SP500, US10y, and the CRB.  If we use the optimal f as our allocation to SP500 and CRB, then the results look like this.

From TimelyPortfolio
From TimelyPortfolio

In general, my biggest problem is applying my systems to an entire portfolio.  The components are easy, but the blending is troublesome.  If I apply a basic method of monthly rebalancing of the SP500 and CRB, I get something like this (not using any leverage).

From TimelyPortfolio

SP500 component improves, but because of my 50% limit on the more mean reverting CRB, the CRB component underperforms the straight optimal f allocation.

From TimelyPortfolio

Since most clients don’t like cash, we can fill the portfolio with bonds when there is room left over after the CRB and SP500 allocation.  The total package look like this.

From TimelyPortfolio

As always, I would like my posts to stimulate discussion and thought.  The drawdown here is much more severe than I would like.  Please let me know how you would improve this system.

R code:

#Please see au.tra.sy blog http://www.automated-trading-system.com/
#for original walkforward/optimize code and http://www.fosstrading.com
#for other techniques

require(PerformanceAnalytics)
require(quantmod)
require(RQuantLib)

#get bond returns to avoid proprietary data problems
#see previous timelyportfolio blogposts for explanation
#probably need to make this a function since I will be using so much
getSymbols("GS10",src="FRED") #load US Treasury 10y from Fed Fred

GS10pricereturn<-GS10  #set this up to hold price returns

GS10pricereturn[1,1]<-0
colnames(GS10pricereturn)<-"PriceReturn"
#I know I need to vectorize this but not qualified enough yet
#Please feel free to comment to show me how to do this
for (i in 1:(NROW(GS10)-1)) {
  GS10pricereturn[i+1,1]<-FixedRateBondPriceByYield(yield=GS10[i+1,1]/100,issueDate=Sys.Date(),
    maturityDate= advance("UnitedStates/GovernmentBond", Sys.Date(), 10, 3),
    rates=GS10[i,1]/100,period=2)[1]/100-1
}

#interest return will be yield/12 for one month
GS10interestreturn<-lag(GS10,k=1)/12/100
colnames(GS10interestreturn)<-"Interest Return"

#total return will be the price return + interest return
GS10totalreturn<-GS10pricereturn+GS10interestreturn
colnames(GS10totalreturn)<-"Bond Total Return"

#get sp500 returns from FRED
getSymbols("SP500",src="FRED") #load SP500

#unfortunately cannot get substitute for proprietary CRB data
#get data series from csv file
CRB<-as.xts(read.csv("spxcrbndrbond.csv",row.names=1))[,2]

#do a little manipulation to get the data lined up on monthly basis
GS10totalreturn<-to.monthly(GS10totalreturn)[,4]
SP500<-to.monthly(SP500)[,4]
#get monthly format to yyyy-mm-dd with the first day of the month
index(SP500)<-as.Date(index(SP500))
#my CRB data is end of month; could change but more fun to do in R
CRB<-to.monthly(CRB)[,4]
index(CRB)<-as.Date(index(CRB))

#now lets merge to get asset class returns
assetROC<-na.omit(merge(ROC(SP500,type="discrete"),CRB,GS10totalreturn))

# Set Walk-Forward parameters (number of periods)
optim<-12 #1 year = 12 monthly returns
wf<-1 #walk forward 1 monthly returns

numsys<-2
# Calculate number of WF cycles
numCycles = floor((nrow(assetROC)-optim)/wf)
# Define JPT function
# this is now part of LSPM package, but fails when no negative returns
# so I still include this where I can force a negative return
jointProbTable <- function(x, n=3, FUN=median, ...) {
  # Load LSPM
  if(!require(LSPM,quietly=TRUE)) stop(warnings())

  # handle case with no negative returns
  for (sys in 1:numsys) {
    if (min(x[,sys])> -1) x[,sys][which.min(x[,sys])]<- -0.03
  }
  # Function to bin data
  quantize <- function(x, n, FUN=median, ...) {
    if(is.character(FUN)) FUN <- get(FUN)
    bins <- cut(x, n, labels=FALSE)
    res <- sapply(1:NROW(x), function(i) FUN(x[bins==bins[i]], ...))
  }
  # Allow for different values of 'n' for each system in 'x'
  if(NROW(n)==1) {
    n <- rep(n,NCOL(x))
  } else
  if(NROW(n)!=NCOL(x)) stop("invalid 'n'")
  # Bin data in 'x'
  qd <- sapply(1:NCOL(x), function(i) quantize(x[,i],n=n[i],FUN=FUN,...))
  # Aggregate probabilities
  probs <- rep(1/NROW(x),NROW(x))
  res <- aggregate(probs, by=lapply(1:NCOL(qd), function(i) qd[,i]), sum)
  # Clean up output, return lsp object
  colnames(res) <- colnames(x)
  res <- lsp(res[,1:NCOL(x)],res[,NCOL(res)])
  return(res)
}

for (i in 0:(numCycles-1)) {
            # Define cycle boundaries
            start<-1+(i*wf)
            end<-optim+(i*wf)
            # Get returns for optimization cycle and create the JPT
        # specify number of bins; does not seem to drastically affect results
        numbins<-6
            jpt <- jointProbTable(assetROC[start:end,1:numsys],n=rep(numbins,numsys))
            outcomes<-jpt[[1]]
            probs<-jpt[[2]]
            port<-lsp(outcomes,probs)
            # DEoptim parameters (see ?DEoptim)
            np=numsys*10       # 10 * number of mktsys
            imax=1000       #maximum number of iterations
            crossover=0.6       #probability of crossover
            NR <- NROW(port$f)
            DEctrl <- list(NP=np, itermax=imax, CR=crossover, trace=TRUE)
            # Optimize f
            res <- optimalf(port, control=DEctrl)
        # use upper to restrict to a level that you might feel comfortable
            #res <- optimalf(port, control=DEctrl, lower=rep(0,13), upper=rep(0.2,13))

    # these are other possibilities but I gave up after 24 hours
        #maxProbProfit from Foss Trading
        #res<-maxProbProfit(port, 1e-6, 6, probDrawdown, 0.1, DD=0.2, control=DEctrl)
        #probDrawdown from Foss Trading
        #res<-optimalf(port,probDrawdown,0.1,DD=0.2,horizon=6,control=DEctrl)

            # Save leverage amounts as optimal f
        # Examples in the book Ralph Vince Leverage Space Trading Model
        # all in dollar terms which confuses me
        # until I resolve I changed lev line to show optimal f output
        lev<-res$f[1:numsys]
            #lev<-c(res$f[1]/(-jpt$maxLoss[1]/10),res$f[2]/(-jpt$maxLoss[2]/10))
            levmat<-c(rep(1,wf)) %o% lev #so that we can multiply with the wfassetROC
            # Get the returns for the next Walk-Forward period
            wfassetROC <- assetROC[(end+1):(end+wf),1:numsys]
            wflevassetROC <- wfassetROC*levmat #apply leverage to the returns
            if (i==0) fullassetROC<-wflevassetROC else fullassetROC<-rbind(fullassetROC,wflevassetROC)
        if (i==0) levered<-levmat else levered<-rbind(levered,levmat)
}

#not super familiar with xts, but this add dates to levered
levered<-xts(levered,order.by=index(fullassetROC) )
colnames(levered)<-c("sp500 optimal f","crb optimal f")
chart.TimeSeries(levered, legend.loc="topleft", cex.legend=0.6)

#review the optimal f values
#I had to fill the window to my screen to avoid a error from R on margins
par(mfrow=c(numsys,1))
for (i in 1:numsys) {
    chart.TimeSeries(levered[,i],xlab=NULL)
}

#charts.PerformanceSummary(fullassetROC, ylog=TRUE, main="Performance Summary with Optimal f Applied as Allocation")

assetROCAnalyze<-merge(assetROC,fullassetROC)
colnames(assetROCAnalyze)<-c("sp500","crb","US10y","sp500 f","crb f")
charts.PerformanceSummary(assetROCAnalyze,ylog=TRUE,main="Performance Summary with Optimal f Applied as Allocation")

#build a portfolio with sp500 and crb
leveredadjust<-levered
#allow up to 50% allocation in CRB
leveredadjust[,2]<-ifelse(levered[,2]<0.25,0,0.5)
#allow up to 100% allocation in SP500 but portfolio constrained to 1 leverage
leveredadjust[,1]<-ifelse(levered[,1]<0.25,0,1-levered[,2])
colnames(leveredadjust)<-c("sp500 portfolio allocation","crb portfolio allocation")
assetROCadjust<-merge(assetROCAnalyze,leveredadjust[,1:2]*assetROC[,1:2])
colnames(assetROCadjust)<-c("sp500","crb","US10y","sp500 f","crb f","sp500 system component","crb system component")
charts.PerformanceSummary(assetROCadjust,ylog=TRUE)

#review the allocations versus optimal f
#I had to fill the window to my screen to avoid a error from R on margins
par(mfrow=c(numsys,1))
for (i in 1:numsys) {
    chart.TimeSeries(merge(levered[,i],leveredadjust[,i]),xlab=NULL,legend.loc="topleft",main="")
}

#add bonds when out of sp500 or crb
assetROCportfolio<-assetROCadjust[,6]+assetROCadjust[,7]+ifelse(leveredadjust[,1]+leveredadjust[,2] >= 1,0,(1-leveredadjust[,1]-leveredadjust[,2])*assetROC[,3])
assetROCadjust<-merge(assetROCAnalyze,assetROCportfolio)
colnames(assetROCadjust)<-c("sp500","crb","US10y","sp500 f","crb f","system portfolio")
charts.PerformanceSummary(assetROCadjust[,c(1:3,6)],ylog=TRUE,main="Optimal f System Portfolio with Bond Filler")

Thursday, April 28, 2011

Slightly Different Use of Ralph Vince’s Leverage Space Trading Model

In honor of the press release Dow Jones Indexes To Develop, Co-Brand Index Family With LSP Partners two days ago, I thought I would show another slightly different use of Ralph Vince’s The Leverage Space Trading Model.

Using the R LSPM package, we can build a monthly system around the probProfit calculation.  This particular system will enter long when the short term (12 month) probProfit exceeds the longer term (36 month) probProfit.  It exits when the short term falls below the longer term.

From TimelyPortfolio
From TimelyPortfolio

Feel free to substitute any index.  Some of my favorites are German Dax GDAXI, Japan Nikkei N225, Korea Kospi KS11, and Signapore Straits Times STI for international testing.  Additional US testing might look at NDX, RUT, CYC, XBD, HGX, REI, DJUSBK, OSX or anything that you can think of to break it.

The results are not fantastic, but the considerable drawdown reductions is nice.  Let me know how you would improve.

R code:

#Please see au.tra.sy blog http://www.automated-trading-system.com/
#for original code and http://www.fosstrading.com for some of the
#other techniques

require(PerformanceAnalytics)
require(PApages)
require(quantmod)
require(LSPM)

tckr<-"^GSPC"

start<-"1929-01-01"
end<- format(Sys.Date(),"%Y-%m-%d") # yyyy-mm-dd

# Pull tckr index data from Yahoo! Finance
getSymbols(tckr, from=start, to=end)

GSPC<-adjustOHLC(GSPC,use.Adjusted=T)

GSPC<-to.monthly(GSPC)
rtn<-monthlyReturn(GSPC[,4])
# Define JPT function
jointProbTable <- function(x, n=3, FUN=median, ...) {
  # handle case with no negative returns; use -0.01
  for (sys in 1:numsys) {
      if (min(x[,sys])> -1) x[,sys][which.min(x[,sys])]<- -0.01
  }
  # Function to bin data
  quantize <- function(x, n, FUN=median, ...) {
    if(is.character(FUN)) FUN <- get(FUN)
    bins <- cut(x, n, labels=FALSE)
    res <- sapply(1:NROW(x), function(i) FUN(x[bins==bins[i]], ...))
  }
  # Allow for different values of 'n' for each system in 'x'
  if(NROW(n)==1) {
    n <- rep(n,NCOL(x))
  } else
  if(NROW(n)!=NCOL(x)) stop("invalid 'n'")
  # Bin data in 'x'
  qd <- sapply(1:NCOL(x), function(i) quantize(x[,i],n=n[i],FUN=FUN,...))
  # Aggregate probabilities
  probs <- rep(1/NROW(x),NROW(x))
  res <- aggregate(probs, by=lapply(1:NCOL(qd), function(i) qd[,i]), sum)
  # Clean up output, return lsp object
  colnames(res) <- colnames(x)
  res <- lsp(res[,1:NCOL(x)],res[,NCOL(res)])
  return(res)
}

# I know there are prettier ways to accomplish
# but I have to live within my limits

numsys<-1
numbins<-12

# Set Walk-Forward parameters (number of periods) for short
optim<-9 # 9 monthly returns
wf<-1 #walk forward 1 month; we'll set horizon separately
# Calculate number of WF cycles
numCycles = floor((nrow(rtn)-optim)/wf)

for (i in 0:(numCycles-1)) {
            # Define cycle boundaries
            start<-1+(i*wf)
            end<-optim+(i*wf)
            # Get returns for optimization cycle and create the JPT
            jpt <- jointProbTable(rtn[start:end,1:numsys],n=rep(numbins,numsys))
            outcomes<-jpt[[1]]
            probs<-jpt[[2]]
            port<-lsp(outcomes,probs)
            profitProb<-probProfit(port,target=0,horizon=6)
        profitProbWF<-c(rep(1,wf)) %o% profitProb
        maxLossWF<-c(rep(1,wf)) %o% jpt$maxLoss
        #make xts
        profitProbWF<-xts(profitProbWF,order.by=index(rtn[(end+1):(end+wf)]))
        maxLossWF<-xts(maxLossWF,order.by=index(rtn[(end+1):(end+wf)]))

            if (i==0) profitProbHistory<-profitProbWF else profitProbHistory<-rbind(profitProbHistory,profitProbWF)
            if (i==0) maxLossHistory<-maxLossWF else maxLossHistory<-rbind(maxLossHistory,maxLossWF)       
}

# Set Walk-Forward parameters (number of periods) for long
optim<-30 # 30 monthly returns
wf<-1 #walk forward 1 month; we'll set horizon separately
# Calculate number of WF cycles
numCycles = floor((nrow(rtn)-optim)/wf)

for (i in 0:(numCycles-1)) {
            # Define cycle boundaries
            start<-1+(i*wf)
            end<-optim+(i*wf)
            # Get returns for optimization cycle and create the JPT
            jpt <- jointProbTable(rtn[start:end,1:numsys],n=rep(numbins,numsys))
            outcomes<-jpt[[1]]
            probs<-jpt[[2]]
            port<-lsp(outcomes,probs)
            profitProb<-probProfit(port,target=0,horizon=3)
        profitProbWF<-c(rep(1,wf)) %o% profitProb
        maxLossWF<-c(rep(1,wf)) %o% jpt$maxLoss
        #make xts
        profitProbWFlong<-xts(profitProbWF,order.by=index(rtn[(end+1):(end+wf)]))
        maxLossWFlong<-xts(maxLossWF,order.by=index(rtn[(end+1):(end+wf)]))

            if (i==0) profitProbHistorylong<-profitProbWFlong else profitProbHistorylong<-rbind(profitProbHistorylong,profitProbWFlong)
            if (i==0) maxLossHistorylong<-maxLossWFlong else maxLossHistorylong<-rbind(maxLossHistorylong,maxLossWFlong)       
}

signalshortterm<-profitProbHistory
#adjust the long term with maxLoss to hopefully reduce drawdown
signallongterm<-profitProbHistorylong - maxLossHistorylong

chartSeries(signalshortterm,TA="addTA(signallongterm,on=1)", theme="white", name="Short and Long Term Probability of Profit")

# Create the signals and enter when long term is < short term
sigup <- ifelse(signallongterm < signalshortterm,1,0)

# no need for lag since signal generated from previous months]
# sigup <- lag(sigup,1) # Note k=1 implies a move *forward*

# Replace missing signals with no position
# (generally just at beginning of series)
sigup[is.na(sigup)] <- 0

#Calculate equity curves
eq_up <- cumprod(1+(rtn)*sigup)

perf_compare<-merge(sigup*rtn,rtn[(optim+1):NROW(rtn)])
colnames(perf_compare)<-c("LSPM probProfit System",tckr)

charts.PerformanceSummary(perf_compare,ylog=TRUE,legend.loc="topleft",main="LSPM probProfit System Performance Comparison")

Tuesday, April 26, 2011

Great FAJ Article on Statistical Measure of Financial Turbulence Part 3

Building on posts Great FAJ Article on Statistical Measure of Financial Turbulence and Great FAJ Article on Statistical Measure of Financial Turbulence Part 2, I will now build a system incorporating a new correlation-based measure of turbulence and a simple relative strength technique.  First, the new correlation-based measure looks like this

From TimelyPortfolio
From TimelyPortfolio

Then if we build a system with just this new correlation-based turbulence, we get fairly good results.  However, if we take one step extra and apply a slope-based relative strength to choose S&P500 or CRB rather than equal-weight both, we get even better results.

From TimelyPortfolio

I always want more data in countries and years, but 1957 is pretty thorough.  Please let me know how you would improve this system.

 

R code:

require(quantmod)
require(PerformanceAnalytics)

#get data from St. Louis Federal Reserve (FRED)
getSymbols("GS20",src="FRED") #load 20yTreasury; 20y has gap 86-93; 30y has gap in early 2000s
getSymbols("GS30",src="FRED") #load 30yTreasury to fill 20y gap 86-93
getSymbols("BAA",src="FRED") #load BAA
getSymbols("SP500",src="FRED") #load SP500

#get CRB data from a csv file
CRB<-as.xts(read.csv("crb.csv",row.names=1))[,1]

#fill 20y gap from discontinued 20y Treasuries with 30y
GS20["1987-01::1993-09"]<-GS30["1987-01::1993-09"]

#do a little manipulation to get the data lined up on monthly basis
SP500<-to.monthly(SP500)[,4]
#get monthly format to yyyy-mm-dd with the first day of the month
index(SP500)<-as.Date(index(SP500))
#my CRB data is end of month; could change but more fun to do in R
CRB<-to.monthly(CRB)[,4]
index(CRB)<-as.Date(index(CRB))

#let's merge all this into one xts object; CRB starts last in 1956
assets<-na.omit(merge(GS20,SP500,CRB))
#use ROC for SP500 and CRB and momentum for yield data
assetROC<-na.omit(merge(momentum(assets[,1])/100,ROC(assets[,2:3],type="discrete")))

#get Correlations
corrBondsSp<-runCor(assetROC[,1],assetROC[,2],n=7)
corrBondsCrb<-runCor(assetROC[,1],assetROC[,3],n=7)
corrSpCrb<-runCor(assetROC[,2],assetROC[,3],n=7)
#composite measure of correlations between asset classes and roc-weighted correlations
assetCorr<-(corrBondsSp+corrBondsCrb+corrSpCrb+
    (corrBondsSp*corrSpCrb*assetROC[,2])+
    (corrBondsCrb*corrSpCrb*assetROC[,3])-
    assetROC[,1])/6
#sum of ROCs of asset classes
assetROCSum<-assetROC[,1]+assetROC[,2]+assetROC[,3]
#finally the turbulence measure
turbulence<-abs(assetCorr*assetROCSum*100)
colnames(turbulence)<-"Turbulence-correlation"

chartSeries(turbulence,theme="white",name="Correlation and % Change as Measure of Financial Turbulence")
abline(h=0.8)

chart.ACF(turbulence,main="Auto-correlation of Turbulence")

#wish I could remember where I got some of this code
#most likely candidate is www.fosstrading.com
#please let me know if you know the source
#so I can give adequate credit

#use turbulence to determine in or out of equal-weighted sp500 and crb
signal<-ifelse(turbulence>0.8,0,1)
#use slope of sp500/crb to determine sp500 or crb

signal<-lag(signal,k=1)
# Replace missing signals with no position
# (generally just at beginning of series)
signal[is.na(signal)] <- 0

#get returns from equal-weighted crb and sp500 position; Return.portfolio was causing problems, so did the hard way
ret<-(assetROC[,2]+assetROC[,3])/2
ret[1] <- 0

#get system performance
system_perf<-ret*signal
system_eq<- cumprod(1+ret*signal)

perf_comparison<-merge(lag((assetROC[,2]+assetROC[,3])/2,k=1),system_perf)
colnames(perf_comparison)<-c("Equal-weighted","System-with-turbulence-filter")

charts.PerformanceSummary(perf_comparison,ylog=TRUE,main="Turbulence-based System vs Equal-Weighted CRB and SP500")

#let's use basic relative strength to pick sp500 or crb
#know I can do this better in R but here is my ugly code
#to calculate 12 month slope of sp500/crb
width=12
for (i in 1:(NROW(assets)-width)) {
    model<-lm(assets[i:(i+width),2]/assets[i:(i+width),3]~index(assets[i:(i+width)]))
    ifelse(i==1,assetSlope<-model$coefficients[2],assetSlope<-rbind(assetSlope,model$coefficients[2]))
}
assetSlope<-xts(cbind(assetSlope),order.by=index(assets)[(width+1):NROW(assets)])
#use turbulence to determine in or out of equal-weighted sp500 and crb
signal<-ifelse(turbulence[(width):NROW(turbulence)]>0.8,0,1)

#use slope of sp500/crb to determine sp500 or crb
signal2<-ifelse(assetSlope>0,1,2)

signal<-lag(signal,k=1)
signal[1]<-0
signal2<-lag(signal2,k=1)
signal2[1]<-0

#get sp500 or crb return based on slope
ret<-ifelse(signal2==1,assetROC[(width):NROW(assetROC),2],ifelse(signal2==2,assetROC[(width):NROW(assetROC),3],0))

#get system performance
system_perf_rs<-ret*signal
system_eq_rs<- cumprod(1+ret*signal)

perf_comparison<-merge((assetROC[,2]+assetROC[,3])/2,assetROC[,2],assetROC[,3],system_perf,system_perf_rs)
colnames(perf_comparison)<-c("Equal-weighted","S&P500","CRB","System-with-turbulence-filter","System-with-turbulence-filter and RS")

charts.PerformanceSummary(perf_comparison,ylog=TRUE,main="Turbulence-based System with RS vs Equal-Weighted, CRB, and SP500")

Great FAJ Article on Statistical Measure of Financial Turbulence Part 2

faj abstract

I did not intend for this to be a multi-part series, but after some clear thinking at the beach over the weekend, I decided that it needed some more analysis.  For those of you that read the article or know Mahalanobis distance, the measure I presented in my post Great FAJ Article on Statistical Measure of Financial Turbulence was a derivative of actual Mahalanobis distance.  To close the gap between the actual calculation and my measure, I thought I should show both calculations.

From TimelyPortfolio
From TimelyPortfolio

The Mahalanobis calculation here is based on the entire dataset.  Of course, we do not have the luxury of hindsight in trading, so I think my real-time measure works well as displayed by the system shown in the first post.

 

R code:

require(quantmod)
require(PerformanceAnalytics)

#get data from St. Louis Federal Reserve (FRED)
getSymbols("GS20",src="FRED") #load 20yTreasury; 20y has gap 86-93; 30y has gap in early 2000s
getSymbols("GS30",src="FRED") #load 30yTreasury to fill 20y gap 86-93
#getSymbols("BAA",src="FRED") #load BAA
getSymbols("SP500",src="FRED") #load SP500

#get CRB data from a csv file
CRB<-as.xts(read.csv("crb.csv",row.names=1))[,1]

#fill 20y gap from discontinued 20y Treasuries with 30y
GS20["1987-01::1993-09"]<-GS30["1987-01::1993-09"]

#do a little manipulation to get the data lined up on monthly basis
SP500<-to.monthly(SP500)[,4]
#get monthly format to yyyy-mm-dd with the first day of the month
index(SP500)<-as.Date(index(SP500))
#my CRB data is end of month; could change but more fun to do in R
CRB<-to.monthly(CRB)[,4]
index(CRB)<-as.Date(index(CRB))

#let's merge all this into one xts object; CRB starts latest in 1956
assets<-na.omit(merge(GS20,SP500,CRB))
#use ROC for SP500 and CRB and momentum for yield data
assetROC<-na.omit(merge(momentum(assets[,1])/100,ROC(assets[,2:3],type="discrete")))
#get Covariances and multiply to by 100000 for 20y to sp500 and crb and 1000 for sp500 to crb to standardize
#don't like this manual intervention; next post will use correlation instead
assetCovar<-runCov(assetROC[,1],assetROC[,2],n=2)*100000+runCov(assetROC[,1],assetROC[,3],n=2)*100000+runCov(assetROC[,2],assetROC[,3],n=2)*1000
assetROCSum<-assetROC[,1]+assetROC[,2]+assetROC[,3]
turbulence<-abs(assetCovar*assetROCSum)

Sx <- cov(assetROC)
mahalanobisturbulence <- mahalanobis(assetROC, colMeans(assetROC), Sx)

mahalanobisturbulence<-as.xts(cbind(mahalanobisturbulence))

turbulencecompare<-merge(mahalanobisturbulence,turbulence)
colnames(turbulencecompare)<-c("Mahalanobis","MyEstimate")

chart.TimeSeries(turbulencecompare,main="Mahalanobis Distance Compared to My Measure",
    colorset=c("cadetblue","darkolivegreen3"),legend.loc="topleft")
chart.Correlation(turbulencecompare)

Thursday, April 21, 2011

Great FAJ Article on Statistical Measure of Financial Turbulence

faj abstract

I particularly liked this well-written paper, since unlike most academic research, I was able to understand it, replicate it, and incorporate it.  I know that the Financial Analyst Journal is not considered by the academic community as a top-tier journal, but its insights have deepened my understanding and stimulated a lot of very beneficial thoughts.

As usual, I have struggled with how I could write this post with public data, but have not found an adequate public replacement for the CRB Index, so I will show how I used R to get similar results and build a basic system around it without including a file with the CRB monthly closes.  Everything else is publicly available, and I hope that even though the data is not included, the method and result will be very clear.

From TimelyPortfolio

When building the system around it, I was hesitant to do an optimization window with a walkforward as the article suggests, so with hindsight unfortunately, I picked a value (1) that works well.  In the next iteration, I achieve much better results with far less hindsight optimization bias by using correlation instead of covariance.  I’ll save that for my next post after a short weekend trip to the beach.

From TimelyPortfolio

R code:

require(quantmod)
require(PerformanceAnalytics)

#get data from St. Louis Federal Reserve (FRED)
getSymbols("GS20",src="FRED") #load 20yTreasury; 20y has gap 86-93; 30y has gap in early 2000s
getSymbols("GS30",src="FRED") #load 30yTreasury to fill 20y gap 86-93
#getSymbols("BAA",src="FRED") #load BAA
getSymbols("SP500",src="FRED") #load SP500

#get CRB data from a csv file
CRB<-as.xts(read.csv("crb.csv",row.names=1))[,1]

#fill 20y gap from discontinued 20y Treasuries with 30y
GS20["1987-01::1993-09"]<-GS30["1987-01::1993-09"]

#do a little manipulation to get the data lined up on monthly basis
SP500<-to.monthly(SP500)[,4]
#get monthly format to yyyy-mm-dd with the first day of the month
index(SP500)<-as.Date(index(SP500))
#my CRB data is end of month; could change but more fun to do in R
CRB<-to.monthly(CRB)[,4]
index(CRB)<-as.Date(index(CRB))

#let's merge all this into one xts object; CRB starts latest in 1956
assets<-na.omit(merge(GS20,SP500,CRB))
#use ROC for SP500 and CRB and momentum for yield data
assetROC<-na.omit(merge(momentum(assets[,1])/100,ROC(assets[,2:3],type="discrete")))
#get Covariances and multiply to by 100000 for 20y to sp500 and crb and 1000 for sp500 to crb to standardize
#don't like this manual intervention; next post will use correlation instead
assetCovar<-runCov(assetROC[,1],assetROC[,2],n=2)*100000+runCov(assetROC[,1],assetROC[,3],n=2)*100000+runCov(assetROC[,2],assetROC[,3],n=2)*1000
assetROCSum<-assetROC[,1]+assetROC[,2]+assetROC[,3]
turbulence<-abs(assetCovar*assetROCSum)
chartSeries(turbulence,theme="white",name="Covariance and % Change as Measure of Financial Turbulence")

#wish I could remember where I got some of this code
#most likely candidate is www.fosstrading.com
#please let me know if you know the source
#so I can give adequate credit

#use turbulence to determine in or out of equal-weighted sp500 and crb
signal<-ifelse(turbulence>1,0,1)
signal<-lag(signal,k=1)
# Replace missing signals with no position
# (generally just at beginning of series)
signal[is.na(signal)] <- 0

#get returns from equal-weighted crb and sp500 position; Return.portfolio was causing problems, so did the hard way
ret<-(assetROC[,2]+assetROC[,3])/2
ret[1] <- 0

#get system performance
system_perf<-ret*signal
system_eq<- cumprod(1+ret*signal)

perf_comparison<-merge(lag((assetROC[,2]+assetROC[,3])/2,k=1),system_perf)
colnames(perf_comparison)<-c("Equal-weighted","System-with-turbulence-filter")

charts.PerformanceSummary(perf_comparison,ylog=TRUE,main="Turbulence-based System vs Equal-Weighted CRB and SP500")

Wednesday, April 20, 2011

New Favorite Test of US Monetary Policy Limits

After a little additional thought, I discovered that my Death Spiral Warning Graph post can be improved through the isolation of the expected inflation component of US 10y yields provided by the US 10y yield – US 10y TIP yield.  Unfortunately, it more accurately reveals the cost of aggressive US Fed monetary policy, and I think offers a more dire assessment of the limits of this policy.

From TimelyPortfolio

The gains we have experienced seem much more transitory and spurious than a simple chart of the S&P 500 reveals.

R code:

require(quantmod)
require(PerformanceAnalytics)

#get data from St. Louis Federal Reserve (FRED)
getSymbols("DGS10",src="FRED") #load 10yTreasury
getSymbols("DFII10",src="FRED") #load 10yTIP for real return
getSymbols("DTWEXB",src="FRED") #load US dollar
getSymbols("SP500",src="FRED") #load SP500

#new FED monetary policy limit monitor
#US10Y TIP Breakeven Inflation divided by the US Dollar
#old favorite is US 10y Nominal divided by the US Dollar

fedpolicytest<-merge((DGS10-DFII10)/DTWEXB,DGS10/DTWEXB)
colnames(fedpolicytest)<-c("US10yTIPInflationDivByUSDollar (new favorite)","US10yDivByDollar (old favorite)")
chart.TimeSeries(fedpolicytest["2003::2011"],legend.loc="topleft",
    main="US Monetary Policy Test-Does Bernanke Have a Limit?",colorset=c("cadetblue","darkolivegreen3"),ylab="",
    event.lines="2008-07-02",event.labels="July 2, 2008")
#add label to denote new recent high
text(NROW(fedpolicytest["2003::2011"])-15,fedpolicytest[NROW(fedpolicytest),1]+.001,"New High")
#label previous high
abline(h=fedpolicytest["2008-07-02",1],col="cadetblue")

Tuesday, April 19, 2011

Barron’s Spring 2008 Big Money Poll

Barron's April 28, 2008, Cover Story "Back in the Pool" offers a great hindsight look at our wonderful foresight:

“AND NOW, FOR SOME GOOD NEWS: THE OTHER SHOE isn't going to drop. After a winter of discontent marked by massive write-offs on Wall Street and a wilting economy on Main, America's portfolio managers have declared that the worst is over. More than half of the institutional investors participating in our latest Big Money poll say they're bullish or very bullish about the prospects for stocks through the end of 2008. Their forecasts suggest they're even more upbeat about the first half of 2009.”

From TimelyPortfolio

Even the world’s greatest investors had been trained by the market from 1980-2000 to buy every 20% dip and ignore fears.  However, as usual, the market trains you over rolling 5-20 year periods to make the wrong decision.  Closing your eyes and buying in 2008 was clearly the wrong decision, and I’m uncertain if closing your eyes and buying in 2009 was the right decision.  Closing your eyes and buying/holding anything now without regard to how we got here is very precarious.  It will be very interesting to see the results of this spring’s Barrons Big Money Poll.

As always, I use this blog as my own record to avoid the errors of hindsight bias and cognitive dissonance.  In no way am I advising you to believe me since I could be as wrong as the money managers in Barrons April 2008.  You are responsible for your own investment gains and losses.

R code:

require(quantmod)

tckr<-"^GSPC"

start<-"1990-01-01"
end<- format(Sys.Date(),"%Y-%m-%d") # yyyy-mm-dd

# Pull tckr index data from Yahoo! Finance
getSymbols(tckr, from=start, to=end, adjust=TRUE)

chartSeries(to.weekly(GSPC["2007::2009"]),theme="white",name="S&P 500-Safe to Get Back in the Pool?",TA=NULL)
#has to be a better way to label; let me know if you know how
text(230, GSPC["2008-05-02",4]+25, labels="Safe and Bullish???")

Sunday, April 17, 2011

Historical Sources of Bond Returns-Comparison of Daily to Monthly

Thanks so much for the comment on my last post Historical Bond Price and Total Returns from 10y Yield Series

“I know this might sound antithetical to a bond guy, but won't the monthly series get you close enough? “

which proved me wrong and allowed me to get nine additional years of history.  The question caused me to go one step further in proving my statement “…but monthly yields from Shiller and Fed are averages of the daily yields and not month-end values.  Applying this method to these month averages does not work.”  This monthly averaging has caused many problems in the past for me, but in this case, it does not matter.  We can see the difference in monthly averaging (GS10) and month-end daily values (DGS10) in this chart.

From TimelyPortfolio

However, if we look at the price and total return series from both, the differences are insignificant.

From TimelyPortfolio

Now that I have cleared that up, I can continue Bond Market as a Casino Game with data back to 1953.

R code:

#do everything twice to compare monthly average to daily

require(RQuantLib)
require(quantmod)
require(PerformanceAnalytics)

getSymbols("DGS10",src="FRED") #load daily US Treasury 10y from Fed Fred
getSymbols("GS10",src="FRED") # load monthly average US 10y from Fed Fred

#Fed monthly series of yields is the monthly average of daily yields
#Shiller series also is monthly average
#to calculate returns monthly average does not work
#so we get daily, take monthend with to.monthly
#and unfortunately the series does not go back as far
DGS10<-to.monthly(DGS10)[,4]
#set index to yyyy-mm-dd format rather than to.monthly mmm yyy for better merging later
index(DGS10)<-as.Date(index(DGS10))
DGS10pricereturn<-DGS10  #set this up to hold price returns
GS10pricereturn<-GS10

DGS10pricereturn[1,1]<-0
colnames(DGS10pricereturn)<-"PriceReturn-daily to monthly DGS10"
#I know I need to vectorize this but not qualified enough yet
#Please feel free to comment to show me how to do this
for (i in 1:(NROW(DGS10)-1)) {
  DGS10pricereturn[i+1,1]<-FixedRateBondPriceByYield(yield=DGS10[i+1,1]/100,issueDate=Sys.Date(),
    maturityDate= advance("UnitedStates/GovernmentBond", Sys.Date(), 10, 3),
    rates=DGS10[i,1]/100,period=2)[1]/100-1
}

GS10pricereturn[1,1]<-0
colnames(GS10pricereturn)<-"PriceReturn-monthly avg GS10"
#I know I need to vectorize this but not qualified enough yet
#Please feel free to comment to show me how to do this
for (i in 1:(NROW(GS10)-1)) {
  GS10pricereturn[i+1,1]<-FixedRateBondPriceByYield(yield=GS10[i+1,1]/100,issueDate=Sys.Date(),
    maturityDate= advance("UnitedStates/GovernmentBond", Sys.Date(), 10, 3),
    rates=GS10[i,1]/100,period=2)[1]/100-1
}

#total return will be the price return + yield/12 for one month
DGS10totalreturn<-DGS10pricereturn+lag(DGS10,k=1)/12/100
colnames(DGS10totalreturn)<-"Total Return-daily to monthly DGS10"
#total return will be the price return + yield/12 for one month
GS10totalreturn<-GS10pricereturn+lag(GS10,k=1)/12/100
colnames(GS10totalreturn)<-"TotalReturn-monthly avg GS10"

chartSeries(DGS10,TA="addTA(GS10,on=1);addTA((DGS10-GS10)/100)",name="Comparison of DGS10 and GS10",theme="white")
mtext("Source: Federal Reserve FRED",side=1,adj=0)

charts.PerformanceSummary(merge(DGS10pricereturn,DGS10totalreturn,GS10pricereturn,GS10totalreturn),ylog=TRUE,colorset=c("cadetblue","darkolivegreen3","goldenrod","gray70"),main="Simulated Returns from US10y Yield")
mtext("Source: Federal Reserve FRED",side=1,adj=0)

Friday, April 15, 2011

Historical Bond Price and Total Returns from 10y Yield Series

Without access to Barclays or Merrill Bond Indicies to the 1970s or Ned Davis to 1950, studying historical bond returns is very difficult.  Here is a way to derive price and total returns on the 10 year US Treasury back to 1962.  I would like to extend to 1871, but monthly yields from Shiller and Fed are averages of the daily yields and not month-end values.  Applying this method to these month averages does not work.

From TimelyPortfolio

R code:

require(RQuantLib)
require(quantmod)
require(PerformanceAnalytics)

getSymbols("DGS10",src="FRED") #load US Treasury 10y from Fed Fred

#Fed monthly series of yields is the monthly average of daily yields
#Shiller series also is monthly average
#to calculate returns monthly average does not work
#so we get daily, take monthend with to.monthly
#and unfortunately the series does not go back as far
DGS10<-to.monthly(DGS10)[,4]
DGS10pricereturn<-DGS10  #set this up to hold price returns

DGS10pricereturn[1,1]<-0
colnames(DGS10pricereturn)<-"PriceReturn"
#I know I need to vectorize this but not qualified enough yet
#Please feel free to comment to show me how to do this
for (i in 1:(NROW(DGS10)-1)) {
  DGS10pricereturn[i+1,1]<-FixedRateBondPriceByYield(yield=DGS10[i+1,1]/100,issueDate=Sys.Date(),
    maturityDate= advance("UnitedStates/GovernmentBond", Sys.Date(), 10, 3),
    rates=DGS10[i,1]/100,period=2)[1]/100-1
}

#total return will be the price return + yield/12 for one month
DGS10totalreturn<-DGS10pricereturn+lag(DGS10,k=1)/12/100
colnames(DGS10totalreturn)<-"Total Return"

charts.PerformanceSummary(merge(DGS10pricereturn,DGS10totalreturn),ylog=TRUE,colorset=c("cadetblue","darkolivegreen3"),main="Simulated Returns from US10y Yield")
mtext("Source: Federal Reserve FRED",side=1,adj=0)

Monday, April 11, 2011

Historical Sources of Bond Returns with Shiller Data 1919-2011

And as usual, I always want a longer data set, so after a little playing with R-Excel, we can extend our historical sources of bond returns to 1919.  If nothing else, maybe you can find other uses for the Shiller Dataset in R.

From TimelyPortfolio

And as an added bonus, we can prepare a boxplot.

From TimelyPortfolio

As always, please criticize or make suggestions.

 

R code:

require(quantmod)
require(PerformanceAnalytics)

getSymbols("BAA",src="FRED") #load Corporate for credit from Fed Fred

#get Shiller data for inflation and US Treasury 10 year
require(gdata)
URL <- "http://www.econ.yale.edu/~shiller/data/ie_data.xls"
shillerData <- read.xls(URL,sheet="Data",pattern="Rate GS10")
#strip out date(1), CPI(5), and GS10(7)
shillerData <- shillerData[,c(1,5,7)]
#get data starting in January 1919 to match with BAA xts series
shillerData <- shillerData[rownames(shillerData[which(shillerData$Date==format(index(BAA)[1],"%Y.%m")):NROW(shillerData),]),2:3]
shillerData <- xts(shillerData[1:NROW(BAA),1:2],order.by=index(BAA))

bondReturnSources<-merge(BAA,shillerData)
bondReturnSources[,1]<-bondReturnSources[,1]-bondReturnSources[,3]
#get 12 month CPI change
bondReturnSources[,2]<-ROC(bondReturnSources[,2],12,type="discrete")*100
bondReturnSources[,3]<-bondReturnSources[,3]-bondReturnSources[,2]
bondReturnSources<-merge(bondReturnSources,bondReturnSources[,1]+bondReturnSources[,2]+bondReturnSources[,3])

colnames(bondReturnSources)<-c("Credit","Inflation","Real","Total")

chart.TimeSeries(bondReturnSources,legend.loc="bottom",main="Historical Sources of Bond Returns",ylab="Yield as %",colorset=c("cadetblue","darkolivegreen3","goldenrod","gray70"))
#add source as caption
mtext("Source: Federal Reserve FRED and Robert Shiller",side=1,adj=0)

chart.Boxplot(bondReturnSources,main="Historical Sources of Bond Returns",colorset=c("cadetblue","darkolivegreen3","goldenrod","gray70"))
mtext("Source: Federal Reserve FRED and Robert Shiller",side=1,adj=0)

Historical Sources of Bond Returns

As promised in Monitoring Sources of Bond Return, we can show more history if we use CPI instead of expected inflation (from the TIP inflation breakeven yield).  Here are the results with history back to 1953.

From TimelyPortfolio

However, more history includes negative inflation, so I think the chart is a little more clear with a PerformanceAnalytics chart.TimeSeries graph.

From TimelyPortfolio

 

R code:

require(quantmod)
require(PerformanceAnalytics)
require(reshape2)
require(ggplot2)

getSymbols("GS10",src="FRED") #load 10yTreasury
getSymbols("BAA",src="FRED") #load Corporate for credit
getSymbols("CPIAUCSL",src="FRED") #load CPI for inflation

bondReturnSources<-na.omit(merge(ROC(CPIAUCSL,12,type="discrete")*100,BAA-GS10,GS10-ROC(CPIAUCSL,12,type="discrete")*100))
bondReturnSources<-merge(bondReturnSources,bondReturnSources[,1]+bondReturnSources[,2]+bondReturnSources[,3]) #add for total
colnames(bondReturnSources)<-c("Inflation","Credit","Real","Total")

chart.TimeSeries(bondReturnSources,legend.loc="bottom",main="Historical Sources of Bond Returns",ylab="Yield as %",colorset=c("darkolivegreen3","cadetblue","goldenrod","gray70"))
#use this for stacked bar, but comes out strange with too many months
#chart.StackedBar(bondReturnSources["1953::"],main="Historical Sources of Bond Returns",ylab="Yield as %",colorset=c("darkgreen","red","blue"))

#use ggplot for consistency with previous plot, but negative makes chart hard to read
bondReturnSourcesToGraph<-data.frame(cbind(as.Date(index(bondReturnSources)),coredata(bondReturnSources[,1:3])))
colnames(bondReturnSourcesToGraph)<-c("Date","Inflation","Credit","Real")
bondReturnSourcesToGraph<-melt(bondReturnSourcesToGraph,id.vars=1)
colnames(bondReturnSourcesToGraph)<-c("Date","ReturnSource","Yield")
rownames(bondReturnSourcesToGraph)<-c(1:NROW(bondReturnSourcesToGraph))
ggplot(bondReturnSourcesToGraph, stat="identity", aes(x=Date,y=Yield,fill=ReturnSource,group=ReturnSource)) + geom_area() +scale_x_date(format = "%Y") + opts(title = "Historical Sources of Bond Return")

Friday, April 8, 2011

Monitoring Sources of Bond Return

Here is a way to monitor bond return sources in R.  In the next iteration, I will use CPI to add history to the series.

From TimelyPortfolio

So right now, you can expect about a 5% return from bonds.  How much of that is real return is unknown.

R code:

require(quantmod)
require(PerformanceAnalytics)
require(reshape2)
require(ggplot2)

getSymbols("WGS10YR",src="FRED") #load 10yTreasury
getSymbols("WFII10",src="FRED") #load 10yTIP for real return
getSymbols("BAMLC0A0CM",src="FRED") #load Corporate for credit
getSymbols("CPIAUCSL",src="FRED") #load Corporate for credit

bondReturnSources<-na.omit(merge(WGS10YR-WFII10,WFII10,to.weekly(BAMLC0A0CM)[,4])["2003::"])
colnames(bondReturnSources)<-c("10yTreasury","10yTIPReal","Credit")
bondReturnSourcesToGraph<-data.frame(cbind(as.Date(index(bondReturnSources)),coredata(bondReturnSources)))
colnames(bondReturnSourcesToGraph)<-c("Date","InflationBreakEven","10yTIPReal","Credit")
bondReturnSourcesToGraph<-melt(bondReturnSourcesToGraph,id.vars=1)
colnames(bondReturnSourcesToGraph)<-c("Date","ReturnSource","Yield")
rownames(bondReturnSourcesToGraph)<-c(1:NROW(bondReturnSourcesToGraph))
ggplot(bondReturnSourcesToGraph, stat="identity", aes(x=Date,y=Yield,fill=ReturnSource,group=ReturnSource)) + geom_area() +scale_x_date(format = "%Y") + opts(title = "Sources of Bond Return")

Guaranteed Failure with Bonds

If we define failure as

the inability to meet expectations

then bonds are guaranteed failures.

If we define failure in money management as

losing money

then bonds are likely failures.

Do Not Ignore These Charts

See previous posts

Death Spiral Warning Graph
Death Spiral of a Country
Nine Lives of the Fed Put
Unsustainable Gift

for more discussion.

From TimelyPortfolio

 


via StockCharts.com


via StockCharts.com

Thursday, April 7, 2011

Money Management Success

In the spirit of thisisindexed

From TimelyPortfolio

Bond Market as a Casino Game Part 3

For additional insight into the recent 30 year bull bond market, please see the well-researched historical references in Credit Suisse Global Investment Returns Yearbook 2011.  My favorite chart from this is

 

For those even more ambitious, see this massive 736 page beast on interest rates with details back to 3,000 B.C. Sumerian documents.

Wednesday, April 6, 2011

Money Management Skill vs Gains

In a simple chart, one of my observations on money management

From TimelyPortfolio

Bond Market as a Casino Game Part 2

Before starting Part 2, please see Bonds as a Casino Game Part 1.   For the Monte Carlo random simulation purists, please ignore some unimportant technicalities in my simulation.

To spoil the fun, here is the conclusion:

Any way you look at it, the US Bond Market has been a wonderfully profitable game, but the game is changing, and past performance in no way predicts similar future returns.  Please lower your expectations for your bond investments.

Given the probabilities and average outcomes that we have experienced in the Barclays Capital U.S. Aggregate Bond Index since 1980, let’s see what happens if we turn the bond market into a game, and play it in 10,000 30 year games with Microsoft Excel (sorry R, we’ll use you later).

From TimelyPortfolio
From TimelyPortfolio

The results demonstrate how attractive this game is to play with the worst case of 10,000 trials still offering a 5.6% annualized return.  The chart below shows 200 of the 10,000 trials (Excel limits to 255), the actual Barclays Aggregate results, and the average of 10,000 trials.

From TimelyPortfolio

Most interesting to me is how the game seems to be changing to less favorable as the actual approaches the average.  This is also clearly visible in the pivot table by decade as expected returns diminish with lower interest rates.

From TimelyPortfolio

Although using just two buckets is helpful and eases the analysis, we can use more buckets/bins with Ralph Vince’s Leverage Space Trading Model and the R LSPM package from Foss Trading.

From TimelyPortfolio
From TimelyPortfolio

Any way you look at it, the US Bond Market has been a wonderfully profitable game, but the game is changing, and past performance in no way predicts similar future returns.  Please lower your expectations for your bond investments.

R code:

#Please see au.tra.sy blog http://www.automated-trading-system.com/
#for original code and http://www.fosstrading.com
#I take no credit for the majority of this code
#I simply changed a couple of things to use xts return series
#and to do simulations of bonds

require(xts)
require(PerformanceAnalytics)

numbins<-10
# get data series from csv file and limit to 1980 to 2010
rtn<-as.xts(read.csv("lbustruu.csv",row.names=1))["1980::"]
# Calculate number of WF cycles
numCycles = floor((nrow(rtn)-optim)/wf)
# Define JPT function
jointProbTable <- function(x, n=3, FUN=median, ...) {
  # Load LSPM
  if(!require(LSPM,quietly=TRUE)) stop(warnings())

  # Function to bin data
  quantize <- function(x, n, FUN=median, ...) {
    if(is.character(FUN)) FUN <- get(FUN)
    bins <- cut(x, n, labels=FALSE)
    res <- sapply(1:NROW(x), function(i) FUN(x[bins==bins[i]], ...))
  }
  # Allow for different values of 'n' for each system in 'x'
  if(NROW(n)==1) {
    n <- rep(n,NCOL(x))
  } else
  if(NROW(n)!=NCOL(x)) stop("invalid 'n'")
  # Bin data in 'x'
  qd <- sapply(1:NCOL(x), function(i) quantize(x[,i],n=n[i],FUN=FUN,...))
  # Aggregate probabilities
  probs <- rep(1/NROW(x),NROW(x))
  res <- aggregate(probs, by=lapply(1:NCOL(qd), function(i) qd[,i]), sum)
  # Clean up output, return lsp object
  colnames(res) <- colnames(x)
  res <- lsp(res[,1:NCOL(x)],res[,NCOL(res)])
  return(res)
}
# Get returns and create the JPT
jpt <- jointProbTable(rtn,n=numbins)
outcomes<-jpt[[1]]
probs<-jpt[[2]]
port<-lsp(outcomes,probs)

# Analyze drawdowns
# Use LSPM to get Probability of Drawdown > 5%
drawdownProb<-probDrawdown(port,DD=.05,horizon=12)
# Display probability of 5% drawdown over 1 year
actualDraw<-sum(table.Drawdowns(rtn)$Length)/NROW(rtn)
#plot bar chart of expected vs actual
barplot(rbind(drawdownProb,actualDraw),beside=TRUE,main="Expected vs. Actual Drawdown > 5%",names.arg=c("Expected", "Actual"))

# Analyze gains
# Use LSPM to get Probability of Profit > 5%
profitProb<-probProfit(port,target=.05,horizon=12)
# Get rolling 12 month period returns
yearReturns<-apply.rolling(rtn,width=12,by=1,FUN="Return.annualized")
# Get actual percentage of 12 month returns > 5%
actualProb<-NROW(yearReturns[which(yearReturns$calcs>.05),1])/NROW(yearReturns)
#plot bar chart of expected vs actual
barplot(rbind(profitProb,actualProb),beside=TRUE,main="Expected vs. Actual 12 month > 5%",names.arg=c("Expected", "Actual"))

Friday, April 1, 2011

Bond Market as a Casino Game Part 1

With this post, I am doing something I try very hard to avoid, especially when communicating to my clients, and that is blurring the line between investing and gambling.  But after reading all of Reuven Brenner’s books and finishing Ralph Vince Leverage Space Trading Model, I think blurring the line can offer additional insights and methods not traditionally available.

As a starting point, let’s apply basic probability methods on the most widely used bond index Barclays Aggregate Total Return to test my belief that the 1980-current bull run in bonds has offered one of the best games ever for any investor in the history of the markets (see my post “Bonds Tumble and Questions Start Getting Asked”).  In the Microsoft Excel pivot table shown below, the Barclays Aggregate has been up 70% of all month since 1980, and the up months generate on average 1.44% return compared to the down months –1.02%.  If we put this in casino terms, it is like winning 7 of every 10 games with a payout of 1.4 to 1.

 

For another look, let’s use R and PerformanceAnalytics to see the histogram and boxplot.

From TimelyPortfolio

source:Barclays Capital and thanks to the fine contributors to R and PerformanceAnalytics

Now that we have the basic probabilities secured, we can have all sorts of fun in later posts with the Leverage Space techniques and Monte Carlo simulations.

 

R code:

require(PerformanceAnalytics)
require(quantstrat)

#load index data given as date and total return value
BarAgg<-as.xts(read.csv("lbustruu.csv",row.names=1))

#convert to monthly return series for PerformanceAnalytics
#use discrete ROC to get simple monthly return
#((value this month-value last month)/value last month)-1
BarAggReturn<-ROC(BarAgg,n=1,"discrete")

par(mfrow=c(2,1)) #2 rows and 1 column
chart.Histogram(BarAggReturn["1980::"],main="Barclays Aggregate Total Return Monthly % Histogram since 1980",cex.main=1,xlab=NULL,breaks=15,methods="add.centered")
chart.Boxplot(BarAggReturn["1980::"],main="Barclays Aggregate Total Return Monthly % Boxplot since 1980",cex.main=1,xlab=NULL,names=FALSE)