# R/Integration.Steps.R In fptdApprox: Approximation of First-Passage-Time Densities for Diffusion Processes

```Integration.Steps <-
function (sfptl, variableStep = TRUE, from.t0 = FALSE, to.T = FALSE,
n = 250, p = 0.2, alpha = 1)
{
Call <- match.call()
if (!is.summary.fptl(sfptl))
stop(paste(sQuote("sfptl"), " object is not of class ",
shQuote("summary.fptl"), ".", sep = ""))
if (n <= 0)
stop("n <= 0.")
if (variableStep) {
if ((p <= 0) | (p > 1))
stop("p is not in (0,1].")
if ((alpha <= 0) | (alpha > 1))
stop("alpha is not in (0,1].")
cotes <- c(n = 50, p = 0.1, alpha = 0.75)
logic <- (c(n, p, alpha) < cotes)
if (any(logic)) {
text <- paste(names(cotes)[logic], cotes[logic],
sep = " >= ")
if (length(text) > 1)
text <- paste(paste(text[-length(text)], collapse = ", "),
text[length(text)], sep = " and ")
text <- paste(text, "is recommended. If not, some integration steps can be too large.")
cat(text, "\n")
repeat {
break
}
stop("\n\n")
else cat("\n")
}
}
else {
if (missing(p)) {
if (!missing(alpha))
cat("alpha argument is ignored.\n")
}
else {
if (missing(alpha))
cat("p argument is ignored.\n")
else cat("p and alpha arguments are ignored.\n")
}
}
if (length(sfptl) == 1L) {
t0 <- attr(sfptl, "FPTLCall")[[1]]\$t0
T <- attr(sfptl, "FPTLCall")[[1]]\$T
SFPTL <- sfptl[[1]]\$instants
e <- c(t(SFPTL[, c(2, 5)]))
m <- 2 * nrow(SFPTL) - 1
skip <- rep(c(FALSE, TRUE), len = m)
A <- diff(e)
a <- rep(SFPTL[, 3] - SFPTL[, 2], each = 2, len = m)
q <- A/a
if (any(skip))
q[skip] <- p * ifelse(q[skip] > 1, q[skip]^alpha,
q[skip])
H <- matrix(, nrow = m, ncol = 3)
H[, 1] <- e[-length(e)]
H[, 2] <- e[-1]
H[, 3] <- A/ceiling(n * q)
if (any(skip)) {
index <- which(skip)
logic <- H[index, 3] < H[index - 1L, 3]
if (any(logic)) {
J1 <- index[logic]
J2 <- J1 - 1L
H[J2, 2] <- H[J1, 2]
w <- H[J2, 2] - H[J2, 1]
H[J2, 3] <- w/ceiling(w/H[J2, 3])
H <- H[-J1, ]
e <- e[-J1]
skip <- skip[-J1]
}
}
if (from.t0 & (t0 < e[1])) {
A <- e[1] - t0
q <- A/a[1]
h <- A/ceiling(n * p * ifelse(q > 1, q^alpha, q))
if (h >= H[1, 3]) {
H <- rbind(c(t0, e[1], h), H)
skip <- c(FALSE, skip)
e <- c(t0, e)
}
else {
e[1] <- t0
H[1, 1] <- t0
w <- e[2] - t0
H[1, 3] <- w/ceiling(n * w/a[1])
}
}
if (to.T & (e[length(e)] < T)) {
A <- T - e[length(e)]
q <- A/a[length(a)]
h <- A/ceiling(n * p * ifelse(q > 1, q^alpha, q))
if (h >= H[nrow(H), 3]) {
H <- rbind(H, c(e[length(e)], T, h))
skip <- c(skip, FALSE)
e <- c(e, T)
}
else {
e[length(e)] <- T
H[nrow(H), 2] <- T
w <- T - e[length(e) - 1L]
H[1, 3] <- w/ceiling(n * w/a[length(a)])
}
}
skip <- as.list(skip)
}
else {
t0 <- attr(sfptl, "FPTLCall")[[1]]\$t0
T <- attr(sfptl, "FPTLCall")[[1]]\$T
SFPTL <- lapply(sfptl, "[[", "instants")
m <- sapply(SFPTL, nrow)
S <- lapply(2 * m - 1, rep_len, x = c(FALSE, TRUE))
E <- lapply(SFPTL, function(x) c(t(x[, c(2, 5)])))
a <- lapply(SFPTL, function(x) rep(x[, 3] - x[, 2], each = 2,
len = 2 * nrow(x) - 1))
e <- sort(unique(unlist(E)))
from.t0 <- from.t0 & (t0 < e[1])
to.T <- to.T & (e[length(e)] < T)
if (from.t0)
e <- c(t0, e)
if (to.T)
e <- c(e, T)
l1 <- (sapply(E, head, n = 1L) > t0)
l2 <- (sapply(E, tail, n = 1L) < T)
if (any(l1)) {
S[l1] <- lapply(S[l1], function(x) c(TRUE, x))
E[l1] <- lapply(E[l1], function(x, y) c(y, x), y = t0)
a[l1] <- lapply(a[l1], function(x) c(x[1], x))
}
if (any(l2)) {
S[l2] <- lapply(S[l2], function(x) c(x, TRUE))
E[l2] <- lapply(E[l2], function(x, y) c(x, y), y = T)
a[l2] <- lapply(a[l2], function(x) c(x, x[length(x)]))
}
A <- lapply(E, diff)
Q <- mapply("/", A, a, SIMPLIFY = FALSE)
logic <- sapply(S, any)
if (any(logic))
Q[logic] <- mapply(function(q, s, P, Alpha) {
q[s] <- P * ifelse(q[s] > 1, q[s]^Alpha, q[s])
q
}, Q[logic], S[logic], MoreArgs = list(P = p, Alpha = alpha),
SIMPLIFY = FALSE)
Q <- lapply(lapply(Q, "*", n), ceiling)
P <- mapply("/", A, Q, SIMPLIFY = FALSE)
J <- mapply(function(h, s) if (any(s)) {
i <- setdiff(which(s), 1L)
i[h[i] < h[i - 1]]
}
else integer(0), P, S, SIMPLIFY = FALSE)
logic <- (sapply(J, length) > 0L)
if (any(logic)) {
J <- J[logic]
P[logic] <- mapply(function(e, h, b, j, N) {
k <- j - 1L
w <- e[j + 1L] - e[k]
h[k] <- w/ceiling(N * w/b[j])
h <- h[-j]
h
}, E[logic], P[logic], a[logic], J, MoreArgs = list(N = n),
SIMPLIFY = FALSE)
E[logic] <- mapply(function(e, j) e[-j], E[logic],
J, SIMPLIFY = FALSE)
S[logic] <- mapply(function(s, j) s[-j], S[logic],
J, SIMPLIFY = FALSE)
e <- intersect(e, sort(unique(unlist(E))))
}
if (any(l1)) {
logic <- sapply(P[l1], function(h) h[1] < h[2])
if (any(logic)) {
J <- which(l1)[logic]
P[J] <- mapply(function(e, h, b, N) {
w <- e[3] - e[1]
h[1] <- w/ceiling(N * w/b[1])
h <- h[-2]
h
}, E[J], P[J], a[J], MoreArgs = list(N = n),
SIMPLIFY = FALSE)
E[J] <- lapply(E[J], function(e) e[-2])
S[J] <- lapply(S[J], function(s) s[-2])
e <- intersect(e, sort(unique(unlist(E))))
}
}
H <- matrix(, nrow = length(e) - 1, ncol = 3)
H[, 1] <- e[-length(e)]
H[, 2] <- e[-1]
index <- lapply(E, findInterval, x = H[, 1])
skip <- as.list(as.data.frame(t(mapply("[", S, index))))
P <- mapply("[", P, index)
H[, 3] <- apply(P, 1, min, na.rm = TRUE)
if (from.t0)
skip[[1]] <- rep(FALSE, length(SFPTL))
if (to.T)
skip[[length(skip)]] <- rep(FALSE, length(SFPTL))
d <- rev(which(c(0, diff(H[, 3])) == 0)[-1])
for (i in d) {
if (all(skip[[i - 1]] == skip[[i]])) {
H[i - 1, 2] <- H[i, 2]
H <- H[-i, ]
skip <- skip[-i]
}
}
}
w <- H[, 2] - H[, 1]
H[, 3] <- w/ceiling(w/H[, 3])
if (!variableStep) {
l1 <- sapply(skip, all)
l2 <- !l1
h <- min(H[l2, 3])
w <- H[nrow(H), 2] - H[1, 1]
h <- w/ceiling(w/h)
upper <- e[-1]
upper[l2] <- e[1] + h * ceiling((upper[l2] - e[1])/h)
upper[l1] <- e[1] + h * floor((upper[l1] - e[1])/h)
lower <- c(e[1], upper[-length(upper)])
H[, 1] <- lower
H[, 2] <- upper
H[, 3] <- h
l1 <- (upper == lower)
if (any(l1)) {
H <- H[!l1, ]
skip <- skip[!l1]
}
}
colnames(H) <- c("lower end", "upper end", "integration step")
rownames(H) <- paste("Subinterval", 1:nrow(H))
if (length(skip) > 1L)
names(skip) <- rownames(H)
out <- list(H = H, skip = skip)
return(out)
}
```

## Try the fptdApprox package in your browser

Any scripts or data that you put into this service are public.

fptdApprox documentation built on Nov. 2, 2023, 5:07 p.m.