Analisis de PCA en Python

La dimensionalidad en los datasets, suele ser un problema muy frecuente en los modelos que utilizan técnicas de machine learning o Deep learning. En este nuevo artículo, vamos a tratar de una forma lo mas simplificada posible, un método para reducir la dimensionalidad, en este caso utilizaremos el análisis de componentes principales (PCA) en Python, mediante la librería de sklearn.

Para poner un poco en contexto, vamos a dejar unos enlaces, con una explicacion matematica sobre lo que estamos haciendo, se puede consultar esta informacion aqui:

https://www.cienciadedatos.net/documentos/35_principal_component_analysis

https://towardsdatascience.com/pca-using-python-scikit-learn-e653f8989e60

Mediante el análisis de los componentes principales (PCA) vamos a transformar un dataset con demasiadas features, algunas relevantes, pero otras altamente correlacionadas, en un nuevo dataset, donde aplicamos PCA y reducimos su dimensionalidad, hasta un punto, que podamos afirmar, que no perdemos demasiada información.

Los modelos de Deep learning aplicados para trading,que podemos programar en Python, de una forma realmente sencilla funcionan muchísimo mejor, con datasets con pocas features, y muy relevantes, además de modelos de redes neuronales, con pocas capas y pocas neuronas por capa, además de nuestros preciados dropouts, que nos eliminan información de forma aleatoria, añadiendo dificultad adicional al modelo, y evitando la sobre optimización o el overfit de dicho modelo.

Este articulo será utilizado de base, para los futuros artículos, donde implementaremos técnicas de machine learning y Deep learning en Python para poder predecir la dirección del precio.

ATENCION! El proceso de formación de los precios es un proceso estocastico, un movimiento browniano, en el cual la aleatoriedad está MUY PRESENTE, y que incluso la propia predicción del precio altera el propio precio, al disponer de predicciones, es información que añadimos a nuestras decisiones o modelos, y, en consecuencia, se opera en función a ello, alterando la formación del precio.

En mi opinión es erróneo intentar predecir el precio, pues tal que retornos pasados no garantizan retornos futuros, precios pasados, no influyen de forma determinante en la formación de precios futuros.

Programando un analisis PCA en Python

Para empezar, lo primero que vamos a hacer, es importar todas las librerías necesarias, como viene siendo clásico en Python…

In [ ]:
import pandas as pd
import numpy as np
import datetime
import yfinance as yf
import quandl
import matplotlib.pyplot as plt
import seaborn as sns
import ta
import warnings
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn import metrics 
from sklearn.decomposition import PCA
from sklearn.utils import class_weight
from sklearn.utils.class_weight import compute_class_weight
from sklearn.metrics import classification_report
import matplotlib.dates as mdates
import quantstats as qt
import warnings
%matplotlib inline
warnings.filterwarnings(action='ignore')
plt.style.use('ggplot')
plt.rcParams["figure.figsize"] = [20,9]

Una vez tenemos todas las librerías descargadas, vamos a crear nuestro dataset, en este caso he creado una combinación de features con diferentes datos, procesándolos en distintos datasets, para su posterior manipulación. Se han utilizado los futuros continuos del Nasdaq, el S&P500 el VIX y los bonos a 10 y 30 años americanos.

Una vez tenemos todas las librerías descargadas, vamos a crear nuestro dataset, en este caso he creado una combinación de features con diferentes datos, procesándolos en distintos datasets, para su posterior manipulación. Se han utilizado los futuros continuos del Nasdaq, el S&P500 el VIX y los bonos a 10 y 30 años americanos. Además, añadimos unos datos alternativos, tales como los nuevos máximos y mínimos, los COT o el índice VXST, que es un VIX a 9 días.

In [ ]:
#Datos RAW de mercado
#indces
nq = pd.read_csv('NQ___CCS.csv',parse_dates=True,index_col=['Date'])
es = pd.read_csv('ES___CCS.csv',parse_dates=True,index_col=['Date'])
vx = pd.read_csv('VX___CCS.csv',parse_dates=True,index_col=['Date'])

