A mashup of financial turbulence and regime switching examples having missing bits into a standalone example without missing bits. Uses sources from Quantivity and Systematic Investor blogs as well as the CRAN RHmm and TTR packages. Uses quantmod and FRED as a data source. The turbulence calculation clearly is not the same as referenced original…

currency.png

equity.png

gdp.png

inflation.png

plot.png

regime_mashup.R

# Turbulence and regimes example mash-up. Standalone example with FRED data.
# Matt Barry <mrb@softisms.com>
#
# Original sources:
# http://quantivity.wordpress.com/2012/11/09/multi-asset-market-regimes/
# http://systematicinvestor.wordpress.com/2012/12/01/financial-turbulence-example/

library("RHmm")
library("TTR")

# load SIT for plota() and legends
# setInternet2(TRUE)
# con = gzcon(url('http://www.systematicportfolio.com/sit.gz','rb'))
con = gzcon(url('file://sit.gz','rb'))
source(con)
close(con)

displayKritzmanRegimes <- function()
{
  # DISPLAY REGIMES FROM KRITZMAN ET AL. (2012), PRINTING REGIME
  # STATISTICS AND PLOTTING LOCAL DECODING.
  
  equityTurbulenceRegime <- getEquityTurbulenceRegime()
  inflationRegime <- getInflationRegime()
  growthRegime <- getGrowthRegime()
  currencyTurbulenceRegime <- getCurrencyTurbulenceRegime()
  
  print(equityTurbulenceRegime)
  print(inflationRegime)
  print(growthRegime)
  print(currencyTurbulenceRegime)
  
  plotMarkovRegimes(equityTurbulenceRegime, "Equity (SPX)", plotDensity=F)
  plotMarkovRegimes(inflationRegime, "Inflation (CPIAUCNS)", plotDensity=F)
  plotMarkovRegimes(growthRegime, "Real GDP (GDPC1)", plotDensity=F)
  plotMarkovRegimes(currencyTurbulenceRegime, "G10 Currency Turbulence", plotDensity=F)
  
  plotLocalDecodings(list(equityTurbulenceRegime, growthRegime, inflationRegime, 
                          currencyTurbulenceRegime),
                     list("US Equity Turbulence (SPX)", "Real GDP (GDPC1)", 
                          "Inflation (CPIAUCNS)","G10 Currency Turbulence"),
                     regimeNums=c(2,2,2,2))
}

getEquityTurbulenceRegime <- function(startDate=as.Date("1977-12-01"), 
                                      endDate=Sys.Date(), numOfStates=2)
{
  # ESTIMATE TWO-STATE MARKOV (SPX-BASED) EQUITY REGIME. IN LIEU OF S&P 500
  # SECTOR INDICES, USE SPX INSTEAD.
  #
  # ARGS:
  #     STARTDATE: DATE WHICH TO BEGIN PANEL FOR REGIME ESTIMATION
  #     ENDDATE: END WHICH TO END PANEL FOR REGIME ESTIMATION
  #     NUMOFSTATES: NUMBER OF HIDDEN STATES IN REGIME
  # 
  # RETURNS: HMMFIT FROM HMMFIT(), SUITABLE FOR DISPLAY WITH PLOTMARKOVREGIME()
  
  #spx <- dROC(getOhlcv(instrumentSymbol="^GSPC", startDate=startDate, 
  #                     endDate=endDate, quote=c("close")))
  spx <- dROC(getFREDData("SP500",startDate=startDate,endDate=endDate)) # changed to use FRED
  
  spx[is.na(spx)] <- 0
  spxTurb <- rollingTurbulence(spx, avgWidth=(250 * 10), covarWidth=(250 * 10))
  meanTurb <- apply.monthly(spxTurb, mean)
  estimateMarkovRegimes(na.omit(meanTurb), numOfStates=numOfStates) # added omit
}

getInflationRegime <- function(startDate=as.Date("1946-01-01"), endDate=Sys.Date(),
                               numOfStates=2)
{
  # ESTIMATE TWO-STATE MARKOV (CPI-BASED) INFLATION REGIME.
  #
  # ARGS:
  #     STARTDATE: DATE WHICH TO BEGIN PANEL FOR REGIME ESTIMATION
  #     ENDDATE: END WHICH TO END PANEL FOR REGIME ESTIMATION
  #     NUMOFSTATES: NUMBER OF HIDDEN STATES IN REGIME
  # 
  # RETURNS: HMMFIT FROM HMMFIT(), SUITABLE FOR DISPLAY WITH PLOTMARKOVREGIME()
  
  val <- 100 *dROC(getFREDData(symbol="CPIAUCNS", startDate=startDate, 
                               endDate=endDate))
  estimateMarkovRegimes(val, numOfStates=numOfStates)
}

