[R] optim function multi variables -newbie

[Hans W Borchers]

michalseneca <michalseneca <at> gmail.com> writes:

> I tried to modify the code,thanks for noticing...
> now i get that the function cannot be evaluated at initial parameters.
> However I do not know what does that mean..
> Should I try to modify parameters. I am still not sure of syntax of
> function. I cannot get that optimize work correctly.
> I am trying to test that..however I cannot get to a result..
> In matlab however that works correctly...
> Did anybody of you tried to succesfully run it ?
>
> Michal


I went through your functions and the Matlab code, below I will append my
versions. Matlab and R now return the same values, at least for a few test
cases. One obstacle, for instance, was that solve.polynom() expects the
coefficients in reverse order.

Here is an example of how to apply optim(). Whether the result is
resonable you'll have to find out for yourself:

    MktVol <- c(20.8, 19.8, 18.1, 16.1, 15.1, 14.5, 14.0, 14.3, 15.0,
                15.9, 16.2, 16.4, 16.6, 17.3, 19.1) / 100;
    MktStrike <- c(4.3, 4.8, 5.3, 5.8, 6.3, 6.8, 7.3, 7.8, 8.3, 8.8,
                   9.3, 9.8, 10.3, 10.8, 11.3) / 100;
    F <- 0.073
    ATMVol <- 0.14
    T <- 1
    b<-1
    parsbox<-c(0.7,0.5)

    f <- function(x) EstimateRhoAndVol(x, MktStrike, MktVol, ATMVol, F, T, b)
    optim(parsbox, f)
    # $par
    # [1] -0.06374008  0.60383201

Of course, there are some more open points, for instance will 'findAlpha'
always find a real, positive cubic root? etc.

--Hans Werner

### --- snip ---------------------------

library(polynom)

EstimateRhoAndVol <- function(params, MktStrike, MktVol, ATMVol, F, T, b)
{
# -------------------------------------------------------------------------
# Returns the following SABR parameters:
# r = rho
# v = vol-of-vol
# Uses ATM volatility to estimate alpha
# Required inputs:
# MktStrike = Vector of Strikes
# MktVol    = Vector of corresponding volatilities
# ATMVol = ATM volatility
# F = spot price
# T = maturity
# b = beta parameter
# -------------------------------------------------------------------------
  r <- params[1]
  v <- params[2]
  a <- findAlpha(F, F, T, ATMVol, b, r, v)
  N <- length(MktVol)

  ModelVol <- numeric(N)
  error <- numeric(N)
  # Define the model volatility and the squared error terms
  for (i in 1:N) {
    ModelVol[i] = SABRvol(a, b, r, v, F, MktStrike[i], T)
    error[i] = (ModelVol[i] - MktVol[i])^2
  }

  # Return the SSE
  y <- sum(error, na.rm=TRUE)

  # Impose the constraint that -1 <= rho <= +1
  # via a penalty on the objective function
  if (abs(r) > 1)
    y <- Inf
  return(y)
}


SABRvol <- function(a, b, r, v, F, K, T)
{
# ---------------------------------------------
# Returns the SABR volatility.
# Required inputs:
# a = alpha parameter
# b = beta parameter
# r = rho parameter
# v = vol of vol parameter
# F = spot price
# K = strike price
# T = maturity
# ---------------------------------------------

  # Separate into ATM case (simplest case) and
  # Non-ATM case (most general case)

  if (abs(F-K) <= 0.001) {  # ATM vol
	Term1 <- a/F^(1-b)
	Term2 <- ((1-b)^2/24*a^2/F^(2-2*b) + r*b*a*v/4/F^(1-b) + (2-3*r^2)*v^2/24)
	y <- Term1 * (1 + Term2*T)

  } else {  # Non-ATM vol
	FK <- F*K
	z <- v/a*(FK)^((1-b)/2)*log(F/K)
	x <- log((sqrt(1 - 2*r*z + z^2) + z - r)/(1-r))
	Term1 <- a / FK^((1-b)/2) / (1 + (1-b)^2/24*log(F/K)^2 + 
	         (1-b)^4/1920*log(F/K)^4)
	if (abs(x-z) < 1e-10) Term2 <- 1
	else                  Term2 <- z / x
	Term3 <- 1 + ((1-b)^2/24*a^2/FK^(1-b) + 
	             r*b*v*a/4/FK^((1-b)/2) + (2-3*r^2)/24*v^2)*T
	y <- Term1 * Term2 * Term3
  }
  return(y)
}


findAlpha <- function(F, K, T, ATMvol, b, r, v)
{
# -------------------------------------------------------------------------
# Finds the SABR "alpha" parameter by solving the cubic equation described
# in West (2005) "Calibration of the SABR Model in Illiquid Markets"
# The function can return multiple roots.  In that case, the program 
# eliminates negative roots, and select the smallest root among the 
# positive roots that remain.
# Required inputs:
# F = spot price
# K = strike price
# T = maturity
# ATMvol = ATM market volatility
# b = beta parameter
# r = rho parameter
# v = vol of vol parameter
# -------------------------------------------------------------------------

  # Find the coefficients of the cubic equation for alpha
  C0 <- -ATMvol * F^(1-b)
  C1 <- (1 + (2-3*r^2) * v^2 * T / 24)
  C2 <- r * b * v * T / 4 / F^(1-b)
  C3 <- (1-b)^2 * T / 24 / F^(2-2*b)

  # Return the roots of the cubic equation (multiple roots)
  AlphaVector <- solve(as.polynomial(c(C0, C1, C2, C3)))

  # Find and return the smallest positive root
  index <- which(Im(AlphaVector) == 0 & Re(AlphaVector) > 0)
  Alpha <- AlphaVector[index]
  if (length(Alpha) == 0) Alpha <- 0
  return(min(Re(Alpha)))
}