Recentemente, tem havido um interesse crescente em modelos de regressão para séries temporais de contagens e um número bastante considerável de publicações sobre este assunto tem surgido na literatura. As séries temporais de contagem aparecem naturalmente em várias áreas sempre que um número de eventos por período de tempo é observado ao longo do tempo. Exemplos que mostram a ampla gama de aplicações são o número diário de internações hospitalares do sistema de saúde pública, o número de transações do mercado de ações por minuto em bolsa de valores por minuto ou o número horário de itens defeituosos do controle de qualidade industrial.

Modelos para séries temporais de contagem devem levar em conta que as observações são números inteiros não negativos e devem capturar adequadamente a dependência entre as observações. Uma abordagem conveniente e flexível é empregar a metodologia do modelo linear generalizado (GLM) (Nelder and Wedderburn, 1972) para modelar as observações condicionalmente às informações passadas, escolhendo uma distribuição adequada para dados de contagem e uma função de ligação apropriada. Essa abordagem é seguida em detalhes por Fahrmeir and Tutz (2001) e Kedem and Fokianos (2002), entre outros.

Outra classe importante de modelos para séries temporais de contagens é baseada no operador de diluição, como os modelos de média móvel autorregressiva inteira (INARMA), que, de certa forma, imitam a estrutura dos modelos comuns de média móvel autorregressiva (ARMA), para um revisão recente ver Wei (2008). Outro tipo de modelos de séries temporais de contagem são os chamados modelos de espaço de estados. Referimo-nos às revisões de Fokianos (2011), Jung and Tremayne (2011), Fokianos (2012), Tjøstheim (2012) e Fokianos (2015) para uma visão geral aprofundada dos modelos para séries temporais de contagem.

Os métodos mencionados aqui estão disponíveis na biblioteca de funções R tscount. Esta implementação, semelhante a sintaxe da função R glm, permite uma estrutura de dependência mais geral que pode ser especificada convenientemente pelo usuário. Assim, temos modelos de séries temporais que incluem um processo latente, semelhante à classe de modelos GARCH. O uso e a saída de nossas funções são inspirados nas funções do R glm e arima.

Modelos

Seja uma série temporal de contagem denotada por \(\{Y_t : t\in \mathbb{N}\}\). Vamos denotar por \(\{X_t:t\in\mathbb{N}\}\) um vetor de covariáveis \(r\)-dimensional variável no tempo, digamos \(X_t=\big(X_{t,1},\cdots,X_{t,r}\big)^\top\). Modelamos a esperança condicional \(\mbox{E}\big(Y_t|\mathcal{F}_{t-1}\big)\) da série temporal de contagem por um processo de média latente, digamos \(\{\lambda_t:t\in\mathbb{N}\}\), tal que \(\mbox{E}\big(Y_t|\mathcal{F}_{t-1}\big)=\lambda_t\).

Denotamos por \(\mathcal{F}_{t-1}\) a história do processo conjunto \(Y_t,\lambda_t,X_{t+1}:t\in\mathbb{N}\) até o tempo \(t\) incluindo a informação da covariável no tempo \(t+1\). A hipótese de distribuição para \(Y_t\) dado \(\mathcal{F}_{t-1}\) é discutida mais adiante. Estamos interessados em modelos da forma geral $$ g(\lambda_t)=\beta_0+\sum_{k=1}^p \beta_k \widetilde{g}(Y_{t-i_k})+\sum_{l=1}^q \alpha_l g(\lambda_{t-j_l})+\eta^\top X_t, $$ onde \(g:\mathbb{R}^+\to \mathbb{R}\) é uma função de ligação e \(\widetilde{g}:\mathbb{R}^+\to \mathbb{R}\) é uma função de transformação.

