通过透视Cholesky分解生成具有秩不足协方差的多元正态r.v. [英] Generate multivariate normal r.v.'s with rank-deficient covariance via Pivoted Cholesky Factorization
问题描述
我只是在脑子里跳动，试图使Cholesky分解起作用，以便模拟相关的价格变动.
I'm just beating my head against the wall trying to get a Cholesky decomposition to work in order to simulate correlated price movements.
我使用以下代码:
cormat <- as.matrix(read.csv("http://pastebin.com/raw/qGbkfiyA"))
cormat <- cormat[,2:ncol(cormat)]
rownames(cormat) <- colnames(cormat)
cormat <- apply(cormat,c(1,2),FUN = function(x) as.numeric(x))
chol(cormat)
#Error in chol.default(cormat) :
# the leading minor of order 8 is not positive definite
cholmat <- chol(cormat, pivot=TRUE)
#Warning message:
# In chol.default(cormat, pivot = TRUE) :
# the matrix is either rank-deficient or indefinite
rands <- array(rnorm(ncol(cholmat)), dim = c(10000,ncol(cholmat)))
V <- t(t(cholmat) %*% t(rands))
#Check for similarity
cor(V) - cormat ## Not all zeros!
#Check the standard deviations
apply(V,2,sd) ## Not all ones!
我不太确定如何正确使用pivot = TRUE语句来生成相关的动作.结果看起来完全是伪造的.
I'm not really sure how to properly use the pivot = TRUE statement to generate my correlated movements. The results look totally bogus.
即使我有一个简单的矩阵并且尝试枢轴"操作，我也会得到虚假结果...
Even if I have a simple matrix and I try out "pivot" then I get bogus results...
cormat <- matrix(c(1,.95,.90,.95,1,.93,.90,.93,1), ncol=3)
cholmat <- chol(cormat)
# No Error
cholmat2 <- chol(cormat, pivot=TRUE)
# No warning... pivot changes column order
rands <- array(rnorm(ncol(cholmat)), dim = c(10000,ncol(cholmat)))
V <- t(t(cholmat2) %*% t(rands))
#Check for similarity
cor(V) - cormat ## Not all zeros!
#Check the standard deviations
apply(V,2,sd) ## Not all ones!
推荐答案
您的代码有两个错误:
-
您没有使用数据透视索引将数据透视还原为Cholesky因子.请注意，半正定矩阵
A的枢轴Cholesky因式分解正在做:
You did not use pivoting index to revert the pivoting done to the Cholesky factor. Note, pivoted Cholesky factorization for a semi-positive definite matrix
Ais doing:
P'AP = R'R
其中，P是列枢轴矩阵，而R是上三角矩阵.要从R恢复A，我们需要应用P的倒数(即P'):
where P is a column pivoting matrix, and R is an upper triangular matrix. To recover A from R, we need apply the inverse of P (i.e., P'):
A = PR'RP' = (RP')'(RP')
协方差矩阵为A的多元法线由以下方式生成:
Multivariate normal with covariance matrix A, is generated by:
XRP'
其中X是具有零均值和同一性协方差的多元正态.
where X is multivariate normal with zero mean and identity covariance.
您生成的X
X <- array(rnorm(ncol(R)), dim = c(10000,ncol(R)))
是错误的.首先，它不应该是ncol(R)，而应该是nrow(R)，即X的等级，由r表示.其次，沿着列循环使用rnorm(ncol(R))，结果矩阵根本不是随机的.因此，cor(X)永远不会接近单位矩阵.正确的代码是:
is wrong. First, it should not be ncol(R) but nrow(R), i.e., the rank of X, denoted by r. Second, you are recycling rnorm(ncol(R)) along columns, and the resulting matrix is not random at all. Therefore, cor(X) is never close to an identity matrix. The correct code is:
X <- matrix(rnorm(10000 * r), 10000, r)
作为上述理论的模型实现，请考虑您的玩具示例:
As a model implementation of the above theory, consider your toy example:
A <- matrix(c(1,.95,.90,.95,1,.93,.90,.93,1), ncol=3)
我们计算较高的三角因子(抑制可能的等级不足警告)，并提取逆向枢轴指标和等级:
We compute the upper triangular factor (suppressing possible rank-deficient warnings) and extract inverse pivoting index and rank:
R <- suppressWarnings(chol(A, pivot = TRUE))
piv <- order(attr(R, "pivot")) ## reverse pivoting index
r <- attr(R, "rank") ## numerical rank
然后我们生成X.为了获得更好的结果，我们将X居中，以使列均值为0.
Then we generate X. For better result we centre X so that column means are 0.
X <- matrix(rnorm(10000 * r), 10000, r)
## for best effect, we centre `X`
X <- sweep(X, 2L, colMeans(X), "-")
然后我们生成目标多元正态:
Then we generate target multivariate normal:
## compute `V = RP'`
V <- R[1:r, piv]
## compute `Y = X %*% V`
Y <- X %*% V
我们可以验证Y具有目标协方差A:
We can verify that Y has target covariance A:
cor(Y)
# [,1] [,2] [,3]
#[1,] 1.0000000 0.9509181 0.9009645
#[2,] 0.9509181 1.0000000 0.9299037
#[3,] 0.9009645 0.9299037 1.0000000
A
# [,1] [,2] [,3]
#[1,] 1.00 0.95 0.90
#[2,] 0.95 1.00 0.93
#[3,] 0.90 0.93 1.00
这篇关于通过透视Cholesky分解生成具有秩不足协方差的多元正态r.v.的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
