proposta de função executável em r para ajuste de modelos arima(p, d, q)

13 Pages·2010·0.31 MB·Portuguese

XXX ENCONTRO NACIONAL DE ENGENHARIA DE PRODUÇÃO
Maturidade e desafios da Engenharia de Produção: competitividade das empresas, condições de trabalho, meio ambiente.

São Carlos, SP, Brasil, 12 a15 de outubro de 2010.

PROPOSTA DE FUNÇÃO EXECUTÁVEL
EM R PARA AJUSTE DE MODELOS
ARIMA(P, D, Q)

Lúcia Adriana dos Santos Gruginskie (UNISINOS)
[email protected]

A proposta deste artigo é construir e testar uma função, executável no
R, que ajuste (p+1)*(q+1)*(d+1) modelos ARIMA(p, d, q) a uma série
temporal univariada, retornando uma matriz cujas colunas apresentem
diversas medidas (como AIC, BIC e mape) sobre os modelos
ajustados. A função é executada sem a necessidade de instalação do
pacote forecast.

Palavras-chaves: Modelos ARIMA, R, previsão
1 Introdução
Encontra-se disponível, hoje, uma gama de artigos e dissertações sobre previsão de demanda
na engenharia de produção.
Esses estudos utilizam-se de dados históricos - ou dados históricos combinados com outros
métodos de previsão, como a opinião de especialistas - para prever a demanda de produtos.
Em relação à análise de séries históricas para previsão, os modelos de Box & Jenkins, também
conhecidos por modelos ARIMA (Autorregressivo Integrado de Média Móvel), são bastante
utilizados em engenharia de produção, de forma individual ou combinada, como, por
exemplo, nos trabalhos de Mileski (2007), Werner e Ribeiro (2006), Werner e Ribeiro (2003),
Dias (2006) e Pellegrini (2000), entre outros.
Na utilização da abordagem de Box & Jenkins, no curso do ciclo iterativo no qual os modelos
são construídos, o passo de identificação pode ser considerado o mais difícil/crítico.
(MORETTIN & TOLOI, 2004).
Pacotes estatísticos comerciais, como o SPSS, sugerem modelos de ajuste, facilitando o
processo de identificação. Softwares livres, como o R, também possuem formas de identificar
o melhor modelo para o ajuste de séries temporais.
A proposta deste artigo é construir e testar uma função, a ser executada no R, que ajuste
vários modelos ARIMA (p, d, q) a uma série e retorne uma matriz contendo as ordens dos
modelos ajustados e uma série de medidas calculadas com o objetivo de facilitar o processo
de identificação do modelo final a ser ajustado.
A matriz resultante não dispensa o conhecimento do analista, mas reduz o tempo de
processamento à medida que ajusta vários modelos de forma consecutiva através da
combinação das ordens “p”, “q” e “d”.
O artigo está estruturado em cinco seções. Após a introdução, segue uma seção com conceitos
básicos sobre a abordagem ARIMA, concentrando-se na fase de identificação, enquanto a
terceira apresenta algumas considerações sobre o R – A language and environment for
statistical computing. A quarta seção testa a função em três séries obtidas na literatura e, por
fim, a quinta seção traz algumas considerações finais.
2 Os modelos de Box e Jenkins
Os modelos foram introduzidos por Box e Jenkins em 1970 e muitos autores ora usam o
termo ARIMA (ou AR, ARI, ARMA, IMA, dependendo das ordens dos modelos), ora
modelos de Box & Jenkins.
A suposição destes modelos é que a variação contida na série pode ser representada em três
componentes: o autorregressivo (AR), o componente integrado (I) ou diferenciado e o
componente média móvel (MA). (SPSS INC, 2006).
A forma geral dos modelos é ARIMA(p, d, q)(P, D, Q), onde:
 “p” refere-se ao processo autorregressivo não sazonal incorporado no modelo ARIMA e
“P” refere-se ao processo autorregressivo sazonal. A ordem da autorregressão (p) refere-se
a diferença de tempo entre o valor atual e os valores que o predizem;
 “d” refere-se à ordem de integração ou diferenciação não sazonal e “D” à ordem de
