diff --git a/exercise223/LognormalCopulaPDF.m b/exercise223/LognormalCopulaPDF.m new file mode 100644 index 0000000..f0b6d4e --- /dev/null +++ b/exercise223/LognormalCopulaPDF.m @@ -0,0 +1,23 @@ +function F_U = LognormalCopulaPDF(u,Mu,Sigma) +% this function computes the pdf of the copula of the lognormal distribution +% at the generic point u in the unit hypercube +% see formula 2.30 +% we need the joint distribution of X as numerator +% and product of marginal distributions as denominator + +N=length(u); % dimension of the hypercube +s=sqrt(diag(Sigma)); % st. deviations + +x=logninv(u,Mu,s); %from uniform we generate a lognormal sample with the + % inverse of lognormal cdf + +%numerator is the joint distribution of a lognormal multivariate variable as defined in 2.213 and 2.156 +Numerator = (2*pi)^(-N/2) * ( (det(Sigma))^(-.5) ) /prod(x) * exp(-.5*(log(x)-Mu)'*inv(Sigma)*(log(x)-Mu)); + +%denominator is the product of marginal distributions +fs=lognpdf(x,Mu,s); %returns values of the lognormal pdf at x +Denominator = prod(fs); + +F_U = Numerator/Denominator; %copula pdf as in 2.30 + + diff --git a/exercise223/exercise223_main.m b/exercise223/exercise223_main.m new file mode 100644 index 0000000..bcb765e --- /dev/null +++ b/exercise223/exercise223_main.m @@ -0,0 +1,42 @@ +% Exercise 2.2.3 - main +% *Author: Giuseppe Mascolo* +% Logonrmal copula pdf + +close all; +clc; clear; + +% input parameters +Mu=[1 -1]'; % exp values +r=-0.7; % correlation +sigmas=[1 1]'; % st. deviations + +%if we want a specific Sigma we just have to change sigmas and r +Sigma=diag(sigmas)*[1 r;r 1]*diag(sigmas); + +%% compute and display results + +%creates grid +GridSide1=[.05:.05:.95]; +GridSide2=GridSide1; +NumGrid=length(GridSide1); + +f_U=zeros(NumGrid); %in order to plot + +%loop inside the grid +for j=1:NumGrid + for k=1:NumGrid + u=[GridSide1(j) + GridSide2(k)]; + f_U(j,k)=LognormalCopulaPDF(u,Mu,Sigma); %function returns copula pdf + end +end + + +[G1,G2]=meshgrid(GridSide1,GridSide2); %creates 2D grid +surf(G1,G2,f_U) %3D surface plot +xlabel('U_1') +ylabel('U_2') +zlabel('Copula pdf') + +%notice that when r is zero we have a flat copula pdf, which means +%independence