Testing Statistico della Mean Reversion – Parte II

Nel precedente articolo abbiamo introdotto il test statistico della mean-reverting. Abbiamo esaminato un paio di tecniche che ci hanno aiutato a determinare se una serie temporale fosse o meno mean reverting. In particolare abbiamo esaminato il test Dickey-Fuller e l’esponente di Hurst. In questo articolo considereremo un altro metodo, ovvero il test Cointegrated Dickey Fuller (CADF).

Innanzitutto, si deve notare che è molto difficile ,in realtà, trovare un asset direttamente negoziabile direttamente che abbia un comportamento mean revering. Ad esempio, le azioni si comportano generalmente come un moto browniano geometrico e quindi rendono relativamente inutili le strategie mean-reverting. Tuttavia, non c’è nulla che ci impedisca di creare un portafoglio di serie di prezzi stazionario. Quindi possiamo applicare strategie di trading a medio termine al portafoglio.

La forma più semplice di strategie di trading di mean-backing è il classico “pairs trade”, che di solito implica una coppia long-short di azioni indipendente dal moneta fiat. Questa tecnica si basa sulla teoria che due società dello stesso settore sono probabilmente esposte a simili fattori di mercato, che influenzano le loro attività. Occasionalmente i loro relativi prezzi azionari divergono a causa di determinati eventi, ma nel lungo periodo tengono a tornare verso la media.

Consideriamo due titoli del settore energetico Approach Resources Inc (AREX) e Whiting Petroleum Corp (WLL). Entrambi sono esposti alle stesse condizioni di mercato e quindi sono collegate in modo stazionario. Si può studiare le relazioni tra questa coppia di titoli usando le librerie Pandas e Matplotlib. Il primo grafico (Figura 1) mostra i rispettivi storici dei prezzi dal 1 ° gennaio 2012 al 1 ° gennaio 2013.

Fig 1 - Grafico dei prezzi di AREX e WLL
Se creiamo il grafico di dispersione dei loro prezzi, si osserva una relazione sostanzialmente lineare (vedi Figura 2) per questo periodo.
Fig 2 - Grafico della dispersione dei prezzi di AREX e WILL

Il pairs trading si basa su un modello lineare per identificare la relazione tra i prezzi di due azioni:

\(\begin{eqnarray}
y(t) = \beta x(t) + \epsilon(t)
\end{eqnarray}\)

Dove y(t) è il prezzo del titolo AREX e x(t) è il prezzo del titolo WLL, entrambi per il giorno t.

Se tracciamo i residui ε(t) = y(t) – βx(t) (per un particolare valore di β che determineremo in seguito) si crea una nuova serie temporale che, a prima vista, sembra relativamente stazionaria. Il grafico è mostrato in Figura 3:

Fig 3 - Grafico dei Residui della Regressione Lineare tra AREX e WLL

Test di Dickey-Fuller Aumetato e Cointegrato

Al fine di confermare statisticamente se questa serie è di tipo mean-reverting, si può utilizzare uno dei metodi descritti nell’articolo precedente, vale a dire il test Dickey-Fuller aumentato o l’esponente di Hurst. Tuttavia, nessuno di questi due test può determinare β, il rapporto di coperta (Hedge Ratio) necessario per formare la combinazione lineare, quindi possono solo verificare se, per un particolare β, la combinazione lineare è stazionaria.

È qui che entra in gioco il test Cointegrated Dickey-Fuller (CADF). Questo metodo determina il rapporto di correlazione ottimale eseguendo una regressione lineare rispetto alle due serie temporali e quindi verifica la stazionarietà di una combinazione lineare.

Implementazione in Python