integração sazonal;
2
 “q” refere-se ao processo de média móvel não sazonal incorporada ao modelo e “Q”
refere-se ao processo de média móvel sazonal.
Este artigo não tratará de modelos sazonais. O termo não estacionariedade de uma série
implica existir inclinação nos dados e/ou flutuações aumentando ou diminuindo com o passar
do tempo (este último pode indicar variância alterada). Para Morettin e Toloi (2004), há duas
razões para transformar a série: estabilizar a variância e tornar o efeito sazonal aditivo.
Uma dada série Z , não estacionária, pode tornar-se estacionária tomando um número finito
t
de diferenças:
dZ
t t
A série , estacionária não sazonal,pode ser escrita como:
t
  ...  a a a ...a
t 1 t1 2 t2 p tp t 1 t1 2 t2 q tq
Onde:
é o valor da série no instante t após a integração (quando for o caso);
t
 são os parâmetros autorregressivos, i = i, 2, ..., p;
i
são os parâmetros das médias móveis, i = 1, 2, ... q;
1
a representa o ruído branco (ruído com média 0 e variância 2).
t a
A estratégia para a construção do modelo está baseada em um ciclo iterativo baseado nos
próprios dados, conforme fluxograma abaixo:
Escolher
uma classe geral de modelo
Identificar o modelo
Estimar os parâmetros
Diagnosticar
(O modelo éadequado?)
Não
Sim
Usar o modelo para previsão
ou controle

Figura 1 – Estágios na abordagem iterativa para construção do modelo
Fonte: Adaptado de Box & Jenkins (1970)
3
Segundo Box e Jenkins (1970), a partir da interação entre teoria e prática, uma classe de
modelos para os fins em questão é considerada. Como essa classe é extensa para ser ajustada
diretamente aos dados de forma conveniente, métodos para a identificação de subclasses são
desenvolvidos. Tais métodos de identificação do modelo empregam o conhecimento do
sistema para sugerir uma subclasse apropriada, parcimoniosa, que pode ser provisoriamente
cogitada.
Usualmente, a forma de identificar o modelo é examinar as ACFs e PACFs, respectivamente,
as funções de autocorrelação e autocorrelação parcial. A ACF e a PACF são ferramentas para
verificar o desempenho e a especificação dos modelos de séries temporais. A autocorrelação é
usada para descrever a correlação entre dois valores da mesma série temporal, em diferentes
períodos de tempo. A autocorrelação parcial é utilizada para descrever a correlação entre dois
valores da mesma série, porém retirando a influência de valores intermediários.
Para auxiliar na identificação, Shumway (1988) relaciona os seguintes comportamentos das
funções de autocorrelação e autocorrelação parcial:
 Não estacionariedade: a ACF decai vagarosamente e a PACF tem um valor significativo
(positivo ou negativo) no lag 1;
 Não estacionariedade sazonal: a ACF é próxima de zero, exceto nos lags s, 2s e assim por
diante e também decai vagarosamente. A série também pode se tornar estacionária através
da diferenciação sazonal, tomando, usualmente, uma única diferença;
 Comportamento autorregressivo: a PACF é diferente de zero para lags 1, 2,..., p e zero
para os demais;
 Comportamento autorregressivo sazonal: A PACF é zero, exceto para os lags s, 2s, ... e
zero para os demais;
 Comportamento média móvel: a ACF é diferente de zero para lags 1, 2, ... q e zero para os
demais;
 Comportamento média móvel sazonal: a ACF é zero exceto para os lags s, 2s...