#bonos
y10 = pd.read_csv('TY2__CCS.csv',parse_dates=True,index_col=['Date'])
y30 = pd.read_csv('US2__CCS.csv',parse_dates=True,index_col=['Date'])

#alternativos

newhighs = quandl.get("URC/NASDAQ_52W_HI", authtoken="xPLyKnHL4wCYWyyziuSs")
newlows = quandl.get("URC/NASDAQ_52W_LO", authtoken="xPLyKnHL4wCYWyyziuSs")
nqhighs =quandl.get("URC/NASDAQ_ADV", authtoken="xPLyKnHL4wCYWyyziuSs")
nqlows =quandl.get("URC/NASDAQ_DEC", authtoken="xPLyKnHL4wCYWyyziuSs")

cot_30 = quandl.get("CFTC/020601_FO_ALL", authtoken="xPLyKnHL4wCYWyyziuSs")
cot_30 = cot_30[['Dealer Longs','Dealer Shorts','Asset Manager Longs','Asset Manager Shorts','Leveraged Funds Longs','Leveraged Funds Shorts',
        'Total Reportable Longs','Total Reportable Shorts','Non Reportable Longs','Non Reportable Shorts']]

cot_vx = quandl.get("CFTC/1170E1_FO_ALL", authtoken="xPLyKnHL4wCYWyyziuSs")
cot_vx = cot_vx[['Dealer Longs','Dealer Shorts','Asset Manager Longs','Asset Manager Shorts','Leveraged Funds Longs','Leveraged Funds Shorts',
        'Total Reportable Longs','Total Reportable Shorts','Non Reportable Longs','Non Reportable Shorts']]

#vx9d
vx9 = pd.read_csv('https://www.cboe.com/publish/scheduledtask/mktdata/datahouse/vix9ddailyprices.csv',skiprows=3,parse_dates=True)
vx9=vx9[['Date','Close']]
vx9['Date']=vx9['Date'].str.strip('*')
vx9['Date'] = pd.to_datetime(vx9['Date'], infer_datetime_format=True)
vx9 = vx9.set_index('Date')
In [ ]:
 

Una vez que ya tenemos toda la información raw (cruda, sin preprocesar) descargada, vamos a hacer lo propio con dichos datos, vamos a obtener features desde ellos, al final estamos añadiendo la misma información, desde diferentes puntos de vista. Creamos un nuevo dataframe llamado features, conde contendremos toda esta información.

In [ ]:
features = pd.DataFrame()
#Features ES
features['es'] = es['Close']
features['es_vol'] = es['Volume']
features['es_oi'] = es['Open Interest']
features['es_delta_vol'] =  features['es_vol'] - features['es_vol'].shift(1)
features['es_delta_vol2'] = features['es_vol'] - features['es_vol'].shift(10)
features['es_delta_vol3'] = features['es_vol'] - features['es_vol'].shift(25)
features['es_range_oh'] = es['Open'] - es['High'].shift(1)
features['es_range_oc'] = es['Open'] - es['Close'].shift(1)
features['es_range_ol'] = es['Open'] - es['Low'].shift(1)
features['es_range_oo'] = es['Open'] - es['Open'].shift(1)
features['es_vol_dolar'] = es['Volume'] * es['Close'] * 50
features['es_dist_max'] = es['Close'].rolling(251).max() - es['Close'] 
features['es_dist_min'] = es['Close'].rolling(251).min() - es['Close']
features['es_autocorr'] = es['Close'].rolling(20).corr() - es['Close'].shift(1).rolling(20).corr()
features['es_autocorr2'] = es['Close'].rolling(20).corr() - es['Close'].shift(5).rolling(20).corr()

