Algoritmos Supervisados: Predicción

En esta sesión seguimos profundizando en las técnicas de Machine Learning.

A lo largo de este cuaderno, usaremos datos de cotización de varias compañias para predecir valores de cierre futuros. Para predecir usaremos una técnica más ``clásica’’ como el método ARIMA y, posteriormente, aplicaremos técnicas propias de ML.

Existen númerosos algoritmos supervisados para la predicción de valores, destacamos:

La gran mayoría están implementados en la librería de python scikit-learn/supervised_learning

|5028afae376a46caa1542fed3bad2e74|

Fuente de la imagen: [5] Machine Learning & Data Science Blueprints for Finance

[1]:
%pip install seaborn
Requirement already satisfied: seaborn in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (0.12.0)
Requirement already satisfied: numpy>=1.17 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from seaborn) (1.22.3)
Requirement already satisfied: pandas>=0.25 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from seaborn) (1.4.2)
Requirement already satisfied: matplotlib>=3.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from seaborn) (3.5.2)
Requirement already satisfied: cycler>=0.10 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from matplotlib>=3.1->seaborn) (0.11.0)
Requirement already satisfied: fonttools>=4.22.0 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from matplotlib>=3.1->seaborn) (4.33.3)
Requirement already satisfied: kiwisolver>=1.0.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from matplotlib>=3.1->seaborn) (1.4.2)
Requirement already satisfied: packaging>=20.0 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from matplotlib>=3.1->seaborn) (21.3)
Requirement already satisfied: pillow>=6.2.0 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from matplotlib>=3.1->seaborn) (9.1.1)
Requirement already satisfied: pyparsing>=2.2.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from matplotlib>=3.1->seaborn) (2.4.7)
Requirement already satisfied: python-dateutil>=2.7 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from matplotlib>=3.1->seaborn) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas>=0.25->seaborn) (2022.6)
Requirement already satisfied: six>=1.5 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from python-dateutil>=2.7->matplotlib>=3.1->seaborn) (1.16.0)

[notice] A new release of pip is available: 23.2.1 -> 23.3.1
[notice] To update, run: pip install --upgrade pip
Note: you may need to restart the kernel to use updated packages.

Obtención de datos

¡Necesitamos datos! Principalmente, datos con índices de tiempo: series temporales. Por ello, vamos a instalar una librería que nos proporcionará datos públicos ya preparados.

[2]:
%pip install pandas_datareader
%pip install yfinance

Requirement already satisfied: pandas_datareader in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (0.10.0)
Requirement already satisfied: lxml in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas_datareader) (4.9.1)
Requirement already satisfied: pandas>=0.23 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas_datareader) (1.4.2)
Requirement already satisfied: requests>=2.19.0 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas_datareader) (2.31.0)
Requirement already satisfied: python-dateutil>=2.8.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas>=0.23->pandas_datareader) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas>=0.23->pandas_datareader) (2022.6)
Requirement already satisfied: numpy>=1.20.0 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas>=0.23->pandas_datareader) (1.22.3)
Requirement already satisfied: charset-normalizer<4,>=2 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.19.0->pandas_datareader) (2.1.1)
Requirement already satisfied: idna<4,>=2.5 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.19.0->pandas_datareader) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.19.0->pandas_datareader) (1.26.6)
Requirement already satisfied: certifi>=2017.4.17 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.19.0->pandas_datareader) (2023.5.7)
Requirement already satisfied: six>=1.5 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from python-dateutil>=2.8.1->pandas>=0.23->pandas_datareader) (1.16.0)

[notice] A new release of pip is available: 23.2.1 -> 23.3.1
[notice] To update, run: pip install --upgrade pip
Note: you may need to restart the kernel to use updated packages.
Requirement already satisfied: yfinance in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (0.2.33)
Requirement already satisfied: pandas>=1.3.0 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (1.4.2)
Requirement already satisfied: numpy>=1.16.5 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (1.22.3)
Requirement already satisfied: requests>=2.31 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (2.31.0)
Requirement already satisfied: multitasking>=0.0.7 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (0.0.11)
Requirement already satisfied: lxml>=4.9.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (4.9.1)
Requirement already satisfied: appdirs>=1.4.4 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (1.4.4)
Requirement already satisfied: pytz>=2022.5 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (2022.6)
Requirement already satisfied: frozendict>=2.3.4 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (2.3.10)
Requirement already satisfied: peewee>=3.16.2 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (3.17.0)
Requirement already satisfied: beautifulsoup4>=4.11.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (4.11.1)
Requirement already satisfied: html5lib>=1.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from yfinance) (1.1)
Requirement already satisfied: soupsieve>1.2 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from beautifulsoup4>=4.11.1->yfinance) (2.3.2.post1)
Requirement already satisfied: six>=1.9 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from html5lib>=1.1->yfinance) (1.16.0)
Requirement already satisfied: webencodings in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from html5lib>=1.1->yfinance) (0.5.1)
Requirement already satisfied: python-dateutil>=2.8.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas>=1.3.0->yfinance) (2.8.2)
Requirement already satisfied: charset-normalizer<4,>=2 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.31->yfinance) (2.1.1)
Requirement already satisfied: idna<4,>=2.5 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.31->yfinance) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.31->yfinance) (1.26.6)
Requirement already satisfied: certifi>=2017.4.17 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from requests>=2.31->yfinance) (2023.5.7)

[notice] A new release of pip is available: 23.2.1 -> 23.3.1
[notice] To update, run: pip install --upgrade pip
Note: you may need to restart the kernel to use updated packages.
[3]:
# código para corregir divergencias entre versiones de la librería pandas_Datareader
import yfinance as yf
yf.pdr_override()
[12]:
import pandas as pd
from pandas_datareader import data as pdr

data = pdr.get_data_yahoo("GOOGL", start="2017-01-01", end="2020-12-31")

#data = yf.Ticker("GOOGL") #, start="2017-01-01", end="2020-12-31")
print(type(data))
data.to_csv("GOOGL.csv")
[*********************100%%**********************]  1 of 1 completed
<class 'pandas.core.frame.DataFrame'>
[14]:
data = pd.read_csv("GOOGL.csv")
[19]:
data.index = pd.DatetimeIndex(data["Date"])
data
[19]:
Date Open High Low Close Adj Close Volume
Date
2017-01-03 2017-01-03 40.030998 40.571999 39.844501 40.400501 40.400501 39180000
2017-01-04 2017-01-04 40.494499 40.671501 40.205502 40.388500 40.388500 30306000
2017-01-05 2017-01-05 40.375000 40.687000 40.296001 40.651001 40.651001 26810000
2017-01-06 2017-01-06 40.749500 41.448002 40.575001 41.260502 41.260502 40342000
2017-01-09 2017-01-09 41.318501 41.521500 41.081001 41.359001 41.359001 28178000
... ... ... ... ... ... ... ...
2020-12-23 2020-12-23 86.196503 87.205498 86.059998 86.411499 86.411499 22974000
2020-12-24 2020-12-24 86.449997 87.120499 86.217499 86.708000 86.708000 9312000
2020-12-28 2020-12-28 87.245499 89.349998 87.091003 88.697998 88.697998 27650000
2020-12-29 2020-12-29 89.361504 89.423500 87.755501 87.888000 87.888000 19726000
2020-12-30 2020-12-30 88.250000 88.388000 86.400002 86.812500 86.812500 21026000