Sobre as ordens dos modelos, Box et al. (apud Pellegrini, 2000), indicam que “p” e “q” são
geralmente menores que 2 para séries temporais estacionárias e Shumway (1988) sugere que
pode-se tornar uma série estacionária através da diferenciação, quando, usualmente, d<=2,
aconselhando a tomar cuidado para evitar a superdiferenciação na série.
O próximo passo é a estimação dos parâmetros. É importante, na prática, que se empregue o
menor número possível de parâmetros para uma adequada representação. A regra central é o
princípio da parcimônia (BOX E JENKINS, 1970).
O ciclo iterativo segue com a fase de verificação, quando pode ser feita a análise dos resíduos
e a realização de testes como o de Box Ljung.
O teste de Box Ljung objetiva verificar se as autocorrelações dos resíduos apresentam valores
“altos”, embora detecte quebras específicas no comportamento de ruído branco (MORETTIN
e TOLOI, 2004).
De acordo com Morettin e Toloi (2004), a comparação dos valores da autocorrelação
estimada dos resíduos com os limites de ±2n fornece uma indicação de possível quebra de
comportamento de ruído branco, com a condição de que, para pequenos valores de k (relações
entre lags próximos), esses limites subestimarão a significância de qualquer discrepância.
O artigo não visa a aprofundar conceitos nos modelos de Box e Jenkins, pois podem ser
4
obtidos em Box e Jenkins (1970), Shumway (1988) e Morettin e Toloi (2004), entre outros.
Mesmo assim, alguns termos utilizados adiante serão conceituados como o teste de Dickey-
Fuller e medidas de bondade de ajuste.
O teste de Dickey-Fuller testa a hipótese nula que a série tem uma raiz unitária contra a
alternativa de estacionariedade e, segundo Werner e Ribeiro (2003), é o teste para raízes
unitárias mais utilizado.
Em relação aos coeficientes estimados para p, q e intercepto, considerando que o tamanho da
amostra é grande, de acordo com Morettin e Toloi (2004), os estimadores de máxima
verossimilhança possuem uma distribuição assintótica normal e, a partir das estimativas das
variâncias, pode-se obter intervalos de confiança para os coeficientes.
Medidas de bondade de ajuste para modelos de séries temporais são usualmente associadas
com o erro das previsões.
Porém, para comparar modelos com diferentes números de parâmetros, os critérios AIC
(Akaike information criteria) e BIC (Bayesian information criterion) podem usados, pois
modelos com muitos parâmetros ou não estacionários obtêm maior penalidade. (DURBIN E
KOOPMAN, 2008).
Para a estimação dos coeficientes, com o R, podem ser utilizados os métodos de máxima
verossimilhança e soma de quadrados condicional.
De acordo com Morettin e Toloi (2004), a fase crítica do procedimento é a identificação,
quando vários pesquisadores podem identificar modelos diferentes para uma mesma série
temporal. Desta forma, pode-se identificar não um único modelo, mas alguns modelos que
serão estimados e verificados.
Neste artigo, o teste de Shapiro-Wilk é usado para testar a normalidade dos resíduos.
3 R – Statistical Computing Environment and Graphics
O R é uma linguagem e ambiente para computação estatística. É um projeto GNU (General
Public License da Free Software Foundation), fornecendo uma ampla variedade de análises
estatísticas e técnicas gráficas, oferecendo o código aberto como uma rota para a participação
nessa atividade. O R pode ser compilado e executado em um grande número de plataformas
UNIX e sistemas semelhantes (como o FreeBSD e Linux), Windows e MacOS. (MILESKI,
2007).
De acordo com Hyndman e Racine (2002), a sintaxe da linguagem R é muito similar à
linguagem S. As principais diferenças são as seguintes: no R, os objetos não são salvos em
arquivos, mas armazenados internamente; os pacotes disponíveis em R são diferentes das
bibliotecas disponíveis no S-Plus.
Entre outras possibilidades, o R possui uma bem desenvolvida, simples e eficaz linguagem de
programação, que inclui condicionantes, loops, funções recursivas definidas pelo usuário e
facilidades de entrada e saída. Para tarefas computacionais intensivas, o C, o C++ e o Fortran
podem ser conectados e acessados durante a execução.
No ambiente R, a análise estatística é geralmente feita como uma série de passos, com os
resultados intermediários sendo armazenados.
Segundo Hyndman e Racine (2002), poderosos softwares comerciais como Matlab, Minitab,
SAS, SPSS e S-Plus, entre outros, podem ser caros, às vezes com restrições de licença e uso.
5
Assim, o R surgiu como uma alternativa a estes softwares comerciais e, como eles, é
dinâmico e evolui continuamente.
4 Proposta do artigo
O Software comercial SPSS, por exemplo, possui o Expert Modeler, que sugere um modelo
de ajuste, entre modelos de alisamento exponencial e ARIMA/SARIMA. O Expert Modeler
automaticamente identifica e estima o melhor ajuste, baseado em uma medida como o BIC
normalizado.
Hyndman e Khandakar (2008) descreveram dois algoritmos de predição implementados no
pacote forecast do R, baseados em modelos de alisamento exponencial e de ARIMA,
aplicados tanto em séries sazonais como não sazonais. Especificamente ao modelo ARIMA, o
algoritmo auto.arima retorna o melhor modelo de acordo com os critérios AIC, AICc e BIC.
O algoritmo realiza uma pesquisa sobre o modelo possível dentro das limitações da ordem
fornecida e retorna o modelo com o menor critério especificado.
Tanto o expert modeler quanto o auto.arima apresentam o melhor modelo segundo algum
critério. Para comparar dois modelos específicos, são necessárias mais de uma execução.
Assim, na tentativa de facilitar o estágio de identificação do modelo final a ser ajustado,
propõe-se a construção de uma função a ser executada no R que ajuste vários modelos
ARIMA(p, d, q), retornando uma matriz que traga nas linhas o número de modelos analisados
e nas colunas as medidas/informações sobre os modelos. A função, no ambiente R, não
necessita de instalação do pacote forecast.
O comando para execução da função requer valores máximos para as ordens “p”, “q” e “d” e
o método de ajuste (máxima verossimilhança, soma de quadrados condicional ou máxima
verssomilhança/soma de quadrados condicional), conforme opções da função arima, no
ambiente R. Caso não seja preenchido nada para “método”, a função ajusta de acordo com a
soma de quadrados condicional.
A função proposta ajusta (p+1)x(q+1)x(d+1) modelos a uma série temporal, sendo que
começa com o modelo ARIMA(0, 0, 0), seguindo com o modelo ARIMA(1, 0, 0) e termina
no modelo (p, d, q). O valor máximo aceitável para “p” e “q” é 6 e para “d” é 2; porém,
acredita-se que são valores exagerados. Poder-se-ia tentar modelos com ordens menores.
A matriz resultante é ordenada segundo o valor de AIC – Akaike information criteria –
quando é possível utilizar o método de máxima verossimilhança, ou pelo mape – média
absoluta percentual dos erros, quando for utilizado o método CSS. Além da matriz, a função
ainda fornece os gráficos da série, da autocorrelação e da autocorrelação parcial da série
original. Fornece, também, os valores da estatística e o p-value do teste de Dickey-Fuller para
a série original, para a 1ª e 2ª diferenças.
Após a obtenção da matriz, o analista pode escolher alguns (poucos) modelos candidatos e
realizar os demais testes e refazer os intervalos de confiança para os coeficientes.
A proposta é a de se obter, a partir da execução de uma linha de comando, uma matriz com
treze colunas:
 “p”: ordem de autorregressão;
 “d”: ordem de integração;
 “q”: ordem de médias móveis;
 sigma2: 2 estimada;
