Yet Another Forecast Dashboard

Recently, I came across quite a few examples of time series forecasting using R. Here are some examples:

• Time series cross-validation 4: forecasting the S&P 500
• Holt-Winters forecast using ggplot2
• Autoplot: Graphical Methods with ggplot2
• Large-Scale Parallel Statistical Forecasting Computations in R (2011) by M. Stokely, F. Rohani, E. Tassone
• Forecasting time series data
• ARIMA Sector Forecasts

So I have decided to roll my own version of Forecast Dashboard as well. First, let’s define some helper functions:

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

# extract forecast info
forecast.helper <- function(fit, h=10, level = c(80,95)) {
	out = try(forecast(fit, h=h, level=level), silent=TRUE)
	if (class(out)[1] != 'try-error') {
		out = data.frame(out)
	} else {
		temp = data.frame(predict(fit, n.ahead=h, doplot=F))
			pred = temp[,1]
			se = temp[,2]
		qq = qnorm(0.5 * (1 + level/100))
		out = matrix(NA, nr=h, nc=1+2*len(qq))
			out[,1] = pred
		for(i in 1:len(qq))
			out[,(2*i):(2*i+1)] = c(pred - qq[i] * se, pred + qq[i] * se)
		colnames(out) = c('Point.Forecast', matrix(c(paste('Lo', level, sep='.'), paste('Hi', level, sep='.')), nr=2, byrow=T))
		out = data.frame(out)
	}	
	return(out)
}

# compute future dates for the forecast
forecast2xts <- function(data, forecast) {
	# length of the forecast
	h = nrow(forecast)
	dates = as.Date(index(data))
 		
	new.dates = seq(last(dates)+1, last(dates) + 2*365, by='day')
	rm.index = date.dayofweek(new.dates) == 6 | date.dayofweek(new.dates) == 0
 	new.dates = new.dates[!rm.index]
 		
 	new.dates = new.dates[1:h] 	
 	return(make.xts(forecast, new.dates))
}

# create forecast plot
forecast.plot <- function(data, forecast, ...) {
	out = forecast2xts(data, forecast)

 	# create plot
 	plota(c(data, out[,1]*NA), type='l', 
 			ylim = range(data,out,na.rm=T), ...) 		
 	
 	# highligh sections
	new.dates = index(out)
		temp = coredata(out)

	n = (ncol(out) %/% 2)
	for(i in n : 1) {
		polygon(c(new.dates,rev(new.dates)), 
			c(temp[,(2*i)], rev(temp[,(2*i+1)])), 
		border=NA, col=col.add.alpha(i+2,150))
	}
	
	plota.lines(out[,1], col='red')
	
	labels = c('Data,Forecast', paste(gsub('Lo.', '', colnames(out)[2*(1:n)]), '%', sep=''))
	plota.legend(labels, fill = c('black,red',col.add.alpha((1:n)+2, 150)))			
}

Now we are ready to fit time series models and create a sample Forecast Dashboard. Below are some examples:

	#*****************************************************************
	# Create models
	#****************************************************************** 
	load.packages('forecast,fGarch,fArma')
	
	sample = last(data$prices$SPY, 200)	
	ts.sample = ts(sample, frequency = 12)
	
	models = list(
		# fGarch		
		garch = garchFit(~arma(1,15)+garch(1,1), data=sample, trace=F),
		# fArma
		arima = armaFit(~ arima(1, 1, 15), data=ts.sample),
				
		# forecast
		arma = Arima(ts.sample, c(1,0,1)),
		arfima = arfima(ts.sample),
		auto.arima = auto.arima(ts.sample),
	
		bats = bats(ts.sample),
		HoltWinters = HoltWinters(ts.sample),
		naive = Arima(ts.sample, c(0,1,0))
	)
	
	
	#*****************************************************************
	# Create Report
	#****************************************************************** 						
	# 30 day forecast with 80% and 95% confidence bands
		layout(matrix(1:9,nr=3))
		for(i in 1:len(models)) {
			out = forecast.helper(models[[i]], 30, level = c(80,95))	 	
			forecast.plot(sample, out, main = names(models)[i]) 	
		}	
	
	# 30 day forecast with 75%,85%,95%,97%, and 99% confidence bands
		layout(matrix(1:9,nr=3))
		for(i in 1:len(models)) {
			out = forecast.helper(models[[i]], 30, level = c(75,85,95,97,99))	 	
			forecast.plot(sample, out, main = names(models)[i]) 	
		}	

I encourage you to read more about various time series models available in R and share your examples of Forecast Dashboards.

