Solve systems of non-linear equations in R / Black-Scholes-Merton Model

Question

I am writing my masters thesis and I got stuck with this problem in my R code. Mathematically, I want to solve this system of non-linear equations with the R-package “nleqslv”:

fnewton <- function(x){

y <- numeric(2)

d1 = (log(x[1]/D1)+(R+x[2]^2/2)*T)/x[2]*sqrt(T)

d2 = d1-x[2]*sqrt(T)

y1 <- SO1 - (x[1]*pnorm(d1) - exp(-R*T)*D1*pnorm(d2))

y2 <- sigmaS*SO1 - pnorm(d1)*x[2]*x[1]

y}

xstart <- c(21623379, 0.526177094846878)

nleqslv(xstart, fnewton, control=list(btol=.01), method="Newton")

I have tried several versions of this code and right now I get the error:

error: error in pnorm(q, mean, sd, lower.tail, log.p): not numerical.

Pnorm is meant to be the cumulative standard Normal distribution of d1and d2 respectively. I really don’t know, what I am doing wrong as I oriented my model on Teterevas slides ( on slide no.5 is her model code), who’s presentation is the first result by googeling

https://www.google.de/search?q=moodys+KMV+in+R&rlz=1C1SVED_enDE401DE401&aq=f&oq=moodys+KMV+in+R&aqs=chrome.0.57.13309j0&sourceid=chrome&ie=UTF-8#q=distance+to+default+in+R

Like me, however more successfull, she calculates the Distance to Default risk measure via the Black-Scholes-Merton approach. In this model, the value of equity (usually represented by the market capitalization, ->SO1) can be written as a European call option – what I labeled y2 in the above code, however, the equation before is set to 0!

The other variables are:

x[1] -> the variable I want to derive, value of total assets

x[2] -> the variable I want to derive, volatility of total assets

D1 -> the book value of debt (1998-2009)

R -> a risk-free interest rate

T -> is set to 1 (time)

sigmaS -> estimated (historical) equity volatility

Thanks already!!! I would be glad, anyone could help me. Caro

Answer 1

I am the author of nleqslv and I'm quite suprised at how you are using it. As mentioned by others you are not returning anything sensible.

y1 should be y[1] and y2 should be y[2]. If you want us to say sensible things you will have to provide numerical values for D1, R, T, sigmaS and SO1. I have tried this:

T <- 1; D1 <- 1000; R <- 0.01; sigmaS <- .1; SO1 <- 1000  

These have been entered before the function definition. See this

library(nleqslv)

T <- 1
D1 <- 1000
R <- 0.01

sigmaS <- .1
SO1 <- 1000

fnewton <- function(x){
    y <- numeric(2)
    d1 <- (log(x[1]/D1)+(R+x[2]^2/2)*T)/x[2]*sqrt(T)
    d2 <- d1-x[2]*sqrt(T)
    y[1] <- SO1 - (x[1]*pnorm(d1) - exp(-R*T)*D1*pnorm(d2))
    y[2] <- sigmaS*SO1 - pnorm(d1)*x[2]*x[1]
    y
}

xstart <- c(21623379, 0.526177094846878)

nleqslv has no problem in finding a solution in this case. Solution found is : c(1990.04983,0.05025). There appears to be no need to set the btol parameter and you can use method Broyden.