1006 rows × 7 columns

[20]:
type(data.index[0])
[20]:
pandas._libs.tslibs.timestamps.Timestamp
[32]:
print(data) # Python tip: una invocación stk_data.__rep__
                 Open       High        Low      Close  Adj Close    Volume
Date
2017-01-03  40.030998  40.571999  39.844501  40.400501  40.400501  39180000
2017-01-04  40.494499  40.671501  40.205502  40.388500  40.388500  30306000
2017-01-05  40.375000  40.687000  40.296001  40.651001  40.651001  26810000
2017-01-06  40.749500  41.448002  40.575001  41.260502  41.260502  40342000
2017-01-09  41.318501  41.521500  41.081001  41.359001  41.359001  28178000
...               ...        ...        ...        ...        ...       ...
2020-12-23  86.196503  87.205498  86.059998  86.411499  86.411499  22974000
2020-12-24  86.449997  87.120499  86.217499  86.708000  86.708000   9312000
2020-12-28  87.245499  89.349998  87.091003  88.697998  88.697998  27650000
2020-12-29  89.361504  89.423500  87.755501  87.888000  87.888000  19726000
2020-12-30  88.250000  88.388000  86.400002  86.812500  86.812500  21026000

[1006 rows x 6 columns]
[21]:
data.columns # ¿Cómo es el índice de acceso a los datos?
[21]:
Index(['Date', 'Open', 'High', 'Low', 'Close', 'Adj Close', 'Volume'], dtype='object')
[22]:
data["Close"] #Hagamos una prueba
[22]:
Date
2017-01-03    40.400501
2017-01-04    40.388500
2017-01-05    40.651001
2017-01-06    41.260502
2017-01-09    41.359001
                ...
2020-12-23    86.411499
2020-12-24    86.708000
2020-12-28    88.697998
2020-12-29    87.888000
2020-12-30    86.812500
Name: Close, Length: 1006, dtype: float64
[23]:
data.index
[23]:
DatetimeIndex(['2017-01-03', '2017-01-04', '2017-01-05', '2017-01-06',
               '2017-01-09', '2017-01-10', '2017-01-11', '2017-01-12',
               '2017-01-13', '2017-01-17',
               ...
               '2020-12-16', '2020-12-17', '2020-12-18', '2020-12-21',
               '2020-12-22', '2020-12-23', '2020-12-24', '2020-12-28',
               '2020-12-29', '2020-12-30'],
              dtype='datetime64[ns]', name='Date', length=1006, freq=None)

Análisis de la serie temporal

Existen númerosos indicadores sobre el comportamiento de una serie temporal. La descomposición permite tener una idea de su “naturaleza”. Hay dos tipos de descomposiciones: aditiva o multiplicativa. - Descomposición aditiva \(Y[t] = T[t] + S[t] + e[t]\) - T_{t}, the trend component at time t, which reflects the long-term progression of the series (secular variation). A trend exists when there is a persistent increasing or decreasing direction in the data. The trend component does not have to be linear. - C_{t}, the cyclical component at time t, which reflects repeated but non-periodic fluctuations. The duration of these fluctuations depend on the nature of the time series. - S_{t}, the seasonal component at time t, reflecting seasonality (seasonal variation). A seasonal pattern exists when a time series is influenced by seasonal factors. Seasonality occurs over a fixed and known period (e.g., the quarter of the year, the month, or day of the week).[1] - I_t, the irregular component (or “noise”) at time t, which describes random, irregular influences. It represents the residuals or remainder of the time series after the other components have been removed.

Necesitaremos una librería especial para facilitar esta tarea.

[24]:
%pip install statsmodels
Requirement already satisfied: statsmodels in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (0.13.2)
Requirement already satisfied: numpy>=1.17 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from statsmodels) (1.22.3)
Requirement already satisfied: scipy>=1.3 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from statsmodels) (1.8.1)
Requirement already satisfied: pandas>=0.25 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from statsmodels) (1.4.2)
Requirement already satisfied: patsy>=0.5.2 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from statsmodels) (0.5.2)
Requirement already satisfied: packaging>=21.3 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from statsmodels) (21.3)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from packaging>=21.3->statsmodels) (2.4.7)
Requirement already satisfied: python-dateutil>=2.8.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas>=0.25->statsmodels) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from pandas>=0.25->statsmodels) (2022.6)
Requirement already satisfied: six in /Users/isaac/.pyenv/versions/3.9.7/envs/my397/lib/python3.9/site-packages (from patsy>=0.5.2->statsmodels) (1.16.0)

[notice] A new release of pip is available: 23.2.1 -> 23.3.1
[notice] To update, run: pip install --upgrade pip
Note: you may need to restart the kernel to use updated packages.
[25]:
# Paso 1
import matplotlib.pyplot as plt
import statsmodels.api as sm
# https://www.statsmodels.org/dev/generated/statsmodels.tsa.seasonal.seasonal_decompose.html

res = sm.tsa.seasonal_decompose(data['Close'],period=30)
#res.plot()

fig, axes = plt.subplots(ncols=1, nrows=4, sharex=True, figsize=(12,8))
res.observed.plot(ax=axes[0], legend=False)
axes[0].set_ylabel('Observed')
res.trend.plot(ax=axes[1], legend=False)
axes[1].set_ylabel('Trend')
res.seasonal.plot(ax=axes[2], legend=False)
axes[2].set_ylabel('Seasonal')
res.resid.plot(ax=axes[3], legend=False)
axes[3].set_ylabel('Residual')
plt.show()
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_15_0.png

Predicción con ARIMA

[26]:
# Podemos trabajar con diferentes periodos para la predicción.
# Vamos a probar meses
data.index = pd.to_datetime(data.index)
month_index = data.index.to_period('M')
month_index
[26]:
PeriodIndex(['2017-01', '2017-01', '2017-01', '2017-01', '2017-01', '2017-01',
             '2017-01', '2017-01', '2017-01', '2017-01',
             ...
             '2020-12', '2020-12', '2020-12', '2020-12', '2020-12', '2020-12',
             '2020-12', '2020-12', '2020-12', '2020-12'],
            dtype='period[M]', name='Date', length=1006)
