# Derivativos de Ações e Commodities
##  Precificação com o Método de Diferenças Finitas
<sub>Uirá Caiado. 29 de Junho, 2016<sub>

#### Resumo

_Neste projeto vou implementar o método de diferenças finitas explícito para precificar diferentes derivativos. Como o preço de opções geralmente é descrito por equações diferenciais de difusão ou parabólicas, o método de diferenças finitas se mostra adequado, uma vez que é utilizado justamente para encontrar soluções numéricas para equações diferencias. Depois de implementar o modelo, vou compará-lo com os resultados obtidos pela solução analítica de cada instrumento._

## 1. Introdução

Nesta seção dou uma breve descrição sobre o método de diferenças finitas e declaro o problema que será abordado.

### 1.1. O Método Utilizado

Como colocado por Wilmott, encontrar soluções fechadas para precificação de muitos derivativos pode ser difícil, ou mesmo inviável. Como o preço destes instrumentos frequentemente são descritos por equações diferencias parciais, frequentemente utilizam-se métodos numéricos para encontrar uma solução, como [árvores binomiais](https://en.wikipedia.org/wiki/Binomial_options_pricing_model), simulação de [monte carlo](https://en.wikipedia.org/wiki/Monte_Carlo_method) e diferenças finitas.

O método de diferenças finitas encontra o valor de um derivativo calculando-o em todo domínio (factivel) de preços do instrumento base, incluíndo a passagem de tempo até seu vencimento.  Assim, são similares às árvores binomiais. Porém, ao invéz e discretizar os preços do ativo e a passagem do tempo em uma estrutura de árvore, discretiza em um grid.

Diferenças finitas também são muitos boas para lidar com problemas com poucas dimensões e equações diferenciais não lineares (o preço e tempo já são duas, e o Wilmott sugere até quatro). Para muitas dimensões, a implementação começa a ficar complicada e pouco eficiente, coisa que o Monte Carlo lida melhor. Porém, se há exercício antecipado e não linearidade, o método de diferenças finitas acaba se mostrando como solução mais viável.

### 1.2. O Problema

Considerando um ativo-objeto cuja dinâmica do preço satisfaz a seguinte EDE:

$$\frac{\mathrm{d} S_t}{S_t}=\mu\cdot \mathrm{d}t + \sigma\cdot \mathrm{d}W_{t}$$

Deseja-se calcular o preço justo de um derivativo um derivativo com característica europeia cujo payoff é descrito por uma função qualquer $V_T=V(T, S_T)$, onde $T$ é o vencimento do derivativo. $S_T$ é o preço do ativo-objeto em $T$ e é possível negociar qualquer quantidade dele em qualquer instante. Não há custo de transação (corretagem, emolumento, bid-ask spread, etc) e posições vendidas a descoberto no subjacente são permitidas, não havendo custos associados.

Pede-se que se implemente o algoritmo de diferenças finitas explícito e se calcule o preço justo, o Delta e o Gamma de cada instrumento abaixo. Os valores encontrados devem ser comparados com o resultado de suas expressões analítica. A simulação deve ser feita para os payoffs abaixo. $K$ é o Strike da opção.
- $V(T, S_T)=ln(S_T)$
- $V(T, S_T)=(ln(S_T))^2$
- $V(T, S_T)=(S_T-K)^2$
- $V(T, S_T)=\mathbf{1}_{S_T > K}$
- $V(T, S_T)=max(S_T-K, 0)$

##  2. Diferenças Finitas

Abaixo, apresentarei os elementos necessários para implementar o método diferenças finitas, sendo eles: a [diferenciação](https://pt.wikipedia.org/wiki/Diferenciação_numérica) das derivadas parciais necessárias usando o Grid; a discretização da condição final para cada derivativo; e suas respectivas condições de contorno (nos limites do grid).

### 2.1. Diferenciação no Grid

De acordo com Wilmott, o grid utilizado pelo método de diferenças finitas tem passos de tempo e de preço (ou log do preço) geralmente homogênios. Porém, **não há** restrição para a forma do Grid, desde que os ajustes necessários sejam feitos. Como exposto por 

Considerando que podemos discretizar $t$ como $t=T-k\delta t$ e do $S$ como $S = i\delta S$, onde $i$ e $k$ são seus respectivos passos no grid, podemos escrever o valor da opção em cada ponto do grig como sendo:

$$V^{k}_{i} = ( i \delta S, T-k \delta t )$$

Seguindo notas de aula, aplicando Taylor ao preço $V(T, S_T)$ de um derivativo genérico, posso escrever a seguinte equação diferencial parabólica:

$$\frac{\partial V}{\partial t} + a(S, t) \cdot \frac{\partial V^2}{\partial S^2} + b(S, t) \cdot \frac{\partial V}{\partial S} + c(S, t)\cdot V = 0 $$

Sendo que $\frac{\partial V}{\partial t}$ também é chamado de theta ($\theta$), $\frac{\partial V}{\partial S}$ de delta ($\Delta$) e $\frac{\partial V^2}{\partial S^2}$ de gamma ($\Gamma$). Como demostrado por Wilmott, como a definição da primeira derivada $V$ em relação a $t$ é dado por

$$\frac{\partial V}{\partial t} = \underset{h \, \rightarrow \, 0}{\lim}\, \frac{V(S, t + h) - V(S, t)}{h}$$

Então podemos aproximar $\theta$ como sendo $\frac{\partial V}{\partial t} (S, t) \approx \frac{V_{i}^{k} - V_{i}^{k+1} }{\delta t}$. Note que $i$, que é relacionado ao passo do ativo, ficou fixo. *Também estou ignorando os erros de aproximação, como ignorarei em todas as outras diferenciações*.

O mesmo raciocínio pode ser utilizado para aproximar o $\Delta$. Porém, Wilmott ainda sugere que se utilize a [diferença centrada](http://math.stackexchange.com/questions/888259/can-someone-explain-in-general-what-a-central-difference-formula-is-and-what-it), que é dada por $\frac{\partial V}{\partial S} (S, t) \approx \frac{V_{i+1}^{k} - V_{i-1}^{k} }{2 \delta S}$. A discretização anterior também poderia ser utilizada, mas esta oferece um erro de aproximação menor. O único problema desta abordagem é que preciso saber os valores $S + \delta S$ e $S - \delta S$ e nas fronteiras do Grid não terei esta informação. Por tanto, nestes casos utilizarei a discretrização utilizando apenas um dos lados.

Por fim, para achar o $\gamma$, Wilmott subtrai o delta forward do delta backward e divide o resultado por $\delta S$, chegando na aproximação  $\frac{\partial V^2}{\partial S^2} (S, t) \approx \frac{V_{i+1}^{k} - 2 V_{i}^{k} + V_{i-1}^{k} }{\delta S^2}$. Como demonstrado em notas de aula, todos estes resultados podem ser checados aplicando a expansão de Taylor.

### 2.2. Condição Terminal e Payoffs

O valor da opção no vencimento é simplemente seu payoff. Assim, não é necessário resolver nada para $T$, apenas para $S$. Usando a notação de diferenças finitas, os payoffs desejados ficam sendo da forma:

- $V_{i}^{0}=ln(i\delta S)$
- $V_{i}^{0}=(ln(i\delta S))^2$
- $V_{i}^{0}=(i\delta S-K)^2$
- $V_{i}^{0}=\mathbf{1}_{i\delta S > K}$
- $V_{i}^{0}=max(i\delta S-K, 0)$

Como o Wilmott explica, o método de diferenças finitas começa de trás para frente, será destes valores que a iteração começará, como se estivéssemos calculando o preço de um derivativo por árvores binomiais, também de trás para frente.

### 2.3. Condições de Contorno

Quando estivermos percorrendo o Grid, necessitaremos definir o preço do derivativo em seus extremos, quando $S=0$ e $S=I \delta S$, onde $I$ é o ponto mais alto do Grid. Esta condição depende do instrumento que está sendo precificado. Para **todos** dos contratos, utilizaremos a seguinte condição para quando $S=0$:

$$V^{k}_{0} = (1 - r \delta t)V^{k-1}_{0}$$

Para a condição de contorno superior, utilizaremos para a **maioria dos contratos**:

$$V^{k}_{I} = 2V^{k}_{I-1} - V^{k}_{I-2}$$

Esta condição é adequada pois, a medida que $S \rightarrow \infty$, o $\Gamma$ da maioria dos contratos tende a zero. Porém, isso não é verdade para o contrato cujo payoff é descrito por $V_{i}^{0}=(i\delta S-K)^2$. Diferenciando novamente o [delta](https://nbviewer.jupyter.org/github/ucaiado/Replicating_Strategy/blob/master/UiraCaiadoEx01.ipynb) deste contrato, chegamos que o gamma será dado por $2 \cdot e^{(r + \sigma^2)(T-t)}$. Neste caso,
partindo da equação discretizada do $\Gamma$ e resolvendo para $V^{k}_{I}$, tenho que:

\begin{equation}
\begin{aligned}
\Gamma &= \frac{V_{I+1}^{k} - 2 V_{I}^{k} + V_{I-1}^{k} }{\delta S^2}\\
2 \cdot e^{(r + \sigma^2)(\delta t)} &= \frac{V_{I+1}^{k} - 2 V_{I}^{k} + V_{I-1}^{k} }{\delta S^2}\\
V_{I}^{k} &= 2 \delta S^2 e^{(r + \sigma^2)(\delta t)} + 2V_{I-1}^{k} - V_{I-2}^{k}\\
V_{I}^{k} &\approx 2 \delta S^2 + 2V_{I-1}^{k} - V_{I-2}^{k}\\
\end{aligned}
\end{equation}



## 3. Implementando o Modelo

Nesta seção vou decsrever o método de diferenças finitas explícito e implementar os códigos necessários.

### 3.1. Método Explícito

Seguindo notas de aula, substituindo as diferenciações encontradas na equação parabólica mencionada anteriormente e reescrevendo os outros termos com a notação de diferenças finitas, temos que:

\begin{equation}
\begin{aligned}
\frac{\partial V}{\partial t} + a(S, t) \cdot \frac{\partial V^2}{\partial S^2} + b(S, t) \cdot \frac{\partial V}{\partial S} + c(S, t)\cdot V &= 0 \\
\frac{V_{i}^{k} - V_{i}^{k+1} }{\delta t} + a_{i}^{k} \cdot \frac{V_{i+1}^{k} - 2 V_{i}^{k} + V_{i-1}^{k} }{\delta S^2} + \\b_{i}^{k} \cdot \frac{V_{i+1}^{k} - V_{i-1}^{k} }{2 \delta S} + c_{i}^{k}\cdot V_{i}^{k} &= 0
\end{aligned}
\end{equation}

Rearranjando equação acima para isolar $V_{i}^{k+1}$ e renomeando alguns termos, ficamos com:

$$V_{i}^{k+1} = A_{i}^{k}V_{i-1}^{k} + \left( 1 + B_{i}^{k} \right)V_{i}^{k} + C_{i}^{k}V_{i+1}^{k}$$

Onde:

\begin{equation}
\begin{aligned}
A_{i}^{k} &= \nu_1 a_{i}^{k} - 0.5\nu_2 b_{i}^{k} \\
B_{i}^{k} &= -2\nu_1 a_{i}^{k} - \delta t c_{i}^{k} \\
C_{i}^{k} &= \nu_1 a_{i}^{k} + 0.5\nu_2 b_{i}^{k}
\end{aligned}
\end{equation}

Sendo que $\nu_1=\frac{\delta t}{\delta S ^2}$ e $\nu_2=\frac{\delta t}{\delta S}$. A equação acima está definida apenas entre $i=1, ..., I-1$. OS pontos restantes necessários para discretização vem das condições de contorno. Como nós conhecemos o valor terminal em $V_{i}^{0}$, podemos calcular o valor de $V_{i}^{1}$ e assim por diante. Como o valor do instrumento em $k+1$ só depende dos valores dele em $k$, chamamos este método de **método explícito**.

### 3.2. Outros Detalhes

#### a. Estabilidade

Para que a solução deste método seja estável, é necessário que satisfaça $\delta t \, \leqslant \frac{\delta S^2}{2\sigma^2 S^2}$ e $\delta S \leqslant \frac{2a}{|b|}$. Vou utilizar a primeira delas explicitamente na implementação, porém vou fazer a seguinte transformação:

$$\delta t \, \leqslant \frac{\delta S^2}{2\sigma^2 S^2} = \frac{1}{\sigma^2} \left(\frac{\delta S}{S}\right)^2 = \frac{1}{\sigma^2 I^2}$$

Onde $I = \delta  S \times S$.

#### b. Interpolação Bilinear

Devido a discretização do modelo, pode ser que seja necessário obter um valor que esteja entre os nós do Grid. Nestes casos, adotarei uma interpolação de duas dimensões (Preço e tempo) chamada [interpolação bilinear](https://en.wikipedia.org/wiki/Bilinear_interpolation). Ela consiste em realizar uma interpolação linear primeiro em uma direção e depois na outra. Considerando um quadrante específico do Grid, em que há quatro preços de opção, basicamente se divide o quadrante em quatro (onde os pontos de divisão são os valores desejados) e se utiliza a área de cada um deles como peso para os valores das opções. Assim, terminamos com a seguinte fórmula:

$$\frac{\sum_{i=1}^{4} A_i V_i}{\sum_{i=1}^{4} A_i}$$

Onde $A$ é a área dos retângulos e $V$ o preço da opção em cada canto do quadrante principal. Cada $A_i$ é oposto ao $V_i$ correspondente (o $A_2$ é oposto ao $V_2$, por exemplo).

### 3.3. Criando o Grid

Primeiro, vou definir as classes para a criação do Grid. Como quero acessar os valores posteriormente, vou manter todos em uma estrutura, ainda que este provavelmente não seja o método mais adequado pensando em velocidade.


```python
import finite_difference; reload(finite_difference);
x = finite_difference.Grid(f_vol=0.5,
                           f_value=100.,
                           f_time=1.,
                           i_nas=10,
                           i_nts=10)
print x
```

           0     1     2     3     4     5     6     7     8     9
    0    0,0   0,1   0,2   0,3   0,4   0,5   0,6   0,7   0,8   0,9
    1    1,0   1,1   1,2   1,3   1,4   1,5   1,6   1,7   1,8   1,9
    2    2,0   2,1   2,2   2,3   2,4   2,5   2,6   2,7   2,8   2,9
    3    3,0   3,1   3,2   3,3   3,4   3,5   3,6   3,7   3,8   3,9
    4    4,0   4,1   4,2   4,3   4,4   4,5   4,6   4,7   4,8   4,9
    5    5,0   5,1   5,2   5,3   5,4   5,5   5,6   5,7   5,8   5,9
    6    6,0   6,1   6,2   6,3   6,4   6,5   6,6   6,7   6,8   6,9
    7    7,0   7,1   7,2   7,3   7,4   7,5   7,6   7,7   7,8   7,9
    8    8,0   8,1   8,2   8,3   8,4   8,5   8,6   8,7   8,8   8,9
    9    9,0   9,1   9,2   9,3   9,4   9,5   9,6   9,7   9,8   9,9
    10  10,0  10,1  10,2  10,3  10,4  10,5  10,6  10,7  10,8  10,9



```python
%timeit x = Grid(f_vol=0.2, f_value=100., f_time=1., i_nas=100)
```

    1000 loops, best of 3: 344 µs per loop


Como comentei, não é uma estrtutura eficiente. Demorou quase 5 segundos para criar uma estrtutura de 10.000 nós. Porém acredito que isso me ajudará adiante.

### 3.3. Testando Método

Seguindo implementação do livro, também há outra opção para se encontrar $V_{i}^{k+1}$. Substituindo a notação da fórmula de precificação e assumindo que possuimos os valores para $\Delta$, $\Gamma$ e $V_{i}^{k}$, podemos isolar o $\theta$, ficando com


\begin{equation}
\begin{aligned}
\theta + a_{i}^{k} \cdot \Gamma + b_{i}^{k} \cdot \Delta + c_{i}^{k}\cdot V_{i}^{k} &= 0\\
-a_{i}^{k} \cdot \Gamma - b_{i}^{k} \cdot \Delta - c_{i}^{k}\cdot V_{i}^{k} &= \theta
\end{aligned}
\end{equation}

Assim, para uma opção que não paga dividendos, que o ativo base segue um movimento brawniano geométrico, e esta sendo avaliado no mundo neutro a risco, temos que o parâmetro $c=-r$, $b=rS$ e $a=0.5 \cdot \sigma^2 S^2$. Assim, fico com o theta sendo da forma:

$$\theta = r V_{i}^{k} - rS_{i}^{k}\Delta - 0.5\sigma^2 S^2\Gamma$$

Então, o preço do derivativo em $k+1$ pode ser dado por:

$$V_{i}^{k+1} = V_{i}^{k} - \delta t \theta$$

Para comparar os resultados obtidos, também serão calculados os preços e as gregas de cada opção analiticamente, seguindo as equações de preficição demostradas em [trabalho](https://nbviewer.jupyter.org/github/ucaiado/Replicating_Strategy/blob/master/UiraCaiadoEx01.ipynb) anterior. Abaixo vou imprimir os resultados obtidos para uma Call Européia para comparação com exemplos do Wilmott.


```python
import finite_difference; reload(finite_difference);
d_param = {"f_St": 100.,  # preco do ativo
           "f_sigma": 0.2, # desvio padra do ativo objeto
           "f_time": 1.,  # tempo para vencimento em anos
           "f_r": 0.05,  # taxa de juros anual
           "i_nas": 20,  # passos que o ativo sera discretizado
           "i_nts": 10,  # passos que o tempo sera discretizado (opcional)
           "f_K": 100.  # strike da opcao
          }

%time self = finite_difference.EuropianCall(**d_param)
```

    Wall time: 582 ms



```python
# imprime apenas a parte da matriz demonstrada no exemplo do livro, pag 1216, tabela 77.1
df = self.df_opt_prices.copy()
df.columns = ["{:.3f}".format(x) for x in df.columns]
df.tail(15).head(9)
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0.000</th>
      <th>0.111</th>
      <th>0.222</th>
      <th>0.333</th>
      <th>0.444</th>
      <th>0.556</th>
      <th>0.667</th>
      <th>0.778</th>
      <th>0.889</th>
      <th>1.000</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>60.0</th>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.001046</td>
      <td>0.004769</td>
      <td>0.012912</td>
      <td>0.027040</td>
      <td>0.048408</td>
    </tr>
    <tr>
      <th>70.0</th>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.010816</td>
      <td>0.040306</td>
      <td>0.092362</td>
      <td>0.168004</td>
      <td>0.266547</td>
      <td>0.386414</td>
    </tr>
    <tr>
      <th>80.0</th>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.084278</td>
      <td>0.248619</td>
      <td>0.475041</td>
      <td>0.746832</td>
      <td>1.051237</td>
      <td>1.379044</td>
      <td>1.723707</td>
    </tr>
    <tr>
      <th>90.0</th>
      <td>0.0</td>
      <td>0.000000</td>
      <td>0.512500</td>
      <td>1.148000</td>
      <td>1.807437</td>
      <td>2.461115</td>
      <td>3.100759</td>
      <td>3.724696</td>
      <td>4.333380</td>
      <td>4.927851</td>
    </tr>
    <tr>
      <th>100.0</th>
      <td>0.0</td>
      <td>2.500000</td>
      <td>4.013889</td>
      <td>5.200154</td>
      <td>6.223712</td>
      <td>7.149693</td>
      <td>8.008182</td>
      <td>8.816223</td>
      <td>9.584443</td>
      <td>10.320129</td>
    </tr>
    <tr>
      <th>110.0</th>
      <td>10.0</td>
      <td>10.555556</td>
      <td>11.571451</td>
      <td>12.561620</td>
      <td>13.500820</td>
      <td>14.389205</td>
      <td>15.235188</td>
      <td>16.045095</td>
      <td>16.824334</td>
      <td>17.576942</td>
    </tr>
    <tr>
      <th>120.0</th>
      <td>20.0</td>
      <td>20.555556</td>
      <td>21.108025</td>
      <td>21.790273</td>
      <td>22.510062</td>
      <td>23.243401</td>
      <td>23.974649</td>
      <td>24.698653</td>
      <td>25.412183</td>
      <td>26.114374</td>
    </tr>
    <tr>
      <th>130.0</th>
      <td>30.0</td>
      <td>30.555556</td>
      <td>31.108025</td>
      <td>31.657425</td>
      <td>32.248867</td>
      <td>32.862026</td>
      <td>33.491183</td>
      <td>34.127231</td>
      <td>34.766578</td>
      <td>35.405282</td>
    </tr>
    <tr>
      <th>140.0</th>
      <td>40.0</td>
      <td>40.555556</td>
      <td>41.108025</td>
      <td>41.657425</td>
      <td>42.203772</td>
      <td>42.764972</td>
      <td>43.335178</td>
      <td>43.915298</td>
      <td>44.500544</td>
      <td>45.090213</td>
    </tr>
  </tbody>
</table>
</div>



O título das colunas se refere ao tempo, começando no vencimento (tempo zero) até $t$ (o dia em que a opção está sendo avaliada). As linhas se referem ao preço do ativo objeto. *Note que começa a aparecer uma diferença em relação ao livro quando o ativo está valendo $R\$ 110.00$  ou menos. Como o método utiliza os três nós subjacentes anteriores para calcular o valor do nó posterior, o nó $(\,90.0 \,; 0.111)$ não deveria apresentar valor algum de fato. Porém, no exemplo do Wilmott, apresentou.*

## 4. Precificando Instrumentos

Nesta seção vou precificar alguns instrumentos e comparar com suas soluções analíticas.

### 4.1. Call Européia

Partindo de implementação anterior, apenas resta imeplementar o gamma analítico da equação, que é dado abaixo:
$$\frac{\partial V^2}{\partial S^2} = \frac{\Delta}{\partial S} = \frac{N'(d1)}{S \sigma \sqrt{t}} $$

Abaixo vou comparar os resultados analíticos com o que obtive do método de diferenças finitas.


```python
import finite_difference; reload(finite_difference);
d_param = {"f_St": 100.,  # preco do ativo
           "f_sigma": 0.5, # desvio padra do ativo objeto
           "f_time": .5,  # tempo para vencimento em anos
           "f_r": 0.10,  # taxa de juros anual
           "i_nas": 20,  # passos que o ativo sera discretizado
           "f_K": 100.  # strike da opcao
          }

%time my_option = finite_difference.EuropianCall(**d_param)
```

    CPU times: user 2.24 s, sys: 23.1 ms, total: 2.26 s
    Wall time: 2.28 s



```python
# plotando graficos do inicio do periodo
l_prices = np.arange(50., 151., 1.)
my_option.compare_to_analytical_solutions(l_prices, d_param['f_time'])
```

### 4.2. Contrato Log

Derivando gamma analítico da equação a paritr do Delta $\Delta = \frac{e^{-r\cdot t}}{S}$, tenho que:

\begin{equation}
\begin{aligned}
\frac{\partial V^2}{\partial S^2} &= \frac{\Delta}{\partial S} = -\frac{e^{-r\cdot t}}{S^2}
\end{aligned}
\end{equation}


Abaixo vou comparar os resultados analíticos com o que obtive do método de diferenças finitas. A única diferença que é necessária tomar nesta implementação é quando o valor do ativo no final é $0$. Nestes casos, fiz o valor da opção também ser $0$.


```python
import finite_difference; reload(finite_difference);
d_param = {"f_St": 100.,  # preco do ativo
           "f_sigma": 0.5, # desvio padra do ativo objeto
           "f_time": .5,  # tempo para vencimento em anos
           "f_r": 0.10,  # taxa de juros anual
           "i_nas": 20,  # passos que o ativo sera discretizado
           "f_K": 100.  # strike da opcao
          }

%time my_option = finite_difference.LogContract(**d_param)
```

    CPU times: user 1.83 s, sys: 11.1 ms, total: 1.84 s
    Wall time: 1.84 s



```python
# plotando graficos do inicio do periodo
l_prices = np.arange(50., 151., 1.)
my_option.compare_to_analytical_solutions(l_prices, d_param['f_time'])
```

### 4.3. Contrato Log Quadrático

Derivando gamma analítico da equação e aplicando a regra de quociente, tenho que, se
$$\Delta = \frac{2e^{-r\cdot t}}{S}\left[ln(S)+ \left(r - \frac{\sigma^2}{2} \right) \cdot t \right]$$

Então

\begin{equation}
\begin{aligned}
\frac{\partial V^2}{\partial S^2} &= \frac{\Delta}{\partial S}= 2e^{-r\cdot t} \frac{\partial}{\partial S} \left [ \frac{ln(S)+ t \left(r - \frac{\sigma^2}{2} \right)}{S} \right ]\\
&= 2e^{-r\cdot t} \frac{S \frac{\partial}{\partial S} \left ( ln(S)+ t \left(r - \frac{\sigma^2}{2} \right) \right) - \left( ln(S)+ t \left(r - \frac{\sigma^2}{2} \right) \right) \frac{\partial}{\partial S} (S)}{S^2}\\
&= 2e^{-r\cdot t} \frac{S \cdot \frac{1}{S} - ln(S)- tr + \frac{t \sigma^2}{2}}{S^2}\\
&= \frac{e^{-r t}}{S^2}\left[2 + \sigma^2 t - 2ln(S) -2rt\right]
\end{aligned}
\end{equation}


Abaixo vou comparar os resultados analíticos com o que obtive do método de diferenças finitas.


```python
import finite_difference; reload(finite_difference);
d_param = {"f_St": 100.,  # preco do ativo
           "f_sigma": 0.5, # desvio padra do ativo objeto
           "f_time": .5,  # tempo para vencimento em anos
           "f_r": 0.10,  # taxa de juros anual
           "i_nas": 20,  # passos que o ativo sera discretizado
           "f_K": 100.  # strike da opcao
          }

%time my_option = finite_difference.SquaredLogContract(**d_param)
```

    CPU times: user 2.02 s, sys: 85.9 ms, total: 2.11 s
    Wall time: 2.09 s



```python
# plotando graficos do inicio do periodo
l_prices = np.arange(50., 151., 1.)
my_option.compare_to_analytical_solutions(l_prices, d_param['f_time'])
```

### 4.4. Exótico Quadrático

Derivando gamma analítico da equação a partir de $\Delta = 2S e^{(r + \sigma^2)} - 2K$, tenho que:

\begin{equation}
\begin{aligned}
\frac{\partial V^2}{\partial S^2} &= 2 e^{(r + \sigma^2)}
\end{aligned}
\end{equation}


Neste caso, também preciso **mudar as condições de contorno**. Abaixo vou comparar os resultados analíticos com o que obtive do método de diferenças finitas.


```python
import finite_difference; reload(finite_difference);
d_param = {"f_St": 100.,  # preco do ativo
           "f_sigma": 0.5, # desvio padra do ativo objeto
           "f_time": .5,  # tempo para vencimento em anos
           "f_r": 0.10,  # taxa de juros anual
           "i_nas": 20,  # passos que o ativo sera discretizado
           "f_K": 100.  # strike da opcao
          }

%time my_option = finite_difference.SquaredExotic(**d_param)
```

    CPU times: user 1.82 s, sys: 34.3 ms, total: 1.86 s
    Wall time: 1.85 s



```python
# plotando graficos do inicio do periodo
l_prices = np.arange(50., 151., 1.)
my_option.compare_to_analytical_solutions(l_prices, d_param['f_time'])
```

### 4.5. Opção Digital

Por último, vou compara o resultado do método explícito com o resultado analítico de uma opção digital. Seguindo a derivação do gamma analítico a [daqui](http://www.stat.nus.edu.sg/~stalimtw/MFE5010/PDF/L2digital.pdf), tenho que:

\begin{equation}
\begin{aligned}
\frac{\partial V^2}{\partial S^2} &= - \frac{e^{-rt} d_1 N(d_2) }{S^2 \sigma^2 t }
\end{aligned}
\end{equation}


Abaixo vou plotar os resultados para comparação


```python
import finite_difference; reload(finite_difference);
d_param = {"f_St": 100.,  # preco do ativo
           "f_sigma": 0.5, # desvio padra do ativo objeto
           "f_time": .5,  # tempo para vencimento em anos
           "f_r": 0.10,  # taxa de juros anual
           "i_nas": 20,  # passos que o ativo sera discretizado
           "f_K": 100.  # strike da opcao
          }

%time my_option = finite_difference.DigitalOption(**d_param)
```

    CPU times: user 2.46 s, sys: 92.1 ms, total: 2.55 s
    Wall time: 2.61 s



```python
# plotando graficos do inicio do periodo
# plotando graficos do inicio do periodo
l_prices = np.arange(50., 151., 1.)
my_option.compare_to_analytical_solutions(l_prices, d_param['f_time'])
```

A solução não aderiu muito bem às soluções analíticas. Vou tentar melhorar a quantidade de discretizaçòes do ativo subjacente para tentar melhorar este resultado


```python
import finite_difference; reload(finite_difference);
d_param = {"f_St": 100.,  # preco do ativo
           "f_sigma": 0.5, # desvio padra do ativo objeto
           "f_time": .5,  # tempo para vencimento em anos
           "f_r": 0.10,  # taxa de juros anual
           "i_nas": 100,  # passos que o ativo sera discretizado
           "f_K": 100.  # strike da opcao
          }

%time my_option = finite_difference.DigitalOption(**d_param)
```

    CPU times: user 4min 26s, sys: 1.62 s, total: 4min 28s
    Wall time: 4min 32s



```python
# plotando graficos do inicio do periodo
# plotando graficos do inicio do periodo
l_prices = np.arange(50., 151., 1.)
my_option.compare_to_analytical_solutions(l_prices, d_param['f_time'])
```

Curiosamente, a melhorar foi marginal. A implementação que fiz também se mostrou inadequada, uma vez que quando criei um grid com $100$ passos no ativo, demorou 4 minutos e meio para realizar todos os cálculos.

## 5. Conclusão

Neste trabalho o método de diferenças finitas explícito foi comparado à precificação de instrumentos com solução analítica. De maneira geral, o método conseguiu aproximar bem o preço das opções e se mostrou bastante flexível, precisando de poucos ajustes para precificar cada instrumento. Apesar de não testado aqui, provavelmente a inclusão de não linearidades nos contratos também não exigiria muitos ajustes, como a utilização de parâmetros incertos, ou mesmo exercício antecipado.


## 6. Últimas Considerações

A decisão de usar classes no lugar de um array dificultou a implementação. A aplicação do método criando um objeto para cada nó do Grid também se mostrou muito lenta. Porém, foi interessante notar que o simples aumento da discretização do grid não melhorou tanto a aderência aos valores calculados pelas soluções analíticas. Talvez a mudança do formato do grid mostrasse melhores resultados.

_Style notebook and change matplotlib defaults_


```python
#loading style sheet
from IPython.core.display import HTML
HTML(open('ipython_style.css').read())
```




<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunss.otf');
    }
    @font-face {
        font-family: "Computer Modern";
        font-weight: bold;
        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunsx.otf');
    }
    @font-face {
        font-family: "Computer Modern";
        font-style: oblique;
        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunsi.otf');
    }
    @font-face {
        font-family: "Computer Modern";
        font-weight: bold;
        font-style: oblique;
        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunso.otf');
    }
    div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    }
    h1 {
        font-family: Helvetica, serif;
    }
    h4{
        margin-top:12px;
        margin-bottom: 3px;
       }
    div.text_cell_render{
        font-family: Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 145%;
        font-size: 130%;
        width:800px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro", source-code-pro,Consolas, monospace;
    }
    .prompt{
        display: None;
    }
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 22pt;
        color: #4057A1;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }
</style>



<style>
    table {
        overflow:hidden;
        font-family: Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        font-size: 12px;
        margin: 10px;
        width: 480px;
        text-align: left;
        border-collapse: collapse;
        border: 1px solid #d3d3d3;
        -moz-border-radius:5px; FF1+;
        -webkit-border-radius:5px; Saf3-4;
        border-radius:5px;
        -moz-box-shadow: 0 0 4px rgba(0, 0, 0, 0.01);
    }
    th
    {
        padding: 12px 17px 12px 17px;
        font-weight: normal;
        font-size: 14px;
        border-bottom: 1px dashed #69c;
    }

    td
    {
        padding: 7px 17px 7px 17px;

    }

    tbody tr:hover th
    {

        background:  #E9E9E9;
    }

    tbody tr:hover td
    {

        background:  #E9E9E9;
    }

</style>







```python
#changing matplotlib defaults
%matplotlib inline
import seaborn as sns
sns.set_palette("deep", desat=.6)
sns.set_context(rc={"figure.figsize": (8, 4)})
sns.set_style("whitegrid")
sns.set_palette(sns.color_palette("Set2", 10))
```


```python

```

dbc717a729a1a2769b1dfbfa71d3fd83e4cfa334
ipynb
Jupyter Notebook
UiraCaiadoEx02.ipynb
ucaiado/FiniteDifference_Pricing
5a7e80eb37077fce90b70d55c80c995f1592e4c8
MIT
2018-02-24T16:23:51.000Z
2022-02-15T09:57:21.000Z
UiraCaiadoEx02.ipynb
ucaiado/FiniteDifference_Pricing
5a7e80eb37077fce90b70d55c80c995f1592e4c8
MIT
UiraCaiadoEx02.ipynb
ucaiado/FiniteDifference_Pricing
5a7e80eb37077fce90b70d55c80c995f1592e4c8
MIT
2018-05-03T14:28:40.000Z
2021-12-21T12:28:55.000Z
Qwen/Qwen-72B
1. YES
2. YES