Regime Detection Update

library(SIT)
load.packages('quantmod')

###############################################################################
# Regime Detection
# http://blogs.mathworks.com/pick/2011/02/25/markov-regime-switching-models-in-matlab/
###############################################################################
#bt.regime.detection.test <- function() 

	#*****************************************************************
	# Generate data as in the post
	#****************************************************************** 
	bull1 = rnorm( 100, 0.10, 0.15 )
	bear  = rnorm( 100, -0.01, 0.20 )
	bull2 = rnorm( 100, 0.10, 0.15 )
	true.states = c(rep(1,100),rep(2,100),rep(1,100))
	returns = c( bull1, bear,  bull2 )

	# find regimes
	load.packages('RHmm', repos ='http://R-Forge.R-project.org')

	y=returns
	ResFit = HMMFit(y, nStates=2)
	VitPath = viterbi(ResFit, y)

	# HMMGraphicDiag(VitPath, ResFit, y)
	# HMMPlotSerie(y, VitPath)

	#Forward-backward procedure, compute probabilities
	fb = forwardBackward(ResFit, y)

	# Plot probabilities and implied states
	layout(1:2)
	plot(VitPath$states, type='s', main='Implied States', xlab='', ylab='State')
	
	matplot(fb$Gamma, type='l', main='Smoothed Probabilities', ylab='Probability')
		legend(x='topright', c('State1','State2'),  fill=1:2, bty='n')

	# http://lipas.uwasa.fi/~bepa/Markov.pdf
	# Expected duration of each regime (1/(1-pii))                
	#1/(1-diag(ResFit$HMM$transMat))
          

	#*****************************************************************
	# It also can be done using depmixS4 package
	#[Getting Started with Hidden Markov Models in R](http://blog.revolutionanalytics.com/2014/03/r-and-hidden-markov-models.html)
	#****************************************************************** 
	load.packages('depmixS4')

	# Construct and fit a regime switching model
	mod = depmix(y ~ 1, family = gaussian(), nstates = 2, data=data.frame(y=y))
	fm2 = fit(mod, verbose = FALSE)

	probs = posterior(fm2)

	layout(1:2)
	plot(probs$state, type='s', main='Implied States', xlab='', ylab='State')
	
	matplot(probs[,-1], type='l', main='Probabilities', ylab='Probability')
		legend(x='topright', c('State1','State2'),  fill=1:2, bty='n')

	#*****************************************************************
	# Add some data and see if the model is able to identify the regimes
	#****************************************************************** 
	bear2  = rnorm( 100, -0.01, 0.20 )
	bull3 = rnorm( 100, 0.10, 0.10 )
	bear3  = rnorm( 100, -0.01, 0.25 )
	true.states = c(true.states, rep(2,100),rep(1,100),rep(2,100))
	y = c( bull1, bear,  bull2, bear2, bull3, bear3 )
	VitPath = RHmm:::viterbi(ResFit, y)$states

	# map states: sometimes HMMFit function does not assign states consistently
	# let's use following formula to rank states
	# i.e. high risk, low returns => state 2 and low risk, high returns => state 1
	map = rank(sqrt(ResFit$HMM$distribution$var) - ResFit$HMM$distribution$mean)
	VitPath = map[VitPath]

	#*****************************************************************
	# Plot regimes
	#****************************************************************** 
	data = xts(y, as.Date(1:len(y)))

	layout(1:3)
		plota.control$col.x.highlight = col.add.alpha(true.states+1, 150)
	plota(data, type='h', plotX=F, x.highlight=T)
		plota.legend('Returns + True Regimes')
	plota(cumprod(1+data/100), type='l', plotX=F, x.highlight=T)
		plota.legend('Equity + True Regimes')
	
		plota.control$col.x.highlight = col.add.alpha(VitPath+1, 150)
	plota(data, type='h', x.highlight=T)
		plota.legend('Returns + Detected Regimes')

#*****************************************************************
# Load historical prices
#****************************************************************** 
data = env()
getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)

price = Cl(data$SPY)
	open = Op(data$SPY)
ret = diff(log(price))
	ret = log(price) - log(open)

atr = ATR(HLC(data$SPY))[,'atr']



#*****************************************************************
# Construct and fit a regime switching model
#****************************************************************** 
load.packages('depmixS4')

index = 14:nrow(price)
temp = data.frame(ret=as.vector(ret), atr=as.vector(atr))
	temp = temp[index,]

# We're setting the LogReturns and ATR as our response variables, 
# want to set 4 different regimes, and setting the response distributions to be gaussian.
mod = depmix(list(ret~1, atr~1),data=temp,nstates=4,family=list(gaussian(),gaussian())) 
fm2 = fit(mod, verbose = FALSE)

print(summary(fm2))

probs = posterior(fm2)

print(head(probs))

layout(1:3)
plota(ret, type='h', col='blue')
	plota.legend('S&P 500 Daily Log Returns', 'blue')
plota(atr, type='l', col='darkgreen')
	plota.legend('Daily ATR(14)', 'darkgreen')

temp = NA * price
	temp[index] = probs$state
plota(temp, type='l', col='darkred')
	plota.legend('Market Regimes', 'darkred')

layout(1:4)
for(i in 2:5) {
	temp = NA * price
		temp[index] = probs[,i]
	plota(temp, type='l', col=i)
		plota.legend(paste('Market Regime', colnames(probs)[i]), i)
}