ProcessoINAR(1)comEstruturaSazonalpara ... · Medeiros, Rodrigo Matheus Rocha de. Processo INAR(1) com estrutura sazonal para séries temporais de valores inteiros com sobredispersão

UNIVERSIDADE FEDERAL DO RIO GRANDE DO NORTE

CENTRO DE CIÊNCIAS EXATAS E DA TERRA

DEPARTAMENTO DE ESTATÍSTICA

Rodrigo Matheus Rocha de Medeiros

Processo INAR(1) com Estrutura Sazonal paraSéries Temporais de Valores Inteiros com

Sobredispersão

Natal - RN

Dezembro de 2017

Rodrigo Matheus Rocha de Medeiros

Processo INAR(1) com Estrutura Sazonal para SériesTemporais de Valores Inteiros com Sobredispersão

Monografia de Graduação apresentada ao De-partamento de Estatística do Centro de Ci-ências Exatas e da Terra da UniversidadeFederal do Rio Grande do Norte como re-quisito parcial para a obtenção do grau deBacharel em Estatística.

UNIVERSIDADE FEDERAL DO RIO GRANDE DO NORTE

CENTRO DE CIÊNCIAS EXATAS E DA TERRA

DEPARTAMENTO DE ESTATÍSTICA

Orientador: Prof. Dr. Marcelo Bourguignon Pereira

Natal - RNDezembro de 2017

Medeiros, Rodrigo Matheus Rocha de. Processo INAR(1) com estrutura sazonal para séries temporaisde valores inteiros com sobredispersão / Rodrigo Matheus Rochade Medeiros. - 2017. 57 f.: il.

Monografia (Bacharelado em Estatística) - UniversidadeFederal do Rio Grande do Norte. Centro de Ciências Exatas e daTerra. Departamento de Estatística. Natal, RN, 2017. Orientador: Marcelo Bourguignon Pereira.

1. Estatística - Monografia. 2. Dados inteiros não-negativos- Monografia. 3. Operador thinning - Monografia. 4. Processoautorregressivo - Monografia. 5. Processos de contagem -Monografia. 6. Sazonalidade - Monografia. I. Pereira, MarceloBourguignon. II. Título.

RN/UF/CCET CDU 519.2

Universidade Federal do Rio Grande do Norte - UFRNSistema de Bibliotecas - SISBI

Catalogação de Publicação na Fonte. UFRN - Biblioteca Setorial Prof. Ronaldo Xavier de Arruda - CCET

Dedicado à Maria José de Medeiros

Agradecimentos

Agradeço primeiramente a Deus por toda força de vontade, perseverança e paciên-cia que tive durante o curso de graduação.

À minha família, minha mãe Isabel, minha irmã Ruthe e meu pai Robson, por teracreditado em mim e me dado todo apoio necessário para concluir esta importante etapa.Agradeço à minha avó Maria José que sempre me incentivou a estudar, mas, infelizmente,perdeu suas principais lembranças tornando o caminho mais tortuoso e difícil para todosda minha família, porém, ainda assim me ensina que nada é fácil nesta vida.

Agradeço à minha namorada Brendda, minha inspiração e fonte de motivaçãodiária, que sempre esteve ao meu lado nos dias bons e também nos maus, tendo umaenorme paciência e sempre entendendo todo meu esforço. Agradeço por cada momentoque passo ao seu lado.

Agradeço ao professor Pledson por todo o incentivo que recebi durante meusprimeiros meses da graduação, pela palestra, sobre o curso de estatística, que ministrouem minha escola do ensino médio que fez com o que eu me interessasse por esta incrívelárea. Por ter ido ao batalhão 17o GAC do exército conversar com o Major Filgueiras paraque eu pudesse continuar no curso sem atrasá-lo, e conseguir. Pelas confraternizações entrealunos e professores.

Agradeço aos professores do departamento de estatística, departamento de ma-temática, departamento de ciência da informação, departamento de educação física, de-partamento de línguas e literaturas estrangeiras modernas e departamento de informáticacom quem tive grande prazer de ser aluno. À professora Dione, que me acompanhou comotutora do PET desde o meu primeiro ano, por todos os conselhos e inspiração. Ao meuorientador Marcelo, por ter aceitado me orientar, pelas conversas, apoio e por todos osconselhos que recebi.

Agradeço a todos os meus amigos que fiz durante a graduação. Aos amigos daminha turma de 2014 que caminharam ao meu lado durante o curso. A todos os talentososamigos que fiz no PET com quem dividi alegrias e preocupações ao longo das tardes, sendoeles das turmas de 2012, 2013, 2014, 2015, 2016 e agora 2017, em especial Erika, Felipe eLucas, com quem passei a maior parte do tempo, obrigado por todo conhecimento queadquiri. Agradeço aos amigos da banda Poty Brigade pela paciência e toda inspiração querecebi.

Agradeço ao PET, pela ajuda financeira, mas principalmente por todo conheci-mento adquirido, e todo incentivo na busca por excelência em um curso tão complicado.

ResumoO estudo de modelos para séries temporais de valores inteiros está cada vez mais presentena literatura. É comum, na prática, que séries contenham componente sazonal, e aocontrário dos modelos contínuos, processos de contagem com estrutura sazonal não têmrecebido muita atenção na literatura até o momento. Os objetivos desta pesquisa são:introduzir um novo modelo autorregressivo para dados inteiros não-negativos com estruturasazonal e que apresentem sobredispersão, definir as principais propriedades do processo,estudar métodos de estimação dos parâmetros do modelo proposto, que são os estimadoresde Yule-Walker, mínimos quadrados condicionais e máxima verossimilhança condicional,compará-los em um estudo de simulação e aplicar o modelo proposto a um conjunto dedados reais. Comparamos o novo modelo com modelos já propostos na literatura e omodelo proposto neste estudo obteve um melhor ajuste.

Palavras-chave: Dados inteiros não-negativos. Operador thinning. Processo autorregres-sivo. Processos de contagem. Sazonalidade.

AbstractThe study of integer-valued time series models is increasingly present in the literature. Itis common in practice for series to contain a seasonal component, and unlike continuousmodels, counting processes with a seasonal structure have not received much attention inthe literature so far. The aims of this paper are: to introduce a new autoregressive modelfor non-negative integer-valued time series with seasonal structure which are overdispersed,to define the main process properties, to study the estimation methods of the parametersof the proposed model, these methods are Yule-Walker, conditional least squares andconditional maximum likelihood estimators, comparing them in a simulation study andfinally applying the proposed model to a real data set. We compared the new model withmodels already proposed in the literature, and, the model proposed in this paper presenteda better fit.

Keywords: Autoregressive process. Counting processes. Non-negative integer-valued.Seasonality. Thinning operators.

Sumário

1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.1 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.2 Desenvolvimento dos capítulos . . . . . . . . . . . . . . . . . . . . . . 12

2 OPERADOR THINNING BINOMIAL NEGATIVO . . . . . . . . . . 132.1 Operador thinning binomial negativo . . . . . . . . . . . . . . . . . . 142.2 Propriedades do operador thinning binomial negativo . . . . . . . . . 15

3 CONSTRUÇÃO DO PROCESSO NGINAR(1)s . . . . . . . . . . . 173.1 Processo NGINAR(1)s . . . . . . . . . . . . . . . . . . . . . . . . . . 173.2 Propriedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.2.1 Momentos ordinários e medidas condicionais . . . . . . . . . . . . . . . . . 203.2.2 Função de autocovariância e autocorrelação . . . . . . . . . . . . . . . . . 213.2.3 Função de Transição . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4 MÉTODOS DE ESTIMAÇÃO PARA OS PARÂMETROS . . . . . . 254.1 Estimador de máxima verossimilhança condicional . . . . . . . . . . . 254.2 Estimador de mínimos quadrados condicionais . . . . . . . . . . . . . 264.3 Estimador de Yule-Walker . . . . . . . . . . . . . . . . . . . . . . . . . 27

5 ESTUDO DE SIMULAÇÃO . . . . . . . . . . . . . . . . . . . . . . . 295.1 Análise dos cenários . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

6 APLICAÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 346.1 Análise descritiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 346.2 Ajustes dos modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

REFERÊNCIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

APÊNDICE A – ESTIMADOR DEMÍNIMOS QUADRADOS CON-DICIONAIS . . . . . . . . . . . . . . . . . . . . . . 43

APÊNDICE B – ROTINAS COMPUTACIONAIS . . . . . . . . . . 44

9

1 Introdução

Existem várias situações práticas em que é necessário recorrer à informações dopassado para obter respostas no presente. Seguindo a noção descrita em Antunes e Cardoso(2015), imagine um empresário que pretende investir em um tipo de mercado em umadeterminada cidade. É essencial que ele realize um estudo a cerca do comportamentohistórico do mercado já existente nessa cidade, para que se possa minimizar eventuaisprejuízos. Muitas vezes também, é necessário antever o futuro e assim antecipar possíveisresultados. Na área de epidemiologia, o esforço para prever futuros cenários é indispensávelpara a redução da carga de doenças na população (ANTUNES; CARDOSO, 2015). Oconjunto de técnicas estatísticas que lidam com esse tipo de problema, se chama análisede séries temporais.

Podemos definir uma série temporal, de maneira informal, por uma coleção deobservações de uma variável aleatória obtidas sequencialmente em instantes no tempo, quepossui como principal característica a presença de dependência entre as observações (BOX;JENKINS; REINSEL, 1994). Ao realizar uma análise de uma série temporal estamosinteressados em descrever, interpretar e realizar previsões acerca de fenômenos de interesseque podem ser aproximados por modelos.

Quando a variável estudada é contínua, i.e., assume valores em um conjunto nãoenumerável, dizemos que a série temporal possui marginal contínua, e caso contrário,dizemos que possui marginal discreta, que muitas vezes será referenciada neste trabalhocomo uma série temporal de valores inteiros, ou tão somente série discreta. Os principaismodelos teóricos desenvolvidos para modelar séries temporais, são baseados em distribuiçõesde probabilidade contínuas (BOX; JENKINS; REINSEL, 1970) que, consequentemente,assumem que a variável de interesse é contínua.

Dados autocorrelacionados discretos podem surgir de contagens de eventos ouindivíduos ao longo do tempo e são bastante comuns na prática, como por exemplo, emfinanças e economia, para modelar o número de transações de seguradoras sul-coreanas(WEIß; KIM, 2013), em estudos sobre trânsito, para análise e prevenção a respeito donúmero de acidentes de trânsito em rodovias do Reino Unido (QUDDUS, 2008), emmedicina para modelar a incidência de um certo tipo de doença contagiosa em Montreal(CARDINAL; ROY; LAMBERT, 1999). Esse tipo de dado é conhecido como processo decontagem.

Como dito anteriormente, em uma modelagem usual de séries temporais, assume-se que a distribuição marginal da série é contínua e, em geral, com distribuição normal.Dados discretos nem sempre podem ser ajustados por modelos contínuos. Em alguns casos

Capítulo 1. Introdução 10

como em que a magnitude dos valores distintos da série temporal discreta são grandes,a modelagem convencional é uma boa aproximação, (PEREIRA, 2012). Contudo, paradados em que não é possível utilizar a metodologia usual, surge a motivação para novaspropostas que sejam capazes de modelar as séries com marginais discretas.