#Features US
features['us'] = y30['Close']
features['us_vol'] = y30['Volume']
features['us_vol2'] = y10['Volume']
features['us_oi'] = y30['Open Interest']
features['us_delta_vol'] =  features['us_vol'] - features['us_vol'].shift(1)
features['us_delta_vol2'] = features['us_vol'] - features['us_vol'].shift(10)
features['us_delta_vol3'] = features['us_vol'] - features['us_vol'].shift(25)
features['us_range_oh'] = y30['Open'] - y30['High'].shift(1)
features['us_range_oc'] = y30['Open'] - y30['Close'].shift(1)
features['us_range_ol'] = y30['Open'] - y30['Low'].shift(1)
features['us_range_oo'] = y30['Open'] - y30['Open'].shift(1)
features['us_vol_dolar'] = y30['Volume'] * y30['Close'] * 31.25
features['us_dist_max'] = y30['Close'].rolling(251).max() - y30['Close'] 
features['us_dist_min'] = y30['Close'].rolling(251).min() - y30['Close']
features['us_autocorr'] = y30['Close'].rolling(20).corr() - y30['Close'].shift(1).rolling(20).corr()
features['us_autocorr2'] = y30['Close'].rolling(20).corr() - y30['Close'].shift(5).rolling(20).corr()

#Features VIX
features['vx'] = vx['Close']
features['vx_vol'] = vx['Volume']
features['vx_oi'] = vx['Open Interest']
features['vx_delta_vol'] =  features['vx_vol'] - features['vx_vol'].shift(1)
features['vx_delta_vol2'] = features['vx_vol'] - features['vx_vol'].shift(10)
features['vx_delta_vol3'] = features['vx_vol'] - features['vx_vol'].shift(25)
features['vx_range_oh'] = vx['Open'] - vx['High'].shift(1)
features['vx_range_oc'] = vx['Open'] - vx['Close'].shift(1)
features['vx_range_ol'] = vx['Open'] - vx['Low'].shift(1)
features['vx_range_oo'] = vx['Open'] - vx['Open'].shift(1)
features['vx_vol_dolar'] = vx['Volume'] * vx['Close'] * 1000
features['vx_dist_max'] = vx['Close'].rolling(251).max() - vx['Close'] 
features['vx_dist_min'] = vx['Close'].rolling(251).min() - vx['Close']
features['vx_autocorr'] = vx['Close'].rolling(20).corr() - vx['Close'].shift(1).rolling(20).corr()
features['vx_autocorr2'] = vx['Close'].rolling(20).corr() - vx['Close'].shift(5).rolling(20).corr()


#Alternative Features

features['yield_curve_30_10'] = (y30['Close'].pct_change() - y10['Close'].pct_change())
features['y30v'] = np.where(features['us_vol'] > features['us_vol2'],1,0)
features['inter_corr1'] = round(features['es'].rolling(20).corr(features['us']) - features['es'].shift(2).rolling(20).corr(features['us']),4)
features['nyse_amp'] = newhighs['Numbers of Stocks'] - newlows['Numbers of Stocks']
features['nq_amp'] = nqhighs['Numbers of Stocks'] - nqlows['Numbers of Stocks']


#Alternative CBOE

features['vxst'] = vx9['Close']
features['vxst_vx'] = features['vxst'].rolling(9).corr(features['vx'])
features['vxst_vx2'] = features['vxst'].rolling(20).corr(features['vx'])
features['delta_vxst_vx'] = features['vxst_vx'] - features['vxst_vx2']
features = features.dropna()

#Alternative COT
features['us_dl'] = cot_30['Dealer Longs']
features['us_ds'] = cot_30['Dealer Shorts']
features['us_aml'] = cot_30['Asset Manager Longs']
features['us_ams'] = cot_30['Asset Manager Shorts']
features['us_lfl'] = cot_30['Leveraged Funds Longs']
features['us_lfl'] = cot_30['Leveraged Funds Shorts']
features['us_longs'] = cot_30['Total Reportable Longs'] + cot_30['Non Reportable Longs']
features['us_shorts'] = cot_30['Total Reportable Shorts'] + cot_30['Non Reportable Shorts']
features['vx_dl'] = cot_vx['Dealer Longs']
features['vx_ds'] = cot_vx['Dealer Shorts']
features['vx_aml'] = cot_vx['Asset Manager Longs']
features['vx_ams'] = cot_vx['Asset Manager Shorts']
features['vx_lfl'] = cot_vx['Leveraged Funds Longs']
features['vx_lfl'] = cot_vx['Leveraged Funds Shorts']
features['vx_longs'] = cot_vx['Total Reportable Longs'] + cot_vx['Non Reportable Longs']
features['vx_shorts'] = cot_vx['Total Reportable Shorts'] + cot_vx['Non Reportable Shorts']
features.fillna(method='ffill', inplace=True)

