Precificação de Opção Européia com Simulação de Monte Carlo

Passos iniciais

Inicialmente será necessário importar os pacotes utilizados na simulação. Definimos também alguns parâmetros de visualização dos resultados.

In [1]:
import pandas_datareader.data as web
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
import pandas as pd
import datetime as dt
from numpy import *
from scipy.stats import norm

%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (15, 9)

Captura dos Preços

Definimos as datas iniciais e finais para capturar os preços, bem como o código de cada ativo que será utilizado.

In [2]:
first_price_date = dt.datetime(2015, 1, 1)
last_price_date = dt.datetime(2016, 11, 14)
symbols = ['VALE5', 'PETR4', 'BVMF3', 'ITUB4']

Com a ajuda do Pandas, faremos a captura dos preços de fechamento de cada ativo utilizando os dados do Google Finance.

In [3]:
data = pd.DataFrame()
for sym in symbols:
    data[sym] = web.DataReader(sym, 'google', first_price_date, last_price_date)['Close']
data = data.dropna()

Retorno e Volatilidade

Faremos o cálculo do retorno diário e da volatilidade histórica de cada ativo.

In [4]:
returns = log(data / data.shift(1))

Apresentaremos o retorno dos últimos 14 dias de negociação.

In [5]:
last_returns = returns.tail(14)
ax = last_returns.plot(kind='bar', title="Daily Returns")
ax.yaxis.set_major_formatter(FuncFormatter(lambda y, _: '{:.0%}'.format(y)))
ax.set_xticklabels([dt.strftime('%d/%b/%y') for dt in last_returns.index.to_pydatetime()])
plt.gcf().autofmt_xdate()
plt.show()

Apresentaremos a volatilidade histórica do período utilizando uma janela de 63 dias de negociação (3 meses).

É importante destacar que a volatilidade deve ser anualizada. Para isso, consideramos um total de 252 dias úteis no ano.

In [6]:
_63_DAYS = 63
_252_DAYS = 252
In [7]:
vol = returns.rolling(window=_63_DAYS).std() * sqrt(_252_DAYS)

ax = vol.plot(title="3 Months Annualized Volatility")
ax.yaxis.set_major_formatter(FuncFormatter(lambda y, _: '{:.0%}'.format(y)))
plt.show()

Matriz de Covariância

Faremos o cálculo da matriz de covariância e do desvio padrão do retorno de cada ativo. Será considerada a janela dos últimos 63 dias de negociação (3 meses). A matriz de covariância também deverá ser anualizada, assim como foi feito para desvio padrão.

In [8]:
covariance_matrix = _252_DAYS * returns.tail(_63_DAYS).cov()
sigmas = pd.Series(diag(sqrt(covariance_matrix)), index=symbols)
In [9]:
days = returns.tail(_63_DAYS).index
print("From %s to %s \n" % (days[0].strftime("%d/%b/%y"), days[-1].strftime("%d/%b/%y")))
print(pd.DataFrame(100 * sigmas, columns=['3 M Vol (%)']).round(2))
From 15/Aug/16 to 14/Nov/16 

       3 M Vol (%)
VALE5        48.26
PETR4        49.80
BVMF3        30.15
ITUB4        26.69
In [10]:
print(covariance_matrix)
          VALE5     PETR4     BVMF3     ITUB4
VALE5  0.232879  0.111141  0.049092  0.034420
PETR4  0.111141  0.248028  0.099472  0.094024
BVMF3  0.049092  0.099472  0.090913  0.056700
ITUB4  0.034420  0.094024  0.056700  0.071255

Características das Opções

Definiremos para cada opção o tipo (call ou put) e o preço de exercício. A data de vencimento será a mesma para todas as opções.

In [11]:
strike = [22.12, 14.00, 17.84, 36.78]
type = ['put', 'put', 'call', 'call']
expire_date = dt.datetime(2016, 12, 19)

time_to_expiration = pd.date_range(start=last_price_date, end=expire_date,
                                   freq='B'  # allow business day only
                                   ).to_pydatetime().size / _252_DAYS  # in years

