Continuing with the list of tasks defined in the previous post, let’s load the edhec data set,fit each of the nine models to each of the funds and store results in a nested list object :

[4] Load data and store results

For each model and fund combination, I store the following information :

  1. The Log-Likelihood value
  2. The Degree of freedom
  3. The AIC criterion
  4. The SBC criterion

The code follows :

########################################################################################
# Load dataset
########################################################################################

orig.csv <- read.csv('edhecfinal.csv',header=T,stringsAsFactors=F)
edhec.xts <- xts(order.by=seq(as.Date('1997-01-31'),as.Date('2014-07-31'),by='1 month'),orig.csv[,-1])
edhec.xts <- log(edhec.xts+1)

########################################################################################
# Nested Lists and Regression results
########################################################################################

n.funds <- ncol(edhec.xts)
fund.names <- names(edhec.xts)
n.obs <- nrow(edhec.xts)

n.models <- 9

edhec.ret <- list()
	for(i in 1:n.funds){
		edhec.ret[[i]] <- list()
		names(edhec.ret)[i] <- fund.names[i]
		edhec.ret[[i]]$Ret.xts <- edhec.xts[,i]
    logret <- Return.clean(edhec.ret[[i]]$Ret.xts,method='boudt')

			#Independent Log-normal model
			edhec.ret[[i]]$ILN <- lm(logret~1)
		  edhec.ret[[i]]$ILN$Stat <- StatComp(edhec.ret[[i]]$ILN,n.obs)

		  edhec.ret[[i]]$ILN$LogLik <- edhec.ret[[i]]$ILN$Stat['model.loglik']
		  edhec.ret[[i]]$ILN$df <- edhec.ret[[i]]$ILN$Stat['model.df']
		  edhec.ret[[i]]$ILN$AIC <- edhec.ret[[i]]$ILN$Stat['model.aic']
		  edhec.ret[[i]]$ILN$SBC <- edhec.ret[[i]]$ILN$Stat['model.sbc']

		  #Autoregressive model
		  edhec.ret[[i]]$AR1 <- arima(x=logret,order=c(1,0,0))
		  edhec.ret[[i]]$AR1$Stat <- StatComp(edhec.ret[[i]]$AR1,n.obs)

		  edhec.ret[[i]]$AR1$LogLik <- edhec.ret[[i]]$AR1$Stat['model.loglik']
		  edhec.ret[[i]]$AR1$df <- edhec.ret[[i]]$AR1$Stat['model.df']
		  edhec.ret[[i]]$AR1$AIC <- edhec.ret[[i]]$AR1$Stat['model.aic']
		  edhec.ret[[i]]$AR1$SBC <- edhec.ret[[i]]$AR1$Stat['model.sbc']

		  #ARCH model
			arch.spec <- ugarchspec(variance.model=list(garchOrder=c(1,0)),mean.model = list(armaOrder=c(0,0)))
		  arch.fit<- ugarchfit(spec=arch.spec, data=logret,solver.control=list(trace = 1))

		  edhec.ret[[i]]$ARCH$LogLik <- likelihood(arch.fit)
		  edhec.ret[[i]]$ARCH$df <- 3
		  edhec.ret[[i]]$ARCH$AIC <- edhec.ret[[i]]$ARCH$LogLik-edhec.ret[[i]]$ARCH$df
		  edhec.ret[[i]]$ARCH$SBC <- edhec.ret[[i]]$ARCH$LogLik-0.5*edhec.ret[[i]]$ARCH$df*log(n.obs)

	    #AR-ARCH model
			arch.spec <- ugarchspec(variance.model=list(garchOrder=c(1,0)),mean.model = list(armaOrder=c(1,0)))
		  arch.fit<- ugarchfit(spec=arch.spec, data=logret,solver.control=list(trace = 1))

		  edhec.ret[[i]]$AR.ARCH$LogLik <- likelihood(arch.fit)
		  edhec.ret[[i]]$AR.ARCH$df <- 4
		  edhec.ret[[i]]$AR.ARCH$AIC <- edhec.ret[[i]]$AR.ARCH$LogLik-edhec.ret[[i]]$AR.ARCH$df
		  edhec.ret[[i]]$AR.ARCH$SBC <- edhec.ret[[i]]$AR.ARCH$LogLik-0.5*edhec.ret[[i]]$AR.ARCH$df*log(n.obs)

	    #GARCH model
			garch.spec <- ugarchspec(variance.model=list(garchOrder=c(1,1)),mean.model = list(armaOrder=c(0,0)))
		  garch.fit<- ugarchfit(spec=garch.spec, data=logret,solver.control=list(trace = 1))

		  edhec.ret[[i]]$GARCH$LogLik <- likelihood(garch.fit)
		  edhec.ret[[i]]$GARCH$df <- 4
		  edhec.ret[[i]]$GARCH$AIC <- edhec.ret[[i]]$GARCH$LogLik-edhec.ret[[i]]$GARCH$df
		  edhec.ret[[i]]$GARCH$SBC <- edhec.ret[[i]]$GARCH$LogLik-0.5*edhec.ret[[i]]$GARCH$df*log(n.obs)

		  #AR-GARCH model
			garch.spec <- ugarchspec(variance.model=list(garchOrder=c(1,1)),mean.model = list(armaOrder=c(1,0)))
		  garch.fit<- ugarchfit(spec=garch.spec, data=logret,solver.control=list(trace = 1))

		  edhec.ret[[i]]$AR.GARCH$LogLik <- likelihood(garch.fit)
		  edhec.ret[[i]]$AR.GARCH$df <- 5
		  edhec.ret[[i]]$AR.GARCH$AIC <- edhec.ret[[i]]$AR.GARCH$LogLik-edhec.ret[[i]]$AR.GARCH$df
		  edhec.ret[[i]]$AR.GARCH$SBC <- edhec.ret[[i]]$AR.GARCH$LogLik-0.5*edhec.ret[[i]]$AR.GARCH$df*log(n.obs)

		  #ILN 2 regimes
		  model.spec <- depmix(eval(parse(text=fund.names[i]))~1,nstates=2,data=logret)
	    model.fit <- fit(model.spec)
		  edhec.ret[[i]]$ILN.Reg2$Stat <- StatComp(model.fit,n.obs)

      edhec.ret[[i]]$ILN.Reg2$LogLik <- edhec.ret[[i]]$ILN.Reg2$Stat['model.loglik']
		  edhec.ret[[i]]$ILN.Reg2$df <- edhec.ret[[i]]$ILN.Reg2$Stat['model.df']
		  edhec.ret[[i]]$ILN.Reg2$AIC <- edhec.ret[[i]]$ILN.Reg2$Stat['model.aic']
		  edhec.ret[[i]]$ILN.Reg2$SBC <- edhec.ret[[i]]$ILN.Reg2$Stat['model.sbc']

	    #ILN 3 regimes
		  model.spec <- depmix(eval(parse(text=fund.names[i]))~1,nstates=3,data=logret)
	    model.fit <- fit(model.spec)
	  	edhec.ret[[i]]$ILN.Reg3$Stat <- StatComp(model.fit,n.obs)

      edhec.ret[[i]]$ILN.Reg3$LogLik <- edhec.ret[[i]]$ILN.Reg3$Stat['model.loglik']
		  edhec.ret[[i]]$ILN.Reg3$df <- edhec.ret[[i]]$ILN.Reg3$Stat['model.df']
		  edhec.ret[[i]]$ILN.Reg3$AIC <- edhec.ret[[i]]$ILN.Reg3$Stat['model.aic']
		  edhec.ret[[i]]$ILN.Reg3$SBC <- edhec.ret[[i]]$ILN.Reg3$Stat['model.sbc']

      #AR-Regime switch model
		  temp.df <- data.frame(logret[2:n.obs],logret[2:n.obs-1])
		  names(temp.df) <- c('left','right')

		  model.spec <- depmix(left~right,nstates=2,data=temp.df)
	    model.fit <- fit(model.spec)

			edhec.ret[[i]]$AR.Reg2$Stat <- StatComp(model.fit,n.obs)

      edhec.ret[[i]]$AR.Reg2$LogLik <- edhec.ret[[i]]$AR.Reg2$Stat['model.loglik']
		  edhec.ret[[i]]$AR.Reg2$df <- edhec.ret[[i]]$AR.Reg2$Stat['model.df']
		  edhec.ret[[i]]$AR.Reg2$AIC <- edhec.ret[[i]]$AR.Reg2$Stat['model.aic']
		  edhec.ret[[i]]$AR.Reg2$SBC <- edhec.ret[[i]]$AR.Reg2$Stat['model.sbc']

	}

