## Not run:
# ## Reproduce Table 1 (p 101) of Loader and Deely (1987)
# ## define a vector of n values
# nLD <- c(8,16,32,64,128)
#
# ## Part 1: c(t) = sqrt(1+t) and tMax=1
# ## define cFct
# cFT1p1 <- function(t) sqrt(1+t)
# ## define the different bFct
# bFT1p1.ii <- function(t) 0.5/sqrt(1+t)
# bFT1p1.iii <- function(t) (cFT1p1(t)-cFT1p1(0))/t
# bFT1p1.iv <- function(t) 0.5*(bFT1p1.ii(t)+bFT1p1.iii(t))
# bFT1p1.v <- function(t) (2*t-4/5*((1+t)^2.5-1))/t^3+3*cFT1p1(t)/2/t
# ## Do the calculations
# round(t(sapply(nLD,
# function(n) {
# c(n=n,
# i=crossGeneral(tMax=1,h=1/n,cFct=cFT1p1)$G[n],
# ii=crossGeneral(tMax=1,h=1/n,cFct=cFT1p1,bFct=bFT1p1.ii)$G[n],
# iii=crossGeneral(tMax=1,h=1/n,cFct=cFT1p1,bFct=bFT1p1.iii)$G[n],
# iv=crossGeneral(tMax=1,h=1/n,cFct=cFT1p1,bFct=bFT1p1.iv)$G[n],
# v=crossGeneral(tMax=1,h=1/n,cFct=cFT1p1,bFct=bFT1p1.v)$G[n])})),
# digits=6)
#
# ## Part 2: c(t) = exp(-t) and tMax=1
# ## define cFct
# cFT1p2 <- function(t) exp(-t)
# ## define the different bFct
# cFT1p2 <- function(t) exp(-t)
# bFT1p2.ii <- function(t) -exp(-t)
# bFT1p2.iii <- function(t) (cFT1p2(t)-cFT1p2(0))/t
# bFT1p2.iv <- function(t) 0.5*(bFT1p2.ii(t)+bFT1p2.iii(t))
# bFT1p2.v <- function(t) 3*(1-t-exp(-t))/t^3+3*cFT1p2(t)/2/t
# ## Do the calculations
# round(t(sapply(nLD,
# function(n) {
# c(n=n,
# i=crossGeneral(tMax=1,h=1/n,cFct=cFT1p2)$G[n],
# ii=crossGeneral(tMax=1,h=1/n,cFct=cFT1p2,bFct=bFT1p2.ii)$G[n],
# iii=crossGeneral(tMax=1,h=1/n,cFct=cFT1p2,bFct=bFT1p2.iii)$G[n],
# iv=crossGeneral(tMax=1,h=1/n,cFct=cFT1p2,bFct=bFT1p2.iv)$G[n],
# v=crossGeneral(tMax=1,h=1/n,cFct=cFT1p2,bFct=bFT1p2.v)$G[n])})),
# digits=6)
#
# ## Part 3: c(t) = t^2 + 3*t + 1 and tMax=1
# ## define cFct
# cFT1p3 <- function(t) t^2+3*t+1
# ## define the different bFct
# bFT1p3.ii <- function(t) 2*t+3
# bFT1p3.iii <- function(t) (cFT1p3(t)-cFT1p3(0))/t
# bFT1p3.v <- function(t) 5*t/4+3
# bFT1p3.vi <- function(t) rep(3,length(t))
# round(t(sapply(nLD,
# function(n) {
# c(n=n,
# i=crossGeneral(tMax=1,h=1/n,cFct=cFT1p3)$G[n],
# ii=crossGeneral(tMax=1,h=1/n,cFct=cFT1p3,bFct=bFT1p3.ii)$G[n],
# iii=crossGeneral(tMax=1,h=1/n,cFct=cFT1p3,bFct=bFT1p3.iii)$G[n],
# v=crossGeneral(tMax=1,h=1/n,cFct=cFT1p3,bFct=bFT1p3.v)$G[n],
# vi=crossGeneral(tMax=1,h=1/n,cFct=cFT1p3,bFct=bFT1p3.vi)$G[n])})),
# digits=6)
#
# ## Part 3: c(t) = t^2 + 3*t + 1 and tMax=1
# ## define cFct
# cFT1p4 <- function(t) 1/t
# ## Here only column (i) and (vii) are reproduced.
# ## If one attempts to reproduce (ii) directly with crossGeneral
# ## a NaN appears (when a -Inf is the correct value) in functions
# ## F() and K(,) (see details) for instance when t=0 in F.
# ## Then as crossGeneral is presently written R attempts to
# ## compute t*b(t), where b(t) is c'(t), that is, t*(-1/t^2) which is
# ## NaN (for R) when t=0.
# bFT1p4.vii <- function(t) rep(-1,length(t))
# round(t(sapply(nLD,
# function(n) {
# c(n=n,
# i=crossGeneral(tMax=1,h=1/n,cFct=cFT1p4)$G[n],
# vii=crossGeneral(tMax=1,h=1/n,cFct=cFT1p4,bFct=bFT1p4.vii)$G[n])})),
# digits=6)
# ## The last 3 rows of column (vii) are not the same as in the paper
#
# ## Reproduce Table 4 (p 104) of Loader and Deely (1987)
# ## As before the probability of first passage between
# ## 0 and 1 is computed. This time only three boundary
# ## functions are considered. The error bounds are
# ## obtained
#
# ## Part 1: c(t) = sqrt(1+t)
# ## Left columns pair: b(t) = c'(t)
# round(t(sapply(nLD,
# function(n) {
# res <- crossGeneral(tMax=1,h=1/n,cFct=cFT1p1,bFct=bFT1p1.ii,withBounds=TRUE,Lplus=TRUE)
# c(Gl=res$Gl[n],Gu=res$Gu[n])
# }
# )
# ),
# digits=5)
#
# ## Right columns pair: b(t) = 0
# round(t(sapply(nLD,
# function(n) {
# res <- crossGeneral(tMax=1,h=1/n,cFct=cFT1p1,withBounds=TRUE,Lplus=TRUE)
# c(n=n,Gl=res$Gl[n],Gu=res$Gu[n])
# }
# )
# ),
# digits=5)
#
# ## Part 2: c(t) = t^2 + 3*t + 1
# ## Left columns pair: b(t) = 3 - 2*t
# round(t(sapply(nLD,
# function(n) {
# res <- crossGeneral(tMax=1,h=1/n,cFct=cFT1p3,bFct=function(t) 3-2*t,withBounds=TRUE,Lplus=TRUE)
# c(n=n,Gl=res$Gl[n],Gu=res$Gu[n])
# }
# )
# ),
# digits=5)
#
# ## Right columns pair: b(t) = 3 - t
# round(t(sapply(nLD,
# function(n) {
# res <- crossGeneral(tMax=1,h=1/n,cFct=cFT1p3,bFct=function(t) 3-2*t,withBounds=TRUE,Lplus=TRUE)
# c(n=n,Gl=res$Gl[n],Gu=res$Gu[n])
# }
# )
# ),
# digits=5)
#
# ## Part 3: c(t) = 1 + sin(t)
# ## Left columns pair: b(t) = c'(t)
# round(t(sapply(nLD,
# function(n) {
# res <- crossGeneral(tMax=1,h=1/n,cFct=function(t) 1+sin(t),bFct=function(t) cos(t),withBounds=TRUE,Lplus=TRUE)
# c(n=n,Gl=res$Gl[n],Gu=res$Gu[n])
# }
# )
# ),
# digits=5)
#
# ## Left columns pair: b(t) = 0.5
# round(t(sapply(nLD,
# function(n) {
# res <- crossGeneral(tMax=1,h=1/n,cFct=function(t) 1+sin(t),bFct=function(t) rep(0.5,length(t)),withBounds=TRUE,Lplus=TRUE)
# c(n=n,Gl=res$Gl[n],Gu=res$Gu[n])
# }
# )
# ),
# digits=5)
#
#
# ## Check crossGeneral against an analytical inverse Gaussian
# ## distribution
# ## Define inverse Gaussian parameters
# mu.true <- 0.075
# sigma2.true <- 3
# ## Define a function transforming the "drift" (mu.true) and
# ## "noise variance" (sigma2.true) of the default inverse
# ## Gaussian parametrization of STAR into a
# ## linear boundary of an equivalent Wiener process first
# ## passage time problem
# star2ld <- function(mu,sigma2) c(sqrt(1/sigma2),-sqrt(1/sigma2)/mu)
# ## Get the "equivalent" boundary parameters (y intercept and slope)
# parB1 <- star2ld(mu.true,sigma2.true)
# ## Plot the "target" inverse Gaussian density
# xx <- seq(0.001,0.3,0.001)
# plot(xx,dinvgauss(xx,mu=mu.true,sigma2=sigma2.true),type="l")
# ## Get the numerical estimate of the density using Loader and
# ## Deely Volterra integral equation method
# igB1 <- crossGeneral(tMax=0.3,h=0.001,cFct=function(t) parB1[1]+parB1[2]*t,withBounds=FALSE)
# ## superpose the numerical estimate to the exact solution
# ## use lines method to do that
# lines(igB1,"density",col=2)
#
# ## Use of crossTight and associated function
# ## Get the paramters, a and b, of the "approximate"
# ## tightest boundary: c(t) = a + b*sqrt(t), giving a
# ## 0.05 probability of exit between 0 and 1
# ## (in fact we are discussing here a pair of symmetrical
# ## bounaries, c(t) and -c(t)). See Kendall et al (2007)
# ## for details
# ## Start by defining the target function
# target95 <- mkTightBMtargetFct(ci=0.95)
# ## get the optimal log(a) and log(b) using
# ## the values of table 1 of Kendall et al as initial
# ## guesses
# p95 <- optim(log(c(0.3,2.35)),target95,method="BFGS")
# ## check the convergence of BFGS
# p95$convergence
# ## check if the parameters changed a lot
# exp(p95$par)
# ## Get the bounds on G(1) for these optimal parameters
# d95 <- crossTight(a=exp(p95$par[1]),b=exp(p95$par[2]),withBound=TRUE,logScale=FALSE)
# ## print out the summary
# summary(d95)
# ## Do the same for the 0.01 probability of first passage
# target99 <- mkTightBMtargetFct(ci=0.99)
# p99 <- optim(p95$par,target99,method="BFGS")
# p99$convergence
# exp(p99$par)
# d99 <- crossTight(a=exp(p99$par[1]),b=exp(p99$par[2]),withBound=TRUE,logScale=FALSE)
# summary(d99)
# ## End(Not run)
Run the code above in your browser using DataLab