[27]:
# Paso 2.a  Cargamos/Fit/Entrenamos ARIMA
from statsmodels.tsa.arima_model import ARIMA
# https://www.statsmodels.org/dev/generated/statsmodels.tsa.arima.model.ARIMA.html
close = data["Close"]
close.index = month_index
print(close.shape)
print(close[-12:])

(1006,)
Date
2020-12    87.612999
2020-12    88.054001
2020-12    87.859497
2020-12    87.025497
2020-12    86.310997
2020-12    86.727997
2020-12    86.011002
2020-12    86.411499
2020-12    86.708000
2020-12    88.697998
2020-12    87.888000
2020-12    86.812500
Freq: M, Name: Close, dtype: float64
[29]:
mod = sm.tsa.arima.ARIMA(endog=close, order=(3, 1, 0))
#U n 3 por autoregression, 1 para hacer la serie estacionaria, y una media móvil de 0.
modfit = mod.fit()
print(modfit.summary())

                               SARIMAX Results
==============================================================================
Dep. Variable:                  Close   No. Observations:                 1006
Model:                 ARIMA(3, 1, 0)   Log Likelihood               -1485.452
Date:                Fri, 15 Dec 2023   AIC                           2978.904
Time:                        17:39:09   BIC                           2998.554
Sample:                    01-31-2017   HQIC                          2986.370
                         - 12-31-2020
Covariance Type:                  opg
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.1338      0.020     -6.863      0.000      -0.172      -0.096
ar.L2          0.0155      0.020      0.777      0.437      -0.024       0.055
ar.L3          0.0296      0.022      1.375      0.169      -0.013       0.072
sigma2         1.1255      0.028     39.852      0.000       1.070       1.181
===================================================================================
Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):              1273.25
Prob(Q):                              0.95   Prob(JB):                         0.00
Heteroskedasticity (H):               4.60   Skew:                            -0.41
Prob(H) (two-sided):                  0.00   Kurtosis:                         8.45
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[30]:
residuals = pd.DataFrame(modfit.resid)
residuals.plot()
plt.show()
print(residuals.describe()) # más o menos valores cercanos al 0 ok
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_20_0.png
                 0
count  1006.000000
mean      0.090508
std       1.655711
min      -6.204212
25%      -0.357341
50%       0.092804
75%       0.539872
max      40.400501
[31]:
# Paso 2.b Predecimos con ARIMA
thres = 0.7
split_samples = int(len(close)*thres)
train, test = close[:split_samples],close[split_samples:]
print(train.shape,test.shape)
(704,) (302,)
[36]:
import math
from sklearn.metrics import mean_squared_error

predictions = []
history = train.to_list()

for i,y in enumerate(test):
    model = sm.tsa.arima.ARIMA(endog=history, order=(3, 1, 0))
    modelfit = model.fit()
    ypred = modelfit.forecast()
    predictions.append(ypred)
    history.append(y)

    if i%20==0: # Cada 20 veces mostramos este mensaje
        print("%i\t Prediction: %.2f - Expected=%.2f"%(i,ypred,y))


rmse = math.sqrt(mean_squared_error(test,predictions))
print("RMSE: %.4f"%rmse)
0        Prediction: 62.17 - Expected=62.21
20       Prediction: 66.63 - Expected=65.99
40       Prediction: 68.06 - Expected=67.74
60       Prediction: 72.03 - Expected=72.51
80       Prediction: 75.64 - Expected=75.94
100      Prediction: 60.50 - Expected=53.65
120      Prediction: 60.56 - Expected=63.26
140      Prediction: 70.11 - Expected=68.76
160      Prediction: 72.65 - Expected=73.24
180      Prediction: 75.14 - Expected=75.93
200      Prediction: 73.90 - Expected=75.25
220      Prediction: 85.37 - Expected=81.48
240      Prediction: 74.24 - Expected=72.78
260      Prediction: 77.42 - Expected=80.81
280      Prediction: 89.25 - Expected=87.72
300      Prediction: 88.45 - Expected=87.89
RMSE: 1.5057
[37]:
fig, ax = plt.subplots()
ax.plot(test.to_list(),color="b")
ax.plot(predictions,color="r")
plt.show()
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_23_0.png

Reflexiones

  • ¿Qué ocurre si el modelo no se entrena con cada nueva observación?

  • ¿Qué ocurre si cambiamos de modelo ARIMA(1,0,0) o ARIMA(5,1,2)?

  • ¿Qué ocurre si reducimos el threshold de entreno?

  • ¿Qué ocurre si el periodo de la serie no es mensual?

Predicción con un algoritmo supervisado

El objetivo de nuestro analisis será predecir el valor del cierre de MSFT en función del valor del cierre de IBM y GOOGLE y el uso de dos indices Dow Jones y S&P500

[13]:
# Cargamos datos
stk_tickers = ['MSFT', 'IBM', 'GOOGL'] # Microsoft, IBM, Google
data = pdr.get_data_yahoo(stk_tickers, start="2017-01-01", end="2020-12-31")