La verifica della relazione di cointegrazione tra AREX e WLL per il periodo dal 1 ° gennaio 2012 al 1 ° gennaio 2013 può essere effettuata facilmente sfruttando le librerie di Python. Come fonte di dati per questo esempio è sufficiente utilizzare Yahoo Finance, mentre per eseguire il test ADF, come sopra, si utilizza la libreria Statsmodels. Il primo passo è creare un nuovo file, cadf.py, e importare le librerie necessarie. Il codice utilizza NumPy, Matplotlib, Pandas e Statsmodels. Per poter etichettare correttamente gli assi e scaricare i dati da Yahoo Finance tramite panda, importiamo il modulo matplotlib.dates e il modulo pandas.io.data. Usiamo anche la funzione dei minimi quadrati ordinari (OLS) presente in pandas:
#!/usr/bin/python
# -*- coding: utf-8 -*-

# cadf.py

import datetime
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import pandas as pd
import pandas.io.data as web
import pprint
import statsmodels.tsa.stattools as ts

from pandas.stats.api import ols

La prima funzione, plot_price_series, accetta un DataFrame panda come input, e due colonne identificate con stringhe “ts1” e “ts2”. Queste saranno le nostre coppie di azioni. Semplicemente la funzione traccia le due serie di prezzi sullo stesso grafico. Questo ci consente di verificare visivamente se esiste una probabile cointegrazione.

Usiamo il modulo dates di Matplotlib per ottenere i mesi dagli oggetti datetime. Quindi creiamo una figura e un insieme di assi su cui applicare l’etichettatura. Infine, stampiamo in output il grafico:

# cadf.py

def plot_price_series(df, ts1, ts2):
    months = mdates.MonthLocator()  # every month
    fig, ax = plt.subplots()
    ax.plot(df.index, df[ts1], label=ts1)
    ax.plot(df.index, df[ts2], label=ts2)
    ax.xaxis.set_major_locator(months)
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
    ax.set_xlim(datetime.datetime(2012, 1, 1), datetime.datetime(2013, 1, 1))
    ax.grid(True)
    fig.autofmt_xdate()

    plt.xlabel('Month/Year')
    plt.ylabel('Price ($)')
    plt.title('%s and %s Daily Prices' % (ts1, ts2))
    plt.legend()
    plt.show()
La seconda funzione, plot_scatter_series, traccia il grafico di dispersione dei due prezzi. Questo ci consente di verificare visivamente se esiste una relazione lineare tra le due serie e quindi se può essere un buon candidato per la procedura OLS e il successivo test ADF:
# cadf.py

def plot_scatter_series(df, ts1, ts2):
    plt.xlabel('%s Price ($)' % ts1)
    plt.ylabel('%s Price ($)' % ts2)
    plt.title('%s and %s Price Scatterplot' % (ts1, ts2))
    plt.scatter(df[ts1], df[ts2])
    plt.show()
La terza funzione, plot_residuals, è progettata per tracciare i valori residui del modello lineare adattato delle due serie di prezzi. Questa funzione richiede che il DataFrame abbia una colonna “res”, che rappresenta i prezzi residui:
# cadf.py

def plot_residuals(df):
    months = mdates.MonthLocator()  # every month
    fig, ax = plt.subplots()
    ax.plot(df.index, df["res"], label="Residuals")
    ax.xaxis.set_major_locator(months)
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
    ax.set_xlim(datetime.datetime(2012, 1, 1), datetime.datetime(2013, 1, 1))
    ax.grid(True)
    fig.autofmt_xdate()

    plt.xlabel('Month/Year')
    plt.ylabel('Price ($)')
    plt.title('Residual Plot')
    plt.legend()

    plt.plot(df["res"])
    plt.show()
Infine, la procedura è racchiusa in una funzione __main__. Il primo compito è scaricare i dati OHLCV per AREX e WLL da Yahoo Finance. Quindi creiamo un DataFrame separato, df, utilizzando lo stesso indice del frame AREX per archiviare entrambi i valori di prezzo di chiusura corretti. Quindi tracciamo la serie di prezzi e il grafico di dispersione. Dopo che i grafici sono completi, i residui vengono calcolati chiamando la funzione pandas ols sulle serie WLL e AREX. Questo ci consente di calcolare il rapporto β. Il rapporto di copertura viene quindi utilizzato per creare una colonna “res” tramite il calcolo della combinazione lineare di WLL e AREX. Infine, i residui vengono tracciati e il test ADF viene eseguito sui residui calcolati. Quindi stampiamo i risultati del test ADF:
# cadf.py