To view the complete source code for this example, please have a look at the bt.forecast.dashboard() function in bt.test.r at github.

Advertisement

More on Factor Attribution to improve performance of the 1-Month Reversal Strategy

In my last post, Factor Attribution to improve performance of the 1-Month Reversal Strategy, I discussed how Factor Attribution can be used to boost performance of the 1-Month Reversal Strategy. Today I want to dig a little dipper and examine this strategy for each sector and also run a sector-neutral back-test.

The initial steps to load historical prices, load Fama/French factors and run Factor Attribution are the same as in the original post, so I will not repeat them here. The complete source code for this example is located in the bt.fa.sector.one.month.test() function in bt.test.r at github for your reference. The new functionality to create sector and sector-neutral back-tests is located in the bt.make.quintiles.sector() function below:

bt.make.quintiles.sector <- function(
	sector,		# sector data
	position.score,	# position.score is a factor to form Quintiles sampled at the period.ends
	data,		# back-test object
	period.ends,	
	n.quantiles = 5,
	start.t = 2,	# first index at which to form Quintiles
	prefix = ''	
) 
{
	#*****************************************************************
	# Re-organize sectors into matrix, assume that sectors are constant in time
	#****************************************************************** 
	temp = factor(sector)
	sector.names = levels(temp)	
		n.sectors = len(sector.names)
	sectors = matrix(unclass(temp), nr=nrow(position.score), nc=ncol(position.score), byrow=T)
	
	#*****************************************************************
	# Create Quantiles
	#****************************************************************** 
	position.score = coredata(position.score)
	quantiles = weights = position.score * NA			
	
	for( s in 1:n.sectors) {
		for( t in start.t:nrow(weights) ) {
			index = sectors[t,] == s
			n = sum(index)
			
			# require at least 3 companies in each quantile
			if(n > 3*n.quantiles) {			
				factor = as.vector(position.score[t, index])
				ranking = ceiling(n.quantiles * rank(factor, na.last = 'keep','first') / count(factor))
			
				quantiles[t, index] = ranking
				weights[t, index] = 1/tapply(rep(1,n), ranking, sum)[ranking]			
			}
		}
	}
	quantiles = ifna(quantiles,0)
	
	#*****************************************************************
	# Create Q1-QN spread for each Sector
	#****************************************************************** 
	long = weights * NA
	short = weights * NA
	models = list()
	
	for( s in 1:n.sectors) {
		long[] = 0
		long[quantiles == 1 & sectors == s] = weights[quantiles == 1 & sectors == s]
		long = long / rowSums(long,na.rm=T)
		
		short[] = 0
		short[quantiles == n.quantiles & sectors == s] = weights[quantiles == n.quantiles & sectors == s]
		short = short / rowSums(short,na.rm=T)
		
		data$weight[] = NA
			data$weight[period.ends,] = long - short
		models[[ paste(prefix,'spread.',sector.names[s], sep='') ]]	= bt.run(data, silent = T)	
	}

	#*****************************************************************
	# Create Sector - Neutral Q1-QN spread
	#****************************************************************** 		
	long[] = 0
	long[quantiles == 1] = weights[quantiles == 1]
	long = long / rowSums(long,na.rm=T)
	
	short[] = 0
	short[quantiles == n.quantiles] = weights[quantiles == n.quantiles]
	short = short / rowSums(short,na.rm=T)
		
	data$weight[] = NA
		data$weight[period.ends,] = long - short
	models$spread.sn = bt.run(data, silent = T)	

	return(models)
}

I have compared below, the overall and sector-neutral versions of the strategy based on the 1 Month returns (one.month) and the Residuals from the Factor Attribution regression (last.e_s). In all cases the Residual strategy outperformed the base one and also the sector-neutral versions have better risk-adjusted coefficients compared to the overall strategy.

Next I looked at the each sector performance for both strategies and surprisingly found that Energy, Materials, and Utilities underperformed in both cases and Financials and Consumer Staples did very well in both cases. Looking at the sector back-test charts, I think the Momentum strategy that selects top few sectors each month would do very well. You can try implementing this idea as a homework.

To view the complete source code for this example, please have a look at the bt.fa.sector.one.month.test() function in bt.test.r at github.

Factor Attribution to improve performance of the 1-Month Reversal Strategy

Today I want to show how to use Factor Attribution to boost performance of the 1-Month Reversal Strategy. The Short-Term Residual Reversal by D. Blitz, J. Huij, S. Lansdorp, M. Verbeek (2011) paper presents the idea and discusses the results as applied to US stock market since 1929. To improve 1-Month Reversal Strategy performance authors investigate an alternative position ranking metric based on the residuals from the rolling Factor Attribution regression.