O vetor de parâmetros \(\eta=\big(\eta_1,\cdots,\eta_r\big)^\top\) corresponde aos efeitos das covariáveis. Na terminologia dos modelos lineares generalizados chamamos \(\nu_t=g(\lambda_t)\) de preditor linear. Para permitir a regressão em observações passadas arbitrárias da resposta, definimos o conjunto \(\mathbb{P}=\{i_1,i_2,\cdots,i_p\}\) de inteiros \(0< i_1< i_2< \cdots < i_p < \infty\) com \(p\in \mathbb{N}_0\). Isso nos permite regredir nas observações defasadas \(Y_{t-i_1},Y_{t-i_2},\cdots,Y_{t-i_p}\). Analogamente, definimos o conjunto \(\mathbb{Q}=\{j_1,j_2,\cdots,j_q\}\) com \(q\in \mathbb{N}\) e inteiros \(0< j_1< j_2< \cdots < j_q< \infty\) para regressão nas médias latentes defasados \(\lambda_{t-j_1},\lambda_{t-j_2},\cdots,\lambda_{t-j_q}\). Este caso mais geral é coberto pela teoria para modelos com \(\mathbb{P}=\{1,2,\cdots,p\}\) e \(\mathbb{Q}=\{1,2,\cdots,q\}\), que geralmente são tratados na literatura, escolhendo \(p\) e \(q\) suficientemes grandes e definindo parâmetros desnecessários para o modelo como zero.

Vejamos vários exemplos do modelo acima. Considere a situação em que \(g\) e \(\widetilde{g}\) são iguais à identidade, ou seja, \(g(x)=\widetilde{g}(x)=x\). Além disso, seja \(\mathbb{P}=\{1,2,\cdots,p\}\), \(\mathbb{Q}=\{1,2,\cdots,q\}\) e \(\eta=0\). Então obtemos que $$ \lambda_t = \beta_0+\sum_{k=1}^p \beta_k Y_{t-k}+\sum_{l=1}^q \alpha_l \lambda_{t-l}\cdot $$

Assumindo ainda que \(Y_t\) dado o passado é distribuído segundo a distribuição Poisson, então obtemos um modelo GARCH de valor inteiro de ordem \(p\) e \(q\), em resumo INGARCH\((p,q)\). Esses modelos foram discutidos por Heinen (2003), Ferland, Latour and Oraichi (2006) e Fokianos, Rahbek and Tjøstheim (2009), entre outros. Quando \(\eta\neq 0\), então o pacote tscount ajusta modelos INGARCH com covariáveis não-negativas; isso ocorre porque precisamos garantir que o processo médio resultante seja positivo.

Considere novamente o modelo geral, mas agora com a função de ligação logarítmico \(g(x)=\log(x)\) e \(\widetilde{g}(x)=\log(x+1)\) e \(\mathbb{P}\) e \(\mathbb{Q}\) como antes. Obtemos então um modelo log-linear de ordem \(p\) e \(q\) para a análise de séries temporais de contagem. De fato, defina \(\nu_t=\log(\lambda_t)\) para obter da expressão geral que $$ \nu_t=\beta_0 + \sum_{k=1}^p \beta_k \log\big( Y_{t-k}+1\big) +\sum_{l=1}^q \alpha_t \nu_{t-l}\cdot $$

Este modelo log-linear é estudado por Fokianos and Tjøstheim (2011), Woodard, Matteson and Henderson (2011) e Douc, Doukhan and Moulines (2013). Seguimos Fokianos and Tjøstheim (2011) na transformação das observações passadas empregando a função \(\widetilde{g}=\log(x+1)\), de modo que estejam na mesma escala que o preditor linear \(\nu_t\) ver Fokianos and e Tjøstheim (2011) para um discussão e por mostrar que a adição de uma constante a cada observação para evitar zeros não afeta a inferência. Observe que o modelo acima permite modelagem de correlação serial negativa, enquanto anterior acomoda apenas correlação serial positiva. Adicionalmente, o modelo acima acomoda covariáveis mais facilmente do que o modelo anteior, uma vez que o modelo log-linear implica positividade do processo médio condicional \(\{\lambda_t\}\).

