Nothing
# Copyright 2014-2020 Steven E. Pav. All Rights Reserved.
# Author: Steven E. Pav
#
# This file is part of MarkowitzR.
#
# MarkowitzR is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# MarkowitzR is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with MarkowitzR. If not, see <http://www.gnu.org/licenses/>.
# Created: 2014.01.31
# Copyright: Steven E. Pav, 2014
# Author: Steven E. Pav
# Comments: Steven E. Pav
#source("utils.r")
#' @title Estimate Markowitz Portfolio
#'
#' @description
#'
#' Estimates the Markowitz Portfolio or Markowitz Coefficient subject
#' to subspace and hedging constraints, and heteroskedasticity.
#'
#' @details
#'
#' Suppose that the expectation of \eqn{p}-vector \eqn{x} is linear
#' in the \eqn{f}-vector \eqn{f}, but the covariance of \eqn{x} is
#' stationary and independent of \eqn{f}. The 'Markowitz Coefficient'
#' is the \eqn{p \times f}{p x f} matrix \eqn{W} such that,
#' conditional on observing \eqn{f}, the portfolio \eqn{Wf} maximizes
#' Sharpe. When \eqn{f} is the constant 1, the Markowitz Coefficient
#' is the traditional Markowitz Portfolio.
#'
#' Given \eqn{n} observations of the returns and features, given
#' as matrices \eqn{X, F}, this code computes the Markowitz Coefficient
#' along with the variance-covariance matrix of the Coefficient and the
#' precision matrix. One may give optional weights, which are inverse
#' conditional volatility. One may also give optional matrix \eqn{J, G}
#' which define subspace and hedging constraints. Briefly, they constrain
#' the portfolio optimization problem to portfolios in the row space of
#' \eqn{J} and with zero covariance with the rows of \eqn{G}. It must
#' be the case that the rows of \eqn{J} span the rows of \eqn{G}.
#' \eqn{J} defaults to the \eqn{p \times p}{p x p} identity matrix,
#' and \eqn{G} defaults to a null matrix.
#'
#' One may use the default method for computing covariance,
#' via the \code{\link{vcov}} function, or via a 'fancy' estimator,
#' like \code{sandwich:vcovHAC}, \code{sandwich:vcovHC}, \emph{etc.}
#'
#' @usage
#'
#' mp_vcov(X,feat=NULL,vcov.func=vcov,fit.intercept=TRUE,weights=NULL,Jmat=NULL,Gmat=NULL)
#'
#' @param X an \eqn{n \times p}{n x p} matrix of observed returns.
#' @param feat an \eqn{n \times f}{n x f} matrix of observed features.
#' defaults to none, in which case \code{fit.intercept} must be
#' \code{TRUE}. If \code{fit.intercept} is true, ones will be prepended
#' to the features.
#' @param weights an optional \eqn{n} vector of the weights. The returns
#' and features will be multiplied by the weights. Weights should be
#' inverse volatility estimates. Defaults to homoskedasticity.
#' @param Jmat an optional \eqn{p_j \times p}{pj x p} matrix of the
#' subspace in which we constrain portfolios. Defaults essentially to the
#' \eqn{p \times p}{p x p} identity matrix.
#' @param Gmat an optional \eqn{p_g \times p}{pg x p} matrix of the
#' subspace to which we constrain portfolios to have zero covariance.
#' The rowspace of \code{Gmat} must be spanned by the rowspace of \code{Jmat}.
#' Defaults essentially to the \eqn{0 \times p}{0 x p} empty matrix.
#'
#' @inheritParams theta_vcov
#' @keywords univar
#'
#' @return a list containing the following components:
#' \item{mu}{Letting \eqn{r = f + p + fit.intercept}, this is a
#' \eqn{q = (r)(r+1)/2} vector...}
#' \item{Ohat}{The \eqn{q \times q}{q x q} estimated variance covariance
#' matrix of \code{mu}.}
#' \item{W}{The estimated Markowitz coefficient, a
#' \eqn{p \times (fit.intercept + f)}{p x (fit.intercept + f)} matrix. The
#' first column corresponds to the intercept term if it is fit. Note that
#' for convenience this function performs the sign flip, which is not performed
#' on \code{mu}.}
#' \item{What}{The estimated variance covariance matrix of \code{vech(W)}.
#' Letting \eqn{s = p(fit.intercept + f)}, this is a \eqn{s \times s}{s x s}
#' matrix.}
#' \item{widxs}{The indices into \code{mu} giving \code{W}, and into
#' \code{Ohat} giving \code{What}.}
#' \item{n}{The number of rows in \code{X}.}
#' \item{ff}{The number of features plus \code{as.numeric(fit.intercept)}.}
#' \item{p}{The number of assets.}
#'
#' @seealso \code{\link{theta_vcov}}, \code{\link{itheta_vcov}}
#' @rdname mp_vcov
#' @export
#' @template etc
#' @template ref-SEP13
#'
#' @note
#'
#' Should also modify to include the theta estimates.
#'
#' @examples
# 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))
#' # subspace constraint
#' Jmat <- matrix(rnorm(6),ncol=3)
#' ism <- mp_vcov(X,fit.intercept=TRUE,Jmat=Jmat)
#' walds <- ism$W / sqrt(diag(ism$What))
#' print(t(walds))
#' # hedging constraint
#' Gmat <- matrix(1,nrow=1,ncol=3)
#' ism <- mp_vcov(X,fit.intercept=TRUE,Gmat=Gmat)
#' walds <- ism$W / sqrt(diag(ism$What))
#'
#' # now 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 <- 2
#' n.ret <- 8
#' n.obs <- 10000
#' set.seed(101)
#' Feat <- matrix(rnorm(n.obs * n.feat),ncol=n.feat)
#' Wtrue <- 10 * 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)
#'
mp_vcov <- function(X,feat=NULL,vcov.func=vcov,fit.intercept=TRUE,
weights=NULL,Jmat=NULL,Gmat=NULL) {
set.coln(X)
if (is.null(vcov.func)) # belt and suspenders defaults
vcov.func=vcov
if (!is.null(feat)) {
set.coln(feat)
if (!is.null(weights)) {
# assume recycling will do this correctly!
feat <- feat * weights
}
f <- dim(feat)[2]
} else {
f <- 0
}
p <- dim(X)[2]
ff <- f + as.numeric(fit.intercept)
twidlize <- function(M) {
rbind(cbind(diag(ff),matrix(0,nrow=ff,ncol=dim(M)[2])),
cbind(matrix(0,nrow=dim(M)[1],ncol=ff),M))
}
if (!is.null(Jmat))
Jtwid <- twidlize(Jmat)
else
Jtwid <- NULL
# compute second moment and variance/covariance;
if (is.null(feat)) {
XY <- as.matrix(X)
} else {
XY <- cbind(as.matrix(feat),as.matrix(X))
}
XY <- na.omit(XY)
asymv <- theta_vcov(XY,fit.intercept=fit.intercept,
vcov.func=vcov.func)
# interpret
pp <- asymv$pp
Theta <- ivech(asymv$mu)
up.t <- .theta_grinder(Theta,M=Jtwid)
PTheta <- up.t$iT
preH <- up.t$iiT
if (!is.null(Gmat)) {
# 2FIX: check that G is spanned by J?
Gtwid <- twidlize(Gmat)
down.t <- .theta_grinder(Theta,M=Gtwid)
PTheta <- PTheta - down.t$iT
preH <- preH - down.t$iiT
}
mu <- matrixcalc::vech(PTheta)
elim.idx <- lower.tri(Theta,diag=TRUE)
dim(elim.idx) <- c(length(elim.idx),1)
Dupp <- matrixcalc::duplication.matrix(pp)
H.mat <- - preH[elim.idx,] %*% Dupp
Ohat <- (H.mat) %*% asymv$Ohat %*% t(H.mat)
# figure out indices of interest and then subselect!
botidx <- 0:(ff-1) * p + sapply(0:(ff-1),function(n)sum((ff-n):ff))
widxs <- outer(1:p,botidx,"+")
dim(widxs) <- c(length(widxs),1)
# sign flip to get Markowitz Coefficient ...
W <- - mu[widxs]
dim(W) <- c(p,ff)
rownames(W) <- colnames(X)
colnames(W)[(1:f) + as.numeric(fit.intercept)] <- colnames(feat)
if (fit.intercept)
colnames(W)[1] <- "Intercept"
What <- Ohat[widxs,widxs]
Wnames <- outer(rownames(W),colnames(W),FUN=function(x,y) { paste(x,"x",y) })
dim(Wnames) <- c(length(Wnames),1)
rownames(What) <- Wnames
colnames(What) <- Wnames
retval <- list(mu=mu,
Ohat=Ohat,
W=W,
What=What,
widxs=widxs,
n=asymv$n,
p=p,
ff=ff)
return(retval)
}
#for vim modeline: (do not edit)
# vim:fdm=marker:fmr=FOLDUP,UNFOLD:cms=#%s:syn=r:ft=r
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.