l
To give an idea of the output,consider the case of CTA.Global fund and the Autoregressive-Regime switching model :

ff

 

The StatComp() function used to calculate these values is simply :

########################################################################################
# Calculating Stats for Comparison across models
########################################################################################

StatComp <- function(fitted.model,n.obs){
	    model.df <- attr(logLik(fitted.model),which='df')
			model.loglik <- as.numeric(logLik(fitted.model))
		  model.aic <- model.loglik-model.df
		  model.sbc <- model.loglik-0.5*model.df*log(n.obs)

	return(c(model.df=model.df,model.loglik=model.loglik,model.aic=model.aic,model.sbc=model.sbc))
}

ll
Once we have fitted each of the nine models to the data for each of the 13 funds in the edhec dataset,let’s summarise this information in a new list object by saving across models.

 

[5] Saving across models and making data ggplot friendly

The purpose of this step in the process is to summarise the data contained in the edhec.ret list object above into something that can be conveniently used in ggplot functions. Problematic funds are subsequently identified and deleted.

# Saving across models & removing problematic funds
model.names <- names(edhec.ret[[1]])[-1]
model.list <- list()
problems <- NULL

for(i in 1:n.models){
	model.list[[i]] <- list()
	names(model.list)[i] <- model.names[i]
	model.list[[i]] <- do.call(rbind,lapply(edhec.ret,function(x) eval(parse(text=paste('x$',model.names[i],sep='')))))[,c('AIC','SBC','df','LogLik')]
  problems <- c(problems,which(is.na(model.list[[i]]==0) | model.list[[i]]<0))
  model.list[[i]] <- cbind(model.list[[i]],Model=rep(model.names[i],n.funds))
}

problem.funds <- unique(problems)[-which(problems>13)][-5]
model.list<-lapply(model.list,function(x) x[-problem.funds,])
subfund.names <- rownames(model.list[[1]])

temp <- apply(do.call(rbind,model.list),2,unlist)

ll
By problematic funds, i mean those for which the following error messages occurred in the fitting process:

gff

 

For me the problematic funds are the ones corresponding to the following indices :

kli

 

Ultimately,our model.list object looks like this :

jjuk

 

This is the first element of the list. It continues this way for all other models listed in the previous post. The last line of code creates a matrix called temp which contains the unlisted version of the data contained in model.list. In the next post we will use this temp object to create three data frames, each of which we will manipulate using the tidyR package to produce the desired data structure for ggplotting purposes.

Advertisements

In my attempt to replicate some of the methodologies of the Kritzman paper, I stumbled across an article by Hardy (2001) which I thought I might try to tinker with as well. To summarise :  the paper itself is concerned with the following issues [1] A comparison of the fit of Regime Switching Log Normal (RSLN) model with other conventional models in common use for the SP500 / TSE 300 indices ; [2] Derivation of the distribution function for the RSLN model ; [3] Derivation of the European option price using the RSLN model ; [4] Application of the model to calculate risk measures.

As usual I only replicate the most basic issues,in this case a comparison of the fit across a host of models on the basis of the edhec data series from previous blog entries (as opposed to the broad market indices used in the paper).

To make this write up easier on me, this and subsequent posts shall deal with the following tasks :

  1. Provide a summary of the models examined by the paper.
  2. Provide a summary of the selection criteria (used to rank fitted results across models within and across fund returns).
  3. Attempt to make a custom ggplot table maker function to simplify text plotting.
  4. Load the edhec data set and store model results in a nested list object such that :
    • We fit 9 models to each of the 13 funds in the data set.
    • We store the estimated parameters and calculate values for chosen selection criteria for each model & fund combination.
    • We can access ILN model in fund CTA.Global by : Nested_List$CTA.Global$ILN.
  5. Save across models & remove problematic funds (funds that had strange values in the fitting process)  such that :
    •  We have a list object of 9 elements :
      • with each element being named after one of the 9 models.
      • with each element containing a data frame of non problematic funds and their selection criteria, formatted in a ggplot friendly way.
  6. Make a facet plot of Selection Criteria across non-problematic funds of the edhec dataset.
  7. Make a snap plot function for the chosen fund and selection criteria such that :
    • We provide a ranking of models for the particular fund and selection criteria combination in tabular as well as graphic form.
  8. Make a plot showing how often each model is the best performer across non-problematic funds.

 

[1] Summary of the models used

(1) The Independent Log Normal Model (ILN)

s1

(2) The First-order autoregressive model (AR1)

s2

(3) The autoregressive conditionally heteroskedastic model (ARCH1)

s3

(4) Combination of AR and ARCH model (AR.ARCH)

s4

(5) The generalized autoregressive conditionally heteroskedastic model (GARCH)

s5

(6) Combination of AR and GARCH model (AR.GARCH)

s7

(7) Regime Switching AR(1) model (AR.Reg2)

s8

(8) Regime Switching ILN model with 2 regimes (ILN.Reg2)

(9) Regime Switching ILN model with 3 regimes (ILN.Reg3)

 

[2] Summary of the Selection Criteria used

While the paper uses the (1) Likelihood ratio test, (2) Akaike Information Criterion and (3) Schwartz Bayes Criterion, I only used the the last two. The AIC uses the model that maximises the difference between the likelihood and the degree of freedom. The SBC uses the model that maximises the following :

lm

where I is the likelihood ; k the degree of freedom ; n the sample size ; j represents the model.

 

[3] Custom ggplot Table Maker

Since it is often useful to be able to plot tabular data on the fly,I have written a simple function for this purpose which allows one to specify typical elements of a basic table :

  • Is there a column and/or row?
  • What is the column title and/or row title?
  • What is the font for these titles?
  • What is the colour of the background of those titles?
  • What is the alpha value?
  • Is there a highlight for the column and/or row?
  • Which column and/or row should be highlighted?
  • What are the colours and alpha values for each highlight?

Mine is of course a fairly naive implementation,hardcoding certain coordinate adjustments instead of dynamically scaling them somehow (don’t know how) and only useful when there are not too many/little rows and columns. In any case,usage is fairly simple :

  1. Specify elements using the ggTableSpec function.
  2. Draw table using the ggTableDrawer function.

l

Example 1 : No highlighting whatsoever

l

  smpl.data <- data.frame(Col1=round(runif(20),3),Col2=round(runif(20),3),Col3=round(runif(20),3),Col4=round(runif(20),3),Col5=round(runif(20),3))
  rownames(smpl.data) <- paste('Row',1:20,sep='')

	smpl.spec <- ggTableSpec(columns.exist=T,columns.txt=colnames(smpl.data),columns.font=rep('bold',ncol(smpl.data)),columns.col='blue',columns.fill='grey',columns.alpha=0.7,
		                      rows.exist=T,rows.txt=rownames(smpl.data),rows.font=rep('bold',10),rows.col='green',rows.fill='blue',rows.alpha=0.7,
		                      data.obj=smpl.data,
		                    	hlt.col.exist=F,hl.col.which=c(1,5,20),hl.col.fill=c('lightgreen','darkgreen','red'),hl.col.alpha=c(0.4,0.4,0.4),
		                    	hlt.row.exist=F,hl.row.which=c(1,2,5),hl.row.fill=c('skyblue','red','yellow'),hl.row.alpha=c(0.4,0.4,0.4)
                )
  ggTableDrawer(smpl.spec)