getGrowthRegime <- function(startDate=as.Date("1946-01-01"), 
                            endDate=as.Date("2012-12-31"), numOfStates=2)
{
  # ESTIMATE TWO-STATE MARKOV (GDP-BASED) GROWTH REGIME.
  #
  # NOTE: GROWTH REGIME APPEARS TO BE BI-MODAL, AND THUS NEED TO ESTIMATE
  # SEVERAL TIMES TO GET CONVERGENCE ON THE REGIME REPORTED BY KRITZMAN.
  #
  # ARGS:
  #     STARTDATE: DATE WHICH TO BEGIN PANEL FOR REGIME ESTIMATION
  #     ENDDATE: END WHICH TO END PANEL FOR REGIME ESTIMATION
  #     NUMOFSTATES: NUMBER OF HIDDEN STATES IN REGIME
  # 
  # RETURNS: HMMFIT FROM HMMFIT(), SUITABLE FOR DISPLAY WITH PLOTMARKOVREGIME()
  
  val <- 100 * dROC(getFREDData(symbol="GDPC1", startDate=startDate, 
                                endDate=endDate))
  estimateMarkovRegimes(val, numOfStates=numOfStates)
}

getCurrencyTurbulenceRegime <- function(startDate=as.Date("1971-01-01"), 
                                        endDate=Sys.Date(), 
                                        numOfStates=2)
{
  # ESTIMATE TWO-STATE MARKOV (G10-BASED) CURRENCY TURBULENCE REGIME.
  #
  # ARGS:
  #     STARTDATE: DATE WHICH TO BEGIN PANEL FOR REGIME ESTIMATION
  #     ENDDATE: END WHICH TO END PANEL FOR REGIME ESTIMATION
  #     NUMOFSTATES: NUMBER OF HIDDEN STATES IN REGIME
  # 
  # RETURNS: HMMFIT FROM HMMFIT(), SUITABLE FOR DISPLAY WITH PLOTMARKOVREGIME()
  
  # g10rates <- getG10Currencies()
  g10rates <- get.G10()
  avgg10rates <- xts(100 * rowMeans(dROC(g10rates), na.rm=T), 
                     order.by=last(index(g10rates), -1))
  turbG10rates <- rollingTurbulence(avgg10rates, avgWidth=(250 * 3), 
                                    covarWidth=(250 * 3))
  meanTurbG10rates <- apply.monthly(turbG10rates, mean)
  
  estimateMarkovRegimes(na.omit(meanTurbG10rates), numOfStates=numOfStates)
}

estimateMarkovRegimes <- function(val, numOfStates=2)
{
  # ESTIMATE N-STATE HIDDEN MARKOV MODEL (HMM) FOR VAL.
  #
  # ARGS:
  #     VAL: SERIES
  #     NUMOFSTATES: NUMBER OF HIDDEN STATES IN HMM
  #
  # RETURNS: HMMFIT FROM HMMFIT(), SUITABLE FOR DISPLAY WITH PLOTMARKOVREGIME()
  
  hmmFit <- HMMFit(val, nStates=numOfStates)
  return (list(val=val, hmmFit=hmmFit))
}

plotLocalDecodings <- function(regimes, 
                               symbols, 
                               plotDateRange="1900::2012", 
                               regimeNums,
                               drawColors=1:9)
{
  # PLOT LOCAL DECODINGS FOR A LIST OF HMM REGIMES, OPTIONALLY OVER A SET
  # DATE RANGE.
  #
  # ARGS:
  #     REGIMES: LIST OF REGIMES, AS RETURNED BY ESTIMATEMARKOVREGIMES()
  #     SYMBOLS: LIST OF HUMAN-READABLE SYMBOLS FOR REGIMES
  #     PLOTDATERANGE: OPTION DATE OVER WHICH TO PLOT REGIME LOCAL DECODINGS
  #     REGIMENUMS: INDEX OF HMM REGIME, INTO REGIMES, TO PLOT
  
  oldpar <- par(mfrow=c(1,1))
  on.exit(par(oldpar))
  
  layout(c(1,2,3,4))
  
  # GENERATE MERGE OF LOCAL DECODINGS
  localList <- lapply(c(1:length(regimes)), function(i) {
    
    regime <- regimes[[i]]
    fb <- forwardBackward(regime$hmmFit, regime$val)
    eu <- exp(fb$Alpha + fb$Beta - fb$LL)
    local <- xts(eu[,regimeNums[i]], index(regime$val))[plotDateRange]

    plota(local, type='l', plotX=T, col=drawColors[i], main=symbols[i])
  })
}