6
 loglik: o log da máxima verossimilhança (dos dados diferenciados ou não);
 Aic: Akaike information criteria;
 n_coef: o número de coeficientes testados no modelo (p+q para modelos com integração e
p+q+1 quando d=0);
 fora: número de coeficientes fora do intervalo coef ± 1.96 . Considera-se que o
coef
número de observações da série seja grande para utilizar a aproximação normal. Na
análise exploratória, o grau de confiança é fixado em 95%;
 Ljung: p-value do teste de Ljung realizado para o lag=1;
 Shapiro: p-value do teste de Shapiro-Wilk realizado com os resíduos;
 c_fora: número de autocorrelações dos resíduos diferentes de zero, segundo o proposto
por Morettin e Toloi (2004). Considerando que a matriz resultante é apenas uma análise
exploratória e não dispensa testes posteriores;
 mape: média absoluta percentual dos resíduos;
 Bic: Bayesian information criteria.
Em relação ao número de coeficientes, considera-se que se a ordem de “d”, número de
diferenças, for maior ou igual a 1, o número de coeficientes é igual a p+q; caso contrário,
p+q+1. Segundo Morettin e Toloi (2004), pode-se supor que, quando d>0, o termo constante
no modelo pode ser igual a zero.
As colunas sigma2, loglik e Aic são obtidas do resultado da execução da função arima no R,
enquanto os testes de Ljung e Shapiro são calculados usando as funções de Box.text e
Shapiro.test, respectivamente. As colunas n_coef, fora, c_fora, mape e Bic são calculadas na
função proposta em Anexo.
A função foi construída na versão 2.10.0.
4.1 Teste da Função
A função pode ser executada em séries univariadas, variáveis regressoras ou valores perdidos
(missing values). A função não analisa os outliers.
A função pode ser chamada da seguinte forma:
ajuste(x=serie, p=p, d=d, q=q, method=(“ML”, “ML-CSS”, “CSS”))
Onde:
 x é a série temporal;
 p, d e q são as ordens máximas;
 method é o método selecionado. A função arima do R dispõe das formas de ajuste de
