% Comparaison des deux méthodes de décomposition d'une matrice de
% variances-covariances: Choleski et PCA (Principal Components Analysis)
% Auteur: A. Hocquart

clear all; close all; clc

% Paramètres de la simulation
nb_simul = [1000 10000 50000 100000];           %Nombres de trajectoires
k = 2:1:25;                                     %Dimensions de la matrice de var-cov
temps = zeros(length(k),2,length(nb_simul));    %Matrice qui entrepose les temps de calcul

for i = 1:1:length(nb_simul)                    %Boucles sur les différentes simulations    
    nb = nb_simul(i);                           %nb de trajectoires   
    for j = 1:1:length(k)
        k_temp = k(j);                          %dimension de la matrice de var-cov
        
        sigma = ones(k_temp,1);                 %Definition du vecteur d'écarts types
        %Construction de la matrice de corrélations
        rho = ones(k_temp,k_temp);              %Préallocation de la mémoire
        for m = 1:k_temp
            for n = m:k_temp
                if m ~= n
                    rho(m,n) = m/n;             %Alors rho(i,j) = i/j;
                    rho(n,m) = m/n;             %Matrice symétrique
                end
            end
        end

        %Construction de la matrice de var-cov
        Varcov = diag(sigma)*rho*diag(sigma);
        % Le temps de calcul étant assez faible, on captera l'influence apportées par
        % les valeurs générées pas le générateur de nombres pseudo-aléatoire, et
        % par la charge initiale du processeur au moment de l'exécution
        % du code
        
        Z = randn(nb,k_temp); % On génere une matrice de nxk de v.a. iid N(0,1)

        %Décomposition de Cholesky
        %-------------------------
        tic;                % Initialisation du chronomètre
        C = chol(Varcov);   % Décomposition de Cholesky de la matrice var-cov
        X = Z*C;            % On obtient une matrice de vecteurs de N(0,Varcov)
        temps(j,1,i) = toc; % Fin du chronomètre
        %On peut vérifier la matrice de corrélation
        %rho_tilde = corr(X)

        % Decomposition en composantes principales
        % ----------------------------------------
        tic;                    % Initialisation du chronomètre
        [V,D] = eig(Varcov);    % Valeurs et vecteurs propres
        A = V*sqrt(D);          
        X = (A*Z')';            % On obtient une matrice de vecteurs de N(0,Varcov)
        temps(j,2,i) = toc;     % Fin du chronomètre
        % On peut vérifier la matrice de corrélation
        % rho_tilde = corr(X)

        clear sigma rho Varcov C Z X rho_tilde
    end
end

%%% ON DECIDE DE NE PAS CONSIDERER LA PREMIERE BOUCLE SUR LES TRAJECTOIRES
%%% (1000 simulations) CAR MATLAB PREND UN TEMPS ADDITIONNEL POUR
%%% INITIALISER SES VARIABLES. LE TEMPS DE CALCUL POUR LA PREMIERE
%%% ITERATION DE LA BOUCLE SUR LE NOMBRE DE TRAJECTOIRES EST DONC BIAISE A
%%% LA HAUSSE

A = permute(temps,[3 1 2]);
figure;
    plot(k,A(2:end,:,1));
    hold on;
    plot(k,A(2:end,:,2),'--');
    hold off;
    xlabel('Nombre d''actifs');
    ylabel('Temps de calcul');
    legend(num2str((nb_simul(2:length(nb_simul)))'),'Location','NorthWest');
    title('Comparaison des Methodes, PCA en pointille, Choleski en trait continu');