El tema principal en este módulo fueron simulaciones Montecarlo. Al finalizar este módulo, se espera que ustedes tengan las siguientes competencias
- Evaluar integrales (o encontrar áreas) numéricamente mendiante métodos Montecarlo.
- Poder replicar fractales aleatorios símples (como los de Barnsley), dadas las características del mismo.
- Realizar evaluaciones de probabilidad precio-umbral.
En clase desarrollamos funciones para la evaluación de integrales con ambos métodos (explicar porqué la segunda se puede ver como una integral). Las funciones son las siguientes:
In [21]:
# Importamos librerías
import numpy as np
import pandas as pd
import random
In [22]:
def int_montecarlo1(f, a, b, N):
# Evaluación numérica de integrales por Montecarlo tipo 1
# f=f(x) es la función a integrar (debe ser declarada previamente) que devuelve para cada x su valor imagen,
# a y b son los límites inferior y superior del intervalo donde se integrará la función, y N es el número
# de puntos con que se aproximará.
return (b-a)/N*np.sum(f(np.random.uniform(a, b, N)))
In [23]:
def int_montecarlo2(region, a1, b1, a2, b2, N):
# Evaluación numérica de integrales por Montecarlo tipo 2
# region=region(x,y) retorna True si la coordenada (x,y) pertenece a la región a integrar y False de lo
# contrario , a1, b1, a2, b2 son los límites del rectángulo que contiene la región, y N es el número de
# puntos con que se aproximará.
A_R = (b1-a1)*(b2-a2)
x = np.random.uniform(a1, b1, N.astype(int))
y = np.random.uniform(a2, b2, N.astype(int))
return A_R*np.sum(region(x,y))/N
Considere las funciones $f_1(x)=\sqrt{1+x^{4}}$, $f_2(x)=\ln(\ln x)$, $f_3(x)=\frac {1}{\ln x}$, $f_4(x)=e^{e^{x}}$, $f_5(x)=e^{-{\frac {x^{2}}{2}}}$ y $f_6(x)=\sin(x^{2})$.
Utilizar las funciones anteriores para realizar la evaluación numérica de las integrales de las funciones anteriores en el intervalo $(4,5)$. Poner los resultados en una tabla cuyas filas correspondan a la cantidad de términos utilizados en la aproximación (usar 10, 100, 1000, 10000 y 100000 términos) y cuyas columnas correspondan a las funciones.
Hacer una tabla por cada método.
¿Se pueden ver diferencias notables en la velocidad de convergencia de los métodos?
In [24]:
def f1(x):
return np.sqrt(x**4+1)
def f2(x):
return np.log(np.log(x))
def f3(x):
return 1/np.log(x)
def f4(x):
return np.exp(np.exp(x))
def f5(x):
return np.exp(-x**2/2)
def f6(x):
return np.sin(x**2)
In [25]:
a, b = 4, 5
N = np.logspace(1,5,5)
tabla_tipo1 = pd.DataFrame(columns=['$f_1(x)$','$f_2(x)$','$f_3(x)$','$f_4(x)$','$f_5(x)$','$f_6(x)$'], index=N)
tabla_tipo1.index.name = 'Cantidad de términos'
for n in N:
n = n.astype('int')
tabla_tipo1.loc[n, '$f_1(x)$'] = int_montecarlo1(f1, a, b, n)
tabla_tipo1.loc[n, '$f_2(x)$'] = int_montecarlo1(f2, a, b, n)
tabla_tipo1.loc[n, '$f_3(x)$'] = int_montecarlo1(f3, a, b, n)
tabla_tipo1.loc[n, '$f_4(x)$'] = int_montecarlo1(f4, a, b, n)
tabla_tipo1.loc[n, '$f_5(x)$'] = int_montecarlo1(f5, a, b, n)
tabla_tipo1.loc[n, '$f_6(x)$'] = int_montecarlo1(f6, a, b, n)
tabla_tipo1
Out[25]:
In [26]:
def region1(x,y):
return y<=f1(x)
def region2(x,y):
return y<=f2(x)
def region3(x,y):
return y<=f3(x)
def region4(x,y):
return y<=f4(x)
def region5(x,y):
return y<=f5(x)
In [27]:
a1, b1 = 4, 5
a2 = 0
b2 = [26, 0.5, 0.8, 2.86e+64, 0.00036]
N = np.logspace(1,5,5)
tabla_tipo2 = pd.DataFrame(columns=['$f_1(x)$','$f_2(x)$','$f_3(x)$','$f_4(x)$','$f_5(x)$'], index=N)
tabla_tipo2.index.name = 'Cantidad de términos'
for n in N:
n = n.astype('int')
tabla_tipo2.loc[n, '$f_1(x)$'] = int_montecarlo2(region1, a1, b1, a2, b2[0], n)
tabla_tipo2.loc[n, '$f_2(x)$'] = int_montecarlo2(region2, a1, b1, a2, b2[1], n)
tabla_tipo2.loc[n, '$f_3(x)$'] = int_montecarlo2(region3, a1, b1, a2, b2[2], n)
tabla_tipo2.loc[n, '$f_4(x)$'] = int_montecarlo2(region4, a1, b1, a2, b2[3], n)
tabla_tipo2.loc[n, '$f_5(x)$'] = int_montecarlo2(region5, a1, b1, a2, b2[4], n)
tabla_tipo2
Out[27]:
In [28]:
tabla_tipo1
Out[28]:
Pues bien, usando la misma idea de transformaciones afines que se escogen con cierta probabilidad, se pueden generar una infinidad inimaginable de fractales. Incluso, se pueden generar fractales aleatorios que poseen un atractor determinístico (¿Qué es esto?).
Como en la clase de fractales, repliquemos el fractal tipo Barnsley descrito por la siguiente tabla...
Referencia:
In [53]:
import pandas as pd
import numpy as np
In [54]:
i = np.arange(4)
df = pd.DataFrame(index=i,columns=['$a_i$', '$b_i$', '$c_i$', '$d_i$', '$e_i$', '$f_i$', '$p_i$'], dtype='float')
df.index.name = "$i$"
df['$a_i$'] = [0.5, 0.5, 0.5, 0.5]
df['$b_i$'] = [0.0, 0.0, 0.0, 0.0]
df['$c_i$'] = [0.0, 0.0, 0.0, 0.0]
df['$d_i$'] = [0.5, 0.5, 0.5, 0.5]
df['$e_i$'] = [1.0, 50.0, 1.0, 50.0]
df['$f_i$'] = [1.0, 1.0, 50.0, 50.0]
df['$p_i$'] = [0.1, 0.2, 0.3, 0.4]
df.round(2)
Out[54]:
In [55]:
mat_barnsley = np.array([[0.5, 0.0, 0.0, 0.5, 1.0, 1.0, 0.1],
[0.5, 0.0, 0.0, 0.5, 50.0, 1.0, 0.2],
[0.5, 0.0, 0.0, 0.5, 1.0, 50.0, 0.3],
[0.5, 0.0, 0.0, 0.5, 50.0, 50.0, 0.4]])
mat_barnsley
Out[55]:
In [56]:
x, y = [0.0], [0.0]
In [62]:
x, y = [0.0], [0.0]
for k in range(100000):
p = random.random()
if p <= mat_barnsley[0,6]:
i = 0
elif p <= mat_barnsley[0,6]+mat_barnsley[1,6]:
i = 1
elif p <= mat_barnsley[0,6]+mat_barnsley[1,6]+mat_barnsley[2,6]:
i = 2
else:
i = 3
x.append(mat_barnsley[i,0]*x[-1]+mat_barnsley[i,1]*y[-1]+mat_barnsley[i,4])
y.append(mat_barnsley[i,2]*x[-2]+mat_barnsley[i,3]*y[-1]+mat_barnsley[i,5])
plt.figure(figsize=(10,10))
plt.scatter(x,y, c='k', s = 0.2)
plt.show()
En las últimas clases vimos una aplicación de simulación montecarlo. Consistía en descargar datos históricos de precio de cierre de acciones de alguna compañía, proyectar esos precios y sacar la probabilidad de que los precios en el siguiente año sobrepasaran cierto precio umbral.
En este ejemplo evaluaremos dos compañías con tendencias más o menos similares (Apple y Microsoft) veremos cuál tiene más probabilidades de darnos un interés deseado.
Además, descargaremos los datos del presente año (2017) para ver si el análisis concuerda.
Las funciones que desarrollamos en la clase del martes son:
In [111]:
# Importamos librerías
import numpy as np
import pandas as pd
import pandas_datareader as data
import matplotlib.pyplot as plt
%matplotlib inline
In [113]:
# Función que creamos en clase para bajar datos
def load_adj_close(ticker, data_source, start_date, end_date):
panel_data = data.DataReader(ticker, data_source, start_date, end_date)
closes = panel_data.loc['Adj Close']
all_weekdays = pd.date_range(start=start_date, end=end_date)
closes = closes.reindex(all_weekdays)
closes = closes.fillna(method='ffill')
return closes
# Función que devuelve rendimientos diarios, media y desviación estándar
def mu_std_daily_ret(closes):
daily_returns = (np.log(closes/closes.shift(1)))[1:]
mu = daily_returns.mean().values
sigma = daily_returns.std().values
return daily_returns, mu, sigma
# Función que simula varios escenarios de rendimientos diarios
def daily_ret_sim(mu, sigma, ndays, ntraj, start_date):
dates = pd.date_range(start=start_date,periods=ndays)
return pd.DataFrame(sigma*np.random.randn(ndays, ntraj)+mu, index = dates)
# Función de proyección de precios
def closes_proj(simret, closes):
return (closes.iloc[-1,:].values[0])*np.exp(simret.cumsum())
Descarguemos datos para Apple y Microsoft en el 2016
In [114]:
# Descargamos datos de microsoft en el 2016
ticker = ['AAPL','MSFT']
data_source = 'yahoo'
start_date = '2016-01-04'
end_date = '2016-12-31'
closes = load_adj_close(ticker, data_source, start_date, end_date)
In [115]:
# Grafiquemos
closes.plot(figsize=(8,6));
Calculamos los rendimientos diarios junto con sus características estadísticas
In [116]:
# Calculamos con la función anterior
daily_returns, mu, sigma = mu_std_daily_ret(closes)
mu, sigma
Out[116]:
In [117]:
# Son similares
daily_returns.plot(figsize=(8,6));
In [118]:
mu_AAPL, sigma_AAPL = mu[0], sigma[0]
mu_MSFT, sigma_MSFT = mu[1], sigma[1]
Simulamos 100 escenarios de rendimientos diarios para el 2017 (para cada una de las empresas)
In [119]:
# Simulamos 100 escenarios para todo el 2017
ndays = 360
ntraj = 100
start_date = '2017-01-01'
simret_AAPL = daily_ret_sim(mu_AAPL, sigma_AAPL, ndays, ntraj, start_date)
simret_MSFT = daily_ret_sim(mu_MSFT, sigma_MSFT, ndays, ntraj, start_date)
Calculamos los precios con base en los rendimientos simulados
In [120]:
# Proyección de precios y concatenación con precios de 2016
simdata_AAPL = closes_proj(simret_AAPL, pd.DataFrame({'AAPL':closes['AAPL']}, index = closes.index))
simdata_MSFT = closes_proj(simret_MSFT, pd.DataFrame({'MSFT':closes['MSFT']}, index = closes.index))
In [121]:
simdata_AAPL.plot(figsize=(8,6), legend=False);
In [122]:
simdata_MSFT.plot(figsize=(8,6), legend=False);
Calculamos las probabilidades con base en una tasa de interés anual deseada
In [62]:
# Cálculo de umbral con base en una tasa de interés anual deseada
tasa_anual = 0.1
K_AAPL, K_MSFT = (1+tasa_anual)*closes.iloc[-1,:]
In [137]:
# Función para calcular probabilidades diarias de que el precio sea mayor que un umbral
def daily_prob(simdata, K):
strike = pd.DataFrame(K*np.ones(simdata.shape).reshape(simdata.shape), index = simdata.index)
return (simdata>strike).T.sum()/simdata.shape[1]
In [138]:
prob = pd.DataFrame({'AAPL':daily_prob(simdata_AAPL, K_AAPL), 'MSFT': daily_prob(simdata_MSFT, K_MSFT)}, index = simdata_AAPL.index)
prob.plot(figsize=(8,6));
In [139]:
prob
Out[139]:
Finalmente, veamos los datos reales del 2017 para ver que tan acertados fueron nuestros análisis...
In [163]:
ticker = ['AAPL','MSFT']
data_source = 'yahoo'
start_date = '2017-01-03'
end_date = '2017-10-15'
closes_actuales = load_adj_close(ticker, data_source, start_date, end_date)
umbral = pd.DataFrame({'Umbral AAPL':K_AAPL*np.ones(closes_actuales.shape[0]), 'Umbral MSFT':K_MSFT*np.ones(closes_actuales.shape[0])}, index= closes_actuales.index)
comparacion = pd.concat((closes_actuales.T, umbral.T)).T
comparacion.plot(figsize=(8,6));