### Archive

Archive for July, 2012

## Yet Another Forecast Dashboard

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

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) != 'try-error') {
out = data.frame(out)
} else {
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)])),
}

plota.lines(out[,1], col='red')

labels = c('Data,Forecast', paste(gsub('Lo.', '', colnames(out)[2*(1:n)]), '%', sep=''))
}
```

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

```	#*****************************************************************
# Create models
#******************************************************************

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.

Categories: R

## 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)

#*****************************************************************
#******************************************************************
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)

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,] = NA
models\$spy = bt.run(data.spy)

# Equal Weight
data\$weight[] = NA
data\$weight[period.ends,] = ntop(prices, n)
data\$weight[1:period.ends,] = 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
#******************************************************************

# 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)

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')

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]

}

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\$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)

#*****************************************************************
#******************************************************************
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)

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,] = NA
models\$spy = bt.run(data.spy)

# Equal Weight
data\$weight[] = NA
data\$weight[period.ends,] = ntop(prices, n)
data\$weight[1:period.ends,] = 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.

```	#*****************************************************************
#******************************************************************
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)

```  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.

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)

#*****************************************************************
#******************************************************************
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)

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'

period.ends = endpoints(data\$prices, periodicity)
period.ends = period.ends[period.ends > 0]

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')

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.

Categories: Backtesting, Factors, R