A volatilidade implícita é uma estimativa da variabilidade futura do ativo subjacente ao contrato de opções. Essa estimativa é derivada do modelo Black-Scholes.

Para calcular a Volatilidade Implícita (VI), primeiro, vamos calcular a VI pela aproximação de Brenner-Subrahmanyam (1988), esse valor será utilizado como 'chute inicial' no Newton-Raphson (publicação baseada nas respostas em Quantitative Finance Stack Exchange; o exemplo é por minha conta):

$$\sigma\approx\sqrt{\frac{2\pi}{T}}\cdot\frac{C}{S_t},$$

no qual:

\(T\): Vencimento da opção,

\(C\): preço da call,

\(S_t\): Preço atual da PETR3.

import numpy as np

def sig0(S,T,C):
    sig = np.sqrt(2 * np.pi / T) * C / S
    return(sig)

T = 32/360
C = 1.19
S = 22.65
sigma0 = sig0(S = S, T = T, C = C)
sigma0
# 0.44171785082871096

Pela aproximação acima, a VI é 44.17%. Agora, vamos calcular a VI usando a equação de BS e Vega - a derivada da função BS com relação a \(\sigma\):

$$\nu = S_t\mathcal{N}'(d_1)\sqrt{T},$$

$$d_1=\frac{\ln\frac{S_t}{K}+(r+\frac{\sigma^2}{2}T)}{\sigma\sqrt{T}},$$

no qual:

\(K\): Strike da call,

\(r\): Taxa de juros livre de risco,

\(\mathcal{N}'(\cdot)\): Derivada da função distribuição acumulada da distribuição normal.

import numdifftools as nd

def vega(S,K,r,T,sigma):
    d1 = (np.log(S / K) + (r + T * 0.5 * sigma**2)) / 
    (sigma * np.sqrt(T))
    fd = nd.Derivative(norm.cdf)
    v = S * fd(d1) * np.sqrt(T)
    return(v)

K = 22.46
r = 0.0265
T = 32/360
sigmaH = 0.3319 # Desvio-padrão da ação a.a.

vega(S = S, K = K, r = r, T = T, sigma = sigmaH)
# 2.4844979800506612

O Vega mede a sensibilidade do preço da opção em relação à volatilidade; quanto varia o preço da opção quando a volalitilidade em termos anualizados varia 1 ponto percentual.

# Definindo a função BS

from scipy.stats import norm

def BS(S,K,r,T,sigma):
    d1 = (np.log(S / K) + (r + T * 0.5 * sigma**2)) / 
    (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    C = norm.cdf(d1) * S - norm.cdf(d2) * K * np.exp(-r * T)
    return(C)

Podemos usar o método de Newton-Raphson para calcular qual seria a VI:

$$\sigma_{n+1}=\sigma_n-\frac{BS(\sigma_n) - S_t}{\nu(\sigma_n)}$$

sigma1 = sigma0 - (( BS(S = S,K = K,r = r,T = T,sigma = sigma0) - C) / 
                   vega(S = S, K = K, r = r, T = T, sigma = sigma0))

sigma1
# 0.4030102497124105

Vamos repetir o passo anterior até a diferença em módulo entre o primeiro termo e o segundo termo seja menor/igual a 0.00001, isto é: \(| \sigma_{n+1} - \sigma_n |\le0.00001 \):

diff = 1
i = 1
sigma0 = 0.44171785082871096
while diff >= 1e-5:
    sigma1 = sigma0 - (( BS(S = S,K = K,r = r,T = T,sigma = sigma0) - C) /
                  vega(S = S, K = K, r = r, T = T, sigma = sigma0))
    diff = abs(sigma1 - sigma0)
    sigma0 = sigma1
    i += 1

[diff,i,sigma0]

# [8.692677371247015e-07, 6, 0.40529296748508575]

Portanto, a VI calculada por esse método é 40.52%. A diferença na aproximação ficou em torno do que havíamos pedido para a função calcular, isto é, menor que 1e-5. Além disso, o algoritmo precisou de 6 iterações até convergir.

Até a próxima!