Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

É possível trabalhar com análise de séries temporais usando Python, mesmo sem dominar muito computação?

0 votos
1,512 visitas
perguntada Fev 27, 2016 em Estatística por danielcajueiro (5,251 pontos)  
Compartilhe

1 Resposta

0 votos
respondida Fev 27, 2016 por danielcajueiro (5,251 pontos)  

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:

spread

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:

python autocorrelation function

python partial autocorrelation funcion

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
...