Os principais modelos para séries temporais de valores inteiros surgiram comoadaptações dos modelos contínuos. O interesse nesses modelos se deu ao final da décadade 70 com o modelo discreto autorregressivo de média móvel (DARMA), proposto emJacobs e Lewis (1977,1978), que é um processo formado por uma combinação linear devariáveis aleatórias discretas independentes e identicamente distribuídas (i.i.d.). Pode-secitar também McKenzie (1986) que propõe modelos autorregressivos de média móvel comdistribuições marginais geométrica e binomial-negativa.

Dentre a nova classe de modelos para séries temporais de valores inteiros quesurgiram, destaca-se o modelo autorregressivo de primeira ordem para dados inteiros,[INAR(1)], proposto por McKenzie (1985) e Al-Osh e Alzaid (1987) de maneira inde-pendente. A estrutura de autocorrelação do processo INAR(1) é a mesma do processoautorregressivo de ordem 1 [AR(1)]. Al-Osh e Alzaid (1987) define o processo INAR(1),bem como suas principais propriedades e métodos de estimação para os parâmetros domodelo.

Em pesquisas recentes, podemos citar Khoo, Ong e Biswas (2017), que propõeum novo modelo de primeira ordem para processos estacionários não negativos discretosbaseados nos operadores Pedram e thinning, Kim e Jun (2017), que propõe o modelo deheterocedasticidade condicional autorregressiva generalizada (GARCH) com distribuiçãomarginal Poisson e binomial, e Bourguignon e Weiß (2017) que propõe um modelo paraséries discretas com equidispersão, sobredispersão ou subdispersão.

Quando um novo modelo estatístico é proposto, é desejável que este capte asprincipais características do fenômeno de estudo. Um dos fatores que influenciam ocomportamento das séries temporais é a sazonalidade. Dados sazonais são comuns naprática, áreas como economia (FRANSES et al., 1996), climatologia (LUTERBACHER etal., 2004), epidemiologia (GRASSLY; FRASER, 2006), entre outras, são exemplos em quesão estudados fenômenos sazonais. Conforme Pereira (2012), a análise da variação sazonalvem sendo analisada desde o final da década de 20, como por exemplo, nos trabalhos deMitchell et al. (1927) e Macaulay et al. (1938).

Uma série temporal apresenta padrões sazonais quando um determinado fenômenose comporta de maneira semelhante em períodos regulares no tempo. Por exemplo, quandodeterminado produto agrícola é influenciado pelas estações ou determinadas épocas doano, seus preços tendem a seguir o padrão sazonal. Chamamos de período sazonal o menorintervalo de tempo entre duas ocorrências de tal fenômeno e o denotaremos neste trabalho


por “s”. O período sazonal resume o comportamento sazonal de uma série temporal, e.g.,dados diários podem apresentar período s = 7, dados mensais podem ter período s = 12,entre outros.

Até o momento, ao contrário dos estudos da componente sazonal em séries commarginais contínuas, poucas pesquisas foram realizadas a fim de estudar os efeitos dasazonalidade em processos de contagem. Em pesquisas recentes podemos citar comoexemplo, Awale (2017) que apresenta o desenvolvimento de um novo modelo com estruturasazonal na área de epidemiologia e Bourguignon et al. (2016) que define e apresenta asprincipais propriedades de um modelo autorregressivo com estrutura sazonal baseado nomodelo INAR(1) com marginal Poisson.

Bourguignon et al. (2016) mostra que o modelo autorregressivo de primeira ordemde valores inteiros e período sazonal s, [INAR(1)s], obtém melhor ajuste que o modeloINAR(1) para um cojunto de dados específico que possui estrutura sazonal. Entretanto,ambos os modelos supõem que a distribuição de probabilidade marginal do processo decontagem segue distribuição de Poisson. Ao utilizar este modelo probabilístico, assume-seque a média e a variância do processo são iguais, fenômeno chamado de equidispersão,algo que não ocorre com frequência em aplicações reais.

Quando a série apresenta sobredispersão, i.e., a variância do processo estudado émaior do que a média, temos indícios de que o modelo de Poisson não seja uma boa escolha,por exemplo, (WEIß, 2008), (WEIß, 2010), (ZHU, 2012) e (BOURGUIGNON; WEIß,2017). Uma das possíveis maneiras de modelar a sobredispersão em processos de contagemé optar por distribuições de probabilidade marginais que assumam que a variância é maiorque a média do processo, como por exemplo a distribuição geométrica. Alguns estudos queconsideraram tal distribuição como a distribuição marginal do processo foram McKenzie(1985,1986), Al-Osh e Aly (1992) e Ristić, Bakouch e Nastić (2009).

Pesquisas que abordam séries temporais de valores inteiros com estrutura sazonale sobredispersão não têm sido discutidas na literatura, em decorrência deste fato, nestetrabalho, estudaremos os efeitos da sazonalidade no processo NGINAR(1) através daconstrução de um novo modelo que considera a componente sazonal em sua estrutura.

1.1 ObjetivosO principal objetivo deste trabalho, é desenvolver um novo modelo de séries

temporais para dados discretos não-negativos, com sobredispersão e que apresentemcomportamento sazonal com base na estrutura desenvolvida em Ristić, Bakouch e Nastić(2009).

Resumimos os objetivos específicos como


• Definir e estudar as principais propriedades do operador thinning binomial negativo;

• Definir o novo modelo e suas principais propriedades tais como função de transição,função de autocovariância e autocorrelação;

• Desenvolver e estudar métodos de estimação para os parâmetros do processo pro-posto, que são os estimadores de Yule-Walker, estimadores de mínimos quadradoscondicionais e os estimadores de máxima verossimilhança condicional;

• Realizar um estudo de Monte Carlo para comparar os estimadores propostos, avali-ando em diferentes cenários suas performances;

• Aplicar o modelo proposto a um conjunto de dados reais e compará-lo com outrosmodelos já existentes na literatura.

1.2 Desenvolvimento dos capítulosEste trabalho se desenvolverá da seguinte forma: no Capítulo 2 realizaremos uma

introdução motivacional a cerca da definição do operador thinning binomial negativo,apresentaremos sua definição e suas principais propriedades conforme apresentadas emRistić, Bakouch e Nastić (2009). No Capítulo 3, realizaremos a construção do processoproposto, estudaremos sua condição de estacionariedade, e demonstraremos suas principaispropriedades como a função de autocorrelação, distribuição do processo de inovação domodelo e a função de transição. No Capítulo 4, propomos métodos de estimação paraos parâmetros do modelo. No Capítulo 5, apresentaremos um estudo de simulação feitocom a finalidade de observar a performance dos estimadores apresentados. Finalmente, noCapítulo 6, aplicaremos o modelo proposto neste trabalho a um conjunto de dados reais,com a finalidade de realizar comparações entre os ajustes de outros modelos presentes naliteratura.

13

2 Operador Thinning Binomial Negativo

A principal adaptação do modelo INAR(1), em relação aos modelos clássicos deséries temporais, é devida a utilização do operador thinning binomial proposto em Steutele Harn (1979), que é definido por uma soma aleatória de variáveis aleatórias i.i.d. comdistribuição Bernoulli, em que essas variáveis são chamadas de série de contagem. Osmodelos de séries temporais para dados contínuos baseiam-se em operações multiplicativasque frequentemente não retornam valores inteiros, como por exemplo o modelo AR(1)

Xt = δ + φXt−1 + εt, t ∈ Z,

em que δ, φ ∈ R são parâmetros e {εt}t∈Z é uma sequência de variáveis aleatórias i.i.d. commédia 0 e variância constante. Mesmo que Xt seja uma variável aleatória discreta, com aevolução do processo, a série não retornará valores inteiros.

Al-Osh e Alzaid (1987) e McKenzie (1985) propuseram utilizar o operador thinningbinomial para que se mantivesse a estrutura de autocorrelação e que os valores retornadosa partir do processo de evolução do modelo fossem inteiros. Gauthier e Latour (1994) apre-sentam uma generalização dos operadores thinning permitindo que as séries de contagemsigam qualquer distribuição discreta, o que possibilita desenvolver modelos com diferentestipos de operadores thinning e a comparação entre eles.

Algumas pesquisas que utilizaram os operadores thinning quando assume-se quedistribuição marginal da série é geométrica, foram McKenzie (1985,1986) que utiliza ooperador thinning binomial e Al-Osh e Aly (1992) que utiliza uma soma de variáveis alea-tórias i.i.d. geométrica para a série de contagem, semelhante ao operador que será definidoneste trabalho, porém se diferencia apenas por considerar uma estrutura condicional auma variável aleatória binomial.

Mais recentemente, Ristić, Bakouch e Nastić (2009) propuseram o operadorthinning binomial negativo na definição do processo NGINAR(1). O operador é definidopor uma soma aleatória de variáveis aleatórias i.i.d. em que cada uma dessas variáveispossui distribuição geométrica. Os objetivos deste capítulo são definir o operador thinningbinomial negativo e apresentar suas principais propriedades.

Capítulo 2. Operador Thinning Binomial Negativo 14

2.1 Operador thinning binomial negativo

Definição 2.1.1. Seja X uma variável aleatória discreta não-negativa e α ∈ [0,1), ooperador thinning binomial negativo, denotado por “∗”, é definido como

α ∗X =

X∑i=1

Wi, se X > 0,

0, se X = 0,

em que a série de contagens {Wi}i∈N é uma sequência de variáveis aleatórias i.i.d. comdistribuição geométrica com média α ∈ [0,1).

Seja W variável aleatória com distribuição geométrica com média α, a função deprobabilidade W é dada por

P (W = w) = αw

(α + 1)w+1 , w = 0,1,2,...,

em que E[W ] = α e V ar[W ] = α(1 +α). Se considerarmos uma sequência finita de ensaiosde Bernoulli expressa em termos de sucessos ou falhas, W representa o número falhasaté o primeiro sucesso com probabilidade de sucesso 1/(α + 1). Utilizamos a notaçãoW ∼ Geo(1/(α + 1)), para denotar a variável aleatória W com distribuição geométricacom média α.

Proposição 2.1.1. Sejam W1, ...,Wk uma amostra aleatória de variáveis aleatórias i.i.d.com Wi ∼ Geo(1/(α + 1)), i = 1,2, ..., k então

V = α ∗X | (X = k) =k∑i=1

Wi ∼ BN (k, α/(1 + α)) ,

em que a notação V ∼ BN(a, b) representa que a variável aleatória V possui distribuiçãobinomial negativa com parâmetros a e b.

Demonstração. Seja φW (t) = E[etW ] a função geradora de momentos da variável aleatóriaW . Seja p = 1/(α + 1), ou seja, W ∼ Geo(p), então

φW (t) = p

1− (1− p)et .

Desejamos obter a distribuição da variável aleatória V =k∑i=1

Wi, como por definição,W1, ...,Wk é uma amostra aleatória i.i.d., segue que

φV (t) =k∏i=1

φWi(t) =

[p

1− (1− p)et

]k,

