Securities

factors_r = ["SP500", "DTWEXAFEGS"] # "SP500" does not contain dividends; note: "DTWEXM" discontinued as of Jan 2020
factors_d = ["DGS10", "BAMLH0A0HYM2"]

Black-Scholes model

from scipy.stats import norm

def level_shock(shock, S, tau, sigma):
    
    result = S * (1 + shock * sigma * np.sqrt(tau))
    
    return result

• https://en.wikipedia.org/wiki/Greeks_(finance)
• https://www.wolframalpha.com/input/?i=option+pricing+formula

factor = "SP500"
types = ["call", "put"]
S = levels_df.ffill()[factor].iloc[-1]
K = S
r = 0 # use "USD3MTD156N"
q = 0 # see https://stackoverflow.com/a/11286679 
tau = 1 # = 252 / 252
sigma = sd_df[factor].iloc[-1] # use "VIXCLS"
shocks = [x / 2 for x in range(-6, 7)]

greeks_df = pd.DataFrame([(x, y) for x in types for y in shocks], 
    columns = ["type", "shock"])
greeks_df["spot"] = level_shock(greeks_df["shock"], S, tau, sigma)

Value

For a given spot price \(S\), strike price \(K\), risk-free rate \(r\), annual dividend yield \(q\), time-to-maturity \(\tau = T - t\), and volatility \(\sigma\):

\[ \begin{aligned} V_{c}&=Se^{-q\tau}\Phi(d_{1})-e^{-r\tau}K\Phi(d_{2}) \\ V_{p}&=e^{-r\tau}K\Phi(-d_{2})-Se^{-q\tau}\Phi(-d_{1}) \end{aligned} \]

def bs_value(type, S, K, r, q, tau, sigma, d1, d2):
    
    if (type == "call"):
        result =  S * np.exp(-q * tau) * Phi(d1) - np.exp(-r * tau) * K * Phi(d2)
    elif (type == "put"):
        result = np.exp(-r * tau) * K * Phi(-d2) - S * np.exp(-q * tau) * Phi(-d1)
        
    return result

def bs_value(type, S, K, r, q, tau, sigma, d1, d2):
    
    r_df = np.exp(-r * tau)
    q_df = np.exp(-q * tau)
    
    call_value = S * q_df * Phi(d1) - r_df * K * Phi(d2)
    put_value = r_df * K * Phi(-d2) - S * q_df * Phi(-d1)
    result = np.where(type == "call", call_value, put_value)
    
    return result

where

\[ \begin{aligned} d_{1}&={\frac{\ln(S/K)+(r-q+\sigma^{2}/2)\tau}{\sigma{\sqrt{\tau}}}} \\ d_{2}&={\frac{\ln(S/K)+(r-q-\sigma^{2}/2)\tau}{\sigma{\sqrt{\tau}}}}=d_{1}-\sigma{\sqrt{\tau}} \\ \phi(x)&={\frac{e^{-{\frac {x^{2}}{2}}}}{\sqrt{2\pi}}} \\ \Phi(x)&={\frac{1}{\sqrt{2\pi}}}\int_{-\infty}^{x}e^{-{\frac{y^{2}}{2}}}dy=1-{\frac{1}{\sqrt{2\pi}}}\int_{x}^{\infty}e^{-{\frac{y^{2}}{2}}dy} \end{aligned} \]

def bs_d1(S, K, r, q, tau, sigma):
    
    result = (np.log(S / K) + (r - q + sigma ** 2 / 2) * tau) / (sigma * np.sqrt(tau))
    
    return result

def bs_d2(S, K, r, q, tau, sigma):
    
    result = (np.log(S / K) + (r - q - sigma ** 2 / 2) * tau) / (sigma * np.sqrt(tau))
    
    return result
    
def phi(x):
    
    result = norm.pdf(x)
    
    return result

def Phi(x):
    
    result = norm.cdf(x)
    
    return result

greeks_df["d1"] = bs_d1(greeks_df["spot"], K, r, q, tau, sigma)
greeks_df["d2"] = bs_d2(greeks_df["spot"], K, r, q, tau, sigma)
greeks_df["value"] = bs_value(greeks_df["type"], greeks_df["spot"], K, r, q, tau, sigma,
                              greeks_df["d1"], greeks_df["d2"]) 