a1

 

Example 2 : Column Highlighting
l

  smpl.data <- data.frame(Col1=round(runif(20),3),Col2=round(runif(20),3),Col3=round(runif(20),3),Col4=round(runif(20),3),Col5=round(runif(20),3))
  rownames(smpl.data) <- paste('Row',1:20,sep='')

	smpl.spec <- ggTableSpec(columns.exist=T,columns.txt=colnames(smpl.data),columns.font=rep('bold',ncol(smpl.data)),columns.col='blue',columns.fill='grey',columns.alpha=0.7,
		                      rows.exist=T,rows.txt=rownames(smpl.data),rows.font=rep('bold',10),rows.col='green',rows.fill='blue',rows.alpha=0.7,
		                      data.obj=smpl.data,
		                    	hlt.col.exist=F,hl.col.which=c(1,5,20),hl.col.fill=c('lightgreen','darkgreen','red'),hl.col.alpha=c(0.4,0.4,0.4),
		                    	hlt.row.exist=T,hl.row.which=c(1,2,5),hl.row.fill=c('skyblue','red','yellow'),hl.row.alpha=c(0.4,0.4,0.4)
                )
  ggTableDrawer(smpl.spec)

a2

Some alignment issues for the last column there.

 

Example 3 : Row Highlighting
l

  smpl.data <- data.frame(Col1=round(runif(20),3),Col2=round(runif(20),3),Col3=round(runif(20),3),Col4=round(runif(20),3),Col5=round(runif(20),3))
  rownames(smpl.data) <- paste('Row',1:20,sep='')

	smpl.spec <- ggTableSpec(columns.exist=T,columns.txt=colnames(smpl.data),columns.font=rep('bold',ncol(smpl.data)),columns.col='blue',columns.fill='grey',columns.alpha=0.7,
		                      rows.exist=T,rows.txt=rownames(smpl.data),rows.font=rep('bold',10),rows.col='green',rows.fill='blue',rows.alpha=0.7,
		                      data.obj=smpl.data,
		                    	hlt.col.exist=T,hl.col.which=c(1,5,20),hl.col.fill=c('lightgreen','darkgreen','red'),hl.col.alpha=c(0.4,0.4,0.4),
		                    	hlt.row.exist=F,hl.row.which=c(1,2,5),hl.row.fill=c('skyblue','red','yellow'),hl.row.alpha=c(0.4,0.4,0.4)
                )
  ggTableDrawer(smpl.spec)

a3

Again some alignment problems at the top and bottom. Maybe I will correct these issues later,for the time being I do not care. My suspicion is that the extreme coordinates for the column and row title panels are -/+ Inf whereas the coordinates for the highlights are not,hence their alignment should not be expected. Changing the Inf values to concrete ones should do the trick (maybe).

 

Example 4 : Row and Column Highlighting
l

 smpl.data <- data.frame(Col1=round(runif(20),3),Col2=round(runif(20),3),Col3=round(runif(20),3),Col4=round(runif(20),3),Col5=round(runif(20),3))
  rownames(smpl.data) <- paste('Row',1:20,sep='')

	smpl.spec <- ggTableSpec(columns.exist=T,columns.txt=colnames(smpl.data),columns.font=rep('bold',ncol(smpl.data)),columns.col='blue',columns.fill='grey',columns.alpha=0.7,
		                      rows.exist=T,rows.txt=rownames(smpl.data),rows.font=rep('bold',10),rows.col='green',rows.fill='blue',rows.alpha=0.7,
		                      data.obj=smpl.data,
		                    	hlt.col.exist=T,hl.col.which=c(1,5,20),hl.col.fill=c('lightgreen','darkgreen','red'),hl.col.alpha=c(0.4,0.4,0.4),
		                    	hlt.row.exist=T,hl.row.which=c(1,2,5),hl.row.fill=c('skyblue','red','yellow'),hl.row.alpha=c(0.4,0.4,0.4)
                )
  ggTableDrawer(smpl.spec)

l
a5Again same issues as above.Do not really see when i would want to highlight so many rows and columns anyways.

 

Example 5 : Problem

As an example of a problem consider the case where columns and rows do not exist but highlight do :
l

	smpl.spec <- ggTableSpec(columns.exist=F,columns.txt=colnames(smpl.data),columns.font=rep('bold',ncol(smpl.data)),columns.col='blue',columns.fill='grey',columns.alpha=0.7,
		                      rows.exist=F,rows.txt=rownames(smpl.data),rows.font=rep('bold',10),rows.col='green',rows.fill='blue',rows.alpha=0.7,
		                      data.obj=smpl.data,
		                    	hlt.col.exist=T,hl.col.which=c(1,5,20),hl.col.fill=c('lightgreen','darkgreen','red'),hl.col.alpha=c(0.4,0.4,0.4),
		                    	hlt.row.exist=T,hl.row.which=c(1,2,5),hl.row.fill=c('skyblue','red','yellow'),hl.row.alpha=c(0.4,0.4,0.4)
                )
  ggTableDrawer(smpl.spec)

l
dd

Clearly there are many issues that still need to be cleared up,but for the purposes of replicating the basic results of the paper I will leave the functions as they are for now. Instead of having me fumble around, I wish that the ggplot2 package creator would add table plotting functions to his package !

The final task we have set ourselves, and which also happens to be a welcome initiation to the excellent ggplot2 package for me, is concerned with summarising the previously saved plots and tables into a single (and hopefully useful) dashboard. I have attempted to reduce clutter as much as possible by omitting axis ticks,values and legends in those cases where interpretation is forthcoming. After exhausting the template themes available from the ggthemes package, I have settled on the economist theme for all the plots in the dashboard,finding it to be the most pleasing on the eyes. The function that collects and arranges all the outputs from previous steps is the custom DashboardPlot() function which accepts as arguments the following :[1] List object of dashboard elements,[2]Fund name of choice ,[3] Regime name of choice.
l

#[7.3.3] Dashboard
dashboard.list <- list(FundRet.list=edhec.ret,
	                     Turbulence.plot=Turbulence.plot,
	                     Turbulence1.tbl=Turbulence1.tbl,
	                     Turbulence2.tbl=Turbulence2.tbl,
	                     plot.list=plot.list,
	                     CompareCumu=CompareCumu,
	                     FacetPlot=FacetPlot.list,
	                     FacetObj=Facet.obj,
	                     MarkovPlots=MarkovEstPlots)

#Dashboard for Global Macro fund in Inflation regime
# DashboardPlot(Inputs=dashboard.list,
# 	            Fund='Global.Macro',
# 	            Regime='Inflation')

#Dashboard for all regimes/fund combinations
for(i in fund.names){
	for(j in regimes.names){
		DashboardPlot(Inputs=dashboard.list,
			            Fund=i,
		              Regime=j
		)
	}
}
########################################################################################

l
This results in n.funds(13)*n.regimes(4)= 52 pdf files, each representing a unique combination of a particular fund and a particular regime.

I have uploaded the pdf files into wordpress and they should be available for viewing by clicking on the links in the following table :

l

