通过 Pivoted Cholesky Factorization 生成具有秩亏协方差的多元正态 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.
即使我有一个简单的矩阵并且我尝试了pivot"，我也会得到虚假的结果...
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
这篇关于通过 Pivoted Cholesky Factorization 生成具有秩亏协方差的多元正态 r.v.的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