que é igual a função geradora de momentos de uma variável aleatória com distribuiçãobinomial negativa com parâmetros k e 1− p, ou seja

V ∼ BN(k,α/(1 + α)),

como queríamos demonstrar.


Sendo assim, a função de probabilidade de V é expressa por

P (V = v) = v + k − 1

v

αv

(α + 1)k+v , v = 0,1,..., (2.1)

com E[V ] = kα e V ar[V ] = kα(α + 1).

A distribuição da variável aleatória V é de grande relevância para a construção doprocesso, uma vez que a distribuição obtida a partir do novo operador thinning binomialnegativo é a distribuição binomial negativa reparametrizada como em (2.1). Além disso,com o conhecimento da distribuição de V , é possível simular valores do processo proposto,e com isso avaliar possíveis estimadores e principais propriedades destes.

2.2 Propriedades do operador thinning binomial negativoAbaixo, no Lema 2.2.1, são apresentadas as propriedades do operador thinning

binomial negativo. As propriedades e demonstrações estão descritas em Ristić, Bakouch eNastić (2009).

Lema 2.2.1. Suponha que a série de contagem de α ∗X é independente de X e de Y e asséries de contagem de αi ∗Xi, i=1,2,...,r, são mutualmente independentes e independentesde Xi, i = 1,2,..., r . Então, o operador thinning binomial negativo ∗ tem as seguintespropriedades

(i) E[r∏i=1

(αi ∗Xi)]=

r∏i=1

αiE[r∏i=1

Xi

], r ≥ 1;

(ii) E [(α ∗X)2] = α2E [X2] + α(1 + α)E [X];

(iii) E [(α ∗X)3] = α3E [X3] + 3α2(1 + α)E [X2] + α(1 + α)(1 + 2α)E [X];

(iv) E [(α ∗X − α ∗ Y )2] = α(1+α)E [|X − Y |]+α2E [(X − Y )2], se os processosde contagem α ∗X e α ∗ Y possuírem a mesma distribuição geométrica com média α.

Levando em consideração o Teorema 2.2.1 abaixo, Ristić, Bakouch e Nastić (2009)mostra que as propriedades do operador thinning binomial negativo não se assemelham aspropriedades do operador thinning binomial, e ao contrário deste, apresenta propriedadesrelativamente menos intuitivas.

Teorema 2.2.1. Seja X variável aleatória em que X ∼ Geo(1/(α + 1)), então

(i)

1 ∗X =

0 com probabilidade 1/(1 + α),X com probabilidade α/(1 + α)2,

X + Y com probabilidade α/(1 + α)2,


em que Y ∼ Geo((1 + α)/(2 + α)) e Y é independente de X.

(ii)

β ∗ γ ∗X =

0 com probabilidade 1 + γ

1 + γ + γα,

(βγ) ∗X + Y1 com probabilidade γ2α2

(1 + γ + γα)(1 + γα) ,

(βγ) ∗X + Y2 com probabilidade γα

(1 + γ + γα)(1 + γα) ,

em que β e γ ∈ [0,1), e as variáveis aleatórias independentes Y1 e Y2 possuem, respectiva-mente, distribuição geométrica com parâmetros βγ/(1+βγ) e β(1+γ+γα)/[1+β(1+γ+γα)]e são independentes de X.

17

3 Construção do Processo NGINAR(1)s

Uma das principais características que, frequentemente, estão presentes em sériestemporais é a dependência entre as observações, que muitas vezes, mesmo com o passardo tempo, continua influenciando os valores da série. Um dos modelos clássicos queconsideram a memória do processo é o modelo autorregressivo proposto em Box, Jenkins eReinsel (1970), que podem ser extremamente úteis em situações práticas (BOX; JENKINS;REINSEL, 1994).

Conforme descrito em (BOX; JENKINS; REINSEL, 1994), no modelo AR(1), aatual observação do processo é expressa em uma estrutura linear em termos das obser-vações anteriores, em razão disto, a estrutura de correlação entre as observações decaiexponencialmente, sugerindo que as observações presentes serão influenciadas cada vezmenos por observações em defasagens de tempos cada vez maiores.

Neste capítulo, definiremos um novo processo autorregressivo de primeira ordempara dados inteiros não-negativos com distribuição marginal geométrica de período sazonals, [NGINAR(1)s]. Assim, consideramos a principal contribuição teórica deste trabalho, naárea de séries temporais de valores inteiros, como sendo o estudo da componente sazonalno modelo NGINAR(1). Os objetivos principais deste capítulo são definir o processo edemonstrar suas principais propriedades como condição de estacionariedade, funções deautocovariância e autocorrelação, assim como a função de transição.

3.1 Processo NGINAR(1)s

Definição 3.1.1. (Processo NGINAR(1)s). Um processo discreto de valores inteirosnão-negativos, {Yt}t∈Z, diz-se um processo sazonal autorregressivo de valores inteiros deordem 1 e período sazonal s ∈ N, NGINAR(1)s, se satisfaz a seguinte equação

Yt = α ∗ Yt−s + εt, t ∈ Z, (3.1)

em que α ∈ [0, 1), {Yt}t∈Z é uma sequência de variáveis aleatórias tal que {Yt}t∈Z possuidistribuição marginal geométrica com média µ > 0. {εt}t∈Z é uma sequência de variáveisaleatórias não negativas i.i.d. independentes da série de contagens {Wi}i∈N ∀i ∈ N e sãoindependentes de Yt−l, ∀ l ≥ s.

Pela construção do processo NGINAR(1)s, {Yt}t∈Z pode ser visto como umacadeia de Markov s-passos à frente, ou seja, a observação em um certo tempo t depende

Capítulo 3. Construção do Processo NGINAR(1)s 18

apenas da observação do tempo t− s, ou seja

P (Yt = yt | Yt−1 = yt−1, Yt−2 = yt−2, ..., Y0 = y0) = P (Yt = yt | Yt−s = yt−s). (3.2)

Além disso, {Yt}t∈Z é um processo homogêneo, i.e., as probabilidades de transiçãosão invariantes no tempo, ver na seção (3.2.3); irredutível, ou seja, cada valor assumidopela variável aleatória {Yt}t∈Z pode ser um valor da cadeia em um determinado tempot, e aperiódica, ou seja, não existe um t0 ∈ Z tal que P (Yt0 = i) = 1 ∀t > t0 para algumi ∈ {0,1,...}.

Quando o período sazonal s = 1, então {Yt}t∈Z resume-se ao processo NGINAR(1)proposto em Ristić, Bakouch e Nastić (2009), ou seja, o processo NGINAR(1) é um casoparticular do processo NGINAR(1)s, caso este que não leva em consideração a componentesazonal na série.

Conforme a Proposição 3.1 em Latour (1998), se {εt}t∈Z é uma sequência devariáveis aleatórias i.i.d. não-negativas e discretas, e α < 1 então o processo {Yt}t∈Zdefinido em (3.1) é estacionário. Em particular, neste trabalho estudaremos o caso emque α obedece a restrição α ∈

[0, µ

1+µ

), então, como εt ≥ 0 ∀t ∈ Z e α < µ/(1 + µ) < 1, o

processo NGINAR(1)s é estacionário.

{εt}t∈Z é chamado de processo de inovação, e acontece por fenômenos aleatóriosque não são explicados no modelo. Em técnicas de diagnóstico, conhecer a distribuiçãodo processo de inovação é necessário para verificar as suposições impostas pelo modeloutilizado. A distribuição de {εt}t∈Z é de interesse também para que se possa calcular afunção de transição do processo, como também para simular valores da série, o que é degrande relevância para estudarmos o comportamento do processo e de seus estimadores.Na Proposição 3.1.1, apresentamos a função de probabilidade de εt.

Proposição 3.1.1. εt é uma mistura de duas variáveis aleatórias, ambas com distribui-ção geométrica com parâmetros 1/(1 + µ) e 1/(1 + α), respectivamente, com função deprobabilidade dada por

P (εt = l) =(

1− αµ

µ− α

)µl

(1 + µ)l+1 +(

αµ

µ− α

)αl

(1 + α)l+1 , µ 6= α, l = 0 ,1 ,...,

em que µε = E[εt] = (1− α)µ e σ2ε = V ar[εt] = (1 + α)µ[(1 + µ)(1− α)− α].

Demonstração. Seja ϕX(r) a função geradora de probabilidades (fgp) da variável aleatóriaX, definida por

ϕX(r) = E[rX].

Se X possui distribuição geométrica de média µ, então sua fgp é dada por

ϕX(r) = 11 + µ− µr

.


Sejam ϕY (r), ϕW (r) e ϕε(r) as fgp’s das variáveis aleatórias Yt, Wi e εt, respecti-vamente. Pela independência entre {εt}t∈Z e a série de contagem, tem-se que

ϕY (r) = E[rYt

]= E

[rα∗Yt−s+εt

]= E

[rα∗Yt−s

]E [rεt ]

⇒ E[rYt

]= E

[E[rα∗Yt−s

∣∣∣Yt−s = yt−s]]E [rεt ] .

Por independência entre as variáveis da série de contagem {Wi}i∈N e a estaciona-riedade do processo

⇒ E[rYt

]= E

[E[rW]Yt−s

]E [rεt ]

⇒ E[rYt

]= E

[E[rW]Yt]E [rεt ] .

Segue queϕY (r) = ϕY (ϕW (r))ϕε(r). (3.3)

Suponha que Yt e Wi possuem distribuição geométrica com parâmetros 1/(1 + µ)e 1/(1 + α), respectivamente, obtemos então, por (3.3),

ϕε(r) = ϕY (r)ϕY (ϕW (r)) = 1 + α(1− r) + µα(1− r)

[1 + µ(1− r)][1 + α(1− r)] .

Utilizando frações parciais, temos que

ϕε(r) =(

1− αµ

µ− α

)1

(1 + µ− µr) +(

αµ

µ− α

)1

(1 + α− αr) ,

Logo, a fgp de εt é uma média ponderada de duas fgp’s de variáveis aleatóriascom distribuição geométrica de média µ e α, respectivamente. Para que P (εt = l) sejauma legítima função de probabilidade, temos a seguinte restrição

0 ≤ αµ/(µ− α) ≤ 1 ⇒ α ≤ µ/(1 + µ).

Assim,

P (εt = l) =(

1− αµ

µ− α

)µl

(1 + µ)l+1 +(

αµ

µ− α

)αl

(1 + α)l+1 , µ 6= α, l = 0 ,1 ,...

Resta demonstrar que µε = E[εt] = (1 − α)µ e σ2ε = V ar[εt] = (1 + α)µ[(1 +


µ)(1− α)− α]. Seja L = εt, com base na função de probabilidade de εt, temos que

E[εt] = E[L] =∞∑l=0

lP (L = l)

=(

1− αµ

µ− α

) ∞∑l=0

lµl

(1 + µ)l+1 +(

αµ

µ− α

) ∞∑l=0

lαl

(1 + α)l+1

=(

1− αµ

µ− α

)µ+

(αµ

µ− α

)α = µ(1− α).

V ar[εt] = V ar[L] = E[L2]− E2[L]

=∞∑l=0

l2P (L = l)− µ2(1− α)2

=(

1− αµ

µ− α

) ∞∑l=0