Funds | Regimes Equity Currency Inflation Growth
Convertible Arbitrage C&EQ C&FX C&Inf C&G
CTA Global CTA&EQ CTA&FX CTA&Inf CTA.&G
Distressed Securities D&EQ D&FX D&Inf D&G
Emerging Markets Em&EQ Em&FX Em&Inf Em&G
Equity Market Neutral Eq&EQ Eq&FX Eq&Inf Eq&G
Event Driven Ev&EQ Ev&FX Ev&Inf Ev&G
Fixed Income Arbitrage Fix&EQ Fix&FX Fix&Inf Fix&G
Funds of Funds Fun&EQ Fun&FX Fun&Inf Fun&G
Global Macro Glo&EQ Glo&FX Glo&Inf Glo&G
Long Short Equity Lon&EQ Lon&FX Lon&Inf Lon&G
Merger Arbitrage Me&EQ Me&FX Me&Inf Me&G
Relative Value Rel&EQ Rel&FX Rel&Inf Rel&G
Short Selling Sh&EQ Sh&FX Sh&Inf Sh&G

l

This html table was generated using The Tables Generator with Compact mode ticked and Do not Generate CSS unticked.

To give an example of the information contained in a dashboard,consider the case of  Global Macro x Inflation and relevant points of interest below :

l

mkl

l

Points of interest :

[Barcharts & Table of estimates in top rows]

  1. The Inflation economic regime variable begins in the normal state,with a high probability of remaining in this state (Persistence : 90.1%) and a low probability of transitioning to the event state (Transition : 9.9%).The estimated mean and sigma for this normal state are 0.2559 and 0.275 respectively.
  2. The Event state for all regime variable is characterised by higher estimated sigmas and means than corresponding values for the normal state.
  3. The persistence of normal states across all economic regime variables is higher than that that of event states. The corollary is found in transition estimates where we are more likely to transition from the event state to the normal state than from the normal state to the event state for all regime variables.
  4. The dodged barchart in the final column seeks to answer the question : What is the average return for fund x when regime y is in state a or state b? For our example, the mean return to the Global Macro strategy when the inflation regime is in the normal state exceeds the corresponding value in the event state.

 

[Selection & Performance plots in bottom rows]

  1. The list of funds and regime variables are given as strips of text in two separate ggplots with a green rectangle beneath the chosen fund and a red rectangle beneath the chosen regime. They are here mainly for eye candy, filling in space and horsing around with the ggplot package!
  2. The monthly,cumulative returns and drawdowns are also drawn.It seems that that for the Global Macro x Inflation case,knowledge of Markov states does not translate to higher cumulative returns.The drawdowns however seem more favourable.
  3. We also have a time series plot of the chosen economic regime variable along with an overlay of event versus non-event states.The thin red bars signify moments in time where the regime variable is in the event state.The blue bars show when that variable is in its normal state.There are some strange discrepancies between the pdf file and this zoomed image above.The pdf file is a better option here even though the colour is strangely off depending on the zoom.Also there is some sort of white taper in the pdf version of this plot which is probably due to the jagged lines I mentioned in the previous post.
  4. The final barchart in the bottom right corner shows the scaled difference in means across states for the chosen regime variable. If I understood this correctly,this should answer the question : By how many standard deviations is a particular fund’s performance higher/lower during the event/normal regime? In our case,the Global Macro strategy performs worse in the event regime than the normal regime by -0.10 standard deviations (?) I suppose this corroborates the previous finding that the mean performance of the chosen fund,when seen in the context of the inflation regime variable,is superior in the normal state versus the event state.

Whether these results are consistent or my interpretations correct I do not know.Any mistakes are of course mine.Although the code for the function is a bit of a mess and for the most part just arranges the plots in the desired layout,I include it here for the sake of completeness :
l

#########################################################################################
# Dashboard
#########################################################################################

DashboardPlot <- function(Inputs,Fund,Regime){

#Set up data
	regime.map <- data.frame(regime=c('Equity','Currency','Inflation','Growth'),idx=c(1:4),stringsAsFactors=F)
	fund.names <- names(Inputs$FundRet.list)
  reg.names <- regime.map[,1]

#regime strip
  lbl.font <- c(rep('bold',4))
	reg.df<-data.frame(lab.x=rep(20,4),lab.y=seq(2,26,length=4),lab.txt=reg.names,lab.font=lbl.font,stringsAsFactors=F)
	empty.df <- data.frame(x=-2:26,y=-2:26)
	ymin <- filter(reg.df,lab.txt==Regime)$lab.y-2
	ymax<-ymin+diff(reg.df$lab.y)[1]-1.5
	rect.df <- data.frame(xmin=-4,xmax=28,ymin=ymin,ymax=ymax)

	gg.reg.list <-  ggplot(empty.df)+geom_blank()+theme_economist(base_size=5)+
                      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
   	 									panel.border = element_blank(), panel.background = element_blank(), axis.title.x = element_blank(),
    									axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
    									axis.ticks = element_blank(),legend.position='none')+
                   		geom_rect(data=rect.df,aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),fill='red',alpha=0.4)+
		                  geom_text(data=reg.df,aes(x=lab.x,y=(lab.y),label=lab.txt,fontface=lab.font,hjust=1,vjust=0),size=rel(3))

#fund strip
  lbl.font <- c(rep('bold',13))
	fund.df<-data.frame(lab.x=rep(20,13),lab.y=seq(0,26,length=13),lab.txt=fund.names,lbl.font=lbl.font,stringsAsFactors=F)
	empty.df <- data.frame(x=-2:26,y=-2:26)
	ymin <- filter(fund.df,lab.txt==Fund)$lab.y-0.5
	ymax<-ymin+diff(fund.df$lab.y)[1]-0.5
	rect.df <- data.frame(xmin=-4,xmax=28,ymin=ymin,ymax=ymax)
	gg.fund.list <-  ggplot(empty.df)+geom_blank()+theme_economist(base_size=5)+
                      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
   	 									panel.border = element_blank(), panel.background = element_blank(), axis.title.x = element_blank(),
    									axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
    									axis.ticks = element_blank(),legend.position='none')+
                   		geom_rect(data=rect.df,aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),fill='green',alpha=0.4)+
		                  geom_text(data=fund.df,aes(x=lab.x,y=(lab.y),label=lab.txt,fontface=lbl.font,hjust=1,vjust=0),size=rel(3))

#Data table
	reg.str <- paste(filter(regime.map,regime==Regime)[2],']]',sep='')
	tbl1 <- Inputs$Turbulence1.tbl
	tbl2 <- Inputs$Turbulence2.tbl
	initial.state <- tbl1[unlist(filter(regime.map,regime==Regime)[2]),1]

	mean.s1 <- round(tbl1[unlist(filter(regime.map,regime==Regime)[2]),4],3)
        mean.s2 <- round(tbl2[unlist(filter(regime.map,regime==Regime)[2]),4],3)
	mean.df <- data.frame(State=c('Normal','Event'),Mean=c(mean.s1,mean.s2),stringsAsFactors=F)

	sd.s1 <- round(tbl1[unlist(filter(regime.map,regime==Regime)[2]),5],3)
        sd.s2 <- round(tbl2[unlist(filter(regime.map,regime==Regime)[2]),5],3)
	sd.df <- data.frame(State=c('Normal','Event'),StDev=c(sd.s1,sd.s2),stringsAsFactors=F)

	#Check State and store fitted information
	if(initial.state=='State 1'){
		temp <-Inputs$Turbulence1.tbl
		state.type <-'Normal'
	}else{
		temp<-Inputs$Turbulence2.tbl
    state.type <- 'Event'
	}

	pers <- temp[unlist(filter(regime.map,regime==Regime)[2]),2]
	trans <- temp[unlist(filter(regime.map,regime==Regime)[2]),3]
	mean.f <- round(temp[unlist(filter(regime.map,regime==Regime)[2]),4],3)
	stddev <- round(temp[unlist(filter(regime.map,regime==Regime)[2]),5],3)