Al final hemos creado un dataset con 69 columnas, algo totalmente inviable, repito lo mencionado anteriormente, “Keep it simple”, Pocas features, y muy relevantes, pocas capas y pocas neuronas.

Al tener cada feature(columna) unos valores máximos y mínimos distintos, debemos normalizar los datos, es una condición fundamental, para ello utilizaremos, la librería sklearn, muy conocida en Python y su método MinMaxScaler()

In [30]:
features.shape
Out[30]:
(2169, 69)
In [18]:
X = features
X2 = X
X = MinMaxScaler().fit_transform(X)
X = pd.DataFrame(X,columns=X2.columns).set_index(X2.index)
X = X.dropna()
In [19]:
X
Out[19]:
eses_voles_oies_delta_voles_delta_vol2es_delta_vol3es_range_ohes_range_oces_range_oles_range_oo...us_lflus_longsus_shortsvx_dlvx_dsvx_amlvx_amsvx_lflvx_longsvx_shorts
Date
2011-02-010.1047160.3064330.7126140.3576900.5185180.6807500.8211620.4738810.2482850.573111...0.2376650.1665150.1665150.3298830.3501670.0000000.0040040.0786010.0692350.069235
2011-02-020.1033860.1918920.7202820.2776340.3669660.6118560.8062590.4664180.2524010.599607...0.2376650.1665150.1665150.3298830.3501670.0000000.0040040.0786010.0692350.069235
2011-02-030.1049580.2633240.7134290.4701700.3919060.6534100.7973170.4776120.1495200.514230...0.2376650.1665150.1665150.3298830.3501670.0000000.0040040.0786010.0692350.069235
2011-02-040.1068920.2232020.7178260.3546800.4233830.6149560.8137110.4738810.2016460.533857...0.2376650.1665150.1665150.3298830.3501670.0000000.0040040.0786010.0692350.069235
2011-02-070.1110040.2092360.7253010.3817590.4563520.6109890.8211620.4738810.1879290.537782...0.2376650.1665150.1665150.3298830.3501670.0000000.0040040.0786010.0692350.069235
..................................................................
2019-12-030.9696490.3159160.6874420.4337890.5970970.6392730.5678090.4738810.1591220.396467...0.7553080.7854940.7854940.6787570.1878680.8758660.8217690.5642180.6171460.617146
2019-12-040.9793230.2067290.6820600.2831770.5042860.5704080.6497760.4850750.2592590.433759...0.7553080.7854940.7854940.6787570.1878680.8758660.8217690.5642180.6171460.617146
2019-12-050.9825880.1791030.6806830.3676170.4292250.5384090.7779430.4738810.2976680.597645...0.7553080.7854940.7854940.6787570.1878680.8758660.8217690.5642180.6171460.617146
2019-12-060.9962520.1760480.6796290.3930560.4564010.4750860.7868850.4738810.2167350.548577...0.7553080.7854940.7854940.6787570.1878680.8758660.8217690.5642180.6171460.617146
2019-12-090.9906890.1504530.6714440.3697190.4884820.5098020.8017880.4813430.2990400.634936...0.7553080.7854940.7854940.6787570.1878680.8758660.8217690.5642180.6171460.617146

2168 rows × 69 columns

In [ ]:
 

Una vez que tenemos los datos homogeneizados, el dataset esta listo para enfrentarse a un proceso de PCA que le reducirá la dimensionalidad. Para ello utilizaremos una vez más, la librería sklearn.

Empezaremos, intentando reducir el dataset con 69 columnas, hasta 3 columnas.

Hemos conseguido un 55% de explicar el dataset con tan solo 3 columnas, es un gran dato, pero no es suficiente, En mi opinión, un proceso de reducción de dimensionalidad mediante PCA necesitamos un 85% de explicación, en el menor número de columnas posibles.

In [21]:
pca = PCA(n_components=3)
principalComponents = pca.fit_transform(X)
principalDf = pd.DataFrame(data = principalComponents
             , columns = ['pca1', 'pca2','pca3']).set_index(X.index)