l2µl

(1 + µ)l+1 +(

αµ

µ− α

) ∞∑l=0

l2αl

(1 + α)l+1 − µ2(1− α)2

= µ2 + µ(1 + µ)− 2αµ2 − 2α2µ− αµ− µ2(1− α)2

= µ(1 + µ)− αµ− α2µ2 − 2α2µ+ αµ2 − αµ2

= (1 + α)µ[(1 + µ)(1− α)− α].

A Figura 1 apresenta seis realizações do processo NGINAR(1)s variando seusparâmetros e o período sazonal. Pode-se perceber a diferença (visual) da magnitudedos dados. Note que para valores pequenos da média, µ, a magnitude dos valores dasérie é pequena e possui muitos valores repetidos, uma das características encontradasem séries temporais de valores inteiros, contudo, quando há um aumento da média, epor consequência na variabilidade dos dados, a magnitude dos valores da série tambémaumenta.

3.2 Propriedades

3.2.1 Momentos ordinários e medidas condicionais

Os momentos ordinários de uma variável aleatória são de importância para quepossam ser calculadas propriedades acerca de sua distribuição de probabilidades, como porexemplo o valor esperado e a variância. Como, por definição, o processo NGINAR(1)spossui distribuição marginal geométrica, para todo t ∈ Z, os momentos ordinais de ordemk, µ′k, são dados por

µ′1 = µ, µ′2 = µ(2µ+ 1), µ′3 = µ[6µ(µ+ 1) + 1], µ′4 = µ(2µ+ 1)[12µ(µ+ 1) + 1].

Note que,V ar[Yt] = µ(µ+ 1), ∀t ∈ Z.


Tempo

Realizações

0 20 40 60 80 100

04

8

(a) Tempo

Realizações

0 20 40 60 80 100

0204060

(b)

Tempo

Realizações

0 20 40 60 80 100

04

8

(c) TempoRealizações

0 20 40 60 80 100

0204060

(d)

Tempo

Realizações

0 20 40 60 80 100

04

8

(e) Tempo

Realizações

0 20 40 60 80 100

0204060

(f)

Figura 1 – Cem realizações do processo NGINAR(1)s em que as variações dos parâmetrosusadas são (a) α = 0,3, µ = 1 e s = 7, (b) α = 0,3, µ = 10 e s = 7, (c) α = 0,3,µ = 1 e s = 12, (d) α = 0,3, µ = 10 e s = 12, (e) α = 0,5, µ = 1,5 e s = 7 e (f)α = 0,5, µ = 10 e s = 7.

Como estamos em um cenário em que admitimos a dependência entre as variáveis,devemos considerar também a estrutura condicional das propriedades estabelecidas. Épossível verificar que

E[Yt | Yt−1, ..., Y0] = E[Yt | Yt−s] = αYt−s + µε,

V ar[Yt | Yt−1, ..., Y0] = V ar[Yt | Yt−s] = α(α + 1)Yt−s + σ2ε ,

em que µε = E[εt] e σ2εt

= V ar[εt].

3.2.2 Função de autocovariância e autocorrelação

Uma das abordagens que é utilizada em análise de séries temporais é a análise nodomínio do tempo, e, em tal abordagem, a ênfase está na estrutura das correlações entreos eventos de interesse em diferentes instantes de tempo. Na metodologia Box-Jenkins,


a principal ferramenta para a identificação de um modelo para uma série temporal é afunção de autocorrelação e a autocorrelação parcial.

As funções de autocovariância e de autocorrelação do processo NGINAR(1)s sãodadas, respectivamente, por

Cov [Yt,Yt+k] = γ(k) =

αks γ(0), se k é múltiplo de s,

0, caso contrário,e

ρ(k) = γ(k)γ(0) =

αks , se k é múltiplo de s,0, caso contrário,

(3.4)

em que k é o número de períodos associados a uma observação precedente chamado dedefasagem.

Demonstração. Suponha que k é múltiplo de s, mais precisamente, ∃ p ∈ N; k = ps, então:

γ(k) = Cov [Yt, Yt+k] = E [YtYt+k]− E [Yt]E [Yt+k]

= E [Yt (α ∗ Yt+k−s + εt+k)]− E [Yt]E [(α ∗ Yt+k−s + εt+k)]

= E [Yt (α ∗ Yt+k−s)] + E [Ytεt+k]− E [Yt]E [α ∗ Yt+k−s]− E [Yt]E [εt+k]

=︸︷︷︸propriedade (i)

αCov [Yt+k−s,Yt] + Cov [εt+k,Yt] =︸︷︷︸independência dos ε′s

αγ(k − s)

Continuando esse processo recursivamente por p− 1 vezes, temos que:

Cov [Yt,Yt+k] = αpγ(k − ps) = αks γ(0), se k = ps.

Suponha agora que k não seja múltiplo se s. Pela propriedade de Markov, descritaem (3.2), a variável aleatória Yt é independente de Yt+k para todo t ∈ Z, o que implica que

Cov(Yt, Yt+k) = 0.

Logo, a função de autocorrelação é dada por

ρ(k) = γ(k)γ(0) =

αks , se k = ps,

0, caso contrário.

A Figura 2 apresenta a função de autocorrelação (FAC) amostral de realizaçõesdo processo NGINAR(1)s em que variamos o período sazonal com os valores s = 3, 4, 7 e12, que podem ser interpretados da seguinte maneira, se a série for tomada mensalmente,os períodos sazonais 3, 4 e 12 representam uma periodicidade trimestral, quadrimestrale anual, respectivamente, enquanto que se a série for tomada diariamente, um períodosazonal s = 7 denota uma periodicidade semanal. Note que dada uma amostra do processoNGINAR(1)s, através da FAC amostral conseguimos encontrar o período sazonal da série.


5 10 15 20 25 30

0.0

0.2

0.4

Defasagem

FAC

(a)

5 10 15 20 25 30

0.0

0.2

0.4

Defasagem

FAC

(b)

5 10 15 20 25 30

0.0

0.2

0.4

Defasagem

FAC

(c)

5 10 15 20 25 300.0

0.2

0.4

Defasagem

FAC

(d)

Figura 2 – Funções de autocorrelação do processo NGINAR(1)s em que somente osvalores do período sazonal são alterados conforme (a) s = 3, (b) s = 4, (c)s = 7 e (d) s = 12.

3.2.3 Função de Transição

Definição 3.2.1. Sejam X e Z duas variáveis aleatórias discretas independentes, em quePX+Z(·), PX(·), PZ(·) são as funções de probabilidade das variáveis aleatórias X +Z, X eZ, respectivamente, a convolução da função de probabilidade da soma de duas variáveisaleatórias discretas independentes é definida por

PX+Z(X + Z = i) =∞∑

k=−∞PX(X = k)PZ(Z = i− k).

Seja {Yt}t∈Z o processo definido em (3.1), a função de transição do processo éobtida através da função de probabilidade da variável aleatória Yt | Yt−s, em que é possívelobservar que esta pode ser escrita como a soma de duas variáveis aleatórias independentesda forma

Yt | Yt−s = α ∗ Yt−s | Yt−s︸︷︷︸Binomial negativa

+ εt︸︷︷︸Mistura de geométricas

.

Assim, podemos encontrar a função de probabilidade de Yt | Yt−s a partir de umaconvolução dada por

P (Yt = j | Yt−s = i) =∑k∈Cj

P (V = k)P (εt = j − k),


em que Cj = {k : k ∈ D(V ) e (j−k) ∈ D(εt)} = {k : 0 ≤ k ≤ j}, com V variável aleatóriacomo definida em (2.1) e D(X) representa o domínio da variável aleatória X.

Proposição 3.2.1. Seja {Yt}t∈Z um processo NGINAR(1)s, a função de transição s-passos à frente de Yt é dada por

P (Yt = j|Yt−s = 0) =(

1− αµ

µ− α

)µj

(1 + µ)j+1 +(

αµ

µ− α

)αj

(1 + α)j+1 ,

(3.5)

P (Yt = j|Yt−s = i) = µαj+1

(µ− α)(1 + α)j+i+1

i+ j

j

+

(1− αµ

µ− α

)µj

(1 + α)i(1 + µ)j+1

j∑k=0

i+ k − 1k

[α(1 + µ)µ(1 + α)

]k,

em que i,j = 0,1,..., e µ 6= α.

Demonstração. Seja pi,j = P (Yt = j|Yt−s = i) e µ 6= α, se i = 0, pela definição do operadorthinning binomial negativo, temos que

p0,j = P (εt = j).

Se i 6= 0, então

pi,j =j∑

k=0

i+ k − 1k

αk

(α + 1)i+k

[(1− αµ

µ− α

)µj−k

(µ+ 1)j−k+1 +(

αµ

µ− α

)αj−k

(α + 1)j−k+1

]

= αj+1µ

(µ− α)(α + 1)j+i+1

j∑k=0

i+ k − 1k

+

(1− αµ

µ− α

)µj

(µ+ 1)j+1(α + 1)ij∑

k=0

i+ k − 1k

(α(µ+ 1)µ(α + 1)

)k

= αj+1µ

(µ− α)(α + 1)j+i+1

i+ j

j

+

(1− αµ

µ− α

)µj

(µ+ 1)j+1(α + 1)ij∑

k=0

i+ k − 1k

(α(µ+ 1)µ(α + 1)

)k

25

4 Métodos de estimação para os parâmetros

Na prática, os parâmetros do modelo NGINAR(1)s não são conhecidos, e paraajustar o modelo a um conjunto de dados se faz necessário estimá-los. Neste capítulo, serãodefinidos os estimadores para os parâmetros do processo NGINAR(1)s. Tais estimadoressão amplamente utilizados como método de estimação para os parâmetros dos processosINAR na literatura, sendo eles: os estimadores de máxima verossimilhança condicional;os estimadores de mínimos quadrados condicionais e os estimadores de Yule-Walker.

4.1 Estimador de máxima verossimilhança condicionalSeja Y1, ..., Yn uma realização do processo NGINAR(1)s e seja P (Yt = i|Yt−s = j)

a função de probabilidade definida em (3.5), com θ = (α, µ)> ∈ Θ = {(α, µ) : (α, µ) ∈[0, µ/(1 + µ))× (0,∞)} ⊂ R2, em que Θ é chamado de espaço paramétrico.

O método de máxima verossimilhança usual que, geralmente, é aplicado à variáveisaleatórias i.i.d., não pode ser empregado diretamente ao processo NGINAR(1)s por razãoda estrutura de dependência entre as variáveis aleatórias do processo. Em razão disto,considera-se a função de probabilidade conjunta de Ys+1, ..., Yn condicionada às variáveisaleatórias Y1, ..., Ys. A função de verossimilhança condicional de θ é definida por

L(θ) =n∏

t=s+1P (Yt = i | Yt−s = j),

em que i é o valor amostral observado da variável aleatória Yt no tempo t, e j é um valorconhecido da variável aleatória Yt−s no tempo t− s.

Por razão da função logaritmo ser uma função monótona, maximizar L(θ) em θ éequivalente a maximizar o logaritmo natural da função de verossimilhança

log[L(θ)] = l(θ) =n∑

