calibração

Superfície SVI

Superfície SVI

Superfície SVI, ou somente SSVI é uma generalização do modelo SVI de (Gatheral 2004) que busca solucionar o problema de restrição dos parâmetros do modelo para evitar a presença de arbitragem do tipo borboleta em um dado smile. Este modelo foi proposto por (Gatheral and Jacquier 2014) e extende o modelo SVI original apresentando duas outras parametrizações equivalentes e então o modelo para superfícies propriamente dito.

Reparametrizações equivalentes

Existem duas outras formas de se apresentar um modelo SVI que são equivalentes a parametrização RAW já apresentada. Estas são as parametrizações “Natural” e “Jump-Wings” que são apresentadas abaixo.

Para um dado conjunto de parâmetros \(\chi_N=\{\Delta, \mu, \rho, \omega, \zeta\}\) a parametrização natural de um SVI é dada por:

\[\begin{equation} w(k; \chi_N)=\Delta+\frac{\omega}{2}\left\lbrace 1+\zeta\rho(k-\mu)+\sqrt{(\zeta(k-\mu)+\rho)^2+(1-\rho^2)} \right\rbrace \tag{1} \end{equation}\]

onde \(\omega\geq 0\), \(\Delta, \mu \in \mathbb R\), \(|\rho|<1\) e \(\zeta>0\). A correspondência entre as parametrizações raw e natural é dada pelo seguinte mapeamento e seu inverso:

\[\begin{equation} (a, b, \rho, m, \sigma)=\left(\Delta+\frac{\omega}{2}(1-\rho^2), \frac{\omega\zeta}{2}, \rho, \mu-\frac{\rho}{\zeta}, \frac{\sqrt{1-\rho^2}}{\zeta}\right) \tag{2} \end{equation}\] \[\begin{equation} (\Delta, \mu, \rho, \omega, \zeta)=\left(a-\frac{\omega}{2}(1-\rho^2), m+\frac{\rho\sigma}{\sqrt{1-\rho^2}}, \rho, \frac{2b\sigma}{\sqrt{1-\rho^2}}, \frac{\sqrt{1-\rho^2}}{\sigma}\right) \tag{3} \end{equation}\]

A desvantagem destas parametrizações é que o valor de seus parâmetros não são intuitivos para os traders, eles não carregam estes valores em sua memória durante a negociação. Valores característicos de uma superfície de volatilidade implícita que traders têm em mente são, por exemplo, volatilidade ATM, skew de volatilidade ATM e assíntotas. Desta forma a parametrização Jump-Wings é útil, pois relaciona estes valores típicos aos parâmetros raw de um SVI.

A parametrização JW é dada em termos da variância implícita (e não da variância total) e portanto existe uma dependência explícita do tempo em sua formulação. Para um dado tempo até a expiração, \(\tau\), o conjunto de parâmetros \(\chi_{J}=\{v_\tau, \psi_\tau, p_\tau, c_\tau, \tilde v_\tau\}\) é definido pelas seguintes equações a partir dos parâmetros raw:

\[\begin{align} v_\tau&=\frac{a+b\{-\rho m + \sqrt{m^2+\sigma^2}\}}{\tau},\\ \psi_\tau&=\frac{1}{\sqrt{w_\tau}}\frac{b}{2}\left(-\frac{m}{\sqrt{m^2+\sigma^2}}+\rho\right),\\ p_\tau&=\frac{1}{\sqrt{w_\tau}}b(1-\rho),\\ c_\tau&=\frac{1}{\sqrt{w_\tau}}b(1+\rho),\\ \tilde v_\tau&=\frac{1}{\tau}\left(a+b\sigma\sqrt{1-\rho^2}\right) \tag{4} \end{align}\]

onde \(w_\tau=v_\tau \tau\) relaciona a variância total ATM com a variância ATM. Os parâmetros possuem as seguintes interpretações: \(v_\tau\) é a variância ATM, \(\psi_\tau\) o skew ATM, \(p_\tau\) a inclinação da asa esquerda (puts), \(c_\tau\) a inclinação da asa direita (calls) e \(\tilde v_\tau\) é a variância implícita mínima.

A figura 1 apresenta uma esquematização destes parâmetros sobre um smile fictício para melhor compreensão.

Interpretação dos parâmetros de um SVI-JW.

Figura 1: Interpretação dos parâmetros de um SVI-JW.

As relações inversas que trazem uma parametrização JW para uma raw, assumindo que \(m \neq 0\) são:

\[\begin{align} b&=\frac{\sqrt{w_\tau}}{2}(c_\tau+p_\tau),\\ \rho&=1-\frac{p_\tau\sqrt{w_\tau}}{b},\\ a&=\tilde v_\tau \tau-b\sigma\sqrt{1-\rho^2},\\ m&=\frac{(v_\tau-\tilde v_\tau)\tau}{b\left\lbrace-\rho+sign(\alpha)\sqrt{1+\alpha^2}-\alpha\sqrt{1-\rho^2}\right\rbrace},\\ \sigma&=\alpha m. \tag{5} \end{align}\]

onde as variáveis auxiliares são definidas da seguinte forma: \(\beta:=\rho-(2\psi_\tau\sqrt{w_\tau})/b\) e \(\alpha:=sign(\beta)\sqrt{1/\beta^2 – 1}\), com \(\beta \in [-1, 1]\) para garantir a convexidade do smile.

Se \(m=0\), então as equações para \(a, b, \rho\) se mantêm, porém, \(\sigma = (v_\tau \tau-a)/b\). Desta forma temos relações entre as três parametrizações SVI, sendo possível navegar entre elas com tranquilidade. Um trader pode verificar no mercado os valores dos parâmetros JW e traduzi-los para raw e simular o smile ou fazer o caminho reverso, calibrar uma fatia da superfície com parâmetros raw, traduzi-los para JW e apresentar para sua mesa, onde todos conseguirão interpretar os valores a que estão habituados.

Superfície SVI

Uma SSVI surge como uma extensão à parametrização natural de um SVI, e fornece em uma única equação, a possibilidade de parametrizar uma superfície de volatilidade implícita por inteiro. É necessário, antes de mais nada, fazer algumas definições preliminares. Defina a variância total implícita no dinheiro (ATM) como \(\theta_\tau:=\sigma_{BS}^2(0, \tau)\tau\) e \(\lim\limits_{\tau\rightarrow 0}\theta_\tau = 0\).

Definição 1 Seja \(\varphi\) uma função suave em \(\mathbb R_+^*\mapsto \mathbb R_+^*\) tal que \(\lim\limits_{\tau\rightarrow 0}\theta_\tau \varphi(\theta_\tau)\) exista em \(\mathbb R\). Uma SSVI é definida por:

\[\begin{equation} w(k, \theta_\tau)=\frac{\theta_\tau}{2}\left\lbrace 1+\rho\varphi(\theta_\tau)k+\sqrt{(\varphi(\theta_\tau)k+\rho)^2+(1-\rho^2)} \right\rbrace \tag{6} \end{equation}\]

Veja que na representação SSVI, engenhosamente os autores substituíram a dimensão de tempo, \(\tau\), típica de superfícies de volatilidde, pela variância total ATM. Com esta representação, eles conseguiram derivar as condições necessárias para a ausência de arbitragem estática na superfície e sempre que necessário, é possível retornar a dimensão de tempo no calendário.

Agora é necessário definir a função \(\varphi\), que então será substituída na equação (6) da definição 1 com seus próprios parâmetros e teremos por fim uma função \(w(k, \theta_\tau; \chi_{ssvi})\) que poderá ser calibrada para o conjunto de parâmetros da SSVI, \(\chi_{ssvi}\) contra os dados observados no mercado. No fundo, qualquer função que obedeça as condições impostas na definição 1 pode ser utilizada, entretanto os autores apresentam dois tipos de função que condizem com as observações empíricas.

Heston

Considere a função \(\varphi\) definida por:

\[\begin{equation} \varphi(\theta)\equiv\frac{1}{\gamma\theta}\left\lbrace 1-\frac{1-e^{-\gamma\theta}}{\gamma\theta}\right\rbrace \tag{7} \end{equation}\]

com \(\gamma > 0\). Esta função recebeu este nome pois, seu skew na variância implícita é compatível com aquele previsto no modelo de Heston.

Lei de potência

