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 Compute variance covariance of 'Unified' Second Moment
#'
#' @description
#'
#' Computes the variance covariance matrix of sample mean and second moment.
#'
#' @details
#'
#' Given \eqn{p}-vector \eqn{x}, the 'unified' sample is the
#' \eqn{(p+1)(p+2)/2} vector of 1, \eqn{x}, and
#' \eqn{\mbox{vech}(x x^{\top})}{vech(x x')} stacked on top
#' of each other.
#' Given \eqn{n} contemporaneous observations of \eqn{p}-vectors,
#' stacked as rows in the \eqn{n \times p}{n x p} matrix \eqn{X},
#' this function computes the mean and the variance-covariance
#' matrix of the 'unified' sample.
#'
#' 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
#'
#' theta_vcov(X,vcov.func=vcov,fit.intercept=TRUE)
#'
#' @param X an \eqn{n \times p}{n x p} matrix of observed returns.
#' @param vcov.func a function which takes an object of class \code{lm},
#' and computes a variance-covariance matrix. If equal to the string
#' \code{"normal"}, we assume multivariate normal returns.
#' @param fit.intercept a boolean controlling whether we add a column
#' of ones to the data, or fit the raw uncentered second moment.
#' For now, must be true when assuming normal returns.
#' @keywords univar
#'
#' @return a list containing the following components:
#' \item{mu}{a \eqn{q = (p+1)(p+2)/2} vector of 1, then the mean,
#' then the vech'd second moment of the sample data.}
#' \item{Ohat}{the \eqn{q \times q}{q x q} estimated variance covariance
#' matrix. When \code{fit.intercept} is true, the left column and top row
#' are all zeros.}
#' \item{n}{the number of rows in \code{X}.}
#' \item{pp}{the number of assets plus \code{as.numeric(fit.intercept)}.}
#' @seealso \code{\link{itheta_vcov}}.
#' @rdname theta_vcov
#' @export
#' @template etc
#' @template ref-SEP13
#'
#' @note
#'
#' Replaces similar functionality from SharpeR package, but with
#' modified API.
#'
#' @examples
#' X <- matrix(rnorm(1000*3),ncol=3)
#' Sigmas <- theta_vcov(X)
#' Sigmas.n <- theta_vcov(X,vcov.func="normal")
#' Sigmas.n <- theta_vcov(X,fit.intercept=FALSE)
#'
#' # make it fat tailed:
#' X <- matrix(rt(1000*3,df=5),ncol=3)
#' Sigmas <- theta_vcov(X)
#' \donttest{
#' if (require(sandwich)) {
#' Sigmas <- theta_vcov(X,vcov.func=vcovHC)
#' }
#' }
#' # add some autocorrelation to X
#' Xf <- filter(X,c(0.2),"recursive")
#' colnames(Xf) <- colnames(X)
#' Sigmas <- theta_vcov(Xf)
#' \donttest{
#' if (require(sandwich)) {
#' Sigmas <- theta_vcov(Xf,vcov.func=vcovHAC)
#' }
#' }
#'
theta_vcov <- function(X,vcov.func=vcov,fit.intercept=TRUE) {
set.coln(X)
X <- na.omit(X)
if (fit.intercept)
X <- cbind(1,X)
p <- dim(X)[2]
TF <- outer(1:p,1:p,FUN=">=")
Y <- X[,col(TF)[TF]] * X[,row(TF)[TF]]
Ynames <- paste(colnames(X)[col(TF)[TF]],"x",colnames(X)[row(TF)[TF]])
colnames(Y) <- Ynames
# model Y and estimate Sigma hat#FOLDUP
if (is.character(vcov.func)) {
n <- dim(X)[1]
if (vcov.func != "normal")
stop("only understand function or 'normal'")
if (!fit.intercept)
stop("unknown form for gaussian zero mean case?")
mu <- colMeans(Y)
rm(Y)
Theta <- ivech(mu)
iTheta <- solve(Theta)
Dupp <- matrixcalc::duplication.matrix(p)
Elim <- matrixcalc::elimination.matrix(p)
# 2FIX there's a faster way to perform elimination, BTW!
LXD <- Elim %*% (iTheta %x% iTheta) %*% Dupp
DLXDU <- Dupp %*% LXD[,2:dim(LXD)[2]]
inner <- (Theta %x% Theta)
H <- t(DLXDU) %*% inner %*% DLXDU
Ohat <- (2 / n) * solve(H)
# now gate with zeros:
Ohat <- cbind(0,rbind(0,Ohat))
} else {
# take first moment of Y, then get standard error
mod2 <- lm(Y ~ 1)
rm(Y)
mu <- mod2$coefficients
Ohat <- vcov.func(mod2)
rm(mod2)
}#UNFOLD
# make nice output
dim(mu) <- c(length(mu),1)
rownames(mu) <- Ynames
dimnames(Ohat) <- list(Ynames,Ynames)
retval <- list(mu=mu,Ohat=Ohat,n=dim(X)[1],pp=p)
return(retval)
}
#' @title Compute variance covariance of Inverse 'Unified' Second Moment
#'
#' @description
#'
#' Computes the variance covariance matrix of the inverse unified
#' second moment matrix.
#'
#' @details
#'
#' Given \eqn{p}-vector \eqn{x} with mean \eqn{\mu}{mu} and
#' covariance, \eqn{\Sigma}{Sigma}, let \eqn{y} be \eqn{x}
#' with a one prepended. Then let
#' \eqn{\Theta = E\left(y y^{\top}\right)}{Theta = E[yy']},
#' the uncentered second moment matrix. The inverse of
#' \eqn{\Theta}{Theta} contains the (negative) Markowitz portfolio
#' and the precision matrix.
#'
#' Given \eqn{n} contemporaneous observations of \eqn{p}-vectors,
#' stacked as rows in the \eqn{n \times p}{n x p} matrix \eqn{X},
#' this function estimates the mean and the asymptotic
#' variance-covariance matrix of \eqn{\Theta^{-1}}{Theta^-1}.
#'
#' 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
#'
#' itheta_vcov(X,vcov.func=vcov,fit.intercept=TRUE)
#'
#'
#' @inheritParams theta_vcov
#' @return a list containing the following components:
#' \item{mu}{a \eqn{q = (p+1)(p+2)/2} vector of 1 + squared maximum
#' Sharpe, the negative Markowitz
#' portfolio, then the vech'd precision matrix of the sample data}
#' \item{Ohat}{the \eqn{q \times q}{q x q} estimated variance
#' covariance matrix.}
#' \item{n}{the number of rows in \code{X}.}
#' \item{pp}{the number of assets plus \code{as.numeric(fit.intercept)}.}
#' @keywords univar
#' @seealso \code{\link{theta_vcov}}, \code{\link{mp_vcov}}.
#' @rdname itheta_vcov
#' @export
#' @template etc
#' @template ref-SEP13
#'
#' @note
#'
#' By flipping the sign of \eqn{X}, the inverse of
#' \eqn{\Theta}{Theta} contains the \emph{positive} Markowitz
#' portfolio and the precision matrix on \eqn{X}. Performing
#' this transform before passing the data to this function
#' should be considered idiomatic.
#'
#' @note
#'
#' A more general form of this function exists as \code{\link{mp_vcov}}.
#'
#' @note
#'
#' Replaces similar functionality from SharpeR package, but with
#' modified API.
#'
#' @examples
#' X <- matrix(rnorm(1000*3),ncol=3)
#' # putting in -X is idiomatic:
#' ism <- itheta_vcov(-X)
#' iSigmas.n <- itheta_vcov(-X,vcov.func="normal")
#' iSigmas.n <- itheta_vcov(-X,fit.intercept=FALSE)
#' # compute the marginal Wald test statistics:
#' qidx <- 2:ism$pp
#' wald.stats <- ism$mu[qidx] / sqrt(diag(ism$Ohat[qidx,qidx]))
#'
#' # make it fat tailed:
#' X <- matrix(rt(1000*3,df=5),ncol=3)
#' ism <- itheta_vcov(X)
#' qidx <- 2:ism$pp
#' wald.stats <- ism$mu[qidx] / sqrt(diag(ism$Ohat[qidx,qidx]))
#' \donttest{
#' if (require(sandwich)) {
#' ism <- itheta_vcov(X,vcov.func=vcovHC)
#' qidx <- 2:ism$pp
#' wald.stats <- ism$mu[qidx] / sqrt(diag(ism$Ohat[qidx,qidx]))
#' }
#' }
#' # add some autocorrelation to X
#' Xf <- filter(X,c(0.2),"recursive")
#' colnames(Xf) <- colnames(X)
#' ism <- itheta_vcov(Xf)
#' qidx <- 2:ism$pp
#' wald.stats <- ism$mu[qidx] / sqrt(diag(ism$Ohat[qidx,qidx]))
#' \donttest{
#' if (require(sandwich)) {
#' ism <- itheta_vcov(Xf,vcov.func=vcovHAC)
#' qidx <- 2:ism$pp
#' wald.stats <- ism$mu[qidx] / sqrt(diag(ism$Ohat[qidx,qidx]))
#' }
#' }
itheta_vcov <- function(X,vcov.func=vcov,fit.intercept=TRUE) {
# delegate
sm_est <- theta_vcov(X,vcov.func,fit.intercept)
# interpret
pp <- sm_est$pp
Theta <- ivech(sm_est$mu)
iTheta <- solve(Theta)
iiTheta <- iTheta %x% iTheta # kron
# elimination & duplication
elim.idx <- lower.tri(Theta,diag=TRUE)
dim(elim.idx) <- c(length(elim.idx),1)
Dupp <- matrixcalc::duplication.matrix(pp)
# delta method outer deriv
H.mat <- iiTheta[elim.idx,] %*% Dupp
Ohat <- (H.mat) %*% sm_est$Ohat %*% t(H.mat)
# make nice output
mu <- iTheta[elim.idx]
Ynames <- rownames(sm_est$mu)
dim(mu) <- c(length(mu),1)
rownames(mu) <- Ynames
dimnames(Ohat) <- list(Ynames,Ynames)
retval <- list(mu=mu,Ohat=Ohat,n=sm_est$n,pp=pp)
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.