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:

Linear regressión / Regresión Lineal: https://es.wikipedia.org/wiki/Regresi%C3%B3n_lineal
Regresión logística: https://es.wikipedia.org/wiki/Regresi%C3%B3n_log%C3%ADstica
Máquinas de vectores de soporte / SVM: https://es.wikipedia.org/wiki/M%C3%A1quinas_de_vectores_de_soporte
CART / Classification And Regression Trees https://www.nature.com/articles/nmeth.4370
Gradient Boosting / Potenciación del gradiente: https://es.wikipedia.org/wiki/Gradient_boosting
Random Forest / Bosques aleatorios: https://es.wikipedia.org/wiki/Random_forest
Artificial Neural Networks (*): https://en.wikipedia.org/wiki/Artificial_neural_network
K-vecinos cercanos / K-nearest neighbors (k-nn): https://es.wikipedia.org/wiki/K_vecinos_m%C3%A1s_pr%C3%B3ximos
Análisis Discriminantel Lineal / LDA: https://es.wikipedia.org/wiki/An%C3%A1lisis_discriminante_lineal
etc.

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

|2326779ae35e42bebb5bd60061c33c65|

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.

Uso de la regularización
Regresión no lineal: un ejemplo
Métricas de regresión
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²

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
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

Isaac Lera and Gabriel Moya Universitat de les Illes Balears isaac.lera@uib.edu, gabriel.moya@uib.edu

Algoritmos Supervisados: Predicción

Obtención de datos

Análisis de la serie temporal

Predicción con ARIMA

Reflexiones

Predicción con un algoritmo supervisado

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

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

1. Uso de la regularización

Ridge Regression

Lasso Regression (L1)

Elastic Net Regression

Regresión polinómica

Regresión no lineal: un ejemplo

Métricas para problemas de regression

3.1 Errores cuadráticos

3.2 Errores absolutos

Error absoluto medio (Mean Absolute Error)

Error absoluto medio porcentual (Mean Average Percentage Error)

3.3 Generalización

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

3.4.1 Consideraciones R²