A parametrização da função \(\varphi\) como uma lei de potência é primeiramente considerada da seguinte forma: \(\varphi(\theta)=\eta\theta^{-\gamma}\) com \(\eta > 0\) e \(0< \gamma<1\). Entretanto, esta parametrização apresenta algumas limitações com relação a arbitragem do tipo borboleta e então é proposta a seguinte forma funcional:

\[\begin{equation} \varphi(\theta)=\frac{\eta}{\theta^\gamma(1+\theta)^{1-\gamma}} \tag{8} \end{equation}\]

que é garantida não possuir arbitragem estática dado que \(\eta(1+|\rho|)\leq 2\).

Condições de não-arbitragem

As condições para ausência de arbitragem estática para uma SSVI são colocadas através de dois teoremas (4.1 e 4.2) e provadas no artigo de (Gatheral and Jacquier 2014), dos quais resulta o seguinte corolário.

Corolário 1 A superfície SVI definida em 1 está livre de arbitragem estática se as seguintes condições são satisfeitas:

  1. \(\partial_\tau\theta_\tau\geq 0, \text{ para todo } \tau > 0\)
  2. \(0\leq \partial_\theta(\theta\varphi(\theta))\leq\frac{1}{\rho^2}\left(1+\sqrt{1-\rho^2}\right)\varphi(\theta), \text{ para todo } \theta>0\)
  3. \(\theta\varphi(\theta)(1+|\rho|)<4, \text{ para todo } \theta>0\)
  4. \(\theta\varphi(\theta)^2(1+|\rho|)\leq 4, \text{ para todo } \theta>0\)

Onde os dois primeiros itens dizem respeito a ausência de arbitragem de calendário, enquanto que os seguintes são exigências para a superfície estar livre de arbitragem do tipo borboleta.

Para uma função \(\varphi\) do tipo Heston, as condições nos seus parâmetros são: \(\gamma>0\) que garante o atendimentos as condições 1 e 2 do corolário 1 e \(\gamma\geq(1+|\rho|)/4\) para garantir a ausência de arbitragem borboleta, a qual subsume a primeira condição.

Para a função do tipo lei de potência dada em (8) são necessárias as condições \(0< \gamma<1\) e \(\eta(1+|\rho|)\leq 2\)

A maneira típica de impor as restrições dos parâmetros no momento da calibração do modelo é inserir uma penalidade na função objetivo, quando a restrição é violada. Por exemplo, consideramos a restrição de inequalidade para a lei de potência, \(\eta(1+|\rho|)-2\leq 0\). No momento da calibração, nós devemos calcular o valor desta expressão e se o seu resultado for maior que zero, uma penalidade é somada a função objetivo da calibração (em geral a soma dos quadrados dos erros).

Densidade neutra ao risco

Já comentamos sobre a distribuição neutra ao risco implícita no preço de opções que pode ser obtida através da fórmula de (Breeden and Litzenberger 1978), assim como já foi introduzida uma certa função \(g(k)\) que desempenha um papel fundamental para garantir a ausência de arbitragem do tipo borboleta em smiles de variância total. O fato é que, para garantir a ausência de arbitragem a função \(g(k)\) deve ser não negativa em todo o seu suporte \(k \in \mathbb R\). Ao mesmo tempo, a definição de ausência de arbitragem borboleta se confunde com o fato que a densidade neutra ao risco deve também ser não negativa, caso contrário não seria uma densidade de probabilidade. Logo percebe-se a estreita relação entre a função \(g(k)\) e a densidade neutra ao risco implícita no preço de opções. De fato, a fórmula de Breeden-Litzenberger nos fornece:

\[\begin{equation} p(k)=\left.\frac{\partial^2C_B(k, w(k))}{\partial K^2}\right|_{K=Fe^k} \end{equation}\]

que após fazer as derivações da equação de Black e as substituições e manipulações algébricas com as reparametrizações do modelo B&S, resulta em:

\[\begin{equation} p(k)=\frac{g(k)}{\sqrt{2\pi w(k)}}\exp\left(-\frac{d_2(k)^2}{2}\right) \tag{9} \end{equation}\]

E para relembrarmos, a função \(g(k)\) é dada por:

\[\begin{equation} g(k)=\left(1-\frac{kw\prime(k)}{2w(k)}\right)^2-\frac{w\prime(k)^2}{4}\left(\frac{1}{w(k)}+\frac{1}{4}\right)+\frac{w\prime\prime(k)}{2} \tag{10} \end{equation}\]

Ou seja, uma vez parametrizado um SVI para um dado prazo de maturidade, \(\tau\), se torna simples a tarefa de extrair a densidade implícita. Possuímos a função \(w(k)\) e suas derivadas e certamente \(d_2(k)\), bastando portanto, aplicar a equação (9) para extrair importante informação do mercado de opções.

Superfície de volatilidade local

De maneira semelhante ao procedimento realizado com a densidade implícita, também é possível, uma vez parametrizada a SSVI, derivar a superfície de volatilidade local através da equação de (Dupire 1994). Esta equação toma a seguinte forma funcional:

\[\begin{equation} \sigma_L^2(K, \tau)=\frac{\partial_\tau C_B(K, \tau)}{\frac{1}{2}K^2\partial_{KK}C_B(K, \tau)} \tag{11} \end{equation}\]

Novamente tomando as derivadas da equação de Black e fazendo as substituições necessárias1 chegamos a relação entre a superfície SVI e variância local.

\[\begin{equation} \sigma^2_L(k, \tau)=\frac{\partial_\tau w(k, \theta_\tau)}{g(k, w(k, \theta_\tau))} \tag{12} \end{equation}\]

De forma bastante simplista, a superfície de volatilidade local pode ser entendida como aquela superfície de volatilidades instantâneas para o ativo subjacente, \(\sigma(S, t)\) que depende tanto do nível de preço deste ativo quanto do tempo, e fornece a previsão do mercado para a volatilidade instantânea de \(S\) dado que este ativo hoje está precificado em \(S_0\) e no tempo futuro \(\tau\) estará em \(K\). Fazendo uma analogia com o mercado de juros, a superfície local está para a curva forward assim como a superfície implícita está para a curva de juros, numa aproximação.

A superfície local é muito utilizada para precificar opções exóticas, aquelas que possuem perfil de payoff distinto das opções europeias e podem ou não terem seus resultados atrelados ao caminho seguido pelo preço do ativo objeto durante a vida da opção. Nestes casos, a precificação se dá através da incorporação da superfície local estipulando um valor de volatilidade instantânea para cada possível combinação \((S_t, t)\) em geral em um ambiente de simulação de Monte Carlo.

No caso da equação (12) temos uma nova derivada a ser computada, \(\partial_\tau w(k, \theta_\tau)\), que só poderia ser feita de forma analítica caso \(\theta_\tau\) fosse parametrizada explicitamente. Nem sempre este é o caso, já que \(\theta_\tau\) são os valores de variância total ATM implícitas, ou seja, observa-se no mercado apenas alguns pontos de \(\theta_\tau\), pontos estes que podem estar espaçados em intervalos diferentes de tempo. A solução para esta derivação é interpolar estes pontos, nem que seja uma simples interpolação linear, e então fazer a derivação de forma numérica através de um método de diferenças finitas. Note que na interpolação deve-se garantir a condição de não arbitragem de calendário, \(\partial_\tau\theta_\tau\geq 0\).

Calibração

Vamos retomar nosso exemplo de superfície de volatilidade do artigo anterior, porém agora com todas as datas de expiração, ou seja, a superfície completa. Os códigos R apresentados abaixo ajudam na compreensão do procedimento.

library(readr)
library(dplyr)
library(purrr)
library(kableExtra)
library(ggplot2)
library(ggthemes)
library(plot3D)
source('svi.R')
source('ssvi.R')

Primeiramente foram carregados os pacotes necessários para a leitura e manipulação dos dados (readr, dplyr e purrr) assim como os pacotes de visualização (kableExtra, ggplot2, ggthemes e plot3D). Em seguida os arquivos svi.R e ssvi.R são implementações do Clube de Finanças para as funções necessárias para calibrar uma SSVI.

Carregados os pacotes e as funções, deve-se carregar os dados da superfície de volatilidade e organizá-los da forma que necessitamos.

ssvi_data <- read_csv("../../static/input/IV_Raw_Delta_surface.csv",
                     col_types = cols(date = col_date(format = "%m/%d/%Y"))) %>% 
  mutate(tau = period / 365,
         theta = theta_vec(moneyness, tau, iv)) %>% 
  rename(k = moneyness)

