通过 Pivoted Cholesky 分解生成具有秩亏协方差的多元正态 r.v.'s

Question

我只是在用头撞墙，试图让 Cholesky 分解起作用，以模拟相关的价格变动。

我使用以下代码：

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 语句来生成我的相关运动。结果看起来完全是假的。

即使我有一个简单的矩阵并尝试 "pivot"，我也会得到错误的结果...

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!

Answer 1

您的代码有两个错误：

您没有使用旋转索引将完成的旋转恢复到 Cholesky 因子。请注意，半正定矩阵 A 的旋转 Cholesky 分解正在做：
```
P'AP = R'R
```
其中P是列旋转矩阵，R是上三角矩阵。要从 R 恢复 A，我们需要应用 P 的逆（即 P'）：
```
A = PR'RP' = (RP')'(RP')
```
具有协方差矩阵 A 的多变量正态分布由：
生成
```
XRP'
```
其中 X 是具有零均值和恒等协方差的多元正态分布。
贵辈X
```
X <- array(rnorm(ncol(R)), dim = c(10000,ncol(R)))
```
错了。首先，它不应该是ncol(R)而是nrow(R)，即X的等级，记为r。其次，您正在沿列回收 rnorm(ncol(R))，并且生成的矩阵根本不是随机的。因此，cor(X) 永远不会接近单位矩阵。正确的代码是：
```
X <- matrix(rnorm(10000 * r), 10000, r)
```

作为上述理论的模型实现，请考虑您的玩具示例：

A <- matrix(c(1,.95,.90,.95,1,.93,.90,.93,1), ncol=3)

我们计算上三角因子（抑制可能的排名不足警告）并提取逆枢轴指数和排名：

R <- suppressWarnings(chol(A, pivot = TRUE))
piv <- order(attr(R, "pivot"))  ## reverse pivoting index
r <- attr(R, "rank")  ## numerical rank

然后我们生成X。为了获得更好的结果，我们将 X 居中，以便列平均值为 0.

X <- matrix(rnorm(10000 * r), 10000, r)
## for best effect, we centre `X`
X <- sweep(X, 2L, colMeans(X), "-")

然后我们生成目标多元正态：

## compute `V = RP'`
V <- R[1:r, piv]

## compute `Y = X %*% V`
Y <- X %*% V

我们可以验证 Y 具有目标协方差 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 分解生成具有秩亏协方差的多元正态 r.v.'s

Generate multivariate normal r.v.'s with rank-deficient covariance via Pivoted Cholesky Factorization

r

matrix

random

statistics

normal-distribution