Сгенерируйте многомерные нормальные с.в. с ковариацией с недостаточным рангом с помощью развернутой факторизации Холецкого

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

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

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 для генерации моих коррелированных движений. Результаты выглядят совершенно фальшивыми.

Даже если у меня есть простая матрица, и я попробую "разворот", я получу фальшивые результаты ...

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!

r random matrix statistics normal-distribution
person JoeBass    schedule 11.02.2016    source источник
comment
Мое первое предположение состоит в том, что cormat действительно p.d. но имеет некоторые собственные значения, которые близки к нулю, что вызывает числовые затруднения. Вы можете доказать или опровергнуть это, вычислив собственные значения. Кстати, а каково происхождение cormat? Если у вас есть некоторый контроль над этим, можете ли вы гарантировать, что cormat будет более определенным? (например, добавьте постоянный множитель к диагонали или сконструируйте его таким образом, чтобы гарантировать p.d.-ness.)   -  person Robert Dodier    schedule 11.02.2016
comment
Некоторые собственные значения действительно отрицательны. Я подозреваю, что это связано с тем, что цены, из которых были взяты корреляции, были сделаны в разные периоды. Например, S1 ‹-› S2 взят из другого периода, как S2 ‹-› S3. Я могу модифицировать кормат по своему усмотрению. Я попытался обнулить отрицательные собственные значения, но безуспешно использовал это в качестве руководства. risklatte.com/Articles/QuantitativeFinance/QF146.php   -  person JoeBass    schedule 12.02.2016
comment
Мой совет на этом этапе - построить меньшую матрицу и убедиться, что ваш метод работает. Я понял, что элементы в матрице являются временными корреляциями, при этом S (i, j) = (корреляция j - i временных шагов друг от друга). Если да, может быть, все элементы на данной поддиагонали должны быть равны, верно? и можно было бы ожидать, что значения будут уменьшаться по мере удаления от диагонали, не так ли? Если вы сделаете вне диагонали достаточно маленькими, матрица гарантированно будет p.d. (так называемое диагональное доминирование). В любом случае, если вы заставите его работать с небольшой построенной матрицей, вы можете вернуться к оригиналу.   -  person Robert Dodier    schedule 12.02.2016

Ответы (1)


arrow_upward
4
arrow_downward

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

  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
person Zheyuan Li    schedule 03.09.2016
comment
Очевидно ... для кого? - person Brandon Bertelsen; 03.09.2016
comment
С моей стороны это было больше похоже на хихиканье. :) - person Brandon Bertelsen; 03.09.2016