print("shape of X_pca", principalDf.shape)
expl = pca.explained_variance_ratio_
print(expl)
print('suma:',sum(expl[0:16]))
shape of X_pca (2168, 3)
[0.28957657 0.17827025 0.08344746]
suma: 0.5512942765472438
In [22]:
principalDf
Out[22]:
pca1pca2pca3
Date
2011-02-01-0.942361-0.047399-0.564197
2011-02-02-0.946661-0.038974-0.574672
2011-02-03-0.938188-0.034657-0.596333
2011-02-04-0.920564-0.030289-0.656189
2011-02-07-0.925990-0.061899-0.661457
............
2019-12-031.2611140.5912680.360982
2019-12-041.2449760.5264960.253212
2019-12-051.2275440.5056290.219236
2019-12-061.2418500.4749190.169300
2019-12-091.1367550.5488480.186680

2168 rows × 3 columns

In [23]:
pca = PCA(n_components=4)
principalComponents = pca.fit_transform(X)
principalDf = pd.DataFrame(data = principalComponents
             , columns = ['pca1', 'pca2','pca3','pca4']).set_index(X.index)
print("shape of X_pca", principalDf.shape)
expl = pca.explained_variance_ratio_
print(expl)
print('suma:',sum(expl[0:16]))
shape of X_pca (2168, 4)
[0.28957657 0.17827025 0.08344746 0.06301154]
suma: 0.6143058188926562
In [25]:
pca = PCA(n_components=5)
principalComponents = pca.fit_transform(X)
principalDf = pd.DataFrame(data = principalComponents
             , columns = ['pca1', 'pca2','pca3','pca4','pca5']).set_index(X.index)
print("shape of X_pca", principalDf.shape)
expl = pca.explained_variance_ratio_
print(expl)
print('suma:',sum(expl[0:16]))
shape of X_pca (2168, 5)
[0.28957657 0.17827025 0.08344746 0.06301154 0.04705533]
suma: 0.6613611473195252
In [26]:
pca = PCA(n_components=8)
principalComponents = pca.fit_transform(X)
principalDf = pd.DataFrame(data = principalComponents
             , columns = ['pca1', 'pca2','pca3','pca4','pca5','pca6','pca7','pca8']).set_index(X.index)
print("shape of X_pca", principalDf.shape)
expl = pca.explained_variance_ratio_
print(expl)
print('suma:',sum(expl[0:16]))
shape of X_pca (2168, 8)
[0.28957657 0.17827025 0.08344746 0.06301154 0.04705533 0.03846761
 0.03172594 0.02789827]
suma: 0.7594529681820242
In [27]:
pca = PCA(n_components=10)
principalComponents = pca.fit_transform(X)
principalDf = pd.DataFrame(data = principalComponents
             , columns = ['pca1', 'pca2','pca3','pca4','pca5','pca6','pca7','pca8','pca9','pca10']).set_index(X.index)
print("shape of X_pca", principalDf.shape)
expl = pca.explained_variance_ratio_
print(expl)
print('suma:',sum(expl[0:16]))
shape of X_pca (2168, 10)
[0.28957657 0.17827025 0.08344746 0.06301154 0.04705533 0.03846764
 0.03172605 0.02789902 0.02613213 0.02342828]
suma: 0.8090142634638263
In [ ]:
 

Vemos como aumentando el numero de outputs que nos devuelve el proceso de PCA, va aumentando la explicación del dataset, al final debemos encontrar el equilibro entre la mayor explicación posible (el numero mas cercano a 100) con el menor numero de columnas (features) posibles.

Conclusiones

Y para concluir, recordamos, que la mejor combinacion es, la que tiene menos columnas, y aporta un mayor grado de explicacion del dataset

Como podemos apreciar en la imagen,a partir de las 6 dimensiones de nuestro PCA, deja de aportar en exceso, y es mas, añades dimensionalidad al dataset y no estas ganando demasiado poder explicativo, por consecuencia, es contraproducente.

Si te ha gustado este articulo,te recompendamos tambien puedes revisar los articulos anterioress seguir aprendiendo! Portfolios1 y Portfolios2