How to integrate over a bivariate gaussian copula using copulapdf?

19 views (last 30 days)
Thibaut
Thibaut on 30 May 2014
Answered: Mike Hosea on 13 Dec 2014
Hi,
I wanted to estimate a conditional tail expectation using a gaussian copula structure of the statistics toolbox in R2012b.
Essentially what I want to integrate is the following formula of a conditional expectation: %E[x|y<C] = integrate_inf^inf { integrate_inf^C { x * f(x,y) dy } dx %where f(x,y) = copula(F(x),F(y)) * f(x) * f(y) I use several handle for the copula and the marginal density distributions, where the copula is gaussian and the marginal densities are from kernel estimation. The codes below reproduce my problem. I get the error message:
"RHO must be a square correlation matrix with size equal to the number of columns in U."
So essentially, it does not succeed in using the handle with the copula (when I remove the copula but keep the rest of the formula, it is fine). What is wrong? did I pass the argument incorrectly to the copula when I do the double integral?
%--------------------------------------------------------------------------
%BIVARIATE GAUSSIAN COPULA INTEGRATION TO ESTIMATE E[x|y<C]
%--------------------------------------------------------------------------
load stockreturns
s1=stocks(:,1);
s2=stocks(:,2);
scatterhist(x,y)
%E[x|y<C] = integrate_inf^inf { integrate_inf^C { x * f(x,y) dy } dx
%where f(x,y) = copula(F(x),F(y)) * f(x) * f(y)
%2) Define kernel density and CDF estimates of x and y
%kernel density
fhatx = @(x) ksdensity(s1, x, 'function','pdf'); %f(x)
fhaty = @(y) ksdensity(s2, y, 'function','pdf'); %f(y)
%kernel CDF:
Fhatx = @(x) ksdensity(s1, x, 'function','cdf'); %F(x)
Fhaty = @(y) ksdensity(s2, y, 'function','cdf'); %F(y)
%3) Get copula estimates
%define copula space
U=zeros(length(s1),2);
U(:,1)=Fhatx(s1);
U(:,2)=Fhaty(s2);
%estimate the gaussian copula and obtain the correlation matrix Rhohat
Rhohat = copulafit('Gaussian',U); % copula(F(x),F(y))
%4) Get conditional tail expectations E[x|y<C]
%set threshold value for conditioning variable
C=-0.02;
%define the integrand : x * copula(F(x),F(y)) * f(x) * f(y)
copfun = @(x,y) copulapdf('Gaussian',[Fhatx(x),Fhaty(y)],Rhohat);
densityfun = @(x,y) fhatx(x).*fhaty(y);
JointDensity = @(x,y) x.*copfun(x,y).*densityfun(x,y);
%evaluate the expectation : double integration over the whole space of x and the truncated space of y
CondTailExpecation = integral2(JointDensity,-inf,+inf,-inf,C); %(fun,xmin,xmax,ymin,ymax)
If I recode the gaussian copula myself, I manage to make it work, although it fails to converge in this example.
%ELSE IF I do it by rewriting the gaussian copula function in the bivariate case copfun2:
Rtemp=(Rhohat^(-1)-eye(2));
copfun2 = @(x,y) (1/(det(Rhohat))^(1/2))*exp((-1/2)* (norminv(Fhatx(x)).*(norminv(Fhatx(x))*Rtemp(1,1)+norminv(Fhaty(y))*Rtemp(2,1)) + norminv(Fhaty(y)).*(norminv(Fhatx(x))*Rtemp(1,2)+norminv(Fhaty(y))*Rtemp(2,2))));
JointDensity2 = @(x,y) x.*copfun2(x,y).*densityfun(x,y);
q = integral2(JointDensity2,-inf,inf,-inf,C); %(fun,xmin,xmax,ymin,ymax)
Thanks in advance for your help.

Answers (2)

Yakup
Yakup on 13 Dec 2014
The first error must be about your copula coefficient Rho, which can be a scalar variable when the size of the copula function is equal to two but when the size of U is greater than 2, you need insert a covariance matrix with all variances equal to 1 and diagonal consisting of correlation variables. Perhaps
copulafit
is a function that returns a vector rather than a matrix which you must be putting into
copulapdf/cdf
functions. I may be wrong but hope this helps.

Mike Hosea
Mike Hosea on 13 Dec 2014
Before you call integral2, make sure your JointDensity2 function satisfies the vectorization requirement.
x = rand(2,3); % or some other values that make sense
y = rand(2,3); % or some other values that make sense
z1 = JointDensity2(x,y);
z2 = arrayfun(JointDensity2,x,y);
isequal(z1,z2)
If these values aren't equal, or of the function errors out, work on that problem first. If I were you, I'd probably just go straight to the arrayfun version
JointDensity2v = @(x,y)arrayfun(JointDensity2,x,y);
and pass that to integral2. That way, just need to make sure JointDensity2 works with scalar inputs. You only need make it work with array inputs if you want to speed things up.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!