Sim, sem dúvida. Veja o exemplo abaixo. Trabalhar com séries temporais em Python é tão simples (ou até mais simples) quanto qualquer outro programa como, por exemplo, R, Eviews, SAS...
Vamos exemplificar usando os dados do arquivo quarterly.xls, que é material de apoio do livro do Walter Enders - Applied Econometric Time Series, que pode ser baixado do site do livro.
Esse arquivo apresenta dados trimestrais da economia Americana do período 1960 a 2012. Usando essa base vamos estudar a variável \(spread\) que é a diferença entre a taxa de juros de títulos do governo americano de longo prazo (5 anos) e da taxa de juros de curto prazo (3 meses) que é usada como um indicador de atividade econômica.
Em qualquer típica estudo de série temporais, o nosso primeiro desafio é carregar a série e plotar ela em um gráfico para verificar o comportamento geral dela. O código abaixa implementa essa tarefa usando Pandas (um pacote do Python)
df = pd.read_excel('quarterly.7775706.xls',sheetname='TB3MS')
# Creation of the variable spread
df['spread']=df['r5']-df['Tbill']
df.plot(x='DATE',y='spread')
A seguinte figura é gerada:
Um procedimento bastante comum é calcular as funções de auto-correlação e as funções de auto-correlação parcial para identificar os modelos que mais se adequam aos dados :
graf.tsaplots.plot_acf(df.spread, lags=12, alpha=0.05)
graf.tsaplots.plot_pacf(df.spread, lags=12, alpha=0.05)
Essas linhas de código geram as seguintes figuras:
Obviamente, em qualquer análise considerando dados reais, existe ambiguidade na avaliação dessas figuras. Normalmente, nessa fase é comum estimar vários modelos e depois efetuar alguns diagnósticos.
Tendo como base essas figuras, parece razoável concluir que:
1) O processo acima não não é um \(MA(q)\) puro, pois não aparenta ter um lag em que depois dele a função de auto-correlação se anula.
2) Se for um \(AR(p)\) puro, ele provavelmente deve ter uma ordem razoavelmente alta [\(AR(6)\) por exemplo] para considerar os termos significativos da função de auto-correlação parcial.
Para exemplificar, vamos estimar os processos \(ARMA(1,1)\) e \(ARMA(2,1)\). O código em Python usado para fazer isso é
dates=sm.tsa.datetools.dates_from_range('1960Q1', length=len(df['spread']))
y = pd.Series(df['spread'].as_matrix(), index=dates,)
# ARMA(1,1)
arma11=sm.tsa.ARMA(y, (1,1),freq='Q').fit()
print arma11.summary
# ARMA(2,1)
arma21=sm.tsa.ARMA(y, (2,1),freq='Q').fit()
print arma21.summary
Esse código gera os seguintes sumários:
ARMA(1,1)
ARMA Model Results
==============================================================================
Dep. Variable: y No. Observations: 212
Model: ARMA(1, 1) Log Likelihood -142.748
Method: css-mle S.D. of innovations 0.473
Date: Sex, 26 Fev 2016 AIC 293.497
Time: 22:52:07 BIC 306.923
Sample: 03-31-1960 HQIC 298.924
- 12-31-2012
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const 1.1932 0.181 6.580 0.000 0.838 1.549
ar.L1.y 0.7572 0.054 13.973 0.000 0.651 0.863
ma.L1.y 0.3771 0.095 3.961 0.000 0.191 0.564
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.3206 +0.0000j 1.3206 0.0000
MA.1 -2.6520 +0.0000j 2.6520 0.5000
-----------------------------------------------------------------------------
ARMA(2,1)
ARMA Model Results
==============================================================================
Dep. Variable: y No. Observations: 212
Model: ARMA(2, 1) Log Likelihood -140.101
Method: css-mle S.D. of innovations 0.467
Date: Sex, 26 Fev 2016 AIC 290.201
Time: 22:53:23 BIC 306.984
Sample: 03-31-1960 HQIC 296.985
- 12-31-2012
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const 1.1896 0.199 5.973 0.000 0.799 1.580
ar.L1.y 0.4131 0.122 3.399 0.001 0.175 0.651
ar.L2.y 0.3188 0.114 2.791 0.006 0.095 0.543
ma.L1.y 0.6979 0.091 7.666 0.000 0.519 0.876
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.2379 +0.0000j 1.2379 0.0000
AR.2 -2.5336 +0.0000j 2.5336 0.5000
MA.1 -1.4328 +0.0000j 1.4328 0.5000
-----------------------------------------------------------------------------
Um critério bem comum para comparar a qualidade dos modelos é comparar os AIC e BIC dos modelos apresentados no sumário acima. Entretanto, enquanto o AIC sugere que o \(ARMA(2,1)\) é melhor, o BIC sugere o contrário. Será que estamos enfrentando uma dúvida cruel?
Não! Existem outros testes que podem nos ajudar. Um bastante importante que nos ajuda a explorar se a estrutura do modelo é adequada é o teste de Ljung-Box que pode ser feito em Python escrevendo essas linhas de código:
[LBvalue,LBpvalue]=sm.stats.diagnostic.acorr_ljungbox(arma11.resid,lags=8)
print LBpvalue
[LBvalue,LBpvalue]=sm.stats.diagnostic.acorr_ljungbox(arma21.resid,lags=8)
print LBpvalue
Eles geram respectivamente os seguintes resultados:
ARMA(1,1): Lags 1 a 8 (pvalores)
[ 0.54311556 0.72990516 0.09439483 0.16795887 0.14216097 0.12066667
0.05890125 0.01988974]
ARMA(2,1) Lags 1 a 8 (pvalores)
[ 0.94877215 0.960222 0.79473837 0.87935057 0.92438263 0.79336478
0.38043862 0.14601597]
Lembre que a hipótese nula do teste de Ljung-Box é que os dados são independentes e claramente essa hipótese nula é rejeitada para o Modelo ARMA(1,1) no caso dos lags mais altos. Logo, o modelo mais adequado considerando esses dois modelos é o ARMA(2,1).
Existem várias outras questões que poderiam ser avaliadas. Por exemplo, previsão fora da amostra...
Se você tem interesse em lidar com problemas de estatística ou econometria usando Python olhe por exemplo:
Python for econometrics - Wooldridge
Como rodar modelos de regressão em Python sem saber computação?
Se você gostaria de melhorar suas habilidade de computação, dê uma olhada aqui:
Como um programador iniciante pode se tornar um programador intermediário?
Como um programador em nível intermediário pode se tornar um programador em nível avançado?
O código completo implementado para esse exercício foi
import numpy as np
import statsmodels.api as sm
import statsmodels.graphics as graf
import pandas as pd
import matplotlib.pyplot as plt
# Data from Walter Enders available at http://time-series.net/data_sets
if __name__ == '__main__':
df = pd.read_excel('quarterly.7775706.xls',sheetname='TB3MS')
# Creation of the variable spread
df['spread']=df['r5']-df['Tbill']
df.plot(x='DATE',y='spread')
#plt.figure()
df['spreadFirstDif']=df['spread'].diff()
df.plot(x='DATE',y='spreadFirstDif')
# Plotting acf and pacf
graf.tsaplots.plot_acf(df.spread, lags=12, alpha=0.05)
graf.tsaplots.plot_pacf(df.spread, lags=12, alpha=0.05)
# Arma(2,1) estimation
dates=sm.tsa.datetools.dates_from_range('1960Q1', length=len(df['spread']))
y = pd.Series(df['spread'].as_matrix(), index=dates,)
arma11=sm.tsa.ARMA(y, (1,1),freq='Q').fit()
print arma11 .summary()
[LBvalue,LBpvalue]=sm.stats.diagnostic.acorr_ljungbox(arma11.resid,lags=8)
print LBpvalue
arma21=sm.tsa.ARMA(y, (2,1),freq='Q').fit()
print arma21.summary()
[LBvalue,LBpvalue]=sm.stats.diagnostic.acorr_ljungbox(arma21.resid,lags=8)
print LBpvalue