First-order

Delta

\[ \begin{aligned} \Delta_{c}&={\frac{\partial V_{c}}{\partial S}}=e^{-q\tau}\Phi(d_{1}) \\ \Delta_{p}&={\frac{\partial V_{p}}{\partial S}}=-e^{-q\tau}\Phi(-d_{1}) \end{aligned} \]

def bs_delta(type, S, K, r, q, tau, sigma, d1, d2):
  
    q_df = np.exp(-q * tau)
    
    call_value = q_df * Phi(d1)
    put_value = -q_df * Phi(-d1)
    result = np.where(type == "call", call_value, put_value)
    
    return result

greeks_df["delta"] = bs_delta(greeks_df["type"], greeks_df["spot"], K, r, q, tau, sigma,
                              greeks_df["d1"], greeks_df["d2"]) 

Delta-beta

Notional market value is the market value of a leveraged position:

\[ \begin{aligned} \text{Equity options }=&\,\#\text{ contracts}\times\text{multiple}\times\text{spot price}\\ \text{Delta-adjusted }=&\,\#\text{ contracts}\times\text{multiple}\times\text{spot price}\times\text{delta} \end{aligned} \]

• https://en.wikipedia.org/wiki/Notional_amount

def bs_delta_diff(type, S, K, r, q, tau, sigma, delta0):
    
    d1 = bs_d1(S, K, r, q, tau, sigma)
    d2 = bs_d2(S, K, r, q, tau, sigma)
    delta = bs_delta(type, S, K, r, q, tau, sigma, d1, d2)
    
    call_value = delta - delta0
    put_value = delta0 - delta
    
    result = np.where(type == "call", call_value, put_value)
        
    return result

beta = 0.35
type = "call"
n = 1
multiple = 100
total = 1000000

d1 = bs_d1(S, K, r, q, tau, sigma)
d2 = bs_d2(S, K, r, q, tau, sigma)
sec = {
    "n": n,
    "multiple": multiple,
    "S": S,
    "delta": bs_delta(type, S, K, r, q, tau, sigma, d1, d2),
    "beta": 1
}

beta_df = pd.DataFrame([(x, y) for x in types for y in shocks], 
    columns = ["type", "shock"])
beta_df["spot"] = level_shock(beta_df["shock"], S, tau, sigma)
beta_df["static"] = beta
beta_df["diff"] = bs_delta_diff(type, beta_df["spot"], K, r, q, tau, sigma, sec["delta"])
beta_df["dynamic"] = beta + sec["n"] * sec["multiple"] * sec["S"] * sec["beta"] * beta_df["diff"] / total

For completeness, duration equivalent is defined as:

\[ \begin{aligned} \text{10-year equivalent }=\,&\frac{\text{security duration}}{\text{10-year OTR duration}} \end{aligned} \]

Vega

\[ \begin{aligned} \nu_{c,p}&={\frac{\partial V_{c,p}}{\partial\sigma}}=Se^{-q\tau}\phi(d_{1}){\sqrt{\tau}}=Ke^{-r\tau}\phi(d_{2}){\sqrt{\tau}} \end{aligned} \]

def bs_vega(type, S, K, r, q, tau, sigma, d1, d2):
    
    q_df = np.exp(-q * tau)
    result = S * q_df * phi(d1) * np.sqrt(tau)
    
    return result

greeks_df["vega"] = bs_vega(greeks_df["type"], greeks_df["spot"], K, r, q, tau, sigma,
                            greeks_df["d1"], greeks_df["d2"]) 

Theta

\[ \begin{aligned} \Theta_{c}&=-{\frac{\partial V_{c}}{\partial \tau}}=-e^{-q\tau}{\frac{S\phi(d_{1})\sigma}{2{\sqrt{\tau}}}}-rKe^{-r\tau}\Phi(d_{2})+qSe^{-q\tau}\Phi(d_{1}) \\ \Theta_{p}&=-{\frac{\partial V_{p}}{\partial \tau}}=-e^{-q\tau}{\frac{S\phi(d_{1})\sigma}{2{\sqrt{\tau}}}}+rKe^{-r\tau}\Phi(-d_{2})-qSe^{-q\tau}\Phi(-d_{1}) \end{aligned} \]

def bs_theta(type, S, K, r, q, tau, sigma, d1, d2):
    
    r_df = np.exp(-r * tau)
    q_df = np.exp(-q * tau)
    
    call_value = -q_df * S * phi(d1) * sigma / (2 * np.sqrt(tau)) - \
        r * K * r_df * Phi(d2) + q * S * q_df * Phi(d1)
        
    put_value = -q_df * S * phi(d1) * sigma / (2 * np.sqrt(tau)) + \
        r * K * r_df * Phi(-d2) - q * S * q_df * Phi(-d1)
        
    result = np.where(type == "call", call_value, put_value)
    
    return result

greeks_df["theta"] = bs_theta(greeks_df["type"], greeks_df["spot"], K, r, q, tau, sigma,
                              greeks_df["d1"], greeks_df["d2"]) 

Second-order

Gamma

\[ \begin{aligned} \Gamma_{c,p}&={\frac{\partial\Delta_{c,p}}{\partial S}}={\frac{\partial^{2}V_{c,p}}{\partial S^{2}}}=e^{-q\tau}{\frac{\phi(d_{1})}{S\sigma{\sqrt{\tau}}}}=Ke^{-r\tau}{\frac{\phi(d_{2})}{S^{2}\sigma{\sqrt{\tau}}}} \end{aligned} \]

def bs_gamma(type, S, K, r, q, tau, sigma, d1, d2):
  
    q_df = np.exp(-q * tau)
    
    result = q_df * phi(d1) / (S * sigma * np.sqrt(tau))
    
    return result

greeks_df["gamma"] = bs_gamma(greeks_df["type"], greeks_df["spot"], K, r, q, tau, sigma,
                              greeks_df["d1"], greeks_df["d2"]) 

Taylor series

First-order

Price-yield formula

For a function of one variable, \(f(x)\), the Taylor series formula is:

\[ \begin{aligned} f(x+\Delta x)&=f(x)+{\frac{f'(x)}{1!}}\Delta x+{\frac{f''(x)}{2!}}(\Delta x)^{2}+{\frac{f^{(3)}(x)}{3!}}(\Delta x)^{3}+\cdots+{\frac{f^{(n)}(x)}{n!}}(\Delta x)^{n}+\cdots\\ f(x+\Delta x)-f(x)&={\frac{f'(x)}{1!}}\Delta x+{\frac{f''(x)}{2!}}(\Delta x)^{2}+{\frac{f^{(3)}(x)}{3!}}(\Delta x)^{3}+\cdots+{\frac{f^{(n)}(x)}{n!}}(\Delta x)^{n}+\cdots \end{aligned} \]

Using the price-yield formula, the estimated percentage change in price for a change in yield is:

\[ \begin{aligned} P(y+\Delta y)-P(y)&\approx{\frac{P'(y)}{1!}}\Delta y+{\frac{P''(y)}{2!}}(\Delta y)^{2}\\ &\approx -D\Delta y +{\frac{C}{2!}}(\Delta y)^{2} \end{aligned} \]

def pnl_bond(duration, convexity, dy):
    
    duration_pnl = -duration * dy
    convexity_pnl = (convexity / 2) * dy ** 2
    income_pnl = dy
    
    result = pd.DataFrame({
        "total": duration_pnl + convexity_pnl + income_pnl,
        "duration": duration_pnl,
        "convexity": convexity_pnl,
        "income": income_pnl
    })
    
    return result

• https://engineering.nyu.edu/sites/default/files/2021-07/CarWuRF2021.pdf
• https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118267967.app1
• https://www.investopedia.com/terms/c/convexity-adjustment.asp

factor = "DGS10"
duration = 6.5
convexity = 0.65
y = levels_df.ffill()[factor].iloc[-width]

bond_df = pd.DataFrame({
    "duration": duration,
    "convexity": convexity,
    "dy": (levels_df.ffill()[factor].iloc[-width:] - y) / 100
})

attrib_df = pnl_bond(bond_df["duration"], bond_df["convexity"], bond_df["dy"])

Duration-yield formula

The derivative of duration with respect to interest rates gives:

\[ \begin{aligned} \text{Drift}&=\frac{\partial D}{\partial y}=\frac{\partial}{\partial y}\left(-\frac{1}{P}\frac{\partial D}{\partial y}\right)\\ &=-\frac{1}{P}\frac{\partial^{2}P}{\partial y^{2}}+\frac{1}{P^{2}}\frac{\partial P}{\partial y}\frac{\partial P}{\partial y}\\ &=-C+D^{2} \end{aligned} \]

Because of market conventions, use the following formula: \(\text{Drift}=\frac{1}{100}\left(-C\times 100+D^{2}\right)=-C+\frac{D^{2}}{100}\). For example, if convexity and yield are percent then \(\text{Drift}=\left(-0.65+\frac{6.5^{2}}{100}\right)\partial y\times100\) or basis points then \(\text{Drift}=\left(-65+6.5^{2}\right)\partial y\).

def yield_shock(shock, tau, sigma):
    
    result = shock * sigma * np.sqrt(tau)
    
    return result

def duration_drift(duration, convexity, dy):
    
    drift = -convexity + duration ** 2 / 100
    change = drift * dy * 100
    
    result = {
        "drift": drift,
        "change": change
    }
    
    return result

# "Risk Management: Approaches for Fixed Income Markets" (page 45)
factor = "DGS10"
sigma = sd_df[factor].iloc[-1]

duration_df = pd.DataFrame(shocks).rename(columns = {0: "shock"})
duration_df["spot"] = yield_shock(duration_df["shock"], tau, sigma)
duration_df["static"] = duration
duration_df["dynamic"] = duration + \
    duration_drift(duration, convexity, duration_df["spot"])["change"]

Second-order

Black’s formula

A similar formula holds for functions of several variables \(f(x_{1},\ldots,x_{n})\). This is usually written as:

\[ \begin{aligned} f(x_{1}+\Delta x_{1},\ldots,x_{n}+\Delta x_{n})&=f(x_{1},\ldots, x_{n})+ \sum _{j=1}^{n}{\frac{\partial f(x_{1},\ldots,x_{n})}{\partial x_{j}}}(\Delta x_{j})\\ &+{\frac {1}{2!}}\sum_{j=1}^{n}\sum_{k=1}^{n}{\frac{\partial^{2}f(x_{1},\ldots,x_{d})}{\partial x_{j}\partial x_{k}}}(\Delta x_{j})(\Delta x_{k})+\cdots \end{aligned} \]

Using Black’s formula, the estimated change of an option price is:

\[ \begin{aligned} V(S+\Delta S,\sigma+\Delta\sigma,t+\Delta t)-V(S,\sigma,t)&\approx{\frac{\partial V}{\partial S}}\Delta S+{\frac{1}{2!}}{\frac{\partial^{2}V}{\partial S^{2}}}(\Delta S)^{2}+{\frac{\partial V}{\partial \sigma}}\Delta\sigma+{\frac{\partial V}{\partial t}}\Delta t\\ &\approx \Delta_{c,p}\Delta S+{\frac{1}{2!}}\Gamma_{c,p}(\Delta S)^{2}+\nu_{c,p}\Delta\sigma+\Theta_{c,p}\Delta t \end{aligned} \]

• https://quant-next.com/option-greeks-and-pl-decomposition-part-1/

def pnl_option(type, S, K, r, q, tau, sigma, dS, dt, dsigma):
    
    d1 = bs_d1(S, K, r, q, tau, sigma)
    d2 = bs_d2(S, K, r, q, tau, sigma)
    value = bs_value(type, S, K, r, q, tau, sigma, d1, d2)
    delta = bs_delta(type, S, K, r, q, tau, sigma, d1, d2)
    vega = bs_vega(type, S, K, r, q, tau, sigma, d1, d2)
    theta = bs_theta(type, S, K, r, q, tau, sigma, d1, d2)
    gamma = bs_gamma(type, S, K, r, q, tau, sigma, d1, d2)
    
    delta_pnl = delta * dS / value
    gamma_pnl = gamma / 2 * dS ** 2 / value
    vega_pnl = vega * dsigma / value
    theta_pnl = theta * dt / value
    
    result = pd.DataFrame({
        "total": delta_pnl + gamma_pnl + vega_pnl + theta_pnl,
        "delta": delta_pnl,
        "gamma": gamma_pnl,
        "vega": vega_pnl,
        "theta": theta_pnl
    })
    
    return result

factor = "SP500"
type = "call"
S = levels_df.ffill()[factor].iloc[-width]
K = S # * (1 + 0.05)
tau = 1 # = 252 / 252
sigma = sd_df[factor].iloc[-width]

options_df = pd.DataFrame({
    "spot": levels_df.ffill()[factor].iloc[-width:],
    "sigma": sd_df[factor].iloc[-width:]
})
options_df["dS"] = options_df["spot"] - S
options_df["dt_diff"] = (options_df.index - options_df.index[0]).days
options_df["dt"] = options_df["dt_diff"] / options_df["dt_diff"].iloc[-1]
options_df["dsigma"] = options_df["sigma"] - sigma

attrib_df = pnl_option(type, S, K, r, q, tau, sigma,
                       options_df["dS"], options_df["dt"], options_df["dsigma"])

Ito’s lemma

For a given diffiusion \(X(t, w)\) driven by:

\[ \begin{aligned} dX_{t}&=\mu_{t}dt+\sigma_{t}dB_{t} \end{aligned} \]

Then proceed with the Taylor series for a function of two variables \(f(t,x)\):

\[ \begin{aligned} df&={\frac{\partial f}{\partial t}}dt+{\frac{\partial f}{\partial x}}dx+{\frac{1}{2}}{\frac{\partial^{2}f}{\partial x^{2}}}dx^{2}\\ &={\frac{\partial f}{\partial t}}dt+{\frac{\partial f}{\partial x}}(\mu_{t}dt+\sigma_{t}dB_{t})+{\frac{1}{2}}{\frac{\partial^{2}f}{\partial x^{2}}}\left(\mu_{t}^{2}dt^{2}+2\mu_{t}\sigma _{t}dtdB_{t}+\sigma_{t}^{2}dB_{t}^{2}\right)\\ &=\left({\frac{\partial f}{\partial t}}+\mu_{t}{\frac{\partial f}{\partial x}}+{\frac{\sigma _{t}^{2}}{2}}{\frac{\partial ^{2}f}{\partial x^{2}}}\right)dt+\sigma_{t}{\frac{\partial f}{\partial x}}dB_{t} \end{aligned} \]

Note: set the \(dt^{2}\) and \(dtdB_{t}\) terms to zero and substitute \(dt\) for \(dB^{2}\).

Geometric Brownian motion

The most common application of Ito’s lemma in finance is to start with the percent change of an asset:

\[ \begin{aligned} \frac{dS}{S}&=\mu_{t}dt+\sigma_{t}dB_{t} \end{aligned} \]

Then apply Ito’s lemma with \(f(S)=log(S)\):

\[ \begin{aligned} d\log(S)&=f^{\prime}(S)dS+{\frac{1}{2}}f^{\prime\prime}(S)S^{2}\sigma^{2}dt\\ &={\frac {1}{S}}\left(\sigma SdB+\mu Sdt\right)-{\frac{1}{2}}\sigma^{2}dt\\ &=\sigma dB+\left(\mu-{\tfrac{\sigma^{2}}{2}}\right)dt \end{aligned} \]

It follows that:

\[ \begin{aligned} \log(S_{t})-\log(S_{0})=\sigma dB+\left(\mu-{\tfrac{\sigma^{2}}{2}}\right)dt \end{aligned} \]

Exponentiating gives the expression for \(S\):

\[ \begin{aligned} S_{t}=S_{0}\exp\left(\sigma B_{t}+\left(\mu-{\tfrac{\sigma^{2}}{2}}\right)t\right) \end{aligned} \]

This provides a recursive procedure for simulating values of \(S\) at \(t_{0}<t_{1}<\cdots<t_{n}\):

\[ \begin{aligned} S(t_{i+1})&=S(t_{i})\exp\left(\sigma\sqrt{t_{i+1}-t_{i}}Z_{i+1}+\left[\mu-{\tfrac{\sigma^{2}}{2}}\right]\left(t_{i+1}-t_{i}\right)\right) \end{aligned} \]

where \(Z_{1},Z_{2},\ldots,Z_{n}\) are independent standard normals.

def sim_gbm(n_sim, S, mu, sigma, dt):
    
    result = S * np.exp(np.cumsum(sigma * np.sqrt(dt) * np.random.normal(size = n_sim)) + \
                        (mu - 0.5 * sigma ** 2) * dt)
    
    return result

This leads to an algorithm for simulating a multidimensional geometric Brownian motion:

\[ \begin{aligned} S_{k}(t_{i+1})&=S_{k}(t_{i})\exp\left(\sqrt{t_{i+1}-t_{i}}\sum_{j=1}^{d}{A_{kj}Z_{i+1,j}}+\left[\mu_{k}-{\tfrac{\sigma_{k}^{2}}{2}}\right]\left(t_{i+1}-t_{i}\right)\right) \end{aligned} \]

where \(A\) is the Cholesky factor of \(\Sigma\), i.e. \(A\) is any matrix for which \(AA^\mathrm{T}=\Sigma\).

def sim_multi_gbm(n_sim, S, mu, sigma, dt):
    
    n_cols = sigma.shape[1]
    
    Z = np.random.normal(size = n_sim * n_cols).reshape((n_sim, n_cols))
    X = np.sqrt(dt) * Z @ np.linalg.cholesky(sigma).T + (mu - 0.5 * np.diag(sigma)) * dt
    
    result = S * np.exp(X.cumsum(axis = 0))
    
    return np.asmatrix(result)

• https://arxiv.org/pdf/0812.4210.pdf
• https://quant.stackexchange.com/questions/15219/calibration-of-a-gbm-what-should-dt-be
• https://stackoverflow.com/questions/36463227/geometrical-brownian-motion-simulation-in-r
• https://quant.stackexchange.com/questions/25219/simulate-correlated-geometric-brownian-motion-in-the-r-programming-language
• https://quant.stackexchange.com/questions/35194/estimating-the-historical-drift-and-volatility/

S = [1] * len(factors)
sigma = np.cov(returns_df["returns"].dropna().T, ddof = 1) * scale["periods"]
mu = np.array(returns_df["returns"].dropna().mean()) * scale["periods"]
mu = mu + np.diag(sigma) / 2 # drift
dt = 1 / scale["periods"]

mu_ls = []
sigma_ls = []

for i in range(10000): # "TypeError: 'float' object cannot be interpreted as an integer"
  
    # assumes underlying stock price follows geometric Brownian motion with constant volatility
    levels_sim = pd.DataFrame(sim_multi_gbm(width + 1, S, mu, sigma, dt))
    returns_sim = np.log(levels_sim).diff()

    mu_sim = returns_sim.mean() * scale["periods"]
    sigma_sim = returns_sim.std() * np.sqrt(scale["periods"])

    mu_ls.append(mu_sim)
    sigma_ls.append(sigma_sim)

mu_df = pd.DataFrame(mu_ls)
sigma_df = pd.DataFrame(sigma_ls)

pd.DataFrame({
    "empirical": np.array(returns_df["returns"].dropna().mean()) * scale["periods"],
    "theoretical": mu_df.mean()
})

   empirical  theoretical
0   0.107832     0.109578
1   0.009068     0.009240
2  -0.000834    -0.000760
3   0.003054     0.002855

pd.DataFrame({
    "empirical": np.sqrt(np.diag(sigma)),
    "theoretical": sigma_df.mean()
})

   empirical  theoretical
0   0.180453     0.180255
1   0.062865     0.062732
2   0.008495     0.008482
3   0.016955     0.016931

Vasicek model

# assumes interest rates follow mean-reverting process with stochastic volatility

Newton’s method

Implied volatility

Newton’s method (main idea is also from a Taylor series) is a method for finding approximations to the roots of a function \(f(x)\):

\[ \begin{aligned} x_{n+1}=x_{n}-{\frac{f(x_{n})}{f'(x_{n})}} \end{aligned} \]

To solve \(V(\sigma_{n})-V=0\) for \(\sigma_{n}\), use Newton’s method and repeat until \(\left|\sigma_{n+1}-\sigma_{n}\right|<\varepsilon\):

\[ \begin{aligned} \sigma_{n+1}=\sigma_{n}-{\frac{V(\sigma_{n})-V}{V'(\sigma_{n})}} \end{aligned} \]

def implied_vol_newton(params, type, S, K, r, q, tau):
    
    target0 = 0
    sigma = params["sigma"]
    sigma0 = sigma
    
    while (abs(target0 - params["target"]) > params["tol"]):
        
        d1 = bs_d1(S, K, r, q, tau, sigma0)
        d2 = bs_d2(S, K, r, q, tau, sigma0)
        
        target0 = bs_value(type, S, K, r, q, tau, sigma0, d1, d2)
        d_target0 = bs_vega(type, S, K, r, q, tau, sigma0, d1, d2)
        
        sigma = sigma0 - (target0 - params["target"]) / d_target0
        sigma0 = sigma
        
    return sigma

• http://www.aspenres.com/documents/help/userguide/help/bopthelp/bopt2Implied_Volatility_Formula.html
• https://books.google.com/books?id=VLi61POD61IC&pg=PA104

S = levels_df.ffill()[factor].iloc[-1]
K = S # * (1 + 0.05)
sigma = sd_df[factor].iloc[-1] # overrides matrix
start1 = 0.2

d1 = bs_d1(S, K, r, q, tau, sigma)
d2 = bs_d2(S, K, r, q, tau, sigma)
target1 = bs_value(type, S, K, r, q, tau, sigma, d1, d2)
params1 = {
    "target": target1,
    "sigma": start1,
    "tol": 1e-4 # np.finfo(float).eps
}

implied_vol_newton(params1, type, S, K, r, q, tau) 

0.12287141700380808

Yield-to-maturity

def yld_newton(params, cash_flows):
    
    target0 = 0
    yld0 = params["cpn"] / params["freq"]
    yld = yld0 # assignment to `yield` variable is not possible
    
    while (abs(target0 - params["target"]) > params["tol"]):
        
        target0 = 0
        d_target0 = 0
        dd_target0 = 0
        
        for i in range(len(cash_flows)):
          
            t = i + 1
          
            # present value of cash flows
            target0 += cash_flows[i] / (1 + yld0) ** t
            
            # first derivative of present value of cash flows
            d_target0 -= t * cash_flows[i] / (1 + yld0) ** (t + 1) # use t for Macaulay duration
            
            # second derivative of present value of cash flows
            dd_target0 -= t * (t + 1) * cash_flows[i] / (1 + yld0) ** (t + 2)
        
        yld = yld0 - (target0 - params["target"]) / d_target0
        yld0 = yld
        
    result = {
        "price": target0,
        "yield": yld * params["freq"],
        "duration": -d_target0 / params["target"] / params["freq"],
        "convexity": -dd_target0 / params["target"] / params["freq"] ** 2
    }
        
    return result

• https://www.bloomberg.com/markets/rates-bonds/government-bonds/us
• https://quant.stackexchange.com/a/61025
• https://pages.stern.nyu.edu/~igiddy/spreadsheets/duration-convexity.xls

target2 = 0.9928 * 1000 # present value
start2 = 0.0438 # coupon
cash_flows = [start2 * 1000 / 2] * 10 * 2
cash_flows[-1] += 1000

params2 = {
    "target": target2,
    "cpn": start2,
    "freq": 2,
    "tol": 1e-4 # np.finfo(float).eps
}

yld_newton(params2, cash_flows)

{'price': 992.8000005704454, 'yield': 0.044700757710159994, 'duration': 8.0165959605043, 'convexity': 76.68109754287481}

Optimization

from scipy.optimize import minimize

Implied volatility

If the derivative is unknown, try optimization:

def implied_vol_obj(param, type, S, K, r, q, tau, target):
    
    d1 = bs_d1(S, K, r, q, tau, param)
    d2 = bs_d2(S, K, r, q, tau, param)
    target0 = bs_value(type, S, K, r, q, tau, param, d1, d2)
    
    result = abs(target0 - target)
    
    return result

def implied_vol_optim(param, type, S, K, r, q, tau, target):
    
    result = minimize(implied_vol_obj, param, args = (type, S, K, r, q, tau, target))
    
    return result.x.item()

implied_vol_optim(start1, type, S, K, r, q, tau, target1)

0.12287141451913787

Yield-to-maturity

def yld_obj(param, cash_flows, target):
    
    target0 = 0
    
    for i in range(len(cash_flows)):
        
        target0 += cash_flows[i] / (1 + param) ** (i + 1)
    
    target0 = abs(target0 - target)
    
    return target0
  
def yld_optim(params, cash_flows, target):
    
    result = minimize(yld_obj, params["cpn"], args = (cash_flows, target))
    
    return result.x.item() * params["freq"]

yld_optim(params2, cash_flows, target2)

0.04470075001603727