Nothing
###############################################################################
# R (https://r-project.org/) Numeric Methods for Optimization of Portfolios
#
# Copyright (c) 2004-2021 Brian G. Peterson, Peter Carl, Ross Bennett, Kris Boudt
#
# This library is distributed under the terms of the GNU Public License (GPL)
# for full details see the file COPYING
#
# $Id: charts.DE.R 3378 2014-04-28 21:43:21Z rossbennett34 $
#
###############################################################################
# Note that many of these functions were provided by Kris Boudt and modified
# only slightly to work with this package
#' Statistical Factor Model
#'
#' Fit a statistical factor model using Principal Component Analysis (PCA)
#'
#' @details
#' The statistical factor model is fitted using \code{prcomp}. The factor
#' loadings, factor realizations, and residuals are computed and returned
#' given the number of factors used for the model.
#'
#' @param R xts of asset returns
#' @param k number of factors to use
#' @param \dots additional arguments passed to \code{prcomp}
#' @return
#' #' \describe{
#' \item{factor_loadings}{ N x k matrix of factor loadings (i.e. betas)}
#' \item{factor_realizations}{ m x k matrix of factor realizations}
#' \item{residuals}{ m x N matrix of model residuals representing idiosyncratic
#' risk factors}
#' }
#' Where N is the number of assets, k is the number of factors, and m is the
#' number of observations.
#' @export
statistical.factor.model <- function(R, k=1, ...){
if(!is.xts(R)){
R <- try(as.xts(R))
if(inherits(R, "try-error")) stop("R must be an xts object or coercible to an xts object")
}
# dimensions of R
m <- nrow(R)
N <- ncol(R)
# checks for R
if(m < N) stop("fewer observations than assets")
x <- coredata(R)
# Make sure k is an integer
if(k <= 0) stop("k must be a positive integer")
k <- as.integer(k)
# Fit a statistical factor model using Principal Component Analysis (PCA)
fit <- prcomp(x, ...=...)
# Extract the betas
# (N x k)
betas <- fit$rotation[, 1:k]
# Compute the estimated factor realizations
# (m x N) %*% (N x k) = (m x k)
f <- x %*% betas
# Compute the residuals
# These can be computed manually or by fitting a linear model
tmp <- x - f %*% t(betas)
b0 <- colMeans(tmp)
res <- tmp - matrix(rep(b0, m), ncol=N, byrow=TRUE)
# Compute residuals via fitting a linear model
# tmp.model <- lm(x ~ f)
# tmp.beta <- coef(tmp.model)[2:(k+1),]
# tmp.resid <- resid(tmp.model)
# all.equal(t(tmp.beta), betas, check.attributes=FALSE)
# all.equal(res, tmp.resid)
# structure and return
# stfm = *st*atistical *f*actor *m*odel
structure(list(factor_loadings=betas,
factor_realizations=f,
residuals=res,
m=m,
k=k,
N=N),
class="stfm")
}
#' Center
#'
#' Center a matrix
#'
#' This function is used primarily to center a time series of asset returns or
#' factors. Each column should represent the returns of an asset or factor
#' realizations. The expected value is taken as the sample mean.
#'
#' x.centered = x - mean(x)
#'
#' @param x matrix
#' @return matrix of centered data
#' @export
center <- function(x){
if(!is.matrix(x)) stop("x must be a matrix")
n <- nrow(x)
p <- ncol(x)
meanx <- colMeans(x)
x.centered <- x - matrix(rep(meanx, n), ncol=p, byrow=TRUE)
x.centered
}
##### Single Factor Model Comoments #####
#' Covariance Matrix Estimate
#'
#' Estimate covariance matrix using a single factor statistical factor model
#'
#' @details
#' This function estimates an (N x N) covariance matrix from a single factor
#' statistical factor model with k=1 factors, where N is the number of assets.
#'
#' @param beta vector of length N or (N x 1) matrix of factor loadings
#' (i.e. the betas) from a single factor statistical factor model
#' @param stockM2 vector of length N of the variance (2nd moment) of the
#' model residuals (i.e. idiosyncratic variance of the stock)
#' @param factorM2 scalar value of the 2nd moment of the factor realizations
#' from a single factor statistical factor model
#' @return (N x N) covariance matrix
covarianceSF <- function(beta, stockM2, factorM2){
# Beta of the stock with the factor index
beta = as.numeric(beta)
N = length(beta)
# Idiosyncratic variance of the stock
stockM2 = as.numeric(stockM2)
if(length(stockM2) != N) stop("dimensions do not match for beta and stockM2")
# Variance of the factor
factorM2 = as.numeric(factorM2)
# Coerce beta to a matrix
beta = matrix(beta, ncol = 1)
# Compute estimate
# S = (beta %*% t(beta)) * factorM2
S = tcrossprod(beta) * factorM2
D = diag(stockM2)
return(S + D)
}
#' Coskewness Matrix Estimate
#'
#' Estimate coskewness matrix using a single factor statistical factor model
#'
#' @details
#' This function estimates an (N x N^2) coskewness matrix from a single factor
#' statistical factor model with k=1 factors, where N is the number of assets.
#'
#' @param beta vector of length N or (N x 1) matrix of factor loadings
#' (i.e. the betas) from a single factor statistical factor model
#' @param stockM3 vector of length N of the 3rd moment of the model residuals
#' @param factorM3 scalar of the 3rd moment of the factor realizations from a
#' single factor statistical factor model
#' @return (N x N^2) coskewness matrix
coskewnessSF <- function(beta, stockM3, factorM3){
# Beta of the stock with the factor index
beta = as.numeric(beta)
N = length(beta)
# Idiosyncratic third moment of the stock
stockM3 = as.numeric(stockM3)
if(length(stockM3) != N) stop("dimensions do not match for beta and stockM3")
# Third moment of the factor
factorM3 = as.numeric(factorM3)
# Coerce beta to a matrix
beta = matrix(beta, ncol = 1)
# Compute estimate
# S = ((beta %*% t(beta)) %x% t(beta)) * factorM3
S = (tcrossprod(beta) %x% t(beta)) * factorM3
D = matrix(0, nrow=N, ncol=N^2)
for(i in 1:N){
col = (i - 1) * N + i
D[i, col] = stockM3[i]
}
return(S + D)
}
#' Cokurtosis Matrix Estimate
#'
#' Estimate cokurtosis matrix using a single factor statistical factor model
#'
#' @details
#' This function estimates an (N x N^3) cokurtosis matrix from a statistical
#' factor model with k factors, where N is the number of assets.
#'
#' @param beta vector of length N or (N x 1) matrix of factor loadings
#' (i.e. the betas) from a single factor statistical factor model
#' @param stockM2 vector of length N of the 2nd moment of the model residuals
#' @param stockM4 vector of length N of the 4th moment of the model residuals
#' @param factorM2 scalar of the 2nd moment of the factor realizations from a
#' single factor statistical factor model
#' @param factorM4 scalar of the 4th moment of the factor realizations from a
#' single factor statistical factor model
#' @return (N x N^3) cokurtosis matrix
cokurtosisSF <- function(beta, stockM2, stockM4, factorM2, factorM4){
# Beta of the stock with the factor index
beta = as.numeric(beta)
N = as.integer(length(beta))
# Idiosyncratic second moment of the stock
stockM2 = as.numeric(stockM2)
if(length(stockM2) != N) stop("dimensions do not match for beta and stockM2")
# Idiosyncratic fourth moment of the stock
stockM4 = as.numeric(stockM4)
if(length(stockM4) != N) stop("dimensions do not match for beta and stockM4")
# Second moment of the factor
factorM2 = as.numeric(factorM2)
# Fourth moment of the factor
factorM4 = as.numeric(factorM4)
# Compute estimate
# S = ((beta %*% t(beta)) %x% t(beta) %x% t(beta)) * factorM4
S = (tcrossprod(beta) %x% t(beta) %x% t(beta)) * factorM4
D = .residualcokurtosisSF(NN=N, sstockM2=stockM2, sstockM4=stockM4, mfactorM2=factorM2, bbeta=beta)
return(S + D)
}
# Wrapper function to compute the residual cokurtosis matrix of a statistical
# factor model with k = 1.
# Note that this function was orignally written in C++ (using Rcpp) by
# Joshua Ulrich and re-written using the C API by Ross Bennett
#' @useDynLib "PortfolioAnalytics"
.residualcokurtosisSF <- function(NN, sstockM2, sstockM4, mfactorM2, bbeta){
# NN : integer
# sstockM2 : vector of length NN
# sstockM4 : vector of length NN
# mfactorM2 : double
# bbeta : vector of length NN
if(!is.integer(NN)) NN <- as.integer(NN)
if(length(sstockM2) != NN) stop("sstockM2 must be a vector of length NN")
if(length(sstockM4) != NN) stop("sstockM4 must be a vector of length NN")
if(!is.double(mfactorM2)) mfactorM2 <- as.double(mfactorM2)
if(length(bbeta) != NN) stop("bbeta must be a vector of length NN")
.Call('residualcokurtosisSF', NN, sstockM2, sstockM4, mfactorM2, bbeta, PACKAGE="PortfolioAnalytics")
}
##### Multiple Factor Model Comoments #####
#' Covariance Matrix Estimate
#'
#' Estimate covariance matrix using a statistical factor model
#'
#' @details
#' This function estimates an (N x N) covariance matrix from a statistical
#' factor model with k factors, where N is the number of assets.
#'
#' @param beta (N x k) matrix of factor loadings (i.e. the betas) from a
#' statistical factor model
#' @param stockM2 vector of length N of the variance (2nd moment) of the
#' model residuals (i.e. idiosyncratic variance of the stock)
#' @param factorM2 (k x k) matrix of the covariance (2nd moment) of the factor
#' realizations from a statistical factor model
#' @return (N x N) covariance matrix
covarianceMF <- function(beta, stockM2, factorM2){
# Formula for covariance matrix estimate
# Sigma = beta %*% factorM2 %*% beta**T + Delta
# Delta is a diagonal matrix with the 2nd moment of residuals on the diagonal
# N = number of assets
# k = number of factors
# Use the dimensions of beta for checks of stockM2 and factorM2
# beta should be an (N x k) matrix
if(!is.matrix(beta)) stop("beta must be a matrix")
N <- nrow(beta)
k <- ncol(beta)
# stockM2 should be a vector of length N
stockM2 <- as.numeric(stockM2)
if(length(stockM2) != N) stop("dimensions do not match for beta and stockM2")
# factorM2 should be a (k x k) matrix
if(!is.matrix(factorM2)) stop("factorM2 must be a matrix")
if((nrow(factorM2) != k) | (ncol(factorM2) != k)){
stop("dimensions do not match for beta and factorM2")
}
# Compute covariance matrix
S <- beta %*% tcrossprod(factorM2, beta)
D <- diag(stockM2)
return(S + D)
}
#' Coskewness Matrix Estimate
#'
#' Estimate coskewness matrix using a statistical factor model
#'
#' @details
#' This function estimates an (N x N^2) coskewness matrix from a statistical
#' factor model with k factors, where N is the number of assets.
#'
#' @param beta (N x k) matrix of factor loadings (i.e. the betas) from a
#' statistical factor model
#' @param stockM3 vector of length N of the 3rd moment of the model residuals
#' @param factorM3 (k x k^2) matrix of the 3rd moment of the factor
#' realizations from a statistical factor model
#' @return (N x N^2) coskewness matrix
coskewnessMF <- function(beta, stockM3, factorM3){
# Formula for coskewness matrix estimate
# Phi = beta %*% factorM3 %*% (beta**T %x% beta**T) + Omega
# %x% is the kronecker product
# Omega is the (N x N^2) matrix matrix of zeros except for the i,j elements
# where j = (i - 1) * N + i, which is corresponding to the expected third
# moment of the idiosyncratic factors
# N = number of assets
# k = number of factors
# Use the dimensions of beta for checks of stockM2 and factorM2
# beta should be an (N x k) matrix
if(!is.matrix(beta)) stop("beta must be a matrix")
N <- nrow(beta)
k <- ncol(beta)
# stockM3 should be a vector of length N
stockM3 <- as.numeric(stockM3)
if(length(stockM3) != N) stop("dimensions do not match for beta and stockM3")
# factorM3 should be an (k x k^2) matrix
if(!is.matrix(factorM3)) stop("factorM3 must be a matrix")
if((nrow(factorM3) != k) | (ncol(factorM3) != k^2)){
stop("dimensions do not match for beta and factorM3")
}
# Compute coskewness matrix
beta.t <- t(beta)
S <- (beta %*% factorM3) %*% (beta.t %x% beta.t)
D <- matrix(0, nrow=N, ncol=N^2)
for(i in 1:N){
col <- (i - 1) * N + i
D[i, col] <- stockM3[i]
}
return(S + D)
}
#' Cokurtosis Matrix Estimate
#'
#' Estimate cokurtosis matrix using a statistical factor model
#'
#' @details
#' This function estimates an (N x N^3) cokurtosis matrix from a statistical
#' factor model with k factors, where N is the number of assets.
#'
#' @param beta (N x k) matrix of factor loadings (i.e. the betas) from a
#' statistical factor model
#' @param stockM2 vector of length N of the 2nd moment of the model residuals
#' @param stockM4 vector of length N of the 4th moment of the model residuals
#' @param factorM2 (k x k) matrix of the 2nd moment of the factor
#' realizations from a statistical factor model
#' @param factorM4 (k x k^3) matrix of the 4th moment of the factor
#' realizations from a statistical factor model
#' @return (N x N^3) cokurtosis matrix
cokurtosisMF <- function(beta , stockM2 , stockM4 , factorM2 , factorM4){
# Formula for cokurtosis matrix estimate
# Psi = beta %*% factorM4 %*% (beta**T %x% beta**T %x% beta**T) + Y
# %x% is the kronecker product
# Y is the residual matrix.
# see Asset allocation with higher order moments and factor models for
# definition of Y
# N = number of assets
# k = number of factors
# Use the dimensions of beta for checks of stockM2 and factorM2
# beta should be an (N x k) matrix
if(!is.matrix(beta)) stop("beta must be a matrix")
N <- nrow(beta)
k <- ncol(beta)
# stockM2 should be a vector of length N
stockM2 <- as.numeric(stockM2)
if(length(stockM2) != N) stop("dimensions do not match for beta and stockM2")
# stockM4 should be a vector of length N
stockM4 <- as.numeric(stockM4)
if(length(stockM4) != N) stop("dimensions do not match for beta and stockM4")
# factorM2 should be a (k x k) matrix
if(!is.matrix(factorM2)) stop("factorM2 must be a matrix")
if((nrow(factorM2) != k) | (ncol(factorM2) != k)){
stop("dimensions do not match for beta and factorM2")
}
# factorM4 should be a (k x k^3) matrix
if(!is.matrix(factorM4)) stop("factorM4 must be a matrix")
if((nrow(factorM4) != k) | (ncol(factorM4) != k^3)){
stop("dimensions do not match for beta and factorM4")
}
# Compute cokurtosis matrix
beta.t <- t(beta)
S <- (beta %*% factorM4) %*% (beta.t %x% beta.t %x% beta.t)
# betacov
betacov <- as.numeric(beta %*% tcrossprod(factorM2, beta))
# Compute the residual cokurtosis matrix
D <- .residualcokurtosisMF(NN=N, sstockM2=stockM2, sstockM4=stockM4, bbetacov=betacov)
return( S + D )
}
# Wrapper function to compute the residual cokurtosis matrix of a statistical
# factor model with k > 1.
# Note that this function was orignally written in C++ (using Rcpp) by
# Joshua Ulrich and re-written using the C API by Ross Bennett
#' @useDynLib "PortfolioAnalytics"
.residualcokurtosisMF <- function(NN, sstockM2, sstockM4, bbetacov){
# NN : integer, number of assets
# sstockM2 : numeric vector of length NN
# ssstockM4 : numeric vector of length NN
# bbetacov : numeric vector of length NN * NN
if(!is.integer(NN)) NN <- as.integer(NN)
if(length(sstockM2) != NN) stop("sstockM2 must be a vector of length NN")
if(length(sstockM4) != NN) stop("sstockM4 must be a vector of length NN")
if(length(bbetacov) != NN*NN) stop("bbetacov must be a vector of length NN*NN")
.Call('residualcokurtosisMF', NN, sstockM2, sstockM4, bbetacov, PACKAGE="PortfolioAnalytics")
}
##### Extract Moments #####
#' Covariance Estimate
#'
#' Extract the covariance matrix estimate from a statistical factor model
#'
#' @param model statistical factor model estimated via
#' \code{\link{statistical.factor.model}}
#' @param \dots not currently used
#' @return covariance matrix estimate
#' @seealso \code{\link{statistical.factor.model}}
#' @author Ross Bennett
#' @export
extractCovariance <- function(model, ...){
if(!inherits(model, "stfm")) stop("model must be of class 'stfm'")
# Extract elements from the model
beta <- model$factor_loadings
f <- model$factor_realizations
res <- model$residuals
m <- model$m
k <- model$k
N <- model$N
# Residual moments
denom <- m - k - 1
stockM2 <- colSums(res^2) / denom
# Factor moments
factorM2 <- cov(f)
# Compute covariance estimate
if(k == 1){
out <- covarianceSF(beta, stockM2, factorM2)
} else if(k > 1){
out <- covarianceMF(beta, stockM2, factorM2)
} else {
# invalid k
message("invalid k, returning NULL")
out <- NULL
}
return(out)
}
#' Coskewness Estimate
#'
#' Extract the coskewness matrix estimate from a statistical factor model
#'
#' @param model statistical factor model estimated via
#' \code{\link{statistical.factor.model}}
#' @param \dots not currently used
#' @return coskewness matrix estimate
#' @seealso \code{\link{statistical.factor.model}}
#' @author Ross Bennett
#' @export
extractCoskewness <- function(model, ...){
if(!inherits(model, "stfm")) stop("model must be of class 'stfm'")
# Extract elements from the model
beta <- model$factor_loadings
f <- model$factor_realizations
res <- model$residuals
m <- model$m
k <- model$k
N <- model$N
# Residual moments
denom <- m - k - 1
stockM3 <- colSums(res^3) / denom
# Factor moments
# f.centered <- center(f)
# factorM3 <- M3.MM(f.centered)
factorM3 <- PerformanceAnalytics::M3.MM(f)
# Compute covariance estimate
if(k == 1){
# Single factor model
out <- coskewnessSF(beta, stockM3, factorM3)
} else if(k > 1){
# Multi-factor model
out <- coskewnessMF(beta, stockM3, factorM3)
} else {
# invalid k
message("invalid k, returning NULL")
out <- NULL
}
return(out)
}
#' Cokurtosis Estimate
#'
#' Extract the cokurtosis matrix estimate from a statistical factor model
#'
#' @param model statistical factor model estimated via
#' \code{\link{statistical.factor.model}}
#' @param \dots not currently used
#' @return cokurtosis matrix estimate
#' @seealso \code{\link{statistical.factor.model}}
#' @author Ross Bennett
#' @export
extractCokurtosis <- function(model, ...){
if(!inherits(model, "stfm")) stop("model must be of class 'stfm'")
# Extract elements from the model
beta <- model$factor_loadings
f <- model$factor_realizations
res <- model$residuals
m <- model$m
k <- model$k
N <- model$N
# Residual moments
denom <- m - k - 1
stockM2 <- colSums(res^2) / denom
stockM4 <- colSums(res^4) / denom
# Factor moments
factorM2 <- cov(f)
# f.centered <- center(f)
# factorM4 <- M4.MM(f.centered)
factorM4 <- PerformanceAnalytics::M4.MM(f)
# Compute covariance estimate
if(k == 1){
# Single factor model
out <- cokurtosisSF(beta, stockM2, stockM4, factorM2, factorM4)
} else if(k > 1){
# Multi-factor model
out <- cokurtosisMF(beta, stockM2, stockM4, factorM2, factorM4)
} else {
# invalid k
message("invalid k, returning NULL")
out <- NULL
}
return(out)
}
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.