Nella precedente serie di articoli (Parti 1 , 2 e 3 ) abbiamo approfondito in modo significativo i modelli di serie temporali lineari AR(p), MA(q) e ARMA(p,q). Abbiamo utilizzato questi modelli per generare set di dati simulati, modelli adattati per recuperare i parametri e quindi applicare questi modelli ai dati delle azioni finanziarie.
In questo articolo introduciamo un’estensione del modello ARMA, ovvero il modello Autoregressive Integrated Moving Average, o modello ARIMA(p,d,q). Vediamo che è necessario considerare il modello ARIMA quando abbiamo serie non stazionarie. Tali serie si verificano in presenza di tendenze stocastiche.
Riepilogo rapido e prossimi passi
Ad oggi abbiamo considerato i seguenti modelli (i link vi porteranno agli articoli appropriati):
- Rumore bianco discreto
- Passeggiata casuale
- Modello autoregressivo di ordine p, AR(p)
- Modello a media mobile di ordine q, MA(q)
- Modello autoregressivo a media mobile di ordine p e q, ARMA(p, q)
Abbiamo costantemente sviluppato lo studio delle serie temporali con concetti come correlazione seriale, stazionarietà, linearità, residui, correlogrammi, simulazione, adattamento, stagionalità, eteroschedasticità condizionale e verifica delle ipotesi.
Finora non abbiamo effettuato alcuna previsione o predizione dai nostri modelli e quindi non abbiamo avuto alcun meccanismo per produrre un sistema di trading o una curva di equity.
Dopo aver studiato ARIMA (in questo articolo), ARCH e GARCH (nei prossimi articoli), saremo in grado di costruire una strategia di trading di base a lungo termine basata sulla previsione dei rendimenti degli indici del mercato azionario.
Nonostante il fatto che abbiamo descritto dettagliatamente i modelli base (AR, MA, ARMA) che non hanno grandi prestazioni nel trading, ora siamo esperti nel processo di modellazione delle serie temporali.
Questo significa che possiamo studiare i modelli più recenti (e anche quelli attualmente presenti nella letteratura accademica) con una significati base di conoscenza a cui attingere per valutare efficacemente questi modelli, invece di trattarli come una “svolta chiave” o “scatola nera”.
Ancora più importante, ci darà la sicurezza di estenderli e modificarli da soli e capire cosa stiamo facendo quando lo facciamo!
Vorrei ringraziarvi per essere stati pazienti finora, poiché potrebbe sembrare che questi articoli siano lontani dalla “vera azione” del trading reale. Tuttavia, la vera ricerca sul trading quantitativo è attenta, misurata e richiede molto tempo per essere corretta. Non esiste una soluzione rapida o uno “schema di arricchimento” nel trading quantitativo.
Siamo quasi pronti per prendere in considerazione il nostro primo modello di trading, che sarà un misto di ARIMA e GARCH, quindi è imperativo dedicare un po’ di tempo alla comprensione del modello ARIMA!
Una volta costruito il nostro primo modello di trading, prenderemo in considerazione modelli più avanzati come i processi a memoria lunga, i modelli dello spazio degli stati (ad esempio il filtro di Kalman) e i modelli Vector Autoregressive (VAR), che ci porteranno verso strategie di trading più sofisticate.
Modelli Autoregressivi Integrati a Media Mobile (ARIMA) di ordine p, d, q
Fondamento logico
I modelli ARIMA vengono utilizzati perché possono ridurre una serie non stazionaria a una serie stazionaria utilizzando una sequenza di passaggi differenziati.
Come descritto nell’articolo sul rumore bianco e le passeggiate casuali, se applichiamo l’operatore differenziale a una serie di passeggiate casuali \(\{x_t \}\) (una serie non stazionaria) rimaniamo con il rumore bianco \(\{w_t \}\) (una serie stazionaria):
\(\begin{eqnarray} \nabla x_t = x_t – x_{t-1} = w_t \end{eqnarray}\)
L’ARIMA esegue essenzialmente questa funzione, ma lo fa ripetutamente, d volte, per ridurre una serie non stazionaria a una stazionaria.
Per gestire altre forme di non stazionarietà, oltre le tendenze stocastiche possono essere utilizzati modelli aggiuntivi.
Gli effetti della stagionalità (come quelli che si verificano nei prezzi delle materie prime) possono essere affrontati con il modello ARIMA stagionale (SARIMA), mentre gli effetti eteroscedastici condizionali (come il raggruppamento della volatilità negli indici azionari) possono essere affrontati con i modelli ARCH/GARCH.
In questo articolo consideriamo le serie non stazionarie con tendenze stocastiche e adattiamo i modelli ARIMA a queste serie. Infine, otteniamo le previsioni per le serie finanziarie.
Definizioni
Prima di definire i processi ARIMA dobbiamo discutere il concetto di serie integrata:
Serie integrata di ordine d
Una serie temporale \(\{x_t \}\) è integrata di ordine \(I\), \(I(d)\), Se:\(\begin{eqnarray}\nabla^d x_t = w_t \end{eqnarray}\)
Cioè, se differenziamo la serie d volte riceviamo una serie discreta di rumore bianco.
In alternativa, utilizzando l’operatore di spostamento all’indietro B una condizione equivalente è:
\(\begin{eqnarray}(1-{\bf B}^d) x_t = w_t \end{eqnarray} \)
Ora che abbiamo definito una serie integrata possiamo definire il processo ARIMA stesso:
Modello Autoregressivo Integrato a Media Mobile di ordine p, d, q
Una serie temporale \(\{x_t \}\) è un modello autoregressivo integrato a media mobile di ordine p, d, q, ARIMA(p,d,q) , se \(\nabla^d x_t\) è una media mobile autoregressiva di ordine p, q, ARMA(p,q).Cioè, se la serie \(\{x_t \}\) è differenziata d volte e segue un processo ARMA(p,q), allora è una serie ARIMA(p,d,q).
Se utilizziamo la notazione polinomiale della Parte 1 e della Parte 2 della serie ARMA, è possibile scrivere un processo ARIMA(p,d,q) in termini dell’ operatore di spostamento all’indietro, B:
\(\begin{eqnarray}\theta_p({\bf B})(1-{\bf B})^d x_t = \phi_q ({\bf B}) w_t\end{eqnarray}\)
Dove \(w_t\) è una serie discreta di rumore bianco.
Ci sono alcuni punti da notare su queste definizioni.
Dato che la passeggiata casuale è data da \(x_t = x_{t-1} + w_t\) si può vedere che \(I(1)\) è un’altra rappresentazione, dal momento che \(\nabla^1 x_t = w_t\).
Se sospettiamo una tendenza non lineare, siamo in grado di utilizzare differenziazioni ripetute (es \(d > 1\)) per ridurre una serie a rumore bianco stazionario.
In Python possiamo utilizzare il comando diff
della libreria Pandas comando con parametri aggiuntivi, ad esempio l’intervallo, per eseguire differenziazioni ripetute.
Simulazione, Correlogramma e Model Fitting
La procedura è simile a quella eseguita per simulare un processo ARMA(p,q), come descritto nella Parte 3 della serie ARMA. La differenza principale consiste nell’impostare \(d=1\), ovvero produrre una serie temporale non stazionaria con una componente stocastica di trend.
Dobbiamo adattare un modello ARIMA ai nostri dati simulati, provare a recuperare i parametri, creare intervalli di confidenza per questi parametri, produrre un correlogramma dei residui del modello adattato e infine eseguire un test di Ljung-Box per stabilire se abbiamo una buona vestibilità.
Simuliamo un modello ARIMA(1,1,1), con il coefficiente autoregressivo \(\alpha=0.6\) e il coefficiente di media mobile \(\beta=-0.5\). Di seguito il codice Python per simulare e tracciare questa tipologia di serie:
import numpy as np
import pandas as pd
from statsmodels.tsa.arima_process import ArmaProcess
np.random.seed(200)
ar_params = np.array([1, 0.6])
ma_params = np.array([1, -0.5])
returns = ArmaProcess(ar_params, ma_params).generate_sample(nsample=1000)
returns = pd.Series(returns)
drift = 100
price = pd.Series(np.cumsum(returns)) + drift
log_return = np.log(price) - np.log(price.shift(1))
log_return = log_return[1:]
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
plt.plot(log_return)
plot_acf(log_return, alpha=1, lags=30)
plt.show()
Ora che abbiamo la serie simulata, proviamo ad adattare un modello ARIMA(1,1,1). Dato che conosciamo l’ordine, possiamo semplicemente specificarlo nell’adattamento:
from statsmodels.tsa.arima_model import ARIMA
def fit_arima(log_returns):
ar_p = 1
ma_q = 1
diff_d = 0
# create tuple : (p, d, q)
order = (ar_p, diff_d, ma_q)
# create an ARIMA model object, passing in the values of the lret pandas series,
# and the tuple containing the (p,d,q) order arguments
arima_model = ARIMA(log_returns.values, order=order)
arima_result = arima_model.fit()
return arima_result
arima_result = fit_arima(log_return)
alpha = arima_result.arparams
beta = arima_result.maparams
# Print out summary information on the fit
print(arima_result.summary())
print(f"fitted AR parameter {alpha[0]:.2f}, MA parameter {beta[0]:.2f}")
ARMA Model Results
==============================================================================
Dep. Variable: y No. Observations: 999
Model: ARMA(1, 1) Log Likelihood 3222.558
Method: css-mle S.D. of innovations 0.010
Date: 10 Oct 2017 AIC -6437.117
Time: 17:13:24 BIC -6417.490
Sample: 0 HQIC -6429.657
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -9.744e-05 9.82e-05 -0.992 0.321 -0.000 9.5e-05
ar.L1.y -0.6311 0.029 -21.891 0.000 -0.688 -0.575
ma.L1.y -0.4736 0.032 -14.695 0.000 -0.537 -0.410
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 -1.5845 +0.0000j 1.5845 0.5000
MA.1 2.1113 +0.0000j 2.1113 0.0000
-----------------------------------------------------------------------------
fitted AR parameter -0.63, MA parameter -0.47
Entrambe le stime dei parametri rientrano negli intervalli di confidenza e sono vicine ai valori reali dei parametri della serie ARIMA simulata. Quindi, non dovremmo essere sorpresi di vedere i residui che sembrano una realizzazione di rumore bianco discreto:
resid = arima_result.resid
plot_acf(resid, alpha=1, lags=30)
plt.show()
Infine, possiamo eseguire un test Ljung-Box per fornire le prove statistiche di un buon adattamento:
from statsmodels.stats.diagnostic import acorr_ljungbox
# Ljung-Box test
ljung_box, p_value = acorr_ljungbox(resid, lags=[20])
print(f'Ljung-Box test: {ljung_box[-1]}')
print(f'p-value: {p_value[-1]}')
Ljung-Box test: 10.23667454613418
p-value: 0.9636744906155854
Possiamo vedere che il valore p è significativamente maggiore di 0,05 e come tale possiamo affermare che ci sono prove evidenti che il rumore bianco discreto si adatti bene ai residui. Quindi, il modello ARIMA(1,1,1) è adatto, come previsto.
Previsione delle Serie Finanziari
In questa sezione adatteremo i modelli ARIMA ad Amazon, Inc. (AMZN) e all’indice S&P500 US Equity (^GPSC, in Yahoo Finance). Usiamo la libreria StatsModels di Python.
Iniziamo usando yfinance per scaricare le serie dei prezzi giornalieri di Amazon dall’inizio del 2013. Inoltre calcoliamo i differenziali di primo ordine delle serie, in questo modo l’adattamento al modello ARIMA non richiede \(d > 0\) per il componente integrato:
import numpy as np
import yfinance as yf
import pandas as pd
symbol = 'AMZN'
data = yf.download(symbol, start="2013-01-01")
data = data.asfreq('b').fillna(method='ffill')
data['Return'] = data['Adj Close'].pct_change()
data['LogRet'] = np.log(data['Adj Close']).diff()
Come nella 3° parte della serie di articoli su modelli ARMA, analizziamo le combinazioni di \(p\), \(d\) e \(q\), per trovare il modello ARIMA(p,d,q) ottimale. Per “ottimale” si intende la combinazione di ordini che riduce al minimo l’Akaike Information Criterion (AIC):
from statsmodels.tsa.arima_model import ARIMA
order_aic = []
# Loop over p values from 0-4
for p in range(5):
# Loop over d values from 0-1
for d in range(2):
# Loop over q values from 0-4
for q in range(5):
try:
# create and fit ARMA(p,q) model
model = ARIMA(LogRet, order=(p, d, q))
results = model.fit()
# Append order and results tuple
order_aic.append((p, d, q, results.aic, results))
except Exception as e:
print(str(e))
pass
# Construct DataFrame from order_aic
order_df = pd.DataFrame(order_aic, columns=['p', 'd', 'q', 'AIC', 'Results'])
# Get the result of best order model
best_result = order_df.sort_values('AIC').iloc[0]
print("p = ", best_result.p)
print("d = ", best_result.d)
print("q = ", best_result.q)
Possiamo vedere che è stato selezionato un ordine di \(p=3\), \(d=0\), \(q=2\). In particolare si conferma \(d=0\) dato che abbiamo già considerato le differenze di primo ordine nella serie in ingresso:
p = 3
d = 0
q = 2
from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
# Get residuals of the model
resid = best_result.Results.resid
# Plot ACF of residuals
plot_acf(resid, alpha=0.05, lags=30)
plt.show()
Notiamo che ci sono due picchi significativi, a k=10 e k=23, nonostante dovremmo aspettarci di vedere picchi statisticamente significativi semplicemente dovuti alla variazione del 5% del periodo di campionamento.
Eseguiamo un test Ljung-Box (vedi articolo precedente) e verifichiamo se abbiamo ottenuto un buon adattamento:
from statsmodels.stats.diagnostic import acorr_ljungbox
# Ljung-Box test
ljung_box, p_value = acorr_ljungbox(resid, lags=[20])
print(f'Ljung-Box test: {ljung_box[-1]}')
print(f'p-value: {p_value[-1]}')
Ljung-Box test: 12.072540376919562
p-value: 0.9135565591469954
Come possiamo vedere, il valore p è maggiore di 0,05 e quindi abbiamo prove di un buon adattamento al livello del 95%.
Ora possiamo utilizzare il comando predict
per prevedere con alcuni giorni di anticipo la serie finanziaria di Amazon:
# Forecast
fc_len = int(len(data.index)*0.1)
fc, se, conf = results.forecast(fc_len, alpha=0.05) # 95% conf
# Make as pandas series
logret_series = pd.Series(LogRet.values)
start=logret_series.index[-1]
fc_series = pd.Series(fc, index=np.arange(start+1, start+len(fc)+1,1))
lower_series = pd.Series(conf[:, 0], index=fc_series.index)
upper_series = pd.Series(conf[:, 1], index=fc_series.index)
# Plot
plt.figure(figsize=(12,5), dpi=100)
plt.plot(logret_series, label='actual')
plt.plot(fc_series, label='forecast')
plt.fill_between(lower_series.index, lower_series, upper_series,
color='k', alpha=.15)
plt.title('Forecast vs Actuals')
plt.legend(loc='upper left', fontsize=8)
plt.show()
Notiamo che le previsioni puntuali per i prossimi giorni con bande di errore del 95%. Utilizzeremo queste previsioni nella nostra prima strategia di trading di serie temporali quando arriveremo a combinare ARIMA e GARCH.
Eseguiamo la stessa procedura per l’S&P500. Per prima cosa otteniamo i dati da yfinance e li convertiamo in un flusso di rendimenti logaritmici:
import numpy as np
import yfinance as yf
import pandas as pd
symbol = '^GSPC'
data = yf.download(symbol, start="2013-01-01")
data = data.asfreq('b').fillna(method='ffill')
Return = data['Adj Close'].pct_change()
LogRet = np.log(data['Adj Close']).diff().dropna()
Fittiamo un modello ARIMA eseguendo un ciclo sui valori di p, d e q:
from statsmodels.tsa.arima_model import ARIMA
order_aic = []
# Loop over p values from 0-4
for p in range(5):
# Loop over d values from 0-1
for d in range(2):
# Loop over q values from 0-4
for q in range(5):
try:
# create and fit ARMA(p,q) model
model = ARIMA(LogRet, order=(p, d, q))
results = model.fit()
# Append order and results tuple
order_aic.append((p, d, q, results.aic, results))
except Exception as e:
print(str(e))
pass
# Construct DataFrame from order_aic
order_df = pd.DataFrame(order_aic, columns=['p', 'd', 'q', 'AIC', 'Results'])
# Get the result of best order model
best_result = order_df.sort_values('AIC').iloc[0]
print("p = ", best_result.p)
print("d = ", best_result.d)
print("q = ", best_result.q)
L’AIC ci dice che il modello “migliore” è il modello ARIMA(4,0,3). Da notare ancora una volta che d=0, poiché abbiamo già considerato le differenze di primo ordine della serie:
p = 4
d = 0
q = 3
Possiamo tracciare i residui del modello adattato per verificare se otteniamo un rumore bianco discreto:
from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
# Get residuals of the model
resid = best_result.Results.resid
# Plot ACF of residuals
plot_acf(resid, alpha=0.05, lags=30)
plt.show()
Il correlogramma sembra promettente, quindi il passo successivo è eseguire il test Ljung-Box e confermare che abbiamo un buon adattamento del modello:
from statsmodels.stats.diagnostic import acorr_ljungbox
# Ljung-Box test
ljung_box, p_value = acorr_ljungbox(resid, lags=[20])
print(f'Ljung-Box test: {ljung_box[-1]}')
print(f'p-value: {p_value[-1]}')
Ljung-Box test: 13.44485001656763
p-value: 0.8574924325615751
Poiché il p-value è maggiore di 0,05, abbiamo verificato il buon adattamento del modello.
Perché nell’articolo precedente il nostro test Ljung-Box per l’S&P500 ha mostrato che l’ARMA(3,3) non era adatto ai rendimenti giornalieri logaritmici?
Si noti che abbiamo deliberatamente troncato i dati dell’S&P500 per iniziare dal 2013 in poi , che esclude convenientemente i periodi volatili intorno al 2007-2008. Quindi abbiamo escluso gran parte dell’S&P500 dove avevamo un eccessivo cluster di volatilità. Ciò influisce sulla correlazione seriale della serie e quindi ha l’effetto di far sembrare la serie “più stazionaria” rispetto al passato.
Questo è un punto molto importante . Quando si analizzano le serie temporali è necessario prestare la massima attenzione alle serie condizionatamente eteroschedastiche, come gli indici di borsa. Nella finanza quantitativa, il tentativo di determinare periodi di diversa volatilità è spesso noto come “rilevamento del regime”. È uno dei compiti più difficili da raggiungere!
Discuteremo questo punto a lungo nel prossimo articolo dove descriviamo i modelli ARCH e GARCH.
Tracciamo ora una previsione per i prossimi giorni dei rendimenti del registro giornaliero dell’S&P500:
Ora che abbiamo la capacità di adattare e prevedere modelli come ARIMA, siamo molto vicini alla possibilità di creare indicatori per una strategia di trading.
Prossimi passi
Nel prossimo articolo daremo un’occhiata al modello Generalized Autoregressive Conditional Heteroscedasticity (GARCH) e lo useremo per spiegare meglio la correlazione seriale in alcune serie di titoli azionari e serie di indici azionari.
Dopo aver descritto il GARCH, saremo in grado di combinarlo con il modello ARIMA e creare indicatori di segnale e quindi una strategia base di trading quantitativo.