1############################################################################### 2# Functions to perform multivariate matrix 3# calculations on portfolios of assets. 4# 5# I've modified these to minimize the number of 6# times the same calculation or statistic is run, and to minimize duplication 7# of code from function to function. This should make things more 8# efficient when running against very large numbers of instruments or portfolios. 9# 10# Copyright (c) 2008 Kris Boudt and Brian G. Peterson 11# Copyright (c) 2004-2020 Peter Carl and Brian G. Peterson for PerformanceAnalytics 12# This R package is distributed under the terms of the GNU Public License (GPL) 13# for full details see the file COPYING 14############################################################################### 15# $Id$ 16############################################################################### 17 18 19M3.MM.old = function(R,...){ 20 cAssets = ncol(R); T = nrow(R); 21 if(!hasArg(mu)) mu = apply(R,2,'mean') else mu=mu=list(...)$mu 22 M3 = matrix(rep(0,cAssets^3),nrow=cAssets,ncol=cAssets^2) 23 for(t in c(1:T)) 24 { 25 centret = as.numeric(matrix(R[t,]-mu,nrow=cAssets,ncol=1)) 26 M3 = M3 + ( centret%*%t(centret) )%x%t(centret) 27 } 28 return( 1/T*M3 ); 29} 30 31#'@useDynLib PerformanceAnalytics 32#'@export 33#'@rdname CoMoments 34M3.MM <- function(R, unbiased = FALSE, as.mat = TRUE, ...) { 35 if(!hasArg(mu)) mu = colMeans(R) else mu=list(...)$mu 36 37 x <- coredata(R) 38 39 NN <- NROW(x) 40 PP <- NCOL(x) 41 42 Xc <- x - matrix(mu, nrow = NN, ncol = PP, byrow = TRUE) 43 44 if (unbiased) { 45 if (NN < 3) stop("R should have at least 3 observations") 46 CC <- NN / ((NN - 1) * (NN - 2)) 47 } else { 48 CC <- 1 / NN 49 } 50 51 M3 <- .Call('M3sample', as.numeric(Xc), NN, PP, CC, PACKAGE="PerformanceAnalytics") 52 if (as.mat) M3 <- M3.vec2mat(M3, PP) 53 54 return( M3 ) 55} 56 57M4.MM.old = function(R,...){ 58 cAssets = ncol(R); T = nrow(R); 59 if(!hasArg(mu)) mu = apply(R,2,'mean') else mu=list(...)$mu 60 M4 = matrix(rep(0,cAssets^4),nrow=cAssets,ncol=cAssets^3); 61 for(t in c(1:T)) 62 { 63 centret = as.numeric(matrix(R[t,]-mu,nrow=cAssets,ncol=1)) 64 M4 = M4 + ( centret%*%t(centret) )%x%t(centret)%x%t(centret) 65 } 66 return( 1/T*M4 ); 67} 68 69#'@useDynLib PerformanceAnalytics 70#'@export 71#'@rdname CoMoments 72M4.MM <- function(R, as.mat = TRUE, ...) { 73 if(!hasArg(mu)) mu = colMeans(R) else mu=list(...)$mu 74 75 x <- coredata(R) 76 77 NN <- NROW(x) 78 PP <- NCOL(x) 79 80 Xc <- x - matrix(mu, nrow = NN, ncol = PP, byrow = TRUE) 81 82 M4 <- .Call('M4sample', as.numeric(Xc), NN, PP, PACKAGE="PerformanceAnalytics") 83 if (as.mat) M4 <- M4.vec2mat(M4, PP) 84 85 return( M4 ) 86} 87 88multivariate_mean <- function(w, mu) { 89 return( t(w) %*% mu ) 90} 91 92StdDev.MM <- function(w, sigma) { 93 return( sqrt(t(w) %*% sigma %*% w) ) 94} 95 96skewness.MM <- function(w, sigma, M3) { 97 w <- as.numeric(w) 98 if (NCOL(M3) != 1) M3 <- M3.mat2vec(M3) 99 m3_univ <- .Call('M3port', w, M3, length(w), PACKAGE="PerformanceAnalytics") 100 return( m3_univ / (StdDev.MM(w, sigma))^3 ) 101} 102 103kurtosis.MM <- function(w, sigma, M4) { 104 w <- as.numeric(w) 105 if (NCOL(M4) != 1) M4 <- M4.mat2vec(M4) 106 m4_univ <- .Call('M4port', w, M4, length(w), PACKAGE="PerformanceAnalytics") 107 return( m4_univ / (StdDev.MM(w, sigma))^4 ) 108} 109 110SR.StdDev.MM <- function(w, mu, sigma) { 111 return( multivariate_mean(w, mu) / StdDev.MM(w, sigma) ) 112} 113 114GVaR.MM <- function(w, mu, sigma, p) { 115 return( -multivariate_mean(w, mu) - qnorm(1 - p) * StdDev.MM(w, sigma) ) 116} 117 118SR.GVaR.MM <- function(w, mu, sigma, p) { 119 return( multivariate_mean(w, mu) / GVaR.MM(w, mu, sigma, p) ) 120} 121 122mVaR.MM <- function(w, mu, sigma, M3, M4, p) { 123 skew <- skewness.MM(w, sigma, M3) 124 exkurt <- kurtosis.MM(w, sigma, M4) - 3 125 z <- qnorm(1 - p) 126 zc <- z + (1 / 6) * (z^2 - 1) * skew 127 Zcf <- zc + (1 / 24) * (z^3 - 3 * z) * exkurt - (1 / 36) * (2 * z^3 - 5 * z) * skew^2 128 return( -multivariate_mean(w, mu) - Zcf * StdDev.MM(w, sigma) ) 129} 130 131SR.mVaR.MM <- function(w, mu, sigma, M3, M4, p) { 132 return( multivariate_mean(w, mu) / mVaR.MM(w, mu, sigma, M3, M4, p) ) 133} 134 135GES.MM <- function(w, mu, sigma, p) { 136 return( -multivariate_mean(w, mu) + dnorm(qnorm(1 - p)) * StdDev.MM(w, sigma) / (1 - p) ) 137} 138 139SR.GES.MM <- function(w, mu, sigma, p) { 140 return( multivariate_mean(w, mu) / GES.MM(w, mu, sigma, p) ) 141} 142 143Ipower = function(power,h){ 144 # probably redundant now 145 fullprod = 1; 146 if( (power%%2)==0 ) #even number: number mod is zero 147 { 148 pstar = power/2; 149 for(j in c(1:pstar)){ 150 fullprod = fullprod*(2*j) } 151 I = fullprod*dnorm(h); 152 153 for(i in c(1:pstar) ) 154 { 155 prod = 1; 156 for(j in c(1:i) ){ 157 prod = prod*(2*j) } 158 I = I + (fullprod/prod)*(h^(2*i))*dnorm(h) 159 } 160 }else{ 161 pstar = (power-1)/2 162 for(j in c(0:pstar) ) { 163 fullprod = fullprod*( (2*j)+1 ) } 164 I = -fullprod*pnorm(h); 165 for(i in c(0:pstar) ){ 166 prod = 1; 167 for(j in c(0:i) ){ 168 prod = prod*( (2*j) + 1 )} 169 I = I + (fullprod/prod)*(h^( (2*i) + 1))*dnorm(h) } 170 } 171 return(I) 172} 173 174mES.MM <- function(w, mu, sigma, M3, M4, p) { 175 skew <- skewness.MM(w, sigma, M3) 176 exkurt <- kurtosis.MM(w, sigma, M4) - 3 177 z <- qnorm(1 - p) 178 h <- z + (1 / 6) * (z^2 - 1) * skew 179 h <- h + (1 / 24) * (z^3 - 3 * z) * exkurt - (1 / 36) * (2 * z^3 - 5 * z) * skew^2 180 181 # E = dnorm(h) 182 # E = E + (1/24)*( Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h) )*exkurt 183 # E = E + (1/6)*( Ipower(3,h) - 3*Ipower(1,h) )*skew; 184 # E = E + (1/72)*( Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2) 185 # E = E/(1-p) 186 E <- dnorm(h) * (1 + h^3 * skew / 6.0 + 187 (h^6 - 9 * h^4 + 9 * h^2 + 3) * skew^2 / 72 + 188 (h^4 - 2 * h^2 - 1) * exkurt / 24) / (1 - p) 189 190 return( -multivariate_mean(w, mu) - StdDev.MM(w, sigma) * min(-E, h) ) 191} 192 193SR.mES.MM <- function(w, mu, sigma, M3, M4, p){ 194 return( multivariate_mean(w, mu) / mES.MM(w, mu, sigma, M3, M4, p) ) 195} 196 197# Computes inner product between M3_1 and M3_2 and otherwise the Frobenius norm of M3_1 198M3.innprod <- function(p, M3_1, M3_2 = NULL) { 199 # p : dimension of the data 200 # M3_1 : numeric vector with unique coskewness elements (p * (p + 1) * (p + 2) / 6) 201 # M3_2 : numeric vector with unique coskewness elements (p * (p + 1) * (p + 2) / 6) 202 203 if (NCOL(M3_1) != 1) stop("M3_1 should only contain the unique coskewness elements (see e.g. output of M3.MM(R, as.mat = FALSE))") 204 205 if (is.null(M3_2)) { 206 M3_2 <- M3_1 207 } else { 208 if (length(M3_1) != length(M3_2)) stop("M3_2 should only contain the unique coskewness elements (see e.g. output of M3.MM(R, as.mat = FALSE))") 209 } 210 211 .Call('M3innprod', M3_1, M3_2, as.integer(p), PACKAGE="PerformanceAnalytics") 212} 213 214# Computes inner product between M4_1 and M4_2 and otherwise the Frobenius norm of M4_1 215M4.innprod <- function(p, M4_1, M4_2 = NULL) { 216 # p : dimension of the data 217 # M4_1 : numeric vector with unique coskewness elements (p * (p + 1) * (p + 2) * (p + 3) / 24) 218 # M4_2 : numeric vector with unique coskewness elements (p * (p + 1) * (p + 2) * (p + 3) / 24) 219 220 if (NCOL(M4_1) != 1) stop("M4_1 should only contain the unique coskewness elements (see e.g. output of M4.MM(R, as.mat = FALSE))") 221 222 if (is.null(M4_2)) { 223 M4_2 <- M4_1 224 } else { 225 if (length(M4_1) != length(M4_2)) stop("M4_2 should only contain the unique coskewness elements (see e.g. output of M4.MM(R, as.mat = FALSE))") 226 } 227 228 .Call('M4innprod', M4_1, M4_2, as.integer(p), PACKAGE="PerformanceAnalytics") 229} 230 231# Wrapper function for casting the vector with unique coskewness elements into the matrix format 232M3.vec2mat <- function(M3, p) { 233 # M3 : numeric vector with unique coskewness elements (p * (p + 1) * (p + 2) / 6) 234 # p : dimension of the data 235 236 if (NCOL(M3) != 1) stop("M3 must be a vector") 237 238 .Call('M3vec2mat', M3, as.integer(p), PACKAGE="PerformanceAnalytics") 239} 240 241# Wrapper function for casting the coskewness matrix into the vector of unique elements 242M3.mat2vec <- function(M3) { 243 # M3 : numeric matrix of dimension p x p^2 244 245 if (is.null(dim(M3))) stop("M3 must be a matrix") 246 247 .Call('M3mat2vec', as.numeric(M3), NROW(M3), PACKAGE="PerformanceAnalytics") 248} 249 250# Wrapper function for casting the vector with unique cokurtosis elements into the matrix format 251M4.vec2mat <- function(M4, p) { 252 # M4 : numeric vector with unique cokurtosis elements (p * (p + 1) * (p + 2) * (p + 3) / 24) 253 # p : dimension of the data 254 255 if (NCOL(M4) != 1) stop("M4 must be a vector") 256 257 .Call('M4vec2mat', M4, as.integer(p), PACKAGE="PerformanceAnalytics") 258} 259 260# Wrapper function for casting the cokurtosis matrix into the vector of unique elements 261M4.mat2vec <- function(M4) { 262 # M4 : numeric matrix of dimension p x p^3 263 264 if (is.null(dim(M4))) stop("M4 must be a matrix") 265 266 .Call('M4mat2vec', as.numeric(M4), NROW(M4), PACKAGE="PerformanceAnalytics") 267} 268 269#' Functions for calculating shrinkage-based comoments of financial time series 270#' 271#' calculates covariance, coskewness and cokurtosis matrices using linear shrinkage 272#' between the sample estimator and a structured estimator 273#' 274#' The coskewness and cokurtosis matrices are defined as the matrices of dimension 275#' p x p^2 and p x p^3 containing the third and fourth order central moments. They 276#' are useful for measuring nonlinear dependence between different assets of the 277#' portfolio and computing modified VaR and modified ES of a portfolio. 278#' 279#' Shrinkage estimation for the covariance matrix was popularized by Ledoit and 280#' Wolf (2003, 2004). An extension to coskewness and cokurtosis matrices by 281#' Martellini and Ziemann (2010) uses the 1-factor and constant-correlation structured 282#' comoment matrices as targets. In Boudt, Cornilly and Verdonck (2017) the framework 283#' of single-target shrinkage for the coskewness and cokurtosis matrices is 284#' extended to a multi-target setting, making it possible to include several target matrices 285#' at once. Also, an option to enhance small sample performance for coskewness estimation 286#' was proposed, resulting in the option 'unbiasedMSE' present in the 'M3.shrink' function. 287#' 288#' The first four target matrices of the 'M2.shrink', 'M3.shrink' and 'M4.shrink' 289#' correspond to the models 'independent marginals', 'independent and identical marginals', 290#' 'observed 1-factor model' and 'constant correlation'. Coskewness shrinkage includes two 291#' more options, target 5 corresponds to the latent 1-factor model proposed in Simaan (1993) 292#' and target 6 is the coskewness matrix under central-symmetry, a matrix full of zeros. 293#' For more details on the targets, we refer to Boudt, Cornilly and Verdonck (2017) and 294#' the supplementary appendix. 295#' 296#' If f is a matrix containing multiple factors, then the shrinkage estimator will 297#' use each factor in a seperate single-factor model and use multi-target shrinkage 298#' to all targets matrices at once. 299#' @name ShrinkageMoments 300#' @concept co-moments 301#' @concept moments 302#' @aliases ShrinkageMoments M2.shrink M3.shrink M4.shrink 303#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of 304#' asset returns 305#' @param targets vector of integers selecting the target matrices to shrink to. The first four 306#' structures are, in order: 'independent marginals', 'independent and identical marginals', 307#' 'observed 1-factor model' and 'constant correlation'. See Details. 308#' @param f vector or matrix with observations of the factor, to be used with target 3. See Details. 309#' @param unbiasedMSE TRUE/FALSE whether to use a correction to have an unbiased 310#' estimator for the marginal skewness values, in case of targets 1 and/or 2, default FALSE 311#' @param as.mat TRUE/FALSE whether to return the full moment matrix or only 312#' the vector with the unique elements (the latter is advised for speed), default 313#' TRUE 314#' @author Dries Cornilly 315#' @seealso \code{\link{CoMoments}} \cr \code{\link{StructuredMoments}} \cr \code{\link{EWMAMoments}} \cr \code{\link{MCA}} 316#' @references Boudt, Kris, Brian G. Peterson, and Christophe Croux. 2008. 317#' Estimation and Decomposition of Downside Risk for Portfolios with Non-Normal 318#' Returns. Journal of Risk. Winter. 319#' 320#' Boudt, Kris, Cornilly, Dries and Verdonck, Tim. 2017. A Coskewness Shrinkage 321#' Approach for Estimating the Skewness of Linear Combinations of Random Variables. 322#' Submitted. Available at SSRN: https://ssrn.com/abstract=2839781 323#' 324#' Ledoit, Olivier and Wolf, Michael. 2003. Improved estimation of the covariance matrix 325#' of stock returns with an application to portfolio selection. Journal of empirical 326#' finance, 10(5), 603-621. 327#' 328#' Ledoit, Olivier and Wolf, Michael. 2004. A well-conditioned estimator for large-dimensional 329#' covariance matrices. Journal of multivariate analysis, 88(2), 365-411. 330#' 331#' Martellini, Lionel and Ziemann, V\"olker. 2010. Improved estimates of higher-order 332#' comoments and implications for portfolio selection. Review of Financial 333#' Studies, 23(4), 1467-1502. 334#' 335#' Simaan, Yusif. 1993. Portfolio selection and asset pricing: three-parameter framework. 336#' Management Science, 39(5), 68-577. 337###keywords ts multivariate distribution models 338#' @examples 339#' 340#' data(edhec) 341#' 342#' # construct an underlying factor (market-factor, observed factor, PCA, ...) 343#' f <- rowSums(edhec) 344#' 345#' # multi-target shrinkage with targets 1, 3 and 4 346#' # as.mat = F' would speed up calculations in higher dimensions 347#' targets <- c(1, 3, 4) 348#' sigma <- M2.shrink(edhec, targets, f)$M2sh 349#' m3 <- M3.shrink(edhec, targets, f)$M3sh 350#' m4 <- M4.shrink(edhec, targets, f)$M4sh 351#' 352#' # compute equal-weighted portfolio modified ES 353#' mu <- colMeans(edhec) 354#' p <- length(mu) 355#' ES(p = 0.95, portfolio_method = "component", weights = rep(1 / p, p), mu = mu, 356#' sigma = sigma, m3 = m3, m4 = m4) 357#' 358#' # compare to sample method 359#' sigma <- cov(edhec) 360#' m3 <- M3.MM(edhec) 361#' m4 <- M4.MM(edhec) 362#' ES(p = 0.95, portfolio_method = "component", weights = rep(1 / p, p), mu = mu, 363#' sigma = sigma, m3 = m3, m4 = m4) 364#' 365#' @export M2.shrink 366M2.shrink <- function(R, targets = 1, f = NULL) { 367 # @author Dries Cornilly 368 # 369 # DESCRIPTION: 370 # computes the shrinkage estimator of the covariance matrix 371 # 372 # Inputs: 373 # R : numeric matrix of dimensions NN x PP 374 # targets : vector of integers indicating which targets to use 375 # : T1 : independent, unequal marginals, Ledoit and Wolf (2003) 376 # : T2 : independent, equal marginals, Ledoit and Wolf (2003) 377 # : T3 : 1-factor model of Ledoit and Wolf (2003) 378 # : if multiple factors are provided, additionall 1-factor structured matrices are added the end 379 # : T4 : constant-correlation model of Ledoit and Wolf (2004) 380 # f : numeric vector with factor observations, needed for 1-factor coskewness matrix of Martellini and Ziemann 381 # : or a numeric matrix with columns as factors 382 # as.mat : output as a matrix or as the vector with only unique coskewness eleements 383 # 384 # Outputs: 385 # M2sh : the shrinkage estimator 386 # lambda : vector with shrinkage intensities 387 # A : A matrix in the QP 388 # b : b vector in QP 389 390 X <- coredata(R) 391 392 # input checking 393 if (NCOL(X) < 2) stop("R must have at least 2 variables") 394 if (length(targets) == 0) stop("No targets selected") 395 if (prod(targets %in% 1:4) == 0) stop("Select valid targets (out of 1, 2, 3, 4)") 396 tt <- rep(FALSE, 4) 397 tt[targets] <- TRUE 398 targets <- tt 399 if (targets[3] && is.null(f)) stop("Provide the factor observations") 400 401 # prepare for additional factors if necessary 402 if (targets[3] && (NCOL(f) != 1)) { 403 nFactors <- NCOL(f) 404 if (nFactors > 1) { 405 f_other <- matrix(f[, 2:nFactors], ncol = nFactors - 1) 406 f <- f[, 1] 407 extraFactors <- TRUE 408 targets <- c(targets, rep(TRUE, nFactors - 1)) 409 } else { 410 f <- c(f) 411 extraFactors <- FALSE 412 } 413 } else { 414 extraFactors <- FALSE 415 } 416 417 # compute useful variables 418 NN <- dim(X)[1] # number of observations 419 PP <- dim(X)[2] # number of assets 420 nT <- sum(targets) # number of targets 421 422 Xc <- X - matrix(colMeans(X), nrow = NN, ncol = PP, byrow = TRUE) # center the observations 423 Xc2 <- Xc^2 424 margvars <- colMeans(Xc2) 425 m11 <- as.numeric(t(Xc) %*% Xc) / NN 426 m22 <- as.numeric(t(Xc2) %*% Xc2) / NN 427 428 ### coskewness estimators 429 M2 <- cov(X) * (NN - 1) / NN 430 T2 <- matrix(NA, nrow = PP * PP, ncol = nT) 431 iter <- 1 432 433 if (targets[1]) { 434 # independent marginals 435 T2[, iter] <- c(diag(margvars)) 436 iter <- iter + 1 437 } 438 if (targets[2]) { 439 # independent and equally distributed marginals 440 T2[, iter] <- c(mean(margvars) * diag(PP)) 441 iter <- iter + 1 442 } 443 if (targets[3]) { 444 # 1-factor model (Ledoit and Wolf (2003)) 445 beta <- apply(Xc, 2, function(a) cov(a, f) / var(f)) 446 fc <- f - mean(f) 447 fvar <- mean(fc^2) 448 T2_1f <- fvar * beta %*% t(beta) 449 diag(T2_1f) <- margvars 450 T2[, iter] <- T2_1f 451 iter <- iter + 1 452 } 453 if (targets[4]) { 454 # constant-correlation (Ledoit and Wolf (2004)) 455 sd_vec <- sqrt(margvars) 456 R2 <- diag(1 / sd_vec) %*% M2 %*% diag(1 / sd_vec) 457 rcoef <- mean(R2[upper.tri(R2)]) 458 R2 <- matrix(rcoef, nrow = PP, ncol = PP) 459 diag(R2) <- 1 460 T2[, iter] <- diag(sd_vec) %*% R2 %*% diag(sd_vec) 461 iter <- iter + 1 462 } 463 if (extraFactors) { 464 # 1-factor model (Ledoit and Wolf (2003)) - extra factors 465 for (ii in 1:(nFactors - 1)) { 466 f_bis <- f_other[, ii] 467 beta_bis <- apply(Xc, 2, function(a) cov(a, f_bis) / var(f_bis)) 468 fc_bis <- f_bis - mean(f_bis) 469 fvar_bis <- mean(fc_bis^2) 470 T2_1f <- fvar_bis * beta_bis %*% t(beta_bis) 471 diag(T2_1f) <- margvars 472 T2[, iter] <- T2_1f 473 iter <- iter + 1 474 } 475 } 476 477 ### build A for the QP 478 A <- matrix(NA, nrow = nT, ncol = nT) 479 for (ii in 1:nT) { 480 for (jj in ii:nT) { 481 A[ii, jj] <- A[jj, ii] <- sum((T2[, ii] - M2) * (T2[, jj] - M2)) 482 } 483 } 484 485 ### build b for the QP 486 VM2vec <- .Call('VM2', m11, m22, NN, PP, PACKAGE="PerformanceAnalytics") 487 488 b <- rep(VM2vec[1], nT) 489 iter <- 1 490 if (targets[1]) { 491 # independent marginals 492 b[iter] <- b[iter] - VM2vec[3] 493 iter <- iter + 1 494 } 495 if (targets[2]) { 496 # independent and equally distributed marginals 497 b[iter] <- b[iter] - VM2vec[2] 498 iter <- iter + 1 499 } 500 if (targets[3]) { 501 # 1-factor model (Ledoit and Wolf (2003)) 502 b[iter] <- b[iter] - .Call('CM2_1F', as.numeric(Xc), fc, fvar, m11, m22, NN, PP, PACKAGE="PerformanceAnalytics") 503 iter <- iter + 1 504 } 505 if (targets[4]) { 506 # constant-correlation (Ledoit and Wolf (2004)) 507 b[iter] <- b[iter] - .Call('CM2_CC', as.numeric(Xc), rcoef, m11, m22, NN, PP, PACKAGE="PerformanceAnalytics") 508 iter <- iter + 1 509 } 510 if (extraFactors) { 511 # 1-factor model (Ledoit and Wolf (2003)) - extra factors 512 for (ii in 1:(nFactors - 1)) { 513 f_bis <- f_other[, ii] 514 fc_bis <- f_bis - mean(f_bis) 515 fvar_bis <- mean(fc_bis^2) 516 b[iter] <- b[iter] - .Call('CM2_1F', as.numeric(Xc), fc_bis, fvar_bis, m11, m22, NN, PP, PACKAGE="PerformanceAnalytics") 517 iter <- iter + 1 518 } 519 } 520 521 ### solve the QP 522 if (nT == 1) { 523 # single-target shrinkage 524 lambda <- b / A # compute optimal shrinkage intensity 525 lambda <- max(0, min(1, lambda)) # must be between 0 and 1 526 M2sh <- (1 - lambda) * M2 + lambda * T2[, 1] # compute shrinkage estimator 527 } else { 528 # multi-target shrinkage 529 Aineq <- rbind(diag(nT), rep(-1, nT)) # A matrix for inequalities quadratic program 530 bineq <- matrix(c(rep(0, nT), -1), ncol = 1) # b vector for inequalities quadratic program 531 lambda <- quadprog::solve.QP(A, b, t(Aineq), bineq, meq = 0)$solution # solve quadratic program 532 M2sh <- (1 - sum(lambda)) * M2 # initialize estimator at percentage of sample estimator 533 for (tt in 1:nT) { 534 M2sh <- M2sh + lambda[tt] * T2[, tt] # add the target matrices 535 } 536 } 537 538 return (list("M2sh" = M2sh, "lambda" = lambda, "A" = A, "b" = b)) 539} 540 541#'@useDynLib PerformanceAnalytics 542#'@export 543#'@rdname ShrinkageMoments 544M3.shrink <- function(R, targets = 1, f = NULL, unbiasedMSE = FALSE, as.mat = TRUE) { 545 # @author Dries Cornilly 546 # 547 # DESCRIPTION: 548 # computes the shrinkage estimator of the coskewness matrix as in Boudt, Cornilly and Verdonck (2017) 549 # 550 # Inputs: 551 # R : numeric matrix of dimensions NN x PP 552 # targets : vector of integers indicating which targets to use 553 # : T1 : independent, unequal marginals 554 # : T2 : independent, equal marginals 555 # : T3 : 1-factor model of Martellini and Ziemann (2010) 556 # : if multiple factors are provided, additionall 1-factor structured matrices are added the end 557 # : T4 : constant-correlation model of Martellini and Ziemann (2010), symmetrized 558 # : T5 : latent 1-factor model of Simaan (1993) 559 # : T6 : central-symmetric coskewness matrix (all zeros) 560 # f : numeric vector with factor observations, needed for 1-factor coskewness matrix of Martellini and Ziemann 561 # : or a numeric matrix with columns as factors 562 # unbiasedMSE : boolean determining if bias is corrected when estimating the MSE loss function 563 # as.mat : output as a matrix or as the vector with only unique coskewness eleements 564 # 565 # Outputs: 566 # M3sh : the shrinkage estimator 567 # lambda : vector with shrinkage intensities 568 # A : A matrix in the QP 569 # b : b vector in QP 570 571 X <- coredata(R) 572 573 # input checking 574 if (NCOL(X) < 2) stop("R must have at least 2 variables") 575 if (length(targets) == 0) stop("No targets selected") 576 if (prod(targets %in% 1:6) == 0) stop("Select valid targets (out of 1, 2, 3, 4, 5, 6)") 577 tt <- rep(FALSE, 6) 578 tt[targets] <- TRUE 579 targets <- tt 580 if (targets[3] && is.null(f)) stop("Provide the factor observations") 581 if (unbiasedMSE && (sum(targets[c(3, 4, 5)]) > 0)) stop("UnbiasedMSE can only be combined with T2, T3 and T6") 582 583 # prepare for additional factors if necessary 584 if (targets[3] && (NCOL(f) != 1)) { 585 nFactors <- NCOL(f) 586 if (nFactors > 1) { 587 f_other <- matrix(f[, 2:nFactors], ncol = nFactors - 1) 588 f <- f[, 1] 589 extraFactors <- TRUE 590 targets <- c(targets, rep(TRUE, nFactors - 1)) 591 } else { 592 f <- c(f) 593 extraFactors <- FALSE 594 } 595 } else { 596 extraFactors <- FALSE 597 } 598 599 # compute useful variables 600 NN <- dim(X)[1] # number of observations 601 if (unbiasedMSE) { 602 if (NN < 6) stop("R should have at least 6 observations") 603 } 604 PP <- dim(X)[2] # number of assets 605 nT <- sum(targets) # number of targets 606 ncosk <- PP * (PP + 1) * (PP + 2) / 6 # number of unique coskewness elements 607 608 Xc <- X - matrix(colMeans(X), nrow = NN, ncol = PP, byrow = TRUE) # center the observations 609 Xc2 <- Xc^2 610 m11 <- as.numeric(t(Xc) %*% Xc) / NN 611 m21 <- as.numeric(t(Xc2) %*% Xc) / NN 612 m22 <- as.numeric(t(Xc2) %*% Xc2) / NN 613 m31 <- as.numeric(t(Xc^3) %*% Xc) / NN 614 m42 <- as.numeric(t(Xc^4) %*% Xc2) / NN 615 m33 <- as.numeric(t(Xc^3) %*% Xc^3) / NN 616 617 ### coskewness estimators 618 M3 <- M3.MM(X, unbiased = unbiasedMSE, as.mat = FALSE) 619 T3 <- matrix(NA, nrow = length(M3), ncol = nT) 620 iter <- 1 621 622 if (targets[1]) { 623 # independent marginals 624 margskews <- colMeans(Xc^3) 625 if (unbiasedMSE) margskews <- margskews * NN^2 / ((NN - 1) * (NN - 2)) 626 T3[, iter] <- .Call('M3_T23', margskews, PP, PACKAGE="PerformanceAnalytics") 627 iter <- iter + 1 628 } 629 if (targets[2]) { 630 # independent and equally distributed marginals 631 margskews <- colMeans(Xc^3) 632 if (unbiasedMSE) margskews <- margskews * NN^2 / ((NN - 1) * (NN - 2)) 633 margskews <- rep(mean(margskews), PP) 634 T3[, iter] <- .Call('M3_T23', margskews, PP, PACKAGE="PerformanceAnalytics") 635 iter <- iter + 1 636 } 637 if (targets[3]) { 638 # 1-factor model (Martellini and Ziemann (2010)) 639 margskews <- colMeans(Xc^3) 640 beta <- apply(Xc, 2, function(a) cov(a, f) / var(f)) 641 fc <- f - mean(f) 642 fskew <- mean(fc^3) 643 T3[, iter] <- .Call('M3_1F', margskews, beta, fskew, PP, PACKAGE="PerformanceAnalytics") 644 iter <- iter + 1 645 } 646 if (targets[4]) { 647 # constant-correlation (Martellini and Ziemann (2010)) (symmetrized version) 648 margvars <- colMeans(Xc^2) 649 margskews <- colMeans(Xc^3) 650 margkurts <- colMeans(Xc^4) 651 r_generalized <- .Call('M3_CCoefficients', margvars, margkurts, m21, 652 m22, as.numeric(Xc), NN, PP, PACKAGE="PerformanceAnalytics") 653 T3[, iter] <- .Call('M3_CC', margvars, margskews, margkurts, 654 r_generalized[1], r_generalized[2], 655 r_generalized[3], PP, PACKAGE="PerformanceAnalytics") 656 iter <- iter + 1 657 } 658 if (targets[5]) { 659 # 1-factor model of Simaan (1993) 660 margskews <- colMeans(Xc^3) 661 margskewsroot <- sign(margskews) * abs(margskews)^(1 / 3) 662 T3[, iter] <- .Call('M3_Simaan', margskewsroot, PP, PACKAGE="PerformanceAnalytics") 663 iter <- iter + 1 664 } 665 if (targets[6]) { 666 # central-symmetric 667 T3[, iter] <- rep(0, ncosk) 668 iter <- iter + 1 669 } 670 if (extraFactors) { 671 # 1-factor model (Martellini and Ziemann (2010)) - extra factors 672 for (ii in 1:(nFactors - 1)) { 673 f_bis <- f_other[, ii] 674 margskews <- colMeans(Xc^3) 675 beta_bis <- apply(Xc, 2, function(a) cov(a, f_bis) / var(f_bis)) 676 fc_bis <- f_bis - mean(f_bis) 677 fskew_bis <- mean(fc_bis^3) 678 T3[, iter] <- .Call('M3_1F', margskews, beta_bis, fskew_bis, PP, PACKAGE="PerformanceAnalytics") 679 iter <- iter + 1 680 } 681 } 682 683 ### build A for the QP 684 A <- matrix(NA, nrow = nT, ncol = nT) 685 for (ii in 1:nT) { 686 for (jj in ii:nT) { 687 A[ii, jj] <- A[jj, ii] <- .Call('M3innprod', T3[, ii] - M3, T3[, jj] - M3, PP, PACKAGE="PerformanceAnalytics") 688 } 689 } 690 691 ### build b for the QP 692 if (unbiasedMSE) { 693 VM3vec <- .Call('VM3kstat', as.numeric(Xc), as.numeric(Xc2), m11 * NN, m21 * NN, m22 * NN, 694 m31 * NN, m42 * NN, m33 * NN, NN, PP, PACKAGE="PerformanceAnalytics") 695 } else { 696 VM3vec <- .Call('VM3', as.numeric(Xc), as.numeric(Xc2), m11, m21, m22, m31, m42, 697 m33, NN, PP, PACKAGE="PerformanceAnalytics") 698 } 699 700 b <- rep(VM3vec[1], nT) 701 iter <- 1 702 if (targets[1]) { 703 # independent marginals 704 b[iter] <- b[iter] - VM3vec[3] 705 iter <- iter + 1 706 } 707 if (targets[2]) { 708 # independent and equally distributed marginals 709 b[iter] <- b[iter] - VM3vec[2] 710 iter <- iter + 1 711 } 712 if (targets[3]) { 713 # 1-factor model (Martellini and Ziemann (2010)) 714 fvar <- mean(fc^2) 715 b[iter] <- b[iter] - .Call('CM3_1F', as.numeric(Xc), as.numeric(Xc2), 716 fc, fvar, fskew, m11, m21, m22, m42, NN, PP, PACKAGE="PerformanceAnalytics") 717 iter <- iter + 1 718 } 719 if (targets[4]) { 720 # constant-correlation (Martellini and Ziemann (2010)) (symmetrized version) 721 marg5s <- colMeans(Xc^5) 722 marg6s <- colMeans(Xc^6) 723 m41 <- as.numeric(t(Xc^4) %*% Xc) / NN 724 m61 <- as.numeric(t(Xc^6) %*% Xc) / NN 725 m32 <- as.numeric(t(Xc^3) %*% Xc^2) / NN 726 b[iter] <- b[iter] - .Call('CM3_CC', as.numeric(Xc), as.numeric(Xc2), margvars, margskews, margkurts, 727 marg5s, marg6s, m11, m21, m31, m32, m41, m61, r_generalized[1], r_generalized[2], 728 r_generalized[3], NN, PP, PACKAGE="PerformanceAnalytics") 729 iter <- iter + 1 730 } 731 if (targets[5]) { 732 # 1-factor model of Simaan (1993) 733 margskewsroot <- margskewsroot^(-2) 734 m51 <- as.numeric(t(Xc^5) %*% Xc) / NN 735 b[iter] <- b[iter] - .Call('CM3_Simaan', as.numeric(Xc), as.numeric(Xc2), margskewsroot, m11, m21, m22, 736 m31, m42, m51, NN, PP, PACKAGE="PerformanceAnalytics") 737 iter <- iter + 1 738 } 739 if (targets[6]) iter <- iter + 1 740 if (extraFactors) { 741 # 1-factor model (Martellini and Ziemann (2010)) - extra factors 742 for (ii in 1:(nFactors - 1)) { 743 f_bis <- f_other[, ii] 744 fc_bis <- f_bis - mean(f_bis) 745 fvar_bis <- mean(fc_bis^2) 746 fskew_bis <- mean(fc_bis^3) 747 b[iter] <- b[iter] - .Call('CM3_1F', as.numeric(Xc), as.numeric(Xc2), fc_bis, fvar_bis, fskew_bis, 748 m11, m21, m22, m42, NN, PP, PACKAGE="PerformanceAnalytics") 749 iter <- iter + 1 750 } 751 } 752 753 ### solve the QP 754 if (nT == 1) { 755 # single-target shrinkage 756 lambda <- b / A # compute optimal shrinkage intensity 757 lambda <- max(0, min(1, lambda)) # must be between 0 and 1 758 M3sh <- (1 - lambda) * M3 + lambda * T3 # compute shrinkage estimator 759 } else { 760 # multi-target shrinkage 761 Aineq <- rbind(diag(nT), rep(-1, nT)) # A matrix for inequalities quadratic program 762 bineq <- matrix(c(rep(0, nT), -1), ncol = 1) # b vector for inequalities quadratic program 763 lambda <- quadprog::solve.QP(A, b, t(Aineq), bineq, meq = 0)$solution # solve quadratic program 764 M3sh <- (1 - sum(lambda)) * M3 # initialize estimator at percentage of sample estimator 765 for (tt in 1:nT) { 766 M3sh <- M3sh + lambda[tt] * T3[, tt] # add the target matrices 767 } 768 } 769 if (as.mat) M3sh <- M3.vec2mat(M3sh, PP) 770 771 return (list("M3sh" = M3sh, "lambda" = lambda, "A" = A, "b" = b)) 772} 773 774 775#'@useDynLib PerformanceAnalytics 776#'@export 777#'@rdname ShrinkageMoments 778M4.shrink <- function(R, targets = 1, f = NULL, as.mat = TRUE) { 779 # @author Dries Cornilly 780 # 781 # DESCRIPTION: 782 # computes the shrinkage estimator of the cokurtosis matrix as in Boudt, Cornilly and Verdonck (2017) 783 # 784 # Inputs: 785 # R : numeric matrix of dimensions NN x PP 786 # targets : vector of integers indicating which targets to use 787 # : T1 : independent, unequal marginals 788 # : T2 : independent, equal marginals 789 # : T3 : 1-factor model of Martellini and Ziemann (2010) 790 # : if multiple factors are provided, additionall 1-factor structured matrices are added the end 791 # : T4 : constant-correlation model of Martellini and Ziemann (2010) 792 # f : numeric vector with factor observations, needed for 1-factor coskewness matrix of Martellini and Ziemann 793 # : or a numeric matrix with columns as factors 794 # as.mat : output as a matrix or as the vector with only unique coskewness eleements 795 # 796 # Outputs: 797 # M4sh : the shrinkage estimator 798 # lambda : vector with shrinkage intensities 799 # A : A matrix in the QP 800 # b : b vector in QP 801 802 X <- coredata(R) 803 804 # input checking 805 if (NCOL(X) < 2) stop("R must have at least 2 variables") 806 if (length(targets) == 0) stop("No targets selected") 807 if (prod(targets %in% 1:4) == 0) stop("Select valid targets (out of 1, 2, 3, 4)") 808 tt <- rep(FALSE, 4) 809 tt[targets] <- TRUE 810 targets <- tt 811 if (targets[3] && is.null(f)) stop("Provide the factor observations") 812 813 # prepare for additional factors if necessary 814 if (targets[3] && (NCOL(f) != 1)) { 815 nFactors <- NCOL(f) 816 if (nFactors > 1) { 817 f_other <- matrix(f[, 2:nFactors], ncol = nFactors - 1) 818 f <- f[, 1] 819 extraFactors <- TRUE 820 targets <- c(targets, rep(TRUE, nFactors - 1)) 821 } else { 822 f <- c(f) 823 extraFactors <- FALSE 824 } 825 } else { 826 extraFactors <- FALSE 827 } 828 829 # compute useful variables 830 NN <- dim(X)[1] # number of observations 831 PP <- dim(X)[2] # number of assets 832 nT <- sum(targets) # number of targets 833 834 Xc <- X - matrix(colMeans(X), nrow = NN, ncol = PP, byrow = TRUE) # center the observations 835 Xc2 <- Xc^2 836 margvars <- colMeans(Xc2) 837 margkurts <- colMeans(Xc^4) 838 m11 <- as.numeric(t(Xc) %*% Xc) / NN 839 m21 <- as.numeric(t(Xc2) %*% Xc) / NN 840 m22 <- as.numeric(t(Xc2) %*% Xc2) / NN 841 m31 <- as.numeric(t(Xc^3) %*% Xc) / NN 842 m32 <- as.numeric(t(Xc^3) %*% Xc2) / NN 843 m41 <- as.numeric(t(Xc^4) %*% Xc) / NN 844 m42 <- as.numeric(t(Xc^4) %*% Xc2) / NN 845 846 ### coskewness estimators 847 M4 <- M4.MM(X, as.mat = FALSE) 848 T4 <- matrix(NA, nrow = length(M4), ncol = nT) 849 iter <- 1 850 851 if (targets[1]) { 852 # independent marginals 853 T4[, iter] <- .Call('M4_T12', margkurts, margvars, PP, PACKAGE="PerformanceAnalytics") 854 iter <- iter + 1 855 } 856 if (targets[2]) { 857 # independent and equally distributed marginals 858 meanmargkurts <- mean(margkurts) 859 meank_iikk <- sqrt(mean(margvars^2)) 860 T4[, iter] <- .Call('M4_T12', rep(meanmargkurts, PP), rep(meank_iikk, PP), PP, PACKAGE="PerformanceAnalytics") 861 iter <- iter + 1 862 } 863 if (targets[3]) { 864 # 1-factor model (Martellini and Ziemann (2010)) 865 beta <- apply(Xc, 2, function(a) cov(a, f) / var(f)) 866 fc <- f - mean(f) 867 fvar <- mean(fc^2) 868 fkurt <- mean(fc^4) 869 epsvars <- margvars - beta^2 * fvar 870 T4[, iter] <- .Call('M4_1f', margkurts, fvar, fkurt, epsvars, beta, PP, PACKAGE="PerformanceAnalytics") 871 iter <- iter + 1 872 } 873 if (targets[4]) { 874 # constant-correlation (Martellini and Ziemann (2010)) 875 marg6s <- colMeans(Xc^6) 876 r_generalized <- .Call('M4_CCoefficients', margvars, margkurts, marg6s, 877 m22, m31, as.numeric(Xc), NN, PP, PACKAGE="PerformanceAnalytics") 878 T4[, iter] <- .Call('M4_CC', margvars, margkurts, marg6s, r_generalized[1], r_generalized[2], 879 r_generalized[3], r_generalized[4], PP, PACKAGE="PerformanceAnalytics") 880 iter <- iter + 1 881 } 882 if (extraFactors) { 883 # 1-factor model (Martellini and Ziemann (2010)) - extra factors 884 for (ii in 1:(nFactors - 1)) { 885 f_bis <- f_other[, ii] 886 beta_bis <- apply(Xc, 2, function(a) cov(a, f_bis) / var(f_bis)) 887 fc_bis <- f_bis - mean(f_bis) 888 fvar_bis <- mean(fc_bis^2) 889 fkurt_bis <- mean(fc_bis^4) 890 epsvars_bis <- margvars - beta_bis^2 * fvar_bis 891 T4[, iter] <- .Call('M4_1f', margkurts, fvar_bis, fkurt_bis, epsvars_bis, 892 beta_bis, PP, PACKAGE="PerformanceAnalytics") 893 iter <- iter + 1 894 } 895 } 896 897 ### build A for the QP 898 A <- matrix(NA, nrow = nT, ncol = nT) 899 for (ii in 1:nT) { 900 for (jj in ii:nT) { 901 A[ii, jj] <- A[jj, ii] <- .Call('M4innprod', T4[, ii] - M4, T4[, jj] - M4, PP, PACKAGE="PerformanceAnalytics") 902 } 903 } 904 905 ### build b for the QP 906 VM4vec <- .Call('VM4', as.numeric(Xc), as.numeric(Xc2), m11, m21, 907 m22, m31, m32, m41, m42, NN, PP, PACKAGE="PerformanceAnalytics") 908 909 b <- rep(VM4vec[1], nT) 910 iter <- 1 911 if (targets[1]) { 912 # independent marginals 913 b[iter] <- b[iter] - VM4vec[3] 914 iter <- iter + 1 915 } 916 if (targets[2]) { 917 # independent and equally distributed marginals 918 b[iter] <- b[iter] - VM4vec[2] 919 iter <- iter + 1 920 } 921 if (targets[3]) { 922 # 1-factor model (Martellini and Ziemann (2010)) 923 fskew <- mean(fc^3) 924 b[iter] <- b[iter] - .Call('CM4_1F', as.numeric(Xc), as.numeric(Xc2), fc, fvar, fskew, fkurt, 925 m11, m21, m22, m31, NN, PP, PACKAGE="PerformanceAnalytics") 926 iter <- iter + 1 927 } 928 if (targets[4]) { 929 # constant-correlation (Martellini and Ziemann (2010)) 930 marg6s <- colMeans(Xc^6) 931 marg7s <- colMeans(Xc^7) 932 m33 <- as.numeric(t(Xc^3) %*% Xc^3) / NN 933 b[iter] <- b[iter] - .Call('CM4_CC', as.numeric(Xc), as.numeric(Xc2), m11, m21, m22, m31, m32, m33, m41, 934 r_generalized[1], r_generalized[2], r_generalized[3], r_generalized[4], 935 marg6s, marg7s, NN, PP, PACKAGE="PerformanceAnalytics") 936 iter <- iter + 1 937 } 938 if (extraFactors) { 939 # 1-factor model (Martellini and Ziemann (2010)) - extra factors 940 for (ii in 1:(nFactors - 1)) { 941 f_bis <- f_other[, ii] 942 fc_bis <- f_bis - mean(f_bis) 943 fvar_bis <- mean(fc_bis^2) 944 fskew_bis <- mean(fc_bis^3) 945 fkurt_bis <- mean(fc_bis^4) 946 b[iter] <- b[iter] - .Call('CM4_1F', as.numeric(Xc), as.numeric(Xc2), fc_bis, fvar_bis, fskew_bis, fkurt_bis, 947 m11, m21, m22, m31, NN, PP, PACKAGE="PerformanceAnalytics") 948 iter <- iter + 1 949 } 950 } 951 952 ### solve the QP 953 if (nT == 1) { 954 # single-target shrinkage 955 lambda <- b / A # compute optimal shrinkage intensity 956 lambda <- max(0, min(1, lambda)) # must be between 0 and 1 957 M4sh <- (1 - lambda) * M4 + lambda * T4 # compute shrinkage estimator 958 } else { 959 # multi-target shrinkage 960 Aineq <- rbind(diag(nT), rep(-1, nT)) # A matrix for inequalities quadratic program 961 bineq <- matrix(c(rep(0, nT), -1), ncol = 1) # b vector for inequalities quadratic program 962 lambda <- quadprog::solve.QP(A, b, t(Aineq), bineq, meq = 0)$solution # solve quadratic program 963 M4sh <- (1 - sum(lambda)) * M4 # initialize estimator at percentage of sample estimator 964 for (tt in 1:nT) { 965 M4sh <- M4sh + lambda[tt] * T4[, tt] # add the target matrices 966 } 967 } 968 if (as.mat) M4sh <- M4.vec2mat(M4sh, PP) 969 970 return (list("M4sh" = M4sh, "lambda" = lambda, "A" = A, "b" = b)) 971} 972 973 974#' Functions for calculating structured comoments of financial time series 975#' 976#' calculates covariance, coskewness and cokurtosis matrices as structured estimators 977#' 978#' The coskewness and cokurtosis matrices are defined as the matrices of dimension 979#' p x p^2 and p x p^3 containing the third and fourth order central moments. They 980#' are useful for measuring nonlinear dependence between different assets of the 981#' portfolio and computing modified VaR and modified ES of a portfolio. 982#' 983#' Structured estimation is based on the assumption that the underlying data-generating 984#' process is known, or at least resembles the assumption. The first four structured estimators correspond to the models 'independent marginals', 985#' 'independent and identical marginals', 'observed multi-factor model' and 'constant correlation'. 986#' Coskewness estimation includes an additional model based on the latent 1-factor model 987#' proposed in Simaan (1993). 988#' 989#' The constant correlation and 1-factor coskewness and cokurtosis matrices can be found in 990#' Martellini and Ziemann (2010). If f is a matrix containing multiple factors, 991#' then the multi-factor model of Boudt, Lu and Peeters (2915) is used. For information 992#' about the other structured matrices, we refer to Boudt, Cornilly and Verdonck (2017) 993#' @name StructuredMoments 994#' @concept co-moments 995#' @concept moments 996#' @aliases StructuredMoments M2.struct M3.struct M4.struct 997#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of 998#' asset returns 999#' @param struct string containing the preferred method. See Details. 1000#' @param f vector or matrix with observations of the factor, to be used with 'observedfactor'. See Details. 1001#' @param unbiasedMarg TRUE/FALSE whether to use a correction to have an unbiased 1002#' estimator for the marginal skewness values, in case of 'Indep' or 'IndepId', default FALSE 1003#' @param as.mat TRUE/FALSE whether to return the full moment matrix or only 1004#' the vector with the unique elements (the latter is advised for speed), default 1005#' TRUE 1006#' @author Dries Cornilly 1007#' @seealso \code{\link{CoMoments}} \cr \code{\link{ShrinkageMoments}} \cr \code{\link{EWMAMoments}} \cr \code{\link{MCA}} 1008#' @references Boudt, Kris, Lu, Wanbo and Peeters, Benedict. 2015. Higher order comoments of multifactor 1009#' models and asset allocation. Finance Research Letters, 13, 225-233. 1010#' 1011#' Boudt, Kris, Brian G. Peterson, and Christophe Croux. 2008. 1012#' Estimation and Decomposition of Downside Risk for Portfolios with Non-Normal 1013#' Returns. Journal of Risk. Winter. 1014#' 1015#' Boudt, Kris, Cornilly, Dries and Verdonck, Tim. 2017. A Coskewness Shrinkage 1016#' Approach for Estimating the Skewness of Linear Combinations of Random Variables. 1017#' Submitted. Available at SSRN: https://ssrn.com/abstract=2839781 1018#' 1019#' Ledoit, Olivier and Wolf, Michael. 2003. Improved estimation of the covariance matrix 1020#' of stock returns with an application to portfolio selection. Journal of empirical 1021#' finance, 10(5), 603-621. 1022#' 1023#' Martellini, Lionel and Ziemann, V\"olker. 2010. Improved estimates of higher-order 1024#' comoments and implications for portfolio selection. Review of Financial 1025#' Studies, 23(4), 1467-1502. 1026#' 1027#' Simaan, Yusif. 1993. Portfolio selection and asset pricing: three-parameter framework. 1028#' Management Science, 39(5), 68-577. 1029###keywords ts multivariate distribution models 1030#' @examples 1031#' 1032#' data(edhec) 1033#' 1034#' # structured estimation with constant correlation model 1035#' # 'as.mat = F' would speed up calculations in higher dimensions 1036#' sigma <- M2.struct(edhec, "CC") 1037#' m3 <- M3.struct(edhec, "CC") 1038#' m4 <- M4.struct(edhec, "CC") 1039#' 1040#' # compute equal-weighted portfolio modified ES 1041#' mu <- colMeans(edhec) 1042#' p <- length(mu) 1043#' ES(p = 0.95, portfolio_method = "component", weights = rep(1 / p, p), mu = mu, 1044#' sigma = sigma, m3 = m3, m4 = m4) 1045#' 1046#' # compare to sample method 1047#' sigma <- cov(edhec) 1048#' m3 <- M3.MM(edhec) 1049#' m4 <- M4.MM(edhec) 1050#' ES(p = 0.95, portfolio_method = "component", weights = rep(1 / p, p), mu = mu, 1051#' sigma = sigma, m3 = m3, m4 = m4) 1052#' 1053#' @export M2.struct 1054M2.struct <- function(R, struct = c("Indep", "IndepId", "observedfactor", "CC"), f = NULL) { 1055 # @author Dries Cornilly 1056 # 1057 # DESCRIPTION: 1058 # computes different strutured covariance estimators 1059 # 1060 # Inputs: 1061 # R : numeric matrix of dimensions NN x PP 1062 # struct : select the structured estimator 1063 # : Indep : independent, unequal marginals (Ledoit and Wolf (2004)) 1064 # : IndepId : independent, equal marginals (Ledoit and Wolf (2003)) 1065 # : observedfactor : 1-factor or multi-factor linear model 1066 # : CC : constant-correlation model of (Ledoit and Wolf (2004)) 1067 # f : numeric vector or matrix with factor observations 1068 # 1069 # Outputs: 1070 # covariance matrix 1071 1072 X <- coredata(R) 1073 1074 # input checking 1075 struct <- match.arg(struct) 1076 if (NCOL(X) < 2) stop("R must have at least 2 variables") 1077 if ((struct == "observedfactor") && is.null(f)) stop("Provide the factor observations") 1078 1079 # compute useful variables 1080 NN <- dim(X)[1] # number of observations 1081 PP <- dim(X)[2] # number of assets 1082 1083 Xc <- X - matrix(colMeans(X), nrow = NN, ncol = PP, byrow = TRUE) # center the observations 1084 margvars <- colMeans(Xc^2) 1085 1086 # compute the coskewness matrix 1087 if (struct == "Indep") { 1088 # independent marginals 1089 T2 <- diag(margvars) 1090 return (T2) 1091 1092 } else if (struct == "IndepId") { 1093 # independent and equally distributed marginals 1094 T2 <- mean(margvars) * diag(PP) 1095 return (T2) 1096 1097 } else if (struct == "observedfactor") { 1098 # linear factor model 1099 if (NCOL(f) == 1) { 1100 beta <- apply(Xc, 2, function(a) cov(a, f) / var(f)) 1101 fc <- f - mean(f) 1102 fvar <- mean(fc^2) 1103 T2 <- fvar * beta %*% t(beta) 1104 diag(T2) <- margvars 1105 } else { 1106 mod <- stats::lm(X ~ f) 1107 beta <- t(mod$coefficients[-1,]) 1108 fcov <- cov(f) * (NN - 1) / NN 1109 T2 <- beta %*% fcov %*% t(beta) 1110 epsvars <- colMeans(mod$residuals^2) 1111 T2 <- T2 + diag(epsvars) 1112 } 1113 return (T2) 1114 1115 } else if (struct == "CC") { 1116 # constant-correlation 1117 sd_vec <- sqrt(margvars) 1118 M2 <- t(Xc) %*% Xc / NN 1119 R2 <- diag(1 / sd_vec) %*% M2 %*% diag(1 / sd_vec) 1120 rcoef <- mean(R2[upper.tri(R2)]) 1121 R2 <- matrix(rcoef, nrow = PP, ncol = PP) 1122 diag(R2) <- 1 1123 T2 <- diag(sd_vec) %*% R2 %*% diag(sd_vec) 1124 return (T2) 1125 1126 } 1127} 1128 1129 1130#'@useDynLib PerformanceAnalytics 1131#'@export 1132#'@rdname StructuredMoments 1133M3.struct <- function(R, struct = c("Indep", "IndepId", "observedfactor", "CC", "latent1factor", "CS"), 1134 f = NULL, unbiasedMarg = FALSE, as.mat = TRUE) { 1135 # @author Dries Cornilly 1136 # 1137 # DESCRIPTION: 1138 # computes different strutured coskewness estimators 1139 # 1140 # Inputs: 1141 # R : numeric matrix of dimensions NN x PP 1142 # struct : select the structured estimator 1143 # : Indep : independent, unequal marginals 1144 # : IndepId : independent, equal marginals 1145 # : observedfactor : observed linear factor model 1146 # : CC : constant-correlation model of Martellini and Ziemann (2010), symmetrized 1147 # : latent1factor : latent 1-factor model of Simaan (1993) 1148 # : CS : central-symmetric coskewness matrix (all zeros) 1149 # f : numeric vector with factor observations, 1150 # unbiasedMarg : boolean determining if bias is corrected when estimating the marginals with 1151 # : methods "Indep" or "INdepID" 1152 # 1153 # Outputs: 1154 # coskewness matrix (as matrix or as vector, depending on as.mat) 1155 1156 X <- coredata(R) 1157 1158 # input checking 1159 struct <- match.arg(struct) 1160 if (NCOL(X) < 2) stop("R must have at least 2 variables") 1161 if ((struct == "observedfactor") && is.null(f)) stop("Provide the factor observations") 1162 if (unbiasedMarg && !((struct == "Indep") || (struct == "IndepId"))) stop("unbiasedMarg can only be combined with T2, T3 and T6") 1163 1164 # compute useful variables 1165 NN <- dim(X)[1] # number of observations 1166 if (unbiasedMarg) { 1167 if (NN < 6) stop("R should have at least 6 observations") 1168 } 1169 PP <- dim(X)[2] # number of assets 1170 ncosk <- PP * (PP + 1) * (PP + 2) / 6 # number of unique coskewness elements 1171 1172 Xc <- X - matrix(colMeans(X), nrow = NN, ncol = PP, byrow = TRUE) # center the observations 1173 margskews <- colMeans(Xc^3) 1174 1175 # compute the coskewness matrix 1176 if (struct == "latent1factor") { 1177 # 1-factor model of Simaan (1993) 1178 margskewsroot <- sign(margskews) * abs(margskews)^(1 / 3) 1179 T3 <- .Call('M3_Simaan', margskewsroot, PP, PACKAGE="PerformanceAnalytics") 1180 if (as.mat) T3 <- M3.vec2mat(T3, PP) 1181 return (T3) 1182 1183 } else if (struct == "Indep") { 1184 # independent marginals 1185 if (unbiasedMarg) margskews <- margskews * NN^2 / ((NN - 1) * (NN - 2)) 1186 T3 <- .Call('M3_T23', margskews, PP, PACKAGE="PerformanceAnalytics") 1187 if (as.mat) T3 <- M3.vec2mat(T3, PP) 1188 return (T3) 1189 1190 } else if (struct == "IndepId") { 1191 # independent and equally distributed marginals 1192 if (unbiasedMarg) margskews <- margskews * NN^2 / ((NN - 1) * (NN - 2)) 1193 margskews <- rep(mean(margskews), PP) 1194 T3 <- .Call('M3_T23', margskews, PP, PACKAGE="PerformanceAnalytics") 1195 if (as.mat) T3 <- M3.vec2mat(T3, PP) 1196 return (T3) 1197 1198 } else if (struct == "observedfactor") { 1199 # observed factor model 1200 if (NCOL(f) == 1) { 1201 beta <- apply(Xc, 2, function(a) cov(a, f) / var(f)) 1202 fc <- f - mean(f) 1203 fskew <- mean(fc^3) 1204 T3 <- .Call('M3_1F', margskews, beta, fskew, PP, PACKAGE="PerformanceAnalytics") 1205 } else { 1206 mod <- stats::lm(X ~ f) 1207 beta <- t(mod$coefficients[-1,]) 1208 M3_factor <- M3.MM(f, as.mat = FALSE) 1209 M3_factor <- .Call('M3timesFull', M3_factor, beta, NCOL(beta), PP, PACKAGE="PerformanceAnalytics") 1210 epsskews <- colMeans(mod$residuals^3) 1211 T3 <- M3_factor + .Call('M3_T23', epsskews, PP, PACKAGE="PerformanceAnalytics") 1212 } 1213 if (as.mat) T3 <- M3.vec2mat(T3, PP) 1214 return (T3) 1215 1216 } else if (struct == "CC") { 1217 # constant-correlation (Martellini and Ziemann (2010)) (symmetrized version) 1218 margvars <- colMeans(Xc^2) 1219 margkurts <- colMeans(Xc^4) 1220 m21 <- as.numeric(t(Xc^2) %*% Xc) / NN 1221 m22 <- as.numeric(t(Xc^2) %*% Xc^2) / NN 1222 r_generalized <- .Call('M3_CCoefficients', margvars, margkurts, m21, 1223 m22, as.numeric(Xc), NN, PP, PACKAGE="PerformanceAnalytics") 1224 T3 <- .Call('M3_CC', margvars, margskews, margkurts, r_generalized[1], r_generalized[2], 1225 r_generalized[3], PP, PACKAGE="PerformanceAnalytics") 1226 if (as.mat) T3 <- M3.vec2mat(T3, PP) 1227 return (T3) 1228 1229 } else if (struct == "CS") { 1230 # coskewness matrix under central-symmetry 1231 return (rep(0, ncosk)) 1232 1233 } 1234} 1235 1236 1237#'@useDynLib PerformanceAnalytics 1238#'@export 1239#'@rdname StructuredMoments 1240M4.struct <- function(R, struct = c("Indep", "IndepId", "observedfactor", "CC"), f = NULL, as.mat = TRUE) { 1241 # @author Dries Cornilly 1242 # 1243 # DESCRIPTION: 1244 # computes different strutured cokurtosis estimators 1245 # 1246 # Inputs: 1247 # R : numeric matrix of dimensions NN x PP 1248 # struct : select the structured estimator 1249 # : Indep : independent, unequal marginals 1250 # : IndepId : independent, equal marginals 1251 # : observedfactor : linear factor model with observed factors 1252 # : CC : constant-correlation model of Martellini and Ziemann (2010), symmetrized 1253 # f : numeric vector with factor observations 1254 # 1255 # Outputs: 1256 # cokurtosis matrix (as matrix or as vector, depending on as.mat) 1257 1258 X <- coredata(R) 1259 1260 # input checking 1261 struct <- match.arg(struct) 1262 if (NCOL(X) < 2) stop("R must have at least 2 variables") 1263 if ((struct == "observedfactor") && is.null(f)) stop("Provide the factor observations") 1264 1265 # compute useful variables 1266 NN <- dim(X)[1] # number of observations 1267 PP <- dim(X)[2] # number of assets 1268 1269 Xc <- X - matrix(colMeans(X), nrow = NN, ncol = PP, byrow = TRUE) # center the observations 1270 margkurts <- colMeans(Xc^4) 1271 margvars <- colMeans(Xc^2) 1272 1273 # compute the coskewness matrix 1274 if (struct == "Indep") { 1275 # independent marginals 1276 T4 <- .Call('M4_T12', margkurts, margvars, PP, PACKAGE="PerformanceAnalytics") 1277 if (as.mat) T4 <- M4.vec2mat(T4, PP) 1278 return (T4) 1279 1280 } else if (struct == "IndepId") { 1281 # independent and equally distributed marginals 1282 meanmargkurts <- mean(margkurts) 1283 meank_iikk <- sqrt(mean(margvars^2)) 1284 T4 <- .Call('M4_T12', rep(meanmargkurts, PP), rep(meank_iikk, PP), PP, PACKAGE="PerformanceAnalytics") 1285 if (as.mat) T4 <- M4.vec2mat(T4, PP) 1286 return (T4) 1287 1288 } else if (struct == "observedfactor") { 1289 # linear factor model with observed factors 1290 if (NCOL(f) == 1) { 1291 beta <- apply(Xc, 2, function(a) cov(a, f) / var(f)) 1292 fc <- f - mean(f) 1293 fvar <- mean(fc^2) 1294 fkurt <- mean(fc^4) 1295 epsvars <- margvars - beta^2 * fvar 1296 T4 <- .Call('M4_1f', margkurts, fvar, fkurt, epsvars, beta, PP, PACKAGE="PerformanceAnalytics") 1297 } else { 1298 mod <- stats::lm(X ~ f) 1299 beta <- t(mod$coefficients[-1,]) 1300 M4_factor <- M4.MM(f, as.mat = FALSE) 1301 M4_factor <- .Call('M4timesFull', M4_factor, beta, NCOL(beta), PP, PACKAGE="PerformanceAnalytics") 1302 epskurts <- colMeans(mod$residuals^4) 1303 epsvars <- colMeans(mod$residuals^2) 1304 Stransf <- beta %*% cov(f) %*% t(beta) * (NN - 1) / NN 1305 M4_factor <- M4_factor + .Call('M4_MFresid', as.numeric(Stransf), epsvars, PP, PACKAGE="PerformanceAnalytics") 1306 T4 <- M4_factor + .Call('M4_T12', epskurts, epsvars, PP, PACKAGE="PerformanceAnalytics") 1307 } 1308 if (as.mat) T4 <- M4.vec2mat(T4, PP) 1309 return (T4) 1310 1311 } else if (struct == "CC") { 1312 # constant-correlation (Martellini and Ziemann (2010)) 1313 marg6s <- colMeans(Xc^6) 1314 m22 <- as.numeric(t(Xc^2) %*% Xc^2) / NN 1315 m31 <- as.numeric(t(Xc^3) %*% Xc) / NN 1316 r_generalized <- .Call('M4_CCoefficients', margvars, margkurts, marg6s, 1317 m22, m31, as.numeric(Xc), NN, PP, PACKAGE="PerformanceAnalytics") 1318 T4 <- .Call('M4_CC', margvars, margkurts, marg6s, r_generalized[1], r_generalized[2], 1319 r_generalized[3], r_generalized[4], PP, PACKAGE="PerformanceAnalytics") 1320 if (as.mat) T4 <- M4.vec2mat(T4, PP) 1321 return (T4) 1322 1323 } 1324} 1325 1326 1327#' Functions for calculating EWMA comoments of financial time series 1328#' 1329#' calculates exponentially weighted moving average covariance, coskewness and cokurtosis matrices 1330#' 1331#' The coskewness and cokurtosis matrices are defined as the matrices of dimension 1332#' p x p^2 and p x p^3 containing the third and fourth order central moments. They 1333#' are useful for measuring nonlinear dependence between different assets of the 1334#' portfolio and computing modified VaR and modified ES of a portfolio. 1335#' 1336#' EWMA estimation of the covariance matrix was popularized by the RiskMetrics report in 1996. 1337#' The M3.ewma and M4.ewma are straightforward extensions to the setting of third and fourth 1338#' order central moments 1339#' @name EWMAMoments 1340#' @concept co-moments 1341#' @concept moments 1342#' @aliases EWMAMoments M2.ewma M3.ewma M4.ewma 1343#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of 1344#' asset returns (with mean zero) 1345#' @param lambda decay coefficient 1346#' @param last.M2 last estimated covariance matrix before the observed returns R 1347#' @param last.M3 last estimated coskewness matrix before the observed returns R 1348#' @param last.M4 last estimated cokurtosis matrix before the observed returns R 1349#' @param as.mat TRUE/FALSE whether to return the full moment matrix or only 1350#' the vector with the unique elements (the latter is advised for speed), default 1351#' TRUE 1352#' @param \dots any other passthru parameters 1353#' @author Dries Cornilly 1354#' @seealso \code{\link{CoMoments}} \cr \code{\link{ShrinkageMoments}} \cr \code{\link{StructuredMoments}} \cr \code{\link{MCA}} 1355#' @references 1356#' JP Morgan. Riskmetrics technical document. 1996. 1357###keywords ts multivariate distribution models 1358#' @examples 1359#' 1360#' data(edhec) 1361#' 1362#' # EWMA estimation 1363#' # 'as.mat = F' would speed up calculations in higher dimensions 1364#' sigma <- M2.ewma(edhec, 0.94) 1365#' m3 <- M3.ewma(edhec, 0.94) 1366#' m4 <- M4.ewma(edhec, 0.94) 1367#' 1368#' # compute equal-weighted portfolio modified ES 1369#' mu <- colMeans(edhec) 1370#' p <- length(mu) 1371#' ES(p = 0.95, portfolio_method = "component", weights = rep(1 / p, p), mu = mu, 1372#' sigma = sigma, m3 = m3, m4 = m4) 1373#' 1374#' # compare to sample method 1375#' sigma <- cov(edhec) 1376#' m3 <- M3.MM(edhec) 1377#' m4 <- M4.MM(edhec) 1378#' ES(p = 0.95, portfolio_method = "component", weights = rep(1 / p, p), mu = mu, 1379#' sigma = sigma, m3 = m3, m4 = m4) 1380#' 1381#' @export M2.ewma 1382M2.ewma <- function(R, lambda = 0.97, last.M2 = NULL, ...) { 1383 # R : numeric matrix of dimensions NN x PP; top of R are oldest observations, bottom are the newest 1384 # : assumes a mean of zero 1385 # lambda : decay parameter for the correlations 1386 if(hasArg(lambda_var)) lambda_var <- list(...)$lambda_var else lambda_var <- NULL 1387 1388 x <- coredata(R) 1389 NN <- NROW(x) 1390 PP <- NCOL(x) 1391 1392 if (is.null(last.M2)) { 1393 if (lambda > 1 - 1e-07) { 1394 lambda_vec <- rep(1 / NN, NN) 1395 } else { 1396 lambda_vec <- (1 - lambda) / (1 - lambda^NN) * lambda^((NN - 1):0) 1397 } 1398 Xc <- x * (matrix(lambda_vec, nrow = NN, ncol = PP, byrow = FALSE))^(1 / 2) 1399 M2 <- t(Xc) %*% Xc 1400 1401 if (!is.null(lambda_var)) { 1402 sd_lambda <- sqrt(diag(M2)) 1403 R2 <- diag(1 / sd_lambda) %*% M2 %*% diag(1 / sd_lambda) 1404 1405 if (lambda_var > 1 - 1e-07) { 1406 lambda_var_vec <- rep(1 / NN, NN) 1407 } else { 1408 lambda_var_vec <- (1 - lambda_var) / (1 - lambda_var^NN) * lambda_var^((NN - 1):0) 1409 } 1410 sd_lambda_var <- sqrt(colSums(x^2 * matrix(lambda_var_vec, nrow = NN, ncol = PP, byrow = FALSE))) 1411 M2 <- diag(sd_lambda_var) %*% R2 %*% diag(sd_lambda_var) 1412 } 1413 } else { 1414 if (PP == 1) { 1415 x <- matrix(x, nrow = 1) 1416 NN <- 1 1417 } 1418 M2 <- last.M2 1419 for (tt in 1:NN) { 1420 xn <- matrix(x[tt,], nrow = 1) 1421 new.M2 <- t(xn) %*% xn 1422 M2 <- (1 - lambda) * new.M2 + lambda * M2 1423 } 1424 } 1425 1426 return( M2 ) 1427} 1428 1429#'@useDynLib PerformanceAnalytics 1430#'@export 1431#'@rdname EWMAMoments 1432M3.ewma <- function(R, lambda = 0.97, last.M3 = NULL, as.mat = TRUE, ...) { 1433 # R : numeric matrix of dimensions NN x PP; top of R are oldest observations, bottom are the newest 1434 # : assumes a mean of zero 1435 # lambda : decay parameter for the standardized coskewnesses 1436 # lambda_var : if not NULL, use lambda_var for the variances 1437 # as.mat : TRUE/FALSE whether or not the output is a matrix or not 1438 if(hasArg(lambda_var)) lambda_var <- list(...)$lambda_var else lambda_var <- NULL 1439 1440 x <- coredata(R) 1441 NN <- NROW(x) 1442 PP <- NCOL(x) 1443 1444 if (is.null(last.M3)) { 1445 if (lambda > 1 - 1e-07) { 1446 lambda_vec <- rep(1 / NN, NN) 1447 } else { 1448 lambda_vec <- (1 - lambda) / (1 - lambda^NN) * lambda^((NN - 1):0) 1449 } 1450 Xc <- x * (matrix(lambda_vec, nrow = NN, ncol = PP, byrow = FALSE))^(1 / 3) 1451 M3 <- .Call('M3sample', as.numeric(Xc), NN, PP, 1.0, PACKAGE="PerformanceAnalytics") 1452 1453 if (!is.null(lambda_var)) { 1454 sd_lambda <- sqrt(colSums(x^2 * matrix(lambda_vec, nrow = NN, ncol = PP, byrow = FALSE))) 1455 R3 <- .Call('M3timesDiag', M3, 1 / sd_lambda, PP, PACKAGE="PerformanceAnalytics") 1456 1457 if (lambda_var > 1 - 1e-07) { 1458 lambda_var_vec <- rep(1 / NN, NN) 1459 } else { 1460 lambda_var_vec <- (1 - lambda_var) / (1 - lambda_var^NN) * lambda_var^((NN - 1):0) 1461 } 1462 sd_lambda_var <- sqrt(colSums(x^2 * matrix(lambda_var_vec, nrow = NN, ncol = PP, byrow = FALSE))) 1463 1464 M3 <- .Call('M3timesDiag', R3, sd_lambda_var, PP, PACKAGE="PerformanceAnalytics") 1465 } 1466 } else { 1467 if (PP == 1) { 1468 x <- matrix(x, nrow = 1) 1469 NN <- 1 1470 PP <- ncol(x) 1471 } 1472 if (NROW(last.M3) == PP) last.M3 <- M3.mat2vec(last.M3) 1473 M3 <- last.M3 1474 for (tt in 1:NN) { 1475 xn <- matrix(x[tt,], nrow = 1) 1476 new.M3 <- M3.MM(xn, as.mat = FALSE, mu = rep(0, PP)) 1477 M3 <- (1 - lambda) * new.M3 + lambda * M3 1478 } 1479 } 1480 1481 if (as.mat) M3 <- M3.vec2mat(M3, PP) 1482 1483 return( M3 ) 1484} 1485 1486#'@useDynLib PerformanceAnalytics 1487#'@export 1488#'@rdname EWMAMoments 1489M4.ewma <- function(R, lambda = 0.97, last.M4 = NULL, as.mat = TRUE, ...) { 1490 # R : numeric matrix of dimensions NN x PP; top of R are oldest observations, bottom are the newest 1491 # : assumes a mean of zero 1492 # lambda : decay parameter for the standardized cokurtosisses 1493 # lambda_var : if not NULL, use lambda_var for the variances 1494 # as.mat : TRUE/FALSE whether or not the output is a matrix or not 1495 if(hasArg(lambda_var)) lambda_var <- list(...)$lambda_var else lambda_var <- NULL 1496 1497 x <- coredata(R) 1498 NN <- NROW(x) 1499 PP <- NCOL(x) 1500 1501 if (is.null(last.M4)) { 1502 if (lambda > 1 - 1e-07) { 1503 lambda_vec <- rep(1 / NN, NN) 1504 } else { 1505 lambda_vec <- (1 - lambda) / (1 - lambda^NN) * lambda^((NN - 1):0) 1506 } 1507 Xc <- x * (matrix(lambda_vec, nrow = NN, ncol = PP, byrow = FALSE))^(1 / 4) 1508 M4 <- .Call('M4sample', as.numeric(Xc), NN, PP, PACKAGE="PerformanceAnalytics") * NN 1509 1510 if (!is.null(lambda_var)) { 1511 sd_lambda <- sqrt(colSums(x^2 * matrix(lambda_vec, nrow = NN, ncol = PP, byrow = FALSE))) 1512 R4 <- .Call('M4timesDiag', M4, 1 / sd_lambda, PP, PACKAGE="PerformanceAnalytics") 1513 1514 if (lambda_var > 1 - 1e-07) { 1515 lambda_var_vec <- rep(1 / NN, NN) 1516 } else { 1517 lambda_var_vec <- (1 - lambda_var) / (1 - lambda_var^NN) * lambda_var^((NN - 1):0) 1518 } 1519 sd_lambda_var <- sqrt(colSums(x^2 * matrix(lambda_var_vec, nrow = NN, ncol = PP, byrow = FALSE))) 1520 1521 M4 <- .Call('M4timesDiag', R4, sd_lambda_var, PP, PACKAGE="PerformanceAnalytics") 1522 } 1523 } else { 1524 if (PP == 1) { 1525 x <- matrix(x, nrow = 1) 1526 NN <- 1 1527 PP <- ncol(x) 1528 } 1529 if (NROW(last.M4) == PP) last.M4 <- M4.mat2vec(last.M4) 1530 M4 <- last.M4 1531 for (tt in 1:NN) { 1532 xn <- matrix(x[tt,], nrow = 1) 1533 new.M4 <- M4.MM(xn, as.mat = FALSE, mu = rep(0, PP)) 1534 M4 <- (1 - lambda) * new.M4 + lambda * M4 1535 } 1536 } 1537 1538 if (as.mat) M4 <- M4.vec2mat(M4, PP) 1539 1540 return( M4 ) 1541} 1542 1543 1544#' Functions for doing Moment Component Analysis (MCA) of financial time series 1545#' 1546#' calculates MCA coskewness and cokurtosis matrices 1547#' 1548#' The coskewness and cokurtosis matrices are defined as the matrices of dimension 1549#' p x p^2 and p x p^3 containing the third and fourth order central moments. They 1550#' are useful for measuring nonlinear dependence between different assets of the 1551#' portfolio and computing modified VaR and modified ES of a portfolio. 1552#' 1553#' MCA is a generalization of PCA to higher moments. The principal components in 1554#' MCA are the ones that maximize the coskewness and cokurtosis present when projecting 1555#' onto these directions. It was introduced by Lim and Morton (2007) and applied to financial returns 1556#' data by Jondeau and Rockinger (2017) 1557#' 1558#' If a coskewness matrix (argument M3) or cokurtosis matrix (argument M4) is passed in using ..., then 1559#' MCA is performed on the given comoment matrix instead of the sample coskewness or cokurtosis matrix. 1560#' @name MCA 1561#' @concept co-moments 1562#' @concept moments 1563#' @aliases MCA M3.MCA M4.MCA 1564#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of 1565#' asset returns 1566#' @param k the number of components to use 1567#' @param as.mat TRUE/FALSE whether to return the full moment matrix or only 1568#' the vector with the unique elements (the latter is advised for speed), default 1569#' TRUE 1570#' @param \dots any other passthru parameters 1571#' @author Dries Cornilly 1572#' @seealso \code{\link{CoMoments}} \cr \code{\link{ShrinkageMoments}} \cr \code{\link{StructuredMoments}} \cr \code{\link{EWMAMoments}} 1573#' @references 1574#' Lim, Hek-Leng and Morton, Jason. 2007. Principal Cumulant Component Analysis. working paper 1575#' 1576#' Jondeau, Eric and Jurczenko, Emmanuel. 2017. Moment Component Analysis: An Illustration 1577#' With International Stock Markets. Journal of Business and Economic Statistics 1578#' @examples 1579#' data(edhec) 1580#' 1581#' # coskewness matrix based on two components 1582#' M3mca <- M3.MCA(edhec, k = 2)$M3mca 1583#' 1584#' # screeplot MCA 1585#' M3dist <- M4dist <- rep(NA, ncol(edhec)) 1586#' M3S <- M3.MM(edhec) # sample coskewness estimator 1587#' M4S <- M4.MM(edhec) # sample cokurtosis estimator 1588#' for (k in 1:ncol(edhec)) { 1589#' M3MCA_list <- M3.MCA(edhec, k) 1590#' M4MCA_list <- M4.MCA(edhec, k) 1591#' 1592#' M3dist[k] <- sqrt(sum((M3S - M3MCA_list$M3mca)^2)) 1593#' M4dist[k] <- sqrt(sum((M4S - M4MCA_list$M4mca)^2)) 1594#' } 1595#' par(mfrow = c(2, 1)) 1596#' plot(1:ncol(edhec), M3dist) 1597#' plot(1:ncol(edhec), M4dist) 1598#' par(mfrow = c(1, 1)) 1599#' 1600#' @export M3.MCA 1601M3.MCA <- function(R, k = 1, as.mat = TRUE, ...) { 1602 # @author Dries Cornilly 1603 # 1604 # DESCRIPTION: 1605 # MCA on coskewness matrix 1606 # 1607 # Inputs: 1608 # R : numeric matrix of dimensions NN x PP 1609 # k : number of components to use 1610 # as.mat : TRUE/FALSE whether or not the output is a matrix or not 1611 # 1612 # Outputs: 1613 # M3mca : coskewness matrix based on k factors 1614 # converged : logical indicating convergence 1615 # iter : number of iterations at convergence / at end 1616 # U : matrix with principal components 1617 1618 x <- coredata(R) 1619 1620 if (hasArg(maxit)) maxit <- list(...)$maxit else maxit <- 1000 1621 if (hasArg(abstol)) abstol <- list(...)$abstol else abstol <- 1e-05 1622 1623 p <- NCOL(x) 1624 n <- NROW(x) 1625 1626 # initialize projection matrix and sample coskewness matrix 1627 if (hasArg(M3)) { 1628 M3 <- list(...)$M3 1629 if (NROW(M3) == p) M3 <- M3.mat2vec(M3) 1630 } else { 1631 M3 <- M3.MM(x, as.mat = FALSE) 1632 } 1633 if (hasArg(U0)) U0 <- list(...)$U0 else U0 <- svd(M3.vec2mat(M3, p), nu = k, nv = 0)$u 1634 1635 # iterate until convergence or maximum number of iterations is reached 1636 iter <- 0 1637 converged <- FALSE 1638 while ((iter < maxit) && !converged) { 1639 # project using the last projection matrix U0 and build new projection matrix U1 1640 Z3 <- .Call('M3HOOIiterator', M3, as.numeric(U0), p, k, PACKAGE="PerformanceAnalytics") 1641 U1 <- svd(Z3, nu = k, nv = 0)$u 1642 1643 # check for convergence - absolute values are because the direction of the eigenvectors might alternate 1644 if (sqrt(sum((abs(U1) - abs(U0))^2)) < abstol) converged <- TRUE 1645 1646 # update U and the iterator 1647 U0 <- U1 1648 iter <- iter + 1 1649 } 1650 1651 # build estimated coskewness matrix from the projection matrix U0 1652 C3 <- .Call('M3timesFull', M3, as.numeric(t(U0)), p, k, PACKAGE="PerformanceAnalytics") 1653 M3mca <- .Call('M3timesFull', C3, as.numeric(U0), k, p, PACKAGE="PerformanceAnalytics") 1654 1655 if (as.mat) M3mca <- M3.vec2mat(M3mca, p) 1656 1657 return ( list("M3mca" = M3mca, "converged" = converged, "iter" = iter, "U" = U0) ) 1658} 1659 1660 1661#'@useDynLib PerformanceAnalytics 1662#'@export 1663#'@rdname MCA 1664M4.MCA <- function(R, k = 1, as.mat = TRUE, ...) { 1665 # @author Dries Cornilly 1666 # 1667 # DESCRIPTION: 1668 # MCA on cokurtosis matrix 1669 # 1670 # Inputs: 1671 # R : numeric matrix of dimensions NN x PP 1672 # k : number of components to use 1673 # as.mat : TRUE/FALSE whether or not the output is a matrix or not 1674 # 1675 # Outputs: 1676 # M4mca : cokurtosis matrix based on k factors 1677 # converged : logical indicating convergence 1678 # iter : number of iterations at convergence / at end 1679 # U : matrix with principal components 1680 1681 x <- coredata(R) 1682 1683 if (hasArg(maxit)) maxit <- list(...)$maxit else maxit <- 1000 1684 if (hasArg(abstol)) abstol <- list(...)$abstol else abstol <- 1e-05 1685 1686 p <- NCOL(x) 1687 n <- NROW(x) 1688 1689 # initialize projection matrix and sample cokurtosis matrix 1690 if (hasArg(M4)) { 1691 M4 <- list(...)$M4 1692 if (NROW(M4) == p) M4 <- M4.mat2vec(M4) 1693 } else { 1694 M4 <- M4.MM(x, as.mat = FALSE) 1695 } 1696 if (hasArg(U0)) U0 <- list(...)$U0 else U0 <- svd(M4.vec2mat(M4, p), nu = k, nv = 0)$u 1697 1698 # iterate until convergence or maximum number of iterations is reached 1699 iter <- 0 1700 converged <- FALSE 1701 while ((iter < maxit) && !converged) { 1702 # project using the last projection matrix U0 and build new projection matrix U1 1703 Z4 <- .Call('M4HOOIiterator', M4, as.numeric(U0), p, k, PACKAGE="PerformanceAnalytics") 1704 U1 <- svd(Z4, nu = k, nv = 0)$u 1705 1706 # check for convergence - absolute values are because the direction of the eigenvectors might alternate 1707 if (sqrt(sum((abs(U1) - abs(U0))^2)) < abstol) converged <- TRUE 1708 1709 # update U and the iterator 1710 U0 <- U1 1711 iter <- iter + 1 1712 } 1713 1714 # build estimated coskewness matrix from the projection matrix U0 1715 C4 <- .Call('M4timesFull', M4, as.numeric(t(U0)), p, k, PACKAGE="PerformanceAnalytics") 1716 M4mca <- .Call('M4timesFull', C4, as.numeric(U0), k, p, PACKAGE="PerformanceAnalytics") 1717 1718 if (as.mat) M4mca <- M4.vec2mat(M4mca, p) 1719 1720 return ( list("M4mca" = M4mca, "converged" = converged, "iter" = iter, "U" = U0) ) 1721} 1722 1723 1724############################################################################### 1725# R (http://r-project.org/) Econometrics for Performance and Risk Analysis 1726# 1727# Copyright (c) 2004-2020 Peter Carl and Brian G. Peterson and Kris Boudt 1728# 1729# This R package is distributed under the terms of the GNU Public License (GPL) 1730# for full details see the file COPYING 1731# 1732# $Id$ 1733# 1734############################################################################### 1735