data.to_csv("Stickers.csv")
idx_tickers = ['^GSPC', '^DJI'] # S&P 500 y Dow Jones index
idx_data = pdr.get_data_yahoo(idx_tickers, start="2017-01-01", end="2020-12-31")
idx_data.to_csv("idx_data.csv")
[                       0%%                      ][*********************100%%**********************]  3 of 3 completed
[*********************100%%**********************]  2 of 2 completed
[69]:
data = pd.read_csv("Stickers.csv")
data
[69]:
Unnamed: 0 Adj Close Adj Close.1 Adj Close.2 Close Close.1 Close.2 High High.1 High.2 Low Low.1 Low.2 Open Open.1 Open.2 Volume Volume.1 Volume.2
0 NaN GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT
1 Date NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 2017-01-03 40.4005012512207 115.34630584716797 57.138736724853516 40.4005012512207 159.8374786376953 62.58000183105469 40.571998596191406 160.48757934570312 62.84000015258789 39.84450149536133 158.70936584472656 62.130001068115234 40.03099822998047 159.6558380126953 62.790000915527344 39180000 3069278 20694100
3 2017-01-04 40.38850021362305 116.77437591552734 56.88307189941406 40.38850021362305 161.81643676757812 62.29999923706055 40.67150115966797 162.3996124267578 62.75 40.205501556396484 160.0 62.119998931884766 40.49449920654297 160.3919677734375 62.47999954223633 30306000 3536944 21340000
4 2017-01-05 40.6510009765625 116.38805389404297 56.88307189941406 40.6510009765625 161.28106689453125 62.29999923706055 40.6870002746582 161.9407196044922 62.65999984741211 40.29600143432617 159.90440368652344 62.029998779296875 40.375 161.806884765625 62.189998626708984 26810000 2805686 24876000
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1003 2020-12-23 86.4114990234375 102.3642578125 215.23953247070312 86.4114990234375 118.45124053955078 221.02000427246094 87.20549774169922 119.7036361694336 223.55999755859375 86.05999755859375 118.29827880859375 220.8000030517578 86.19650268554688 118.43212127685547 223.11000061035156 22974000 2817819 18699600
1004 2020-12-24 86.70800018310547 103.01693725585938 216.9242706298828 86.70800018310547 119.20649719238281 222.75 87.12049865722656 119.59847259521484 223.61000061035156 86.21749877929688 118.74761199951172 221.1999969482422 86.44999694824219 119.50286865234375 221.4199981689453 9312000 1842111 10550600
1005 2020-12-28 88.697998046875 103.12435150146484 219.07647705078125 88.697998046875 119.3307876586914 224.9600067138672 89.3499984741211 121.03250122070312 226.02999877929688 87.09100341796875 118.98661804199219 223.02000427246094 87.24549865722656 119.59847259521484 224.4499969482422 27650000 3781499 17933500
1006 2020-12-29 87.88800048828125 102.28163146972656 218.28765869140625 87.88800048828125 118.35564422607422 224.14999389648438 89.42350006103516 119.96176147460938 227.17999267578125 87.75550079345703 117.82026672363281 223.5800018310547 89.36150360107422 119.83747863769531 226.30999755859375 19726000 3647402 17403200
1007 2020-12-30 86.8125 102.7277603149414 215.8822784423828 86.8125 118.87189483642578 221.67999267578125 88.38800048828125 119.35946655273438 225.6300048828125 86.4000015258789 118.193115234375 221.47000122070312 88.25 118.35564422607422 225.22999572753906 21026000 3535794 20272300

1008 rows × 19 columns

[71]:
data = pd.read_csv("Stickers.csv",header=[0, 1],skiprows=[2])
data.index = pd.DatetimeIndex(data[data.columns[0]])
data

[71]:
Unnamed: 0_level_0 Adj Close Close High Low Open Volume
Unnamed: 0_level_1 GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03 2017-01-03 40.400501 115.346306 57.138737 40.400501 159.837479 62.580002 40.571999 160.487579 62.840000 39.844501 158.709366 62.130001 40.030998 159.655838 62.790001 39180000 3069278 20694100
2017-01-04 2017-01-04 40.388500 116.774376 56.883072 40.388500 161.816437 62.299999 40.671501 162.399612 62.750000 40.205502 160.000000 62.119999 40.494499 160.391968 62.480000 30306000 3536944 21340000
2017-01-05 2017-01-05 40.651001 116.388054 56.883072 40.651001 161.281067 62.299999 40.687000 161.940720 62.660000 40.296001 159.904404 62.029999 40.375000 161.806885 62.189999 26810000 2805686 24876000
2017-01-06 2017-01-06 41.260502 116.960693 57.376114 41.260502 162.074570 62.840000 41.448002 162.447418 63.150002 40.575001 160.152969 62.040001 40.749500 161.271515 62.299999 40342000 3080993 19922900
2017-01-09 2017-01-09 41.359001 115.663658 57.193520 41.359001 160.277252 62.639999 41.521500 162.332703 63.080002 41.081001 160.248566 62.540001 41.318501 162.017212 62.759998 28178000 3336635 20382700
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2020-12-23 2020-12-23 86.411499 102.364258 215.239532 86.411499 118.451241 221.020004 87.205498 119.703636 223.559998 86.059998 118.298279 220.800003 86.196503 118.432121 223.110001 22974000 2817819 18699600
2020-12-24 2020-12-24 86.708000 103.016937 216.924271 86.708000 119.206497 222.750000 87.120499 119.598473 223.610001 86.217499 118.747612 221.199997 86.449997 119.502869 221.419998 9312000 1842111 10550600
2020-12-28 2020-12-28 88.697998 103.124352 219.076477 88.697998 119.330788 224.960007 89.349998 121.032501 226.029999 87.091003 118.986618 223.020004 87.245499 119.598473 224.449997 27650000 3781499 17933500
2020-12-29 2020-12-29 87.888000 102.281631 218.287659 87.888000 118.355644 224.149994 89.423500 119.961761 227.179993 87.755501 117.820267 223.580002 89.361504 119.837479 226.309998 19726000 3647402 17403200
2020-12-30 2020-12-30 86.812500 102.727760 215.882278 86.812500 118.871895 221.679993 88.388000 119.359467 225.630005 86.400002 118.193115 221.470001 88.250000 118.355644 225.229996 21026000 3535794 20272300

1006 rows × 19 columns

[55]:
idx_data = pd.read_csv("idx_data.csv",header=[0, 1],skiprows=[2])
idx_data.index = pd.DatetimeIndex(idx_data[idx_data.columns[0]])
[72]:
idx_data[:5]
[72]:
Unnamed: 0_level_0 Adj Close Close High Low Open Volume
Unnamed: 0_level_1 ^DJI ^GSPC ^DJI ^GSPC ^DJI ^GSPC ^DJI ^GSPC ^DJI ^GSPC ^DJI ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03 2017-01-03 19881.759766 2257.830078 19881.759766 2257.830078 19938.529297 2263.879883 19775.929688 2245.129883 19872.859375 2251.570068 339180000 3773010000
2017-01-04 2017-01-04 19942.160156 2270.750000 19942.160156 2270.750000 19956.140625 2272.820068 19878.830078 2261.600098 19890.939453 2261.600098 280010000 3768890000
2017-01-05 2017-01-05 19899.289062 2269.000000 19899.289062 2269.000000 19948.599609 2271.500000 19811.119141 2260.449951 19924.560547 2268.179932 269920000 3785080000
2017-01-06 2017-01-06 19963.800781 2276.979980 19963.800781 2276.979980 19999.630859 2282.100098 19834.080078 2264.060059 19906.960938 2271.139893 277700000 3342080000
2017-01-09 2017-01-09 19887.380859 2268.899902 19887.380859 2268.899902 19943.779297 2275.489990 19887.380859 2268.899902 19931.410156 2273.590088 287510000 3219730000
[57]:
# Tomamos un periodo de 5 dias, de valores medios sobre los indices y cierres de  L - V.
# Primero, unimos lo datos en un DF
# print(type(data))
# print(type(idx_data))
# print(data.index)
# print(idx_data.index)
# Ambos son DF. Entonces es una operación de Merge
dataset = data.merge(idx_data, how="inner",left_index=True,right_index=True)
print(dataset[:5])
                                         Unnamed: 0_level_0_x  Adj Close  \
                                           Unnamed: 0_level_1      GOOGL
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                         2017-01-03  40.400501
2017-01-04                                         2017-01-04  40.388500
2017-01-05                                         2017-01-05  40.651001
2017-01-06                                         2017-01-06  41.260502
2017-01-09                                         2017-01-09  41.359001

                                                                     Close  \
                                                 IBM       MSFT      GOOGL
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                115.346306  57.138737  40.400501
2017-01-04                                116.774376  56.883072  40.388500
2017-01-05                                116.388054  56.883072  40.651001
2017-01-06                                116.960693  57.376114  41.260502
2017-01-09                                115.663658  57.193520  41.359001

                                                                      High  \
                                                 IBM       MSFT      GOOGL
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                159.837479  62.580002  40.571999
2017-01-04                                161.816437  62.299999  40.671501
2017-01-05                                161.281067  62.299999  40.687000
2017-01-06                                162.074570  62.840000  41.448002
2017-01-09                                160.277252  62.639999  41.521500

                                                                 ...  \
                                                 IBM       MSFT  ...
