-
Notifications
You must be signed in to change notification settings - Fork 5
/
coverage.R
296 lines (278 loc) · 12.4 KB
/
coverage.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
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
### actuar: Actuarial Functions and Heavy Tailed Distributions
###
### Create modified density and modified cumulative distribution
### function for data with deductible, limit, coinsurance and
### inflation.
###
### See Chapter 8 of Klugman, Panjer & Willmot, Loss Models, Fourth
### Edition, Wiley, 2012.
###
### AUTHORS: Mathieu Pigeon, Vincent Goulet <[email protected]>
coverage <- function(pdf, cdf, deductible = 0, franchise = FALSE,
limit = Inf, coinsurance = 1, inflation = 0,
per.loss = FALSE)
{
Call <- match.call()
## First determine if the cdf is needed or not. It is needed when
## there is a deductible or a limit and, of course, if the output
## function should compute the cdf.
is.cdf <- missing(pdf) || is.null(pdf) # return cdf?
has.limit <- limit < Inf # used often
needs.cdf <- any(deductible > 0, has.limit, is.cdf) # cdf needed?
## Sanity check of arguments
if (any(deductible < 0, limit < 0, coinsurance < 0, inflation < 0))
stop("coverage modifications must be positive")
if (limit <= deductible)
stop("deductible must be smaller than the limit")
if (coinsurance > 1)
stop("coinsurance must be between 0 and 1")
if (missing(cdf) & needs.cdf)
stop(sprintf("%s must be supplied", sQuote("cdf")))
## Quantites often used. Leave as expressions for the output
## function to preserve accuracy.
r <- 1 + inflation
d <- if (inflation) substitute(d/r, list(d = deductible, r = r)) else deductible
u <- if (inflation) substitute(u/r, list(u = limit, r = r)) else limit
## The modified cdf or pdf are usually defined in branches. To
## avoid using nested ifelse(), we will rather rely on sets of
## expressions to make the required calculations for each branch
## separately. This is actually much faster.
##
## The output function will have varying number of said
## expressions depending on the case that is dealt with. We will
## build the body of the output function piece by piece as we go
## along.
e <- expression(Call <- match.call())
## One main discriminating factor is whether the cdf is needed for
## the output function of not.
if (needs.cdf)
{
## Get argument list of 'cdf' to transfert them to the output
## function.
argv <- formals(cdf) # arguments as list
argn <- names(argv) # arguments names as strings
## Remember if argument 'lower.tail' is available, so we can
## use it later. Then, drop unsupported arguments 'lower.tail'
## and 'log.p'.
has.lower <- "lower.tail" %in% argn
argn <- setdiff(argn, c("lower.tail", "log.p"))
## Calculations in the output function are done by evaluating
## function calls built upon the call to the output function
## itself. This is convenient as we do not need to fiddle with
## the values of the formal arguments.
if (is.cdf) # output function computes F(y)
{
## The output function will have the same formal arguments
## as the cdf. Object 'x' holds the symbol of the first
## argument.
argsFUN <- argv[argn] # arguments of output function
x <- as.name(argn[1]) # symbol
## Calls needed in this case:
## 1. one to compute F(y);
## 2. one to compute F(d) if there is a deductible;
## 3. one to compute 1 - F(d) for the per payment cases.
## Never need to compute F(u).
##
## If 'lower.tail' is available in 'cdf', then set it to
## FALSE to compute S(d) = 1 - F(d) more accurately.
e <- c(e,
quote(F <- Call), # 1.
substitute(F[[1L]] <- as.name(fun),
list(fun = as.character(Call$cdf))))
if (deductible)
{
e <- c(e,
quote(Fd <- F), # 2.
substitute(Fd[[2L]] <- a, list(a = d)))
if (!per.loss & has.lower)
e <- c(e,
quote(Sd <- Fd), # 3.
quote(Sd$lower.tail <- FALSE))
}
}
else # output function computes f(y)
{
## When there is a limit, we will need to compute 1 - F(u)
## as is or using 'lower.tail = FALSE' to improve
## accuracy. For clarity in the output function, we will
## use Fu as object name in the first case and Su in the
## second.
if (has.limit)
{
if (has.lower)
{
Fu.name <- as.name("Su")
Su.quote <- quote(eval.parent(Su))
}
else
{
Fu.name <- as.name("Fu")
Su.quote <- quote((1 - eval.parent(Fu)))
}
}
## Calls needed in this case:
## 1. one to compute F(d) if there is a deductible for the
## per loss cases;
## 2. one to compute 1 - F(d) if there is a deductible for the
## per payment cases;
## 3. one to compute 1 - F(u) when there is a limit.
## No function to compute F(y) needed.
##
## If 'lower.tail' is available in 'cdf', then set it to
## FALSE to compute S(d) = 1 - F(d) and S(u) = 1 - F(u)
## more accurately.
if (deductible) # f(y) with deductible
{
Fd.name <- as.name(if (!per.loss & has.lower) "Sd" else "Fd")
e <- c(e,
substitute(G <- Call, list(G = Fd.name)), # 1. or 2.
if (!per.loss & has.lower) quote(Sd$lower.tail <- FALSE),
substitute(G[[1L]] <- as.name(fun),
list(G = Fd.name, fun = as.character(Call$cdf))),
substitute(names(G)[2L] <- q, list(G = Fd.name, q = argn[1])),
substitute(G[[2L]] <- a, list(G = Fd.name, a = d)))
if (has.limit)
e <- c(e,
substitute(H <- G, list(H = Fu.name, G = Fd.name)), # 3.
if (per.loss & has.lower) quote(Su$lower.tail <- FALSE),
substitute(H[[2L]] <- a, list(H = Fu.name, a = u)))
}
else # f(y) with limit only
{
## Since 'needs.cdf == TRUE', then this case
## necessarily has 'limit < Inf'. Only call needed is
## one to compute 1 - F(u).
e <- c(e,
substitute(G <- Call, list(G = Fu.name)),
if (has.lower) quote(Su$lower.tail <- FALSE),
substitute(G[[1L]] <- as.name(fun),
list(G = Fu.name, fun = as.character(Call$cdf))),
substitute(names(G)[2L] <- q, list(G = Fu.name, q = argn[1])),
substitute(G[[2L]] <- a, list(G = Fu.name, a = u)))
}
}
}
## Repeat same steps as above for case needing the pdf. The output
## function is a pdf and in this case the arguments of the output
## function are those of 'pdf'.
if (!is.cdf)
{
argv <- formals(pdf) # arguments as list
argn <- setdiff(names(argv), "log") # drop argument 'log'
argsFUN <- argv[argn] # arguments of output function
x <- as.name(argn[1]) # symbol
e <- c(e,
quote(f <- Call),
substitute(f[[1L]] <- as.name(fun),
list(fun = as.character(Call$pdf))))
}
## Build the value at which the underlying pdf/cdf will be called
## for non special case values of 'x'. We need to index 'x' to
## only compute for the correct values of a given branch.
x.mod <- as.call(c(as.name("["), x, as.name("w")))
if (coinsurance < 1)
x.mod <- substitute(x/alpha, list(x = x.mod, alpha = coinsurance))
if (deductible & !franchise)
x.mod <- substitute(x + d, list(x = x.mod, d = deductible))
if (inflation)
x.mod <- substitute((x)/r, list(x = x.mod, r = r))
## Each pdf/cdf is defined in three branches. Define the
## boundaries and conditions for the first two branches. Those for
## the third branch are defined further down.
if (franchise)
{
bound1 <- coinsurance * deductible
bound2 <- coinsurance * limit
cond1 <- if (is.cdf)
substitute(0 <= x & x <= b1, list(x = x, b1 = bound1))
else
quote(x == 0)
cond2 <- substitute(b1 < x & x < b2,
list(x = x, b1 = bound1, b2 = bound2))
}
else
{
bound1 <- 0
bound2 <- coinsurance * (limit - deductible)
cond1 <- substitute(x == 0, list(x = x))
cond2 <- substitute(0 < x & x < b, list(x = x, b = bound2))
}
## Initialization of the results vector in the output function
## with 0s.
e <- c(e,
substitute(res <- numeric(length(x)), list(x = x)))
## Definition of the output function for the first branch. There
## is a computation to make only if there is a deductible with the
## payment per loss random variable. For all other cases, the
## value in the first branch is 0 and we rely on the
## initialization with numeric() done at the previous step.
if (per.loss & deductible)
e <- c(e,
substitute(res[which(cond1)] <- eval.parent(Fd),
list(cond1 = cond1)))
## Definition of the output function for the second and third
## branches. The 'is.cdf = TRUE' and 'is.cdf = FALSE' cases must
## be treated separately.
if (is.cdf)
{
cond3 <- substitute(x >= b, list(x = x, b = bound2))
f2 <- quote(eval.parent(F))
if (!per.loss & deductible)
f2 <- if (has.lower)
substitute((f - F)/S,
list(f = f2,
F = quote(eval.parent(Fd)),
S = quote(eval.parent(Sd))))
else
substitute((f - F)/S,
list(f = f2,
F = quote((p <- eval.parent(Fd))),
S = quote((1 - p))))
e <- c(e,
substitute(w <- which(cond), list(cond = cond2)),
substitute(F[[2L]] <- x, list(x = x.mod)),
substitute(res[w] <- f, list(f = f2)),
if (has.limit) substitute(res[cond] <- 1, list(cond = cond3)))
}
else
{
cond3 <- substitute(x == b, list(x = x, b = bound2))
f2 <- quote(eval.parent(f))
if (has.limit) f3 <- Su.quote
if (!per.loss & deductible)
{
if (has.limit)
{
f2 <- if (has.lower)
substitute(f/(p <- S),
list(f = f2, S = quote(eval.parent(Sd))))
else
substitute(f/(p <- S),
list(f = f2, S = quote(1 - eval.parent(Fd))))
f3 <- substitute(f/p, list(f = f3))
}
else
f2 <- if (has.lower)
substitute(f/S,
list(f = f2, S = quote(eval.parent(Sd))))
else
substitute(f/S,
list(f = f2, S = quote((1 - eval.parent(Fd)))))
}
if (inflation | coinsurance < 1)
f2 <- substitute(f/k, list(f = f2, k = coinsurance * r))
e <- c(e,
substitute(w <- which(cond), list(cond = cond2)),
substitute(f[[2L]] <- x, list(x = x.mod)),
substitute(res[w] <- f, list(f = f2)),
if (has.limit) substitute(res[cond] <- f, list(cond = cond3, f = f3)))
}
## Last expression of the output function.
e <- c(e, quote(res))
## Wrap up the output function.
FUN <- function() {}
body(FUN) <- as.call(c(as.name("{"), e)) # taken from help(body)
formals(FUN) <- argsFUN # set arguments
environment(FUN) <- new.env() # new, empty environment
FUN
}