Ya hablamos del contenido de memoria y las propiedades estadísticas de las series temporales fraccionalmente diferenciadas (fracdiff en corto) en esta otra publicación. Ahora vamos a comparar el desempeño de las series de tiempo fracdiff como predictoras del futuro con predicciones de línea base usando series de precio limpias, sin diferenciar, o diferenciadas cuando d=1 (la primera derivada de toda la vida, generalmente los retornos diarios).
Vamos a utilizar la premisa en esta otra publicación, esto es, el comportamiento pasado del índice SPY puede predecir o adelantar al valor del indicador VIX. Vamos a asumir también que podemos obtener estas predicciones mediante regresión y utilizaremos métodos comúnmente utilizados para generar señales de "arriba" o "abajo" como predicción varios días hacia el futuro.
El primer paso es obtener los precios requeridos y el valor de VIX utilizando los métodos de agregado de instrumentos y datos de Quantconnect:
qb = QuantBook()
spy = qb.AddEquity('SPY').Symbol
vix = qb.AddData(CBOE, 'VIX').Symbol
symbol_list = [spy, vix]
#Dates for training the model:
start_date, end_date = datetime(2005, 1, 1), datetime(2015, 1, 1)
history = qb.History(qb.Securities.Keys, start_date, end_date, Resolution.Daily)El dataframe que se genera ha de ser reindexado y pivotado, solo vamos a utilizar el precio de cierre ("close") esta vez, ignoraremos el volumen y los diferenciales de compraventa que podrían, porqué no, ser utilizados en un modelo más complejo:
reindex_history = history.reset_index()
reindex_history['time'] = reindex_history['time'].dt.date
close_prices = reindex_history.pivot(index = 'time', columns ='symbol', values = 'close').dropna()Sabemos por análisis anteriores que el precio de SPY no es estacionario durante este periodo histórico. En cambio, sabemos que los valores de cierre del VIX en este periodo sí que son estacionarios. En cualquier caso vamos a verificar esta estacionaridad para las diferenciales fraccionarias entre 0 y 1, desde la serie original hasta la primera derivada. Utilizamos la biblioteca de MLFinlab de nuevo. MLFinlab ha cambiado su estrategia en cuanto al desarrollo de la biblioteca, seguimos creyendo que se trata de una biblioteca muy buena y útil y mantenemos nuestro patronaje dentro de nuestras (modestas) posibilidades:
symbol_list = [spy, vix]
fracdiff_series_dictionary = {}
fracdiff_series_adf_tests = {}
fracdiff_values = np.linspace(0,1,11)
for frac in fracdiff_values:
for symbol in symbol_list:
frac = round(frac,1)
dic_index = (symbol, frac)
fracdiff_series_dictionary[dic_index] = frac_diff_ffd(close_prices[[symbol]], frac, thresh=0.0001)
fracdiff_series_dictionary[dic_index].dropna(inplace=True)
fracdiff_series_adf_tests[dic_index] = adfuller(fracdiff_series_dictionary[dic_index])Nos quedaremos solo con aquellas series que son estacionarias. Las consignamos a diccionarios en caso de querer, en el futuro, utilizar símbolos adicionales o datos, o simplemente para facilitar la inspección de estas series.
stationary_fracdiffs = [key for key in fracdiff_series_adf_tests.keys() if fracdiff_series_adf_tests[key][1] < 0.05 ]
first_stationary_fracdiff = {}
diff_series = {}
for symbol, d in stationary_fracdiffs:
ds = [kvp for kvp in stationary_fracdiffs if kvp[0]==symbol]
first_stationary_fracdiff[symbol] = sorted(ds , key = lambda x: x[1], reverse=False)[0][1]
diff_series[symbol] = fracdiff_series_dictionary[(symbol,first_stationary_fracdiff[symbol])]
diff_series[symbol].columns =['diff_close_'+str(symbol)]
for symbol in symbol_list:
print("The lowest stationary fractionally differentiated series for ", symbol," is for d=", first_stationary_fracdiff[symbol],".")The lowest stationary fractionally differentiated series for SPY is for d= 0.5 .
The lowest stationary fractionally differentiated series for VIX.CBOE is for d= 0.0 .Comprobamos que, efectivamente, para este periodo el VIX es estacionario en d=0 y SPY no es estacionario hasta d=0.5. Ahora estas dos series deberían exhibir propiedades estadísticas decentes que permitirán su uso como herramientas de predicción en modelos de aprendizaje automático y como generadores de predicciones. Construiremos nuestra predicción objetivo como movimiento para N "días en el futuro", Cierto hacia arriba, Falso hacia abajo. De esta forma nuestro modelo de predicción se reduce a nada más que (como si fuera poco) una clasificación binaria:
prediction_target_days = [1,3,5,7,13,23]
predicted_symbol = vix
predicting_symbol = spy
data = close_prices
data['SPY_1D_Return'] = fracdiff_series_dictionary[(spy,1)]
data['VIX_1D_Return'] = fracdiff_series_dictionary[(vix,1)]
for symbol in symbol_list:
data = data.merge(diff_series[symbol],left_index=True, right_index=True)
for future in prediction_target_days:
data[('fut_dir_'+str(future)+'_'+str(symbol))] = data[[symbol]].shift(-future) > data[[symbol]]
data.dropna(inplace=True)Separamos estos datos en conjuntos de aprendizaje y prueba (train y test). Utilizamos un conjunto de prueba relativamente pequeño confiando en que el uso de validación cruzada durante el entrenamiento y el trabajo en el futuro (siguiente publicación) nos sirvan de validación efectiva del modelo:
# Get sizes of each of the datasets
test_size = 0.1
num_test = int(test_size*len(data))
num_train = len(data)- num_test
print("Training data points: %d " %num_train)
print("Test data points: %d" % num_test)
# Split into train, cv, and test
train = data[:num_train]
test = data[num_train:]
print("Training DataFrame Shape: ",train.shape)
print("Testing DataFrame Shape: ",test.shape)Training data points: 2086
Test data points: 231
Training DataFrame Shape: (2086, 17)
Testing DataFrame Shape: (231, 17)Ahora podemos escalar los conjuntos de aprendizaje y prueba. Recordemos que se deben ajustar los escaladores a los datos de aprendizaje primero y aplicar esta transformación a los datos de prueba, de otro modo estaríamos filtrando datos hacia el conjunto de pruebas. Utilizaremos un escalado de maximos y minimos y otro escalado de estandarización en secuencia. La función se vuelve monstruosa al tener que preservar los factores de escalado para el futuro en forma de diccionario al que poder acudir en caso de necesidad:
from sklearn.preprocessing import StandardScaler, MinMaxScaler
def scale_data(selected_features, selected_target, predict_window, train, test):
min_max_scaler_parameters = {}
standard_scaler_parameters = {}
scaled_train_features = pd.DataFrame()
scaled_test_features = pd.DataFrame()
for feature in selected_features:
min_max_scaler = MinMaxScaler(feature_range=(-1,1))
min_max_scaler_parameters[feature] = min_max_scaler.fit(np.array(train[feature]).reshape(-1,1))
scaled_train_features[str(feature)+"_scaled"] = min_max_scaler.transform(np.array(train[feature]).reshape(-1,1)).flatten()
standard_scaler = StandardScaler()
standard_scaler_parameters[feature] = standard_scaler.fit(np.array(scaled_train_features[str(feature)+"_scaled"]).reshape(-1,1))
scaled_train_features[str(feature)+"_scaled"] = standard_scaler.transform(np.array(scaled_train_features[str(feature)+"_scaled"]).reshape(-1,1)).flatten()
if test_size == 0.0: continue
scaled_test_features[str(feature)+"_scaled"] = min_max_scaler.transform(np.array(test[feature]).reshape(-1,1)).flatten()
scaled_test_features[str(feature)+"_scaled"] = standard_scaler.transform(np.array(scaled_test_features[str(feature)+"_scaled"]).reshape(-1,1)).flatten()
X_train = np.array(scaled_train_features)
X_test = np.array(scaled_test_features)
y_train = train['fut_dir_'+str(predict_window)+'_'+str(selected_target)]
y_test = test['fut_dir_'+str(predict_window)+'_'+str(selected_target)]
return X_train, X_test, y_train, y_test
Tenemos ahora las entradas para cada dia con la derivada fraccional del cierre de SPY a d=0.5, el cambio diario de SPY y el valor actual del VIX. Podemos utilizar estos datos para predecir la dirección futura y para ello alimentaremos nuestro modelo automático con los puntos de datos de manera secuencial conformando ventanas de tiempo de varios días. Nuestros datos de aprendizaje incluyen una serie de datos de "mirada atrás" de nuestros valores de predicción, conteniendo tantos puntos previos como queramos. Podemos generar estas entradas de predicción de manera flexible para obtener un análisis más robusto:
def input_features(data, lookback_window, predict_window):
nb_samples = len(data) - lookback_window - predict_window
input_list = [np.expand_dims(data[i:lookback_window+i,:], axis=0) for i in range(nb_samples)] #here nb_samples comes in handy
input_mat = np.concatenate(input_list, axis=0)
return input_mat
def output_targets(data, lookback_window, predict_window):
nb_samples = len(data) - lookback_window - predict_window
input_list = [np.expand_dims(data[i:lookback_window+i][-1], axis=0) for i in range(nb_samples)] #here nb_samples comes in handy
input_mat = np.concatenate(input_list, axis=0)
return input_mat
def generate_windows(X_train, X_test, y_train, y_test, lookback_window, predict_window, predict_symbol):
#format the features into the batch size
X_train_windows = input_features(X_train, lookback_window, predict_window).reshape(-1,lookback_window*len(selected_features))
X_test_windows = input_features(X_test, lookback_window, predict_window).reshape(-1,lookback_window*len(selected_features))
y_train = np.array(train['fut_dir_'+str(predict_window)+"_"+str(predict_symbol)])
y_test = np.array(test['fut_dir_'+str(predict_window)+"_"+str(predict_symbol)])
y_train_windows = output_targets(y_train, lookback_window, predict_window).astype(int)
y_test_windows = output_targets(y_test, lookback_window, predict_window).astype(int)
return X_train_windows, X_test_windows, y_train_windows, y_test_windowsExiste una doble transformación de forma (.reshape) en la definición de las matrices de mirada atrás. Esta ventana hacia atrás puede ser utilizada para alimentar otro modelo LSTM que requiere matrices tridimensionales, sería una entrada válida para LSTM con pocos cambios. Con los datos en posición comprobamos las características diferenciadas con una batería de modelos de ojeo para conseguir una indicación rápida de qué clasificador (dentro de la oferta básica en sklearn) ofrece el mejor rendimiento inicial:
from sklearn.linear_model import LogisticRegression
from sklearn.neural_network import MLPClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn import svm
from sklearn.model_selection import cross_val_score
def spot_models(X_train_matrix, y_train_matrix):
models = []
models.append(('Logistic Regression', LogisticRegression(solver = 'liblinear')))
models.append(('K-Nearest Neighbours', KNeighborsClassifier()))
models.append(('Random Forest', RandomForestClassifier(n_estimators = 10)))
models.append(('Support Vector Classifier', svm.SVC(gamma='auto')))
models.append(('Multi-layer Perceptron', MLPClassifier(max_iter = 1000, solver = 'sgd', learning_rate='adaptive', hidden_layer_sizes=(3))))
results = []
names = []
tscv = TimeSeriesSplit(n_splits=10)
for name, model in models:
cv_results = cross_val_score(model, X_train_matrix, y_train_matrix, cv=tscv, scoring='roc_auc')
results.append(cv_results)
names.append(name)
print('%s: Mean Acc.: %.2f std: %.2f - min: %.2f - max: %.2f' % (name, cv_results.mean(), cv_results.std(),cv_results.min(), cv_results.max()))
# Compare Algorithms
plt.boxplot(results, labels=names)
plt.title('Algorithm Accuracy Comparison')
plt.show()Utilizamos validación cruzada del tipo "series temporales" ya que encaja bien con el tipo de datos.
Ajustamos nuestros modelos en secuencia a una regresión logística, un perceptrón multicapa, un clasificador de K vecinos más próximos, un bosque aleatorio y una máquina clasificadora de vectores de soporte. No nos vamos a preocupar por los detalles de los modelos, solamente nos interesa una visión global de su comportamiento y una dirección general que seguir. Ajustamos después los modelos y evaluamos su puntuación con todos los datos disponibles; precios estacionarios de SPY, precio actual de SPY, valor actual de VIX y los retornos de un dia de ambos valores:
selected_features = ['diff_close_SPY', vix, 'SPY_1D_Return', spy, 'VIX_1D_Return']
lookback_window = 22
predict_window = 5
selected_target = vix
X_train, X_test, y_train, y_test = scale_data(selected_features, selected_target, predict_window, train, test)
X_train_windows, X_test_windows, y_train_windows, y_test_windows = generate_windows(X_train, X_test, y_train, y_test, lookback_window, predict_window, predict_symbol)
spot_models(X_train_windows, y_train_windows)La puntuación de precisión para los modelos es (utilizando la curva de característica operativa del receptor, roc_auc, nombre tan largo por venir del desarrollo de radares militares):
Regresión Logística: Prec. Promedio: 0.56, dev. est.: 0.07 - min: 0.47 - max: 0.68 K-Vecinos: Prec. Promedio: 0.52, dev. est.: 0.04 - min: 0.46 - max: 0.62 Bosque Aleatorio: Prec. Promedio: 0.55, dev. est.: 0.07 - min: 0.42 - max: 0.67 Clas. Soporte-Vector: Prec. Promedio: 0.56, dev. est.: 0.06 - min: 0.46 - max: 0.70 Perceptrón Multicapa: Prec. Promedio: 0.56, dev. est.: 0.07 - min: 0.45 - max: 0.71
En un formato de gráficos de caja, que resulta más accesible:
Para precisiones promedio del 55% el modelo de vectores de soporte ofrece el mayor potencial entre los modelos clasificadores ojeados. Seleccionaremos el modelo SVC como nuestra línea base e intentaremos mejorar el rendimiento del modelo mediante una busca de parámetros en parrilla, básicamente un método de fuerza bruta sin más complicaciones:
from sklearn.model_selection import GridSearchCV
model = svm.SVC()
param_search = {
'kernel': ['linear', 'rbf', 'sigmoid'],
'gamma': [1e-2, 1e-3],
'class_weight' : ['balanced'],
'C': [1, 10, 100],
'decision_function_shape': ['ovo', 'ovr']
}
tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = 'roc_auc', verbose=2)
gsearch.fit(X_train_windows, y_train_windows)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_La opción de verbosidad (una palabra muy sonora, verbose=2 en el código) en la búsqueda de parrilla de sklearn genera una buena cantidad de datos y nos da una indicación de los posibles tiempos de aprendizaje del modelo para cada uno de los parámetros en el caso de que queramos profundizar en la selección de manera más refinada. Este es el resultado parcial de la búsqueda de parrilla:
Fitting 10 folds for each of 36 candidates, totalling 360 fits
[CV] C=1, class_weight=balanced, decision_function_shape=ovo, gamma=0.01, kernel=linear
[CV] C=1, class_weight=balanced, decision_function_shape=ovo, gamma=0.01, kernel=linear, total= 0.0s
[CV] C=1, class_weight=balanced, decision_function_shape=ovo, gamma=0.01, kernel=linear
[CV] C=1, class_weight=balanced, decision_function_shape=ovo, gamma=0.01, kernel=linear, total= 0.0s
[CV] C=1, class_weight=balanced, decision_function_shape=ovo, gamma=0.01, kernel=linear
[CV] C=1, class_weight=balanced, decision_function_shape=ovo, gamma=0.01, kernel=linear, total= 0.1s
[CV] C=1, class_weight=balanced, decision_function_shape=ovo, gamma=0.01, kernel=linear Después de encontrar el mejor modelo posible de la parrilla de parámetros es momento de comprobar el modelo con los datos de prueba reservados:
print('Best Model Score:', best_score.round(2))
print('Best Model Parameters:', best_model)
validation_score = best_model.score(X_test_windows, y_test_windows).round(2)
print('Best Model Validation Score:', validation_score)Puntuación Mejor Modelo: 0.58 Parámetros Mejor Modelo: SVC(C=1, cache_size=200, class_weight='balanced', coef0=0.0, decision_function_shape='ovo', degree=3, gamma=0.001, kernel='rbf', max_iter=-1, probability=False, random_state=None, shrinking=True, tol=0.001, verbose=False) Puntuación de Validación: 0.55
Hemos sido capaces de clasificar correctamente el 55% de los puntos de prueba. Los datos de prueba son 104 Verdaderos y 100 Falsos, estamos ligeramente por encima de un clasificador ingenuo. Investiguemos finalmente si nuestras características diferenciadas han sido importantes a la hora de obtener estos resultados eliminandolas de los datos de aprendizaje y repitiendo el proceso completo:
Puntuación Mejor Modelo: 0.58 Parámetros Mejor Modelo: SVC(C=1, cache_size=200, class_weight='balanced', coef0=0.0, decision_function_shape='ovo', degree=3, gamma=0.001, kernel='sigmoid', max_iter=-1, probability=False, random_state=None, shrinking=True, tol=0.001, verbose=False) Puntuación Validación Mejor Modelo: 0.49
Las características diferenciadas definitivamente impulsan la puntuación del modelo. El núcleo favorito en ausencia de las series diferenciadas es el sigmoide en vez de rbf, indicando que la característica diferenciada puede contener la mayor parte de la información útil en el modelo anterior constituyendo una distribución gaussiana que se pierde sin esta.
Parece que utilizar diferenciación fraccional para generar los modelos, que sospechábamos útil, mejora los resultados de nuestra predicción. En nuestra próxima publicación implementaremos el mejor modelo de esta investigación en una validación de prueba retrospectiva y evaluaremos la rentabilidad del poder predictivo del modelo.
Por favor recuerde que las publicaciones en Ostirion.net no son consejos financieros. Ostirion.net no mantiene posiciones en ninguno de los instrumentos financieros que se mencionan en esta publicación en el momento de la publicación. Si necesita más información, apoyo con la gestión de activos financieros, desarrollo de estrategias de trading automatizados o despliegue táctico de estrategias existentes no dude en contactar con nosotros aquí.