-
Notifications
You must be signed in to change notification settings - Fork 33
/
BiCopTau2Par.R
230 lines (216 loc) · 8.12 KB
/
BiCopTau2Par.R
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
#' Parameter of a Bivariate Copula for a given Kendall's Tau Value
#'
#' This function computes the parameter of a (one parameter) bivariate copula
#' for a given value of Kendall's tau.
#'
#'
#' @param family integer; single number or vector of size `n`; defines the
#' bivariate copula family: \cr `0` = independence copula \cr `1` =
#' Gaussian copula \cr `2` = Student t copula (Here only the first
#' parameter can be computed) \cr `3` = Clayton copula \cr `4` =
#' Gumbel copula \cr `5` = Frank copula \cr `6` = Joe copula \cr
#' `13` = rotated Clayton copula (180 degrees; ``survival Clayton'') \cr
#' `14` = rotated Gumbel copula (180 degrees; ``survival Gumbel'') \cr
#' `16` = rotated Joe copula (180 degrees; ``survival Joe'') \cr `23`
#' = rotated Clayton copula (90 degrees) \cr `24` = rotated Gumbel copula
#' (90 degrees) \cr `26` = rotated Joe copula (90 degrees) \cr `33` =
#' rotated Clayton copula (270 degrees) \cr `34` = rotated Gumbel copula
#' (270 degrees) \cr `36` = rotated Joe copula (270 degrees)\cr Note that
#' (with exception of the t-copula) two parameter bivariate copula families
#' cannot be used.
#' @param tau numeric; single number or vector of size `n`; Kendall's tau
#' value (vector with elements in \eqn{[-1,1]}).
#' @param check.taus logical; default is `TRUE`; if `FALSE`, checks
#' for family/tau-consistency are omitted (should only be used with care).
#'
#' @return Parameter (vector) corresponding to the bivariate copula family and
#' the value(s) of Kendall's tau (\eqn{\tau}). \tabular{ll}{ No.
#' (`family`) \tab Parameter (`par`) \cr `1, 2` \tab
#' \eqn{\sin(\tau \frac{\pi}{2})}{sin(\tau \pi/2)} \cr `3, 13` \tab
#' \eqn{2\frac{\tau}{1-\tau}}{2\tau/(1-\tau)} \cr `4, 14` \tab
#' \eqn{\frac{1}{1-\tau}}{1/(1-\tau)} \cr `5` \tab no closed form
#' expression (numerical inversion) \cr `6, 16` \tab no closed form
#' expression (numerical inversion) \cr `23, 33` \tab
#' \eqn{2\frac{\tau}{1+\tau}}{2\tau/(1+\tau)} \cr `24, 34` \tab
#' \eqn{-\frac{1}{1+\tau}}{-1/(1+\tau)} \cr `26, 36` \tab no closed form
#' expression (numerical inversion) }
#'
#' @note The number `n` can be chosen arbitrarily, but must agree across
#' arguments.
#'
#' @author Jakob Stoeber, Eike Brechmann, Tobias Erhardt
#'
#' @seealso [BiCopPar2Tau()]
#'
#' @references Joe, H. (1997). Multivariate Models and Dependence Concepts.
#' Chapman and Hall, London.
#'
#' Czado, C., U. Schepsmeier, and A. Min (2012). Maximum likelihood estimation
#' of mixed C-vines with application to exchange rates. Statistical Modelling,
#' 12(3), 229-255.
#'
#' @examples
#' ## Example 1: Gaussian copula
#' tau0 <- 0.5
#' rho <- BiCopTau2Par(family = 1, tau = tau0)
#' BiCop(1, tau = tau0)$par # alternative
#'
#' ## Example 2:
#' vtau <- seq(from = 0.1, to = 0.8, length.out = 100)
#' thetaC <- BiCopTau2Par(family = 3, tau = vtau)
#' thetaG <- BiCopTau2Par(family = 4, tau = vtau)
#' thetaF <- BiCopTau2Par(family = 5, tau = vtau)
#' thetaJ <- BiCopTau2Par(family = 6, tau = vtau)
#' plot(thetaC ~ vtau, type = "l", ylim = range(thetaF))
#' lines(thetaG ~ vtau, col = 2)
#' lines(thetaF ~ vtau, col = 3)
#' lines(thetaJ ~ vtau, col = 4)
#'
#' ## Example 3: different copula families
#' theta <- BiCopTau2Par(family = c(3,4,6), tau = c(0.4, 0.5, 0.6))
#' BiCopPar2Tau(family = c(3,4,6), par = theta)
#'
#' \dontshow{
#' # Test BiCopTau2Par
#' BiCopTau2Par(family = 0, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 1, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 2, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 3, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 4, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 5, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 6, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 13, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 14, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 16, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 23, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 24, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 26, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 33, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 34, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 36, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 41, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 51, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 61, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 71, tau = c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 41, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 51, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 61, tau = -c(0.4,0.5,0.6))
#' BiCopTau2Par(family = 71, tau = -c(0.4,0.5,0.6))
#' }
#'
BiCopTau2Par <- function(family, tau, check.taus = TRUE) {
## sanity check
if (any(family %in% setdiff(allfams[twopar], 2)))
stop("For two parameter copulas (except t) Kendall's tau cannot be inverted.")
if (check.taus && any(abs(tau) > 0.99999))
stop("some tau is too close to -1 or 1")
# fix for SemiParBIVProbit package
dims <- set_dims(family, tau = tau)
## adjust length for input vectors; stop if not matching
family <- c(family)
tau <- c(tau)
n <- max(length(family), length(tau))
if (length(family) == 1)
family <- rep(family, n)
if (length(tau) == 1)
par <- rep(tau, n)
if (!all(c(length(family), length(tau)) %in% c(1, n)))
stop("Input lenghts don't match")
## check for family/tau consistency
if (check.taus)
BiCopCheckTaus(family, tau)
## calculate the parameter
par <- vapply(seq_along(tau),
function(i) calcPar(family[i], tau[i]),
numeric(1))
par <- vapply(seq_along(par),
function(i) adjustPars(family[i], par[i], 0)[1],
numeric(1))
## return result
if (length(dims) > 1)
par <- array(par, dim = dims)
par
}
calcPar <- function(family, tau) {
## calculation of parameter(s) depending on pair-copula family
if (family == 0) {
par <- rep(0, times = length(tau))
} else if (family %in% 1:2) {
par <- sin(pi * tau/2)
} else if (family %in% c(3, 13)) {
par <- 2 * tau/(1 - tau)
} else if (family %in% c(4, 14)) {
par <- 1/(1 - tau)
} else if (family == 5) {
par <- if (tau == 0) 0 else Frank.itau.JJ(tau)
} else if (family %in% c(6, 16)) {
par <- Joe.itau.JJ(tau)
} else if (family %in% c(23, 33)) {
par <- 2 * tau/(1 + tau)
} else if (family %in% c(24, 34)) {
par <- -(1/(1 + tau))
} else if (family %in% c(26, 36)) {
par <- -Joe.itau.JJ(-tau)
} else if (family %in% c(41, 51)) {
par <- ipsA.tau2cpar(tau)
} else if (family %in% c(61, 71)) {
par <- -ipsA.tau2cpar(-tau)
}
## return result
par
}
Frank.itau.JJ <- function(tau) {
if (abs(tau) > 0.99999) return(Inf)
a <- 1
if (tau < 0) {
a <- -1
tau <- -tau
}
v <- uniroot(function(x) tau - frankTau(x),
lower = 0 + .Machine$double.eps^0.5, upper = 100,
tol = .Machine$double.eps^0.5)$root
return(a*v)
}
Joe.itau.JJ <- function(tau) {
if (abs(tau) > 0.99999) return(Inf)
if (tau < 0) {
return(1.000001)
} else {
tauF <- function(par) {
param1 <- 2/par + 1
tem <- digamma(2) - digamma(param1)
tau <- 1 + tem * 2/(2 - par)
tau[par == 2] <- 1 - trigamma(2)
tau
}
v <- uniroot(function(x) tau - tauF(x),
lower = 1,
upper = 5e5,
tol = .Machine$double.eps^0.5)$root
return(v)
}
}
ipsA.tau2cpar <- function(tau, mxiter = 20, eps = 1e-06, dstart = 0, iprint = FALSE) {
con <- log((1 - tau) * sqrt(pi)/2)
de <- dstart
if (dstart <= 0)
de <- tau + 1
iter <- 0
diff <- 1
while (iter < mxiter & max(abs(diff)) > eps) {
g <- con + lgamma(1 + de) - lgamma(de + 0.5)
gp <- digamma(1 + de) - digamma(de + 0.5)
iter <- iter + 1
diff <- g/gp
de <- de - diff
while (min(de) <= 0) {
diff <- diff/2
de <- de + diff
}
if (iprint)
cat(iter, " ", de, " ", diff, "\n")
}
if (iter >= mxiter)
cat("did not converge\n")
de
}