Esta questão foi resolvida no Python.
A equação de Black-Scholes é:
\( C = N(d_1)S_0 - N(d_2)X e^{-rT} \)
Onde C é o valor real da opção, \( S_0 \) é o fluxo de caixa futuro descontado (DFC), \( X \) o custo do investimento, \( r \) a taxa de juros contínua e \( T \) o tempo do projeto. \( d_1 \) e \( d_2 \) são definidos abaixo:
import math
def d1 (S0,X, r, vol,T):
return (math.log(S0/X) + (r+0.5*(vol**2))*T)/(vol*(T**(1/2)))
def d2(d1,vol,T):
return d1(S0,X, r, vol,T)-(vol*(T**(1/2)))
def Nd1 (d1, xbarra, dp): #normalizando d1
return (d1(S0,X, r, vol,T)-xbarra)/dp
def Nd2 (d2, d1, xbarra,dp): #normalizando d2
return (d2(d1,vol,T)-xbarra)/dp
def C(Nd1, S0, Nd2, X, r, T):
return Nd1(d1, xbarra, dp)*S0 -Nd2(d2, d1, xbarra,dp)*X*math.exp(-r*T)
Substituindo para os valores dados na questão (os valores estão em US$ milhões):
if __name__=='__main__':
#BlackScholes Parameters
S0=300 #fluxo de caixa futuro descontado (DFC)
X=300 #custo inicial
vol=0.3 #fator de volatilidade anual
r=0.05 #taxa de juros contínua
T=5 #tempo do projeto
xbarra=-1.277
dp=2.557
Resultado:
print(C(Nd1, S0, Nd2, X, r, T))
112.8121315637216
Neste caso, a equação principal é dada por:
\( S_t = S_{t-1} -S_{t-1}(r\delta_t + \sigma \epsilon \sqrt{\delta_t}) \)
Onde \( S_t\) é o valor corrente do ativo no tempo t ( \( S_0 \) é o valor no tempo zero, ou seja, é o fluxo de caixa futuro descontado (DFC)), \( \sigma \) é o fator de volatilidade anual, \( r \) é a taxa de juros contínua, \( \delta_t \) é o fator de incremento no tempo e \( \epsilon \sim N(0,1) \) que será simulado. Para simular \( \epsilon \) criei uma matriz com 1000 linhas, onde cada linha representa uma das 1000 simulações do modelo de Monte Carlo (os valores estão em US$ milhões).
#MonteCarlo Parameters
S0=300 #DFC
X=300 #custo inicial
vol=0.3 #fator de volatilidade anual
r=0.05 #taxa de juros contínua
T=5 #tempo do projeto
ts_bs=0.5 #incremento no tempo
e1=np.random.normal(0, 1, 1000)
e2=np.random.normal(0, 1, 1000)
e3=np.random.normal(0, 1, 1000)
e4=np.random.normal(0, 1, 1000)
e5=np.random.normal(0, 1, 1000)
colunas = ['S0','S1','S2','S3','S4','S5']
df=pd.DataFrame(columns=colunas, index=list(range(1000)))
df['S0']=np.ones((1000,1))*S0
df['S1']=df['S0'] + df['S0']*(r*ts_bs + vol*e1*((ts_bs*1)**(1/2))) #T=1
df['S2']=df['S1'] + df['S1']*(r*ts_bs + vol*e2*((ts_bs*2)**(1/2))) #T=2
df['S3']=df['S2'] + df['S2']*(r*ts_bs + vol*e3*((ts_bs*3)**(1/2))) #T=3
df['S4']=df['S3'] + df['S3']*(r*ts_bs + vol*e4*((ts_bs*4)**(1/2))) #T=4
df['S5']=df['S4'] + df['S4']*(r*ts_bs + vol*e5*((ts_bs*5)**(1/2))) #T=5
df['S']=df['S5']-X #Descontar do preço inicial
df['S'][df['S'] < 0] = 0 #Retirar os valores negativos
df['SS']=df['S']/ math.exp(T*r) #Trazer para o valor presente
print (df.head())
S0 S1 S2 S3
0 300.0 293.113987 274.047409 273.703754
1 300.0 371.690300 345.183854 330.803967
2 300.0 365.988123 402.431531 529.311190
3 300.0 226.522438 258.380246 185.515708
4 300.0 323.448469 586.806771 845.507553
S4 S5 S SS
0 352.657403 410.919094 110.919094 86.383877
1 443.009509 223.060508 0.000000 0.000000
2 655.303765 837.086725 537.086725 418.283562
3 167.183016 330.285038 30.285038 23.586011
4 104.210120 112.377667 0.000000 0.000000
Resultado:
print (df['SS'].mean())
110.58149678565597
Agora, no método binomial, é preciso calcular todas as bifurcações da árvore de probabilidade, e depois fazer o método de indução reversa para achar o ROV. As probabilidades de subida e descida de cada estado, é respectivamente:
\( up= e^{\sigma\sqrt{\delta_t}} \)
\( down =\frac{1}{up}=\frac{1}{e^{\sigma\sqrt{\delta_t}}} \)
Onde \( \sigma \) é o fator de volatilidade anual e \( \delta_t \) é o fator de incremento no tempo.
Para resolver o problema de forma geral, pode-se fazer uma função que calcula cada árvore de probabilidade, como fiz abaixo:
def binom_tree_call (T,N,t, r, vol, r, X, array_out=False):
dt=T/N
u=np.exp(vol*(dt)**(1/2))
u=1/u
p=(np.exp(r*dt)-d/(u-d)
price_tree = np.zeros([N+1,N+1])
for i in range (N+1):
for j in range (i+1):
price_tree[j,i]=S0*(d**j)*(u**(i-j))
option = np.zeros([N+1,N+1])
option[:,N]=np.maximum(np.zeros(N+1),price_tree[:,N]-X)
for i in np.arrange(N-1, -1, -1):
for j in np.arrange(0,i+1):
option[j,i]=np.exp(-r*dt)*(p*option[j,i+1]+(1-p)*option[j+i,i+1])
if array_out:
return [option[0,0], price_tree, option]
else:
return option[0,0]
#Binomial Parameters
S0=300 #DFC
X=300 #custo inicial
vol=0.3 #fator de volatilidade anual
r=0.05 #taxa de juros contínua
T=5 #tempo do projeto
N=5
ts_bm=1 #incremento no tempo
Ou ainda, se preferir pode calcular o valor de cada período passo a passo:
import math
def u(t):
return math.exp(vol*(ts_bm*t)**(1/2))
def d(t):
return 1/u(t)
def p(t):
return (math.exp(r*ts_bm*t)-d(t))/(u(t)-d(t))
S0u = S0*u(1)
S0d = S0*d(1)
S1= [S0u,S0d]
S0uu=S0u*u(1)
S0ud=S0u*d(1)
S0dd=S0d*d(1)
S2 = [S0uu, S0ud, S0dd]
S0uuu=S0uu*u(1)
S0uud=S0ud*u(1)
S0udd=S0ud*d(1)
S0ddd=S0dd*d(1)
S3 = [S0uuu, S0uud, S0udd, S0ddd]
S0uuuu=S0uuu*u(1)
S0uuud=S0uuu*d(1)
S0uudd=S0uud*d(1)
S0uddd=S0udd*d(1)
S0dddd=S0ddd*d(1)
S4 = [S0uuuu,S0uuud,S0uudd,S0uddd,S0dddd]
S0uuuuu=S0uuuu*u(1)
S0uuuud=S0uuuu*d(1)
S0uuudd=S0uuud*d(1)
S0uuddd=S0uudd*d(1)
S0udddd=S0uddd*d(1)
S0ddddd=S0dddd*d(1)
S5 = [S0uuuuu,S0uuuud,S0uuudd,S0uuddd,S0udddd, S0ddddd]
#Indução Reversa
I5 = [S0uuuuu - X,S0uuuud -X,S0uuudd - X,S0uuddd-X,S0udddd-X, S0ddddd-X]
def replace_neg(old_list):
new_list = []
for i in old_list:
if i < 0:
new_list.append(0)
else:
new_list.append(i)
return new_list
I5 = replace_neg(I5)
I0uuuuu=I5[0]
I0uuuud=I5[1]
I0uuudd=I5[2]
I0uuddd=I5[3]
I0udddd=I5[4]
I0ddddd=I5[5]
I0uuuu=(p(1)*I0uuuuu + (1-p(1))*I0uuuud)*math.exp(-r*ts_bm*1)
I0uuud=(p(1)*I0uuuud + (1-p(1))*I0uuudd)*math.exp(-r*ts_bm*1)
I0uudd=(p(1)*I0uuudd + (1-p(1))*I0uuddd)*math.exp(-r*ts_bm*1)
I0uddd=(p(1)*I0uuddd + (1-p(1))*I0udddd)*math.exp(-r*ts_bm*1)
I0dddd=(p(1)*I0udddd + (1-p(1))*I0ddddd)*math.exp(-r*ts_bm*1)
I4 = [I0uuuu, I0uuud, I0uudd, I0uddd, I0dddd]
I0uuu=(p(1)*I0uuuu + (1-p(1))*I0uuud)*math.exp(-r*ts_bm*1)
I0uud=(p(1)*I0uuud + (1-p(1))*I0uudd)*math.exp(-r*ts_bm*1)
I0udd=(p(1)*I0uudd + (1-p(1))*I0uddd)*math.exp(-r*ts_bm*1)
I0ddd=(p(1)*I0uddd + (1-p(1))*I0dddd)*math.exp(-r*ts_bm*1)
I3= [I0uuu, I0uud, I0udd, I0ddd]
I0uu=(p(1)*I0uuu + (1-p(1))*I0uud)*math.exp(-r*ts_bm*1)
I0ud=(p(1)*I0uud + (1-p(1))*I0udd)*math.exp(-r*ts_bm*1)
I0dd=(p(1)*I0udd + (1-p(1))*I0ddd)*math.exp(-r*ts_bm*1)
I2 = [I0uu, I0ud, I0dd]
I0u = (p(1)*I0uu + (1-p(1))*I0ud)*math.exp(-r*ts_bm*1)
I0d = (p(1)*I0ud + (1-p(1))*I0dd)*math.exp(-r*ts_bm*1)
I1 = [I0u, I0d]
I0 = (p(1)*I0u + (1-p(1))*I0d)*math.exp(-r*ts_bm*1)
Resultado:
print (binom_tree_call (T,N,t, r, vol, r, X, array_out=False))
print (I0)
110.45581453294359
Assim, os valores reais das opções são:
- Black-Scholes: 112 milhões
- Monte Carlo: 110 milhões
- Binomial: 110 milhões
Assim, fica fácil de observar que os ROV's dos três métodos são muito parecidos. A equação de Black-Scholes chega em um resultado muito aproximado, através de uma equação bem simplificada. Porém, os outros métodos podem nos oferecer uma análise mais detalhada da opção. Na simulação de Monte Carlo podemos ver vários cenários possíveis de qual seria o valor da opção. No método Binomial também é possível ver os cenários extremamentes positivos e negativos do valor real da opcão.