simva.m

Gaussiennes, Gaussienne standard and Lognormal型随机变量的数值模拟

%==========================================================================
%=== SIMULATION NUMERIQUE DE VARIABLES ALEATOIRES VECTORIELLES          ===
%=== GAUSSIENNES ET LOGNORMALES                                         ===
%=== Auteurs : P. Grange, M. Fogli (Juin 2006)                          ===
%==========================================================================

function[X]=simva(va,mx,Cx,N)

% Donn閑s
%==========================================================================
% va : type de v.a. ?simuler (gaussienne standard, gaussienne ou lognormale). Il faut saisir
% le num閞o correspondant au type de v.a. ?simuler : 
% 1 - gaussienne standard 
% 2 - gaussienne
% 3 - lognormale

% n : dimension de la v.a. simul閑
% mx : moyenne de la v.a. ?simuler (de la forme n*1)
% Cx : matrice de covariance de la v.a. ?simuler (de la forme n*n)
% N : nombre de copies ind閜endantes de la v.a.

% NB : La fonction "boxmuller" est une fonction
% permettant de simuler une v.a. gaussienne standard n-dimensionnelle.

% Exemple d'utilisation du programme
%==========================================================================
% mx=[10;17.25;11]; moyenne th閛rique de la v.a.
% Cx=[0.25,0.03,0.22;0.03,2.6,0.47;0.22,0.47,1.14]; matrice de covariance
% th閛rique de la v.a.
% [X]=simva(3,mx,Cx,100000); simulation de 100000 copies ind閜endantes 
% d'une v.a. lognormale vectorielle (de dimension 3 -d'apr鑣 dimension de 
% mx et Cx-), de moyenne mx et de matrice de covariance Cx.

n=size(mx,1);

% Algorithme de simulation de v.a. gaussienne standard vectorielle 
% n-dimensionnelle
%==========================================================================
if va==1
    [X]=boxmuller(n,N);

% Algorithme de simulation de v.a. gaussienne vectorielle n-dimensionnelle
%==========================================================================
elseif va==2
    if n==1
        [GS]=boxmuller(n,N);
        X=sqrt(Cx)*GS;
        for j=1:N
            X(:,j)=mx+X(:,j);
        end
    elseif n==2
        [GS]=boxmuller(n,N);
        Lx(1,1)=sqrt(Cx(1,1));
        Lx(1,2)=0;
        Lx(2,1)=Cx(1,2)/sqrt(Cx(1,1));
        Lx(2,2)=sqrt(Cx(2,2)-Cx(1,2)^2/Cx(1,1));
        X=Lx*GS;
        for j=1:N
            X(:,j)=mx+X(:,j);
        end
    else
    [GS]=boxmuller(n,N);
    Lx=chol(Cx);
    Lx=Lx';
    X=Lx*GS;
        for j=1:N
            X(:,j)=mx+X(:,j);
        end
    end

% Algorithme de simulation de v.a. lognormale vectorielle n-dimensionnelle
%==========================================================================
else
    if n==1
        my=log(mx/sqrt(1+Cx/mx^2));
        sigmay=sqrt(log(1+Cx/mx^2));
        [GS]=boxmuller(n,N);
        X=exp(my)*exp(sigmay*GS);
    elseif n==2
        my1=log(mx(1,1)/sqrt(1+Cx(1,1)/mx(1,1)^2));
        my2=log(mx(2,1)/sqrt(1+Cx(2,2)/mx(2,1)^2));
        sigmay1=sqrt(log(1+Cx(1,1)/mx(1,1)^2));
        sigmay2=sqrt(log(1+Cx(2,2)/mx(2,1)^2));
        rhoy1y2=log(1+Cx(1,2)/(mx(1,1)*mx(2,1)))/(sigmay1*sigmay2);
        [GS]=boxmuller(n,N);
        Y(1,:)=sigmay1*GS(1,:);
        Y(2,:)=rhoy1y2*sigmay2*GS(1,:)+sqrt(1-rhoy1y2^2)*sigmay2*GS(2,:);
        for j=1:N
            Y(1,j)=my1+Y(1,j);
            Y(2,j)=my2+Y(2,j);
        end
        X=exp(Y);
    else
        for p=1:n
            cvx(p,1)=sqrt(Cx(p,p))/mx(p,1); % coeff de variation de X
            my(p,1)=log(mx(p,1)/sqrt(1+cvx(p,1)^2)); % moyenne de la v.a. Y gaussienne associ閑 ?X lognormale
            for q=1:n
                Cy(p,q)=log(1+Cx(p,q)/(mx(p,1)*mx(q,1))); % matrice de covariance de Y
            end
        end
        [GS]=boxmuller(n,N);
        Ly=chol(Cy);
        Ly=Ly';
        Y=Ly*GS;
        for j=1:N
            Y(:,j)=my+Y(:,j);
        end
        X=exp(Y);
    end
end

% Sauvegarde des r閍lisations de la v.a. dans un fichier Excel
%==========================================================================
%xlswrite('Simulation v.a.',X);
%save('test',)
X=X';
save('test.xls','X','-ASCII','-double') %si nombre de ligne<65 536
%save('test.mat','X','-ASCII','-double') %si un fichier de plus de 65536 lignes ou 256 colonnes.