Comment générer efficacement des matrices de corrélation positives-semi-définies aléatoires?

J'aimerais pouvoir générer efficacement des matrices de corrélation positive semi-définie (PSD). Ma méthode ralentit considérablement lorsque j'augmente la taille des matrices à générer.

1. Pourriez-vous suggérer des solutions efficaces? Si vous connaissez des exemples dans Matlab, je vous en serais très reconnaissant.
2. Lors de la génération d'une matrice de corrélation PSD, comment choisir les paramètres pour décrire les matrices à générer? Une corrélation moyenne, écart type de corrélations, valeurs propres?

Vous pouvez le faire en arrière: chaque matrice (l’ensemble de toutes les matrices PSD p × p symétriques ) peut être décomposée en tant que$C \in \mathbb{R}_{++}^p$$p \times p$

où O est une matrice orthonormée$C=O^{T}DO$$O$

Pour obtenir , d' abord générer une base aléatoire ( v 1 , . . . , V p ) (où v i sont des vecteurs aléatoires, typiquement dans ( - 1 , 1 ) ). De là, en utilisant le procédé d'orthogonalisation de Gram-Schmidt pour obtenir ( u 1 , . . . . , U p ) = $O$$(v_1,...,v_p)$$v_i$$(-1,1)$$(u_1,....,u_p)=O$

possède un certain nombre de packages qui peuvent effectuer l'orthogonalisation GS d'une base aléatoire de manière efficace, même pour les grandes dimensions, par exemple le package "far". Bien que vous trouviez l'algorithme GS sur wiki, il vaut probablement mieux ne pas réinventer la roue et opter pour une implémentation de matlab (il en existe sûrement, je ne peux en recommander aucune).$R$

Enfin, est une matrice diagonale dont tous les éléments sont positifs (encore une fois, c’est facile à générer: générez p nombres aléatoires, positionnez-les, triez-les et placez-les à la diagonale d’une matrice identité p sur p ).$D$$p$$p$$p$

Un article Génération de matrices de corrélation aléatoire à base de vignes et de méthode d'oignon étendu par Lewandowski, Kurowicka et Joe (LKJ), 2009, fournit un traitement et une présentation unifiés des deux méthodes efficaces de génération de matrices de corrélation aléatoire. Les deux méthodes permettent de générer des matrices à partir d'une distribution uniforme dans un certain sens défini ci-dessous, sont simples à mettre en œuvre, rapides et présentent l'avantage supplémentaire d'avoir des noms amusants.