if __name__ == "__main__":
    start = datetime.datetime(2012, 1, 1)
    end = datetime.datetime(2013, 1, 1)

    arex = web.DataReader("AREX", "yahoo", start, end)
    wll = web.DataReader("WLL", "yahoo", start, end)

    df = pd.DataFrame(index=arex.index)
    df["AREX"] = arex["Adj Close"]
    df["WLL"] = wll["Adj Close"]

    # Plot the two time series
    plot_price_series(df, "AREX", "WLL")

    # Display a scatter plot of the two time series
    plot_scatter_series(df, "AREX", "WLL")

    # Calculate optimal hedge ratio "beta"
    res = ols(y=df['WLL'], x=df["AREX"])
    beta_hr = res.beta.x

    # Calculate the residuals of the linear combination
    df["res"] = df["WLL"] - beta_hr*df["AREX"]

    # Plot the residuals
    plot_residuals(df)

    # Calculate and output the CADF test on the residuals
    cadf = ts.adfuller(df["res"])
    pprint.pprint(cadf)
L’output del codice (ad esclusione dei grafici Matplotlib) è il seguente:
(-2.9607012342275936,
 0.038730981052330332,
 0,
 249,
 {'1%': -3.4568881317725864,
  '10%': -2.5729936189738876,
  '5%': -2.8732185133016057},
 601.96849256295991)
Si può vedere che il test statistico di -2.96 è inferiore al valore critico del 5% di -2.87, il che significa che possiamo rifiutare l’ipotesi che non ci sia una relazione di cointegrazione al livello del 5%. Quindi possiamo concludere, con un ragionevole grado di certezza, che AREX e WLL possiedono una relazione di cointegrazione, almeno nel il campione temporale considerato.

Perchè fare il test statistico?

Fondamentalmente, per quanto riguarda il trading algoritmico, i test statistici descritti sopra sono utili per analizzare i profitti che generano quando sono applicati a strategie di trading. Quindi, ha sicuramente senso valutare semplicemente le prestazioni a livello di strategia, al contrario del livello di serie prezzo / tempo? Perché prendersi il disturbo di calcolare tutte le metriche di cui sopra quando possiamo semplicemente utilizzare l’analisi a livello di trade, le metriche di rischio/rendimento e le valutazioni del drawdown?

In primo luogo, qualsiasi strategia di trading implementata basata su una misura statistica delle serie temporali avrà un campione molto più ampio su cui lavorare, dato che quando si calcolano questi test statistici, stiamo facendo uso delle informazioni contenute in ogni candela OHLCV, invece che di ogni trade. Ci saranno molti meno trade round-trip rispetto alle barre e quindi la rilevanza statistica di qualsiasi metrica del livello di trading sarà molto più piccola.

In secondo luogo, qualsiasi strategia che implementiamo dipenderà da alcuni parametri, come i periodi di osservazione per il calcolo dell medie o le misure z-score per entrare / uscire da una operazione, in un contesto di mean-reverting. Quindi le metriche del livello di strategia sono più appropriate per questi parametri, mentre i test statistici sono validi per il campione di serie temporali sottostante.

In pratica, vogliamo calcolare entrambe le serie di statistiche. Python, tramite le librerie pandas e statsmodels, rende questo estremamente semplice. Lo sforzo aggiuntivo è in realtà piuttosto minimo!

Testing Statistico della Mean Reversion