t=s+1log [P (Yt = i | Yt−s = j)] ,

que na grande maioria das vezes se torna mais atrativo devido a facilidade do processo deotimização da função l(θ) em relação a função L(θ). Definimos como estimativa de máximaverossimilhança condicional, denotada por θmvc = (αmvc, µmvc)>, o vetor θmvc ∈ Θ tal quel(θ) ≤ l(θmvc) ∀θ ∈ Θ. Assim, θmvc é tal que U

(θmvc

)= 0, em que

U (θ) = ∂l(θ)∂θ

=

∂l(θ)∂α

∂l(θ)∂µ

, (4.1)

Capítulo 4. Métodos de estimação para os parâmetros 26

e a matriz de segundas derivadas de l(θ) é negativa definida.

Para o processo NGINAR(1)s, o sistema de equações em (4.1) não possui soluçãoanalítica, sendo assim, as estimativas de máxima verossimilhança condicional (MVC)devem ser obtidas via métodos numéricos.

Para estimar os parâmetros do modelo NGINAR(1)s numericamente, através dosestimadores MVC, foi utilizado o método de Nelder-Mead, que é um método de otimização,proposto por Nelder e Mead (1965). Uma das vantagens do método de Nelder-Mead sobreoutros métodos de otimização, é que não é necessário o uso de derivadas nem da matrizde segundas derivadas, mas somente da função que se deseja minimizar. Observe quemaximizar a função l(θ) é equivalente a minimizar a função h(θ) = −l(θ) em relação a θ,ou seja, o método de Nelder-Mead deve ser implementado para minimizar a função h(θ).

4.2 Estimador de mínimos quadrados condicionaisO estimador de mínimos quadrados condicionais (MQC), proposto por Klimko e

Nelson (1978), é um procedimento de estimação para processos estocásticos baseado naminimização da soma dos quadrados dos desvios das observações em relação a esperançacondicional do processo. O estimador MQC é frequentemente abordado como método deestimação para os processos INAR (AL-OSH; ALZAID, 1987), (RISTIĆ; BAKOUCH;NASTIĆ, 2009), (BOURGUIGNON et al., 2016).

Seja Y1, ..., Yn uma amostra do processo NGINAR(1)s, em que a distribuição deprobabilidade de Yi, i ∈ {1,2, ..., n}, depende do vetor de parâmetros θ = (α, µ)>. Definag(θ) = E[Yt | Y1, ..., Yn], função de θ. Para o modelo NGINAR(1)s, a função g(θ) éexpressa por

g(θ) = E[Yt | Yt−s = yt−s] = E [α ∗ Yt−s + εt | Yt−s = yt−s] = αyt−s + µ(1− α).

Definimos como estimador de MQC, denotado por θmqc, o vetor θmqc = (αmqc, µmqc)>

que minimiza a funçãoQ(θ) =

n∑t=s+1

[Yt − g(θ)]2 ,

ou seja, θmvc é tal que∂Q (θ)∂θ

∣∣∣∣∣θ=θmqc

= 0, (4.2)

e a matriz de segundas derivadas de Q(θ) é positiva definida.

Proposição 4.2.1. Dada a amostra Y1, ..., Yn, os estimadores MQC, θmvc, que resolvem


o sistema de equações (4.2) são dados por

αmqc =(n− s)

n∑t=s+1

YtYt−s −(

n∑t=s+1

Yt

)(n∑

t=s+1Yt−s

)

(n− s)n∑

t=s+1Y 2t−s −

(n∑

t=s+1Yt−s

)2 ,

e

µmqc =

n∑t=s+1

Yt − αmqcn∑

t=s+1Yt−s

(n− s)(1− αmqc).

Demonstração. Ver no apêndice A.

As propriedades assintóticas do estimador MQC foram primeiramente demons-tradas por Klimko e Nelson (1978) e abordadas sob o contexto de estimação em modelosde séries temporais não lineares por Tjøstheim (1986). Conforme demonstrado em Ris-tić, Bakouch e Nastić (2009), os estimadores de MQC para os parâmetros do processoNGINAR(1) possuem distribuição assintótica normal.

Como os processos NGINAR(1) e NGINAR(1)s possuem as mesmas estruturasmarginais e probabilísticas, segue que a distribuição assintótica do estimador de MQC paraos parâmetros do processo NGINAR(1)s é a mesma verificada em Ristić, Bakouch e Nastić(2009). Portanto, satisfeitas as condições necessárias, pelo Teorema 3.2 em Tjøstheim(1986), temos que

√n

αmqc

µmqc

− α

µ

a∼ N

0

0

,

(1+α)(µ+µ2+α+αµ−αµ2)µ(1+µ)

α(1+α)1−α

α(1+α)1−α

µ(1+µ)(1+α)1−α

, (4.3)

em que a notação Xna∼ N(a, b) significa que, quando n aumenta indefinidamente, a

variável aleatória Xn possui distribuição normal com parâmetros a e b.

4.3 Estimador de Yule-WalkerA estimação dos parâmetros pelo método dos momentos geralmente é empregada

pela fácil obtenção dos estimadores pontuais. Frequentemente, as estimativas pelo métododos momentos são utilizadas como chute inicial quando se tem a necessidade de estimar osparâmetros através da implementação de um algoritmo de otimização. O método consisteem equiparar os momentos populacionais aos momentos amostrais, e então encontrar umafunção da amostra em que seja possível estimar o parâmetro desejado.

Considere uma amostra Y1, ..., Yn do processo NGINAR(1)s. Com base em (3.4),temos que ρ(k) = αk/s. Utilizamos a função de autocorrelação amostral, para k = s, como


estimador de Yule-Walker (YW) para o parâmetro α. Então, o estimador de Yule-Walker,αyw, para o parâmetro α é dado por

αyw = ρ(s) =

n∑t=s+1

(Yt − Y )(Yt−s − Y )n∑t=1

(Yt − Y )2,

em que Y = (1/n)∑ni=1 Yi. Com base no primeiro momento, o estimador de YW para

E [Yt] = µ, denotado por µyw, é dado por

µyw = Y .

Em Freeland e McCabe (2005) foi estudado o comportamento dos estimadoresde YW e MQC para os parâmetros do processo INAR(1) com marginais Poisson, efoi demonstrado que a distribuição do estimador MQC é assintoticamente equivalente adistribuição dos estimadores de YW. Tal propriedade também foi verificada para o processoNGINAR(1) em Ristić, Bakouch e Nastić (2009), o que implica que os estimadores de YWpara os parâmetros do processo NGINAR(1)s possuem distribuição assintótica equivalentea distribuição assintótica do estimador de MQC expressa em (4.3).

29

5 Estudo de Simulação

Neste capítulo, é apresentado um estudo de simulação realizado com a finalidadede comparação entre os estimadores definidos no Capítulo 4. Para isso, realizaremos umestudo de Monte Carlo com 5000 réplicas, com tamanhos amostrais n = 100, 200, 400 e800 utilizando os seguintes parâmetros α = 0,3 e 0,5, µ = 5 e 10 com período sazonal s= 7 e 12. Para este estudo, utilizamos o software R (Core Team, 2017) versão 3.3.3 naplataforma windows.

As informações apresentadas nas tabelas da Seção 5.1 constam da estimativamédia obtida em cada réplica e o erro quadrático médio do estimador que pode ser estimadopor

EQM(θi) = 15000

5000∑i=1

(θi − θi

)2,

em que θi ∈ {α, µ} é o estimador de θi ∈ {α, µ}, i ∈ {1,2, ..., 5000}.

Consideramos como melhor estimador aquele que apresentar menor erro quadráticomédio. Tal critério leva em conta a variabilidade do estimador como também o seu viés,de acordo com a seguinte expressão

EQM(θ) = V ar[θ] +[E(θ − θ)

]2.

5.1 Análise dos cenáriosA Tabela 1 apresenta os resultados do estudo de simulação para o processo

NGINAR(1)7 com média µ = 5, tamanho amostral variando em 100, 200, 400 e 800, comα = 0,3 e 0,5. As estimativas estão próximas de seus respectivos parâmetros estimados,além disso, o erro quadrático médio diminui quando o tamanho da amostra aumenta.Quanto ao aumento do parâmetro α, os erros quadráticos médios de todos os estimadoressão maiores do que ou iguais quando aumentamos de α = 0,3 para α = 0,5. Acreditamosque tal comportamento é observado devido à relativa proximidade, do valor 0,5 escolhidopara α, à restrição, que neste cenário, é de 5/6 ≈ 0,83. Além disso, percebe-se que adiferença entre tais erros é maior para os estimadores de µ do que para os de α. Nota-setambém que na sua grande maioria, os estimadores subestimam o verdadeiro valor doparâmetro.

Quanto a análise do melhor estimador, para o parâmetro α, o estimador demáxima verossimilhança condicional obteve um menor erro quadrático médio em todasas combinações estudadas na Tabela 1, já para o parâmetro µ, tal comportamento só foi

Capítulo 5. Estudo de Simulação 30

verificado para tamanhos amostrais maiores que n = 100, em que para este, o estimadorde Yule-Walker obteve menor erro quadrático médio.

Tabela 1 – Estimativa média e erro quadrático médio (em parênteses) para os estimadoresdo processo NGINAR(1)7 com µ = 5.

Estimador YW Estimador MQC Estimador MVCn α αyw µyw αmqc µmqc αmvc µmvc100 0,3 0,2612 5,0001 0,2817 4,9940 0,2996 4,9831

(0,0105) (0,5236) (0,0106) (0,5898) (0,0044) (0,5485)0,5 0,4368 4,9775 0,4719 4,9931 0,4931 4,9711

(0,0132) (0,7975) (0,0108) (0,9772) (0,0050) (0,8892)

200 0,3 0,2800 4,9951 0,2904 4,9933 0,2996 4,9893(0,0054) (0,2605) (0,0054) (0,2757) (0,0022) (0,2572)

0,5 0,4675 5,0215 0,4854 5,0223 0,4967 5,0090(0,0058) (0,4317) (0,0050) (0,4655) (0,0022) (0,4252)

400 0,3 0,2904 4,9929 0,2957 4,9930 0,3000 4,9919(0,0027) (0,1323) (0,0027) (0,1365) (0,0010) (0,1277)

0,5 0,4832 4,9951 0,4921 4,9947 0,4977 4,9896(0,0028) (0,2220) (0,0026) (0,2323) (0,0012) (0,2178)

800 0,3 0,2939 4,9973 0,2965 4,9971 0,2999 4,9949(0,0013) (0,0711) (0,0013) (0,0722) (0,0005) (0,0671)

0,5 0,4929 4,9951 0,4972 4,9946 0,4994 4,9907(0,0013) (0,1133) (0,0013) (0,1160) (0,0006) (0,1063)

Com relação as combinações dos parâmetros utilizadas neste experimento, a Tabela2 apresenta como diferença da Tabela 1 um aumento na média do processo de 5 para10, que reflete no aumento da variabilidade e no aumento da magnitude dos valores doprocesso que levam a menos observações repetidas.