The base 1-Month Reversal Strategy forms portfolios each month by buying 20% of loosers and short selling 20% of winners from the S&P 500 index based on the prior 1-Month returns. The 1-Month Residual Reversal Strategy forms portfolios each month by buying 20% of loosers and short selling 20% of winners from the S&P 500 index based on the residuals from the rolling Factor Attribution regression. Following are the steps to form 1-Month Residual Reversal Strategy portfolio each month:

• 1. for each company in the S&P 500 index, run a rolling Factor Attribution regression on the prior 36 months and compute residuals: e.1, e.2, …, e.t, …, e.T for t in 1:36
• 2. alternative position ranking metric = e.T / standard deviation of (e)
• 3. form portfolios by buying 20% of loosers and short selling 20% of winners from the S&P 500 index based on the alternative position ranking metric

Let’s start by loading historical prices for all companies in the S&P 500 and create SPY and Equal Weight benchmarks using the Systematic Investor Toolbox:

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = sp500.components()$tickers
	
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
		# remove companies with less than 5 years of data
		rm.index = which( sapply(ls(data), function(x) nrow(data[[x]])) < 1000 )	
		rm(list=names(rm.index), envir=data)
		
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)		
	bt.prep(data, align='keep.all', dates='1994::')
		tickers = data$symbolnames
	
	
	data.spy <- new.env()
	getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data.spy, auto.assign = T)
	bt.prep(data.spy, align='keep.all', dates='1994::')
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices
		n = ncol(prices)
					
	#*****************************************************************
	# Setup monthly periods
	#****************************************************************** 
	periodicity = 'months'
	
	period.ends = endpoints(data$prices, periodicity)
		period.ends = period.ends[period.ends > 0]
	
	prices = prices[period.ends, ]		
		
	#*****************************************************************
	# Create Benchmarks, omit results for the first 36 months - to be consistent with Factor Attribution
	#****************************************************************** 	
	models = list()
	
	# SPY
	data.spy$weight[] = NA
		data.spy$weight[] = 1
		data.spy$weight[1:period.ends[36],] = NA
	models$spy = bt.run(data.spy)
	
	# Equal Weight
	data$weight[] = NA
		data$weight[period.ends,] = ntop(prices, n)
		data$weight[1:period.ends[36],] = NA		
	models$equal.weight = bt.run(data)

Next let’s run Factor Attribution on each the stocks in the S&P 500 index to determine it’s alternative position ranking metric. I will save both e.T and e.T / standard deviation of (e) metrics.

# function to compute additional custom stats for factor.rolling.regression
factor.rolling.regression.custom.stats <- function(x,y,fit) {
	n = len(y)
	e = y - x %*% fit$coefficients
	se = sd(e)
	return(c(e[n], e[n]/se))
}

	#*****************************************************************
	# Load factors and align them with prices
	#****************************************************************** 	
	# load Fama/French factors
	factors = get.fama.french.data('F-F_Research_Data_Factors', periodicity = periodicity,download = T, clean = F)
	
	# align monthly dates
	map = match(format(index(factors$data), '%Y%m'), format(index(prices), '%Y%m'))
		dates = index(factors$data)
		dates[!is.na(map)] = index(prices)[na.omit(map)]
	index(factors$data) = as.Date(dates)
		
	
	# add factors and align
	data.fa <- new.env()
		for(i in tickers) data.fa[[i]] = data[[i]][period.ends, ]
		data.fa$factors = factors$data / 100
	bt.prep(data.fa, align='remove.na')

	
	index = match( index(data.fa$prices), index(data$prices) )
		prices = data$prices[index, ]
		
	#*****************************************************************
	# Compute Factor Attribution for each ticker
	#****************************************************************** 	

	temp = NA * prices
	factors	= list()
		factors$last.e = temp
		factors$last.e_s = temp
	
	for(i in tickers) {
		cat(i, '\n')
		
		# Facto Loadings Regression
		obj = factor.rolling.regression(data.fa, i, 36, silent=T,
			factor.rolling.regression.custom.stats)

		for(j in 1:len(factors))		
			factors[[j]][,i] = obj$fl$custom[,j]
			
	}

	# add base strategy
	factors$one.month = coredata(prices / mlag(prices))

