R/functions.R

In FMStable: Finite Moment Stable Distributions

# 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(!inherits(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(!inherits(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(!inherits(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(!inherits(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(!inherits(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(!inherits(stableParamObj, "stableParameters")) stop(paste("qFMstable:",
    "Parameter stableParamObj must be of class stableParameters"))
  return(exp(-qEstable(p, stableParamObj, lower.tail=!lower.tail, log=FALSE)))
}