Conforme a Tabela 2, semelhante ao ocorrido com a Tabela 1, as estimativasestão próximas de seus respectivos parâmetros estimados, com o erro quadrático médiodiminuindo quando o tamanho da amostra aumenta. Comparando as estimativas obtidasem relação ao aumento de α, o erro quadrático médio de todos os estimadores de µ émaior quando α = 0,5 do que para α = 0,3, mas desta vez, o estimador de mínimosquadrados condicionais de α teve seu erro quadrático médio reduzido. Em sua grandemaioria, como observado na Tabela 1, os estimadores tendem a subestimar o verdadeirovalor do parâmetro.

Nota-se como principal diferença entre a Tabela 1 e a Tabela 2, a magnitudedo erro quadrático médio dos estimadores de µ. Enquanto que para as estimativas de αpraticamente não houve diferenças, o erro quadrático médio dos estimadores de µ sãotodos maiores do que 3 quando n = 100 e chega a ser maior do que 4 para o estimador de


mínimos quadrados condicionais.

Novamente, para o parâmetro α, o estimador de máxima verossimilhança condicio-nal obteve um menor erro quadrático médio em todas as combinações estudadas na Tabela2, para o parâmetro µ, o estimador de Yule-Walker apresentou menor erro quadráticomédio quando n = 100, e a partir de n = 200, o estimador de máxima verossimilhançacondicional passa a ser o melhor estimador de µ.



(0,0102) (1,9669) (0,0101) (2,1764) (0,0024) (1,9723)0,5 0,4414 10,0412 0,4765 10,0641 0,4986 10,0105

(0,0117) (3,0939) (0,0094) (4,8880) (0,0026) (3,2367)

200 0,3 0,2792 9,9983 0,2896 9,9985 0,2994 9,9961(0,0049) (1,0040) (0,0048) (1,0767) (0,0011) (0,9738)

0,5 0,4694 10,0114 0,4866 10,0125 0,4990 9,9983(0,0053) (1,5557) (0,0045) (1,6977) (0,0012) (1,4795)

400 0,3 0,2913 10,0063 0,2965 10,0065 0,3003 10,0057(0,0025) (0,4955) (0,0025) (0,5115) (0,0005) (0,4682)

0,5 0,4855 10,0435 0,4943 10,0462 0,5001 10,0324(0,0025) (0,7970) (0,0023) (0,8340) (0,0006) (0,7456)

800 0,3 0,2958 9,9976 0,2985 9,9984 0,3003 10,0012(0,0012) (0,2533) (0,0012) (0,2576) (0,0003) (0,2327)

0,5 0,4926 10,0097 0,4970 10,0078 0,4997 10,0080(0,0011) (0,3937) (0,0011) (0,4046) (0,0003) (0,3614)

Ao comparar a Tabela 3 com a Tabela 1, é possível verificar o comportamento dosestimadores em amostras com períodos sazonais diferentes. Análogo as tabelas anteriores,as estimativas estão próximas de seus respectivos parâmetros estimados, além disso, oerro quadrático médio de todos os estimadores diminuem quando aumentamos o tama-nho da amostra. Quando o parâmetro α aumenta de 0,3 para 0,5, percebe-se o mesmocomportamento presente nas análises anteriores, com o aumento do erro quadrático médio.

Com respeito ao aumento do período sazonal, não verificamos alterações signifi-cantes no erro quadrático médio dos estimadores, poderíamos usar como exceção quandon = 100, em que nota-se um aumento do erro quadrático médio para os estimadores de µ.

Quanto ao melhor estimador, mais uma vez verificou-se que para o parâmetro α, oestimador de máxima verossimilhança condicional obteve um menor erro quadrático médioem todas as combinações estudadas na Tabela 3, já para o parâmetro µ, o estimador de


Yule-Walker se mostrou mais eficiente para n = 100, enquanto que a partir de n = 200, oestimador de máxima verossimilhança condicional obteve menor erro.



(0,0118) (0,5072) (0,0114) (0,6340) (0,0046) (0,6021)0,5 0,4155 5,0010 0,4762 5,0244 0,4957 5,0091

(0,0165) (0,7725) (0,0115) (1,1272) (0,0054) (1,0781)

200 0,3 0,2732 4,9994 0,2906 4,9982 0,3003 4,9967(0,0056) (0,2653) (0,0054) (0,2971) (0,0022) (0,2816)

0,5 0,4587 5,0118 0,4887 5,0156 0,4985 5,0035(0,0065) (0,4193) (0,0051) (0,4928) (0,0023) (0,4540)

400 0,3 0,2858 5,0046 0,2947 5,0059 0,2996 5,0022(0,0027) (0,1364) (0,0027) (0,1440) (0,0010) (0,1348)

0,5 0,4771 4,9949 0,4921 4,9946 0,4978 4,9881(0,0030) (0,2123) (0,0026) (0,2265) (0,0011) (0,2116)

800 0,3 0,2927 5,0004 0,2972 5,0007 0,2996 4.9988(0,0013) (0,0678) (0,0013) (0,0693) (0,0005) (0,0647)

0,5 0,4889 4,9930 0,4963 4,9923 0,4992 4,9901(0,0014) (0,1150) (0,0013) (0,1195) (0,0006) (0,1093)

Em relação a Tabela 3, a Tabela 4 se representa um aumento da média do processoµ de 5 para 10, ocasionando um aumento da variabilidade e da magnitude dos valores dasérie, como já comentado anteriormente.

Conforme a Tabela 4, é possível observar a proximidade entre as estimativas médiascom seus respectivos parâmetros, além disso, o erro quadrático médio dos estimadoresdiminui conforme o tamanho da amostra aumenta. Com relação ao aumento do parâmetroα, fica mais evidente que os estimadores do parâmetro µ sofrem mais influência de talaumento, principalmente em amostras pequenas. Conforme a Tabela 4, assim como aTabela 2, observa-se que o erro quadrático médio dos estimadores do parâmetro µ é altoquando os valores de µ são altos.

Para o parâmetro α, o estimador de máxima verossimilhança condicional obteveum menor erro quadrático médio em todas as combinações estudadas na Tabela 4, jápara o parâmetro µ, o estimador de Yule-Walker apresentou menor erro quadrático médioquando n = 100 e n = 200, e a partir de n = 400, o estimador de máxima verossimilhançacondicional passa a ser o melhor estimador de µ.

De modo geral, percebe-se que as estimativas médias estão próximas de seus




(0,0112) (1,9184) (0,0108) (2,3992) (0,0026) (2,1816)0,5 0,4182 9,9752 0,4781 10,0074 0,4983 9,9714

(0,0149) (2,7999) (0,0100) (4,0417) (0,0027) (3,4048)

200 0,3 0,2736 10,0089 0,2913 10,0049 0,3007 9,9975(0,0053) (0,9660) (0,0052) (1,0762) (0,0012) (0,9789)

0,5 0,4571 9,9984 0,4878 10,0080 0,4985 9,9830(0,0062) (1,5077) (0,0047) (1,7186) (0,0012) (1,5171)

400 0,3 0,2869 10,0015 0,2958 10,0009 0,3000 10,0024(0,0026) (0,5053) (0,0025) (0,5313) (0,0005) (0,4864)

0,5 0,4783 10,0019 0,4933 10,0020 0,4994 9,9956(0,0027) (0,7816) (0,0023) (0,8484) (0,0006) (0,7440)

800 0,3 0,2926 10,0144 0,2971 10,0165 0,3001 10,0115(0,0013) (0,2578) (0,0013) (0,2627) (0,0003) (0,2399)

0,5 0,4893 10,0053 0,4968 10,0056 0,4998 10,0050(0,0012) (0,4018) (0,0011) (0,4199) (0,0003) (0,3791)

respectivos valores dos parâmetros, e pode-se dizer que, na maioria das vezes, os estimadoressubestimam o verdadeiro valor do parâmetro. Quando o tamanho da amostra aumenta,é possível notar que o erro quadrático médio, de todos os estimadores, diminui, umindicativo de que tais estimadores são assintoticamente consistentes, isto é, quando naumenta indefinidamente, os estimadores ficam mais próximos do verdadeiro valor doparâmetro com variância tendendo para 0.

É notável o aumento do erro quadrático médio dos estimadores quando o valorde µ é grande, denotando uma diminuição na precisão das estimativas. Além disso, suasestimativas também são afetadas quando o parâmetro α está próximo à restrição. Nãofoi verificado diferenças significantes nas estimativas com relação ao aumento do períodosazonal.

Ainda que o estimador de máxima verossimilhança condicional exija um trabalhocomputacional considerável, acreditamos que este obteve um desempenho satisfatório, emrelação aos outros estimadores abordados neste capítulo, na estimação do parâmetro α.Sugerimos utilizar o estimador de máxima verossimilhança condicional para estimar α, poistal estimador se mostrou mais preciso do que os outros em qualquer cenário, enquanto quesugerimos utilizá-lo para estimar µ apenas com n > 100. Quando n ≤ 100, recomendamosutilizar o estimador de Yule-Walker para estimar o parâmetro µ.

34

6 Aplicação

A aplicabilidade do modelo teórico a dados reais se faz necessária para justificar apesquisa. Os dados utilizados neste trabalho se referem a um recorte mensal, no períodode janeiro de 1985 a dezembro de 1994, de 120 observações do número de reivindicações debenefícios de invalidez a curto prazo, em que todos os requerentes são do sexo masculino comidade entre 35 e 54 anos, trabalham na indústria madeireira e relataram sua reivindicaçãoao Conselho de Compensação de Trabalhadores Work Safe BC em Richmond, Canadá.Somente reclamantes cujos ferimentos foram devidos a cortes e lacerações foram incluídosno conjunto de dados.

Conforme Work Safe BC (2017), a Work Safe BC é uma empresa com enfoqueem segurança do trabalho, que, possui como principais objetivos promover a prevenção delesões e doenças no local de trabalho, reabilitar aqueles que são feridos e providenciar oretorno ao trabalho, fornecer uma compensação justa para substituir a perda de saláriosdos trabalhadores ao se recuperar de lesões e garantir uma boa gestão financeira para umsistema viável de compensação dos trabalhadores.

As observações fazem parte do conjunto de dados analisados previamente porFreeland (1998), em que o modelo INAR(1) com marginal Poisson é ajustado aos dadoscom a finalidade de modelar e realizar previsões a cerca do número de reivindicaçõesdos trabalhadores. O principal objetivo deste capítulo é ajustar o modelo NGINAR(1)sao conjunto de dados e comparar tal ajuste com o dos modelos INAR(1), INAR(1)s eNGINAR(1).

6.1 Análise descritivaA Tabela 5 apresenta algumas medidas descritivas. É possível notar que a série

possui valores de baixa magnitude, em que 75% das observações são de até 4 reivindicaçõesmensais, o valor máximo de 15 reivindicações em um mês destaca-se como ponto discrepante.Note que a série possui sobredispersão, com índice de dispersão, razão entre a variância ea média, de 1,86, i.e., a variância é 86% maior do que a média.

A Figura 3 apresenta o gráfico de barras e boxplot das observações. Neste primeirográfico, é possível observar uma leve assimetria, com coeficiente de assimetria de 8,05,em que a maioria das observações se concentram entre nenhuma e seis reivindicaçõesmensais. Já no boxplot, é possível observar a presença do ponto discrepante já mencionadoanteriormente. Tal observação é referente ao número de reivindicações no mês de julho de1987.