#Create Plots

	gg.ret.plot <- eval(parse(text=paste('Inputs$FundRet.list$',Fund,sep='')))
	gg.turb.plot <- eval(parse(text=paste('Inputs$Turbulence.plot[[',reg.str,sep='')))

	gg.turb.map1 <- eval(parse(text=paste('Inputs$Turbulence.map2[[',reg.str,'[[1]]',sep='')))
	gg.turb.map2 <- eval(parse(text=paste('Inputs$Turbulence.map2[[',reg.str,'[[2]]',sep='')))

	gg.avg.plot <- eval(parse(text=paste('Inputs$plot.list[[',reg.str,sep='')))
	gg.cumu.plot1 <- eval(parse(text=paste('Inputs$CompareCumu$',Fund,'[[',reg.str,'$plot[[1]]',sep='')))
	gg.cumu.plot2 <- eval(parse(text=paste('Inputs$CompareCumu$',Fund,'[[',reg.str,'$plot[[2]]',sep=''))) 

	gg.facet.plot <- eval(parse(text=paste('Inputs$FacetPlot$',Fund,sep='')))
	gg.mkv.mean <- Inputs$MarkovPlots[[1]]
	gg.mkv.sd <- Inputs$MarkovPlots[[3]]
	gg.mkv.p <- Inputs$MarkovPlots[[5]]
	gg.mkv.t <- Inputs$MarkovPlots[[6]]

#Dashboard table
	lab1 <- paste("~Fund   :  ",Fund,sep=' ')
	lab2 <- paste('~Regime   :  ',Regime,sep=' ')
        lab3 <- paste('       |__Initial State :  ',state.type,sep=' ')
        lab4 <- paste('       |__Persistence :  ',pers,sep=' ')
        lab5 <- paste('       |__Transition :  ',trans,sep=' ')
        lab6 <- paste('       |__Mean  :  ',mean.f,sep=' ')
        lab7 <- paste('       |__StDev  :  ',stddev,sep=' ')
        lbl.cols <- c('green','red','white','white','white','white','white')
        lbl.font <- c('bold','bold','bold','bold','bold','bold','bold')
        lbl.df <-data.frame(lab.x=rep(-2,7),lab.y=seq(-1,8,length=7),lab.txt=c(lab1,lab2,lab3,lab4,lab5,lab6,lab7),lbl.font=lbl.font,lbl.cols=lbl.cols,stringsAsFactors=F)

	empty.df <- data.frame(x=-2:12,y=-2:12)
	gg.empty.plot <-  ggplot(empty.df,aes(x=x,y=y))+geom_blank()+geom_hline(y=11,colour='white',size=1)+theme_tufte(base_size=5)+
                      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
   	 									panel.border = element_blank(), panel.background = element_blank(), axis.title.x = element_blank(),
    									axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
    									axis.ticks = element_blank(),plot.title = element_text(size = rel(4),colour='white'),plot.background = element_rect(colour = "skyblue4", fill = "skyblue4"))+
		                  labs(title='\n.:: Dashboard ::.')+
		                  geom_text(data=lbl.df,aes(x=lab.x,y=rev(lab.y),label=lab.txt,hjust=0,vjust=1,fontface=lbl.font),size=rel(3),colour=lbl.cols)

	gg.empty2.plot <-  ggplot(empty.df,aes(x=x,y=y))+geom_blank()+theme_tufte(base_size=5)+
                      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
   	 									panel.border = element_blank(), panel.background = element_blank(), axis.title.x = element_blank(),
    									axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
    									axis.ticks = element_blank(),plot.title = element_text(size = rel(4)),plot.background = element_rect(colour = "skyblue4", fill = "skyblue4"))

#Data table
  temptbl <- eval(parse(text=paste('Inputs$FacetObj$',Fund,'$df',sep='')))
	gg.tbl.plot <- QuickTbl(temptbl,title='Averaging Returns')

#Markov estimate tables
	gg.mkv.mean.tbl <- Inputs$MarkovPlots[[2]]
   	gg.mkv.sd.tbl <- Inputs$MarkovPlots[[4]]
	gg.mkv.p.tbl <- Inputs$MarkovPlots[[7]]
        gg.mkv.t.tbl <- Inputs$MarkovPlots[[8]]

#Arrange multiple plots and save as pdf
pdf(file = paste(Fund,'&',Regime,'.pdf',sep=''), width = 30, height = 17)
  layOut(
  	     list(gg.empty.plot,1:2,1),
  	     list(gg.empty2.plot,3:6,1),
  	     list(gg.reg.list,3,2),
  	     list(gg.fund.list,4:6,2),
   	  	 list(gg.mkv.mean,1,2),
  	     list(gg.mkv.sd,1,3),
  	     list(gg.mkv.p,1,4),
  	     list(gg.mkv.t,1,5),
  	  	 list(gg.facet.plot,1,6),
  	     list(gg.mkv.mean.tbl,2,2),
  	     list(gg.mkv.sd.tbl,2,3),
         list(gg.mkv.p.tbl,2,4),
  	     list(gg.mkv.t.tbl,2,5),
  	     list(gg.tbl.plot,2,6),
  	     list(gg.ret.plot,3:4,3:4),
  	     list(gg.cumu.plot1,5,3:4),
   	     list(gg.cumu.plot2,6,3:4),
  	     list(gg.avg.plot,5:6,5:6),
  	  	 list(gg.turb.plot,3:4,5:6)
   	)
dev.off()
}
#########################################################################################

l
Wish I could make this collapsable but wordpress seems to need a plugin for that.

Now that we have (Gods willing) aligned fund returns,regime indices and state probabilities by time,we can finally turn to the fifth issue in the previously defined list and examine in-sample performance along the following lines :

  1. Standard-deviation-Scaled difference across state means (as per the paper)
  2. Cumulative returns and drawdowns when states are known versus unknown
  3. Event-state mean vs non-event-state mean for the chosen fund across all economic regimes.

 

[5] In-Sample Performance

To every element in the Aligned list object from the end of the previous post,
l

#[7.3] Fund Performance : In-sample
 Regimes <- list()
 Regimes <- lapply(Aligned,function(x) FundPerf(x))

 plot.list <- list()
 plot.list <- lapply(Regimes,function(x) FundMeanPlot(x))

l
I applied the custom FundPerf() function to calculate among other things, the scaled differences in means.
l

#########################################################################################
# Fund Performance
#########################################################################################
FundPerf <- function(aligned.list.element){

	col.turb <- names(aligned.list.element)[1]
	aligned.sub <- select(aligned.list.element,eval(parse(text=col.turb)),Normal.Class:Funds.Of.Funds)
 	std.smpl <- apply(aligned.sub[,-c(1:3)],2,sd)

	event.ret <- select(aligned.sub,Convertible.Arbitrage:Funds.Of.Funds)*aligned.sub[,'Event.Class']
  event.mean <- colMeans(event.ret[which(aligned.sub[,'Event.Class']==1),])

	nonevent.ret <- select(aligned.sub,Convertible.Arbitrage:Funds.Of.Funds)*aligned.sub[,'Normal.Class']
  nonevent.mean <- colMeans(nonevent.ret[which(aligned.sub[,'Normal.Class']==1),])

	result.df <- data.frame(Fund=colnames(aligned.list.element[-c(1:5)]),EventMean=event.mean,NonEventMean=nonevent.mean,Stdev=std.smpl)
	fund.result <- mutate(result.df,Perf=(event.mean-nonevent.mean)/std.smpl)

	return(list(Aligned.df = aligned.sub,
		     			Event.df = event.ret,
		     			NonEvent.df = nonevent.ret,
		     			FundPerf.df = fund.result
		     )
	)
}
#########################################################################################

l
To give a sense of the output that this function returns. the following shows a snapshot of the returned list object for the Regimes$Growth case :

  • A subset of the original aligned data frame now containing only the time series of the regime variable index,state maps and fund returns.

d1

  •  A data frame which contains the time series of fund returns multiplied by a vector (of zeros and ones) corresponding to the event state.