máxima verossimilhança e soma de quadrados condicional. O padrão, a não ser que haja
valores perdidos, é a utilização da soma de quadrados condicional para localizar os
valores iniciais e então usar a máxima verossimilhança.
A seguir o teste da função em 3 séries pesquisadas na literatura. Não é objetivo do artigo
propor um modelo para as séries, mas fornecer informações para a escolha de um.
4.1.1 Série Nile
São 100 medições da vazão anual do rio Nilo no período de 1871–1970. A série está no banco
de dados do R.
7
Comando:
fit_A<-ajuste(st=Nile, p=6, d=1, q=6, method="ML")
Ou seja, o ajuste de (6+1)*(1+1)*(6+1) modelos de Box & Jenkins aplicados à mesma série.
A função retornou a estatística do teste de Dickey-Fuller e o p-value, a seguir.
 estatística e p-value do Teste de Dickey-Fuller - série original: -6.690051 0.01
 estatística e p-value do Teste de Dickey-Fuller - 1ª diferença: -16.93959 0.01
 estatística e p-value do Teste de Dickey-Fuller - 2ª diferença: -32.16724 0.01
A realização do teste é para auxiliar na tomada de decisão de quantas diferenças tomar.
Os gráficos da série original a ACF e a PACF estão a seguir.

Figura 2 – Gráfico da série, da função de autocorrelação e autocorrelação parcial
Abaixo, os três modelos com menores AIC para a série. São as três primeiras linhas da matriz
resultante, porém, transposta. As demais linhas foram suprimidas.

