-
Notifications
You must be signed in to change notification settings - Fork 6
/
copulapdfadv.m
77 lines (60 loc) · 2.49 KB
/
copulapdfadv.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
function pdf = copulapdfadv(family,u,theta)
% Computes the pdf of a copula at u.
%
% call: pdf = copulapdfadv(family,u,theta)
%
% input family - the copula family: 'gumbel', 'clayton', 'frank', 't',
% 'gauss', 'ind', 'amhaq', 'tawn',
% 'fgm', 'plackett', 'joe',
% 'surclayton', 'surgumbel',
% 'surjoe'
% u - nx2 matrix of points to be evaluated
% theta - copula parameter: for t-copula [rho, nu]
%
% output pdf - the pdf at u
%
%
% Copyright 2020, Maximilian Coblenz
% This code is released under the 3-clause BSD license.
%
% some parsing
p = inputParser;
p.addRequired('family',@isstr);
p.addRequired('u',@ismatrix);
p.addRequired('theta',@isvector);
p.parse(family,u,theta);
% sanity checks
if ~cpcheck(family,theta)
error('copulapdfadv:InvalidParameter',['invalid parameter for ',family,' copula']);
end
% compute pdf
switch lower(family)
case{'gauss','gumbel','clayton','frank'}
pdf = copulapdf(family,u,theta);
case('t')
pdf = copulapdf(family,u,theta(1),theta(2));
case('ind')
pdf = ones(size(u,1),1);
case('amhaq')
pdf = (1-theta+2*theta*u(:,1).*u(:,2)-theta*(1-theta)*(1-u(:,1)).*(1-u(:,2)))./...
(1-theta*(1-u(:,1)).*(1-u(:,2))).^3;
case('tawn')
pdf = exp(-theta*log(u(:,1)).*log(u(:,2))./log(u(:,1).*u(:,2))).*...
((1-theta*(log(u(:,2))).^2./(log(u(:,1).*u(:,2))).^2).*(1-theta*(log(u(:,1))).^2./(log(u(:,1).*u(:,2))).^2)-...
theta*(2*log(u(:,1)).*log(u(:,2)))./(log(u(:,1).*u(:,2))).^3);
case('joe')
pdf = ((1-u(:,1)).^theta+(1-u(:,2)).^theta-(1-u(:,1)).^theta.*(1-u(:,2)).^theta).^(1/theta-2).*...
(1-u(:,1)).^(theta-1).*(1-u(:,2)).^(theta-1).*(theta-1+(1-u(:,1)).^theta+(1-u(:,2)).^theta-...
(1-u(:,1)).^theta.*(1-u(:,2)).^theta);
case('fgm')
pdf = 1+theta*(1-2*u(:,1)).*(1-2*u(:,2));
case('plackett')
pdf = theta*(1+(theta-1)*(u(:,1)+u(:,2)-2*u(:,1).*u(:,2)))./...
((1+(theta-1)*(u(:,1)+u(:,2))).^2-4*theta*(theta-1)*u(:,1).*u(:,2)).^(3/2);
case{'surclayton','surgumbel','surjoe'}
pdf = copulapdfadv(erase(lower(family),'sur'),1-u,theta);
% add new copulas here
% case('XXX')
% pdf = formula;
end % switch
end