d2

  • A data frame which contains the time series of fund returns multiplied by a vector (of zeros and ones) corresponding to the normal state.

d3

  • A data frame which contains the fund names,the event mean,the non event mean and the final scaled value (given in the paper as event mean minus non event mean divided by standard deviation)

d4

 

The FundMeanPlot() function is then used to generate a ggplot on the basis of this output  :
l

#########################################################################################
# Fund Perormance event/non event mean plot
#########################################################################################
FundMeanPlot <- function(regime.obj){

	localenv <- environment()

	dfdata <- regime.obj$FundPerf.df
	regime.name <- names(regime.obj$Aligned.df)[1]

	min.ret <- round(min(dfdata[,'Perf']),1)
	max.ret <- round(max(dfdata[,'Perf']),1)

	ggp <- ggplot(dfdata,aes(x=Fund,y=Perf),environment=localenv)+
				    theme_economist(base_size=5)+
						geom_bar(stat='identity',position=position_dodge(),aes(fill=(dfdata$Perf)))+
						scale_y_continuous(limits=c(min.ret-0.2,max.ret+0.2),breaks=seq(min.ret-0.2,max.ret+0.2,length=5))+
						labs(x='',title=paste('Scaled Difference in Means for',regime.name,'regime',sep=' '),y='Values')+
						theme(legend.position = "none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
						coord_flip()

return(ggp)
}
#########################################################################################

l
In addition to the in-sample scaled performance score calculated above,I also tried to look at how cumulative returns varied depending on one’s knowledge of the estimated states in each regime. I simply multiplied the chosen fund return by the state map for each economic regime variable (so that fund returns are multiplied by 1 in the case of a normal state and 0 of an event state),calculated the cumulative return as well as drawdowns for that series (using the PerformanceAnalytics package) and compared the results to the case where states are unknown.

I don’t know if these steps make sense or if my implementation is correct but that is what i did in any case. Although the data for this is already available from previous code segments,I purposefully started afresh so I can toy with it more easily (eg.multiplying by -1 in the event state as opposed to 0). Here is the code to create the data frame to be passed to the custom CumuCalc() function.
l

#[7.3.2] Comparing Cumulative returns when states are known vs unknown

n.regimes <- num.turbs
regimes.names <- names.turbs
fund.names <- names(edhec.xts)

CompareCumu <- list()
for(i in 1:n.funds){

	CompareCumu[[i]] <- list()
	names(CompareCumu)[i] <- fund.names[i]

	for(j in 1:n.regimes){
		CompareCumu[[i]][[j]] <- list()
		names(CompareCumu[[i]])[[j]] <- regimes.names[j]

		CompareCumu[[i]][[j]]$name.str <- paste(regimes.names[j],'::',fund.names[i],'Fund')
		CompareCumu[[i]][[j]]$data.orig <- data.frame(TimeIdx   = as.Date(rownames(Regimes[[j]]$Aligned.df)),
																						 FundReturn     = Regimes[[j]]$Aligned.df[,fund.names[i]],
																						 State			    = Regimes[[j]]$Aligned.df$Normal.Class
																				)
		temp.store <- CompareCumu[[i]][[j]]$data.orig

		CompareCumu[[i]][[j]]$data.new <- mutate(temp.store,State.new=ifelse(State==1,1,0),FundReturn.new=State.new*FundReturn)
		CompareCumu[[i]][[j]]$plot <- CumuCalc(CompareCumu[[i]][[j]])
	}
}

l
The most important issue here is the call to the CumuCalc() function which returns a ggplot of cumulative returns and drawdowns for the chosen fund/regime combination.The nested list structure of the CompareCumu variable alows us to [1] loop over all funds and [2] within each fund loop over all regime variables,simplifying data access to a command like : CompareCumu$Short.Selling$Inflation for example.

For the sake of completeness,the code for that function is given here :
l

#########################################################################################
# Cumulative  Calculation and ggplots saved
#########################################################################################

CumuCalc <- function(fund.df){

	tindex <- as.Date(fund.df$data.new$TimeIdx)
	fund.names <- fund.df$name.str
	orig.ret <- fund.df$data.new[,'FundReturn']
	orig.cumu.ret <- cumprod(orig.ret+1)-1

	state.ret <- fund.df$data.new[,'FundReturn.new']
	state.cumu.ret <- cumprod(state.ret+1)-1

	orig.dd <- Drawdowns(orig.ret)
	state.dd <- Drawdowns(state.ret)

	ggplot.df <- data.frame(TimeIdx  = tindex,
		                      Ret.Original = orig.cumu.ret,
		                      Ret.Markov   = state.cumu.ret,
													DD.Original  = orig.dd,
													DD.Markov    = state.dd,
		                      stringsAsFactors=F
	             )
	lbl.df <- data.frame(lbl.x=rep(max(ggplot.df$TimeIdx)-300,2),lbl.y=c(last(ggplot.df$Ret.Original),last(ggplot.df$Ret.Markov)),lbl.txt=c('Original','Markov'),stringsAsFactors=F)
	ggp1 <- ggplot(data=ggplot.df,aes(x=TimeIdx))+
 				  theme_economist(base_size=5)+
		      geom_line(aes(x=TimeIdx,y=Ret.Original,colour='Original'),size = 1)+
				  geom_line(aes(x=TimeIdx,y=Ret.Markov,colour='Markov'),size = 1)+
          scale_colour_manual("",breaks = c('Markov','Original'),values=c('skyblue4','cyan'))+
				  labs(x='',y='Values')+
		      theme(legend.position=c(0.3,1),legend.key.size=unit(0.2, "cm"),legend.key.width=unit(2,'cm'),legend.direction='horizontal',panel.grid.major = element_blank(), panel.grid.minor = element_blank())

	ggp2 <- ggplot(data=ggplot.df,aes(x=TimeIdx))+
 				  theme_economist(base_size=5)+
		      geom_line(aes(x=TimeIdx,y=DD.Original,colour='Original'),size=1)+
				  geom_line(aes(x=TimeIdx,y=DD.Markov,colour='Markov'),size=1)+
 		      scale_colour_manual("Legend : ",breaks = c('Markov','Original'),values=c('skyblue4','cyan'))+
				  labs(x='Time',y='Values')+
		      theme(legend.position='none',panel.grid.major = element_blank(), panel.grid.minor = element_blank())

	return(list(ggp1,ggp2))

}
####################################################################################################

l
This function simply returns a list object with two elements, ggplots of cumulative returns and drawdowns for the chosen fund and regime combination when states are known versus when they are unknown.

A final issue I tried to address,so as to make the dashboard less empty, was to compare the event-state mean with its non event counterpart for the chosen fund and economic regime variable. Once again,although all the necessary data frames are already in place,I dealt with this issue separately to ease my pain :
l

#[7.3.4] Facet Plot Date frame and graph

Facet.obj <- list()
for(i in 1:n.funds){
	Facet.obj[[i]] <- list()
	names(Facet.obj)[i] <- fund.names[i]

	values.matrix <- NULL
	for(j in 1:n.regimes){
		temp <- select(Regimes[[j]]$FundPerf.df,c(Fund,EventMean:NonEventMean))
		filtered <- filter(temp,Fund==fund.names[i])[,-1]
		values.matrix <- rbind(values.matrix,filtered)
	}

	Facet.obj[[i]]$df <- data.frame(Regime=regimes.names,Event=values.matrix[,1],Normal=values.matrix[,2])

	Facet.obj[[i]]$reshaped.df <- gather(Facet.obj[[i]]$df,variable,value,-Regime)
	colnames(Facet.obj[[i]]$reshaped.df)[-1] <- c('State','Mean')
}

FacetPlot.list <- lapply(Facet.obj,function (x) FacetPlot(x))

l
The FacetPlot() function called at the end there simply returns a ggplot of ‘dodged’ barcharts of average returns across regimes and states for the chosen fund.
l

#########################################################################################
# Facet Plotter
#########################################################################################

FacetPlot <- function(plot.obj){
	plotdf <- plot.obj$reshaped.df

		ggbb <- ggplot(plotdf)+
        			theme_economist(base_size=5)+
						  geom_bar(data = plotdf,aes(x=Regime,y=Mean,fill=State),stat = "identity",position='dodge') +
		  				scale_fill_brewer(palette='Set1')+
							labs(title='Average returns across\nRegimes & States')+
      				coord_flip()+
      			  theme(axis.ticks=element_blank(),legend.position='none',axis.title=element_blank(),axis.text.x=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

	return(ggbb)
}
#########################################################################################

ll
Across these three posts,we have gathered all the inputs necessary for our final dashboard.Let’s save these into a single list and pass it to the DashboardPlot() function in the next post.
l

#[7.3.3] Dashboard
dashboard.list <- list(FundRet.list=edhec.ret,
	                     Turbulence.plot=Turbulence.plot,
	                     Turbulence1.tbl=Turbulence1.tbl,
	                     Turbulence2.tbl=Turbulence2.tbl,
	                     plot.list=plot.list,
	                     CompareCumu=CompareCumu,
	                     FacetPlot=FacetPlot.list,
	                     FacetObj=Facet.obj,
	                     MarkovPlots=MarkovEstPlots)

l

Continuing from the previous post where we downloaded,transformed (where necessary) and stored the data series for each of the two categories of variables (economic regime variables & risk premia),the remainder of this series of posts shall be concerned with the following tasks :

  1. Fit a two state hidden markov model to economic regime variables (Equity turbulence,Currency Turbulence,Inflation,Economic Growth).
  2. For each regime variable,overlay the resulting probability state map on a time series plot of the regime variable and save final graph as an element of a list.
  3. For each regime variable,extract from the fitted object a table of estimation results and store the table plot as an element of a list.
  4. For each hedge fund strategy,ensure time series are aligned across economic regime variables.
  5. Explore In-Sample Performance across funds :
    • Event-state Mean vs Normal-state Mean.
    • Scaled difference in state means.
    • Cumulative returns & Drawdowns when states are known vs unknown
  6. Create Dashboard of results for chosen Fund/Regime combination.

Although the code was not initially created with a dashboard in mind,it is probably for the best to just show the dashboard of plots and tables at the end as opposed to displaying the outputs of intermediate steps. 


 

[1] Fitting the hidden markov model

The Turbulence.list object holds our regime variables as elements of a list. To fit a two-state hidden markov model to each of the variables,I used the lapply function to pass each element of that list to a custom HMfit() function.
g

#[3] Fit the turbulence indices where available
	Turbulence.fit <- lapply(Turbulence.list,function(x) HMfit(x))

k
The HMfit() function calls the fitting functions in the depmixS4 package and returns a list object as given below :
j

#########################################################################################
# Hidden Markov Fit
#########################################################################################
HMfit <- function(inp.ret){

	ret.xts <- inp.ret
	states <- 2

	model.spec <- depmix(eval(parse(text=names(ret.xts)))~1,nstates=states,data=ret.xts)
	model.fit <- fit(model.spec)
	model.posterior <- posterior(model.fit)
        sum.frame <- summary(model.fit)
	normal.idx <- which(sum.frame[,1]>0 & sum.frame[,1]==min(sum.frame[,1]))
	event.idx <- index(sum.frame)[-normal.idx]

	HM.fit <- list(
		normal.state = normal.idx,
		event.state = event.idx,
		turb =ret.xts,
		spec =model.spec,
                fit = model.fit,
                summary = summary(model.fit),
		normal.prob = model.posterior[,normal.idx+1],
                normal.class = model.posterior[,'state']==normal.idx,
		event.prob = model.posterior[,event.idx+1],
		event.class = model.posterior[,'state']==event.idx
	)

	HM.fit$df = data.frame(Time.Idx=index(HM.fit$turb),
			       Turbulence=HM.fit$turb,
                               Normal.Prob=HM.fit$normal.prob,
			       Event.Prob=HM.fit$event.prob,
			       Normal.Class=ifelse(HM.fit$normal.class==T,1,0),
			       Event.Class=ifelse(HM.fit$event.class==T,1,0),
			       row.names=NULL)
  return(HM.fit)
}

l
The most useful structure this function returns is a data frame containing :a time stamp,the turbulence index,state probabilities and state maps.

 

[2] Plotting economic regime variables and overlay with state maps

For a reason I cannot understand, the state map which shows when a particular economic regime variable is in the event (or normal) state and which therefore should be a line plot of either 0 or 1, becomes jagged when it should be straight. To overlay the state map on top of a timeseries plot of an economic regime variable,I simply multiplied the event state map ( a time series of 1’s and 0s corresponding to when the regime variable is in an event or non-event state ) by the highest value of the regime variable and used the resulting index as an input to the geom_ribbon layer. As before I use the lapply function to pass the previously fitted objects to a custom function TurbPlot().

l

#[4] Plot the Turbulence indices and overlay state maps
	Turbulence.plot <- lapply(Turbulence.fit,function(x) TurbPlot(x))

j
The TurbPlot() function creates one ‘overlay’ data frame for each state (event,normal) from the fitted object and then passes these to the ggplot function ( and most importantly the geom_ribbon layer ):
l

#########################################################################################
# Turbulence and regime Plotting
#########################################################################################
TurbPlot <- function(fitted.obj){

  plot.df <- fitted.obj$df
  ymax <- max(plot.df[,names(fitted.obj$turb)])+1

  overlay.df <- mutate(select(plot.df,Time.Idx,Event.Class),ID=Event.Class*ymax)
    overlay.df$yvals <- fitted.obj$turb
    overlay2.df <- mutate(select(plot.df,Time.Idx,Normal.Class),ID=Normal.Class*ymax)
    overlay2.df$yvals <- fitted.obj$turb

  gg.s1 <- ggplot()+
            theme_tufte(base_size=11)+
            geom_ribbon(data=overlay.df,aes(x=Time.Idx,ymin=0,ymax=ID),alpha=0.6,fill='red')+
            geom_ribbon(data=overlay2.df,aes(x=Time.Idx,ymin=0,ymax=ID),alpha=0.6,fill='lightblue')+
            geom_area(data=plot.df,aes_string(x='Time.Idx',y=names(fitted.obj$turb)),fill='skyblue4',colour='darkblue',size=0.5)+
            labs(x='Time',y='Values',title=paste('Turbulence Index For',names(fitted.obj$turb),'Regime',sep=' '))

return(gg.s1)
}
#########################################################################################

l
The function returns a ggplot to be used in the final dashboard.

 

[3] Construct table of estimation results

Continuing with our list of fitted objects,let’s extract pertinent information from it across the event and normal states and make a simple text plot of mean and sigma estimates across all four regime variables.
;l

#[6] Table of Estimation results
  Turbulence1.tbl <- Extract.est(Turbulence.fit,state='Normal')
  Turbulence2.tbl <- Extract.est(Turbulence.fit,state='Event')
  MarkovEstPlots <- MarkovEstPlot(Turbulence1.tbl,Turbulence2.tbl)

ll
Since we want a table that collects information across regimes,we do not need lapply to act on every element of the list and instead pass the entire list to a custom Extract.est() function.

j

#########################################################################################
#Extracting Estimates from list object
#########################################################################################
Extract.est <- function(fitted.list,state){
 list.store <- list()
 idx.str<-paste('ifelse(state==','"Normal"',',','x$normal.state',',','x$event.state)',sep='')
 list.store <- lapply(fitted.list,function(x) data.frame(
                      Initial_State = paste('State',which(x$fit@init==1)),
  		      Persistence   = paste(round(x$fit@transition[[eval(parse(text=idx.str))]]@parameters$coefficients[eval(parse(text=idx.str))],3)*100,'%',sep=''),
  		      Transition    = paste(round(x$fit@transition[[eval(parse(text=idx.str))]]@parameters$coefficients[-eval(parse(text=idx.str))],3)*100,'%',sep=''),
  		      Mean          = x$summary[eval(parse(text=idx.str)),1],
                      Std_Dev	    = x$summary[eval(parse(text=idx.str)),2],
                      stringsAsFactors=F
                    ))
  est.tbl <- do.call(rbind,list.store)
  return(est.tbl)
}
#########################################################################################

f
This function returns a data frame that contains information on (1)The initial state of the regime variable,(2)The probability of staying in one state,(3) The probability of transitioning to another state,(3) The mean estimate and (4) The sigma estimate across all regimes.

Once we have extracted the desired data from our fitted objects,I use the MarkovEstPlot() function to make relevant stacked bar charts and associated value tables to be used in the final dashboard.

l

#########################################################################################
# Markov Estimates Plotter
#########################################################################################

MarkovEstPlot <- function (tbl.norm,tbl.event){

	tbl.norm <- coredata(tbl.norm)
	tbl.event <- coredata(tbl.event)
	
	mean.df <- data.frame(c(rownames(tbl.norm)),select(tbl.event,Mean),select(tbl.norm,Mean),stringsAsFactors=F)
	colnames(mean.df) <- c('Regime','Event','Normal')
	mean.tbl <- mean.df
	
	sd.df <-data.frame(c(rownames(tbl.event)),select(tbl.event,Std_Dev),select(tbl.norm,Std_Dev),stringsAsFactors=F)
	colnames(sd.df) <- c('Regime','Event','Normal')
	sd.tbl <- sd.df
	
	pers.df <- data.frame(c(rownames(tbl.norm)),select(tbl.event,Persistence),select(tbl.norm,Persistence),stringsAsFactors=F)
  colnames(pers.df) <- colnames(mean.df)
	pers.df <-data.frame(Regime=c(rownames(tbl.norm)),apply(select(pers.df,Event,Normal),2,function(x) round(as.numeric(sub('%','',x)),4)),stringsAsFactors=F)
	pers.tbl <- pers.df
	
	trans.df <- data.frame(c(rownames(tbl.norm)),select(tbl.event,Transition),select(tbl.norm,Transition),stringsAsFactors=F)
  colnames(trans.df) <- colnames(mean.df)
	trans.df <- data.frame(Regime=c(rownames(tbl.norm)),apply(select(trans.df,Event,Normal),2,function(x) round(as.numeric(sub('%','',x)),4)),stringsAsFactors=F)
	trans.tbl <- trans.df
	
	mean.df <- gather(mean.df,variable,value,-Regime)
	sd.df<-gather(sd.df,variable,value,-Regime)
	pers.df<-gather(pers.df,variable,value,-Regime)
	trans.df<-gather(trans.df,variable,value,-Regime)
	colnames(mean.df)[2] <- colnames(sd.df)[2] <-colnames(pers.df)[2] <- colnames(trans.df)[2] <-'State' 
	
	
	gg.m <- ggplot()+
		        theme_economist(base_size=5)+
		  			geom_bar(data = mean.df,aes(x=Regime,y=value,fill=State),stat = "identity",position='dodge')+
		  			scale_fill_brewer(palette='Set1')+
						labs(title='Mean Estimates Across\nRegimes & States')+
		  			coord_flip()+
      			theme(axis.ticks=element_blank(),legend.position='none',axis.title=element_blank(),axis.text.x=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

	
	gg.sd <- ggplot()+
				    theme_economist(base_size=5)+
		  			geom_bar(data = sd.df,aes(x=Regime,y=value,fill=State),stat = "identity",position='dodge')+
		  			scale_fill_brewer(palette='Set1')+
						labs(title='Sigma Estimates Across\nRegimes & States')+
      			coord_flip()+
      			theme(axis.ticks=element_blank(),legend.position='none',axis.title=element_blank(),axis.text.x=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

	  gg.m.tbl <- QuickTbl(mean.tbl,title='Estimates of Mean')
	  gg.sd.tbl <-QuickTbl(sd.tbl,title='Estimates of Sigma')
	
		gg.p <- ggplot()+
				    theme_economist(base_size=5)+
		  			geom_bar(data = pers.df,aes(x=Regime,y=value,fill=State),stat = "identity",position='dodge')+
		  			scale_fill_brewer(palette='Set1')+
					  labs(title='Persistence Estimates Across\nRegimes & States')+
      			coord_flip()+
      		  theme(axis.ticks=element_blank(),legend.position='none',axis.title=element_blank(),axis.text.x=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank())


		gg.t <- ggplot()+
				    theme_economist(base_size=5)+
		  			geom_bar(data = trans.df,aes(x=Regime,y=value,fill=State),stat = "identity",position='dodge')+
		  			scale_fill_brewer(palette='Set1')+
						labs(title='Transition Estimates Across\nRegimes & States')+
      		  coord_flip()+
      			theme(axis.ticks=element_blank(),legend.position='none',axis.title=element_blank(),axis.text.x=element_blank(),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

		gg.p.tbl <- QuickTbl(pers.tbl,title='Estimates of Persistence (%)')
	  gg.t.tbl <-QuickTbl(trans.tbl,title='Estimates of Transition (%)')
	
						
	return(list(gg.m,gg.m.tbl,gg.sd,gg.sd.tbl,gg.p,gg.t,gg.p.tbl,gg.t.tbl))
}
##########################################################################################

hk
This function returns a list object with 4 stacked bar chart plots and 4 text plots corresponding to mean,sigma,persistence and transition estimates across regimes.

 

[4] Aligning risk premia and regime variable time series

Since returns on various hedge fund strategy returns may have different time indices than those of our regime variables,it is first necessary to ensure that we are looking at data that is correctly aligned in time. This is where mistakes may easily creep up (at least for me) and if there is anything wrong in subsequent outputs,the culprit is likely to lurk in this step. The Idea is to align returns across 13 funds with economic regime variables and store the output in a list object named Aligned.

l

#[7.2] Create list object to hold Turbulence aligned edhec return data
num.turbs <- length(Turbulence.fit)
names.turbs <- names(Turbulence.fit)
wind.beg <- min(index(edhec.xts))
wind.end <- max(index(edhec.xts))

Aligned <- list()
  for(i in 1:num.turbs){
    	Aligned[[i]] <- list()
        	names(Aligned)[i] <- names.turbs[i]
		sub.turb <- filter(Turbulence.fit[[i]]$df,Time.Idx>=wind.beg & Time.Idx<=wind.end)
		min.t <- paste(format(min(sub.turb$Time.Idx),'%Y-%m'),'-01',sep='')
		max.t <- paste(format(max(sub.turb$Time.Idx),'%Y-%m'),'-01',sep='')
		if(names.turbs[i]=='Growth'){
         		edhec.sub <- edhec.xts[sub.turb$Time.Idx,]
		}else{
			edhec.sub <- window(edhec.xts,start=min.t,end=max.t)
		}
	Aligned[[i]] <- cbind(sub.turb[,-1],as.data.frame(edhec.sub))
	}

l

Each element of the Aligned list is named after an economic regime variable and contains  a data frame of aligned fund returns. To get an idea as to the contents of the Aligned list above, let us look at  head(Aligned$Equity)

 sdsd

 

The columns continue for the remaining strategies.Data frames for other aligned time series are accessed by Aligned$Currency,Aligned$Inflation,Aligned$Growth. Now that we have aligned fund returns,regime indices and state probabilities,we are in the position to explore the in-sample performance of each fund.