Next let’s group stocks into Quantiles based on 1-Month Reversal factors and create reports.

	#*****************************************************************
	# Create Quantiles
	#****************************************************************** 
	quantiles = list()
	
	for(name in names(factors)) {
		cat(name, '\n')
		quantiles[[name]] = bt.make.quintiles(factors[[name]], data, index, start.t =  1+36, prefix=paste(name,'_',sep=''))
	}	

	#*****************************************************************
	# Create Report
	#****************************************************************** 					
	plotbt.custom.report.part1(quantiles$one.month$spread,quantiles$last.e$spread,quantiles$last.e_s$spread)
	
	plotbt.strategy.sidebyside(quantiles$one.month$spread,quantiles$last.e$spread,quantiles$last.e_s$spread)
		
	plotbt.custom.report.part1(quantiles$last.e_s)

The 1-Month Residual Reversal Strategy have done well over the last 10 years and handsomely outperformed the base 1-Month Reversal Strategy.

To view the complete source code for this example, please have a look at the bt.fa.one.month.test() function in bt.test.r at github.

1-Month Reversal Strategy

Today I want to show a simple example of the 1-Month Reversal Strategy. Each month we will buy 20% of loosers and short sell 20% of winners from the S&P 500 index. The loosers and winners are measured by prior 1-Month returns. I will use this post to set the stage for my next post that will show how Factor Attribution can boost performance of the 1-Month Reversal Strategy. Following is the references for my next post, in case you want to get a flavor, Short-Term Residual Reversal by D. Blitz, J. Huij, S. Lansdorp, M. Verbeek (2011) paper.

Let’s start by loading historical prices for all companies in the S&P 500 and create SPY and Equal Weight benchmarks using the Systematic Investor Toolbox:

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = sp500.components()$tickers
	
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
		# remove companies with less than 5 years of data
		rm.index = which( sapply(ls(data), function(x) nrow(data[[x]])) < 1000 )	
		rm(list=names(rm.index), envir=data)
		
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)		
	bt.prep(data, align='keep.all', dates='1994::')
		tickers = data$symbolnames
	
	
	data.spy <- new.env()
	getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data.spy, auto.assign = T)
	bt.prep(data.spy, align='keep.all', dates='1994::')
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices
		n = ncol(prices)
					
	#*****************************************************************
	# Setup monthly periods
	#****************************************************************** 
	periodicity = 'months'
	
	period.ends = endpoints(data$prices, periodicity)
		period.ends = period.ends[period.ends > 0]
	
	prices = prices[period.ends, ]		
		
	#*****************************************************************
	# Create Benchmarks, omit results for the first 36 months - to be consistent with Factor Attribution
	#****************************************************************** 	
	models = list()
	
	# SPY
	data.spy$weight[] = NA
		data.spy$weight[] = 1
		data.spy$weight[1:period.ends[36],] = NA
	models$spy = bt.run(data.spy)
	
	# Equal Weight
	data$weight[] = NA
		data$weight[period.ends,] = ntop(prices, n)
		data$weight[1:period.ends[36],] = NA		
	models$equal.weight = bt.run(data)

Next let’s group stocks into Quantiles based on 1-Month returns and create back-test for each Quantile. I will rely on the code in the Volatility Quantiles post to create Quantiles.

	#*****************************************************************
	# Create Reversal Quantiles
	#****************************************************************** 
	n.quantiles = 5
	start.t = 1 + 36
	quantiles = weights = coredata(prices) * NA			
	
	one.month = coredata(prices / mlag(prices))
	
	for( t in start.t:nrow(weights) ) {
		factor = as.vector(one.month[t,])
		ranking = ceiling(n.quantiles * rank(factor, na.last = 'keep','first') / count(factor))
		
		quantiles[t,] = ranking
		weights[t,] = 1/tapply(rep(1,n), ranking, sum)[ranking]			
	}
	
	quantiles = ifna(quantiles,0)
	
	#*****************************************************************
	# Create backtest for each Quintile
	#****************************************************************** 
	temp = weights * NA

	for( i in 1:n.quantiles) {
		temp[] = 0
		temp[quantiles == i] = weights[quantiles == i]
	
		data$weight[] = NA
			data$weight[period.ends,] = temp
		models[[ paste('M1_Q',i,sep='') ]] = bt.run(data, silent = T)
	}

Finally, let’s construct Q1/Q5 spread and create summary performance report.

	#*****************************************************************
	# Create Q1-Q5 spread
	#****************************************************************** 
	temp[] = 0
	temp[quantiles == 1] = weights[quantiles == 1]
	temp[quantiles == n.quantiles] = -weights[quantiles == n.quantiles]
	
	data$weight[] = NA
		data$weight[period.ends,] = temp
	models$spread = bt.run(data, silent = T)	

	#*****************************************************************
	# Create Report
	#****************************************************************** 	
	plotbt.custom.report.part1(models)

	plotbt.custom.report.part1(models[spl('spy,equal.weight,spread')])