Os efeitos das covariáveis na resposta são multiplicativos para o modelo acima; é aditivo para o modelo anterior a este. Para uma discussão sobre a inclusão de covariáveis dependentes do tempo, veja Fokianos and Tjøstheim (2011).

O modelo geral juntamente com a suposição de Poisson, ou seja, \(Y_t|\mathcal{F}_{t-1}\sim \mbox{Poisson}(\lambda_t)\), implica que $$ P\big(Y_t=y|\mathcal{F}_{t-1} \big)=\dfrac{\lambda_t^y \exp(-\lambda_t)}{y!}, \qquad y=0,1,\cdots\cdot $$

Obviamente, \(\mbox{Var}\big(Y_t=y|\mathcal{F}_{t-1} \big)=\mbox{E}\big(Y_t=y|\mathcal{F}_{t-1} \big)=\lambda_t\). Portanto, no caso de um modelo de resposta de Poisson condicional, o processo médio latente é idêntico à variância condicional do processo observado.

A distribuição Binomial Negativa permite uma variância condicional maior que \(\lambda_t\). Seguindo Christou and Fokianos (2014), assume-se que \(Y_t|\mathcal{F}_{t-1}\sim \mbox{Binomial Negativa}(\lambda_t,\phi)\), onde a distribuição Binomial Negativa é parametrizada em função de sua esperança com um parâmetro de dispersão adicional \(\phi\in (0,\infty)\), ou seja, $$ P\big(Y_t=y|\mathcal{F}_{t-1} \big)=\dfrac{\Gamma(\phi+y)}{\Gamma(y+1)\Gamma(\phi)}\left(\dfrac{\phi}{\phi+\lambda_t}\right)^\phi\left(\dfrac{\lambda_t}{\phi+\lambda_t}\right)^y, \qquad y=0,1,\cdots\cdot $$

Neste caso, \(\mbox{Var}\big(Y_t=y|\mathcal{F}_{t-1} \big)=\lambda_t+\lambda_t^2/\phi\), ou seja, a variância condicional aumenta quadraticamente com \(\lambda_t\). A distribuição de Poisson é um caso limite da Binomial Negativo quando \(\phi\to \infty\). Note que a distribuição Binomial Negativa pertence à classe dos processos mistos de Poisson. Um processo misto de Poisson é especificado definindo \(Y_t\sim N_t(0,Z_t\lambda_t)\), onde \(\{N_t\}\) representa variáveis aleatórias igualmente distribuídas (i.i.d). Os processos de Poisson com intensidade unitária e \(\{Z_t\}\) são variáveis aleatórias i.i.d. com esperança 1 e variância \(\sigma^2\). Quando \(\{Z_t\}\) é um processo i.i.d. de variáveis aleatórias Gamma, então obtemos o processo Binomial Negativo com \(\sigma^2=1/\phi\). Referimo-nos a \(\sigma^2\) como o coeficiente de superdispersão porque é proporcional à extensão da superdispersão da distribuição condicional. O caso limite de \(\sigma^2=0\) corresponde à distribuição de Poisson, ou seja, sem superdispersão. O procedimento de estimação não se limita ao caso Binomial Negativo, mas a qualquer distribuição mista de Poisson. No entanto, a suposição Binomial Negativa é necessária para intervalos de previsão e avaliação do modelo.

No modelo geral o efeito de uma covariável entra plenamente na dinâmica do processo e se propaga para observações futuras tanto pela regressão em observações passadas quanto pela regressão em médias latentes passadas. O efeito de tais covariáveis pode ser visto como uma influência interna no processo de geração de dados, razão pela qual nos referimos a ele como um efeito covariável interno. Também é permitido incluir covariáveis de forma que seu efeito apenas se propague para observações futuras pela regressão em observações passadas, mas não diretamente pela regressão em médias latentes passadas.