Finora su DataTrading abbiamo parlato dell’identificazione di strategie di trading algoritmico, il backtesting, i securities master database dei titoli e come creare un ambiente di sviluppo software. È giunto il momento di rivolgere la nostra attenzione alla creazione di strategie di trading e di come implementarle.

Uno dei concetti chiave nella cassetta degli attrezzi del trader quantitativo è quello della mean reversion. Questo concetto si riferisce ad una serie temporale che mostra una tendenza a ritornare verso valore medio. Matematicamente, una tale serie temporale (continua) viene definita processo Ornstein-Uhlenbeck. Questo è in contrasto con random walk (moto Browniano), dove non si prevede “memoria” dei valori pregressi in ogni specifica istanza di tempo. La proprietà mean-reverting di una serie temporale può essere sfruttata al fine di produrre strategie di trading redditizie.

In questo articolo verranno illustrati i test statistici necessari per identificare la mean reversion. In particolare, studieremo il concetto di stazionarietà e come testarla.

Testing per la Mean Reversion

Una serie temporale continua di mean-reverting può essere rappresentata da un’equazione differenziale stocastica di Ornstein-Uhlenbeck:
\(\begin{eqnarray}
d x_t = \theta (\mu – x_t) dt + \sigma dW_t
\end{eqnarray}\)

Dove θ è la velocità del processo di reversion, μ è il prezzo medio di lungo periodo verso cui avviene la reversione, σ è la varianza del processo e Wt è un processo Wiener o moto browniano.

L’equazione afferma che la “tendenza” delle serie di prezzi nel prossimo periodo di tempo è proporzionale alla differenza tra il prezzo medio e il prezzo corrente, con l’aggiunta del rumore gaussiano.

Questa proprietà è dimostata del test Augmented Dickey-Fuller, che descriveremo di seguito.

Il Test Dickey-Fuller Aumentato (ADF)

Matematicamente, l’ADF si basa sulla verifica della presenza di un trend o radici unitarie a root unit in un campione di serie temporali autoregressive. Sfrutta il fatto che se una serie di prezzi possiede una mean revertion, allora il livello di prezzo successivo sarà proporzionale al livello di prezzo corrente. 

Queste serie storiche sono rappresentate matimaticamente da un modello di regressione lineare di ordine p:

\(\begin{eqnarray}
\Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \delta_1 \Delta y_{t-1} + \cdots + \delta_{p-1} \Delta y_{t-p+1} + \epsilon_t
\end{eqnarray}\)

Dove α è una costante, β rappresenta il coefficiente di andamento temporale e Δyt = y (t) -y (t-1).

Il test ADF prevede che  il processo autogressivo di ordine p=1 abbia una media nulla, cioè γ = 0, che si traduce con α = β = 0 e quindi identifica il processo come casuale e non mean-reverting.

Se si può eliminare l’ipotesi che γ = 0, allora il movimento della serie di prezzi è proporzionale al prezzo corrente e quindi è improbabile che sia casuale.

Quindi, come viene eseguito il test ADF? Il primo compito è calcolare la statistica di prova \(DF_{\tau}\), che è data dalla costante di proporzionalità del campione \(\hat{\gamma}\) divisa per l’errore standard della costante di proporzionalità del campione:

\(\begin{eqnarray} DF_{\tau} = \frac{\hat{\gamma}}{SE(\hat{\gamma})} \end{eqnarray}\)

Dickey e Fuller hanno calcolato la distribuzione di questa statistica di test, che ci consente di scartare l’ipotesi per qualsiasi specifico valore percentuale scelto. La statistica del test è un numero negativo e quindi, per essere significativo, il numero deve essere più negativo di questi valori, cioè inferiore ai valori critici.

In sostanza, per i trader è che qualsiasi deriva a lungo termine in termini di  prezzo è molto probabile rispetto a qualsiasi fluttuazione a breve termine e quindi il modello approssima a zero la deriva (β = 0).