Capítulo 6. Aplicação 35

Tabela 5 – Medidas descritivas.Medidas Valores AmostraisMínimo 01o Quartil 1,75Mediana 3Média 3,243o Quartil 4Máximo 15Variância 6,02

0 1 2 3 4 5 6 7 8 9

Nº de Reivindicações

Fre

quên

cia

05

1015

20

05

1015

Nº

de R

eivi

ndic

açõe

s

Figura 3 – Gráfico de barras e boxplot das observações.

A Figura 4 apresenta o gráfico da série, bem como a função de autocorrelação efunção de autocorrelação parcial. Ao observar o gráfico da série, nota-se que não existemtendências de crescimento, ou decrescimento, e que a evolução da série se dá em torno deuma média constante ao longo do tempo. Além disso, a menos do ponto discrepante, avariação das observações permanecem constantes. Por tais razões, as evidências sugeremque a série é estacionária. No gráfico da função de autocorrelação é possível observar umcomportamento senoidal, porém, sem o decaimento exponencial característico de processosautorregressivos. Entretanto, tal disposição pode ser fruto de uma componente sazonal,em que, ponderamos que a FAC evidencia que o fenômeno se repete de forma acentuadade 12 em 12 meses. A função de autocorrelação parcial sugere um modelo de segundaordem, entretanto, neste trabalho, nos restringiremos aos modelos de primeira ordem.


Mês

Nº

de R

eivi

ndic

açõe

s

0 20 40 60 80 100 120

05

1015

5 10 15 20 25 30

-0.2

0.2

0.6

1.0

Defasagem

FA

C

5 10 15 20 25 30

-0.1

0.1

DefasagemF

AC

P

Figura 4 – Gráfico da série, função de autocorrelação e função de autocorrelação parcialamostral.

6.2 Ajustes dos modelosAjustamos o modelo proposto neste trabalho, bem como os modelos NGINAR(1)

(RISTIĆ; BAKOUCH; NASTIĆ, 2009), INAR(1)12 (BOURGUIGNON et al., 2016) eINAR(1) (AL-OSH; ALZAID, 1987) e (MCKENZIE, 1985) à série temporal obtida.Utilizamos os estimadores recomendados pelos autores dos modelos ajustados ao conjuntode dados para a estimação dos parâmetros. A Tabela 6 apresenta as estimativas dosparâmetros do modelo NGINAR(1)12. Note que as estimativas para o parâmetro α

encontram-se dentro do espaço paramétrico α ∈ (0, µ/(1 + µ)).

Tabela 6 – Estimativas para os parâmetros.

EstimadoresParâmetro MVC MQC YWα 0,56 0,32 0,29µ 2,72 3,12 3,24

Conforme sugerido no Capítulo 5, utilizaremos apenas o estimador de MVC paraestimar os parâmetros do modelo, e como chutes iniciais, as estimativas de YW. Sendoassim, estima-se que a média do processo é de 2,72 reivindicações mensais. Estima-setambém, que a autocorrelação do número mensal de reivindicações com defasagem 12 sejade ρ(12) = 0,56, i.e., o número de reivindicações se comporta de maneira semelhante paramesmos meses ao longo dos anos.

Semelhantemente ao processo NGINAR(1)s, os modelos INAR(1), INAR(1)s


e NGINAR(1) possuem dois parâmetros, que também chamaremos neste trabalho de αe µ, em que α representa a autocorrelação em defasagem 1 para os modelos INAR(1) eNGINAR(1), e o parâmetro µ representa a média do processo NGINAR(1), porém, paraos modelos INAR(1) e INAR(1)s, representa a média do processo de inovação destes. ATabela 7 apresenta as estimativas de máxima verossimilhança condicional dos parâmetrosdos modelos.

Tabela 7 – Estimativas para os parâmetros.

EstimativasModelo α µNGINAR(1) 0,51 2,88INAR(1) 0,19 2,64INAR(1)12 0,22 2,45

Finalmente, iremos comparar o ajuste dos modelos discutidos nesta seção. Paraisto, utilizamos como critérios de comparação o critério de informação de Akaike (AKAIKE,1970), conhecido como AIC, e o critério de informação Bayesiano, ou de Schwarz (SCHWARZet al., 1978), conhecido como BIC. Estes são os critérios que são comumente utilizadospara comparação entre o ajuste de modelos a um conjunto de dados.

Nas expressões do AIC e BIC, são levadas em consideração a qualidade de ajuste,em termos de verossimilhança, e a penalização por inclusão de parâmetros. Sendo assim,entre os modelos utilizados, o que possuir menor AIC e BIC é considerado como tendo omelhor ajuste aos dados com a menor quantidade de parâmetros.

Conforme a Tabela 8, para ambos os critérios considerados, o modeloNGINAR(1)12

obteve um melhor ajuste, em termos de AIC e BIC, portanto, para o conjunto de dadosdo número de reivindicações mensais estudado, o modelo NGINAR(1)12 é o modelo maisindicado dentre os abordados neste capítulo.

Tabela 8 – Critérios de comparação para os modelos NGINAR(1)12, NGINAR(1),INAR(1) e INAR(1)12.

Critérios de InformaçãoModelos AIC BICNGINAR(1)12 482,51 497,67NGINAR(1) 540,41 555,56INAR(1) 536,79 551,94INAR(1)12 487,47 502,62

Chamamos atenção para o fato de que os modelos NGINAR(1)12 e INAR(1)12,que consideram a estrutura sazonal dos dados, obtiveram melhor ajuste, evidenciando queos modelos sazonais para valores inteiros melhoram a qualidade de ajuste. Acreditamosque a diferença entre os ajustes dos modelos NGINAR(1)12 e INAR(1)12 se deu pela


presença de sobredispersão nos dados e que o modelo NGINAR(1)12 obteve um melhorajuste por considerá-la entre suas suposições.

Como método de diagnóstico, utilizamos os resíduos definidos por

rt = yt − E[Yt|Yt−s], t = 13,..., 120,

A Figura 5 apresenta a análise dos resíduos para o ajuste do modeloNGINAR(1)12

aos dados. Conforme a função de autocorrelação dos resíduos, a menos de uma defasagem,existem evidências de que os resíduos se comportam de forma independente, que é umadas suposições fundamentais do processo.

5 10 15 20 25 30

-0.2

0.2

0.6

1.0

Defasagem

FAC

Figura 5 – Análise de resíduos para o ajuste do modelo NGINAR(1)12.

39

Considerações Finais

Neste trabalho, consideramos um novo modelo para processos de contagem queapresentem estrutura sazonal e que possuem sobredispersão. Com relação ao desenvolvi-mento teórico, estudamos o operador thinning binomial negativo, definindo-o e apresen-tando suas principais propriedades como descritas em Ristić, Bakouch e Nastić (2009).Propomos métodos de estimação, baseados nos principais estimadores que são encontradosna literatura para os parâmetros dos modelos de séries temporais de valores inteiros, ealém disso, avaliamos a performance desses estimadores em diferentes cenários por meiode um estudo de simulação.

Concluímos que, para tamanhos amostrais maiores que 100, o estimador de máximaverossimilhança condicional é o melhor estimador, em termos de viés e erro quadráticomédio, dentre os apresentados nesta pesquisa, para estimar os parâmetros do modeloNGINAR(1)s. Embora sua estrutura seja complicada e envolva um certo conhecimentocomputacional, acreditamos que seu desempenho se mostrou relevante ao ponto de serescolhido.

Com relação a aplicabilidade do modelo proposto, dados sazonais estão ligadosa repetições regulares dos fenômenos de interesse e são comuns em diversas áreas. Nestetrabalho, analisamos um conjunto de dados referente ao número de reivindicações debenefícios de invalidez a curto prazo, no período de janeiro de 1985 a dezembro de 1994,e percebemos, segundo os critérios de informação de Akaike e bayesiano, que o modeloproposto possui um melhor ajuste do que outros modelos propostos anteriormente naliteratura.

Como principais considerações acerca deste trabalho, conforme já discutidos noCapítulo 6, os modelos que consideraram a estrutura sazonal dos dados obtiveram ummelhor ajuste do que os modelos que não consideraram tal estrutura, algo que era esperado.Além disso, o modelo NGINAR(1)s obteve melhor ajuste do que o modelo INAR(1)s,em que supomos que tal resultado se deva a sobredispersão nos dados observados.

40

Referências

AKAIKE, H. Statistical predictor identification. Annals of the Institute of StatisticalMathematics, Springer, v. 22, n. 1, p. 203–217, 1970. 37

AL-OSH, M.; ALZAID, A. A. First-order integer-valued autoregressive INAR(1) process.Journal of Time Series Analysis, Wiley Online Library, v. 8, n. 3, p. 261–275, 1987. 10,13, 26, 36

AL-OSH, M. A.; ALY, E.-E. A. First order autoregressive time series with negativebinomial and geometric marginals. Communications in Statistics-Theory and Methods,Taylor & Francis, v. 21, n. 9, p. 2483–2492, 1992. 11, 13

ANTUNES, J. L. F.; CARDOSO, M. R. A. Uso da análise de séries temporais emestudos epidemiológicos. Epidemiologia e Serviços de Saúde, Coordenação-Geralde Desenvolvimento da Epidemiologia em Serviços/Secretaria de Vigilância emSaúde/Ministério da Saúde, v. 24, n. 3, p. 565–576, 2015. 9

AWALE, M. Modeling seasonality in epidemic surveillance data using count time seriesmodels. In: 2nd ISI Regional Statistics Conference. [S.l.: s.n.], 2017. 11

BOURGUIGNON, M.; VASCONCELLOS, K. L.; REISEN, V. A.; ISPÁNY, M. A poissonINAR(1) process with a seasonal structure. Journal of Statistical Computation andSimulation, Taylor and Francis, v. 86, n. 2, p. 373–387, 2016. 11, 26, 36

BOURGUIGNON, M.; WEIß, C. H. An INAR(1) process for modeling count time serieswith equidispersion, underdispersion and overdispersion. TEST, Springer, p. 1–22, 2017.10, 11

BOX, G. E.; JENKINS, G. M.; REINSEL, G. Forecasting and control. Time SeriesAnalysis, v. 3, p. 75, 1970. 9, 17

BOX, G. E.; JENKINS, G. M.; REINSEL, G. C. Time series analysis: forecasting andcontrol. New Jersey: Prentice Hall, 1994. 9, 17

CARDINAL, M.; ROY, R.; LAMBERT, J. On the application of integer-valued timeseries models for the analysis of disease incidence. Statistics in Medicine, Wiley OnlineLibrary, v. 18, n. 15, p. 2025–2039, 1999. 9

Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria,2017. Disponível em: <https://www.R-project.org/>. 29

FRANSES, P. H. et al. Periodicity and stochastic trends in economic time series. OUPCatalogue, Oxford University Press, 1996. 10

FREELAND, R. K. Statistical analysis of discrete time series with application to theanalysis of workers’ compensation claims data. Tese (Doutorado) — University of BritishColumbia, 1998. 34

