How to integrate over a bivariate gaussian copula using copulapdf?

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.

Answer 1

[Yakup]

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.

Answer 2

[Mike Hosea]

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.