(Unnamed: 0_level_0, Unnamed: 0_level_1)                         ...
2017-01-03                                160.487579  62.840000  ...
2017-01-04                                162.399612  62.750000  ...
2017-01-05                                161.940720  62.660000  ...
2017-01-06                                162.447418  63.150002  ...
2017-01-09                                162.332703  63.080002  ...

                                                 Close               \
                                                  ^DJI        ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                19881.759766  2257.830078
2017-01-04                                19942.160156  2270.750000
2017-01-05                                19899.289062  2269.000000
2017-01-06                                19963.800781  2276.979980
2017-01-09                                19887.380859  2268.899902

                                                  High               \
                                                  ^DJI        ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                19938.529297  2263.879883
2017-01-04                                19956.140625  2272.820068
2017-01-05                                19948.599609  2271.500000
2017-01-06                                19999.630859  2282.100098
2017-01-09                                19943.779297  2275.489990

                                                   Low               \
                                                  ^DJI        ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                19775.929688  2245.129883
2017-01-04                                19878.830078  2261.600098
2017-01-05                                19811.119141  2260.449951
2017-01-06                                19834.080078  2264.060059
2017-01-09                                19887.380859  2268.899902

                                                  Open               \
                                                  ^DJI        ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                19872.859375  2251.570068
2017-01-04                                19890.939453  2261.600098
2017-01-05                                19924.560547  2268.179932
2017-01-06                                19906.960938  2271.139893
2017-01-09                                19931.410156  2273.590088

                                             Volume
                                               ^DJI       ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-03                                339180000  3773010000
2017-01-04                                280010000  3768890000
2017-01-05                                269920000  3785080000
2017-01-06                                277700000  3342080000
2017-01-09                                287510000  3219730000

[5 rows x 32 columns]

Reconvertimos las series a una frecuencia de 5 días laborales

[73]:
# Reconvertimos la frecuencia agrupando los valores.
# https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.resample.html
data5days = dataset.resample("1W-MON").mean()  # averaging all values along 5 days starting on Mondays
data5days.head()

[73]:
Adj Close Close High Low ... Close High Low Open Volume
GOOGL IBM MSFT GOOGL IBM MSFT GOOGL IBM MSFT GOOGL ... ^DJI ^GSPC ^DJI ^GSPC ^DJI ^GSPC ^DJI ^GSPC ^DJI ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-09 40.811901 116.226617 57.094903 40.811901 161.057361 62.532000 40.980000 161.921606 62.896001 40.400401 ... 19914.878125 2268.691992 19957.335938 2273.158008 19837.467969 2260.027979 19905.346094 2265.216016 290864000.0 3.577758e+09
2017-01-16 41.454250 115.311779 57.321336 41.454250 159.789677 62.780000 41.554250 160.693115 63.142500 41.189251 ... 19896.634766 2272.324951 19952.964355 2276.262512 19822.260254 2262.965027 19900.620117 2270.549988 301407500.0 3.455635e+09
2017-01-23 41.534500 116.330086 57.162468 41.534500 161.200766 62.606000 41.677701 161.780115 62.864000 41.245400 ... 19798.198047 2267.995947 19842.650000 2273.432031 19736.757812 2261.320020 19814.990234 2269.583936 337072000.0 3.353020e+09
2017-01-30 42.338200 122.245384 58.869874 42.338200 169.397702 64.475999 42.755399 170.248566 64.816001 42.054700 ... 20029.408203 2290.141992 20060.282031 2294.039990 19958.350391 2281.860010 19999.477734 2288.083936 352768000.0 3.602052e+09
2017-02-06 40.954400 120.744125 58.201519 40.954400 167.317398 63.744000 41.201500 167.986615 64.106001 40.713400 ... 19952.764062 2285.850049 19997.016016 2289.337988 19885.774219 2277.549951 19937.026172 2283.824023 354716000.0 3.707408e+09

5 rows × 30 columns

[74]:
# Limpiamos series innecesarias
for ix,col in enumerate(data5days.columns):
    print(ix,col)
0 ('Adj Close', 'GOOGL')
1 ('Adj Close', 'IBM')
2 ('Adj Close', 'MSFT')
3 ('Close', 'GOOGL')
4 ('Close', 'IBM')
5 ('Close', 'MSFT')
6 ('High', 'GOOGL')
7 ('High', 'IBM')
8 ('High', 'MSFT')
9 ('Low', 'GOOGL')
10 ('Low', 'IBM')
11 ('Low', 'MSFT')
12 ('Open', 'GOOGL')
13 ('Open', 'IBM')
14 ('Open', 'MSFT')
15 ('Volume', 'GOOGL')
16 ('Volume', 'IBM')
17 ('Volume', 'MSFT')
18 ('Adj Close', '^DJI')
19 ('Adj Close', '^GSPC')
20 ('Close', '^DJI')
21 ('Close', '^GSPC')
22 ('High', '^DJI')
23 ('High', '^GSPC')
24 ('Low', '^DJI')
25 ('Low', '^GSPC')
26 ('Open', '^DJI')
27 ('Open', '^GSPC')
28 ('Volume', '^DJI')
29 ('Volume', '^GSPC')
[75]:
columns_si = [3,4,5,15,16,17,20,21]
columns_to_remove = set(range(len(data5days.columns)))-set(columns_si)
columns_to_remove = list(columns_to_remove)
columns_to_remove
[75]:
[0,
 1,
 2,
 6,
 7,
 8,
 9,
 10,
 11,
 12,
 13,
 14,
 18,
 19,
 22,
 23,
 24,
 25,
 26,
 27,
 28,
 29]