FREELAND, R. K.; MCCABE, B. Asymptotic properties of cls estimators in the PoissonAR(1) model. Statistics & probability letters, Elsevier, v. 73, n. 2, p. 147–153, 2005. 28

https://www.R-project.org/

Referências 41

GAUTHIER, G.; LATOUR, A. Convergence forte des estimateurs des paramètresd’un processus GENAR(p). In: UNIVERSITÉ DU QUÉBEC À MONTRÉAL,DÉPARTEMENT DE MATHÉMATIQUES ET INFORMATIQUE. Annales des SciencesMathématiques du Québec. [S.l.], 1994. v. 18, n. 1, p. 49–71. 13

GRASSLY, N. C.; FRASER, C. Seasonal infectious disease epidemiology. Proceedings ofthe Royal Society of London B: Biological Sciences, The Royal Society, v. 273, n. 1600, p.2541–2550, 2006. 10

JACOBS, P. A.; LEWIS, P. A. Stationary discrete autoregressive-moving average timeseries generated by mixtures. Journal of Time Series Analysis, Wiley Online Library, v. 4,n. 1, p. 19–36, 1977. 10

JACOBS, P. A.; LEWIS, P. A. Discrete time series generated by mixtures II: asymptoticproperties. Journal of the Royal Statistical Society. Series B (Methodological), JSTOR, p.222–228, 1978. 10

KHOO, W. C.; ONG, S. H.; BISWAS, A. Modeling time series of counts with a new classof inar (1) model. Statistical Papers, Springer, p. 1–24, 2017. 10

KIM, J.-M.; JUN, S. Integer-valued garch processes for apple technology analysis.Industrial Management & Data Systems, Emerald Publishing Limited, n. just-accepted, p.00–00, 2017. 10

KLIMKO, L. A.; NELSON, P. I. On conditional least squares estimation for stochasticprocesses. The Annals of statistics, JSTOR, p. 629–642, 1978. 26, 27

LATOUR, A. Existence and stochastic structure of a non-negative integer-valuedautoregressive process. Journal of Time Series Analysis, Wiley Online Library, v. 19, n. 4,p. 439–455, 1998. 18

LUTERBACHER, J.; DIETRICH, D.; XOPLAKI, E.; GROSJEAN, M.; WANNER,H. European seasonal and annual temperature variability, trends, and extremes since1500. Science, American Association for the Advancement of Science, v. 303, n. 5663, p.1499–1503, 2004. 10

MACAULAY, F. R. et al. Some theoretical problems suggested by the movements ofinterest rates, bond yields and stock prices in the united states since 1856. NBER Books,National Bureau of Economic Research, Inc, 1938. 10

MCKENZIE, E. Some simple models for discrete variate time series. JAWRA Journal ofthe American Water Resources Association, Wiley Online Library, v. 21, n. 4, p. 645–650,1985. 10, 11, 13, 36

MCKENZIE, E. Autoregressive moving-average processes with negative-binomial andgeometric marginal distributions. Advances in Applied Probability, Cambridge UniversityPress, v. 18, n. 3, p. 679–705, 1986. 10, 11, 13

MITCHELL, W. C. et al. Business cycles: The problem and its setting. NBER Books,National Bureau of Economic Research, Inc, 1927. 10

NELDER, J. A.; MEAD, R. A simplex method for function minimization. The computerjournal, Oxford University Press, v. 7, n. 4, p. 308–313, 1965. 26

Referências 42

PEREIRA, M. B. Modelos inar sazonais e de raízes unitárias. Dissertação (Mestrado) —Universidade Federal de Pernambuco, Recife - PE, 2012. 10

QUDDUS, M. A. Time series count data models: an empirical application to trafficaccidents. Accident Analysis & Prevention, Elsevier, v. 40, n. 5, p. 1732–1741, 2008. 9

RISTIĆ, M. M.; BAKOUCH, H. S.; NASTIĆ, A. S. A new geometric first-orderinteger-valued autoregressive (NGINAR(1)) process. Journal of Statistical Planning andInference, Elsevier, v. 139, n. 7, p. 2218–2226, 2009. 11, 12, 13, 15, 18, 26, 27, 28, 36, 39

SCHWARZ, G. et al. Estimating the dimension of a model. The annals of statistics,Institute of Mathematical Statistics, v. 6, n. 2, p. 461–464, 1978. 37

STEUTEL, F.; HARN, K. V. Discrete analogues of self-decomposability and stability. TheAnnals of Probability, JSTOR, p. 893–899, 1979. 13

TJØSTHEIM, D. Estimation in nonlinear time series models. Stochastic Processes andtheir Applications, Elsevier, v. 21, n. 2, p. 251–273, 1986. 27

WEIß, C. H. Modelling time series of counts with overdispersion. Statistical Methods &Applications, Springer, v. 18, n. 4, p. 507–519, 2008. 11

WEIß, C. H. The INARCH(1) model for overdispersed time series of counts.Communications in Statistics-Simulation and Computation, Taylor & Francis, v. 39, n. 6,p. 1269–1291, 2010. 11

WEIß, C. H.; KIM, H.-Y. Parameter estimation for the binomial AR(1) models withapplications in finance and industry. Statistical Papers, Springer, v. 54, n. 3, p. 563–590,2013. 9

Work Safe BC. Work Safe BC: Our mandate, vision, mission, goals & values. 2017.<https://www.worksafebc.com/en/about-us/who-we-are/mission-vision-values>. Accessoem: 21 de novembro de 2017. 34

ZHU, F. Modeling overdispersed or underdispersed count data with generalized Poissoninteger-valued GARCH models. Journal of Mathematical Analysis and Applications,Elsevier, v. 389, n. 1, p. 58–71, 2012. 11

https://www.worksafebc.com/en/about-us/who-we-are/mission-vision-values

43

APÊNDICE A – Estimador de mínimosquadrados condicionais

Demonstraremos neste apêndice que os estimadores de MQC são pontos críticos

da função Q(θ), ou seja, θmqc é tal que ∂Q(θ)∂θ

∣∣∣∣∣θ=θmqc

= 0. As derivadas parciais da

função Q(θ) são dadas por

∂Q(θ)∂µ

=− 2n∑

t=s+1[Yt − αYt−s − µ(1− α)] (1− α),

∂Q(θ)∂α

=− 2n∑

t=s+1[Yt − αYt−s − µ(1− α)] (Yt−s − µ) .

Sejam µmqc e αmqc tais que

∂Q(θ)∂µ

∣∣∣∣∣µ=µmqc

= 0 e ∂Q(θ)∂α

∣∣∣∣∣α=αmqc

= 0, (A.1)

temos que, pela primeira das equações em (A.1)n∑

t=s+1Yt − αmqc

n∑t=s+1

Yt−s − (n− s)µmqc(1− αmqc) = 0

⇒ µmqc =

n∑t=s+1

Yt − αmqcn∑

t=s+1Yt−s

(n− s)(1− αmqc). (A.2)

Por outro lado,n∑

t=s+1YtYt−s − αmqc

n∑t=s+1

Y 2t−s − µmqc

n∑t=s+1

Yt−s(1− αmqc) = 0, (A.3)

já que,

− 2µmqcn∑

t=s+1[Yt − αmqcYt−s − µmqc(1− αmqc)]

=− 2µmqc[

n∑t=s+1

Yt − αmqcn∑

t=s+1Yt−s − (n− s)(1− αmqc)µmqc

]= 0.

Substituindo µmqc em (A.3), e realizando as devidas manipulações algébricas,temos que

αmqc =(n− s)

n∑t=s+1

YtYt−s −(

n∑t=s+1

Yt

)(n∑

t=s+1Yt−s

)

(n− s)n∑

t=s+1Y 2t−s −

(n∑

t=s+1Yt−s

)2 .

44

APÊNDICE B – Rotinas computacionais

####################################################################### Processo NGIGINAR(1)s #######################################################################

# Geração de valores aleatórios de EtEt = function(n,alpha,mu){

e = rep(NA,n)a = (alpha*mu)/(mu-alpha)u = runif(n,0,1)n1=length(e[as.numeric(u<a)==1])n2=length(e[as.numeric(u<a)==0])e[as.numeric(u<a)==1] = rgeom(n1, 1/(1 + alpha))e[as.numeric(u<a)==0] = rgeom(n2, 1/(1 + mu))return(e)

}

# Função de probabilidade de Etpet = function(l,alpha,mu){

a = (alpha*mu)/(mu-alpha)return( (1-a)*dgeom(l,1/(1+mu)) + a*dgeom(l,1/(1+alpha)) )

}

# Geração de valores aleatórios da sérieNGINAR1s = function(n,alpha,mu,s){y = vector();n = n + 150; y[1:s] = rgeom(s,1/(1+mu)) ;e=Et(n,alpha,mu)for (t in (s+1):n){

if(y[t-s]==0){y[t] = e[t]

}else{y[t] = rnbinom(1,y[t-s],1/(1+alpha)) + e[t]

}}

return(ts(y[151:n]))

APÊNDICE B. Rotinas computacionais 45

}

####################################################################### Estimador de Yule Walker #######################################################################YW = function(y,s){

return(c(acf(y,plot=F)$acf[s+1],mean(y)))}

####################################################################### Estimador de MQC #######################################################################MQC = function(y,s){

n = length(y);yt = y[(s+1):n];yt_s = y[1:(n-s)]a = ((n-s)*sum(as.numeric(yt*yt_s)) -

(sum(as.numeric(yt))*sum(as.numeric(yt_s))))/((n-s)*sum(as.numeric(yt_s*yt_s)) - (sum(as.numeric(yt_s))*sum(as.numeric(yt_s))))

m = (sum(yt) - a*sum(yt_s))/((n-s)*(1-a))return(c(a,m))

}

####################################################################### Estimador de Máxima Verossimilhança Condicional #######################################################################mvc = function(y,s){

n = length(y)llik<-function(theta){

alpha <- theta[1];mu <- theta[2];p = vector(); p[1:s] = 0

for(t in (s+1):n){

if(y[t-s]==0){p[t] = log(pet(y[t],alpha,mu))

}else{

APÊNDICE B. Rotinas computacionais 46

i=0:y[t]Soma = sum( choose(y[t-s]+i-1,i)*( ( (alpha*(mu+1))/(mu*(alpha+1)) )^i ) )

p[t] = log( mu*alpha^(y[t]+1)/((mu - alpha)*(1 + alpha)^(y[t]+y[t-s]+1))*choose(y[t]+y[t-s],y[t]) +

(1 - mu*alpha/(mu - alpha))*mu^y[t]/(((1+alpha)^y[t-s])*((1+mu)^(y[t]+1)))*Soma)

}

}

-sum(p)}

optim(c(YW(y,s)[1]+0.1,YW(y,s)[2]+0.5),llik,method="Nelder-Mead",hessian = T)

}

## Resíduosres = function(x,alpha,mu,s){

E = alpha*lag(x,-s) + (1 - alpha)*mue = x-Ereturn(e)

}

Documents

ProcessoINAR(1)comEstruturaSazonalpara ... · Medeiros, Rodrigo Matheus Rocha de. Processo INAR(1) com estrutura sazonal para séries temporais de valores inteiros com sobredispersão