Poiché stiamo considerando un modello di regressione di ordine p, è necessario definire p con uno specifico valore. Di solito è sufficiente, per le strategie di trading, impostare p = 1 per consentirci scartare l’ipotesi nulla.

Per calcolare il test Dickey-Fuller aumentato possiamo utilizzare le librerie pandas e statsmodels. La prima ci fornisce un metodo diretto per ottenere i dati Open-High-Low-Close-Volume (OHLCV) da Yahoo Finance, mentre il secondo implementa il test ADF in una funzione facile da richiamare.

Effettueremo il test ADF su i dati storici delle azioni di Google, dal 1 ° gennaio 2000 al 1 ° gennaio 2013.

Serie dei prezzi Google dal 01/01/2000 to 01/01/2013
Di seguito il codice Python per eseguire il test:
# Import the Time Series library
import statsmodels.tsa.stattools as ts

# Import Datetime and the Pandas DataReader
from datetime import datetime
from pandas.io.data import DataReader

# Download the Google OHLCV data from 1/1/2000 to 1/1/2013
goog = DataReader("GOOG", "yahoo", datetime(2000,1,1), datetime(2013,1,1))

# Output the results of the Augmented Dickey-Fuller test for Google
# with a lag order value of 1
ts.adfuller(goog['Adj Close'], 1)
Ecco il risultato del test di Dickey-Fuller Aumentato per le azioni Google nel periodo considerato. Il primo valore è il valore della statistica test, il secondo valore è il p-value. Il quarto è il numero di punti nel campione. Il quinto valore, il dizionario, contiene i valori critici della statistica del test rispettivamente ai valori 1, 5 e 10%.
(-2.1900105430326064,
 0.20989101040060731,
 0,
 2106,
 {'1%': -3.4334588739173006,
  '10%': -2.5675011176676956,
  '5%': -2.8629133710702983},
 15436.871010333041)
Poiché il valore della statistica test è maggiore di uno qualsiasi dei valori critici, sui livelli 1, 5 o 10 percento, non possiamo respingere l’ipotesi nulla di γ = 0 e quindi è improbabile che abbiamo trovato una serie temporale mean reverting. Un metodo alternativo per identificare una serie temporale mean reverting sfrutta il concetto di stazionarietà, di cui parleremo ora.

Testing per la Stazionarità

Una serie temporale (o processo stocastico) è definita fortemente stazionaria se la sua distribuzione di probabilità congiunta non cambia se viene traslata nel tempo o nello spazio. In particolare, e di fondamentale importanza per i trader, la media e la varianza del processo non cambiano nel tempo o nello spazio e nessuno di essi stia seguendo un trend. Una caratteristica fondamentale delle serie di prezzi stazionarie è che i prezzi all’interno della serie si allontanano dal loro valore iniziale più lentamente rispetto al moto browniano geometrico. Misurando il tasso di questo comportamento possiamo identificare la natura delle serie temporali. Descriveremo ora un parametro, ovvero l’Esponente di Hurst, che ci aiuta a caratterizzare la stazionarietà di una serie temporale.

Esponente Hurst

L’Espondente di Hurst Exponent rappresenta un valore scalare che può essere d’aiuto ad identificare (entro i limiti della stima statistica) se una serie è mean reverting, random walking o trending. L’idea alla base del calcolo dell’esponente consiste nell’usare la varianza dei logaritmi di una serie di prezzi per valutare il tipo di comportamento della struttura. Per un intervallo di tempo arbitrario τ, la varianza è data da:
\(\begin{eqnarray} {\rm Var}(\tau) = \langle |\log(t+\tau)-\log(t)|^2 \rangle \end{eqnarray}\)
Dal momento che si sta confrontando il tasso di diffusione con quello di un moto browniano geometrico, si può dire che τ è abbastanza grande e quindi la varianza è proporzionale a τ, come nel caso di un GBM:
\(\begin{eqnarray} \langle |\log(t+\tau)-\log(t)|^2 \rangle \sim \tau \end{eqnarray} \)
L’idea chiave è che se esistono alcune autocorrelazioni (ad esempio un qualsiasi movimento sequenziale del prezzo con una correlazione diversa da zero), la relazione di cui sopra non è valida. Ciò nonostate, può essere modificato per includere il valore esponenziale “2H”, da cui calcolare il valore H dell’Esponente di Hurst:
\(\begin{eqnarray} \langle |\log(t+\tau)-\log(t)|^2 \rangle \sim \tau^{2H} \end{eqnarray} \)

