Calibrating an Implied Volatility Model with Nested Optimization
Summary
The document explains how to fit a five-parameter implied volatility model when its pricing relation gives strike as a function of call price, rather than directly giving call price for a strike. The proposed structure uses an inner numerical search to recover each model-implied call price, then an outer optimization to minimize squared differences from observed prices. It illustrates this with bounded scalar minimization for the inner step and a brute-force parameter search for the outer step.
The example is a rough implementation guide, not a validated calibration procedure. The answer cautions that global optimization over five parameters is difficult and may produce unreliable results, and expresses doubt about practical usability. It does not report empirical calibration results or demonstrate that the fitted parameters produce a reliable risk-neutral density. The search bounds and numerical precision in the sample are implementation choices that may need adjustment for actual data and model behavior.
Key ideas
- The inverse pricing relation requires numerical recovery of each model-implied call price.
- An inner scalar optimization can estimate call price for each observed strike.
- An outer optimization can then minimize squared errors across observed prices.
- The example uses bounded scalar minimization inside a brute-force parameter search.
- The answer warns that five-parameter global optimization may yield unreliable fits.
Tags
Full text
# Fitting parameters given an inverse function. (Orosi, 2015)
# Fitting parameters given an inverse function. (Orosi, 2015)
In trying to replicate Orosi's (2015) 5-parameter implied volatility model, but I can't wrap my head around the parameter fitting procedure Orosi proposes. My main goal is to calibrate the model to my set of observed call prices, so that I can derive the Risk-neutral-density function.
Consider that I have empirical observations of call prices, $c$ for a range of exercise prices, $x$, for a given time and maturity. I then have to minimise the squared errors:
$\sum_{i=1}^{n}(c(x_i)-c_i)^2$
Where $c(x_i)$ has to be numerically determined from the equation:
$x_i=1-c(x_i)+G\frac{(1-c(x_i))^\beta}{c(x_i)^\alpha}*(Ac(x_i)+1)^\gamma$
With the following constraints:
$G>0,\alpha>0,\beta>2,\gamma>1,0>A>-1$
Now I am pretty comfortable with minimising the objective function in Python using e.g. scipy's fmin function etc. However, as noted by the author, as $c(x)$ is provided in an inverted manner "$c(x)$ has to be numerically derived using a golden section search or something similar". I have no idea how to go about this and to implement it in a "usual" non-linear least squares minimisation method?
If anyone could point me in the right direction or provide Python or pseudo-code so I can figure it out myself it would be highly appreciated
## Answer by Trevor Hansen (score 3, accepted)
https://quant.stackexchange.com/a/65659
You need to use some form of optimization technique twice. Once to solve for $c(x_i)$ and another one to minimize the MSE. I have added below a quick hack of the problem where minimize_scalar is used to solve for $c(x_i)$ and the brute force algorithm is used to solve for the parameters that minimize the MSE. There are definitely better optimizations to be used, but this should provide the structure you are looking for.
Words of caution: as mentioned in the comments I am skeptical if this is at all usable in practice. Also a 5 variable global minimization problem is very difficult to solve and I suspect the solutions will be unreliable and unusable in practice.
Good luck!
```
from scipy import optimize
PRECISION = 1e-3
def error_c(c, params, x):
error = abs(1-c+params[0]*((1-c)**params[2]/(c**params[1]))*(params[4]*c+1)**params[3]-x)
return error
def get_c(x, params):
c = optimize.minimize_scalar(error_c,
bounds=(PRECISION, 1/PRECISION),
method='bounded',
args=(params, x),
tol=PRECISION)
return c.x
def mse(params, x, c):
mse = 0
for x_i, c_i in zip(x, c):
mse += (get_c(x_i, params)-c_i)**2
return(mse)
x = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
c = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
bounds = ((PRECISION, 1/PRECISION),
(PRECISION, 1/PRECISION),
(2+PRECISION, 1/PRECISION),
(1+PRECISION, 1/PRECISION),
(-1+PRECISION, -PRECISION))
result = optimize.brute(mse, ranges=bounds, args=(x, c), Ns=10)
print(result)
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.