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

Don't call a variable `T`. That's an alias for `TRUE` in R and overwriting that could have unforeseen consequences. However, you need to check if all variables used for calculating `d1` and `d2` are numeric. — Roland, Sep 23 '13 at 12:13
Great, .... except that doesn't help people who are trying to help you. — IRTFM, Sep 23 '13 at 12:46
As an example: sigmaS <- read.csv2("Equity Volatility BoA.csv"). However, I think this is not the problem. I struggle more with the formula: does it make sense like this? Is the start-command correct? Are I am right, that these are the starting values for the optimization process? — schloni, Sep 23 '13 at 12:46
You do not seem to assign anything to `y`: your function is just returning `c(0,0)` each time. — Vincent Zoonekynd, Sep 23 '13 at 12:48
Don't depend on variables in the parent environment showing up properly inside your function. It's always a much safer route to require **all** variables to be arguments to your function (D1,R,T, SO1 in this case). Then do some sanity checks such as graphing your function to make sure it's got the shape you expect and that all variables are of appropriate class. Only then are you ready to start working with nonlinear fitting tools. — Carl Witthoft, Sep 23 '13 at 12:50
I am the author of nleqslv and I'm quite suprised at how you are using it. — Bhas, Sep 23 '13 at 17:36
My previous comment got lost because of the time limit and my use of the ENTER key. See my answer. — Bhas, Sep 23 '13 at 17:45
You can also consider other methods such as maximum likelihood. This is implemented in [this package](https://github.com/boennecd/DtD). It tends to more stable but requires that you have more than two periods of observations. — Benjamin Christoffersen, Apr 29 '18 at 12:53

score 11 · Answer 1 · edited Sep 13 '17 at 14:33

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.

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

1 Answers1

Linked