Seguindo Liboschik, Kerschke, Fokianos and Fried (2014), que fazem essa distinção para o caso de efeitos de intervenção descritos por covariáveis determinísticas, nos referimos ao efeito de tais covariáveis como um efeito covariável externo. Seja \(e=(e_1,\cdots,e_r)^\top\) um vetor especificado pelo usuário com \(e_i=1\) se o \(i\)-ésimo componente do vetor covariável tiver um efeito externo e \(e_i=0\) caso contrário, \(i=1,\cdots,r\). Denotando por \(\mbox{diag}(e)\) uma matriz diagonal com elementos diagonais dados por \(e\).

A generalização, permitindo efeitos covariáveis internos e externos então seria $$ g(\lambda_t)=\beta_0+\sum_{k=1}^p \beta_k \widetilde{g}(Y_{t-i_k})+\sum_{l=1}^q \alpha_l\Big( g(\lambda_{t-j_l})-\eta^\top\mbox{diag}(e) X_{t-j_l} \Big)+\eta^\top X_t\cdot $$

Basicamente, o efeito de todas as covariáveis com efeito externo é subtraído nos termos de feedback de tal forma que seu efeito entra na dinâmica do processo apenas através das observações. Referimo-nos a Liboschik et al. (2014) para uma extensa discussão de efeitos internos versus externos. É nossa experiência com esses modelos que uma discriminação empírica entre efeitos covariáveis internos e externos é difícil e que não é crucial qual tipo de efeito covariável se encaixa em aplicações práticas.

Estimação e inferência

Estimam-se os parâmetros dos modelos por máxima verossimilhança quase condicional (ML), função tsglm. Se a suposição de Poisson for verdadeira, então obtemos uma estimador de máxima verossimilhança ordinário. No entanto, sob a suposição mista de Poisson, obtemos um estimador de quase-máxima verossimilhnaça.

Vamos denotar o vetor de parâmetros da regressão como \(\theta=\big(\beta_0,\beta_1,\cdots,\beta_p,\alpha_1,\cdots,\alpha_q,\eta_1,\cdots,\eta_r\big)^\top\). Independentemente da suposição de distribuição, o espaço paramétrico para o modelo INGARCH com covariáveis é dado por \begin{equation*} \Theta=\left\{ \theta\in\mathbb{R}^{p+q+r+1}: \beta_0>0, \beta_1,\cdots,\beta_p,\alpha_1,\cdots,\alpha_q,\eta_1,\cdots,\eta_r\geq 0, \sum_{k=1}^p \beta_k+\sum_{l=1}^q \alpha_l<1\right\}\cdot \end{equation*}

O intercepto \(\beta_0\) deve ser positivo e todos os outros parâmetros devem ser não negativos para garantir a positividade do processo de média latente. A outra condição garante que o modelo ajustado tenha uma solução estacionária (Ferland et al. 2006). Para o modelo log-linear com covariáveis, o espaço paramétrico é considerado \begin{equation*} \Theta=\left\{ \theta\in\mathbb{R}^{p+q+r+1}: \beta_0>0, |\beta_1|,\cdots,|\beta_p|,|\alpha_1|,\cdots,|\alpha_q|<1, \left|\sum_{k=1}^p \beta_k+\sum_{l=1}^q \alpha_l\right| <1\right\}\cdot \end{equation*}

Esta pretende ser a generalização para a ordem do modelo \(p\), \(q\) das condições \(\beta_1|<1\), \(\alpha_1<1\) e \(\beta_1+\alpha_1<1\) que Douc et al. (2013) derivam para o modelo de primeira ordem. Christou and Fokianos (2014) apontam que com a parametrização da distribuição Binomial Negativa a estimação dos parâmetros de regressão \(\theta\) não depende do parâmetro de dispersão adicional \(\phi\). Isso permite empregar uma abordagem de quase máxima verossimilhança baseada na verossimilhança de Poisson para estimar os parâmetros de regressão \(\theta\), que é descrita abaixo. O parâmetro de incômodo ou perturbação \(\phi\) é então estimado separadamente em uma segunda etapa.

