R/functions.R

Defines functions qFMstable qEstable standardStableQuantile tailsGstable lnorm.param ImpliedVol BSOptionValue optionsFMstable callFMstable putFMstable pFMstable.alpha0 moments setParam iidcombine matchQuartiles fitGivenQuantile pEstable dEstable pFMstable dFMstable tailsFMstable tailsEstable bigAsInf print.stableParameters setMomentsFMstable setMomentsFMstable

Documented in BSOptionValue callFMstable dEstable dFMstable fitGivenQuantile iidcombine ImpliedVol lnorm.param matchQuartiles moments optionsFMstable pEstable pFMstable pFMstable.alpha0 print.stableParameters putFMstable qEstable qFMstable setMomentsFMstable setParam tailsEstable tailsFMstable tailsGstable

# Find parameters corresponding to mean and sd, given alpha
setMomentsFMstable <- function(mean=1, sd=1,
    alpha, oneminusalpha, twominusalpha){
  if(!missing(alpha)){
    # Case where alpha is specified
    if(length(alpha) != 1) stop(
      "setMomentsFMstable: alpha must be of length 1")
    if(!missing(oneminusalpha) || !missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    oneminusalpha <- 1 - alpha
    twominusalpha <- 2 - alpha
  } else if(!missing(oneminusalpha)){
    # Case where oneminusalpha is specified
    if(length(oneminusalpha) != 1) stop(
      "setMomentsFMstable: oneminusalpha must be of length 1")
    if(!missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    alpha <- 1 - oneminusalpha
    twominusalpha <- 1 + oneminusalpha
  } else {
    if(missing(twominusalpha))stop(
      "Specify one of alpha, oneminusalpha, twominusalpha")
    # Case where twominusalpha is specified
    if(length(twominusalpha) != 1) stop(
      "setMomentsFMstable: twominusalpha must be of length 1")
    alpha <- 2 - twominusalpha
    oneminusalpha <- twominusalpha - 1
  }
  if(length(mean) != 1) stop("setMomentsFMstable: mean must be of length 1")
  if(length(sd) != 1) stop("setMomentsFMstable: sd must be of length 1")
  if(alpha <= 0 || twominusalpha < 0) stop(
    "setMomentsFMstable: alpha must be >= 0 and <= 2")

  cv <- sd/mean
  logmr <- log1p(cv*cv)

  # Case when alpha is less than 0.5 
  if(alpha < 0.5){
    scale.to.alpha <- logmr/(2- 2^alpha)
    logscale <- log(scale.to.alpha)/alpha
    location <- -log(mean) - scale.to.alpha
  } else if(oneminusalpha == 0){
    # Case when alpha is precisely 1
    scale <- pi*logmr/(4*log(2))
    logscale <- log(scale)
    location <- logscale*2*scale/pi - log(mean)
  } else{
    # Case when alpha exceeds 0.5 but is not precisely 1
    s <- if(alpha < 1.5) sin(oneminusalpha*.5*pi) else -cos(twominusalpha*.5*pi)
    sinalpha <- if(alpha < 1) sin(alpha*.5*pi) else sin(twominusalpha*.5*pi)
    scale.to.alpha <- -.5*logmr*s/expm1(-oneminusalpha*log(2))
    logscale <- log(scale.to.alpha)/alpha
    scale <- exp(logscale)
    location <- (scale*sinalpha - scale.to.alpha)/s - log(mean)
  }

  res <- list(alpha=alpha, oneminusalpha=oneminusalpha,
    twominusalpha=twominusalpha, location=location, logscale=logscale,
    created.by="setMomentsFMstable")
  class(res) <- "stableParameters"
  return(res)
}
# Find parameters corresponding to mean and sd, given alpha
setMomentsFMstable <- function(mean=1, sd=1,
    alpha, oneminusalpha, twominusalpha){
  if(!missing(alpha)){
    # Case where alpha is specified
    if(length(alpha) != 1) stop(
      "setMomentsFMstable: alpha must be of length 1")
    if(!missing(oneminusalpha) || !missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    oneminusalpha <- 1 - alpha
    twominusalpha <- 2 - alpha
  } else if(!missing(oneminusalpha)){
    # Case where oneminusalpha is specified
    if(length(oneminusalpha) != 1) stop(
      "setMomentsFMstable: oneminusalpha must be of length 1")
    if(!missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    alpha <- 1 - oneminusalpha
    twominusalpha <- 1 + oneminusalpha
  } else {
    if(missing(twominusalpha))stop(
      "Specify one of alpha, oneminusalpha, twominusalpha")
    # Case where twominusalpha is specified
    if(length(twominusalpha) != 1) stop(
      "setMomentsFMstable: twominusalpha must be of length 1")
    alpha <- 2 - twominusalpha
    oneminusalpha <- twominusalpha - 1
  }
  if(length(mean) != 1) stop("setMomentsFMstable: mean must be of length 1")
  if(length(sd) != 1) stop("setMomentsFMstable: sd must be of length 1")
  if(alpha <= 0 || twominusalpha < 0) stop(
    "setMomentsFMstable: alpha must be >= 0 and <= 2")

  cv <- sd/mean
  logmr <- log1p(cv*cv)

  # Case when alpha is less than 0.5 
  if(alpha < 0.5){
    scale.to.alpha <- logmr/(2- 2^alpha)
    logscale <- log(scale.to.alpha)/alpha
    location <- -log(mean) - scale.to.alpha
  } else if(oneminusalpha == 0){
    # Case when alpha is precisely 1
    scale <- pi*logmr/(4*log(2))
    logscale <- log(scale)
    location <- logscale*2*scale/pi - log(mean)
  } else{
    # Case when alpha exceeds 0.5 but is not precisely 1
    s <- if(alpha < 1.5) sin(oneminusalpha*.5*pi) else -cos(twominusalpha*.5*pi)
    sinalpha <- if(alpha < 1) sin(alpha*.5*pi) else sin(twominusalpha*.5*pi)
    scale.to.alpha <- -.5*logmr*s/expm1(-oneminusalpha*log(2))
    logscale <- log(scale.to.alpha)/alpha
    scale <- exp(logscale)
    location <- (scale*sinalpha - scale.to.alpha)/s - log(mean)
  }

  res <- list(alpha=alpha, oneminusalpha=oneminusalpha,
    twominusalpha=twominusalpha, location=location, logscale=logscale,
    created.by="setMomentsFMstable")
  class(res) <- "stableParameters"
  return(res)
}
if(FALSE){
yearDsn <- setMomentsFMstable(mean=10, sd=5, alpha=1.7)
# Compute mean by integration
f <- function(x) dFMstable(x, yearDsn) * x
integrate(f, lower=-Inf, upper=Inf, rel.tol=1.e-12)$value -10
f <- function(x) dFMstable(x, yearDsn) * x*x
integrate(f, lower=-Inf, upper=Inf, rel.tol=1.e-6)$value -125

yearDsn <- setMomentsFMstable(mean=10, sd=5, alpha=1.000001)
yearDsn
yearDsn <- setMomentsFMstable(mean=10, sd=5, alpha=.999999)
yearDsn
yearDsn <- setMomentsFMstable(mean=10, sd=5, alpha=1)
yearDsn
}

print.stableParameters <- function(x, ...){
  cat(paste(" ******** Stable parameter object created by", x$created.by))
  cat(paste("\n ** Alpha =", x$alpha,"  location =", x$location,
    "  logscale =", x$logscale))
  if(abs(x$oneminusalpha) < .1) cat(paste("\n ** oneminusalpha =",
    x$oneminusalpha))
  if(abs(x$twominusalpha) < .1) cat(paste("\n ** twominusalpha =",
    x$twominusalpha))
  if(max(abs(1:2 -x$alpha - c(x$oneminusalpha, x$twominusalpha))) >
    .Machine$double.eps) cat(paste("\n Versions of alpha not consistent"))
  m <- moments(1:2, x)
  cat(paste("\n ** Logstable distribution has mean=", m[1],
    "  sd=", sqrt(m[2] - m[1]^2),"\n ********\n"))
}

# Replace numbers outside range -1.e-308 to +1.e308 with -Inf or Inf
bigAsInf <- function(x){
  x[x > 1.e308] <- Inf
  x[x < -1.e308] <- -Inf
  return(x)
}

tailsEstable <- function(x, stableParamObj){
  if(class(stableParamObj) != "stableParameters")stop(paste("tailsEstable:",
    "Parameter stableParamObj must be of class stableParameters"))
  storage.mode(x) <- "double"
  n <- as.double(length(x))
  temp <- .C("RtailsMSS", stableParamObj$alpha, stableParamObj$oneminusalpha,
    stableParamObj$twominusalpha, stableParamObj$location,
    stableParamObj$logscale, n, x, out1=double(n), out2=double(n),
    out3=double(n), out4=double(n), out5=double(n), out6=double(n),
    COPY=rep(c(FALSE, TRUE), c(7, 6)), PACKAGE="FMStable")
  list(density=bigAsInf(temp$out1), F=bigAsInf(temp$out3),
    righttail=bigAsInf(temp$out5), logdensity=bigAsInf(temp$out2),
    logF=bigAsInf(temp$out4), logrighttail=bigAsInf(temp$out6))
}

# Given a set of x values for a log
#  maximally skew stable distribution, find density,
#  distribution function, and right tail probability.
tailsFMstable <- function(x, stableParamObj){
  if(class(stableParamObj) != "stableParameters")stop(paste("tailsFMstable:",
    "Parameter stableParamObj must be of class stableParameters"))
  storage.mode(x) <- "double"
  # Send finite positive x values to C function
  forC <- is.finite(x) & x>0
  n <- as.double(sum(forC))
  if(n > 0) temp <- .C("RtailslogMSS", stableParamObj$alpha,
    stableParamObj$oneminusalpha, stableParamObj$twominusalpha,
    stableParamObj$location, stableParamObj$logscale, n, x[forC],
    out1=double(n), out2=double(n), out3=double(n), out4=double(n),
    out5=double(n), out6=double(n), COPY=rep(c(FALSE, TRUE), c(7, 6)),
    PACKAGE="FMStable")
  nx <- length(x)
  righttail <- F <- density <- double(nx)
  logrighttail <- logF <- logdensity <- rep(-Inf, nx)
  nas <- !is.finite(x)
  if(any(nas)){
    righttail[nas] <- F[nas] <- density[nas] <- NA
    logrighttail[nas] <- logF[nas] <- logdensity[nas] <- NA
  }
  le0 <- is.finite(x) & x <= 0
  if(any(le0)){
    righttail[le0] <- 1
    logrighttail[le0] <- F[le0] <- density[le0] <- 0
    logF[le0] <- logdensity[le0] <- -Inf
  }
  if(n > 0){
    density[forC] <- temp$out1
    F[forC] <- temp$out3
    righttail[forC] <- temp$out5
    logdensity[forC] <- temp$out2
    logF[forC] <- temp$out4
    logrighttail[forC] <- temp$out6
  }
  list(density=density, F=F, righttail=righttail, logdensity=logdensity,
    logF=logF, logrighttail=logrighttail)
}

dFMstable <- function(x, stableParamObj, log=FALSE){
  tls <- tailsFMstable(x, stableParamObj)
  if(log) tls$logdensity else tls$density
}

pFMstable <- function(x, stableParamObj, log=FALSE, lower.tail=TRUE){
  tls <- tailsFMstable(x, stableParamObj)
  if(log){ if(lower.tail) tls$logF else tls$logrighttail } else {
    if(lower.tail) tls$F else tls$righttail }
}

dEstable <- function(x, stableParamObj, log=FALSE){
  tls <- tailsEstable(x, stableParamObj)
  if(log) tls$logdensity else tls$density
}

pEstable <- function(x, stableParamObj, log=FALSE, lower.tail=TRUE){
  tls <- tailsEstable(x, stableParamObj)
  if(log){ if(lower.tail) tls$logF else tls$logrighttail } else {
    if(lower.tail) tls$F else tls$righttail }
}

# Find alpha for FMStable distribution by specifying mean,
#  standard deviation and probability of exceeding a value
fitGivenQuantile <- function(mean, sd, prob, value, tol=1.e-10){
  if(length(mean) != 1) stop("fitGivenQuantile: mean must be of length 1")
  if(length(sd) != 1) stop("fitGivenQuantile: sd must be of length 1")
  if(length(prob) != 1) stop("fitGivenQuantile: prob must be of length 1")
  if(length(value) != 1) stop("fitGivenQuantile: value must be of length 1")
  if(sd <= 0) stop("alpha.given.quantile: sd must be > 0")
  if(prob <= 0 || prob >= 1) stop("
    alpha.given.quantile: prob must be > 0 and < 1")
  if(value <= 0) stop("alpha.given.quantile: value must be > 0")
  if(mean <= 0) stop("alpha.given.quantile: mean must be > 0")
  # Check whether alpha=0 distribution would have sd larger than specified
  q <- 1-prob
  A <- mean/q
  minvar <- prob*mean*mean + q*(A-mean)^2
  if(sd*sd <= minvar) stop(paste("fitGivenQuantile: Quantile cannot be fitted.",
    "\nDistribution with probability", prob, "at zero and", q, "at", A,
    "\nhas standard deviation", sqrt(minvar), "which exceeds specified sd."))
  cv <- sd/mean
  Pdiscrepancy <- function(x) pFMstable(value/mean,
    setMomentsFMstable(mean=1, sd=cv, twominusalpha=x)) - prob
  xx <- uniroot(Pdiscrepancy, lower=0, upper=2-1.e-15, tol=tol)$root
  setMomentsFMstable(mean=mean, sd=mean*cv, twominusalpha=xx)
}

# Given alpha, find parameters of logstable distribution to matchQuartiles
matchQuartiles <- function(quartiles, alpha, oneminusalpha, twominusalpha,
    tol=1.e-10){
  if(!missing(alpha)){
    # Case where alpha is specified
    if(length(alpha) != 1) stop(
      "setMomentsFMstable: alpha must be of length 1")
    if(!missing(oneminusalpha) || !missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    case <- 1
  } else if(!missing(oneminusalpha)){
    # Case where oneminusalpha is specified
    if(length(oneminusalpha) != 1) stop(
      "setMomentsFMstable: oneminusalpha must be of length 1")
    if(!missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    case <- 2
  } else {
    if(missing(twominusalpha))stop(
      "Specify one of alpha, oneminusalpha, twominusalpha")
    # Case where twominusalpha is specified
    if(length(twominusalpha) != 1) stop(
      "setMomentsFMstable: twominusalpha must be of length 1")
    case <- 3
  }
  if(length(quartiles) != 2) stop(
    "matchQuartiles: quartiles must be of length 2")
  if(quartiles[1] <= 0 || quartiles[2] <= quartiles[1]) stop(
    "matchQuartiles: Require 0 < quartile[1] < quartile[2]")

  # First approximation to mean and sd
  av <- mean(quartiles)
  sd <- diff(quartiles)*.74
  switch(case, {
    disc <- function(param) sum((pFMstable(quartiles, setMomentsFMstable(
      exp(param[1]), exp(param[2]), alpha=alpha)) - c(.25,.75))^2)
  },{
    disc <- function(param) sum((pFMstable(quartiles, setMomentsFMstable(
      exp(param[1]), exp(param[2]), oneminusalpha=oneminusalpha)) -
      c(.25,.75))^2)
  },{
    disc <- function(param) sum((pFMstable(quartiles, setMomentsFMstable(
      exp(param[1]), exp(param[2]), twominusalpha=twominusalpha)) -
      c(.25,.75))^2)
  })

  # Find mean and sd that give match to specified quartiles
  opt <- optim(log(c(av, sd)), disc, control=list(reltol=tol*tol))
  if(opt$convergence > 0) print("matchQuartiles: Poor convergence")
  av <- exp(opt$par[1])
  sd <- exp(opt$par[2])
  switch(case,
    setMomentsFMstable(mean=av, sd, alpha=alpha),
    setMomentsFMstable(mean=av, sd, oneminusalpha=oneminusalpha),
    setMomentsFMstable(mean=av, sd, twominusalpha=twominusalpha))
}

# Find location and scale for convolution of n (not necessarily integral)
#  log maximally skew stable distributions.
# New moment generating function is original to power n
iidcombine <- function(n, stableParamObj){
  if(length(n) != 1) stop("iidcombine: n must be of length 1")
  if(class(stableParamObj) != "stableParameters") stop(paste("iidcombine:",
    "Parameter stableParamObj must be of class stableParameters"))
  if(n <=0) stop("iidcombine: n must be > 0")

  logscale.conv <- stableParamObj$logscale +log(n)/stableParamObj$alpha
  # For C parametrization, location parameters add and scales^(alpha) add
  if(stableParamObj$alpha < .5){
    logscale.conv <- stableParamObj$logscale +log(n)/stableParamObj$alpha
    location.conv <- n*stableParamObj$location
  } else {
    scale <- exp(stableParamObj$logscale)
    if (stableParamObj$oneminusalpha == 0){
      location.conv <- n*stableParamObj$location - 2/pi*(
         n*stableParamObj$logscale*scale -
              exp(logscale.conv)*logscale.conv)
    } else if (abs(stableParamObj$oneminusalpha) < 0.1) {
      location.conv <- n*stableParamObj$location +
        (exp(logscale.conv)-n*scale)/tan(.5*pi*stableParamObj$oneminusalpha)
    } else {
      tana <- if(stableParamObj$oneminusalpha > 0){
            tan(.5*pi*stableParamObj$alpha)
          } else { -tan(.5*pi*stableParamObj$twominusalpha)}
      location.conv <- n*stableParamObj$location +
        (exp(logscale.conv)-n*scale)*tana
    }
  }
  res <- list(alpha=stableParamObj$alpha,
    oneminusalpha=stableParamObj$oneminusalpha,
    twominusalpha=stableParamObj$twominusalpha, location=location.conv,
    logscale=logscale.conv, created.by="iidcombine")
  class(res) <- "stableParameters"
  return(res)
}

# Set parameters of stable distribution using any of three parametrizations
setParam <- function(alpha, oneminusalpha, twominusalpha,
    location, logscale, pm){
  if(!missing(alpha)){
    # Case where alpha is specified
    if(length(alpha) != 1) stop(
      "setParam: alpha must be of length 1")
    if(!missing(oneminusalpha) || !missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    oneminusalpha <- 1 - alpha
    twominusalpha <- 2 - alpha
  } else if(!missing(oneminusalpha)){
    # Case where oneminusalpha is specified
    if(length(oneminusalpha) != 1) stop(
      "setParam: oneminusalpha must be of length 1")
    if(!missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    alpha <- 1 - oneminusalpha
    twominusalpha <- 1 + oneminusalpha
  } else {
    if(missing(twominusalpha))stop(
      "Specify one of alpha, oneminusalpha, twominusalpha")
    # Case where twominusalpha is specified
    if(length(twominusalpha) != 1) stop(
      "setParam: twominusalpha must be of length 1")
    alpha <- 2 - twominusalpha
    oneminusalpha <- twominusalpha - 1
  }
  if(length(location) != 1) stop("setParam: location must be of length 1")
  if(length(logscale) != 1) stop("setParam: logscale must be of length 1")
  if(length(pm) != 1) stop("setParam: pm must be of length 1")
  pmcode <- rep(0:2, each=3)[ match(as.character(pm),
    c("0","S0","M", "1","S1","A", "2","CMS","C"))]
  if(is.na(pmcode)) stop(paste("setParam:  Parameterization must be",
    "specified as 0, 1 or 2; S0, S1 or CMS; or M, A or C"))
  if(alpha <= 0 || alpha > 2) stop(
    "setParam: Must have 0 < alpha <= 2")
  if(pmcode != 0)if(abs(alpha - 1) < .01)stop(
    "setParam: Only S0=M parametrization suitable near alpha=1")
  if(pmcode == 0)if(alpha < .1)stop(
    "setParam: S0=M parametrization not suitable for small alpha")
  angle <- 0.5*pi*alpha
  if(alpha < .5){
    if(pmcode == 0) {
      location <- location-tan(angle)*exp(logscale)
      logscale <- logscale-log(cos(angle))/alpha
    }
    if(pmcode == 1) logscale <- logscale-log(cos(angle))/alpha
  } else {
    tanangle <- if(alpha > 1) -tan(.5*pi*twominusalpha) else tan(angle)
    if(pmcode == 1) location <- location+tanangle*exp(logscale)
    if(pmcode == 2) {
      logscale <- logscale+log(abs(sin(.5*pi*oneminusalpha)))/alpha
      location <- location+tanangle*exp(logscale)
    }
  }
  res <- list(alpha=alpha, oneminusalpha=oneminusalpha,
    twominusalpha=twominusalpha, location=location, logscale=logscale,
    created.by=paste("setParam with pmcode =",pmcode))
  class(res) <- "stableParameters"
  return(res)
}

moments <- function(powers, stableParamObj, log=FALSE){
  if(class(stableParamObj) != "stableParameters") stop(paste("moments:",
    "Parameter stableParamObj must be of class stableParameters"))
  if(any(powers <= 0)) stop("moments: Must have powers > 0")
  logpowers <- log(powers)
  logprod <- logpowers + stableParamObj$logscale
  if(stableParamObj$alpha < 0.5){	# C parametrization
    # powers*location-(scale*powers)^alpha
    # If alpha is very small (e.g. .01) then scale might be 10^400?
    logmoment <- -powers * stableParamObj$location -
      exp(stableParamObj$alpha*logprod)
  } else {				# M = S0 parametrization
    prod <- exp(logprod)
    if(stableParamObj$oneminusalpha ==0){
      logmoment <- -powers * stableParamObj$location + logprod * 2*prod/pi
    } else {
      # Case when alpha > 0.5 but not precisely 1.
      sineEpsTerm <- sin(.5*pi*stableParamObj$oneminusalpha)
      oneminussineAlphaTerm <- 2 * sin(.25*pi*stableParamObj$oneminusalpha)^2
      logmoment <- -powers*stableParamObj$location - prod *
        (expm1(-logprod*stableParamObj$oneminusalpha)+oneminussineAlphaTerm)/
        sineEpsTerm
    }
  }
  if(log) logmoment else exp(logmoment)
}

# Compute distribution function of logstable distribution for alpha=0
pFMstable.alpha0 <- function(x, mean=1, sd=1, lower.tail=TRUE){
  if(length(mean) != 1) stop("pFMstable.alpha0: mean must be of length 1")
  if(length(sd) != 1) stop("pFMstable.alpha0: sd must be of length 1")
  # Probability of value A=mean+sd^2/mean is pA=mean/A.  Otherwise zero.
  A <- mean + sd^2/mean
  pA <- mean/A
  res <- double(length(x))
  if(lower.tail){
    res[x > 0] <- 1-pA
    res[x >= A] <- 1
  } else {
    res[x < A] <- pA
    res[x <= 0] <- 1
  }
  return(res)
}

# Value of put options under finite moment logstable distribution
putFMstable <- function(strike, paramObj, rel.tol=1.e-10){
  if(!all(is.finite(strike) & strike>0)) stop(
    "putFMstable: All values of strike must be finite and strictly positive")
  result <- double(length(strike))
  f <- function(x) pFMstable(x, paramObj)
  o <- order(strike)
  sum <- 0
  upto <- 0
  for (i in 1:length(strike)){
    s <- strike[o[i]]
    sum <- sum + integrate(f, lower=upto, upper=s, rel.tol=rel.tol)$value
    result[o[i]] <- sum
    upto <- s
  }
  lastput <- result[o[length(strike)]]
  #if(lastput+moments(1, paramObj)-s < .01 * lastput) cat(paste(
  #  "putFMstable: For large strikes, consider using callFMstable",
  #  "and the relationship\nput = call + strike - mean\n"))
  result
}

# Value of call options under finite moment logstable distribution
callFMstable <- function(strike, paramObj, rel.tol=1.e-10){
  if(!all(is.finite(strike) & strike>0)) stop(
    "callFMstable: All values of strike must be finite and strictly positive")
  result <- double(length(strike))
  f <- function(x) pFMstable(x, paramObj, lower.tail=FALSE)
  o <- order(strike, decreasing=TRUE)
  sum <- 0

  # Find where right tail is rel.tol^2 smaller than for largest strike
  pright <- pFMstable(strike[o[1]], paramObj, lower.tail=FALSE)
  if(pright <= 0) upto <- strike[o[1]] else upto <- qFMstable(
    max(min(.5, pright)*rel.tol^2, 1.e-300), paramObj, lower.tail=FALSE)

  # Work through strikes in descending order
  for (i in 1:length(strike)){
    s <- strike[o[i]]
    sum <- sum + integrate(f, lower=s, upper=upto, rel.tol=rel.tol)$value
    result[o[i]] <- sum
    upto <- s
  }
  lastcall <- result[o[length(strike)]]
  #if(lastcall-moments(1, paramObj)+s < .01 * lastcall) cat(paste(
  #  "callFMstable: For small strikes, consider using putFMstable",
  #  "and the relationship\ncall = put + mean - strike\n"))
  result
}
# Value of call and put options under finite moment logstable distribution
optionsFMstable <- function(strike, paramObj, rel.tol=1.e-10){
  if(!all(is.finite(strike) & strike>0)) stop(
    "callFMstable: All values of strike must be finite and strictly positive")
  call <- put <- double(length(strike))
  mean <- moments(1, paramObj)
  low <- strike < mean
  high <- !low
  if(any(low)){
    results <- putFMstable(strike[low], paramObj, rel.tol)
    put[low] <- results
    call[low] <- results + mean - strike[low]
  }
  if(any(high)){
    results <- callFMstable(strike[high], paramObj, rel.tol)
    call[high] <- results
    put[high] <- results + strike[high] - mean
  }
  return(list(put=put, call=call))
}

# Black-Scholes model for value of option
BSOptionValue <- function(spot, strike, expiry, volatility,
    intRate=0, carryCost=0, Call=TRUE){
  d1 <- (log(spot/strike)+(carryCost+.5*volatility*volatility)*expiry)/
    (volatility*sqrt(expiry))
  d2 <- d1 - volatility*sqrt(expiry)
  price <- if(Call){ spot*exp((carryCost-intRate)*expiry)*pnorm(d1) -
      strike*exp(-intRate*expiry)* pnorm(d2)
    } else {strike*exp(-intRate*expiry)*pnorm(-d2) -
      spot*exp((carryCost - intRate)*expiry)*pnorm(-d1)}
  return(price)
}
# tests
if(FALSE){
BSOptionValue(5, 5.1, .3, .2, Call=FALSE)
BSOptionValue(5, 5.1, .3, .2)
BSOptionValue(5, 5.2, .3, .2)
BSOptionValue(5, 5.3, .3, .2)
BSOptionValue(5, 5, .4, .3)
}

# Find implied volatilities to make prices consistent with Black-Scholes model
# May give vectors of strike, expiry and price
ImpliedVol <- function(spot, strike, expiry, price, intRate=0, carryCost=0,
    Call=TRUE, ImpliedVolLowerBound=.01, ImpliedVolUpperBound=1, tol=1.e-9){
  Nprice <- length(price)
  if(!(length(strike) %in% c(1,Nprice)))stop(
    "ImpliedVol: strike must be a scalar or of same length as price")
  if(!(length(expiry) %in% c(1,Nprice)))stop(
    "ImpliedVol: expiry must be a scalar or of same length as price")
  f <- function(vol){
    BSOptionValue(spot,Si,Ei,vol,intRate,carryCost,Call) - Pi
  }
  impVol <- rep(NA,Nprice)
  Si <- strike[1]
  Ei <- expiry[1]
  for (i in 1:Nprice){
    if(length(strike) > 1) Si <- strike[i]
    if(length(expiry) > 1) Ei <- expiry[i]
    Pi <- price[i]
    inrange <- f(ImpliedVolLowerBound)*f(ImpliedVolUpperBound) < 0
    if(is.na(inrange)) inrange <- FALSE
    if(inrange){impVol[i] <- uniroot(f,lower=ImpliedVolLowerBound,
      upper=ImpliedVolUpperBound,tol=tol)$root}
  }
  return(impVol)
}
# tests
if(FALSE){
  ImpliedVol(5,5.1,.3,0.1741752)
  ImpliedVol(5,5.1,.3,0.274152,Call=FALSE)
  ImpliedVol(5,rep(5.1,3),.3,c(0.1741752,.17,-.1))
  ImpliedVol(5,c(5.1,5.3),c(.3,.4),c(0.1741752,.17))
  ImpliedVol(5,5.2,c(.3,.4),c(0.1741752,.17))
  ImpliedVol(5,50,.7,0.001)
}

# Find parameters of lognormal distibution to match given mean and sd
lnorm.param <- function(mean, sd){
  cv <- sd/mean
  varlog <- log1p(cv*cv)
  return(list(meanlog=log(mean)-.5*varlog, sdlog=sqrt(varlog)))
}

######################################################
#Compute distribution function of stable distribution by numerical integration
# Do not bother with alpha=1, since can set oneminusalpha to be, say, 1.e-20.
# Leave out special cases alpha=2 and where distribution function is 0 or 1.

# Tamper with calculation of integrand to give good precision
# Note that the functions which compute integrands are always called with
#  a vector of length 21 in R.
# This tampering is based on the values which will be used when integrating
#  over a range -1 to 1.  These were obtained using the commands:
#    integrand <- function(theta){print(theta, digits=16);theta}
#    integrate(integrand, -1, 1)

####### Start of function to compute tail probabilities ########
tailsGstable <- function(x, logabsx, alpha, oneminusalpha, twominusalpha,
     beta, betaplus1, betaminus1, parametrization, lower.tail=TRUE){

  if(length(x) != 1) stop("tailsGstable: x must be of length 1")
  if(length(logabsx) != 1) stop("tailsGstable: logabsx must be of length 1")

  if(!missing(alpha)){
    # Case where alpha is specified
    if(length(alpha) != 1) stop("tailsGstable: alpha must be of length 1")
    if(!missing(oneminusalpha) || !missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    oneminusalpha <- 1 - alpha
    twominusalpha <- 2 - alpha
  } else if(!missing(oneminusalpha)){
    # Case where oneminusalpha is specified
    if(length(oneminusalpha) != 1) stop(
      "tailsGstable: oneminusalpha must be of length 1")
    if(!missing(twominusalpha)) stop(
      "Only specify one of alpha, oneminusalpha, twominusalpha")
    alpha <- 1 - oneminusalpha
    twominusalpha <- 1 + oneminusalpha
  } else {
    if(missing(twominusalpha))stop(
      "Specify one of alpha, oneminusalpha, twominusalpha")
    # Case where twominusalpha is specified
    if(length(twominusalpha) != 1) stop(
      "tailsGstable: twominusalpha must be of length 1")
    alpha <- 2 - twominusalpha
    oneminusalpha <- twominusalpha - 1
  }

  if(!missing(beta)){
    # Case where beta is specified
    if(length(beta) != 1) stop("tailsGstable: beta must be of length 1")
    if(!missing(betaplus1) || !missing(betaminus1)) stop(
      "Only specify one of beta, betaplus1, betaminus1")
    betaplus1 <- beta + 1
    betaminus1 <- beta - 1
  } else if(!missing(betaplus1)){
    # Case where betaplus1 is specified
    if(length(betaplus1) != 1) stop(
      "tailsGstable: betaplus1 must be of length 1")
    if(!missing(betaminus1)) stop(
      "Only specify one of beta, betaplus1, betaminus1")
    beta <- betaplus1 - 1
    betaminus1 <- betaplus1 - 2
  } else {
    if(missing(betaminus1))stop("Specify one of beta, betaplus1, betaminus1")
    # Case where betaminus1 is specified
    if(length(betaminus1) != 1) stop(
      "tailsGstable: betaminus1 must be of length 1")
    beta <- betaminus1 + 1
    betaplus1 <- betaminus1 + 2
  }

  INT.PTS <- c(0, -0.9739065285171717, 0.9739065285171717, -0.8650633666889845, 
    0.8650633666889845, -0.6794095682990244, 0.6794095682990244, 
    -0.4333953941292472, 0.4333953941292472, -0.1488743389816312, 
    0.1488743389816312, -0.9956571630258081, 0.9956571630258081, 
    -0.9301574913557082, 0.9301574913557082, -0.7808177265864169, 
    0.7808177265864169, -0.5627571346686047, 0.5627571346686047, 
    -0.2943928627014602, 0.2943928627014602)

  POWERSOFHALF <- c(1, cumprod(rep(.5, 60)))
  HPI <- pi/2
  LARGE <- 999

  integrand <- function(p){
    # p is the proportion of the integration range
    # Convention: global variables have names in capital letters
    # Test input data sometimes useful:  p <- INT.PTS*.5+.5
    # Find n such that integration range is 2^(-n) of full range
    n <- round(log(p[3]-p[2])/log(.5))
    if(n>53) stop("Integration interval too small")
    # Mid-point will have no rounding error, provided n < about 54
    # Check by printing:   1-POWERSOFHALF-1
    pmid <- p[1]
    qmid <- 1-pmid
    factor <- POWERSOFHALF[n+2]
    ps <- pmid+INT.PTS*factor
    qs <- qmid-INT.PTS*factor
  
    # Compute linear combinations of theta to high precision
    theta <- qs*THETA.LOW +ps*THETA.HIGH
    pctheta <- qs*PCTHETA.LOW +ps*PCTHETA.HIGH
    nctheta <- qs*NCTHETA.LOW +ps*NCTHETA.HIGH
    arg1 <- qs*ARG1.LOW +ps*ARG1.HIGH
    pcarg1 <- qs*PCARG1.LOW +ps*PCARG1.HIGH
    ncarg1 <- qs*NCARG1.LOW +ps*NCARG1.HIGH
    arg2 <- qs*ARG2.LOW +ps*ARG2.HIGH
    pcarg2 <- qs*PCARG2.LOW +ps*PCARG2.HIGH
    ncarg2 <- qs*NCARG2.LOW +ps*NCARG2.HIGH
  
    # Compute integrand to high precision
    # Integrand is exp(-w) where w=cos(theta-ALPHA*(theta-PHI0))*
    #    sin(ALPHA*(theta-PHI0))^(ALPHA/(1-ALPHA))*
    #    cos(theta)^(-1/(1-ALPHA))*X^(-ALPHA/(1-ALPHA)) 
    # Sometimes both sin(ALPHA*(theta-PHI0)) and X are negative. This is OK.
  
    # pctheta is theta+pi/2 
    # nctheta is theta-pi/2 
    # arg1 is alpha*(theta-phi0) 
    # pcarg1 is alpha*(theta-phi0) + pi/2
    # ncarg1 is alpha*(theta-phi0) - pi/2
    # arg2 is theta-alpha*(theta-phi0) 
    # pcarg2 is theta-alpha*(theta-phi0)+pi/2 
    # ncarg2 is theta-alpha*(theta-phi0)-pi/2 
  
    term1 <- l2 <- l3 <- rep(0, length(p))
    # Compute cos(arg2) as sine of one of
    #  theta+ALPHA*(theta-PHI0)+HPI and theta+ALPHA*(theta-PHI0)-HPI
    first <- abs(pcarg2) < abs(ncarg2)
    term1[first] <- sin(pcarg2[first])
    term1[!first] <- sin(-ncarg2[!first])
    # When arg1 is not near to unity, (say abs(arg1) > 1), 
    #  use sin(arg1)=1-2*sin(a)^2 where a is 0.5*(PI-arg1) or 0.5(arg1+PI)
    first <- abs(pcarg1) < abs(ncarg1)
    a <- ncarg1
    a[first] <- pcarg1[first]
    first <- abs(arg1) > 1
    l2[first] <- log1p(-2*sin(.5*a[first])^2)
    l2[!first] <- log(abs(sin(arg1[!first])))
    # When abs(theta) < 0.7, compute cos(theta) as 1-2*sin(0.5*theta)^5
    #  otherwise use sin(b) where b is theta+HPI or theta-HPI
    first <- abs(theta) < 0.7
    secondtest <- abs(pctheta) < abs(nctheta)
    second <- !first & secondtest
    third <- !first & !secondtest
    l3[first] <- log1p(-2*sin(.5*theta[first])^2)
    l3[second] <- log(sin(pctheta[second]))
    l3[third] <- log(sin(-nctheta[third]))
    w <- term1*exp((l3+ALPHA*(LOGABSX-l2))/ALPHAMINUS1)
    # Note that if w is infinite, we still want exp(-w) to give zero.
    # The next line is a safety precaution, knowing that exp(-LARGE) will give 0
    w[!is.finite(w)] <- LARGE
    # If subinterval is at an extreme of the integration region, 
    #  then ensure that the value for W nearest to the limit of the range
    #  of integration is within 1 of the value at the limit.
    # This is meant to stop the numerical procedure from completely missing
    #  the region where the integrand > 0, or is < 1.
    if(pmid==factor) w[2] <- if(LEFTLIMIT < RIGHTLIMIT){
      min(w[2], LEFTLIMIT+1) } else { max(w[2], LEFTLIMIT-1)}
    if(qmid==factor) w[3] <- if(LEFTLIMIT > RIGHTLIMIT){
      min(w[3], RIGHTLIMIT+1) } else { max(w[3], RIGHTLIMIT-1)}
    # Use complementary function if this will give more precise tail prob.
    if(COMPLEMENT){-expm1(-w)} else {exp(-w)}
  }

  # If x is outside the allowed exponents of real numbers, use +1 or -1 being
  #  careful to give correct sign.  Always input logabsx to full precision.
#  if(alpha > 1.5){alpha <- 2-twominusalpha; oneminusalpha <- twominusalpha-1} else
#   {if(alpha < .5){twominusalpha <- 2-alpha;oneminusalpha <- 1-alpha} else
#      {twominusalpha <- oneminusalpha+1; alpha <- 1-oneminusalpha}}
  k <- if(oneminusalpha>0){alpha}else{twominusalpha}
  # Change parametrization, if necessary
  if(parametrization < 2){
    a <- HPI*k
    # Compute tan(a) to high precision
    if(a<.8){tana <- tan(a)} else {tana <- 1/tan(HPI*abs(oneminusalpha))}
    btn <- beta*tana
    if(beta > 0.8){
      betaminus1 <- atan(betaminus1*tana/(1+btn*tana))/a
      beta <- betaminus1+1
      betaplus1 <- betaminus1+2
    } else { if(beta < -0.8){
      betaplus1 <- atan(betaplus1*tana/(1-btn*tana))/a
      beta <- betaplus1-1
      betaminus1 <- betaplus1-2
    } else {
      beta <- atan(btn)/a
      betaplus1 <- beta+1
      betaminus1 <- beta-1
      }
    }
    # If alpha>1 then change sign of beta
    if(oneminusalpha < 0 && parametrization < 2){
      beta <- -beta
      btn <- -btn
      temp <- betaminus1
      betaminus1 <- -betaplus1
      betaplus1 <- -temp
      }
    scale <- exp(-.5/alpha*log1p(btn^2)) # =(1+btn**2)**(-.5/alpha)
  }

  # Zolotarev's M = Nolan S0 (which is not used near alpha=0 or for extreme x)
  if(parametrization == 0){
    if(abs(btn) > 10*abs(x)){
      logabsx <- log1p(x/btn)-.5/alpha * log1p(1/(btn*btn))-
        oneminusalpha/alpha * log(abs(btn))
      x <- sign(btn)*exp(logabsx)
    } else {
      x <- (x+btn)*scale
      logabsx <- log(abs(x))
    }
  }

  # Zolotarev's A = Nolan S1
  if(parametrization == 1){
    x <- x*scale
    logabsx <- logabsx +log(scale)
  }

  # Zolotarev's C = Chambers, Mallows and Stuck
  if(parametrization == 2){
    scale <- 1.
  }

  # Code is not designed to work if x=0; alpha >= 2; alpha <=0; alpha=1;
  #  alpha<1, beta=1 & x<0; or alpha<1, beta=-1 & x>0
  if( x==0 || twominusalpha <= 0 || alpha <= 0 || oneminusalpha ==0 ||
    oneminusalpha > 0 && ((betaminus1==0 && x<0) || (betaplus1==0 && x>0))){
    print("Check parameters input to function tailsGstable")}
  phi0 <- -HPI*beta*k/alpha
  possibleWlimit <- abs(oneminusalpha)*abs(x/alpha)^(-alpha/oneminusalpha)
  if(oneminusalpha>0){  #Case when alpha<1
    if(x < 0){
  # THETA.LOW is the lower integration limit
  # PCTHETA.LOW is theta+pi/2 at the lower integration limit
  # NCTHETA.LOW is theta-pi/2 at the lower integration limit
  # ARG1.LOW is alpha*(theta-phi0) at the lower integration limit
  # PCARG1.LOW is alpha*(theta-phi0) + pi/2 at the lower integration limit
  # NCARG1.LOW is alpha*(theta-phi0) - pi/2 at the lower integration limit
  # ARG2.LOW is theta-alpha*(theta-phi0) at the lower integration limit
  # PCARG2.LOW is theta-alpha*(theta-phi0)+pi/2 at the lower integration limit
  # NCARG2.LOW is theta-alpha*(theta-phi0)-pi/2 at the lower integration limit

      THETA.LOW <- -HPI
      ARG1.LOW <- HPI*alpha*betaminus1
      THETA.HIGH <- phi0
      ARG1.HIGH <- 0
      LEFTLIMIT <- LARGE
      RIGHTLIMIT <- if(betaplus1==0){possibleWlimit}else{0}
    } else {
      THETA.LOW <- phi0
      ARG1.LOW <- 0
      THETA.HIGH <- HPI
      ARG1.HIGH <- HPI*alpha*betaplus1
      LEFTLIMIT <- if(betaminus1==0){possibleWlimit}else{0}
      RIGHTLIMIT <- LARGE
    }
    # Use complementary integrand whenever right tail
    # required or integrating over less than -HPI to HPI
    singleregion <- (phi0+HPI) == 0 || (phi0-HPI) == 0
    COMPLEMENT <- !(lower.tail && singleregion)
  } else {    #Case when alpha>1
    if(x < 0){
      THETA.LOW <- -HPI
      ARG1.LOW <- HPI*(-alpha+beta*twominusalpha)
      THETA.HIGH <- phi0
      ARG1.HIGH <- 0
      LEFTLIMIT <- if(betaplus1>0){0}else{possibleWlimit}
      RIGHTLIMIT <- LARGE
    } else {
      THETA.LOW <- phi0
      ARG1.LOW <- 0
       THETA.HIGH <-  HPI
      ARG1.HIGH <- HPI*(alpha+beta*twominusalpha)
      LEFTLIMIT <- LARGE
      RIGHTLIMIT <- if(betaminus1<0){0}else{possibleWlimit}
    }
    COMPLEMENT <- FALSE
  }

  PCTHETA.LOW <- THETA.LOW+HPI  
  NCTHETA.LOW <- THETA.LOW-HPI
  PCTHETA.HIGH <- THETA.HIGH+HPI  
  NCTHETA.HIGH <- THETA.HIGH-HPI
  PCARG1.LOW <- ARG1.LOW+HPI
  NCARG1.LOW <- ARG1.LOW-HPI
  PCARG1.HIGH <- ARG1.HIGH+HPI
  NCARG1.HIGH <- ARG1.HIGH-HPI
  ARG2.LOW <- THETA.LOW-ARG1.LOW
  PCARG2.LOW <- PCTHETA.LOW-ARG1.LOW
  NCARG2.LOW <- NCTHETA.LOW-ARG1.LOW
  ARG2.HIGH <- THETA.HIGH-ARG1.HIGH
  PCARG2.HIGH <- PCTHETA.HIGH-ARG1.HIGH
  NCARG2.HIGH <- NCTHETA.HIGH-ARG1.HIGH

  LOGABSX <- logabsx
  ALPHA <- alpha
  ALPHAMINUS1 <- -oneminusalpha

  result <- integrate(integrand, 0, 1, subdivisions=1000, stop.on.error=F, 
    abs.tol=1.e-50, rel.tol=1.e-14)

  contribution <- result$value*(THETA.HIGH-THETA.LOW)/pi

  if(x<0){
    lefttail <- contribution; righttail <- 1-lefttail
  }else{
    righttail <- contribution; lefttail <- 1-righttail}

  # If alpha<1 & lower.tail & singleregion then swap tails around
  swap <- if(oneminusalpha < 0) {FALSE} else {lower.tail && singleregion}
  if(swap){temp <- lefttail; lefttail <- righttail; righttail <- temp}
  list(left.tail.prob=lefttail, right.tail.prob=righttail, 
    est.error=result$abs.error, message=result$message)
}

# Compute quantiles of standard stable distributions
# Prototype code to be converted into C later.
# Logarithms of both tail probabilities assumed to be provided.
# Return both z and logz, with logz being meaningless if z <= 0
# Deals with only one quantile at a time.
# The stable distribution object passed must have standard location and scale
standardStableQuantile <- function(logp, logq, obj, browse=FALSE){
  if(length(logp) > 1 || length(logq) > 1) stop(
    "standardStableQuantile: logp and logq must both be scalars")
  if(abs(1-exp(logp)-exp(logq)) > 1.e-14) stop(
    "standardStableQuantile: logp and logq must be consistent")
  lower.tail <- logp < logq
  if(lower.tail){ if(logp == -Inf) stop(
    "standardStableQuantile: logp input as -Inf")
  } else { if(logq == -Inf) stop(
    "standardStableQuantile: logq input as -Inf")
  }

  if(browse)browser()

  # If alpha==2, return quantile immediately
  if(obj$twominusalpha ==0){
    z <- sqrt(2)* if(lower.tail) qnorm(logp, log.p=TRUE) else qnorm(logq,
      log.p=TRUE, lower.tail=FALSE)
    return(z)
  }

  # First stage: Find reasonable first approximation
  if(obj$alpha < .5){
    uselogz <- TRUE
    if(lower.tail){
      # First approximation to xi
      xi <- -logp
      # One Newton-Raphson iteration
      xi <- xi -(logp + .5*log(2*pi*obj$alpha*xi) +xi)/(.5/xi + 1)
      # Find logz for this xi
      logz <- log(obj$alpha)-obj$oneminusalpha/obj$alpha*log(xi/obj$oneminusalpha)
    } else {
      logCalpha <- log(gamma(obj$alpha)*sin(pi*obj$alpha)/pi)
      logz <- (logq-logCalpha) *(-1/obj$alpha)
    }  
  } else if (obj$alpha > 1.7 && min(logp, logq) > -20){
    # Case where quantile of normal distribution gives a useful starting point

# Speed would be improved by giving a better approximate root for
#    moderate alpha, moderate p

    uselogz <- FALSE
    z <- sqrt(2)* qnorm(logp, log.p=TRUE)
  } else if (!lower.tail){
    uselogz <- TRUE
    # Case of right tail for alpha >= .5
    logy <- (log(2*gamma(obj$alpha + 1)* sin(.5 * pi *obj$twominusalpha)/
      (pi*obj$alpha)) - logq)/obj$alpha
    if(obj$oneminusalpha ==0 || logy > 100) logz <- logy else{
      logz <- logy + log1p(expm1(obj$oneminusalpha * logy)/
       (tan(pi*obj$oneminusalpha/2) * exp(logy)))
    }
  } else {
    # Case of left tail for alpha >= 0.5
    uselogz <- FALSE
    # First approximation to xi
    xi <- -logp
    # One Newton-Raphson iteration
    xi <- xi -(logp + .5*log(2*pi*obj$alpha*xi) +xi)/(.5/xi + 1)
    # Find z for this xi
    if(obj$oneminusalpha == 0){ # i.e. alpha = 1
      z <- -(1+log(.5*pi*xi))/(.5*pi)
    } else if(obj$alpha < 1.9){ # i.e. not near alpha==2
      ang <- .5*pi*obj$oneminusalpha
      cosa <- if(obj$alpha < 1.5) cos(ang) else sin(.5*pi*obj$twominusalpha)
      z <- expm1(  log(obj$alpha/cosa)+
        (obj$oneminusalpha/obj$alpha)*log(obj$oneminusalpha/
        (xi* sin(ang)))  )/tan(ang)
    } else  {	#i.e. when alpha is near 2
      ang <- .5*pi*obj$twominusalpha
      z <- -expm1(  log(obj$alpha/sin(ang))+
        (obj$oneminusalpha/obj$alpha)*log(-obj$oneminusalpha/
        (xi* cos(ang)))  )*tan(ang)
    }
  }

  # Temporary intermediate stage:  Check that approximation is OK
  if(browse){
    if(uselogz){
      if(logz > 700) return(Inf)
      z <- exp(logz)
    }
    if(!is.finite(z)) return(Inf)
    frac.error <- if(lower.tail){
      pEstable(z, obj, lower.tail=TRUE)/exp(logp) - 1
    } else {
      pEstable(z, obj, lower.tail=FALSE)/exp(logq) - 1
    }
    print(paste("Fractional error:", frac.error), quote=FALSE)
  }

  # Second stage: Refine approximation by Newton-Raphson
  if(uselogz){
    if(lower.tail){
      adjustmentlog <- function(logz){
        tls <- tailsEstable(exp(logz), obj)
        (logp - tls$logF)*exp(-logz + tls$logF - tls$logdensity)
      }
    } else {
      adjustmentlog <- function(logz){
        tls <- tailsEstable(exp(logz), obj)
        (tls$logrighttail - logq)*exp(-logz + tls$logrighttail - tls$logdensity)
      }
    }
    for (i in 1:10){
      adj <- adjustmentlog(logz)
      logz <- logz + adj
      if(logz > 700) return(Inf)
      if(abs(adj) < 1.e-8* max(1, abs(logz))) break
    }
    if(browse) print(paste("Log z scale:",i,"iterations used"), quote=FALSE)
    z <- exp(logz)
  } else {
    if(lower.tail){
      adjustment <- function(z){
        tls <- tailsEstable(z, obj)
        (logp - tls$logF)*exp(tls$logF - tls$logdensity)
      }
    } else {
      adjustment <- function(z){
        tls <- tailsEstable(z, obj)
        (tls$logrighttail - logq)*exp(tls$logrighttail - tls$logdensity)
      }
    }
    for (i in 1:10){
      adj <- adjustment(z)
      z <- z + adj
      if(abs(adj) < 1.e-8 * max(1, abs(z))) break
      if(!is.finite(z) || abs(z) > 1.e300) return(NA)
    }
    if(browse) print(paste("Using z scale:",i,"iterations used"), quote=FALSE)
    logz <- if(z>0) log(z) else NA
  }
  return(z)
}

# Quantile-finding using R function standardStableQuantile
# Find quantiles of a stable distribution
qEstable <- function(p, stableParamObj, log=FALSE, lower.tail=TRUE){
  if(class(stableParamObj) != "stableParameters")stop(paste("qEstable:",
    "Parameter stableParamObj must be of class stableParameters"))
  if(!all(is.finite(p))) stop("qEstable: Components of p must all be finite")
  if(log){
    logp <- p
    if(any(logp >= 0)) stop("qEstable: Logs of probabilities must be < 0")
    complement <- log1p(-exp(p))
  } else {
    if(any(p >= 1)) stop("qEstable: Require p < 1")
    logp <- log(p)
    complement <- log1p(-p)
  }
  scale <- exp(stableParamObj$logscale)
  results <- double(length(p))
  standard <- stableParamObj
  standard$location <- 0
  standard$logscale <- 0
  standard$created.by <- "function qEstable"
  for (i in 1:length(p)){
    if(lower.tail){
      r <- standardStableQuantile(logp[i], complement[i], standard)
    } else {
      r <- standardStableQuantile(complement[i], logp[i], standard)
    }
    results[i] <- r * scale + stableParamObj$location
  }
  return(results)
}
# Find quantiles of a log stable distribution
qFMstable <- function(p, stableParamObj, lower.tail=TRUE){
  if(!all(is.finite(p))) stop("qFMstable: Components of p must all be finite")
  if(class(stableParamObj) != "stableParameters")stop(paste("qFMstable:",
    "Parameter stableParamObj must be of class stableParameters"))
  return(exp(-qEstable(p, stableParamObj, lower.tail=!lower.tail, log=FALSE)))
}

Try the FMStable package in your browser

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

FMStable documentation built on May 29, 2017, 7:20 p.m.