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)))

}

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.