[76]:
data5days = data5days.drop(data5days.columns[columns_to_remove],axis=1)
data5days.head()
[76]:
Close Volume Close
GOOGL IBM MSFT GOOGL IBM MSFT ^DJI ^GSPC
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-09 40.811901 161.057361 62.532000 32963200.0 3165907.2 21443140.0 19914.878125 2268.691992
2017-01-16 41.454250 159.789677 62.780000 25812500.0 3535767.5 20125200.0 19896.634766 2272.324951
2017-01-23 41.534500 161.200766 62.606000 29210800.0 6789606.8 22419380.0 19798.198047 2267.995947
2017-01-30 42.338200 169.397702 64.475999 56426000.0 4919526.2 33673920.0 20029.408203 2290.141992
2017-02-06 40.954400 167.317398 63.744000 35357600.0 3245800.8 32173440.0 19952.764062 2285.850049
[77]:
# Facilitamosel manejo de columnas, eliminando un nivel superior (son : close o volume)
data5days.columns = data5days.columns.droplevel(0)
data5days.columns
[77]:
Index(['GOOGL', 'IBM', 'MSFT', 'GOOGL', 'IBM', 'MSFT', '^DJI', '^GSPC'], dtype='object')
[78]:
# Y renombramos columnas: MUY IMPORTANTE EL ORDEN!!
columns_names = ["CloseGOOGL","CloseIBM","CloseMSFT","VolGOOGL","VolIBM","VolMSFT","DJI","SP500"]
data5days.columns = columns_names

data5days.head()

[78]:
CloseGOOGL CloseIBM CloseMSFT VolGOOGL VolIBM VolMSFT DJI SP500
(Unnamed: 0_level_0, Unnamed: 0_level_1)
2017-01-09 40.811901 161.057361 62.532000 32963200.0 3165907.2 21443140.0 19914.878125 2268.691992
2017-01-16 41.454250 159.789677 62.780000 25812500.0 3535767.5 20125200.0 19896.634766 2272.324951
2017-01-23 41.534500 161.200766 62.606000 29210800.0 6789606.8 22419380.0 19798.198047 2267.995947
2017-01-30 42.338200 169.397702 64.475999 56426000.0 4919526.2 33673920.0 20029.408203 2290.141992
2017-02-06 40.954400 167.317398 63.744000 35357600.0 3245800.8 32173440.0 19952.764062 2285.850049

Con los datos preparados, podemos analizar los valores de estas series y sus posibles relaciones

[79]:
# Su desomposición
res = sm.tsa.seasonal_decompose(data5days.CloseMSFT,period=8) # que implica este 8?
fig = res.plot()
plt.show()
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_40_0.png

o incluso podemos mirar la matriz de correlación para detectar algún patrón común entre los datos

[103]:
#con matplotlib
import matplotlib.pyplot as plt
correlation = data5days.corr()
plt.matshow(correlation)
[103]:
<matplotlib.image.AxesImage at 0x29eca61f0>
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_42_1.png
[104]:
#con seaborn
import seaborn as sns
sns.heatmap(correlation, vmax=1, square=True,annot=True,cmap='cubehelix')
[104]:
<AxesSubplot:>
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_43_1.png

Con los datos preparados podemos aplicar modelos supervisados de predicción

Nuestra variable objetivo será el valor de cierre de MSFT

[81]:
# Aplicación de modelos de supervisado

from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split

# Vamos aplicar diferentes modelos al mismo tiempo!
models = [LinearRegression,SVR,RandomForestRegressor]

# Preparamos las variables: features y target
X = data5days.drop(labels=["CloseMSFT"],axis=1)
Y = data5days.CloseMSFT

# Particionamos los datos de entreno y test al 80%
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.20, random_state=33)

[82]:
# Variables de salida
model_names = []
train_results = []
test_results = []

# Entrenamos y testeamos
for model in models:
    print("Model: %s"%model.__name__)
    model_names.append(model.__name__)
    res = model().fit(X=X_train,y=y_train)
    train_result = mean_squared_error(res.predict(X_train), y_train)
    train_results.append(train_result)

    test_result = mean_squared_error(res.predict(X_test), y_test)
    test_results.append(test_result)

Model: LinearRegression
Model: SVR
Model: RandomForestRegressor
[109]:
import numpy as np
ind = np.arange(len(models))
width = 0.35 # the width of the bars

fig,ax = plt.subplots()
plt.bar(ind - width/2, train_results, width=width, label='Train Error')
plt.bar(ind + width/2, test_results, width=width, label='Test Error')
plt.legend()
ax.set_xticks(ind)
ax.set_xticklabels(model_names)
plt.show()
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_47_0.png

Más allá de la regresión lineal simple

En la sesión de hoy seguiremos trabajando con las técnicas de regresión, profundizaremos en las técnicas de regresión lineal y veremos otras técnicas de regresión no lineal. Por otra parte revisaremos las métricas más usadas para este problema. Finalmente trabajaremos en la automatización de los procesos de aprendizaje automático.

  1. Uso de la regularización

  2. Regresión no lineal: un ejemplo

  3. Métricas de regresión

  4. Buenas prácticas I: Conjuntos de entrenamiento y test.

1. Uso de la regularización

Una forma de encontrar una buena relación entre el sesgo y la varianza es ajustar la complejidad del modelo a través de la regularización. La regularización es un método muy útil para manejar colinealidad (alta correlación entre características), filtrar el ruido de los datos y eventualmente evitará el overfitting. El concepto detrás de la regularización es introducir información adicional (sesgo) para penalizar los valores extremas de los parámetros.

Si tenemos que la expresión para una regresión lineal es:

\[y = w_0x_0+w_1x_1+...+w_nx_n=\sum^m_{i=0} w_i X_{i}\]

Y la función que queremos minimizar es:

\[J(w) = \sum^n_{i=1}(y^{(i)}-\hat{y}^{(i)})^2\]

Ridge Regression

La Ridge Regression es un modelo que usa una penalización aplicando la norma L2 donde simplemente agregamos la suma al cuadrado de los pesos de nuestra función de coste:

\[J(w)_{Ridge} = \sum^n_{i=1}(y^{(i)}-\hat{y}^{(i)})^2 + \lambda\| w \|_2\]

El uso de esta regresión sería el siguiente:

[110]:
from sklearn.linear_model import Ridge

ridge = Ridge(alpha=1.0)

Al aumentar el valor del hiperparámetro \(\lambda\), aumentamos la fuerza de la regularización y encogemos los pesos de nuestro modelo. Hay que tener en cuenta que no regularizamos el término del intercepto (\(w_0\)).

Lasso Regression (L1)

