In shabbychef/MarkowitzR: Statistical Significance of the Markowitz Portfolio

# set the knitr options ... for everyone!
# if you unset this, then vignette build bonks. oh, joy.
#opts_knit$set(progress=TRUE)
opts_knit$set(eval.after='fig.cap')
# for a package vignette, you do want to echo.
# opts_chunk$set(echo=FALSE,warning=FALSE,message=FALSE)
opts_chunk$set(warning=FALSE,message=FALSE)
#opts_chunk$set(results="asis")
opts_chunk$set(cache=TRUE,cache.path="cache/")

#opts_chunk$set(fig.path="github_extra/figure/",dev=c("pdf","cairo_ps"))
#opts_chunk$set(fig.path="github_extra/figure/",dev=c("png","pdf"))
#opts_chunk$set(fig.path="github_extra/figure/",dev=c("png"))
opts_chunk$set(fig.path="man/figures/",dev=c("png"))
opts_chunk$set(fig.width=5,fig.height=4,dpi=64)

# doing this means that png files are made of figures;
# the savings is small, and it looks like shit:
#opts_chunk$set(fig.path="figure/",dev=c("png","pdf","cairo_ps"))
#opts_chunk$set(fig.width=4,fig.height=4)
# for figures? this is sweave-specific?
#opts_knit$set(eps=TRUE)

# this would be for figures:
#opts_chunk$set(out.width='.8\\textwidth')
# for text wrapping:
options(width=64,digits=2)
opts_chunk$set(size="small")
opts_chunk$set(tidy=TRUE,tidy.opts=list(width.cutoff=50,keep.blank.line=TRUE))

#MarkowitzR.meta <- packageDescription('MarkowitzR')
library(MarkowitzR)

MarkowitzR

A number of utilities for dealing with the Markowitz portfolio.

-- Steven E. Pav, shabbychef@gmail.com

Installation

This package may be installed from CRAN; the latest version may be found on github via devtools:

if (require(devtools)) {
    # latest greatest
    install_github(repo='MarkowitzR',username='shabbychef',ref='master')
}

Basic Usage

Inference on the Markowitz Portfolio

The (negative) Markowitz portfolio appears in the inverse of the uncentered second moment matrix of the 'augmented' vector of returns. Via the Central Limit Theorem and the delta method the asymptotic distribution of the Markowitz portfolio can be found. From this, Wald statistics on the individual portfolio weights can be computed.

Fake Data

First for unconditional returns:

set.seed(1001)
X <- matrix(rnorm(1000*3),ncol=3)
ism <- mp_vcov(X,fit.intercept=TRUE)
walds <- ism$W / sqrt(diag(ism$What))
print(t(walds))

Now for conditional expectation:

# generate data with given W, Sigma
Xgen <- function(W,Sigma,Feat) {
 Btrue <- Sigma %*% W
 Xmean <- Feat %*% t(Btrue)
 Shalf <- chol(Sigma)
 X <- Xmean + matrix(rnorm(prod(dim(Xmean))),ncol=dim(Xmean)[2]) %*% Shalf
}

n.feat <- 3
n.ret <- 5 
n.obs <- 2000
set.seed(101)
Feat <- matrix(rnorm(n.obs * n.feat),ncol=n.feat)
Wtrue <- 5 * matrix(rnorm(n.feat * n.ret),ncol=n.feat)
Sigma <- cov(matrix(rnorm(100*n.ret),ncol=n.ret))
Sigma <- Sigma + diag(seq(from=1,to=3,length.out=n.ret))
X <- Xgen(Wtrue,Sigma,Feat)
ism <- mp_vcov(X,feat=Feat,fit.intercept=TRUE)

walds <- ism$W / sqrt(diag(ism$What))
print(t(walds))
# results are not much changed when using robust s.e.
library(sandwich)

ism.rse <- mp_vcov(X,feat=Feat,vcov.func=sandwich::vcovHAC,fit.intercept=TRUE)
walds.rse <- ism.rse$W / sqrt(diag(ism.rse$What))
print(t(walds.rse))

# errors should be asymptotically normal with the given covariance.
n.feat <- 5
n.ret <- 15 
n.obs <- 3000
set.seed(101)

Feat <- matrix(rnorm(n.obs * n.feat),ncol=n.feat)
Wtrue <- 5 * matrix(rnorm(n.feat * n.ret),ncol=n.feat)
Sigma <- cov(matrix(rnorm(100*n.ret),ncol=n.ret))
Sigma <- Sigma + diag(seq(from=1,to=3,length.out=n.ret))
X <- Xgen(Wtrue,Sigma,Feat)
ism <- mp_vcov(X,feat=Feat,fit.intercept=TRUE)

Wcomp <- cbind(0,Wtrue)
errs <- ism$W - Wcomp
dim(errs) <- c(length(errs),1)
Zerr <- solve(t(chol(ism$What)),errs)
print(summary(Zerr))
library(ggplot2)
ph <- ggplot(data.frame(Ze=Zerr),aes(sample=Ze)) + stat_qq() + geom_abline(slope=1,intercept=0,colour='red')
print(ph)
#qqnorm(Zerr)
#qqline(Zerr,col=2)

Fama French data

Now load the Fama French 3 factor portfolios.

if (!require(aqfb.data,quietly=TRUE) && require(devtools)) {
  # get the 10 industry data
  devtools::install_github('shabbychef/aqfb_data')
}
library(aqfb.data)
# fama 
data(mff4)

# will not matter, but convert pcts:
ff.data <- 1e-2 * mff4 

# risk free rate:
rfr <- ff.data[,'RF']

# subtract risk free from Mkt, HML and SMB:
ff.ret <- ff.data[,c('Mkt','HML','SMB')] - rep(rfr,2)

Now analyze the Markowitz portfolio on them.

ism <- mp_vcov(ff.ret,fit.intercept=TRUE)
walds <- ism$W / sqrt(diag(ism$What))
print(t(walds))
# now consider the hedging constraint: no covariance
# with the market:
Gmat <- matrix(c(1,0,0),nrow=1)
ism <- mp_vcov(ff.ret,fit.intercept=TRUE,Gmat=Gmat)
walds <- ism$W / sqrt(diag(ism$What))
print(t(walds))