A função de log-verossimilhança, o vetor score e a matriz de informação são derivados condicionalmente nos valores pré-amostrais da série temporal e do processo médio latente \(\{\lambda_t\}\). Uma inicialização apropriada é necessária para sua avaliação, que é discutida na próxima subseção. Para um trecho de observações \(y = (y_1,\cdots,y_n)^\top\), a função de quase -log-verossimilhança condicional, até uma constante, é dada por \begin{equation*} \ell(\theta)=\sum_{t=1}^n \log\big(p_t(y_t;\theta)\big) \approx \sum_{t=1}^n \left(y_t\ln\Big(\lambda_t(\theta)\Big)-\lambda_t(\theta) \right), \end{equation*} onde \(p_t(y_t;\theta)=P(Y_t=y|\mathcal{F}_{t-1})\) é a função de probabilidadae da distribuição Poisson conforme foi apresentada.

O processo médio latente é considerado uma função \(\lambda_t : \Theta\to\mathbb{R}^+\) e, portanto, é denotado por \(\lambda_t(\theta)\) para todo \(t\). A função escore condicional é um vetor de dimensão \(p+q+r+1)\) dado por \begin{equation*} S_n(\theta)=\dfrac{\partial \ell(\theta)}{\partial\theta}=\sum_{t=1}^n \left( \dfrac{y_t}{\lambda_t(\theta)}-1\right)\dfrac{\partial \lambda_t(\theta)}{\partial \theta}\cdot \end{equation*}

O vetor de derivadas parciais \(\lambda_t(\theta)/\partial \theta\) pode ser calculado recursivamente. Finalmente, a matriz de informação condicional é dada por \begin{equation*} \mathcal{I}_n(\theta,\sigma^2)=\sum_{t=1}^n \mbox{Cov}\left(\dfrac{\partial\ell(\theta;y_t)}{\partial\theta} \Big| \mathcal{F}_{t-1}\right)=\sum_{t=1}^n \left( \dfrac{1}{\lambda_t(\theta)}+\sigma^2\right) \left( \dfrac{\partial\lambda_t(\theta)}{\partial\theta}\right) \left( \dfrac{\partial\lambda_t(\theta)}{\partial\theta}\right)^\top\cdot \end{equation*}

O estimador de quase máxima verossimilhança (QML) \(\widehat{\theta}_n\) é, supondo que exista, a solução do problema de otimização restrita não linear \begin{equation*} \widehat{\theta}_n=\mbox{arg}\max_{\theta\in\Theta} \ell(\theta)\cdot \end{equation*}

Conforme proposto por Christou and Fokianos (2014), o parâmetro de dispersão \(\phi\) da distribuição Binomial Negativa é estimado resolvendo a equação \begin{equation*} \sum_{t=1}^n\dfrac{\big(y_t-\widehat{\lambda}_t\big)^2}{\widehat{\lambda}_t(1+\widehat{\lambda}_t/\widehat{\phi})}=n-m, \end{equation*} que é baseado na estatística \(\chi^2\) de Pearson.

O parâmetro de variância \(\sigma^2\) é estimado por \(\widehat{\sigma}^2=1/\widehat{\phi}\). Para a distribuição Poisson definimos \(\widehat{\sigma}^2=0\). A rigor, o modelo log-linear não se enquadra na classe de modelos considerados por Christou and Fokianos (2014). No entanto, os resultados obtidos por Douc et al. (2013) permitem utilizar este estimador também para o modelo log-linear, especificamente no caso em que \(p=q=1\).