Un enfoque alternativo que puede conducir a modelos dispersos es la regresión LASSO. Según sea el valor del término de regularización, ciertos pesos pueden coger el valor cero, lo que hace que este tipo de regresión LASSO también sea útil como técnica de selección de características supervisada.

La función a minimizar en este caso es:

\[J(w)_{Lasso} = \sum^n_{i=1}(y^{(i)}-\hat{y}^{(i)})^2 + \lambda\| w \|_1\]

El uso de esta regresión sería el siguiente:

[111]:
from sklearn.linear_model import Lasso
lasso = Lasso(alpha=1.0)

Elastic Net Regression

Sin embargo, una limitación de LASSO es que selecciona como máximo \(n\) variables si $ m > n $. Elastic Net ofrece un buen compromiso entre la Ridge regression y LASSO, que tiene una Penalización L1 para generar una selección de características y penalización L2 para superar algunas de las limitaciones de la regresión LASSO, como es el número de variables seleccionadas.

La función a minimizar en este caso es:

\[J(w)_{Elasticnet} = \sum^n_{i=1}(y^{(i)}-\hat{y}^{(i)})^2 + \lambda_1 \sum^n_{i=1} w^2_j + \lambda_2 \sum^n_{i=1} w_j\]

El uso de esta regresión sería el siguiente:

[112]:
from sklearn.linear_model import ElasticNet
elastic = ElasticNet(alpha=1.0, l1_ratio=0.5)

Regresión polinómica

La regresión polinomial es un caso especial de regresión lineal en el que se construyen nuevas características en función del grado del polinomio que se quiera construir. En las secciones anteriores, hemos asumido una relación lineal entre las variables del modelo y la variable objetivo. Una forma de no tener en cuenta el supuesto de linealidad es usar un modelo de regresión polinomial, es decir agregando términos polinomiales.

\[y = w_0 + w_1x + w_2x^2 + ... + w_dx^d\]

donde \(d\) es el grado del polinomio. Aunque podemos usar un polinomio para modelar una relación no lineal, todavía se considera una modelo de regresión lineal debido a los coeficientes \(w\).

[113]:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

def true_fun(X):
    return np.cos(1.5 * np.pi * X)

np.random.seed(0)
n_samples = 30
degrees = [1, 4, 15]

X = np.sort(np.random.rand(n_samples))
y = true_fun(X) + np.random.randn(n_samples) * 0.1

plt.figure(figsize=(14, 5))

for i in range(len(degrees)):
    ax = plt.subplot(1, len(degrees), i + 1)
    plt.setp(ax, xticks=(), yticks=())


    linear_regression = LinearRegression()
    polynomial_features = PolynomialFeatures(degree=degrees[i])
    polynomial = polynomial_features.fit_transform(X[:, np.newaxis])

    linear_regression.fit(polynomial, y)

    X_test = np.linspace(0, 1, 100)
    plt.plot(X, linear_regression.predict(polynomial), label="Model")
    plt.plot(X_test, true_fun(X_test), label="True function")
    plt.scatter(X, y, edgecolor="b", s=20, label="Samples")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.xlim((0, 1))
    plt.ylim((-2, 2))
    plt.legend(loc="best")

    plt.title( "Degree: " + str(degrees[i]))

plt.show()
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_55_0.png

Regresión no lineal: un ejemplo

Cuando queremos resolver un problema usando un modelo más complejo podemos elegir alguno de los siguientes:

  • Árboles de regresión

  • Bosques de regresión

  • SVM para la regresión

  • Y todos los modelos que existen…

Realizaremos un ejercicio con los árboles de decisión debido a que son fáciles de entender y que en sesiones futuras trabajaremos con los análogos de los Bosques de regresión y las SVM para problemas de clasificación.

Los árboles de decisión (DT) son un método de aprendizaje supervisado no paramétrico que se utiliza para clasificación y regresión. El objetivo es crear un modelo que prediga el valor de una variable objetivo mediante el aprendizaje de reglas de decisión simples inferidas de las características de los datos. Un árbol puede verse como una aproximación de una función a trozos.

En sci-kit la clase que modela este tipo de árboles se llama DecisionTreeRegressor enlace

Veamos un ejemplo sencillo:

[114]:
# Importamos las librerias necesarias
import numpy as np
from sklearn.tree import DecisionTreeRegressor
import matplotlib.pyplot as plt

# Generamos un conjunto de datos (función seno)
rng = np.random.RandomState(1)
X = np.sort(5 * rng.rand(80, 1), axis=0)
y = np.sin(X).ravel()
y[::5] += 3 * (0.5 - rng.rand(16))

# Entrenamos el modelo

regr_1 = DecisionTreeRegressor(max_depth=2)
regr_1.fit(X, y)

# Realizamos una predicción
X_test = np.expand_dims(np.arange(0.0, 5.0, 0.01), 1)
y_1 = regr_1.predict(X_test)

# Mostramos los resultados
plt.figure()
plt.scatter(X, y, s=20, edgecolor="black", c="darkorange", label="data")
plt.plot(X_test, y_1, color="cornflowerblue", label="max_depth=2", linewidth=2)
plt.xlabel("data")
plt.ylabel("target")
plt.title("Decision Tree Regression")
plt.legend()
plt.show()
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_57_0.png

Vamos a ver la estructura del primer árbol que hemos entrenado para entender mejor el proceso que ha seguido:

[115]:
from sklearn import tree
plt.figure(figsize=(9,9))
tree.plot_tree(regr_1)
[115]:
[Text(0.5, 0.8333333333333334, 'X[0] <= 3.133\nsquared_error = 0.547\nsamples = 80\nvalue = 0.122'),
 Text(0.25, 0.5, 'X[0] <= 0.514\nsquared_error = 0.231\nsamples = 51\nvalue = 0.571'),
 Text(0.125, 0.16666666666666666, 'squared_error = 0.192\nsamples = 11\nvalue = 0.052'),
 Text(0.375, 0.16666666666666666, 'squared_error = 0.148\nsamples = 40\nvalue = 0.714'),
 Text(0.75, 0.5, 'X[0] <= 3.85\nsquared_error = 0.124\nsamples = 29\nvalue = -0.667'),
 Text(0.625, 0.16666666666666666, 'squared_error = 0.124\nsamples = 14\nvalue = -0.452'),
 Text(0.875, 0.16666666666666666, 'squared_error = 0.041\nsamples = 15\nvalue = -0.869')]
../../../_images/notebooks_Part2_03_MachineLearningNotes_04_Regression_59_1.png

