Questo è il terzo e ultimo post della miniserie sui modelli Autoregressive Moving Average (ARMA) per l’analisi delle serie temporali. Abbiamo introdotto i modelli autoregressivi e i modelli a media mobile nei due articoli precedenti. Ora è il momento di combinarli per produrre un modello più sofisticato.
In definitiva, questo ci porterà ai modelli ARIMA e GARCH che ci permetteranno di prevedere i rendimenti degli asset e prevedere la volatilità. Questi modelli costituiranno la base per i segnali di trading e le tecniche di gestione del rischio.
Se hai letto la Parte 1 e la Parte 2 avrai visto che tendiamo a seguire uno schema per l’analisi di un modello di serie temporali. Lo ripeto brevemente qui:
- Motivazione – Perché ci interessa questo particolare modello?
- Definizione – Una definizione matematica per ridurre l’ambiguità.
- Correlogramma – Tracciare un correlogramma di esempio per visualizzare il comportamento di un modello.
- Simulazione e adattamento – Adattamento del modello alle simulazioni, al fine di garantire la corretta comprensione del modello.
- Dati finanziari reali : applica il modello ai prezzi storici reali degli asset.
- Predizione – Prevedi i valori successivi per creare segnali o filtri di trading.
Per seguire questo articolo si consiglia di dare un’occhiata ai precedenti articoli sull’analisi delle serie temporali. Si possono trovare tutti qui .
Criterio di informazione bayesiano
Nella parte 1 di questa serie di articoli abbiamo descritto il Criterio di Informazione di Akaike (AIC) come un mezzo per aiutarci a scegliere tra i “migliori” modelli di serie temporali.
Uno strumento strettamente correlato è il Bayesian Information Criterion (BIC). In sostanza ha un comportamento simile all’AIC in quanto penalizza i modelli che hanno troppi parametri. Ciò può portare a un sovra-adattamento. La differenza tra BIC e AIC è che il BIC è più rigoroso con la sua penalizzazione dei parametri aggiuntivi.
Criterio di informazione bayesiano
Se prendiamo la funzione di verosimiglianza per un modello statistico, che ha K parametri e L massimizza la verosimiglianza , allora il criterio di informazione bayesiano è dato da:\(\begin{eqnarray}BIC = -2 \text{log} (L) + k \text{log}(n)\end{eqnarray}\)
In cui \(n\) è il numero di punti dati nella serie temporale.
Usiamo l’AIC e il BIC quando dobbiamo scegliere gli appropriati modelli ARMA(p,q).
Test Ljung-Box
Il test di Ljung-Box è un classico test di ipotesi progettato per verificare se un insieme di autocorrelazioni di un adattato modello di serie temporale differisce significativamente da zero. Il test non verifica la casualità di ogni singolo ritardo, ma piuttosto verifica la casualità su un gruppo di ritardi.
Formalmente:
Test Ljung-Box
Definiamo l’ipotesi nulla \({\bf H_0}\) come: I dati delle serie temporali ad ogni ritardo sono iid , ovvero le correlazioni tra i valori delle serie della popolazione sono zero.Definiamo l’ipotesi alternativa \({\bf H_a}\) come: I dati delle serie temporali non sono iid e possiedono una correlazione seriale.
Calcoliamo il seguente test statistico, \(Q\):
\(\begin{eqnarray}Q = n (n+2) \sum_{k=1}^{h} \frac{\hat{\rho}^2_k}{n-k}\end{eqnarray}\)
In cui si n è la lunghezza del campione della serie storica, \(\hat{\rho}_k\) è l’autocorrelazione del campione per il ritardo k e h è il numero di ritardi nel test.
La regola per decidere se rifiutare l’ipotesi nulla \({\bf H_0}\)
consiste nel verificare se \(Q > \chi^2_{\alpha,h}\), per una distribuzione del chi quadrato con h gradi di libertà al \(100(1-\alpha)\)-esimo percentile.
Sebbene i calcoli del test possano sembrare leggermente complessi, possiamo utilizzare Python per calcolare il test, semplificando un po’ la procedura.
Modelli di media mobile autogressiva (ARMA) di ordine p, q
Ora che abbiamo discusso del BIC e del test Ljung-Box, siamo pronti per discutere il nostro primo modello misto, ovvero la media mobile autoregressiva di ordine p, q o ARMA(p,q).
Fondamento logico
Ad oggi abbiamo considerato i processi autoregressivi e i processi a media mobile.
Il primo modello considera il proprio comportamento passato come input per il modello e come tale tenta di catturare gli effetti dei partecipanti al mercato, come lo slancio e l’inversione alla media nel trading di azioni.
Il secondo modello viene utilizzato per caratterizzare le informazioni “shock” in una serie, come un annuncio di guadagni a sorpresa o un evento inaspettato (come la fuoriuscita di petrolio della BP Deepwater Horizon ).
Pertanto, un modello ARMA tenta di catturare entrambi questi aspetti durante la modellazione di serie temporali finanziarie.
Si noti che un modello ARMA non tiene conto del clustering della volatilità, un fenomeno empirico chiave di molte serie temporali finanziarie. Non è un modello condizionalmente eteroscedastico. Per questo dovremo aspettare i modelli ARCH e GARCH.
Definizione
Il modello ARMA(p,q) è una combinazione lineare di due modelli lineari e quindi è esso stesso lineare:
Modello autoregressivo di media mobile di ordine p, q
\(\begin{eqnarray}x_t = \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + \ldots + w_t + \beta_1 w_{t-1} + \beta_2 w_{t-2} \ldots + \beta_q w_{t-q}\end{eqnarray}\)
Un modello di serie storica, \(\{ x_t \}\), è un modello a media mobile autoregressiva di ordine p, q, ARMA(p,q), se:Dove \(\{ w_t \}\) è il rumore bianco con \(E(w_t) = 0\) e varianza \(\sigma^2\).
Se consideriamo l’ operatore di spostamento all’indietro B,(vedi articolo precedente) allora possiamo riscrivere quanto sopra come una funzione \(\theta\) e \(\phi\) di B
\(\begin{eqnarray}\theta_p ({\bf B}) x_t = \phi_q ({\bf B}) w_t\end{eqnarray}\)
Vediamo direttamente che impostando \(p \neq 0\) e \(q=0\) ricaviamo il modello AR(p). Allo stesso modo se impostiamo \(p=0\) e \(q \neq 0\) ricaviamo il modello MA(q).
Le caratteristiche chiave del modello ARMA sono moderazione e ridondanza nei suoi parametri. Cioè, un modello ARMA richiederà spesso meno parametri rispetto al solo modello AR(p) o solo modello MA(q). Inoltre, se riscriviamo l’equazione in termini di BSO, allora polinomi \(\theta\) e \(\phi\) possono a volte condividere un fattore comune, portando così a un modello più semplice.
Simulazioni e Correlogrammi
Come per i modelli autoregressivi e a media mobile, ora simuliamo varie serie ARMA e tentiamo di adattare i modelli ARMA a queste realizzazioni. In questo modo possiamo discutere e comprendere la procedura di adattamento, incluso il calcolo degli intervalli di confidenza per i modelli, oltre a garantire che la procedura recuperi effettivamente stime ragionevoli per i parametri ARMA originali.
Nella parte 1 e nella parte 2 abbiamo costruito manualmente le serie AR e MA disegnando \(N\) campioni da una distribuzione normale e quindi abbiamo creato lo specifico modello di serie temporali utilizzando i ritardi di questi campioni.
ARMA(1,1)
Iniziamo con il modello ARMA non banale più semplice possibile, ovvero il modello ARMA(1,1). Questo è un modello autoregressivo di ordine uno combinato con un modello di media mobile di ordine uno. Tale modello ha solo due coefficienti \(\alpha\),e \(\beta\), che rappresentano i primi ritardi della stessa serie storica e i termini “shock” del rumore bianco. Tale modello è dato da:
\(\begin{eqnarray}x_t = \alpha x_{t-1} + w_t + \beta w_{t-1}\end{eqnarray}\)
Dobbiamo specificare i coefficienti prima della simulazione. Considerazioni \(\alpha = 0.5\) e \(\beta = 0.5 \):
import matplotlib.pyplot as plt
import numpy as np
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.graphics.tsaplots import plot_acf
ar1 = np.array([1, 0.5])
ma1 = np.array([1, 0.5])
simulated_ARMA_data = ArmaProcess(ar1, ma1).generate_sample(nsample=1000)
plt.plot(simulated_ARMA_data)
plot_acf(simulated_ARMA_data)
plt.show()
ARMA della libreria Statsmodels:
from statsmodels.tsa.arima_model import ARMA
# Fit an ARMA(1,1) model to the first simulated data
mod = ARMA(simulated_ARMA_data, order=(1, 1))
res = mod.fit()
# Print out summary information on the fit
print(res.summary())
# Print out the estimate for the constant and for alpha and beta
print("When the estimate values of alpha and beta (and the constant) are:")
print(res.params)
Il report prodotto dalla funzione ARMA contiene molte informazioni utili. In primo luogo, possiamo vedere come i parametri stimati \(\hat{\alpha_1} = 0.504\) e \(\hat{\beta_1} = 0.516\), che sono molto vicini ai valore reali di \(\hat{\alpha_1} = 0.5\) e \(\hat{\beta_1} = 0.5\). In secondo luogo, gli errori standard sono già calcolati per noi, rendendo semplice il calcolo degli intervalli di confidenza.
Possiamo vedere che l’intervallo di confidenza al 95% contiene i valori dei parametri veri e quindi possiamo giudicare il modello adatto.
ARMA Model Results
==============================================================================
Dep. Variable: y No. Observations: 1000
Model: ARMA(1, 1) Log Likelihood -1436.502
Method: css-mle S.D. of innovations 1.017
Date: 04 Oct 2017 AIC 2881.004
Time: 16:14:55 BIC 2900.635
Sample: 0 HQIC 2888.465
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -0.0940 0.098 -0.957 0.339 -0.287 0.099
ar.L1.y 0.5043 0.033 15.125 0.000 0.439 0.570
ma.L1.y 0.5164 0.032 16.083 0.000 0.453 0.579
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.9830 +0.0000j 1.9830 0.0000
MA.1 -1.9366 +0.0000j 1.9366 0.5000
-----------------------------------------------------------------------------
When the estimate values of alpha and beta (and the constant) are:
[-0.09400128 0.50428429 0.51637292]
ARMA(2,2)
Proviamo ora un modello ARMA(2,2). Ovvero, un modello AR(2) combinato con un modello MA(2). Dobbiamo specificare quattro parametri per questo modello \(\alpha_1\), \(\alpha_2\), \(\beta_1\), \(\beta_2\).
Consideriamo \(\alpha_1 = 0.5\), \(\alpha_1 = -0.25\), \(\beta_1 = 0.5\) e \(\beta_2 = -0.3\):
import matplotlib.pyplot as plt
import numpy as np
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.graphics.tsaplots import plot_acf
ar1 = np.array([1, -0.5, 0.25])
ma1 = np.array([1, 0.5, -0.3])
simulated_ARMA_data = ArmaProcess(ar1, ma1).generate_sample(nsample=1000)
plt.plot(simulated_ARMA_data)
plot_acf(simulated_ARMA_data, alpha=1, lags=30)
plt.show()
from statsmodels.tsa.arima_model import ARMA
# Fit an ARMA(2,2) model to the first simulated data
mod = ARMA(simulated_ARMA_data, order=(2, 2, 0))
res = mod.fit()
# Print out summary information on the fit
print(res.summary())
# Print out the estimate for the constant and for alpha and beta
print("When the estimate values of alpha and beta (and the constant) are:")
print(res.params)
ARMA Model Results
==============================================================================
Dep. Variable: y No. Observations: 1000
Model: ARMA(2, 2) Log Likelihood -1393.495
Method: css-mle S.D. of innovations 0.974
Date: 04 Oct 2017 AIC 2798.990
Time: 18:31:36 BIC 2828.436
Sample: 0 HQIC 2810.182
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -0.0026 0.048 -0.053 0.958 -0.097 0.092
ar.L1.y 0.4874 0.124 3.920 0.000 0.244 0.731
ar.L2.y -0.2470 0.034 -7.240 0.000 -0.314 -0.180
ma.L1.y 0.4947 0.127 3.896 0.000 0.246 0.744
ma.L2.y -0.3050 0.117 -2.596 0.010 -0.535 -0.075
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 0.9867 -1.7537j 2.0122 -0.1684
AR.2 0.9867 +1.7537j 2.0122 0.1684
MA.1 -1.1731 +0.0000j 1.1731 0.5000
MA.2 2.7953 +0.0000j 2.7953 0.0000
-----------------------------------------------------------------------------
When the estimate values of alpha and beta (and the constant) are:
[-0.00255363 0.4873853 -0.24698618 0.49470116 -0.30496004]
Esiste la possibilità reale che gli intervalli di confidenza per i coefficienti delle due componenti possono effettivamente non contenere il valore dei parametri originale. Questo delinea il pericolo di tentare di adattare i modelli ai dati, anche quando conosciamo i veri valori dei parametri!
Tuttavia, ai fini del trading, abbiamo solo bisogno di avere un potere predittivo che superi la casualità (lancio di una moneta) e produca un profitto sufficiente al netto dei costi di transazione, per essere redditizio nel lungo periodo.
Ora che abbiamo visto alcuni esempi di modelli ARMA simulati, abbiamo bisogno di definire un metodo che ci permetta di scegliere i valori di p e q quando si adattano i modelli di dati finanziari reali.
Scelta del miglior modello ARMA(p,q).
import numpy as np
from statsmodels.tsa.arima_process import ArmaProcess
ar1 = np.array([1, -0.5, 0.25, -0.4])
ma1 = np.array([1, 0.5, -0.3])
simulated_ARMA_data = ArmaProcess(ar1, ma1).generate_sample(nsample=1000)
from statsmodels.tsa.arima_model import ARMA
import pandas as pd
# Create empty list to store search results
order_aic = []
# Loop over p values from 0-4
for p in range(5):
# Loop over q values from 0-4
for q in range(5):
try:
# create and fit ARMA(p,q) model
model = ARMA(simulated_ARMA_data, order=(p, q))
results = model.fit()
# Append order and results tuple
order_aic.append((p, q, results.aic, results))
except:
pass
# Construct DataFrame from order_aic
order_df = pd.DataFrame(order_aic_bic, columns=['p', 'q', 'AIC', 'Results'])
Dopo aver eseguito il fitting dei modelli su diversi ordini p e q, possiamo valutare i risultati e ricavare l’ordine del modello migliore in base all’AIC minore.
from statsmodels.tsa.arima_model import ARMA
import pandas as pd
# Create empty list to store search results
order_aic = []
# Loop over p values from 0-4
for p in range(5):
# Loop over q values from 0-4
for q in range(5):
try:
# create and fit ARMA(p,q) model
model = ARMA(simulated_ARMA_data, order=(p, q))
results = model.fit()
# Append order and results tuple
order_aic.append((p, q, results.aic, results))
except:
pass
# Construct DataFrame from order_aic
order_df = pd.DataFrame(order_aic_bic, columns=['p', 'q', 'AIC', 'Results'])
Possiamo vedere che è l’ordine originale del modello ARMA simulato è stato selezionato come ordine migliore, cioè \(p=3\) e \(q=2\).
# Print order_df in order of increasing AIC
print(order_df.sort_values('AIC'))
# Get the result of best order model
best_result = order_df.sort_values('AIC').iloc[0]
print("p = ", best_result.p)
print("q = ", best_result.q)
res = best_result.Results
ARMA dove si evidenzia che parametri stimati sono molto vicini ai valore reali. Inoltre l’intervallo di confidenza al 95% contiene i valori dei parametri veri e quindi possiamo giudicare il modello adatto.
ARMA Model Results
==============================================================================
Dep. Variable: y No. Observations: 1000
Model: ARMA(3, 2) Log Likelihood -1437.719
Method: css-mle S.D. of innovations 1.018
Date: 04 Oct 2017 AIC 2889.438
Time: 19:34:48 BIC 2923.793
Sample: 0 HQIC 2902.495
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -0.0207 0.094 -0.220 0.826 -0.206 0.164
ar.L1.y 0.5027 0.087 5.750 0.000 0.331 0.674
ar.L2.y -0.2483 0.037 -6.635 0.000 -0.322 -0.175
ar.L3.y 0.3333 0.032 10.427 0.000 0.271 0.396
ma.L1.y 0.5143 0.091 5.642 0.000 0.336 0.693
ma.L2.y -0.3023 0.086 -3.506 0.000 -0.471 -0.133
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.3215 -0.0000j 1.3215 -0.0000
AR.2 -0.2883 -1.4789j 1.5067 -0.2806
AR.3 -0.2883 +1.4789j 1.5067 0.2806
MA.1 -1.1572 +0.0000j 1.1572 0.5000
MA.2 2.8587 +0.0000j 2.8587 0.0000
-----------------------------------------------------------------------------
The estimate of alpha and beta parameter (and the constant) are:
[-0.02072297 0.50269967 -0.24828968 0.33331434 0.51434838 -0.30228735]
from statsmodels.graphics.tsaplots import plot_acf
# Get residuals of the model
resid = res.resid
# Plot ACF of residuals
plot_acf(resid, alpha=1, lags=30)
# 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.287236878999057
p-value: 0.9626594519387415
Analisi Dati Finanziari dell'SP500
Ora che abbiamo delineato la procedura per scegliere il modello di serie temporale ottimale per una serie simulata, è piuttosto semplice applicarlo ai dati finanziari. Per questo esempio scegliamo l’indice S&P500 dell’azionariato statunitenze.
Scarichiamo i prezzi di chiusura giornalieri utilizzando yfinance e quindi creiamo il flusso dei ritorni logaritmici:
import numpy as np
import pandas as pd
import yfinance as yf
data = yf.download("^GSPC", period="max")
sp500 = np.log(data['Adj Close'] / data['Adj Close'].shift(1)).dropna()
from statsmodels.tsa.arima_model import ARMA
# Create empty list to store search results
order_aic = []
# Loop over p values from 0-4
for p in range(5):
# Loop over q values from 0-4
for q in range(5):
try:
# create and fit ARMA(p,q) model
model = ARMA(sp500, order=(p, q))
results = model.fit()
# Append order and results tuple
order_aic.append((p, q, results.aic, results))
except:
pass
# Construct DataFrame from order_aic
order_df = pd.DataFrame(order_aic, columns=['p', '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("q = ", best_result.q)
p = 4
q = 3
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=1, lags=30);
plt.show()
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: 40.390954304955045
p-value: 0.0044565024322041746
Prossimi Passi
In questa serie di articoli abbiamo visto prove di eteroschedasticità condizionale (cluster di volatilità) nella serie S&P500, specialmente nei periodi intorno al 2007-2008. Quando introdurremo un modello GARCH, più avanti nella serie di articoli, vedremo come eliminare queste autocorrelazioni.
In pratica, i modelli ARMA non sono mai generalmente adatti per i rendimenti logaritmici delle azioni. Dobbiamo prendere in considerazione l’eteroscedasticità condizionale e utilizzare una combinazione di ARIMA e GARCH. Il prossimo articolo descriverà il modello ARIMA e mostrerà come il componente “Integrato” differisca dal modello ARMA che abbiamo considerato in questo articolo.