plotMarkovRegimes <- function(regime, 
                              symbol, 
                              plotDateRange="1900::2012",
                              plotDensity=T, 
                              plotTimeSeries=T,
                              drawColors=1:9)
{
  # PLOT MARKOV REGIMES FROM HMM: KERNEL DENSITIES AND PER-REGIME LOCAL DECODINGS.
  #
  # ARGS:
  #     HMMFIT: FIT FOR HMM, AS GENERATED BY ESTIMATEMARKOVREGIMES()
  #     SYMBOL: HUMAN-READABLE DESCRIPTION OF SERIES WITH MARKOV REGIMES
  #     PLOTDATERANGE: CONTIGUOUS RANGE OF TIME WHICH TO PLOT

  
  val <- regime$val
  hmmFit <- regime$hmmFit
  
  # CALCULATE LOCAL DECODING
  fb <- forwardBackward(hmmFit, val)
  eu <- exp(fb$Alpha + fb$Beta - fb$LL)
  hmmMeans <- hmmFit$HMM$distribution$mean
  hmmSD <- sqrt(hmmFit$HMM$distribution$var)
  
  # PLOT KERNEL DENSITY WITH REGIME MEANS
  oldpar <- par(mfrow=c(1,1))
  on.exit(par(oldpar))
  
  if (plotDensity)
  {
    plot(density(val), main=paste("Density with Regime Means:", symbol))
    abline(v=mean(val), lty=2)
    
    sapply(c(1:length(hmmMeans)), function(i) {
      abline(v=hmmMeans[i], lty=2, col=drawColors[(i+1)])
      curve(dnorm(x, hmmMeans[i], hmmSD[i]), add=T, lty=3,
            col=drawColors[(i+1)])
    })
  }
  
  # PLOT TIME SERIES OF PERCENT CHANGE AND LOCAL DECODING FOR EACH REGIME
  if (plotTimeSeries)
  {
    merged <- merge(val, eu)
    layout(c(1:(1+ncol(eu))))
    
    plota(merged[,1][plotDateRange], type='l', paste("Regime:", symbol), plotX=F)
    sapply(c(1:length(hmmMeans)), function(i) {
      abline(h=hmmMeans[i], lty=2, col=drawColors[(i+1)])
    })
    plota.legend("Percent Change:", drawColors[1], last(merged[,1]))
    
    sapply(c(1:ncol(eu)), function(i) {
      
      plota(xts(merged[,(i+1)], index(val))[plotDateRange], type='l', 
            plotX=(i==(ncol(eu))),
            col=drawColors[(i+1)])
      plota.legend(paste0("Event Regime ", i, ":"), drawColors[(i+1)], 
                   last(merged[,(i+1)]))
    })
  }
}


dROC <- function(x, n=1)
{
  # RETURN DISCRETE RATE-OF-CHANGE (ROC) FOR A SERIES, WITHOUT PADDING
  ROC(x, n, type="discrete", na.pad=F)
}


#added 
rollingTurbulence <- function(ret,avgWidth=250,covarWidth=250) {
  turbulence = ret[,1] * NA
  nperiods = nrow(ret)
  look.back = avgWidth
  
  for( i in (look.back+1) : nperiods ) {
    temp = ret[(i - look.back + 1):(i-1), ]
    
    # measures turbulence for the current observation
    turbulence[i] = mahalanobis(ret[i,], 
                                colMeans(temp,na.rm=TRUE), 
                                cov(temp,use="pairwise.complete.obs"))
    
    if( i %% 200 == 0) cat(i, 'out of', nperiods, '\n')
  }
  
  return(turbulence)
}

# added
getFREDData <- function(symbol,startDate,endDate) {
  dataEnv = new.env()
  getSymbols(symbol,src="FRED",env=dataEnv)
  s = get(symbol,env=dataEnv)
  s = s[paste(startDate,endDate,'/')]
  return(s)
}

# get G10 from Timely Portfolio
get.G10 <- function() {
  # FRED acronyms for daily FX rates
  map = '
  FX          FX.NAME        
  DEXUSAL     U.S./Australia 
  DEXUSUK     U.S./U.K.      
  DEXCAUS     Canada/U.S.    
  DEXNOUS     Norway/U.S.    
  DEXUSEU     U.S./Euro      
  DEXJPUS     Japan/U.S.     
  DEXUSNZ     U.S./NewZealand
  DEXSDUS     Sweden/U.S.    
  DEXSZUS     Switzerland/U.S.
  '
  
  map = matrix(scan(text = map, what='', quiet=T), nc=2, byrow=T)
  colnames(map) = map[1,]
  map = data.frame(map[-1,], stringsAsFactors=F)
  
  # convert all quotes to be vs U.S.
  convert.index = grep('DEXUS',map$FX, value=T)   
  
  #*****************************************************************
  # Load historical data
  #****************************************************************** 
  # load.packages('quantmod')
  
  # load fx from fred
  data.fx <- new.env()
  getSymbols(map$FX, src = 'FRED', from = '1970-01-01', env = data.fx, auto.assign = T)       
  for(i in convert.index) data.fx[[i]] = 1 / data.fx[[i]]
  
  # extract fx where all currencies are available
  # bt.prep(data.fx, align='remove.na')
  # fx = bt.apply(data.fx, '[')
  fx = get(map$FX[1],envir=data.fx)
  for ( i in 2:length(map$FX)) {
    f = get(map$FX[i],envir=data.fx)
    fx = cbind(fx,f)
  }
  fx = na.omit(fx)
  
  return(fx)
}

Pretty good!