Ejercicio

  • Cambia el parámetro max_depth del modelo. Que pasa si es 1? Que pasa si es 5 o mayor?

  • Visualiza el árbol resultante.

  • Cuál es la medida de Score de este método (mirar documentación)? Como cambia al permitir árboles más profundos?

  • Busca en la documentación de Sci-kit de los Bosques de Regresión ( RandomForestRegressor ) compara el score obtenido cuando incrementamos el número de árboles, manteniendo la misma profundidad.

Métricas para problemas de regression

Las métricas más comunes en los problemas de regresión son el error cuadrático y el error absoluto, y sus distintas modificaciones.

3.1 Errores cuadráticos

El error cuadrático (Squared Error) de un valor predicho con respecto al valor real, se calcula cómo:

\[SE = \sum_j\left[\hat{y}_j - y_j\right]^2\]

Error cuadrático medio (Mean Squared Error) Da una idea del error de nuestras predicciones dando más peso a los errores grandes.

\[MSE = \frac{1}{m}\sum_j^m \left[\hat{y}_j - y_j\right]^2\]

Raiz del error cuadrático medio (Root Mean Square Error) La raíz cuadrada del MSE produce el error de la raíz cuadrada de la media o la desviación de la raíz cuadrada media (RMSE o RMSD). Tiene las mismas unidades que la cantidad que estima. Para un estimador sin sesgo (bies), el RMSE es la raíz cuadrada de la varianza, es decir la desviación estandar.

\[RMSE = \sqrt{\frac{1}{m}\sum_j^m \left[\hat{y}_j - y_j\right]^2}\]

A pesar de ser una de las métricas más utilizadas, tiene el inconveniente de ser sensible a los valores extremos (outliers). Cuando este comportamiento pueda suponer un problema, los errores absolutos pueden darnos una mejor medida de rendimiento.

[33]:
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

sns.set_style('darkgrid')
sns.set_context('notebook')
sns.set(rc={'figure.figsize':(12, 8)})
[34]:
## Escribe una función que devuelva el MSE y el RMSE dados dos arrays de numpy

def MSE(x1, x2):
    """
    Returns the Root Mean Squared Error of the two input vectors
    """
    sq_error = ##
    mean_sq_error = np.##
    return mean_sq_error

def RMSE(x1, x2):
    """
    Returns the Root Mean Squared Error of two vectors. Depends on the MSE function
    """
    mse = MSE(x1, x2)
    root_mse = np.##
    return root_mse
  Input In [34]
    sq_error = ##
               ^
SyntaxError: invalid syntax

3.2 Errores absolutos

El error absoluto (Absolute Error) se define cómo:

\[AE = \sum_j \left|\hat{y}_j - y_j\right|\]

Error absoluto medio (Mean Absolute Error)

Es más robusto a los valores extremos y su interpretabilidad es más alta que la del RMSE ya que también está en las unidades de la variable a predecir con la ventaja de que el dato no ha sufrido ninguna transformación.

\[MAE = \frac{1}{m} \sum_j^m \left|\hat{y}_j - y_j\right|\]

Error absoluto medio porcentual (Mean Average Percentage Error)

A pesar de su simpleza, presenta varios inconvenientes a la hora de usarlo de forma práctica. Por ejemplo, no puede usarse cuando el valor de referencia es 0. Además, si se usa para elegir métodos predictivos seleccionará de forma sistemática un metodo que prediga valores bajos. wiki

\[MAPE = \frac{1}{m} \sum_j^m \left|\frac{\hat{y}_j - y_j}{y_j}\right|\]

3.3 Generalización

Las dos métricas expuestas anteriormente pueden considerarse como distancias entre el vector de valores reales y el predicho. De hecho, el RMSE corresponde a la distancia euclidiana, también conocida como norma \(l_2\) o \(\lVert{v}\rVert_2\).

Por otro lado, el MAE corresponde a la norma \(l_1\) o \(\lVert{v}\rVert_1\). A esta distancia se la conoce como distancia de manhattan, porque sólo se puede viajar de un bloque a otro de la ciudad a traves de calles ortogonales.

De forma general, una norma \(l_k\) o \(\lVert{v}\rVert_k\) se calcula:

\[\lVert{v}\rVert_k = \left(|v_0|^k + ...+ |v_m|^k \right)^\frac{1}{k}\]

El concepto de distancia será particularmente útil en los problemas de segmentación (aprendizaje no supervisado)

3.4 Coeficiente de determinación (\(R^2\))

El coeficiente determina la calidad del modelo para replicar los resultados, y la proporción de variación de los resultados que puede explicarse por el modelo. El valor más alto obtenible será 1, aunque hay casos en los que puede presentar valores negativos. De forma intuitiva, \(R^2\) compara el “fit” de nuestro modelo al de una linea recta horizontal. Dada una regresión lineal simple, un \(R^2\) negativo sólo es posible cuando la ordenada en el origen o la pendiente están restringidas de forma que el mejor modelo es peor que una linea horizontal.

Si representamos la varianza de la variable dependiente por \(\sigma^{2}\) y la varianza residual por \(\sigma _{r}^{2}\), el coeficiente de determinación viene dado por la siguiente ecuación:

\[R^{2}=1- \frac{\sigma_{r}^{2}}{\sigma ^{2}}\]

Siendo \(\hat{y}_i\) el valor predicho de la muestra i y \(y_i\) el valor real, el \(R^2\) estimado sobre \(n_{\text{muestras}}\) se define como:

\[R^2(y, \hat{y}) = 1 - \frac{\sum_{i=0}^{n_{\text{samples}} - 1} (y_i - \hat{y}_i)^2}{\sum_{i=0}^{n_\text{samples} - 1} (y_i - \bar{y})^2}\]

donde

\[\bar{y} = \frac{1}{n_{\text{samples}}} \sum_{i=0}^{n_{\text{samples}} - 1} y_i\]

3.4.1 Consideraciones R²

  1. R² no puede determinar si los coeficientes y las predicciones tienen bies: Hay que checkear los residos –> Si observamos patrones en los plots de residuos es indicativo de un mal ajuste a pesar de un R2 elevado

  2. Cada vez que añadimos un predictor a un modelo, el R² aumenta aunque sea por suerte, pero nunca decrece. Por consiguiente, un modelo con muchos terminos puede parecer mejor simplemente por el hecho de tener más terminos. Para prevenir este efecto, podemos usar el adjusted R², una versión modificada que se ajusta al número de predictores en el modelo. De ésta forma, el R² solo aumenta si el nuevo término mejora el modelo más que por mera suerte. Siempre es más bajo que el R²

\[\bar{R}^2 = 1 - \frac{N-1}{N-k-1}(1-R^2)\]

n – numero de observaciones

k – numero de parametros

License: CC BY 4.0 Isaac Lera and Gabriel Moya Universitat de les Illes Balears isaac.lera@uib.edu, gabriel.moya@uib.edu