R

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
    Simulação de Monte Carlo

    Simulação de Monte Carlo

    Em artigo anterior, sobre processos estocásticos, fizemos uso de uma poderosa ferramenta computacional, frequentemente utilizada em finanças para fazer simulações. Naquele artigo simulamos cinco realizações de caminhos de um processo estocástico, cada um com 500 passos a frente. Esta técnica é conhecida como Simulação de Monte Carlo – SMC e será abordada no presente artigo.

    Neste artigo também iremos introduzir, no corpo do texto, os códigos em linguagem R utilizados para fazer as simulações, aumentando a didática de nossos artigos. A linguagem R é uma das preferidas para a modelagem estatística, é uma das linguagens de ciência de dados que vem ganhando muitos adeptos e, por conseguinte, é amplamente utilizada pelo mercado financeiro. E claro, é uma das preferidas aqui do CF também.

    Nosso problema será simular a posição de um portfólio composto por uma posição comprada em uma ação PETR4 e uma put PETRV17. A opção de venda (put) tem preço de exercício em R$ 16,92, data de expiração em 15/10/2018 e é do tipo europeia. Ao final do pregão do dia 21/09/2018 a ação PETR4 fechou cotada a R$ 20,14 e a put em R$ 0,12. A partir desta data até o dia da expiração da opção haverão 16 dias de pregão, que será nosso horizonte de simulação.

    Para melhoria da didática do texto e também para simplificação do problema, manteremos algumas variáveis necessárias para a precificação de opções constantes ao longo do período de análise, são elas:

    • Volatilidade: Será calculada a volatilidade implícita da opção da data de compra do portfólio, 21/09/2018, e será mantida constante a partir daí para fins de precificação na SMC;

    • Taxa de juros: constante no valor de 6,5 %a.a. tem termos contínuos;

    • Taxa de dividendos: suposto igual a zero.

    Simulação de Monte Carlo

    Para realizar uma SMC de um ativo financeiro deve-se primeiramente estabelecer uma distribuição de probabilidades que os retornos deste ativo deve seguir. Em nosso exemplo, utilizaremos a distrbuição normal para os retornos logarítimicos da ação PETR4, em linha com o clássico modelo Black & Scholes, certamente existem diversas variantes que se ajustam melhor a realidade dos mercados, entretanto este é o modelo mais conhecido e base de todos os demais.

    Uma vez escolhida a distribuição dos (log) retornos, tem-se de escolher valores para a média e variância desta distribuição. A média dos retornos iremos tirar do histórico da ação, o retorno médio diário dos último ano. A variância da distribuição será encontrada a partir da volatilidade implícita da opção na data de compra do portfólio. A função utilizada para encontrar esta volatilidade retorna um valor em termos anuais, portanto, conforme visto no artigo sobre processos estocásticos devemos reescalar uma volatilidade anual para diária, e isto é obtido fazendo a divisão por (\sqrt{252}), onde temos 252 dias úteis em 1 ano.

    Desta forma é possível fazer a simulação dos log-retornos da ação para cada um dos dias a frente, até a data de exercício da opção, 15/10/2018. Estes retornos são acumulados e o preço simulado da ação PETR4 em uma data intermediária é o preço de compra vezes o retorno acumulado até então.

    Faremos 1.000 simulações destas, gerando caminhos possíveis de preços para a ação. É necessário fazer muitas simulações para termos uma boa ideia da distribuição dos preços na data final, não é raro serem feitas mais de mil simulações, as vezes até dez mil podem ser necessárias.

    Uma vez gerados todos os caminhos simulados do preço do ativo objeto, podemos então precificar a put com base nestes preços simulados e as outras variáveis necessárias para se precificar uma opção europeia. Assim teremos também todos os caminhos de preço para a opção até sua data de exercício.

    O valor de nosso portfólio, em qualquer ponto do intervalo de tempo em análise, será a soma do preço da ação com o preço da opção e será possível verificar o efeito de proteção contra quedas do preço do ativo objeto a compra da put tem no portfólio.

    Cabe ressaltar aqui que o preço da opção não é simulado, não diretamente. Como a opção é um instrumento derivativo o seu preço “deriva” do preço do ativo objeto, este sim que é simulado. Uma vez que tenhamos o preço da ação, dadas nossas premissas de precificação, podemos calcular o prêmio da opção com base no modelo Black & Scholes.

    Implementação em R

    Conforme comentado, utilizamos aqui no CF a linguagem de programação R para realizar nossas atividades que envolvam métodos quantitativos em finanças. Abaixo irei apresentar o código utilizado, trecho a trecho e o resultado obtido ao final.

    Primeiramente, no R, devemos carregar em nossa sessão de trabalho os pacotes que serão utilizados ao longo do código. Os pacotes funcionam como extensões ao R base, nestes pacotes encontramos diversas funções já programadas por outras pessoas que facilitam (e muito!) a nossa codificação.

    library(tidyverse)
    library(ggthemes)
    library(tidyquant)
    library(RQuantLib)
    

    O pacote RQuantLib, por exemplo, possui já implementado dentro dele funções para fazer a precificação de opções europeias, sem que se tenha que implementar o modelo manualmente. Como a intenção deste artigo não é explicar o modelo Black & Scholes, vamos abstrair esta parte e simplesmente chamar uma função que nos retorna o valor da opção dadas as variáveis de entrada.

    Em seguida iremos definir algumas de nossas variáveis, como o ticker da ação para buscar seus dados históricos através da função tq_get() do pacote tidyquant e calcular os retornos logarítimicos e tirar sua média, o preço e data de exercício da opção e também iremos relacionar os dias de negócio entre a data de compra e vencimento.

    acao <- "PETR4.SA"
    p_exer <- 16.92
    d_exer <- as.Date("2018-10-15")
    d_atual <- as.Date("2018-09-21")
    dias <- seq(d_atual, d_exer, by = 1)
    dias <- dias[isBusinessDay("Brazil", dias)]
    nsims <- 1000
    ndias <- length(dias) - 1
    sim_nomes <- paste0("sim", 1:nsims)
    
    # Carregar os precos historicos da acao
    p_hist <- tq_get(acao, from = d_atual - years(1), to = d_atual + days(1)) %>% 
      filter(volume != 0.0)
    ret_hist <- p_hist %>% 
      tq_mutate(select = adjusted,
                mutate_fun = periodReturn,
                period = "daily",
                type = "log",
                leading = FALSE,
                col_rename = "log_ret") %>% 
      na.omit()
    rf <- log(1 + 0.065)
    div <- 0
    S0 <- last(ret_hist$adjusted)
    P0 <- 0.12
    mi <- 252 * mean(ret_hist$log_ret) # retorno medio em termos anuais
    sigma <- EuropeanOptionImpliedVolatility("put", P0, S0, p_exer, div, rf, 
                                             (ndias + 1) / 252, 0.30)
    

    Com o código acima obtemos basicamente todos os dados com os quais poderemos implementar a simulação de Monte Carlo. Entretanto, para realizar as simulações, necessitamos especificar mais algumas funções customizadas para nossas necessidades.

    Primeiro iremos especificar uma função que retorna uma única simulação de log-retornos acumulados em uma coluna de dados, esta função é chamada de mc_sim_fun. A segunda função necessária é a função de precificação da opção europeia. Por padrão, a função do pacote RQuantLib EuropeanOption() retorna uma lista com o valor da opção, mas também todas as suas gregas. Também de forma um tanto quanto estranha, esta função retorna o valor zero na data de exercício, mesmo que a opção esteja In-The-Money, portanto é necessário modificar este comportamento.

    # Funcao para realizar uma simulacao
    mc_sim_fun <- function(valor_i, N, media, volat){
      med_d <- media / 252
      volat_d <- volat / sqrt(252)
      ans <- tibble(c(valor_i, rnorm(N, med_d - (volat_d^2 / 2), volat_d))) %>% 
        `colnames<-`("log_ret") %>%
        mutate(ret_ac = cumsum(log_ret)) %>% 
        select(ret_ac)
    
      return(ans)
    }
    
    # Funcao para precificar uma opcao europeia
    eur_option <- function(type, underlying, strike, dividendYield, riskFreeRate, 
                           maturity, volatility) {
      if (maturity == 0.0) {
        ans <- switch(type,
                      put = max(strike - underlying, 0),
                      call = max(underlying - strike, 0))
        return(ans)
      }
    
      ans <- EuropeanOption(type = type,
                            underlying = underlying,
                            strike = strike,
                            dividendYield = dividendYield,
                            riskFreeRate = riskFreeRate,
                            maturity = maturity,
                            volatility = volatility)$value
      return(ans)
    }
    

    Uma vez com os dados obtidos e as funções auxiliares programadas, podemos passar a SMC propriamente dita. Aqui vamos estabelecer o número de simulações (1.000), calcular um data frame com os log-retornos acumulados e então calcular o preço da ação para cada dia e simulação realizados. O preço da ação na data (t) será (S_t=S_0 e^{r_t}), onde (r_t) é o log-retorno acumulado até a data (t).

    Após termos todos os preços do ativo objeto, passamos a computar qual seria o preço da opção, (P_t), naquelas condições. O valor do portfólio é dado pela soma destes dois preços (lembre-se, nosso portfólio é composto por uma ação e uma opção de venda).

    # Simulacao de Monte Carlo
    # Valores Iniciais
    inic <- rep(0, nsims) 
    set.seed(12345)
    ret_ac_mc <- map_dfc(inic,
                         mc_sim_fun,
                         N = ndias,
                         media = mi,
                         volat = sigma)
    
    precos_acao <- (S0 * exp(ret_ac_mc)) %>% 
      set_names(sim_nomes) %>% 
      mutate(anos_exp = (ndias:0) / 252) %>% 
      gather(key = sims, value = St, -anos_exp)
    
    # Evolucao do Portfolio
    port_mc <- precos_acao %>% 
      mutate(Pt = map2_dbl(St, anos_exp, 
                           ~eur_option(type = "put",
                                       underlying = .x,
                                       strike = p_exer,
                                       dividendYield = div,
                                       riskFreeRate = rf,
                                       maturity = .y,
                                       volatility = sigma)),
             port_valor = Pt + St,
             data = rep(dias, nsims))
    head(port_mc)
    
    ##     anos_exp sims       St         Pt port_valor       data
    ## 1 0.05952381 sim1 20.14000 0.10666291   20.24666 2018-09-21
    ## 2 0.05555556 sim1 20.56213 0.06610065   20.62823 2018-09-24
    ## 3 0.05158730 sim1 21.08354 0.03486500   21.11841 2018-09-25
    ## 4 0.04761905 sim1 21.01296 0.03001788   21.04297 2018-09-26
    ## 5 0.04365079 sim1 20.69410 0.03266826   20.72677 2018-09-27
    ## 6 0.03968254 sim1 21.14278 0.01516756   21.15794 2018-09-28
    

    O data frame port_mc contém todas as informações da SMC de nosso portfólio. Contém as datas desde o dia da compra até a data de vencimento da opção e contém todos os caminhos de (S_t), (P_t) e do portfólio. Vamos plotar os resultados obtidos para a evolução apenas da ação, primeiramente.

    brk <- round(sort(c(p_exer, seq(min(port_mc$St),
                                    max(port_mc$St),
                                    length.out = 5))),
                 digits = 2)
    ggplot(port_mc, aes(x = data, y = St)) + 
      geom_line(aes(color = sims)) +
      geom_hline(yintercept = p_exer, color = "red") +
      guides(color = FALSE) +
      labs(title = "Simulações do Valor da Ação",
           x = "data",
           y = "Valor (R$)") +
      scale_y_continuous(breaks = brk) +
      scale_x_date(date_breaks = "2 days", date_labels = "%d") +
      scale_color_viridis_d() +
      theme_economist_white()
    

    plot of chunk gr_acao

    Podemos verificar pela figura acima que a ação, pela nossa SMC, deve fechar na maioria dos caminhos simulados acima do preço de exercício da put (linha vermelha). Entretanto existe uma menor probabilidade de, até a data de vencimento, o preço da ação cair abaixo do strike desta opção.

    Podemos inferir esta probabilidade através do número de caminhos que terminaram em preço da ação abaixo do valor de referência. O custo de proteção contra este risco é o prêmio por nós ao comprarmos a put. O código para esta inferência está abaixo.

    p_baixo <- port_mc %>% 
      filter(data == d_exer) %>% 
      summarise(num_baixo = sum(St < p_exer)) %>% 
      as.double()
    prob <- p_baixo / nsims
    

    Este cálculo nos mostra que em 82 caminhos simulados do preço de PETR4, este terminou abaixo do preço de exercío da opção PETRV17, ou seja, uma probabilidade de 8.2%.

    Para nos precavermos desta possível queda de preço e garantir um valor mínimo de nosso portfólio até a data de 15/10/2018, podemos comprar uma opção de venda, com preço de exercício no valor que desejamos e então o portfólio passa a ser composto pela ação e também pela opção. Caso na data de vencimento o preço da ação seja menor que o preço de exercício da put, esta opção estará ITM e pode ser exercida pelo valor da diferença entre os preços, ou seja, nos garantindo que nosso portfólio estará avaliado em R$ 16,92.

    Esta dinâmica pode ser verificada pela figura abaixo, que agora apresenta o valor do portfólio completo, ação mais opção. Verificamos que, de fato, no dia 15/10/2018 nosso investimento não estará em situação pior que o preço garantido pela compra da put.

    brk <- round(sort(c(p_exer, seq(min(port_mc$port_valor),
                                    max(port_mc$port_valor),
                                    length.out = 5)[-1])),
                 digits = 2)
    ggplot(port_mc, aes(x = data, y = port_valor)) + 
      geom_line(aes(color = sims)) +
      geom_hline(yintercept = p_exer, color = "red") +
      guides(color = FALSE) +
      labs(title = "Simulações do Valor do Portfolio",
           x = "data",
           y = "Valor (R$)") +
      scale_y_continuous(breaks = brk) +
      scale_x_date(date_breaks = "2 days", date_labels = "%d") +
      scale_color_viridis_d() +
      theme_economist_white()
    

    plot of chunk gr_port

    Ou seja, ao custo de 0.6% do preço da ação, compramos uma proteção contra uma queda de preços com probabilidade de 8.2%.

    Esta é apenas uma (simples) aplicação das inúmeras possíveis que a Simulação de Monte Carlo possui no mundo das finanças. A SMC é uma poderosa ferramenta para avaliação e controle de risco de grandes portfólios, com centenas ou milhares de ativos, onde nem sempre consegue-se aferir medidas de retorno esperado ou de risco de mercado de forma analítica.

    Posted by Rafael F. Bressan in Derivativos & Riscos, 1 comment
    Introdução a funções básicas financeiras no R

    Introdução a funções básicas financeiras no R

    **Se você é um estudante de economia, ou um interessado no assunto, e já parou para pensar “Será que vou me ferrar no futuro por não saber programar?”, com certeza este artigo é para você. **

    Afinal, existe uma pitada em você de mistura de medo com curiosidade e, com tempo, pode ter certeza, vai ficando maior. Para isso, estamos aqui para desmistificar sobre os primeiros passos, em específico, para programar em R, um software tão potente quanto amplo, com vasta aplicação no mundo de finanças.

    O que é R

    Direto ao ponto, R é uma linguagem e ambiente para computação estatística e gráficos desenvolvida por John Chambers e colegas na Bell Laboratories, tendo como diferencial a facilidade para plotar gráficos, além de ser de graça e rodar em diversas plataformas, tais como Windows e MacOS.

    Em específico, para o mercado financeiro existem diversas aplicações para análises financeiras, sendo muito mais intuitivas para compreender conceitos da economia, a exemplo, os de econometria.

    A propósito, com o R é mais fácil baixar dados da internet, cortando o tempo perdido e a chatice de ficar entrando nos sites-fonte para baixar planilhas. Enfim, basta explorá-lo que você poderá fazer qualquer coisa – até criar jogos, apesar de não ser o objetivo do programa.

    Para saber mais, indico fortemente fuçar o site do R, além de ler o que está escrito aqui. Neste último texto algumas coisas interessante são citadas como packages (pacotes) e formatação de documentos. Você vai se surpreender quanta coisa dá pra fazer com o R.

    No mais, muitos que já tenham ouvido falar sobre o R podem também se perguntar sobre a importância do seu aprendizado para aplicação no mercado financeiro, uma vez que já estão familiarizados com o Excel, e, às vezes, também com o VBA. Para isso, te respondo.

    R x Excel

    A batalha do século. De um lado, uma plataforma potente e mais complicada, e por outro, uma mundialmente aplicada e aceita, mas com limites. O que escolher?

    De início, podemos refletir que programar é tornar o que é complicado simples e rápido. Ou seja, se ater somente ao Excel não vai ajudá-lo a resolver problemas mais complexos de forma eficiente. E nessa, incluo VBA, que apesar de ser uma maneira de tornar o Excel mais divertido, não chega perto do R, que tem o potencial computacional muito maior.

    Então, se você quer ser mais eficiente com manipulação de dados, sugiro fortemente o R. No mais, caso não conheça o Excel, também sugiro dar uma olhada nele, já que a grande maioria do mundo corporativo o usa como base (existe uma infinidade de curso de graça na internet para Excel).

    Um fator interessante é que o R é bastante intuitivo (não mais que o Excel), com linguagem simples e direta, o que pode ser uma boa maneira para iniciar nos estudos de linguagens de programação.

    De qualquer forma, é uma discussão muito extensa e que alguns até me condenariam por essa comparação, mas que em resumo, explica-se nesse gráfico:

    **Os primeiRos passos **

    Apesar do trocadilho bobo no título, o que virá a seguir poderá ter muito impacto na sua vida e de repente até despertar uma paixão, então siga-os sem medo de errar. Se houver alguma dificuldade, só mandar nos comentários que faremos o possível para ajudá-lo no processo.

    E como tudo começa? Instalando o R…

    Para muitos pode parecer uma coisa simples, mas não é. Não queremos vírus e nem baixar coisas inúteis, que também podem ser vírus.

    É válido comentar que os passos a seguir incluirão baixar a plataforma R Studio que nos facilitará como meio de escrever códigos e resolver problemas mais rápido.

    “RStudio é um conjunto de ferramentas integradas projetadas para ajudá-lo a ser mais produtivo com R. Ele inclui um console, editor de destaque de sintaxe que suporta execução de código direto e uma variedade de ferramentas robustas para traçar, visualizar histórico, depurar e gerenciar seu espaço de trabalho.” – R Studio

    **Etapas para Instalação do R **

  • Baixar o R: https://www.r-project.org/
  • CRAN > Escolha o servidor da UFPR (Brazil) > Seu sistema operacional (no caso Windows) > Install R for the first time > Download R … for Windows

  • Baixar o R Studio: https://www.rstudio.com/products/rstudio/rstudio/download
  • Installers for Supported Plataforms > Seu sistema operacional (no caso Windows)

    Quando você abrir o R Studio, você verá isto abaixo. Confirma? Show!

    Introdução ao ambiente R

    Vamos começar com códigos simples, e aproveite para ir reaplicando no seu R. Antes, devo-lhes explicar uma coisa importante: todo texto ou código que inicie com o símbolo “#” no R não é lido, ou seja, os textos escritos para a situação abaixo com “#” servirão somente como guia, não tendo impacto na leitura do R caso você copie e cole na plataforma.

    Outro ponto importante é que a interface do R está divida em quatro partes as quais vamos nos atentar somente a esquerda-superior(Script) e a esquerda-inferior(Console), sendo que todo código na área Script é salvo, enquanto no Console, é o Script “rodando”, em termos gerais.

    Os códigos abaixo serão aplicados no Script no intuito de se observar o comportamento do R e para “rodá-los”, basta clicar Ctrl + R.

    Aos primeiros códigos

    Vou ensiná-los cinco códigos básicos para aplicação no R e que explicam o funcionamento da plataforma.

    # 1: Adição 
    3 + 5

    Clique Ctrl + R e veja que resulta em 8. Teste com subtração e outras funções matemáticas e veja o resultado. Símbolos: multiplicação “*”, divisão “/”, exponenciação “^”.

    A propósito, o Ctrl + R será necessário sempre que quiser que rode um código, como já mencionado.

    # 2: Guardando informação 
    poupança <- 200

    A palavra poupança agora está vinculada ao valor 200. Desta forma criei uma variável chamada poupança que quando quiser, posso requisitá-la na busca de seus valores ou objetos relacionados.:

    # 3: Requisitando valores/objetos da variável
    print(poupança)

    Veja que uma vez que clicou Ctrl + R tanto para a criação da variável, quanto para o “print”, você perceberá que o valor vinculado a palavra será chamado.

    # 4: Guardando mais de um dado
    Ibm_açoes <-  c(159.82, 160.02, 159.84)

    Agora, se você “printar” a variável Ibm_açoes, perceberá que o conjunto de valores aparecerá.

    #5: Plotando um gráfico
    plot(Ibm_açoes)

    Gráfico de pontinho é feio? Temos uma solução: caso você queira uma linha, acrescente o argumento type=l (l de line), ficando: plot(Ibm_açoes, type=”l”).

    Fim

    Com estas cinco funções espero ter contribuído e atiçado sua vontade de aprender mais sobre R e programação, deixando claro que há uma vastidão de funções, sendo que algumas delas já estão prontas (criadas por outras pessoas membros da comunidade R) que nos facilitam e muito a vida.

    #Plus1: não resisti e coloquei uma aplicação com uso de pacotes
    install.packages(“tidyquant”) #baixando os pacotes da internet
    library(tidyquant) #colocando disponível as funções do pacote no R
    apple <- tq_get("AAPL", get = "stock.prices") #pegando os valores da ação Apple
    plot(apple$date, apple$adjusted, type="l")

    Veja o que acontece. O legal é que você nem precisou ir atrás dos dados. O próprio R fez isso para você através do pacote.

    #Plus2: se ficou em dúvida quanto a alguma fórmula, faça, por exemplo
    ?library

    Para aqueles que se interessaram e querem aprender mais sobre a linguagem, fica de sugestão o Curso de Introdução ao R com aplicabilidade em economia do Vitor Wilher, além da plataforma de estudos em programação Data Camp que tem muita coisa boa para R. Ambos pagos, mas que valem a pena considerar o investimento na troca de menos festinhas no final de semana.

    Para o próximo mês estarei trabalhando em um artigo dedicado somente aos códigos iniciais não abordados aqui para que vocês possam praticar e desenvolver no R, pelo menos um pouco, sem custo algum a não ser seu “gostei” ou “compartilhar” no facebook. Barato não? 

    Seja bem-vindo ao mundo do R.

     

    Posted by Henrique Rosa in Diversos, 5 comments