A inferência para os parâmetros de regressão é baseada na normalidade assintótica do estimador QML, que foi mostrada por Fokianos et al. (2009) e Christou and Fokianos (2014) para modelos sem covariáveis. Para um processo covariado bem comportado \(\{X_t\}\) temos que \begin{equation*} \sqrt{n}\big( \widehat{\theta}_n-\theta_0\big) \xrightarrow{D} N_{p+q+r+1}\Big(0,\mathcal{I}_n^{-1}(\widehat{\theta},\widehat{\sigma}^2) \mathcal{I}_n^*(\widehat{\theta})\mathcal{I}_n^{-1}(\widehat{\theta},\widehat{\sigma}^2)\Big), \end{equation*} quando \(n\to\infty\), onde \(\theta_0\) denota o verdadeiro valor do parâmetro e \(\widehat{\sigma}^2\) é um estimador consistente de \(\sigma^2\).

Suponhamos que isso se aplica sob as mesmas suposições geralmente feitas para o modelo de regressão linear ordinária (veja por exemplo, Demidenko 2013). Para covariáveis deterministas, essas suposições são \(|| X_t || < c\), isto é, o processo covariado é limitado e \(lim_{n\to\infty} \frac{1}{n}\sum_{t=1}^n X_t X_t^\top = A\), onde \(c\) é uma constante e \(A\) é uma matriz não singular. Para covariáveis estocásticos, presume-se que as esperanças \(\mbox{E}(X_t)\) e \(\mbox{E}(X_tX_t^\top)\) existam e que \(\mbox{E}(X_tX_t^\top)\) seja não-singular. As suposições implicam que as informações sobre cada covariável crescem linearmente com o tamanho da amostra e que as covariáveis não são linearmente dependentes. Fuller (1996) mostra a normalidade assintótica do estimador de mínimos quadrados para um modelo de regressão com erros de série temporal sob condições ainda mais gerais que permitem a presença de certos tipos de tendências nas covariáveis.

A normalidade assintótica do estimador QML é suportada por simulações. Uma prova formal requer mais pesquisas. Para evitar instabilidades numéricas ao inverter \(\mathcal{I}_n(\widehat{\theta},\widehat{\sigma}^2)\) aplicamos um algoritmo que faz uso do fato de que é uma verdadeira matriz simétrica e definida positiva.

Uma alternativa à aproximação normal acima para obter erros padrão é um procedimento de bootstrap paramétrico, que faz parte do pacote, função se. Por conseguinte, \(B\) séries temporais são simuladas a partir do modelo estimado para os dados originais. Os erros padrão empíricos das estimativas de parâmetros para essas \(B\) séries temporais são os erros padrão do Bootstrap. Este procedimento pode calcular os erros padrão não apenas para os parâmetros de regressão estimados, mas também para o coeficiente de dispersão \(\sigma^2\).

Implementação

As restrições dos parâmetros que são impostas pela condição \(\theta\in\Theta\) podem ser formuladas como \(d\) desigualdades lineares. Isso significa que existe uma matriz \(U\) de dimensão \(d\times (p+q+r+1)\) e um vetor \(c\) de comprimento \(d\), tal que \(\Theta=\{\theta \, |\, U\theta\geq c\}\). Para o modelo linear necessitamos \(d=p+q+r+2\) restrições para garantir a nonegatividade da média latente \(\lambda_t\) e estacionariedade do processo resultante. Para o modelo log-linear, não há restrições no intercepto nem nos coeficientes das covariáveis e o número total de restrições é \(d=2(p+q+1)\). Para reforçar as desigualdades estritas as respectivas restrições são apertadas por uma constante arbitrariamente pequena \(\xi>0\); esta constante é definida como \(\xi=10^{-6}\) por padrão, argumento slackvar.