Una serie temporale può quindi essere caratterizzata nel seguente modo:

  • H <0,5 – Le serie temporali sono mean reverting;
  • H = 0.5 – Le serie temporali sono in un moto browniano geometric;
  • H> 0,5 – Le serie temporali sono in trend

Oltre alla caratterizzazione delle serie temporali, l’esponente di Hurst descrive anche il grado in cui una serie si comporta nel modo categorizzato. Ad esempio, un valore di H vicino allo 0 vuol dire che una serie ha un’elevata mean-reverting, mentre per H vicino a 1 la serie è fortemente in trend.

Per calcolare l’esponente di Hurst per la serie di prezzi di Google, gli stessi ustati come esempio per l’ADF, possiamo usare il seguente codice Python:

from numpy import cumsum, log, polyfit, sqrt, std, subtract
from numpy.random import randn

def hurst(ts):
	"""Returns the Hurst Exponent of the time series vector ts"""
	# Create the range of lag values
	lags = range(2, 100)

	# Calculate the array of the variances of the lagged differences
	tau = [sqrt(std(subtract(ts[lag:], ts[:-lag]))) for lag in lags]

	# Use a linear fit to estimate the Hurst Exponent
	poly = polyfit(log(lags), log(tau), 1)

	# Return the Hurst exponent from the polyfit output
	return poly[0]*2.0

# Create a Gometric Brownian Motion, Mean-Reverting and Trending Series
gbm = log(cumsum(randn(100000))+1000)
mr = log(randn(100000)+1000)
tr = log(cumsum(randn(100000)+1)+1000)

# Output the Hurst Exponent for each of the above series
# and the price of Google (the Adjusted Close price) for 
# the ADF test given above in the article
print "Hurst(GBM):   %s" % hurst(gbm)
print "Hurst(MR):    %s" % hurst(mr)
print "Hurst(TR):    %s" % hurst(tr)

# Assuming you have run the above code to obtain 'goog'!
print "Hurst(GOOG):  %s" % hurst(goog['Adj Close'])
Mandando in esecuzione il code Python si ottiene il seguente output:
Hurst(GBM):   0.500606209426
Hurst(MR):    0.000313348900533
Hurst(TR):    0.947502376783
Hurst(GOOG):  0.50788012261

Da questa output si nota come il movimento browniano geometrico possiede un esponente di Hurst, H, che è quasi esattamente 0,5. La serie di mean reverting ha H quasi uguale a zero, mentre la serie in trend ha H vicino a 1.

È interessante notare che Google ha anche H quasi uguale a 0,5 che indica che è estremamente simile ad una random walk geometrica (almeno nel periodo considerato).

Ora che abbiamo un metodo per caratterizzare la natura di una serie temporale di prezzi, dobbiamo discutere quanto il valore H sia statisticamente significativo. Dobbiamo essere in grado di determinare se sia possiamo rifiutare l’ipotesi nulla che H = 0.5 in modo da valutare un comportamento mean reverting o di trend.

Negli articoli successivi descriveremo come calcolare che H è statisticamente significativa. Inoltre, prenderemo in considerazione il concetto di cointegrazione, che ci consentirà di creare le nostre serie temporali di mean-reverting da più serie di prezzi differenti. Infine, uniremo insieme queste tecniche statistiche al fine di formare una strategia di trading mean reverting di base.