In the next post I will show how Factor Attribution can boost performance of the 1-Month Reversal Strategy using the methodology presented in the Short-Term Residual Reversal by D. Blitz, J. Huij, S. Lansdorp, M. Verbeek (2011) paper.

To view the complete source code for this example, please have a look at the bt.one.month.test() function in bt.test.r at github.

Example of Factor Attribution

In the prior post, Factor Attribution 2, I have shown how Factor Attribution can be applied to decompose fund’s returns in to Market, Capitalization, and Value factors, the “three-factor model” of Fama and French. Today, I want to show you a different application of Factor Attribution. First, let’s run Factor Attribution on each the stocks in the S&P 500 to determine it’s Value exposure. Next let’s group stocks into Quantiles based on Value exposure and create back-test for each Quantile. I will rely on the code in the Volatility Quantiles post to create Quantiles.

Let’s start by loading historical prices for all current components of the S&P 500 index.

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = sp500.components()$tickers
	
	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
		# remove companies with less than 5 years of data
		rm.index = which( sapply(ls(data), function(x) nrow(data[[x]])) < 1000 )	
		rm(list=names(rm.index), envir=data)
		
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)		
	bt.prep(data, align='keep.all', dates='1994::')
		tickers = data$symbolnames
	
	
	data.spy <- new.env()
	getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data.spy, auto.assign = T)
	bt.prep(data.spy, align='keep.all', dates='1994::')
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 
	prices = data$prices
		nperiods = nrow(prices)
		n = ncol(prices)
			
	models = list()
	
	# SPY
	data.spy$weight[] = NA
		data.spy$weight[] = 1
	models$spy = bt.run(data.spy)
	
	# Equal Weight
	data$weight[] = NA
		data$weight[] = ntop(prices, n)
	models$equal.weight = bt.run(data)

Next let’s run Factor Attribution on each the stocks in the S&P 500 to determine it’s Value exposure.

	#*****************************************************************
	# Compute Factor Attribution for each ticker
	#****************************************************************** 
	periodicity = 'weeks'
	
	# load Fama/French factors
	factors = get.fama.french.data('F-F_Research_Data_Factors', periodicity = periodicity,download = T, clean = F)
	
	period.ends = endpoints(data$prices, periodicity)
		period.ends = period.ends[period.ends > 0]
	
	# add factors and align
	data.fa <- new.env()
		for(i in tickers) data.fa[[i]] = data[[i]][period.ends,]
	data.fa$factors = factors$data / 100
	bt.prep(data.fa, align='remove.na')

	
	index = match( index(data.fa$prices), index(data$prices) )
	measure = data$prices[ index, ]	
	for(i in tickers) {
		cat(i, '\n')
		
		# Facto Loadings Regression
		obj = factor.rolling.regression(data.fa, i, 36, silent=T)
		
		measure[,i] = coredata(obj$fl$estimate$HML)
	}

Finally, let’s group stocks into Quantiles based on Value exposure and create back-test for each Quantile.

	#*****************************************************************
	# Create Value Quantiles
	#****************************************************************** 
	n.quantiles=5
	start.t = 1+36
	quantiles = weights = coredata(measure) * NA			
	
	for( t in start.t:nrow(weights) ) {
		factor = as.vector(coredata(measure[t,]))
		ranking = ceiling(n.quantiles * rank(factor, na.last = 'keep','first') / count(factor))
		
		quantiles[t,] = ranking
		weights[t,] = 1/tapply(rep(1,n), ranking, sum)[ranking]			
	}

	quantiles = ifna(quantiles,0)
	
	#*****************************************************************
	# Create backtest for each Quintile
	#****************************************************************** 
	for( i in 1:n.quantiles) {
		temp = weights * NA
			temp[] = 0
		temp[quantiles == i] = weights[quantiles == i]
	
		data$weight[] = NA
			data$weight[index,] = temp
		models[[ paste('Q',i,sep='_') ]] = bt.run(data, silent = T)
	}
	
	#*****************************************************************
	# Create Report
	#****************************************************************** 					
	plotbt.custom.report.part1(models)		
	
	plotbt.strategy.sidebyside(models)

There is no linear relationship between Value Quantiles and historical performance. I’m also suspecting that that implied Value exposure might be quite different than the real Price/Book ratio for each stock. Let me know what do you think about this approach.

In the next post I will show another example of Factor Attribution.

To view the complete source code for this example, please have a look at the bt.fa.value.quantiles.test() function in bt.test.r at github.