Une matrice symétrique réelle de taille avec des unités sur la diagonale a d ( d - 1 ) / 2 éléments uniques hors diagonale et peut donc être paramétrée comme un point dans R d ( d - 1 ) / 2 . Chaque point de cet espace correspond à une matrice symétrique, mais tous ne sont pas définis positivement (comme doivent l'être les matrices de corrélation). Les matrices de corrélation forment donc un sous-ensemble de R d ( d - 1 ) /$d \times d$$d(d-1)/2$$\mathbb R^{d(d-1)/2}$$\mathbb R^{d(d-1)/2}$ (en réalité un sous-ensemble convexe connecté) et les deux méthodes peuvent générer des points à partir d’une distribution uniforme sur ce sous-ensemble.

Je vais fournir ma propre implémentation MATLAB de chaque méthode et les illustrer avec .$d=100$

Méthode de l'oignon

La méthode de l'oignon provient d'un autre document (référence n ° 3 dans LKJ) et tire son nom du fait que les matrices de corrélation sont générées en commençant par la matrice et en le développant colonne par colonne et ligne par ligne. La distribution résultante est uniforme. Je ne comprends pas vraiment le calcul derrière la méthode (et préfère la deuxième méthode de toute façon), mais voici le résultat:$1\times 1$

Ci-dessous et en dessous, le titre de chaque sous-parcelle indique les valeurs propres les plus petites et les plus grandes, ainsi que le déterminant (produit de toutes les valeurs propres). Voici le code:

%// ONION METHOD to generate random correlation matrices distributed randomly
function S = onion(d)
    S = 1;
    for k = 2:d
        y = betarnd((k-1)/2, (d-k)/2); %// sampling from beta distribution
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';             %// R is a square root of S
        q = R*w;
        S = [S q; q' 1];               %// increasing the matrix size
    end
end

Méthode de l'oignon étendu

Les LKJ modifient légèrement cette méthode afin de pouvoir échantillonner les matrices de corrélation partir d’une distribution proportionnelle à [ d e $\mathbf C$ . Plus le η est grand , plus le déterminant sera grand, ce qui signifie que les matrices de corrélation générées s'approcheront de plus en plus de la matrice d'identité. La valeur η = 1 correspond à une distribution uniforme. Sur la figure ci-dessous, les matrices sont générées avec η = 1 , 10 , 100 , 1000 , $[\mathrm{det}\:\mathbf C]^{\eta-1}$$\eta$$\eta=1$ .$\eta={1, 10, 100, 1000, 10\:000, 100\:000}$

Pour une raison quelconque d'obtenir un déterminant du même ordre de grandeur que dans la méthode de l'oignon à la vanille, je dois mettre et non pas η = 1 (comme l'affirme LKJ). Je ne sais pas où est l'erreur.$\eta=0$$\eta=1$

%// EXTENDED ONION METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = extendedOnion(d, eta)
    beta = eta + (d-2)/2;
    u = betarnd(beta, beta);
    r12 = 2*u - 1;
    S = [1 r12; r12 1];  

    for k = 3:d
        beta = beta - 1/2;
        y = betarnd((k-1)/2, beta);
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';
        q = R*w;
        S = [S q; q' 1];
    end
end

Méthode de la vigne

La méthode de la vigne a été proposée à l'origine par Joe (J en LKJ) et améliorée par LKJ. Je l’aime plus, parce que c’est conceptuellement plus facile et aussi plus facile à modifier. L'idée est de générer corrélations partielles (elles sont indépendantes et peuvent avoir n'importe quelle valeur de [ - 1 , 1 $d(d-1)/2$$[-1, 1]$sans contrainte) puis convertissez-les en corrélations brutes via une formule récursive. Il est pratique d'organiser le calcul dans un certain ordre. Ce graphe est appelé "vigne". De manière importante, si des corrélations partielles sont échantillonnées à partir de distributions bêta particulières (différentes pour différentes cellules de la matrice), la matrice résultante sera distribuée de manière uniforme. Là encore, LKJ introduit un paramètre supplémentaire à échantillonner à partir d’une distribution proportionnelle à [ d e $\eta$ . Le résultat est identique à l'oignon étendu:$[\mathrm{det}\:\mathbf C]^{\eta-1}$

%// VINE METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = vine(d, eta)
    beta = eta + (d-1)/2;   
    P = zeros(d);           %// storing partial correlations
    S = eye(d);

    for k = 1:d-1
        beta = beta - 1/2;
        for i = k+1:d
            P(k,i) = betarnd(beta,beta); %// sampling from beta
            P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
            p = P(k,i);
            for l = (k-1):-1:1 %// converting partial correlation to raw correlation
                p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
            end
            S(k,i) = p;
            S(i,k) = p;
        end
    end
end

Méthode de la vigne avec échantillonnage manuel des corrélations partielles

$\pm 1$$[0,1]$$[-1, 1]$$\alpha=\beta={50, 20, 10, 5, 2, 1}$. Plus les paramètres de la distribution bêta sont petits, plus elle est concentrée près des bords.

Notez que dans ce cas, la distribution n'est pas garantie d'être invariante à la permutation, aussi permutez-vous de manière aléatoire de manière aléatoire les lignes et les colonnes après la génération.

%// VINE METHOD to generate random correlation matrices
%// with all partial correlations distributed ~ beta(betaparam,betaparam)
%// rescaled to [-1, 1]
function S = vineBeta(d, betaparam)
    P = zeros(d);           %// storing partial correlations
    S = eye(d);

    for k = 1:d-1
        for i = k+1:d
            P(k,i) = betarnd(betaparam,betaparam); %// sampling from beta
            P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
            p = P(k,i);
            for l = (k-1):-1:1 %// converting partial correlation to raw correlation
                p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
            end
            S(k,i) = p;
            S(i,k) = p;
        end
    end

    %// permuting the variables to make the distribution permutation-invariant
    permutation = randperm(d);
    S = S(permutation, permutation);
end

Voici comment les histogrammes des éléments hors diagonale recherchent les matrices ci-dessus (la variance de la distribution augmente de façon monotone):

Mise à jour: utilisation de facteurs aléatoires

$k<d$$\mathbf W$$k \times d$$\mathbf W \mathbf W^\top$$\mathbf D$$\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$$\mathbf C = \mathbf E^{-1/2}\mathbf B \mathbf E^{-1/2}$$\mathbf E$$\mathbf B$$k={100, 50, 20, 10, 5, 1}$

Et le code:

%// FACTOR method
function S = factor(d,k)
    W = randn(d,k);
    S = W*W' + diag(rand(1,d));
    S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));
end

Voici le code d'emballage utilisé pour générer les chiffres:

d = 100; %// size of the correlation matrix

figure('Position', [100 100 1100 600])
for repetition = 1:6
    S = onion(d);

    %// etas = [1 10 100 1000 1e+4 1e+5];
    %// S = extendedOnion(d, etas(repetition));

    %// S = vine(d, etas(repetition));

    %// betaparams = [50 20 10 5 2 1];
    %// S = vineBeta(d, betaparams(repetition));

    subplot(2,3,repetition)

    %// use this to plot colormaps of S
    imagesc(S, [-1 1])
    axis square
    title(['Eigs: ' num2str(min(eig(S)),2) '...' num2str(max(eig(S)),2) ', det=' num2str(det(S),2)])

    %// use this to plot histograms of the off-diagonal elements
    %// offd = S(logical(ones(size(S))-eye(size(S))));
    %// hist(offd)
    %// xlim([-1 1])
end

Une caractérisation encore plus simple est celle de la matrice réelle $A$, $A^TA$est positif semi-défini. Pour voir pourquoi c'est le cas, il suffit de prouver que$y^T (A^TA) y \ge 0$ pour tous les vecteurs $y$(de la bonne taille, bien sûr). C'est trivial:$y^T (A^TA) y = (Ay)^T Ay = ||Ay||$ce qui est non négatif. Donc, dans Matlab, essayez simplement

A = randn(m,n);   %here n is the desired size of the final matrix, and m > n
X = A' * A;

Selon l’application, il se peut que cela ne vous donne pas la distribution des valeurs propres souhaitée; La réponse de Kwak est bien meilleure à cet égard. Les valeurs propres de Xcet extrait de code doivent correspondre à la distribution de Marchenko-Pastur.

Pour simuler les matrices de corrélation des actions, supposons une approche légèrement différente:

k = 7;      % # of latent dimensions;
n = 100;    % # of stocks;
A = 0.01 * randn(k,n);  % 'hedgeable risk'
D = diag(0.001 * randn(n,1));   % 'idiosyncratic risk'
X = A'*A + D;
ascii_hist(eig(X));    % this is my own function, you do a hist(eig(X));
-Inf <= x <  -0.001 : **************** (17)
-0.001 <= x <   0.001 : ************************************************** (53)
 0.001 <= x <   0.002 : ******************** (21)
 0.002 <= x <   0.004 : ** (2)
 0.004 <= x <   0.005 :  (0)
 0.005 <= x <   0.007 : * (1)
 0.007 <= x <   0.008 : * (1)
 0.008 <= x <   0.009 : *** (3)
 0.009 <= x <   0.011 : * (1)
 0.011 <= x <     Inf : * (1)

En variante à la réponse du kwak: générer une matrice diagonale $\mathbf{D}$ avec des valeurs propres non négatives aléatoires à partir d'une distribution de votre choix, puis effectuez une transformation de similarité $\mathbf{A}=\mathbf{Q}\mathbf{D}\mathbf{Q}^T$ avec $\mathbf{Q}$une matrice orthogonale pseudo-aléatoire distribuée par Haar .

Vous n'avez pas spécifié de distribution pour les matrices. Deux distributions communes sont les distributions de Wishart et inverses de Wishart. La décomposition de Bartlett donne une factorisation de Cholesky d'une matrice de Wishart aléatoire (qui peut également être résolue efficacement pour obtenir une matrice de Wishart inverse aléatoire).

En fait, l'espace de Cholesky est un moyen pratique de générer d'autres types de matrices de PSD aléatoires, car il vous suffit de vous assurer que la diagonale est non négative.

La méthode la plus simple est celle ci-dessus, qui consiste à simuler un ensemble de données aléatoires et à calculer le gramian . Un mot de prudence: la matrice résultante ne sera pas uniformément aléatoire, en ce sens que sa décomposition, disons$U^TSU$Les rotations ne seront pas réparties conformément à la mesure de Haar. Si vous souhaitez disposer de matrices PSD "uniformément réparties", vous pouvez utiliser l'une des approches décrites ici .

Si vous souhaitez mieux contrôler votre matrice PSD symétrique générée, par exemple générer un jeu de données de validation synthétique, vous disposez d'un certain nombre de paramètres. Une matrice PSD symétrique correspond à une hyper-ellipse dans l'espace à N dimensions, avec tous les degrés de liberté correspondants:

1. Rotations
2. Longueurs des axes.

Ainsi, pour une matrice à 2 dimensions (c.-à-d. Ellipse 2d), vous aurez 1 rotation + 2 axes = 3 paramètres.

Si les rotations rappellent les matrices orthogonales, c’est un train correct, car la construction est à nouveau $\Sigma=ODO^T$, avec $\Sigma$ étant la matrice Sym.PSD produite, $O$ la matrice de rotation (qui est orthogonale), et $D$la matrice diagonale, dont les éléments diagonaux contrôleront la longueur des axes de l'ellipse.

Le code Matlab ci-après représente 16 ensembles de données distribués gaussiens bidimensionnels basés sur $\Sigma$, avec un angle croissant. Le code pour la génération aléatoire de paramètres est dans les commentaires.

figure;
mu = [0,0];
for i=1:16
    subplot(4,4,i)
    theta = (i/16)*2*pi;   % theta = rand*2*pi;
    U=[cos(theta), -sin(theta); sin(theta) cos(theta)];
    % The diagonal's elements control the lengths of the axes
    D = [10, 0; 0, 1]; % D = diag(rand(2,1));    
    sigma = U*D*U';
    data = mvnrnd(mu,sigma,1000);
    plot(data(:,1),data(:,2),'+'); axis([-6 6 -6 6]); hold on;
end

Pour plus de dimensions, la matrice Diagonal est simple (comme ci-dessus), et le $U$ découler de la multiplication des matrices de rotation.

Une méthode peu coûteuse et amusante que j’ai utilisée pour tester consiste à générer m N (0,1) n-vecteurs V [k], puis à utiliser P = d * I + Somme {V [k] * V [k] '} en tant que matrice nxn psd. Avec m <n, cela sera singulier pour d = 0, et pour petit d aura un nombre de condition élevé.