Inicialmente será necessário importar os pacotes utilizados na simulação. Definimos também alguns parâmetros de visualização dos resultados.
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)
Definimos as datas iniciais e finais para capturar os preços, bem como o código de cada ativo que será utilizado.
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.
data = pd.DataFrame()
for sym in symbols:
data[sym] = web.DataReader(sym, 'google', first_price_date, last_price_date)['Close']
data = data.dropna()
Faremos o cálculo do retorno diário e da volatilidade histórica de cada ativo.
returns = log(data / data.shift(1))
Apresentaremos o retorno dos últimos 14 dias de negociação.
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.
_63_DAYS = 63
_252_DAYS = 252
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()
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.
covariance_matrix = _252_DAYS * returns.tail(_63_DAYS).cov()
sigmas = pd.Series(diag(sqrt(covariance_matrix)), index=symbols)
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
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
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.
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
Será utilizada a Taxa Selic do período considerado (14% ao ano).
risk_free_rate = 0.14
Definição do método que executará a simulação de Monte Carlo.
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.
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
Foi implementada a equação de Back-Scholes para conferir os valores obtidos com a simulação.
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
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.
Abaixo um exemplo de caminho simulado. É possível observar no gráfico a correlação entre os retornos dos ativos.
(_, path) = monte_carlo_simulation(1000, 1)
path.plot(title="A Random Walk")
plt.show()
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