Генерация многомерных нормальных rv с ранг-дефицитной ковариацией с помощью Pivoted Cholesky Factorization

Question

Я просто бьюсь головой об стену, пытаясь заставить разложение Холецкого работать, чтобы имитировать коррелированные ценовые движения.

Я использую следующий код:

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

В вашем коде есть две ошибки:

1. Вы не использовали индекс поворота, чтобы вернуть значение поворота к коэффициенту Холецкого. Обратите внимание, развернутая холеристская факторизация для полуположительной определенной матрицы A делается:

   P'AP = R'R
   

   где P является матрицей поворота столбца, и R верхняя треугольная матрица. Восстановить A от Rнам нужно применить обратное P (То есть, P'):

   A = PR'RP' = (RP')'(RP')
   

   Многомерный нормаль с ковариационной матрицей A, генерируется:

   XRP'
   

   где X многомерная нормальная с нулевым средним и единичной ковариацией.

2. Ваше поколение 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