Para calibrar os parâmetros de uma SSVI, com a função \(\varphi\) do tipo lei de potência, necessitaremos dos dados de moneyness (\(k\)), vencimento (\(\tau\)) e volatilidade implícita (\(iv\)) que será internamente convertida em variância total implícita (\(w\)). Dentro da função de ajuste dos parâmetros fit_ssvi() a variância total ATM implícita \(\theta_\tau\) é calculada através de interpolação spline para cada uma das fatias da superfície, pois, nossos dados não necessariamente possuem esta observação. Cabe ressaltar também que os dados estão organizados de forma tidy, sendo portanto, todos os argumentos passados para as funções na forma de vetores, não sendo necessário gerar uma matriz contendo o grid \((k, \tau)\).

k <- ssvi_data$k
tau <- ssvi_data$tau
iv <- ssvi_data$iv
theta <- ssvi_data$theta

set.seed(12345)
powerlaw_par <- fit_ssvi(tau, k, iv, "powerlaw")
kable(powerlaw_par,
      caption = "Parâmetros da SSVI-power-law estimados.",
      col.names = "SSVI-PL")
Tabela 1: Parâmetros da SSVI-power-law estimados.
SSVI-PL
rho -0.6479238
gamma 0.4926757
eta 0.8607807

Podemos rapidamente checar se estes parâmetros estimados geram uma superfície livre de arbitragem.

paste("Parametrização livre de arbitragem borboleta? :", 
      ssvi_butterfly_cons(powerlaw_par, "powerlaw") <= 0)
## [1] "Parametrização livre de arbitragem borboleta? : TRUE"

Com os parâmetros estimados, é simples plotar todos os smiles que compõe a superfície e verificar visualmente o ajuste. Vamos plotar a variância total em função do moneyness e vencimentos, pois neste gráfico se verifica a condição de arbitragem de calendário simplesmente através do fato que as curvas geradas não devem se cruzar.

plt_df <- ssvi_data %>% 
  mutate(w_pl = ssvi_fun(powerlaw_par, theta, k)) %>% 
  select(k, tau, w_pl, iv) %>% 
  filter(tau < 0.5)

ggplot(plt_df, aes(x = k, y = w_pl)) + 
  geom_line(aes(color = as.factor(format(tau, digits = 3)))) +
  geom_point(aes(y = iv^2 * tau)) +
  guides(color = guide_legend(title = "Vencimento")) +
  labs(title = "SSVI Lei de Potência",
       x = "Forward log-moneyness (k)",
       y = "Variância total implícita (w)",
       caption = "Elaborado por Rafael Bressan para o Clube de Finanças.") +
#  scale_y_continuous(labels = scales::percent) +
  scale_color_viridis_d() +
  theme_economist_white()

Foram retirados os vencimentos mais longos apenas para uma melhor visualização da parte mais curta da superfície, que em geral é a mais complicada de se calibrar através de SVI. O resultado parece interessante, com as variâncias ao redor do dinheiro tendo boa aderência aos dados observados, porém para puts muito fora do dinheiro, a parametrização aparenta subestimar o valor da variância total e portanto estaria subprecificando as opções.

Iremos agora verificar a densidade da distribuição implícita para o período de vencimento de 90 dias. Utiliza-se para tanto a equação (9) onde \(w(k)\) (e suas derivadas e o próprio \(d_2\)) será agora calculado com base nos parâmetros da tabela 1 para um vetor de moneyness mais amplo e denso. O resultado pode ser observado através da figura 2

thetadens <- ssvi_data %>% 
  filter(period == 90) %>% 
  pull(theta) %>% 
  `[`(1)
kdens <- seq(-0.5, 0.3, length.out = 100)
dens <- ssvi_density(powerlaw_par, thetadens, kdens, "powerlaw")
dens_tbl <- tibble(kdens = kdens, dens = dens)

ggplot(dens_tbl, aes(x = kdens, y = dens)) + 
  geom_line() +
  labs(title = "Densidade neutra ao risco SSVI",
     x = "Forward log-moneyness (k)",
     y = "Densidade",
     caption = "Elaborado por Rafael Bressan para o Clube de Finanças.") +
  scale_color_viridis_d() +
  theme_economist_white()
Densidade implícita estimada a partir da SSVI. Presença de assimetria e leptocurtose a esquerda.

Figura 2: Densidade implícita estimada. Presença de assimetria e leptocurtose a esquerda.

Este é um gráfico interessante, que nos mostra exatamente o que se espera de um smile típico de equities que possui skew negativo. Percebemos como a densidade de probabilidades é assimétrica, com a cauda esquerda muito mais longa, refletindo o sentimento de mercado de uma maior probabilidade de grandes quedas no preço do ativo objeto que altas equivalentes. Sempre que verificamos um smile com skew negativo, como é o presente caso, a distribuição de probabilidades é assimétrica a esquerda.

Vamos conferir se a área sob esta curva de densidade integra aproximadamente 1, como deve ser.

area <- integrate(function(x) ssvi_density(powerlaw_par, thetadens, x, "powerlaw"),
                  lower = kdens[1],
                  upper = kdens[length(kdens)])
paste("Área sob a curva de densidade é: ", area$value)

[1] “Área sob a curva de densidade é: 0.999302879362191”

Chegou o momento de inferirmos a superfície de volatilidade local a partir de nossa parametrização SSVI. O método será a aplicação direta da equação (12) pois, com a parametrização realizada, dispomos de todos os dados necessários para seu cômputo.

Antes porém, devemos observar uma peculiaridade dos nossos dados. O site ivolatility.com fornece dados alinhados por Delta, ou seja, para cada período de vencimento temos um mesmo conjunto de Deltas, mas não de moneyness. Os valores deste último irá variar conforme o vencimento, uma vez que pela definição de forward log-moneyness que vimos utilizando é dependente do tempo para maturidade da opção. Desta forma precisamos gerar um novo grid \((k, \tau)\) para plotar nossa superfície de volatilidade local.

Para tanto, criaremos uma sequência uniforme de valores de \(k\) e tomaremos os valores que já possuímos de \(\tau\) e \(\theta_\tau\) recombinando-os através da função expand.grid() para gerar o data frame que precisamos.

Feito isso, é apenas questão de aplicar a equação (12) e criar uma matriz (pois assim pede a função surf3D()) com os valores da volatilidade local.

kloc <- seq(-.4, 0.4, length.out = 17)
utau <- unique(tau)
utheta <- unique(theta)
names(utheta) <- utau  # utheta will be a lookup vector
grid_df <- expand.grid(kloc, utau) %>% 
  rename(kloc = Var1,
         tau = Var2) %>% 
  mutate(theta = utheta[as.character(tau)])

loc_vol_vec <- ssvi_local_vol(powerlaw_par, 
                              grid_df$tau, 
                              grid_df$theta, 
                              grid_df$kloc, 
                              "powerlaw")
# Matrix where k is row
loc_vol_m <- matrix(loc_vol_vec, nrow = length(kloc))  

# Plot Local Volatility Surface
# x is k, y is tau
M <- mesh(kloc, utau)
surf3D(M$x, M$y, loc_vol_m, colkey = FALSE, bty = "b2", 
       phi = 20, ticktype = "detailed",
       xlab = "k",
       ylab = "\u03c4",
       zlab = "Volatilidade local \u03c3")

Conclusão

Apresentamos o modelo de superfícies SVI, que faz uma generalização dos smiles SVI e apresenta vantagens sobre a parametrização fatia-a-fatia em virtude dos teoremas sobre arbitragem estática, apresentando restrições para os parâmetros do modelo que garantem a ausência deste tipo de arbitragem.

Uma vez parametrizada toda uma SSVI, torna-se simples, uma mera aplicação de fórmulas a obtenção tanto da densidade da distribuição neutra ao risco implícita no preços das opções, como da superfície de volatilidade local através da equação de Dupire.

Referências

Breeden, Douglas T, and Robert H Litzenberger. 1978. “Prices of State-Contingent Claims Implicit in Option Prices.” Journal of Business. JSTOR, 621–51.

Dupire, Bruno. 1994. “Pricing with a Smile.” Risk 7 (1): 18–20.

Gatheral, Jim. 2004. “A Parsimonious Arbitrage-Free Implied Volatility Parameterization with Application to the Valuation of Volatility Derivatives.” Presentation at Global Derivatives & Risk Management, Madrid.