Taxa Livre de Risco

Será utilizada a Taxa Selic do período considerado (14% ao ano).

In [12]:
risk_free_rate = 0.14

Simulação de Monte Carlo

Definição do método que executará a simulação de Monte Carlo.

In [13]:
def monte_carlo_simulation(number_steps, number_paths):
    dt = time_to_expiration / number_steps
    Z = random.multivariate_normal(zeros(len(symbols)), covariance_matrix, (number_steps + 1, number_paths))
    option_montecarlo = pd.DataFrame()
    path = pd.DataFrame()
    for idx, sym in enumerate(symbols):
        S0 = data[sym][-1] # Last price is the first price in simulation
        sigma = sigmas[sym] # Standard deviation
        z = Z[:, :, idx] # random values to Brownian Motion
        # Simulating "number_paths" paths with "number_steps" time steps
        S = S0 * exp(cumsum((risk_free_rate - 0.5 * sigma ** 2) * dt + sqrt(dt) * z, axis=0))
        path[sym] = S[:, -1] / S[0, -1]
        # Call or Put Value
        intrinsic = S[-1] - strike[idx] \
            if type[idx] == 'call' \
            else strike[idx] - S[-1]
        option_montecarlo[sym] = exp(-risk_free_rate * time_to_expiration) * maximum(intrinsic, 0) # Present Value
    return (option_montecarlo.sum(axis=0) / number_paths), path

Execução da simulação para três valores diferentes de quantidade de caminhos para o ativo subjacente.

In [14]:
option_value = pd.DataFrame()
for k in [10, 100, 10000]:
    (option_montecarlo, _) = monte_carlo_simulation(1000, k)
    option_value['MC %d' % k] = option_montecarlo.round(2)

print(option_value)
       MC 10  MC 100  MC 10000
VALE5   0.27    0.36      0.43
PETR4   0.23    0.82      0.80
BVMF3   0.14    0.27      0.32
ITUB4   0.36    0.42      0.49

Equação de Black-Scholes

Foi implementada a equação de Back-Scholes para conferir os valores obtidos com a simulação.

In [15]:
prices = pd.Series(index=symbols)
T = time_to_expiration
for idx, sym in enumerate(symbols):
    S = data[sym][-1]
    K = strike[idx]
    sigma = sigmas[sym]

    d1 = (log(S / K) + (risk_free_rate + 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
    d2 = d1 - sigma * sqrt(T)

    price = S * norm.cdf(d1) - K * exp(-risk_free_rate * T) * norm.cdf(d2) \
        if type[idx] == 'call' \
        else K * exp(-risk_free_rate * T) * norm.cdf(-d2) - S * norm.cdf(-d1)

    prices[sym] = price
In [16]:
option_value['BS'] = prices.round(2)
option_value['Type'] = type
print(option_value)
       MC 10  MC 100  MC 10000    BS  Type
VALE5   0.27    0.36      0.43  0.42   put
PETR4   0.23    0.82      0.80  0.79   put
BVMF3   0.14    0.27      0.32  0.33  call
ITUB4   0.36    0.42      0.49  0.51  call

Os preços obtidos na simulação se aproximam do preço calculado pela equação de Black-Scholes à medida que aumentamos o número de caminhos simulados.

Um Exemplo de Caminho do Preço dos Ativos

Abaixo um exemplo de caminho simulado. É possível observar no gráfico a correlação entre os retornos dos ativos.

In [17]:
(_, path) = monte_carlo_simulation(1000, 1)
path.plot(title="A Random Walk")
plt.show()
In [18]:
print(returns.tail(_63_DAYS).corr().round(2))
       VALE5  PETR4  BVMF3  ITUB4
VALE5   1.00   0.46   0.34   0.27
PETR4   0.46   1.00   0.66   0.71
BVMF3   0.34   0.66   1.00   0.70
ITUB4   0.27   0.71   0.70   1.00