Para resolver numericamente o problema de maximização emprega-se a função constrOptim. Esta função aplica um algoritmo descrito por Lange (1999), o que essencialmente aplica as restrições, adicionando um valor de barreira à função objetivo e, em seguida, emprega um algoritmo para a otimização irrestrita desta nova função objetivo, iterando esses dois passos, se necessário. Por padrão, o algoritmo quase-Newton Broyden-Fletcher-Goldfarb-Shanno (BFGS) é empregado para a última tarefa de otimização irrestrita, que adicionalmente faz uso do vetor escore.

Observe que as funções de log-verossimilhança e escore são condicionais em valores pré-amostrais não observados. Elas dependem do preditor linear e suas derivadas parciais, que podem ser computados recursivamente usando qualquer inicialização. Para visualizarmos recursões e várias estratégias para sua inicialização dispomos dos argumentos init.Method e init.drop. Christou and Fokianos (2014) mostram que o efeito da inicialização desaparece assintoticamente. No entanto, do ponto de vista prático, a inicialização das recursões é crucial. Especialmente na presença de forte dependência serial, as estimativas resultantes podem diferir substancialmente mesmo por longas séries de tempo com 1000 observações.

Resolver o problema de otimização não linear requer um valor inicial para o vetor de parâmetros \(\theta\). Este valor inicial pode ser obtido a partir de um modelo mais simples para o qual um procedimento de estimação está prontamente disponível. Consideramos ou o ajuste de modelo GLM ou o ajuste de um modelo autoregressivo de médias móveis (ARMA). Uma terceira possibilidade é ajustar um mdelo ingênuo (naive) i.i.d. sem covariates. Como última escolha, observe que poderíamos usar valores fixos que precisam ser fornecidos pelo estatístico. Todas as possibilidades estão disponíveis no pacote tscount, argumento start.control.

Acontece que o algoritmo de otimização converge de forma confiável, mesmo que os valores iniciais não sejam próximos do ótimo global da função de verossimilhança. Claro, um valor inicial mais próximo do máximo global geralmente requer menos iterações até a convergência. No entanto, foram encontrados alguns exemplos de dados em que os valores iniciais próximos a um ótimo local, obtido por um dos dois métodos de estimativa, podem até impedir encontrar o ótimo global. Consequentemente, recomendamos que ajustem o modelo ingênuo i.i.d. sem covariáveis para obter valores iniciais.

Exemplo

Série temporal com o número de casos de infecções por Campylobacter no norte da província de Quebec (Canadá) em intervalos de quatro semanas de janeiro de 1990 até o final de outubro de 2000. Possui 13 observações por ano e 140 observações no total. A campilobacterose é uma doença infecciosa bacteriana aguda que ataca o sistema digestivo. 


library(tscount)
data("campy")
par(mfrow=c(1,1), mar = c(4,5,3,1), pch=19)
plot(campy, type="b", xlab="", ylab="", main="Número de casos de infecções por Campylobacter")
grid()


campyfit <- tsglm(ts=campy, model=list(past_obs=1, past_mean=c(7,13)))
summary(campyfit)
## 
## Call:
## tsglm(ts = campy, model = list(past_obs = 1, past_mean = c(7, 
##     13)))
## 
## Coefficients:
##              Estimate  Std.Error  CI(lower)  CI(upper)
## (Intercept)    1.4965     0.7101     0.1048      2.888
## beta_1         0.5856     0.0527     0.4823      0.689
## alpha_7        0.0873     0.0640    -0.0381      0.213
## alpha_13       0.1838     0.0750     0.0369      0.331
## Standard errors and confidence intervals (level =  95 %) obtained
## by normal approximation.
## 
## Link function: identity 
## Distribution family: poisson 
## Number of coefficients: 4 
## Log-likelihood: -434.3842 
## AIC: 876.7685 
## BIC: 888.535 
## QIC: 877.9338


par(mfrow=c(1,1), mar = c(4,5,4,1), pch=19)
plot(campy, type="b", xlab="", ylab="", main="Número de casos de infecções por Campylobacter \n 
     observados vs estimativas (vermelho)")
points(fitted(campyfit), col="red")
grid()

Referências