Gatheral, Jim, and Antoine Jacquier. 2014. “Arbitrage-Free Svi Volatility Surfaces.” Quantitative Finance 14 (1). Taylor & Francis: 59–71.

  • Referimos o leitor ao material do prof. Antoine Jacquier para uma prova desta relação aqui

    Posted by Rafael F. Bressan in Derivativos & Riscos, 0 comments
    Calibrando um SVI

    Calibrando um SVI

    Neste post iremos mostrar como fazer uma calibração de um smile SVI baseado nos trabalhos de (Gatheral 2004) e (De Marco and Martini 2009). Escolheremos apenas uma fatia da superfície de volatilidade, fixando o tempo para expiração (maturidade) e coletando as volatilidades implícitas para diversos strikes.

    Como já apresentado em posts anteriores, existem diversas formas de interpolar, extrapolar, parametrizar e calibrar smiles de volatilidade. Exsitem vantagens e desvantagens para cada método. Neste post iremos fixar nossa atenção no modelo paramétrico de smile chamado SVI – Stochastic Volatility Inspired – uma forma que une a “simplicidade” de um modelo paramétrico com o poder de adesão aos dados de mercado dos modelos de volatilidade estocástica (i.e. Heston, SABR e afins).

    Chamar o modelo SVI de simples é puro eufemismo, ele é um modelo poderoso, com fundamento teórico avançado e diversos detalhes para sua calibração.

    Modelo SVI

    Este modelo foi apresentado por Jim Gatheral na conferência Global Derivatives & Risk Management 2004 e foi bem recebido pelos profissionais de mercado interessados em superfícies de volatilidade para equities, principalmente.

    O modelo possui duas propriedades que são as razões para sua popularidade. Ele satisfaz a fórmula do momento de Lee (2004), que é um resultado independente de modelo que especifica os limites assintóticos para um smile de volatilidade implícita. Portanto, o modelo SVI é válido para extrapolação além da região central dos dados disponíveis. Além disso, afirma-se que o modelo SVI é relativamente fácil de calibrar para dados de mercado, de modo que a superfície de volatilidade implícita correspondente é livre de arbitragem de calendário. As condições que garantem a ausência de arbitragem de borboleta forma resolvidas em um segundo artigo por Gatheral and Jacquier (2014).

    No SVI é possível se estabelecer condições explícitas em seus parâmetros, de modo que o modelo não gere preços onde oportunidades de arbitragem estáticas possam ocorrer. A calibração para dados reais de mercado requer algoritmos de otimização não-linear e pode ser bastante demorada. Mais recentemente, um método para calibração que usa a estrutura inerente do modelo para reduzir as dimensões do problema de otimização foi desenvolvido em De Marco and Martini (2009).

    A parametrização conhecida como RAW do SVI é apresentada na equação (1), seguindo a notação já introduzida anteriormente, portanto, estamos modelando a variância total implícita para um determinado prazo. Para diferentes maturidades, teremos diferentes conjuntos de parâmetros.

    \[\begin{equation} w(k) = a + b\left(\rho(k-m)+\sqrt{(k-m)^2 + \sigma^2}\right) \tag{1} \end{equation}\]

    onde: \(a \in \mathbb R\), \(b \geq 0\), \(|\rho| < 1\), \(m \in \mathbb R\), \(\sigma > 0\), e \(a+b \sigma\sqrt{1 − \rho^2} \geq 0\) para garantir que \(\min w(k)>0, \, \forall k \in \mathbb R\).

    Restrições de não-arbitragem

    Antes de demonstrar a restrição imposta aos parâmetros \(b\) e \(\rho\) em função dos limites de inclinação das asas do smile, vamos derivar as expressões para \(w\prime(k)\) e \(w\prime\prime(k)\) que nos serão úteis na demonstração.

    A expressão para \(w\prime(k)\) é bastante simples:

    \[\begin{equation} w\prime(k) = b \left[\rho + \frac{(k-m)}{\sqrt{(k-m)^2+\sigma^2}}\right] \tag{2} \end{equation}\]

    Derivando novamente a equação (2) em relação a \(k\) teremos uma expressão ainda mais simples, mesmo que após alguma manipulação algébrica um tanto tediosa1, e resulta em:

    \[\begin{equation} w\prime\prime(k)=\frac{b\sigma^2}{[(k-m)^2+\sigma^2]^{3/2}} \tag{3} \end{equation}\]

    onde, se considerarmos \(b>0\) temos que \(w\prime\prime(k)>0, \,\forall k\in \mathbb R\), ou seja, o smile de volatilidade definido pela equação (1) é estritamente convexo.

    Rogers and Tehranchi (2010) definiram os limites possíveis para a inclinação das asas em função do tempo para expiração, provando que o smile tende a ficar mais horizontal a medida que o prazo aumenta. Este limite pode ser escrito da seguinte forma e é uma condição necessária para a ausência de arbitragem:

    \[\begin{equation} |w\prime(k)|\leq \frac{4}{\tau} \qquad \forall k \in \mathbb R, \quad \forall \tau \in (0, \infty) \tag{4} \end{equation}\]

    Sendo o smile convexo, suas máximas inclinações ocorrem quando \(k\rightarrow \pm \infty\). Portanto, deve-se avaliar a restrição dada pela equação (4) nestes limites da seguinte maneira:

    \[\begin{align} \lim\limits_{k\rightarrow\infty}w\prime(k)&=b(1+\rho)\geq 0\\ \lim\limits_{k\rightarrow-\infty}w\prime(k)&=-b(1-\rho)\leq 0 \end{align}\]

    que satisfazendo estas duas relações ao mesmo tempo em que se restringe os parâmetros \(b\) e \(\rho\) através da inequalidade de Rogers e Tehranchi nos garante o seguinte resultado para um SVI livre de arbitragem de travas.

    \[\begin{equation} b(1+|\rho|)\leq\frac{4}{\tau} \tag{5} \end{equation}\]

    Para garantir que a superfície gerada está livre de arbitragem do tipo borboleta deve-se primeiramente definir uma função2 \(g: \mathbb R\rightarrow \mathbb R\), tal que:

    \[\begin{equation} g(k)=\left(1-\frac{kw\prime(k)}{2w(k)}\right)^2-\frac{w\prime(k)^2}{4}\left(\frac{1}{w(k)}+\frac{1}{4}\right)+\frac{w\prime\prime(k)}{2} \tag{6} \end{equation}\]

    e seguir o lema :

    Lema 1 Uma fatia da superfície de volatilidade está livre de arbitragem do tipo borboleta se, e somente se, \(g(k) \geq 0\) para todo \(k \in \mathbb R\) e \(\lim\limits_{k\rightarrow+\infty}d_1(k)=-\infty\).

    Infelizmente, a natureza altamente não-linear da função \(g(k)\) impossibilita a derivação de restrições gerais aos parâmetros do SVI. A forma mais simples de eliminar arbitragem do tipo borbobleta é incluir a restrição \(g(k) \geq 0\) na função perda e proceder com a calibração dos parâmetros.

    Reparametrização Quasi-explicit

    Um dos problemas mais marcantes com a calibração do SVI dado pela equação (1) é sua natureza altamente não-linear que gera inúmeros pontos de mínimo locais. Mesmo em um ambiente simulado, um típico otimizador de mínimos quadrados como Levenberg-Marquardt não consegue chegar ao mínimo global, onde a função perda é igual a zero. A solução encontrada é dependente dos valores iniciais inputados ao otimizador e a robustez do conjunto de parâmetros encontrados não é garantida.

    Uma forma de contornar este problema pode ser a utilização de otimizadores globais, como algortimos genéticos, em um primeiro estágio e então o refinamento desta solução através de um otimizador local (LM, por exemplo).

    Outra forma, adotada em De Marco and Martini (2009) é a reparametrização da equação (1) de forma que esta possa ser tratada como um prolema linear. Para tanto, considere a seguinte troca de variáveis:

    \[\begin{equation} y = \frac{k-m}{\sigma} \end{equation}\]

    então a parametrização RAW do SVI se torna.

    \[\begin{equation} w(y) = a + b\sigma\left(\rho y + \sqrt{y^2 + 1}\right) \end{equation}\]

    Definindo agora as seguintes variáveis reparametrizadas é possível reduzir a dimensão de parâmetros de um SVI de 5 para apenas 3:

    \[\begin{align} c = &b\sigma\\ d = &\rho b \sigma \end{align}\] \[\begin{equation} w(y)=a+dy+c\sqrt{y^2+1} \tag{7} \end{equation}\]

    Portanto, para um par fixo de \((m, \sigma)\) nosso problema reduzido é:

    \[\begin{equation} P_{m, \sigma}:=\min\limits_{a, c, d \in D}f_y(a, c, d) \tag{8} \end{equation}\]

    onde \(f_y(\cdot)\) é a função objetivo da reparametrização, e é dada pela seguinte equação:

    \[\begin{equation} f_y(a, c, d)=\sum_{i=1}^{n}\left[w(y_i)-\tilde w_i\right]^2 \tag{9} \end{equation}\]

    onde \(\tilde w_i\) é a variância total observada correspondente ao moneyness \(k_i\).

    O domínio \(D\) dos parâmetros \(\{a, c, d\}\) é encontrado a partir do limite imposto por (5).

    \[\begin{equation} D = \begin{cases} 0 \leq c \leq 4\sigma\\ |d| \leq c \quad \text{e}\quad |d| \leq 4\sigma – c\\ 0 \leq a \leq \max\{\tilde w_i\}\\ \end{cases} \tag{10} \end{equation}\]

    O problema reduzido dado pela equação (8), é um típico problema de mínimos quadrados com restrições lineares. Este problema, por ser convexo, admite uma única solução interior (se existente) que será o mínimo global para este problema e é encontrada através do gradiente da função objetivo igualando-o a zero, \(\nabla f_y = 0\). Esta equação gera um sistema linear nos parâmetros \(a, c, d\) que pode ser explicitamente resolvido. Caso a solução encontrada para este problema esteja contida no domínio \(D\), esta solução é interior e é o mínimo desejado, caso contrário, deve-se percorrer o perímetro do domínio e encontrar o menor valor da função objetivo que será uma solução de canto.

    Seja \((a^*, c^*, d^*)\) a solução de (8) e \((a^*, b^*, \rho^*)\) os correspondentes parâmetros originais recuperados, então o problema completo de calibração é:

    \[\begin{equation} P:=\min\limits_{m, \sigma}\sum_{i=1}^n (w_*(k_i)-\tilde w_i)^2 \tag{11} \end{equation}\]

    onde \(w_*(k)=a^*+b^*\left(\rho^*(k-m)+\sqrt{(k-m)^2 + \sigma^2}\right)\).

    O problema completo, (11) é um problema em apenas duas dimensões, \((m, \sigma)\) e não-linear, que deve ser abordado através de algum tipo de otimizador global.

    Solução explícita do problema reduzido

    Nesta seção apresentaremos a solução, sem considerar as restrições impostas em (10) para o problema reduzido em (8), algo omitido em De Marco and Martini (2009). Esta seção é opcional para o leitor atento que já percebeu a semelhança entre o problema reduzido e um típico problema de regressão linear múltipla.

    Para encontrar o conjunto de parâmetros \((a^*, c^*, d^*)\) que representam os valores ótimos na equação (8), devemos resolver o seguinte sistema de equações:

    \[\begin{equation} \nabla f_y = \left[ \begin{array}{c} \partial f_y / \partial a\\ \partial f_y / \partial d\\ \partial f_y / \partial c \end{array} \right] = \boldsymbol{0} \tag{12} \end{equation}\]

    Cada uma das derivadas parciais da equação acima quando igualadas a zero dão origem ao sistema linear apresentado abaixo:

    \[\begin{equation} \scriptsize \begin{bmatrix} &n &\sum y_i &\sum\sqrt{y_i^2+1}\\ &\sum y_i &\sum y_i^2 &\sum(y_i\sqrt{y_i^2+1})\\ &\sum\sqrt{y_i^2+1} &\sum(y_i\sqrt{y_i^2+1}) &\sum(y_i^2+1)\\ \end{bmatrix} \cdot \begin{bmatrix} a \\ d \\ c \end{bmatrix} = \begin{bmatrix} \sum\tilde w_i \\ \sum \tilde w_i y_i \\ \sum(\tilde w_i\sqrt{y_i^2+1}) \end{bmatrix} \tag{13} \end{equation}\]

    Portanto, o problema reduzido pode ser resolvido através de um sistema linear de ordem 3 sob restrições também lineares.

    Algoritmo

    A otimização para encontrar os parâmetros ótimos de um SVI dadas observações de mercado e a técnica de calibração Quasi-explicit de De Marco and Martini (2009) pode ser resumida nos seguintes passos:

    1. Definir valores iniciais para os parâmetros \((m, \sigma)\),

    2. Iniciar algum otimizador global com estes parâmetros e resolver o problema completo (11)

      2.1 Dentro da otimização global, resolver o problema reduzido (8) para os parâmetros \((m, \sigma)\) dados,

    3. Na convergência do problema completo do passo 2, otimizar uma última vez o problema reduzido,

    4. Recuperar os parâmetros \((a, b, \rho, m, \sigma)\)

    A escolha dos otimizadores fica a cargo pessoal, sendo sugerido testar vários para o mesmo problema. Eventualmente, para um determinado smile um otimizador pode se mostrar melhor que outro que vinha sendo utilizado em outras ocasiões.

    Nos testes realizados pelo Clube de Finanças, entre os otimizadores globais para o problema completo utilizamos Algoritmos Genéticos, Nelder-Mead restrito e um Método não-linear generalizado. Para o problema reduzido, apesar de ser linear, o método de Nelder-Mead restrito se mostrou tão eficiente quanto e de mais fácil implementação. Se o objetivo for fazer uma calibração direta, dos cinco parâmetros ao mesmo tempo, uma combinação de otimizador global em primeiro estágio e o método de Levenberg-Marquardt restrito para refinamento da solução é o ideal.

    Resultados

    A seguir apresentamos um smile de referência para a calibração, obtido de ivolatility.com e então partimos para diferentes técnicas de calibração de um SVI. Os códigos em R também estão apresentados ao longo do texto para melhor compreensão e estudo do leitor.

    Os dados utilizados neste exemplo estão apresentados na tabela 1 abaixo. Esta é uma típica apresentação de um slice de superfície, ou seja, dados para um smile apenas. As principais variáveis são: a data em que os dados foram coletados (date), o preço de fechamento do ativo (stock_price), o prazo para expiração em dias (period), e medidas de moneyness como delta, strike e o próprio moneyness, além é claro da volatilidade implícita (iv) retirada do mercado.

    Tabela 1: Dados reais para exemplo de calibração de uma SVI.
    date symbol exchange stock_price_for_iv period delta moneyness strike iv
    2017-09-21 IWM NYSEArca 143.73 30 10 0.03 148.41 0.09
    2017-09-21 IWM NYSEArca 143.73 30 15 0.03 147.49 0.09
    2017-09-21 IWM NYSEArca 143.73 30 20 0.02 146.80 0.09
    2017-09-21 IWM NYSEArca 143.73 30 25 0.02 146.21 0.09
    2017-09-21 IWM NYSEArca 143.73 30 30 0.01 145.69 0.09
    2017-09-21 IWM NYSEArca 143.73 30 35 0.01 145.19 0.10
    2017-09-21 IWM NYSEArca 143.73 30 40 0.01 144.69 0.10
    2017-09-21 IWM NYSEArca 143.73 30 45 0.00 144.18 0.10
    2017-09-21 IWM NYSEArca 143.73 30 50 0.00 143.66 0.10
    2017-09-21 IWM NYSEArca 143.73 30 55 0.00 143.12 0.11
    2017-09-21 IWM NYSEArca 143.73 30 60 -0.01 142.53 0.11
    2017-09-21 IWM NYSEArca 143.73 30 65 -0.01 141.88 0.11
    2017-09-21 IWM NYSEArca 143.73 30 70 -0.02 141.13 0.12
    2017-09-21 IWM NYSEArca 143.73 30 75 -0.02 140.26 0.13
    2017-09-21 IWM NYSEArca 143.73 30 80 -0.03 139.16 0.13
    2017-09-21 IWM NYSEArca 143.73 30 85 -0.04 137.66 0.14
    2017-09-21 IWM NYSEArca 143.73 30 90 -0.06 135.32 0.16

    Esta tabela poderia conter (de fato contém no arquivo original) outros períodos de expiração, e neste caso uma das colunas de moneyness começa a se repetir, no caso seria o delta pois baixamos uma tabela de volatilidades implícitas por delta. Assim, em uma tabela simples em formato tidy é possível armazenar informações de uma superfície inteira, a qual de outra forma necessitaria de um arranjo em 3 dimensões.

    Ressaltamos aqui que a unidade de volatilidade implícita está em percentuais ao ano, equanto que nosso período é de dias corridos. É necessário harmonizar estas medidas de forma que, para volatiliades dadas em percentual ao ano, o período também seja dado em anos. Logo nosso \(\tau = 30/365\), ou seja, 0.08219.

    Demonstraremos aqui os resultados para a calibração de uma RAW SVI pelos métodos “Direto”, “GA”, “Quasi-NM” e “Quasi-PQ”, abaixo explicados.

    O método “Direto” é uma calibração direta através de um algoritmo de Levenberg-Marquardt da equação (1), ou seja, não existe reparametrização Quasi-explicit e o problema resolvido é não-linear em 5 dimensões. São realizadas 10 calibrações com estimativas iniciais dos parâmetros aleatórias, mas dentro de seus respectivos domínios. A melhor solução, aquela com o menor valor para a função objetivo, é selecionada. Todos os outros métodos utilizam a reparametrização, ocorrendo variações apenas nos algoritmos de otimização utilizados nos problemas reduzido e completo.

    A calibração “GA” faz uso do otimizador global de algoritmos genéticos para o problema completo, ou seja, para estimar o par \((m, \sigma)\) que corresponde ao mínimo global. Após, o problema reduzido é resolvido através do algoritmo de Nelder-Mead. Este método é robusto, pois o algoritmo genético tem grande probabilidade de encontrar a região onde se encontra o mínimo global e não ficar preso localmente. Entretanto a robustez ocorre as expensas do tempo de computação.

    Os métodos ditos “Quasi” diferem entre si na resolução do problema reduzido. Enquanto “PQ” remete a programação quadrática e faz uso da resolução do sistema linear apresentado na equação (13) com as restrições impostas por (10), o método “Quasi-NM” utiliza o método de Nelder-Mead com restrições para a resolução deste mesmo problema reduzido. Em ambos os métodos, o problema completo é resolvido com um algoritmo de Nelder-Mead com 50 reinicializações das estimativas iniciais dos parâmetros \((m, \sigma)\), o que causa algum impacto no tempo de computação destes métodos.

    smile <- dados %>% 
      mutate(tau = period / 365) %>% 
      select(moneyness, iv, tau)
    
    par_names <- factor(c("a", "b", "$\\rho$", "m", "$\\sigma$"),
                        levels = c("a", "b", "$\\rho$", "m", "$\\sigma$"))
    k <- smile$moneyness
    w <- smile$iv^2 * smile$tau
    
    init_direct <- proc.time()
    par_direct <- svi_fit_direct(k, w)
    end_direct <- proc.time()
    time_direct <- end_direct - init_direct 
    
    init_ga <- proc.time()
    par_ga <- svi_fit_ga(k, w)
    end_ga <- proc.time()
    time_ga <- end_ga - init_ga 
    
    init_quasipq <- proc.time()
    par_quasipq <- svi_fit_quasi(k, w, inner = "quadprog")
    end_quasipq <- proc.time()
    time_quasipq <- end_quasipq - init_quasipq 
    
    init_quasinm <- proc.time()
    par_quasinm <- svi_fit_quasi(k, w)
    end_quasinm <- proc.time()
    time_quasinm <- end_quasinm - init_quasinm 
    
    iv_direct <- sqrt(svi_fun(par_direct$par[[1]], k) / smile$tau)
    iv_ga <- sqrt(svi_fun(par_ga$par[[1]], k) / smile$tau)
    iv_quasipq <- sqrt(svi_fun(par_quasipq$par[[1]], k) / smile$tau)
    iv_quasinm <- sqrt(svi_fun(par_quasinm$par[[1]], k) / smile$tau)
    
    plot_tbl <- tibble(k = k,
                  Direct = iv_direct,
                  GA = iv_ga,
                  QuasiPQ = iv_quasipq,
                  QuasiNM = iv_quasinm,
                  observed = smile$iv) %>% 
      gather(key = method, value = iv, -c(k, observed))
    
    par_tbl <- bind_rows(par_direct, par_ga, par_quasipq, par_quasinm) %>% 
      select(method, par) %>% 
      unnest() %>% 
      mutate(names = rep(par_names, 4)) %>% 
      spread(method, par) %>% 
      select(names, Direct, GA, QuasiPQ, QuasiNM) %>% 
      mutate(names = as.character(names)) %>% 
      mutate_at(vars(Direct:QuasiNM), arred)
    
    rmse_tbl <- bind_rows(par_direct, par_ga, par_quasipq, par_quasinm) %>% 
      select(method, par) %>% 
      unnest() %>% 
      group_by(method) %>% 
      summarise(RMSE = rmse(par, k, w)) %>% 
      spread(method, RMSE) %>% 
      mutate(names = "RMSE") %>% 
      select(names, Direct, GA, QuasiPQ, QuasiNM) %>% 
      mutate_at(vars(Direct:QuasiNM), format, digits = 3, scientific = TRUE)
    
    time_tbl <- tibble(method = c("Direct", "GA", "QuasiPQ", "QuasiNM"),
                       time = rbind(time_direct, time_ga, 
                                    time_quasipq, time_quasinm)[, 3]) %>% 
      spread(method, time) %>% 
      mutate(names = "Tempo") %>% 
      select(names, Direct, GA, QuasiPQ, QuasiNM) %>% 
      mutate_at(vars(Direct:QuasiNM), arred)
    
    frame_tbl <- bind_rows(par_tbl, rmse_tbl, time_tbl)

    Abaixo é apresetanda uma tabela com os valores estimados para os parâmetros da SVI, o RMSE (root mean square error) e o tempo total em segundos para a calibração. Aqui o RMSE é definido como \(\sqrt{1/n\sum(w(k_i)-\tilde w_i)^2}\) e nos fornece um valor típico de erro na variância.

    kable(frame_tbl,
          col.names = c("Estimativa", "Direto", "GA", 
                        "QuasiPQ", "QuasiNM"),
          caption = "Parâmetros estimados da calibração, RMSE e tempo de computação em segundos.") %>% 
      kable_styling(bootstrap_options = "striped",
                    font_size = 20,
                    full_width = FALSE)
    Tabela 2: Parâmetros estimados da calibração, RMSE e tempo de computação em segundos.
    Estimativa Direto GA QuasiPQ QuasiNM
    a 0.00000 0.00000 0.00000 0.00011
    b 0.01964 0.01952 0.01651 0.01870
    \(\rho\) -0.81157 -0.80220 -1.00000 -0.90090
    m -0.00861 -0.00773 -0.01131 -0.01145
    \(\sigma\) 0.05101 0.05039 0.07100 0.05027
    RMSE 8.71e-06 8.68e-06 9.62e-05 1.02e-05
    Tempo 0.11800 24.93800 0.14100 5.94600

    O método Direto, com algoritmo de Levenberg-Marquardt se mostrou muito mais rápido que os demais, principalmente com relação ao algoritmo genético, e com um bom ajuste dado o baixo valor de RMSE. O algoritmo genético é consideravelmente mais lento, entretanto durante as várias calibrações realizadas em testes (e que não estão apresentadas na tabela 2), este algoritmo sempre se mostrou robusto, com baixo RMSE, diferentemente dos outros métodos que por vezes, denpendendo das estimativas iniciais, podem convergir para um mínimo local.

    O gráfico com os ajustes realizados pode ser observado abaixo.

    ggplot(plot_tbl, aes(x = k)) + 
      geom_point(aes(y = observed)) +
      geom_line(aes(y = iv, color = method)) +
      guides(color = guide_legend(title = "")) +
      labs(title = "",
           x = "Forward log-moneyness",
           y = "Volatility",
           caption = "") +
      scale_y_continuous(labels = scales::percent) +
      scale_color_viridis_d() +
      theme_economist_white()
    Comparação entre diferentes métodos de calibração de uma SVI.

    Figura 1: Comparação entre diferentes métodos de calibração de uma SVI.

    De fato o método direto e o método Quasi-explicit com otimizador global do tipo algoritmo genético se mostram mais adequados para a calibração de um SVI. Enquanto o método direto é muito mais eficiente em termos computacionais, o método Quasi-explicit com GA é mais robusto. Desta forma, deve-se salientar que é necessário que o usuário, ao fazer uma calibração de smile de volatilidade, deve dispor de diferentes métodos de fazê-lo, e a inspeção visual do resultado é obrigatória para determinar qual método foi mais eficiente em ajustar a curva aos dados.

    Conclusão

    Apesar de neste exemplo ter se mostrado um método efetivo, com bom ajuste e baixo tempo de calibração, o método direto é altamente dependente dos valores iniciais dos parâmetros. Para tornar este método mais robusto, um número maior de reinicializações deve ser feita o que penaliza o tempo de calibração. O método Quasi-explicit com algoritmo genético para encontrar a região de \((m, \sigma)\) onde se encontra o mínimo global se mostrou bastante robusta, entretanto, de convergência lenta. Para ajustar apenas um smile alguns segundos a mais não representam problema. Porém, se imaginarmos que em uma grande instituição financeira são necessárias calibrações de, talvez, milhares de smiles representando inúmeras superfícies de diversos instrumentos, este método pode se mostrar computacionalmente caro.

    Já os métodos Quasi-explicit que utilizam um algoritmo de Nelder-Mead para a resolução do problema completo se mostraram muito sensíveis às estimativas iniciais dos parâmetros. Mesmo utilizando 50 reinicialzações do método, diversas vezes o ajuste realizado foi insatisfatório. A resolução através de programação quadrática é rápida, se comparada com NM, entretanto, quando as restrições impostas pela equação (10) se tornam ativas, este método parece sofrer com algum viés em sua solução.

    Referências

    De Marco, S, and C Martini. 2009. “Quasi-Explicit Calibration of Gatheral’s Svi Model’.” Zeliade White Paper, 1–15.

    Gatheral, Jim. 2004. “A Parsimonious Arbitrage-Free Implied Volatility Parameterization with Application to the Valuation of Volatility Derivatives.” Presentation at Global Derivatives & Risk Management, Madrid.

    Gatheral, Jim, and Antoine Jacquier. 2014. “Arbitrage-Free Svi Volatility Surfaces.” Quantitative Finance 14 (1). Taylor & Francis: 59–71.

    Lee, Roger W. 2004. “The Moment Formula for Implied Volatility at Extreme Strikes.” Mathematical Finance: An International Journal of Mathematics, Statistics and Financial Economics 14 (3). Wiley Online Library: 469–80.

    Rogers, Leonard CG, and MR Tehranchi. 2010. “Can the Implied Volatility Surface Move by Parallel Shifts?” Finance and Stochastics 14 (2). Springer: 235–48.

  • A resolução desta derivada é uma simples regra da divisão, entretanto a simplificação do resultado pede alguma manipulação algébrica. É possível utilizar sistemas de computação simbólica, o qual recomendamos o SymPy

  • Condições para ausência de arbitragem do tipo borboleta em um SVI estão detalhadas na seção 2.2 do artigo de Gatheral and Jacquier (2014).

    Posted by Rafael F. Bressan in Derivativos & Riscos, 0 comments
    Métodos de calibração de superfícies de volatilidade

    Métodos de calibração de superfícies de volatilidade

    2019/02/08

    Métodos de calibração são as diferentes formas existentes entre “interpolação”, “suavização” e “parametrização” que podem ser utilizadas para fazer o ajustes dos dados obtidos do mercado de opções às superfícies de volatilidade.

    Como já apresentado em artigos anteriores, existem diversas formas de interpolar, extrapolar, parametrizar e calibrar smiles de volatilidade. Existem vantagens e desvantagens para cada método.

    Calibração versus Interpolação

    Uma forma simples de gerar um smile de volatilidade a partir de dados observados no mercado é a interpolação destes dados. Diversas formas de interpolação existem, sendo talvez a mais conhecida a spline cúbica. Não é a proposta deste artigo detalhar os procedimentos de interpolação, restando saber que em tal procedimento é gerada uma função contínua em partes (piecewise) que passa por todos os pontos observados.

    Uma interpolação força a passagem da função interpolada em todos os seus pontos de referência, como se estivesse ligando estes pontos em um desenho a mão livre. Portanto, nestes pontos o erro da interpolação é zero por definição, entretanto em pontos intermediários podem surgir erros, inclusive aqueles que possibilitam a existência de arbitragem entre strikes de um mesmo smile1.

    Em contraposição a métodos de interpolação, podem ser derivados métodos de suavização (smoothing) ou então a parametrização do smile de volatilidade. Seja suavização, ou parametrização, estes métodos não forçam a passagem da função que representa o smile pelos pontos de mercado, mas buscam minimizar alguma função perda dos desvios em relação a estes pontos ao mesmo tempo em que buscam “suavizar” o smile, para que este não apresente variações bruscas entre os strikes ou alterações de convexidade na curva de preços, que não são condizentes com a teoria de precificação de derivativos.

    Um método paramétrico, como o SVI, Heston ou SABR, busca ajustar às volatilidades implícitas observadas através dos preços das opções sendo praticados no mercado a uma determinada função, que possui parâmetros em sua definição que por sua vez determinam a forma desta. Ao se ajustar os parâmetros, pode-se adequar a função para ficar “o mais próxima possível” dos dados observados, sem necessariamente, no entanto, passar por todos estes pontos.

    A figura abaixo tenta mostrar as diferenças entre uma interpolação spline cúbica, uma suavização e uma parametrização SVI. Enquanto que a interpolação liga todos os pontos marcados, a suavização e a parametrização não necessariamente passam sobre estes mas fornecem uma curva mais “suave”, sem trocas de convexidade, o que geraria oportunidades de arbitragem e probabilidades negativas de ocorrência de determinados preços para o ativo subjacente, que ferem os princípios de precificação de opções. Os dados utilizados neste e nos próximos artigos sobre superfícies de volatilidade foram obtidos do site ivolatility.com na forma de amostra gratuita fornecida livremente. O ativo subjacente é o ETF IWM para a data de 21/09/2017.

    Diferentes métodos de ajuste de dados a um smile.

    Figura 1: Diferentes métodos de ajuste de dados a um smile.

    Pode-se verificar como os métodos SVI e a suavização não passam sobre todos os pontos marcados, com a suavização tendo dificuldade com a curvatura nos valores mais altos de moneyness e a SVI possuindo uma inclinação mais branda na asa esquerda do smile.

    Spline cúbica

    Este método é possivelmente um dos mais flexíveis e conhecidos de interpolação de dados univariados existente, embora também exista sua versão bi-dimensional. Uma spline nada mais é que “uma curva definida matematicamente por dois ou mais pontos de controle”2.

    No caso da spline cúbica, esta é uma função polinomial de grau 3 definida em cada subintervalo demarcados pelos pontos de controle, no caso de interpolação são todos nós. Ou seja, considere um segmento entre dois pontos consecutivos \([c, d]\in S\) a spline é uma função cúbica com seus parâmetros calculados pelo algoritmo de ajuste. Para o próximo intervalo de pontos dentro do domínio da função, um novo polinômio de grau 3 é ajustado, sendo que nos pontos de nós uma restrição de igualdade entre as derivadas nos dois segmentos é aplicada para garantir a suavidade da função interpolada como um todo.

    Assim, uma spline cúbica é uma função contínua, suave e diferenciável até a segunda ordem. Entretanto, suas derivadas, apesar de contínuas, podem não ser suaves, especialmente aquela de segunda ordem que pode apresentar pontos de “ruptura”. Esta característica de uma spline cúbica a torna pouco atrativa para a inferência de distribuições de probabilidade a partir de dados de volatilidade ou mesmo dos preços de opções.

    Cada segmento de uma spline cúbica é um polinômio de grau 3 diferente.

    Figura 2: Cada segmento de uma spline cúbica é um polinômio de grau 3 diferente.

    Suavização

    A técnica de suavização é muito semelhante a interpolação, inclusive o método spline também é aplicado, com algumas modificações de forma que nem todos os pontos fornecidos serão nós.

    Na spline de suavização (ou aproximação), os pontos fornecidos são separados entre os nós, onde a função deve passar e pontos de controle, que são utilizados para controlar a curvatura da função nestes pontos.

    Estas suavizações são principalmente utilizadas quando se possui muitas observações sujeitas a ruídos, de forma que uma interpolação entre todos os pontos seria tanto impraticável quanto sem sentido. O que se deseja, portanto, é uma função aproximada que melhor descreva o processo sob análise.

    Um ponto em comum entre estas técnicas é o parâmetro de suavização, ausente, na interpolação, que controla a “suavidade” da função estimada.

    Menor parâmetro de suavização gera granularidade na curva.

    Figura 3: Menor parâmetro de suavização gera granularidade na curva.

    Parametrização

    E por fim as técnicas de parametrização. Nesta categoria estão diversos conhecidos modelos de superfícies de volatilidade implícita, dentre eles os modelos de Heston (1993) e SVI de Gatheral (2004).

    Em comum, estes modelos tentam parametrizar a superfície, e por conseguinte o smile de volatilidade, de acordo com alguma função, em geral não-linear, que possui características condizentes com a teoria de precificão de derivativos e também a observação empírica das superfícies.

    Por exemplo, a parametrização raw da SVI possui a seguinte forma para a variância total3 :

    \[ w(k) = a + b\left(\rho(k-m)+\sqrt{(k-m)^2 + \sigma^2}\right)\]

    que fornece um espaço de cinco parâmetros \(\chi_B=\{a, b, \rho, m, \sigma\}\) que definem o smile e devem, portanto, serem calibrados a partir de dados observados no mercado.

    O procedimento de calibração consiste em encontrar o conjunto de parâmetros que minimizam uma função perda entre a volatilidade prevista pelo modelo e os dados de mercado, enquanto satisfazem algumas restrições adicionais, como “ausência de arbitragem”, suavidade, etc. Trata-se, via de regra, de problemas de otimização não-linear com restrições de inequalidade também não-lineares.

    Função perda

    A função perda, ou função de calibração pode ser definida de diversas maneiras, de forma geral, para uma determinada maturidade, ela toma a forma:

    \[L=\sum\limits_{i=1}^n\lambda_i||\hat w(k_i)-w_{imp}(k_i)||\] onde \(||\cdot||\) é alguma medida de norma, sendo a mais conhecida o quadrado das diferenças, dando origem a minimização do erro quadrático médio (MSE). Para o presente smile sendo calibrado existem \(n\) strikes (\(k_i\)) e suas volatilidades implícitas observadas são \(w_{imp}(k_i)\). A resposta do modelo para um determinado strike é \(\hat w(k_i)\) e \(\lambda_i\) são os pesos dados na função perda para cada um destes strikes.

    Os pesos \(\lambda_i\) são utilizados para ponderar as observações das volatilidades mais importantes para o cálculo, onde se deseja que a curva ajustada possua um menor erro. Em geral, estes pesos são calculado como inversamente proporcionais:

    • ao quadrado dos spreads bid-ask, para dar mais importância às opções mais líquidas
    • ao quadrado da grega vega calculada a partir do modelo BSM

    Otimizadores

    Os otimizadores são os algoritmos pelos quais o problema de minimização posto é resolvido. Se a função perda é convexa, e ela deve ser construída de forma a ser, mesmo que não estritamente, então ela possui um ou mais pontos de mínimo onde o gradiente desta função é igual a zero. O que os otimizadores fazem é buscar o conjunto de parâmetros que minimizam a função perda e atendem as restrições impostas simultaneamente. Os otimizadores podem ser classificados em dois grandes grupos, globais e locais.

    Algoritmos locais dependem de uma estimativa inicial dos parâmetros para começarem a busca pelo mínimo. Seguindo uma regra utilizando-se o gradiente da função ou alguma heurística, estes otimizadores caminham em direção ao ponto de mínimo mais próximo da estimativa inicial, daí o nome “local”. Como desvantagem destes otimizadores é a mais evidente é que se a função perda for altamente não-linear, com diversos pontos de mínimo local, este otimizador pode ficar preso em um destes pontos sem nunca, no entanto, encontrar o mínimo global. Eles são, portanto muito sensíveis à estimativa inicial dos parâmetros.

    Por sua vez, otimizadores globais buscam mapear todo o espaço factível para os parâmetros e encontrar o ponto mínimo da função perda dentro deste espaço. Estes algoritmos não dependem de estimativas iniciais, uma vez que tentarão avaliar o espaço completo. São utilizados quando o problema de minimização é não-linear e possui múltiplos pontos de mínimo local. Estes algoritmos usam alguma forma de heurística para encontrar a região onde o mínimo global está localizado, mas são, em geral, ineficientes em apontar rapidamente onde este ponto de mínimo se encontra com precisão. Por esta razão, é frequente a utilização de otimizadores globais com um posterior refinamento de sua solução por algum algoritmo local.

    Abaixo apresentamos alguns exemplos mais comuns de otimizadores, tanto locais quanto globais:

    • Gauss-Newton: Este método é utilizado para encontrar as raízes de alguma função. Para encontrar o ponto de mínimo da função perda, precisa-se encontrar as raízes do gradiente desta função, portanto o método de Newton em otimização faz uso da função gradiente. Este é um método de otimização local.

    • Levenberg-Marquardt: Método muito utilizado para problemas não-lineares, ele parte de uma modificação ao método de Gauss-Newton ao introduzir um fator de amortecimento calculado iterativamente.

    • L-BFGS-B: BFGS é um método conhecido como quasi-Newton, onde não é necessário calcular a Hessiana do problema, ela é aproximada a partir do próprio gradiente. É bastante utilizado para resolver problemas não-lineares e em sua versão L-BFGS-B pode lidar com restrições do tipo box, intervalo dos parâmetros é fixo.

    • Nelder-Mead: Este é um método livre do uso de gradiente, já que usa uma heurística para construir um simplex e a partir deste “procurar” por um mínimo. Bastante utilizado quando a função objetivo pode não ser diferenciável. Faz uso de um simplex inicial, que pode ser grande o suficiente para encampar o mínimo global, entretanto, não se classifica como um otimizador global.

    • Algoritmo Genético: Este método utiliza conceitos da seleção natural para gerar os resultados da otimização. É um otimizador global, no sentido que independe de uma estimativa inicial de parâmetros e faz uma busca por todo o espaço factível. Em um algoritmo genético, uma população aleatória inicial de parâmetros é criada e a partir desta, as gerações evoluem conforme mutações e cross-over de características e é avaliado o fitness de cada conjunto de parâmetros até um deles ser considerado adequado.

    • Evolução Diferencial: É um método de otimização global, assim como o Algoritmo Genético e o Enxame de Partículas. Sua diferença reside no fato de que sua população inicial é constantemente avaliada e deslocada de posição. Se o agente obtiver uma situação melhor (menor valor para a função perda) na nova posição, esta agora faz parte da população. Desta forma os agentes, antes espalhados pelo espaço factível dos parâmetros, tendem a convergir para um ponto com o menor valor da função perda.

    • Enxame de Partículas: Do inglês, Particle Swarm Optimization – PSO este método é semelhante ao DE (Differential Evolution) porém as partículas (o equivalente dos agentes no DE) matém informações sobre a posição da melhor partícula até então, de forma a fazer com que as partículas tendam para a melhor solução.

    Conclusão

    Dependendo do objetivo da aplicação, superfícies de volatilidade podem ser interpoladas, suavizadas ou parametrizadas. A parametrização tem recebido especial interesse pois pode, ao mesmo tempo que garante uma superfície livre de arbitragem estática se devidamente construída, ajustar-se muito bem aos dados observados e gerar distribuições neutras ao risco implícitas factíveis.

    Para gerar uma superfície parametrizada, primeiramente é necessário um modelo teórico com propriedades desejáveis e que se ajuste aos dados de mercado quando calibrado. Escolhido este modelo paramétrico, passa-se a calibração do mesmo onde existem diversas opções de escolha entre otimizadores. Ao final do processo teremos um modelo de superfície devidamente parametrizado com valores que melhor se ajustam segundo alguma função perda escolhida.

    Com a superfície de volatilidade calibrada, as aplicações possíveis incluem a precificação de derivativos, gerenciamento de risco, simulações de Monte Carlo, análises de stress, entre outras.

    Referências

    Gatheral, Jim. 2004. “A Parsimonious Arbitrage-Free Implied Volatility Parameterization with Application to the Valuation of Volatility Derivatives.” Presentation at Global Derivatives & Risk Management, Madrid.

    Heston, Steven L. 1993. “A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options.” The Review of Financial Studies 6 (2). Oxford University Press: 327–43.

  • Veja mais detalhes no artigo anterior: Smile de volatilidade – parte 2
  • Definição retirada de https://pt.wikipedia.org/wiki/Spline
  • A variância total é o produto entre a variância implícita e o tempo para expiração, (w=\sigma^2_{imp}\cdot\tau).
  • Posted by Rafael F. Bressan in Derivativos & Riscos, 0 comments