Linha
8
Coluna 1 2 3
p 2 1 0
d 1 1 1
q 4 1 2
Sigma2 17480,38 19769,29 19912,62
Loglik -626,61 -630,63 -630,98
Aic 1267,23 1267,26 1267,96
Coeficientes estimados 6 2 2
Coeficientes estimados igual a zero 0 0 1
Box Ljung (p-value) 0,94 0,74 0,99
Shapiro test (p-value) 0,87 0,73 0,69
Número de autocorrelações residuais diferentes de zero 0 0 0
Mape 12,29 12,79 12,81
BIC 1285,39 1275,04 1275,74
Fonte: Resultado da aplicação da função
Tabela 1 – Resultado da aplicação da função (Matriz transposta das 3 primeiras linhas)
Com os resultados da função executada, é possível eliminar modelos que não são
interessantes por apresentar autocorrelação dos resíduos no lag=1, verificado pelo teste de
Box Ljung ou por apresentar autocorrelações dos resíduos diferentes de zero, por exemplo.
Embora a função tenha sido executada com sucesso no R, gerando uma matriz com 98 linhas
(98 modelos ajustados), ela gerou 9 avisos: 7 problemas de convergência e 2 notificando pelo
desvio padrão igual a NaN (Not a Number). Quando acontece o problema de convergência, a
função retorna NA (missing) para todas as medidas daquele modelo, exceto para as ordens
“p”, “d” e “q”. Quando a variância de algum dos coeficientes for NaN, o valor retornado para
coeficientes estimados igual a zero é NA, na matriz de resultados.
As demais séries testadas terão as saídas omitidas.
4.1.2 Série B
São 114 observações do preço médio mensal (corrente) recebido pelos produtores de café,
exposto em Morettin e Toloi (1981).
Comando:
fit_b<-ajuste(st=B, p=6, d=2, q=6, method="ML")
A função também foi executada com sucesso, gerando 147 linhas (modelos). Como na série
anterior, a função gerou 25 avisos acusando problemas de convergência.
4.1.3 Série S
Série obtida de Makridakis et al. (apud Pellegrini, 2000), referente a usuários conectados a
um servidor da internet durante um período de 100 minutos.
Comando:
fit_c<-ajuste(st=S, p=6, d=2, q=6, method="ML")
A função retornou os gráficos e as saídas dos testes de Dickey-Fuller, porém acusou erro. As
alternativas foram ou reduzir as ordens de “p”, “d” e “q” ou mudar o método para “CSS”-
soma de quadrados condicionais. O comando utilizado novamente foi:
fit_d<-ajuste(st=S, p=6, d=2, q=6, method="CSS")
9
A função foi executada com sucesso, retornando 29 avisos acusando problemas de
convergência. A desvantagem de utilizar o método “CSS” em vez de “ML” é que aquele
retorna NA para os critérios AIC e BIC.
Quando a série original apresenta zeros, o mape retorna NA.
5 Considerações finais
A função executada no R apresenta as seguintes limitações:
Pode ser aplicada apenas em séries univariadas, sem componentes sazonais, variáveis
regressoras e valores perdidos.
A função ainda pode incorporar análises complementares, como a detecção de valores
extremos e formas de comparar as previsões com os modelos ajustados, além de outras
medidas importantes na comparação de modelos.
Além das implementações elencadas nos parágrafos precedentes, é necessário resolver os
avisos (warnings) posteriormente, sendo esses os desafios futuros.
Por outro lado, a função proposta pode auxiliar na identificação quando se ajusta um modelo
de Box & Jenkins a uma série temporal, à medida que permite comparar vários modelos
simultaneamente, permitindo descartar aqueles que não são interessantes.

Referências
BOX, G.E.P. & JENKINS, G.M. Time series analysis: forecasting and control. San Francisco: Holden-Day,
1970.
DIAS, G.J.C. Planejamento estratégico de um centro de distribuição: uma aplicação de redes neurais
artificiais de funções de bases radiais para previsão de séries temporais. Dissertação (Métodos Numéricos em
Engenharia). Universidade Federal do Paraná, 2006.
DURBIN J. & KOOPMAN S.J. Time series analysis by state space methods. Oxford: Oxford University Press,
2001. Reimpreted 2008.
HYNDMAN, R. & RACINE, J. Using R to teach econometrics. Journal of Applied Econometrics. Vol. 17, p.
175-189, 2002.
HYNDMAN, R. & KHANDAKAR, Y. Automatic time series forecasting: the forecast package for R. Journal
of Statistical Software. Vol. 27, 2008.
MILESKI, A.J. Análise de métodos de previsão de demanda baseados em séries temporais de uma empresa do
setor de perfumes e cosméticos. Dissertação (Engenharia de Produção e Sistemas). Pontifícia Universidade
Católica do Paraná, 2007.
MORETTIN, P.A. TOLOI, C.M.C. Modelos para previsão de séries temporais. Rio de Janeiro: IMPA/CNPQ,
1981. v. 2.
MORETTIN, P.A. TOLOI, C.M.C. Análise de séries temporais. São Paulo: Blucher, 2004.
PELLEGRINI, F.R. Metodologia para implementação de sistemas de previsão de demanda. Dissertação
(Engenharia de Produção). Universidade Federal do Rio Grande do Sul, 2000.
R Development Core Team (2009). R: A language and environment for statistical computing. R Foundation
for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.
SHUMWAY, R.H. Applied statistical time series analysis. New Jersey: Prentice-Hall, 1988.
SPSS Inc. TIME SERIES ANALYSIS AND FORECASTING WITH SPSS 14 FOR WINDOWS. Copyright 2006
by SPSS Inc. Printed in the United States of America.
WERNER, L. & RIBEIRO, J.L.D. Modelo composto para prever demanda através da integração de
10