2.1 Pressuposição 1: especificação

O modelo de regressão é linear nos parâmetros

A relação entre Y e X é linear

Esta pressuposição em princípio implica na consideração de uma reta estimada, ou seja, uma função linear nas variáveis do tipo

\[{Y_i} = {\beta _0} + {\beta _1}{X_{1i}} + {\beta _2}{X_{2i}} + \ldots + {\beta _k}{X_{ki}} + {\varepsilon _i}\]

ou pela forma matricial:

\[ Y = X\beta + \varepsilon , \]
em que \(Y\) é o vetor de variáveis explicadas, \(X\) é uma matriz de variáveis explicativas (incluindo uma coluna de uns para o intercepto) e \(\varepsilon\) é um vetor de resíduos aleatórios. Os resíduos são a parte inexplicada da regressão (das variações em \(Y\)), enquanto a parte de \(X\beta\) é a parte explicada das variações em \(Y\).
Entretanto, deve-se atentar para outros tipos de linearidades implícitas na pressuposição. Têm-se os seguintes tipos de linearidades: linearidade das variáveis explicativas (\(X\)) e linearidade dos parâmetros (\(\beta\)). A não linearidade nas variáveis às vezes pode ser contornada por transformações nas variáveis, mas a não linearidade dos parâmetros é mais complicada e requer outros métodos de estimação não lineares.
É fácil imaginar que o comportamento de um fenômeno econômico não siga a relação retilínea, como por exemplo, as tradicionais relações de oferta e demanda não necessariamente serão retas que se cruzam. É muito mais fácil admitir que o comportamento de variáveis econômicas seja curvilíneo. Na Figura 2.1, dos retornos das ações das Lojas Americanas em função de uma variável \(Z\) qualquer, observa-se que as possibilidades de ajustamentos em reta ou em parábola apresentam diferentes resultados em termos de melhor representar a nuvem de pontos.

Figure 2.1: Retornos das ações das Lojas Americanas em função de uma variável Z qualquer.

Fonte: Elaboração própria com dados do Yahoo Finance.

Quando as variáveis explicativas são elevadas a alguma potência diferente de um, a função que relaciona o comportamento dessas variáveis com a variável explicada será diferente de uma reta e os estimadores tradicionais de Mínimos Quadrados Ordinários (MQO) não mais serão válidos. O estimador de um parâmetro é uma “regra” ou “expressão” que resultará de uma amostra selecionada. Por exemplo, um estimador simples comumente utilizado é o da média aritmética de uma amostra.

Outro exemplo de modelo linear, agora nos parâmetros, refere-se ao que pode-se chamar de “intrinsecamente lineares”, ou que podem se tornar lineares por transformação das variáveis.

\[ Y = \alpha + \left( \beta_1 * X_1 \right) + \left( \beta_2 * X_2^2 \right) \]

\[{Y_i} = {\beta _0} + {\beta _1}{X_{1i}} + {\beta _2}{X_{2i}^2} + \ldots + {\beta _k}{X_{ki}} + {\varepsilon _i}\]

Neste caso, a variável \(X_{2i}\) está ao quadrado, e a equação é dita linear nos parâmetros mas não linear nas variáveis. Ainda é satisfatória nesse caso, se houver uma transformação da variável não-linear (\(X_{2i}\)).

O caso mais comum na literatura econômica é o de funções do tipo Cobb-Douglas, ou seja,

\[Y = AX_1^{\beta 1}X_2^{\beta 2}X_3^{\beta 3}{e^\varepsilon }\]

em que os parâmetros podem assumir valores diferentes de um e, ainda, tem-se a multiplicação de variáveis explicativas. A função acima pode ser linearizada transformando-se as variáveis em logaritmos, obtendo:

\[\ln Y = \ln A + {\beta _1}\ln {X_1} + {\beta _2}\ln {X_2} + {\beta _3}\ln {X_3} + \varepsilon \]

ou, simbolizando o \(ln\) por \(*\):

\[{Y^*} = {\beta _0} + {\beta _1}X_1^* + {\beta _2}X_2^* + {\beta _3}X_3^* + \varepsilon \]

A função linearizada pode ser estimada da forma tradicional lembrando que os parâmetros estimados serão agora da função transformada, que no caso log-log (Cobb-Douglas), equivalem às elasticidades. A função transformada pode ser vista como linear nos parâmetros (os parâmetros \(\beta\) são todos em primeira potência) e nas variáveis transformadas (\(X^*=\ln X\)).
Outros modelos não podem ser transformados e são os chamados intrinsecamente não lineares. Por exemplo, é possível perceber que a função abaixo não pode ser linearizada:

\[Y = A{\beta _1} + {\beta _2}{e^{{\beta _3}{X_1}}} + {\beta _4}{e^{{\beta _5}{X_2}}} + {e^\varepsilon }\]

Esses modelos devem ser estimados por Mínimos Quadrados não lineares ou Máxima Verossimilhança não linear. Algumas formas funcionais utilizadas em economia da produção podem ser:

• Cobb-Douglas logaritmizada: \[\log y = {a_0} + \sum\limits_{i = 1}^n {{a_i}\log {x_i}} \]
• Elasticidade Constante de Substituição ou CES: \[{y^\rho } = {a_0} + \sum\limits_{i = 1}^n {{a_i}x_i^\rho } \]
• Generalizada Leontief: \[y = {a_0} + \sum\limits_{i = 1}^n {{a_i}\sqrt {{x_i}} } + \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{a_{ij}}\sqrt {{x_i}{x_j}} } } \]
• Transcendental Logaritmica ou Translog: \[\log y = {a_0} + \sum\limits_{i = 1}^n {{a_i}\log {x_i}} + \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{a_{ij}}\log {x_i}\log {x_j}} } \]
  
• Quadrática: \[y = {a_0} + \sum\limits_{i = 1}^n {{a_i}{x_i}} + \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{a_{ij}}{x_i}{x_j}} } \]

A utilização de uma forma mais complexa em detrimento de uma mais simples dependerá da disposição dos dados e do rigor científico desejado. A função Cobb-Douglas de modo geral oferece um ajustamento satisfatório e é fácil de executar. As funções elasticidade de substituição constante (CES), Generalizada Leontief, Transcendental Logarítmica e Quadrática são generalizações da função Cobb-Douglas para contornar pressuposições econômicas de substitutibilidade dos fatores e produtos ou ainda de concorrência perfeita, entre outras situações.

Juntamente ao problema da forma funcional (linearidade dos parâmetros e variáveis), quando se especifica um modelo, automaticamente estão sendo cometidos outros dois tipos de erros que poderão ou não comprometer a análise. Um está associado à “omissão de uma variável relevante” e outro associado à “inclusão de variável irrelevante”.
- Omissão de variável relevante

Imagine que a revisão de literatura, revisão teórica, indique que a quantidade demandada (\(Q\)) de um produto seja função do preço do produto (\(P\)) e da renda (\(R\)), e que o comportamento da demanda do produto analisado na realidade está em conformidade com a teoria. O modelo “correto” seria o modelo A:

\[(A) \space Q_t = \beta_0 + \beta_1.P_t + \beta_2.R_t + \varepsilon_t^* \]

em que os \(\beta\) são parâmetros estimados e \(\epsilon\) é o resíduo aleatório.

Imagine agora que, por algum motivo, estimou-se a demanda em função apenas do preço do produto (modelo B), fazendo:

\[(B) \space Q_t = \alpha_0 + \alpha_1.P_t + \varepsilon_t. \]

em que \(\alpha\) são parâmetros e as demais variáveis como anteriormente citadas.

A questão é: quais as consequências sobre os estimadores de MQO (ou sobre os \(\beta\) estimados)? Qual o efeito sobre \(\alpha_0\) e \(\alpha_1\) em razão da exclusão(ou omissão) de \(R\) do modelo?

Se \(P_t\) for altamente correlacionado com \(R_t\), a retirada de \(R_t\) trará um alto viés (alta tendenciosidade) e os parâmetros estimados serão muito diferentes do valor esperado:

\[\beta_{estimado} ≠ E(\beta), \]

ou seja, os parâmetros estimados serão inconsistentes e no limite \(E(\beta) ≠ \beta\). Os testes de hipóteses não serão válidos e as estimativas de variâncias também serão tendenciosas.

• Inclusão de variável irrelevante

Imagine agora a situação inversa: o modelo estimado contempla mais variáveis explicativas do que as que deveriam estar no modelo “correto”. Imagine que o modelo deveria ter apenas \(P\) (modelo A) e que foi estimado com \(P\) e \(Z\) (modelo B), sendo \(Z\) uma variável irrelevante no modelo.

\[(A) \space Q_t = \beta_0 + \beta_1.P_t + \varepsilon_t, \] - modelo correto

\[(B) \space Q_t = \alpha_0 + \alpha_1.P_t + \alpha_2.Z_t + \varepsilon_t^*. \] - modelo estimado

e que Z não tem relevância teórica.

A questão é: quais as consequências de \(\alpha\), em razão da inclusão de \(Z_t\), sobre \(\beta\)?

As consequências da inclusão de uma variável irrelevante serão menos problemáticas que no caso da omissão de uma variável relevante. Primeiro, a presença das variáveis “irrelevantes” não viesa as outras estimativas. Segundo, aumentam-se a variância dos parâmetros e o desvio-padrão. Tende, portanto, a fazer com que o parâmetro “\(\alpha\)” seja não significativo, mas aumenta o coeficiente de determinação (a proporção explicada da variância total - medida usualmente chamada de \(R^2\)).

2.2 Pressuposição 2: média residual

O resíduo aleatório tem média zero

Significa que o resíduo (\(\varepsilon_i\)) tem uma distribuição de probabilidade centralizada em zero (com média zero). O resíduo é o efeito das variáveis que não consigo explicar no modelo. A média pode ser considerada como o valor esperado do erro, ou seja,

\(E(\varepsilon_i) = 0\), para \(i = 1,2,...,n\).

Ou na forma matricial,
\[E(\varepsilon) = 0. \] Dado que \(E(\epsilon) = 0\), então
\[E(Y) = E[ X\beta + \varepsilon] = E[X\beta] +E[\varepsilon] = E[X\beta] + 0 \]

Portanto, \(E(Y) = X\beta\) e o modelo fornece soluções adequadas estatisticamente. Essa pressuposição é importante para ter confiança na estimação por mínimos quadrados ordinários (\(MQO\)) em que (veremos mais adiante) o estimador dos parâmetros será \(\hat\beta_{MQO} = (X^\prime X)^{-1} X^\prime Y\). Caso os erros não tenham média zero, o estimador \(\hat\beta_{MQO}\) será tendencioso. Mais a frente serão discutidas as propriedades dos estimadores de \(MQO\).

Observe na Figura 2.1 (anteriormente mostrada) que traz as taxas de retorno observadas e estimadas para a ação das Lojas Americanas S.A. (LAME4), e verifique que existem momentos em que os pontos vermelhos (com marcador quadrado) estão acima que os verdes (com marcador de x) e em outros momentos estão abaixo. O gráfico dos resíduos (neste gráfico representados por \(u_i\)) obtidos fazendo resíduo igual a diferença entre o retorno observado e o estimado, têm-se valores positivos e negativos. A pressuposição prevê que estes, na média, sejam nulos. Ainda, no gráfico de dispersão dos retornos mensais dos ativos das Lojas Americanas (LAME4) e Lojas Renner (LREN3), Jan/2005 a Mar/2012 (aqui denominados RLAME x RREN, Figura 2.2), pode-se observar que existem resíduos \(u_i\) positivos e negativos e que a reta de regressão estimada como a reta de tendência passa aproximadamente no meio da nuvem de pontos.

Figure 2.2: Gráfico de dispersão dos retornos mensais dos ativos das Lojas Americanas (LAME4) e Lojas Renner (LREN3), Jan/2005 a Mar/2012.

Fonte: Elaboração própria com dados do Yahoo Finance.


Outra situação pode ser observada na Figura 2.3, com os resultados de uma estimação dos retornos mensais do ativo Lojas Americanas (LAME4) contra os do ativo Lojas Renner (LREN3) e os do índice da Bolsa de Valores e Mercadorias e Futuros de São Paulo (BMF-BOVESPA), Jan/2005 a Mar/2012. São mostrados os resíduos (linha azul), resultado da diferença entre os valores observados (linha vermelha) e os valores estimados (linha verde).

Figure 2.3: Representação dos resultados de uma estimação dos retornos mensais do ativo Lojas Americanas (LAME4) contra os do ativo Lojas Renner (LREN3) e os do índice da Bolsa de Valores e Mercadorias e Futuros de São Paulo (B3), Jan/2005 a Mar/2012.

Fonte: Elaboração própria com dados do Yahoo Finance.

Mais a frente traremos o capítulo específico onde faremos a checagem dessa pressuposição e como tratar o modelo em caso de violação do pressuposto.

2.3 Pressuposição 3: homocedasticidade

O resíduo aleatório tem variância constante (presença de homocedasticidade)

A variância é calculada com base no valor esperado do quadrado da diferença entre a média da variável e o valor esperado da média desta variável. Ou seja, a definição estatística para um resíduo \(\varepsilon_i\) é

\[V(\varepsilon_i) = E[\varepsilon_i – E(\varepsilon_i)]^2\],

\[V(\varepsilon_i) = E(\varepsilon_i^2) = \sigma^2\] (populacional) para todo \(i\),

ou seja, presença de Homocedasticidade nos resíduos. A presença da homocedasticidade implica que a variância para todos os resíduos é a mesma.

O caso contrário será:

\[V(e_i) = E(e_i^2) = \sigma^2_i \] - presença de Heterocedasticidade, ou uma variância para cada \(i\).

O problema de heterocedasticidade é típico de dados de seção cruzada. Pode significar, por exemplo, uma heterogeneidade da amostra. A dispersão dos valores para cada observação é diferente entre as observações. A amostra vem de uma população onde os erros não são homogêneos. Na Figura 2.4, observa-se que a nuvem de pontos não tem uma dispersão constante em torno da reta estimada, o que caracteriza a variabilidade distinta ao longo da amostra de \(X\). Com o aumento de \(X\), a variância dos resíduos (a dispersão em torno da reta estimada) reduz no valores menores e volta a aumentar para os vaores mais elevados de \(X\), caracterizando a presença de heterocedasticidade.

Figure 2.4: Representação da dispersão dos pontos em torno de uma reta estimada – ilustração da heterocedasticidade.

Fonte: Elaboração própria com dados aleatórios.


Mais a frente, em capítulo próprio, serão discutidas as formas de checar a presença de homocedasticidade e a correção no caso de violação desse pressuposto.

2.4 Pressuposição 4: autocorrelação

Os resíduos aleatórios são independentes (ou não autocorrelacionados)

Também chamada de ausência de correlação serial, neste caso pressupõe-se que os resíduos \(\varepsilon_i\) de uma observação \(i\) não afetam outros resíduos. Ou seja, matematicamente, isso equivale a

\[ COV\left( {{\varepsilon_i},{\varepsilon_j}} \right){\rm{ }} = {\rm{ }}E{\rm{ }}\left\{ {{\rm{ }}\left[ {{\varepsilon_i}-{\rm{ }}E\left( {{\varepsilon_i}} \right)} \right]{\rm{ }}\left[ {{\varepsilon_j}-{\rm{ }}E\left( {{\varepsilon_j}} \right)} \right]{\rm{ }}} \right\}\]

\[E\left( {{\varepsilon_i},{\rm{ }}{\varepsilon_j}} \right){\rm{ }} = {\rm{ }}0,\;\;{\rm{ }}i \ne j.\]

Esta pressuposição é denominada “ausência de autocorrelação”. A violação desta pressuposição é um problema típico de séries temporais. A correlação é, de modo simplificado, “a relação entre duas séries”. Pode ser positiva (diretamentte proporcional) ou negativa (inversa), ou mesmo nula. O valor da correlação varia então entre -1 e +1. Diz-se “auto” correlação porque estamos relacionando a série dos resíduos com ela mesma, para observações distintas. No caso de séries temporais, serão resíduos em um período \(t\) relacionados com resíduos em alguma defasagem temporal (valor passado) ou \(t-p\) em que \(p\) seria o número de períodos defasados. Ou seja, existe alguma dependência temporal entre o valor presente e valores passados, gerando um padrão inexplicado na variável explicada \(Y\) e que aparecerá nos resíduos.

Quando se trabalha com ajustamentos de séries temporais, essa pressuposição em geral não é obedecida, visto que nas séries temporais como, por exemplo, as séries de preços, de salários e de produção têm no seu comportamento o reflexo de movimentos cíclicos e/ou sazonais. Observa-se na Figura 2.5 que existem relações entre os resíduos das observações à medida que X aumenta – correlação negativa (caso a) e positiva (caso b).

Figure 2.5: Representações da dispersão dos pontos em torno de uma reta estimada – ilustração da correlação serial.

Fonte: Figura 3.3 de Pindyck e Rubinfeld (2004).

Algumas causas da autocorrelação nos resíduos estão relacionadas a variáveis não especificadas no modelo, forma funcional inadequada e inércia temporal no fenômeno.

A principal consequência da violação desta pressuposição é a ineficiência dos estimadores de \(MQO\), mas continuam não tendenciosos. Nesta situação, da mesma forma que para a heterocedasticidade, é melhor utilizar o método de Mínimos Quadrados Generalizados (\(MQG\)).

Uma forma usual é olhar os gráficos de dispersão entre os resíduos da regressão. Padrões geométricos podem indicar o tipo de correlação, como na Figura 2.6.

Figure 2.6: Padrões de correlação entre os distúrbios. (a) correlação serial positiva; (b) correlação serial negativa; e (c) correlação nula.

Fonte: Adaptado da Figura 3.6 de Gujarati (2006).

2.5 Pressuposição 5: exogeneidade

As variáveis explicativas são não aleatórias (são fixas)

Também pode ser expressa como as variáveis \(X\) e os resíduos são não correlacionados, ou ausência de correlação entre \(X\) e \(\varepsilon\). Neste caso, pressupõem-se fixos os valores da variável explicativa e observa-se o que ocorre com a variável dependente. Se o \(X\) é aleatório, mas independente do resíduo, pode-se mostrar que os parâmetros estimados serão não tendenciosos. Assim, a confirmação desta pressuposição significa dizer que as variáveis explicativas são distribuídas independentemente dos resíduos.

Em linguagem matemática, pode-se dizer que a covariância entre os resíduos \(\varepsilon_i\) e \(X_i\) é igual a zero (\(cov(\varepsilon_i,X_i)=0\)). Formalmente,

\(cov(\varepsilon_i,X_i)=E[\varepsilon_i-E(\varepsilon_i )][X_i-E(X_i )]\);

\(cov(\varepsilon_i,X_i)=E[\varepsilon_i (X_i-E(X_i ))]\) pois \(E(\varepsilon_i )=0\);

\(cov(\varepsilon_i,X_i)=E(\varepsilon_i X_i )-E(X_i )E(\varepsilon_i )\) pois \(E(X_i )\) é não estocástica;

\(cov(\varepsilon_i,X_i)=E(\varepsilon_i X_i )\) pois \(E(\varepsilon_i )=0\);

\(cov(\varepsilon_i,X_i)=0\) por pressuposição.

Entretanto, se as variáveis explicativas e os termos aleatórios forem correlacionados, haverá inconsistência dos estimadores de mínimos quadrados ordinários. Deve-se utilizar o estimador de variáveis instrumentais. O método de Variáveis instrumentais (VI) prevê que \(\beta_{VI} = (Z^\prime X)^{-1} Z^\prime Y\), e \(Z\) é uma matriz de instrumentos independentes dos resíduos aleatórios.

2.6 Pressuposição 6: normalidade

Os resíduos têm distribuição normal, com média zero e variância constante:

As vezes expressa simplesmente por os resíduos têm distribuição normal, usualmente se especifica essa pressuposição da forma:

\(\varepsilon_i ∼ N (0, \sigma^2)\), para toda observação \(i = 1, 2, ..., n\).

As consequências associadas à não normalidade dos resíduos são parâmetros estimados não normais e não será possível fazer os testes de hipóteses com distribuições baseadas na normal, como os usuais testes \(t\) e \(F\) para avaliar a qualidade dos ajustamentos, e para construir intervalos de confiança para os parâmetros conforme exposto ao longo do capítulo de regressão. Em termos gráficos, pode-se plotar o histograma da série de resíduos e comparar com a distribuição normal teórica como na Figura 2.7.

Figure 2.7: Histograma dos resíduos e a curva normal.

Fonte: Elaboração própria.

Os estimadores continuam sendo os Melhores Estimadores Lineares Não-Tendenciosos (\(MELNT\)), como será melhor exemplificado na seção das propriedades dos estimadores de \(MQO\).

2.7 Pressuposição 7: multicolinearidade

Ausência de multicolinearidade prejudicial entre as variáveis explicativas

Neste caso, pode-se interpretar como ausência de multicolinearidade perfeita assim como multicolinearidade expressiva entre as variáveis explicativas (\(X\)), a ponto de prejudicar a devida interpretação da significância dos parâmetros, como explicado a seguir.

A multicolinearidade é um problema relacionado com fortes relações entre as variáveis explicativas (\(X\)) no modelo de regressão. O que espera-se é que a matriz \(X\) explique \(Y\), mas não que \(X_k\) esteja relacionado com outro(s) \(X_l\) (\(k\ne l\)). Considere a matriz de variáveis explicativas como composta por colunas das variáveis \(X_1, X_2, ... , X_k\) e ainda uma coluna de uns (1) para incluir o intercepto.

A pressuposição prevê a inexistência de qualquer relação linear entre as variáveis explicativas \(X\), como por exemplo,

\[X_1 = 2.X_2 \\ X_1 + 3.X_2 = X_5\]

No método de mínimos quadrados ordinários, a existência de uma relação linear entre colunas de \(X\) representa uma redução no posto da matriz (\(X\)) e o determinante de \(X^\prime X\) será próximo de zero. No caso de uma relação linear exata, haverá uma singularidade perfeita na matriz (\(X^\prime X\)) e seu determinante será zero. Como o método de mínimos quadrados ordinários prevê a inversão da matriz \(X^\prime X\), o determinante próximo de zero fará com que os parâmetros sejam indeterminados. A matriz \((X^\prime X)^{-1}\) não existirá e não será possível estimar o modelo. No caso exato (multicolinearidade perfeita), o software interrompe e acusará erro. Em muitos casos aplicados, o que se observa são valores de determinantes muito próximo de zeros (por Ex.: 1e-5 = 0,00001), e muitas vezes o sistema não interrompe, requerendo a observação das estatísticas de teste.

O problema da correlação entre as variáveis explicativas pode ser visto da seguinte maneira:

1. ausência de correlação ou ausência de multicolinearidade: a regressão múltipla dá o mesmo resultado que as regressões simples quando as correlações parciais entre as variáveis explicativas forem nulas;

2. correlação perfeita ou multicolinearidade perfeita: a relação linear perfeita entre os X´s causa a indeterminação de \(\hat \beta = (X^\prime X)^{-1}X^\prime Y\) pois \((X^\prime X)^{-1}\) é singular; e,

3. alto grau de correlação entre as colunas de \(X\) ou multicolinearidade imperfeita: multicolinearidade que prejudica a significância dos parâmetros estimados..

O enfoque é diferente das outras pressuposições: trata-se de um problema da amostra. Enquanto as outras pressuposições se referiam mais ao resíduo e à população, esta se refere mais à amostra, ao número de variáveis relativo ao tamanho da amostra. Por isto, muitas vezes se fala do problema como micronumerosidade, ou que o número de observações é pequeno relativamente ao número de variáveis. Não se trata, portanto, de testar a pressuposição, mas sim de pensar como lidar com o problema.

2.8 Resumo das Pressuposições:

Apresentadas as pressuposições, o Quadro 2.1 tem um resumo com a expressão matemática em forma escalar e matricial, assim como o problema que se tem caso as pressuposições sejam violadas ou não atendidas. Em geral, pode-se dizer que se testará o modelo e, em caso de violação, se “tratará” ou corrigirá adequadamente.

Fonte: Elaboração própria. Notas:\(^*\) Em que \(Y = [Y_i]\) é um vetor (n x 1) das observações da variável dependente; \(X = [X_{ij}]\) é uma matriz (n x p) das observações das variáveis independentes; \(\varepsilon = [\varepsilon_i]\) é um vetor (nx1) dos resíduos aleatórios; \(\beta = [\beta_j]\), \(j = 0, 1, 2, ..., k\) é um vetor (p x 1) de parâmetros a serem estimados; \(\sigma^2\) é a variância do erro, também a ser estimada; \(I\) é uma matriz identidade de ordem (m x n); \(k\) é o número de variáveis independentes; \(p= (k + 1)\) é o número de parâmetros; \(n\) é o número de observações; \(E\) significa valor esperado ou esperança matemática.

\(^+\) PROBLEMA (o que acontece se as pressuposições não forem atendidas).

3.1 Introdução

A estimação dos parâmetros do modelo linear pressupõe a satisfação aos pressupostos básicos anteriormente mencionados. O princípio que norteia os cálculos é obter valores de parâmetros que minimizem a Soma dos Quadrados dos Resíduos - SQRes, ou comumente chamado de Mínimos Quadrados Ordinários - MQO.

Ou seja, para o modelo \[Y=X\beta+\varepsilon\]

em que, \(\beta\) é um vetor de parâmetros com um \(\beta\) associado a cada variável da matriz de variáveis explicativas \(X\). O conjunto dos parâmetros \(\beta\) será o conjunto de coeficientes da regressão. O \(\varepsilon\) é o termo de resíduos (muitas vezes chamado de termo de erros), ou a parte inexplicada de \(Y\) no modelo de regressão.
A estimação requer a minimização conforme a seguir:

• Forma algébrica: \[ \min \left( SQRes = \sum {\varepsilon _i^2} = \sum {\left( {Y_i-\beta _0 -\beta _1 X_{1i}-\beta _2 X_{2i}-...-\beta _k X_{ki}} \right)}^2 \right) \]

ou

• Forma matricial: \[ \min \left( SQRes = {\varepsilon}' \varepsilon \right) \]

O problema matemático é de otimizar, ou seja, minimizar um produto de um vetor linha por um vetor coluna. Portanto, deriva-se e iguala a zero obtendo a solução para o vetor de parâmetros. Segue abaixo:

\[\varepsilon '\varepsilon = {\left( {Y - X\hat \beta } \right)^\prime }\left( {Y - X\hat \beta } \right)\]

\[\varepsilon '\varepsilon = Y'Y - Y'X\hat \beta - \hat \beta 'X'Y + \hat \beta 'X'X\hat \beta \]

\[\frac{{\partial (\varepsilon '\varepsilon )}}{{\partial \hat \beta }} = - 2X'Y + 2X'X\hat \beta = 0\]

A partir dessa solução tem-se o sistema de equações normais dos mínimos quadrados e o estimador dos parâmetros (\(\hat\beta\)) por MQO. O estimador é válido para uma matriz não multicolinear de \(X\).

\[X'X\hat \beta = X'Y\]

\[\hat \beta = {\left( {X'X} \right)^{ - 1}}X'Y\]

O vetor de \(\hat\beta\) terá dimensão \((k+1)\) x \(1\) sendo \(k\) o número de variáveis explicativas. Assim, com as matrizes \(X\) e \(Y\) obtêm-se os parâmetros estimados. O estimador da variância dos resíduos será \(s^2\), para os \((n-p)\) graus de liberdade (simbolizado por \(GL\), igual ao número de observações, \(n\), menos o número de parâmetros, \(p\)):

\[ s^2 = \frac{e'e}{n - p} = \frac{SQRes}{n - p} = \frac{SQRes}{G.L.} \]

A matriz de variância-covariância dos parâmetros será:

\[\widehat {Var - Cov\left( {\hat \beta } \right)} = E\left[ {\left( {\hat \beta - \beta } \right){{\left( {\hat \beta - \beta } \right)}^\prime }} \right]\]

mas, se

\[\hat \beta = {\left( {X'X} \right)^{ - 1}}X'Y\]

então,

\[ \hat \beta = {\left( {X'X} \right)^{ - 1}}X'\left( {X\beta + \varepsilon } \right) = {\left( {X'X} \right)^{ - 1}}X'X\beta + {\left( {X'X} \right)^{ - 1}}X'\varepsilon \]

\[ \hat \beta = I.\beta + {\left( {X'X} \right)^{ - 1}}X'\varepsilon \]

\[ \hat \beta - \beta = {\left( {X'X} \right)^{ - 1}}X'\varepsilon \]

\[ Var - Cov(\hat \beta ) = E\left[ {\left( {{{\left( {X'X} \right)}^{ - 1}}X'\varepsilon } \right){{\left( {{{\left( {X'X} \right)}^{ - 1}}X'\varepsilon } \right)}^\prime }} \right]\]

\[Var - Cov(\hat \beta ) = E\left[ {{{\left( {X'X} \right)}^{ - 1}}X'\varepsilon \varepsilon 'X{{\left( {X'X} \right)}^{ - 1}}} \right]\]

Mas como as variáveis \(X\) são fixas, independentes dos resíduos (pressuposição clássica do modelo de regressão linear), o valor esperado se reduz a:

\[ Var - Cov(\hat \beta ) = {\left( {X'X} \right)^{ - 1}}X'E\left[ {\varepsilon \varepsilon '} \right]X{\left( {X'X} \right)^{ - 1}}\]

\[ Var - Cov(\hat \beta ) = {\left( {X'X} \right)^{ - 1}}X'{\sigma ^2}IX{\left( {X'X} \right)^{ - 1}}\]

Ou seja,

\[ Var - Cov(\hat \beta ) = {\sigma ^2}{\left( {X'X} \right)^{ - 1}}X'X{\left( {X'X} \right)^{ - 1}}\]

\[ Var - Cov(\hat \beta ) = {\sigma ^2}I{\left( {X'X} \right)^{ - 1}}\]

\[ Var - Cov(\hat \beta ) = {\sigma ^2}{\left( {X'X} \right)^{ - 1}}\]

ou

\[ Var - Cov(\hat \beta ) = {s^2}{\left( {X'X} \right)^{ - 1}}\]

Desta forma, têm-se as equações essenciais para a estimação. Segue a Figura 3.1 com um resumo dos estimadores de MQO.

Figura 3.1 Estimadores de Mínimos Quadrados Ordinários.

Fonte: Elaboração própria.

O valor dos erros padrões dos parâmetros será obtido a partir da raiz da variância dos parâmetros, ou seja, tirando-se a raiz da diagonal principal da \(Var-Cov(\hat \beta )\). Os parâmetros devem ter análise de significância, por meio de um teste de hipótese do tipo t:

\[ H_0:{\beta _j} = 0 \] \[ H_1:{\beta _j} \ne 0 (bilateral) \]

\[ t_{calculado} = \frac{{{{\hat \beta }_j}}}{{{s_{{{\hat \beta }_j}}}}}\sim\mathop {{t_{n - p}}}\limits_{G.L.} \]

O teste t bicaudal prevê a área de rejeição de \({H_0}\) para valores em módulo maiores que \(t_c = t_{n-p}\) para o nível de significância escolhido, que geralmente é 1%, 5% ou 10% (Figura 3.2).

Figura 3.2 Ilustração da área de rejeição do teste de hipótese bicaudal da distribuição t-Student.

Fonte: Elaboração própria.

Na Figura 3.3, têm-se as áreas para a distribuição de t para um exemplo genérico com 20 graus de liberdade e três graus de confiança: 90% (área vermelha), 95% (área vermelha mais azul), e 99% (áreas vermelha mais azul mais verde).

As áreas coloridas são de não-rejeição de \({H_0}:{\beta} = 0\), ou seja, nos bordos se têm as áreas de rejeição e no centro as áreas de não-rejeição. O que os softwares reportam são as os valores das probabilidades para valores maiores que os limites críticos, ou seja, os softwares como R, Stata e Eviews reportam os valores das probabilidades de +infinito (\(+\infty\)) ou –infinito (\(-\infty\)) até os limites das áreas pintadas. Os limites para 90% (área vermelha), 95%(área vermelha mais azul), e 99% (áreas vermelha mais azul mais verde), para por exemplo, 20 graus de liberdade serão, respectivamente: \({\pm}\) 1,724718;\({\pm}\) 2,085963; \({\pm}\) 2,845340.

Figura 3.3 Áreas para a distribuição de t para um exemplo genérico com 20 graus de liberdade e três graus de confiança: 90%, 95% e 99%.

Fonte: Elaboração própria.

Para maiores graus de liberdade, a distribuição de \(t\) converge para a distribuição normal, como pode ser observado na Figura 3.4. A distribuição teórica para df = 30 já é praticamente sobreposta à normal. Por esse motivo muitos livros falam de ter mais de 30 observações para uma regressão confiável, mas isto dependerá do número de parâmetros a estimar.

x <- seq(-4, 4, length=100)
hx <- dnorm(x)
degf <- c(1, 3, 8, 30)
colors <- c("red", "blue", "darkgreen", "gold", "black")
labels <- c("df=1", "df=3", "df=8", "df=30", "normal")
plot(x, hx, type="l", lty=2, xlab="valor da estatística",
     ylab="Densidade", main="Comparação de Distribuições t")
for (i in 1:4){
  lines(x, dt(x,degf[i]), lwd=2, col=colors[i])
}
legend("topright", inset=.05, title="Distribuições",
       labels, lwd=2, lty=c(1, 1, 1, 1, 2), col=colors)

Figura 3.4 Comparação de distribuições t-Student com a normal, em diferentes graus de liberdade (df).

Fonte: Elaboração própria a partir do script.

Os softwares econométricos em geral disponibilizam o valor da probabilidade (p-value) associado ao valor de t calculado. Desta forma, pode-se comparar com níveis predeterminados de significância para rejeitar ou não a hipótese nula. Em geral, costuma-se observar os valores das probabilidades comparando a 10%, 5% ou 1% para concluir a respeito da hipótese nula. Espera-se, para que a variável X tenha efeito não nulo sobre Y, que se rejeite a hipótese nula e que assim, os valores calculados dos parâmetros permitam uma interpretação econômica deste efeito.

Pr(>|t|) ou p-value ou valor-p é a probabilidade associada ao valor de \(t\) calculado ser maior que \(|t|\). Portanto, se Pr(>|t|) é pequeno, os coeficientes são ditos estatisticamente significativos diferentes de zero. Se Pr(>|t|), a probabilidade de \(t\) é elevada (leia-se \(>0.10\)), então os coeficientes são ditos não significativos.
> Resumo: um valor de \(t\) grande indica ser mais provável que o coeficiente estimado seja não nulo. Portanto, se desejar que haja efeitos de \(X\) sobre \(Y\), deseja-se maiores valores de t. Por sua vez, maiores valores de t indicam menores valores da probabilidade de t.

3.2 Decomposição da Variação de \(Y\)

Para auxiliar o entendimento, é possível decompor a variação de \(Y\) como na Figura.

Figura 3.5 Decomposição da variação de Y em função de X.

Fonte: Elaboração própria.

Então pode-se escrever como

\[{y_i} = {Y_i} - \bar Y\]

\[{y_i} = {\hat Y_i} + {\hat e_i} - \bar Y\]

A variação total (\(y_i\) ) será a variação explicada por \(X\) (\(\hat y_i\)) [em que \(\hat y_i=\hat Y_i-\bar Y\) é a variação devida à regressão] mais a variação não explicada (decorrente do resíduo), de modo que \(y_i=(\hat y_i+ \hat e_i)\). Assim,

\[ SQTot=SQReg+SQRes \]

em que \(SQTot\) é a soma dos quadrados totais (relativa à variação total, \(y_i^2\)); \(SQRes\) é a soma dos quadrados dos resíduos (relativa à variação não explicada) e \(SQReg\) é a soma dos quadrados da regressão (relativa à variação explicada por \(X\)). Portanto,

\[{y_i}^2 = {({\hat y_i} + {\hat e_i})^2}\]

\[\begin{array}{l} {\rm{SQTot}} = \sum\limits_{}^{} {y_i^2} = \sum {\hat y_i^2 + 2\sum\limits_{}^{} {\hat y_i^{}\hat e_i^{}} + } \sum\limits_{}^{} {\hat e_i^2} = \sum\limits_{}^{} {\left( {{Y_i} - \bar Y} \right)_{}^2} = Y'Y - n{{\bar Y}^2}\\ {\rm{SQRes}} = \sum\limits_{}^{} {\hat e_i^2} = e'e = Y'Y - \hat \beta 'X'Y\\ {\rm{SQReg}} = \sum {\hat y_i^2 = } \sum\limits_{}^{} {\left( {{{\hat Y}_i} - \bar Y} \right)_{}^2} = \hat Y'\hat Y - n{{\bar Y}^2}\\ {\rm{SQTot}} = {\rm{SQReg + SQRes}} \end{array}\].

O coeficiente de determinação (\(R^2\) = R-squared ou R quadrado) é utilizado para avaliar quanto da variação total é explicada. Define-se como:

\[{R^2} = \frac{{SQReg}}{{SQTot}} = 1 - \frac{{SQRes}}{{SQTot}}\]

Seu intervalo de variação é de zero a um em condições normais: \(0 < R^2 < 1\). Não necessariamente se descarta modelos com baixos valores de \(R^2\), mas deve-se olhar os critérios de informação (no próximo capítulo), as significâncias dos parâmetros estimados e outros indicadores de testes do modelo. Se \(SQRes=SQTot\), então \(R^2=0\); Se \(SQRes\approx 0\) então \(R^2=1\).

Ou seja, mede quanto da variação de \(Y\) está sendo explicada por variações de \(X\), ou seja, mede a qualidade do ajustamento. Procura-se estimar um modelo com o maior \(R^2\) possível. Em geral, acredita-se ter um modelo bem ajustado para valores maiores que 0,8, mas sempre se deve ter cautela quanto a esses indicadores usualmente aceitos.

Na forma matricial, o cálculo será;

\[{R^2} = \frac{{\hat \beta 'X'Y - n{{\bar Y}^2}}}{{Y'Y - n{{\bar Y}^2}}} = 1 - \frac{{Y'Y - \hat \beta 'X'Y}}{{Y'Y - n{{\bar Y}^2}}}\]
Outro indicador útil, principalmente para comparações entre modelos é o \(R^2\) ajustado ( adjusted R-squared). Ele recebe este nome, pois se faz um ajustamento de \(SQRes\) e de \(SQTot\) quanto aos graus de liberdade da respectiva variação. Assim, tem-se:

\[{\bar R^2} = 1 - \frac{{\frac{{{\rm{SQRes}}}}{{\left( {{\rm{n - p}}} \right)}}}}{{\frac{{{\rm{SQTot}}}}{{\left( {{\rm{n - 1}}} \right)}}}}\]

Em geral, quanto maior o número de variáveis \(X\), maior é o valor de \(R^2\), mas para o \(R^2\) ajustado esta regra não vale. Justamente para evitar a inclusão equivocada de variáveis explicativas é que se usa o \(R^2\) ajustado. Assim, a inclusão de uma variável irrelevante poderá elevar o valor de \(R^2\), mas não necessariamente elevará o valor de \(R^2\) ajustado.

Se \(n\) for grande e o número de parâmetros (\(p\)) pequeno em relação a \(n\), a diferença entre \(\bar R^2\) (R quadrado ajustado) e \(R^2\) será pequena. Se \(n\) for pequeno e \(p\) grande em relação a \(n\), a diferença entre ambos pode ser grande e o valor ajustado será mais importante.

Outro indicador é o Teste F da regressão ( F-statistic). Procura-se saber se o modelo tem suporte estatístico. É o Teste de significância global da regressão: os \(X^\prime s\) em conjunto explicam \(Y\) de forma significativa. A hipótese nula é de que todos os parâmetros em conjunto são nulos. A Hipótese alternativa prevê pelo menos um parâmetro não nulo:

\({H_0}:{\beta _1} = 0,{\beta _2} = 0,...,{\beta _k} = 0\);

\({H_1}:\) pelo menos um \({\beta _i} \ne 0\).

Define-se a estatística de teste \(F\) como:

\[F_{calculado} = \frac{{\frac{{{\rm{SQReg}}}}{{{\rm{p - 1}}}}}}{{\frac{{{\rm{SQRes}}}}{{{\rm{n - p}}}}}} \]

e \(F_{calculado}\) é distribuído como \(F_{(p - 1,n - p)}\) graus de liberdade.

Se \(F_{calculado} > F_{tabelado}\) , então rejeita-se H0 e concluo pela existência de ao menos um \(X\) explicando \(Y\). Deseja-se um P-value (F de significação) menor que 10%, 5% ou 1%, similarmente ao teste de \(t\) dos parâmetros.

O teste \(F\) muitas vezes é usado para outras hipóteses como para testar a especificação do modelo, ou a omissão de variáveis. Estes testes alternativos serão apresentados na seção de teste da especificação do modelo.

Esses indicadores em geral são obtidos em todos os softwares econométricos ou estatísticos. Podem-se mencionar alguns: R/RStudio, Python, Excel, Eviews, Stata, Gretl, SAS, SPSS, Gauss, e MatLab.

3.3 Propriedades dos estimadores de \(MQO\)

Os estimadores do \(MQO\) são os Melhores Estimadores Lineares Não-Tendenciosos (\(MELNT\)), ou no inglês, \(BLUE\) ( Best Linear Unbiased Estimator). Como os estimadores do modelo do Quadro 3.1 são obtidos a partir de uma amostra aleatória, os estimadores dos parâmetros (\(\hat\beta\) ), da variância-covariância dos resíduos (\(s^2\)) e dos parâmetros (\(\rm {Var-Cov}(\hat\beta)\) ) terão distribuições de probabilidade associadas, e seus valores normalmente se alteram conforme o tamanho da amostra, ou seja, com o número de observações.

Tendenciosidade

Uma propriedade desejável é a da não-tendenciosidade, ou seja, os valores esperados dos parâmetros estimados por meio da expressão \(\hat \beta=(X^\prime X)^{-1}(X^\prime Y)\) são iguais ao verdadeiro valor de \(\beta\). Ou seja, para cada \(\beta_i\), o valor esperado do estimador \(E(\hat \beta_i)=\beta_i\). Intuitivamente, ao coletar várias amostras para um fenômeno, e ao calcular os parâmetros pela expressão de \(MQO\) ( \(\hat \beta_{MQO}=(X^\prime X)^{-1}(X^\prime Y)\)), então o valor médio das estimativas das várias amostras será \(\beta\) , para um modelo estatisticamente correto, ou que suas pressuposições clássicas sejam válidas. De outro modo, se não convergir para \(\beta\), então a estimativa obtida por \(\hat \beta_{MQO}\) (pelo \(MQO\)) será dita tendenciosa ou viesada (Figura 3.6).

Figura 3.6 Estimadores: com tendenciosidade (\(\beta_{MQ2}\)) e sem tendenciosidade (\(\beta_{MQO}\)).

Fonte: Adaptado de Wooldridge (2016).

O viés de um estimador será exatamente a medida da diferença entre o valor esperado e seu valor observado. Ou seja, para uma variável \(Y\), o viés será a diferença entre o valor estimado \(\hat Y\) e o próprio \(Y\) (ou o \(Y_{observado}\)). Portanto, \({\rm viés} = E(Y)-Y=\hat Y-Y\).

Consistência

Embora tenha sido falado sobre estimadores não-tendenciosos, há casos em séries temporais, nos quais o estimador de \(MQO\) não é não-tendencioso. Neste caso, em geral os pesquisadores costumam observar a propriedade da Consistência dos estimadores, dentro do espectro chamado de propriedades assintóticas dos estimadores, ou propriedades dos estimadores para grandes amostras, ou quando \(n\) tender ao infinito (\(n\) = número de observações). Faz-se uma explicação intuitiva ao leitor. Demonstrações formais podem ser obtidas, entre outros lugares, no Apêndice C de Wooldridge (2016). Se aumentar o tamanho da amostra e o estimador \(\hat \beta\) convergir para o valor verdadeiro, então se tem matematicamente que \(\hat \beta\) é o limite de probabilidade de \(\beta\) , ou seja,
\[\mathop {{\rm{plim}}(\hat \beta )}\limits_{n \to \infty } = \beta \].

A distribuição de \(\hat \beta\) está cada vez mais concentrada em torno de seu valor verdadeiro \(\beta\) . Ou seja, para grandes amostras, a estimativa está ficando cada vez mais próxima do valor verdadeiro (Figura 3.7).

Figura 3.7 Consistência: Distribuições amostrais do parâmetro para tamanhos de amostras diferentes.

Fonte: Adaptado de Wooldridge (2016).

Eficiência

Entre os estimadores da classe dos estimadores lineares não tendenciosos, o de \(MQO\) tem a menor variância, e a isto se deve a ideia de melhor estimador. Espera-se que quanto maior a amostra, tendendo ao infinito, os estimadores de \(MQO\) tenderão aos valores verdadeiros, aos que se teria para a população, não tendenciosos, de variância mínima. Esta ideia está explicada pelo Teorema de Gauss-Markov: Em um modelo estatisticamente correto, ou seja, cujas pressuposições clássicas sejam válidas, os estimadores de Mínimos Quadrados são os melhores, os de variância mínima, quando comparados a outros estimadores lineares não-tendenciosos. Este Teorema independe da satisfação ao pressuposto da normalidade dos resíduos. Será mostrado nos capítulos de violação dos pressupostos clássicos que as variâncias estimadas, na presença daqueles problemas de violação, não serão as mínimas.

Seja a Figura 3.8. Se tiver para a distribuição de \(\beta_2\) , dois estimadores \(\hat \beta_2\) e \(\hat \beta^*_2\) , o estimador eficiente será o de menor variância. Por exemplo, sejam duas distribuições de \(\beta_2\) : a distribuição da Figura 3.8 (fdp de \(\hat \beta_2\)) ; e a distribuição da Figura 3.8 (fdp de \(\hat \beta^*_2\)). A distribuição de variância mínima será a de \(\hat \beta_2\) , a mais concentrada em torno da média \(E(\hat \beta_2)=\beta_2\) . Portanto, \(\hat \beta_2\) será dito o estimador eficiente de \(\beta_2\) , posto que \(Var({\hat \beta _2}) \le Var(\hat \beta _2^*)\) . No caso de estimadores em que algum deles é viesado, uma forma de comparar é olhando o erro quadrático médio (\(EQM\)), pois este pode ser demonstrado como a soma da variância com o quadrado do viés. O \(EQM\) medirá o quanto o estimador está longe de seu valor esperado.

Figura 3.8 Exemplificação do melhor estimador não-tendencioso dos parâmetros.

Fonte: Adaptado de Wooldridge (2016) e Gujarati (2011, p.92).

O modelo clássico de análise de regressão é construído com base numa série de pressuposições referentes ao comportamento da população. Estas pressuposições foram descritas nas seções anteriores e no próximo capítulo se inicia a discussão das violações das pressuposições.

3.4 Avaliando o(s) modelo(s)

Ante o exposto anteriormente, as principais medidas para avaliar e decidir pelo modelo são:

3.5 Apêndice - Estimação com matrizes

Neste arquivo utilizo uma série fictícia da Tabela 2.4 do livro de Gujarati e Porter (2011: p.67) das Despesas familiares de consumo semanal Y e renda familiar semanal X. É um exercício rápido de regressão com uso de matrizes e comparativamente à saída da função lm() nativa do R. O leitor pode carregar direto na Tabela 2.4 do livro de Gujarati e Porter (2011: p.67) como abaixo ou usar a saída do dput(dados) conforme colocado abaixo.

# os dados também podem ser chamados fazendo a rotina abaixo `dput(dados)`:
# A execução da linha abaixo fornece exatamente o mesmo que a linha 2
dados<-structure(list(y = c(70, 65, 90, 95, 110, 115, 120, 140, 155, 150), 
                      x = c(80, 100, 120, 140, 160, 180, 200, 220, 240, 260), 
                      obs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)), 
                 row.names = c(NA, -10L), 
                 class = c("tbl_df", "tbl", "data.frame"))
print(dados)      # funcao para visualizar os dados na tela do RStudio

##      y   x obs
## 1   70  80   1
## 2   65 100   2
## 3   90 120   3
## 4   95 140   4
## 5  110 160   5
## 6  115 180   6
## 7  120 200   7
## 8  140 220   8
## 9  155 240   9
## 10 150 260  10

library(knitr)   # pacote para a funcao kable
attach(dados)    # funcao para que o RStudio entenda os rótulos de 'dados'

# Estatisticas Descritivas

summary(dados)

##        y                x            obs       
##  Min.   : 65.00   Min.   : 80   Min.   : 1.00  
##  1st Qu.: 91.25   1st Qu.:125   1st Qu.: 3.25  
##  Median :112.50   Median :170   Median : 5.50  
##  Mean   :111.00   Mean   :170   Mean   : 5.50  
##  3rd Qu.:135.00   3rd Qu.:215   3rd Qu.: 7.75  
##  Max.   :155.00   Max.   :260   Max.   :10.00

reg1<-lm(y~x,data=dados)    
# sumario de resultados do objeto reg1: 
# saída da regressao
summary(reg1)    

## 
## Call:
## lm(formula = y ~ x, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.364  -4.977   1.409   4.364   8.364 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 24.45455    6.41382   3.813  0.00514 ** 
## x            0.50909    0.03574  14.243 5.75e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.493 on 8 degrees of freedom
## Multiple R-squared:  0.9621, Adjusted R-squared:  0.9573 
## F-statistic: 202.9 on 1 and 8 DF,  p-value: 5.753e-07

# criterios de informacao
reg1$AIC <- AIC(reg1)
reg1$BIC <- BIC(reg1)
suppressMessages(library(stargazer))
stargazer(reg1, 
          title = "Título: Resultado da Regressão OLS ou 
          Mínimos Quadrados Ordinários", 
          align = TRUE,type = "text", 
          style = "all", 
          keep.stat = c("AIC", "BIC", "rsq", "adj.rsq","n"))

## 
## Título: Resultado da Regressão OLS o
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                  y             
## -----------------------------------------------
## x                            0.509***          
##                               (0.036)          
##                             t = 14.243         
##                             p = 0.00000        
## Constant                     24.455***         
##                               (6.414)          
##                              t = 3.813         
##                              p = 0.006         
## -----------------------------------------------
## Observations                    10             
## R2                             0.962           
## Adjusted R2                    0.957           
## Akaike Inf. Crit.             69.562           
## Bayesian Inf. Crit.           70.470           
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

# Para que o R entenda os rótulos, usarei attach(dados)
attach(dados)
# Para criar a matriz X 
# coluna 1 com valores unitários e x1
X<-as.matrix.data.frame(cbind(1,dados[,2]))
print(X)

##       [,1] [,2]
##  [1,]    1   80
##  [2,]    1  100
##  [3,]    1  120
##  [4,]    1  140
##  [5,]    1  160
##  [6,]    1  180
##  [7,]    1  200
##  [8,]    1  220
##  [9,]    1  240
## [10,]    1  260

# Para transformar a variavel (y) em vetor:
y1<-as.matrix(dados[,1])
print(y1)

##       [,1]
##  [1,]   70
##  [2,]   65
##  [3,]   90
##  [4,]   95
##  [5,]  110
##  [6,]  115
##  [7,]  120
##  [8,]  140
##  [9,]  155
## [10,]  150

# Para estimar o vetor (X'X) da equacao, 
# primeiro se obtem o parametro pelo seguintes passos:
# 1) transposta de X: (X')
trX<-(t(X))
# 2) Produto da transposta de X por X, 
#    com o codigo %*%: (X'X)
X_X<-trX %*% X
X_X

##      [,1]   [,2]
## [1,]   10   1700
## [2,] 1700 322000

# Para obter a inversa (X'X)-1 se deve primeiro ativar 
# o pacote library(MASS), e usar ginv() para inverter 
# a matriz
det(X_X)     

## [1] 330000

library(MASS)
invX_X<-(ginv(X_X))
invX_X                 # (X'X)-1

##              [,1]          [,2]
## [1,]  0.975757576 -5.151515e-03
## [2,] -0.005151515  3.030303e-05

# Uma vez que se tem a inversa (X'X)-1, 
# se procede o produto X'y:
Xy<-trX %*% y1
Xy                    # X'y

##        [,1]
## [1,]   1110
## [2,] 205500

#Para calcular o vetor beta:
beta<-invX_X %*% Xy   # beta = (X'X)-1 X'y
beta

##            [,1]
## [1,] 24.4545455
## [2,]  0.5090909

summary(reg1) # comparar beta com coefficients do summary

## 
## Call:
## lm(formula = y ~ x, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.364  -4.977   1.409   4.364   8.364 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 24.45455    6.41382   3.813  0.00514 ** 
## x            0.50909    0.03574  14.243 5.75e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.493 on 8 degrees of freedom
## Multiple R-squared:  0.9621, Adjusted R-squared:  0.9573 
## F-statistic: 202.9 on 1 and 8 DF,  p-value: 5.753e-07

# agora calculamos os desvios-padrões
# obtencao de sigma quadrado estimado
yly<-t(y)%*%y
blXy<-t(beta)%*%Xy
ulu<-yly - blXy
ulu                  # Soma dos Quadrados dos Residuos

##          [,1]
## [1,] 337.2727

uhat<-y-X%*%(beta)
uhat                 # residuos estimados

##              [,1]
##  [1,]   4.8181818
##  [2,] -10.3636364
##  [3,]   4.4545455
##  [4,]  -0.7272727
##  [5,]   4.0909091
##  [6,]  -1.0909091
##  [7,]  -6.2727273
##  [8,]   3.5454545
##  [9,]   8.3636364
## [10,]  -6.8181818

# obter sigma quadrado para gl=n-p=8
n<-nrow(dados)       # numero de observacoes
k<-ncol(X)          # numero de parametros a estimar
sigsqhat<-as.numeric(ulu/(n-k))
sigsqhat            # s2 - variancia estimada do residuo

## [1] 42.15909

# variancia-covariancia de beta
varcovbeta<-sigsqhat*invX_X     
varcovbeta 

##            [,1]         [,2]
## [1,] 41.1370523 -0.217183196
## [2,] -0.2171832  0.001277548

# obtendo a raiz da variancia da diagonal de varcovbeta
#desvio padrao de beta
sebeta<-sqrt(diag(varcovbeta))  
sebeta

## [1] 6.41381730 0.03574281

# parametros beta
beta

##            [,1]
## [1,] 24.4545455
## [2,]  0.5090909

# estatistica t dos parametros
tbeta<-beta/sebeta
tbeta                         

##           [,1]
## [1,]  3.812791
## [2,] 14.243171

#obtendo a probabilidade de tbeta para 5% e df=n-k
pvalue<-2*pt(-abs(tbeta),n-k)
pvalue        # probabilidade de t dos parametros

##              [,1]
## [1,] 5.142172e-03
## [2,] 5.752746e-07

# soma dos quadrados totais 
sqtot<- yly-n*mean(y1)^2
sqtot     

##      [,1]
## [1,] 8890

# soma dos quadrados dos residuos
sqres<- ulu
sqres      

##          [,1]
## [1,] 337.2727

# soma dos quadrados da regressao
sqreg<- sqtot-sqres
sqreg    

##          [,1]
## [1,] 8552.727

# coeficiente de determinação R2
r2<-sqreg/sqtot
r2               

##           [,1]
## [1,] 0.9620616

# coeficiente de determinação ajustado R2-ajustado
r2aj<-1-(sqres/(n-k))/(sqtot/(n-1))
r2aj                          

##           [,1]
## [1,] 0.9573193

F<-(sqreg/(k-1))/(sqres/(n-k))
F                        

##          [,1]
## [1,] 202.8679

# probabilidade de F para H0: betas iguais a zero
probF<-1-pf(F,k-1,n-k)
probF                         

##              [,1]
## [1,] 5.752746e-07

4.1 Teste RESET

Entre outros testes, o teste \(RESET\) de @ramsey1969tests é um dos mais aplicados na literatura. O nome vem do pesquisador Ramsey para o Regression Specification Error Test ou Teste de Erro de Especificação da Regressão. O teste é baseado na regressão aumentada

\[Y = X\beta + Z\alpha + \varepsilon ,\]

em que \(X\) é uma matriz das variáveis explicativas e \(Z\) é uma matriz das variáveis dependentes estimadas e elevadas a uma potência tal que, \(Z = [ \hat Y^2, \hat Y^3, \hat Y^4]\), (neste exemplo, para três fitted terms, três termos acrescentados na regressão aumentada, em que \(\hat Y\) é a variável \(Y\) prevista pelo modelo estimado originalmente: \(Y = X\beta + \varepsilon\)).

A ideia é olhar a significância do vetor de parâmetros \(\alpha\) para ver se os termos acrescentados são relevantes no modelo, indicando erro de especificação.

O procedimento do teste será:

1. Estima-se \(Y = X\beta + \varepsilon\); nesta regressão inicial, será obtido o \(R_{velho}^2\);

2. Obtêm-se os valores previstos de \(Y\) e gera-se \(Z = [ \hat Y^2, \hat Y^3, \hat Y^4]\) ou mais se desejar. Recomenda-se no máximo até 3 termos, ou seja, até \(\hat Y^4]\);

3. Ajusta-se a regressão aumentada, colocando-se os \(X\) e as variáveis do item 2 : \(Y = f(X, \hat Y^2, \hat Y^3\) ); nesta regressão aumentada, será obtido o \(R_{novo}^2\);

4. Com as regressões de 1 e de 3, observam-se os valores de \(R_{novo}^2\) (de 3) e \(R_{velho}^2\) (de 1) e calcula-se a estatística de teste;

5. Estatística de Teste;

\[\begin{equation} $$F = \frac{{\frac{{R_{novo}^2 - R_{velho}^2}}{{{\rm{número\ de \ novos\ regressores }}(m)}}}}{{\frac{{1 - R_{novo}^2}}{{n - {\rm{número \ de \ parâmetros\ no\ novo\ modelo (p)}}}}}};$$ \end{equation}\]

em que \(F ∼ F_{(m, n-p)}\) ; \(m\) é o número de novos regressores (potências incluídas); \((n-p)\) é o número de observações menos o número de parâmetros no novo modelo.

6. Comparar o \(F\) do item 5 com o \(F\) da tabela teórica de distribuição de \(F\) , para o nível de significância, e graus de liberdade do numerador \(m\) e denominador \((n-p)\). Como a hipótese nula (\(H_0\)) é de que não há erro de especificação, espera-se que a hipótese nula não seja rejeitada, ou seja, que \(F-statistic\) seja muito pequena e sua respectiva probabilidade (\(p-value\)) seja elevada.

O teste \(RESET\) indica apenas se o modelo está especificado incorretamente, mas não diz qual seria a solução. A solução para um problema envolverá: incluir outras variáveis relevantes no modelo (variáveis novas ou as \(X\) nas potências); retirar as irrelevantes; ou mudar a forma funcional (por exemplo, alterando para log, CES, Translog ou outra). Portanto, o bom senso indica que é melhor incluir variáveis do que excluir, pois a exclusão pode causar viés, enquanto a inclusão tende a melhorar o modelo, a não ser pela possibilidade de não-significância dos parâmetros.

A inclusão de variáveis pode envolver o uso de dummies para melhorar o detalhamento e/ou discriminar grupos ou categorias de interesse, ou ainda o uso de variável de tendência em séries temporais, ou mesmo para destacar possíveis outliers. A exclusão ou transformação de variáveis pode ser uma opção caso esteja havendo multicolinearidade ou outra violação de pressuposto como resíduos não normais.

No software \(R\), o teste \(RESET\) pode ser feito manualmente, por meio de ajustamento da regressão, obtenção dos valores ‘ajustados’ ou ‘estimados’ e estimação da regressão de teste como nos procedimentos 1 a 6 mencionados nesta seção. Outra forma é usar o comando resettest do pacote (library) lmtest. Serão mostrados no script os dois procedimentos.

library(readxl)
# library(foreign)
dados <- read_excel("soja_apostila.xlsx", 
                    sheet = "dados")
summary(dados)

##       OBS          QSOJA        FERTILIZANTE        TRATOR     
##  Min.   :  1   Min.   :204.6   Min.   : 8.974   Min.   :2.218  
##  1st Qu.: 30   1st Qu.:294.9   1st Qu.:16.592   1st Qu.:3.269  
##  Median : 59   Median :316.7   Median :18.450   Median :4.059  
##  Mean   : 59   Mean   :322.3   Mean   :18.404   Mean   :4.317  
##  3rd Qu.: 88   3rd Qu.:350.2   3rd Qu.:21.392   3rd Qu.:5.455  
##  Max.   :117   Max.   :506.5   Max.   :27.079   Max.   :6.677  
##        MO         
##  Min.   :0.04820  
##  1st Qu.:0.05539  
##  Median :0.07085  
##  Mean   :0.08165  
##  3rd Qu.:0.09346  
##  Max.   :0.26126

attach(dados)
# QSOJA = quantidade produzida de soja; 
# FERTILIZANTE = quantidade utilizada de fertilizantes, 
# TRATOR = número de horas-máquina utilizadas, e 
# MO = quantidade de mão-de-obra em número de pessoas 

## RESET de Ramsey
### Execucao manual do RESET de Ramsey para especificação

regressao1<-lm(QSOJA~FERTILIZANTE+TRATOR+MO)
summary(regressao1)

## 
## Call:
## lm(formula = QSOJA ~ FERTILIZANTE + TRATOR + MO)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.714 -33.171   1.768  24.894 149.637 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   494.9657    25.5723  19.356  < 2e-16 ***
## FERTILIZANTE   -0.5535     1.0589  -0.523   0.6022    
## TRATOR        -33.6899     3.7410  -9.006 6.09e-15 ***
## MO           -209.1407   107.8926  -1.938   0.0551 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41.51 on 113 degrees of freedom
## Multiple R-squared:  0.4651, Adjusted R-squared:  0.4509 
## F-statistic: 32.75 on 3 and 113 DF,  p-value: 2.608e-15

### fazendo os criterios de informacao de Akaike e Schwarz
regressao1$AIC <- AIC(regressao1)
regressao1$BIC <- BIC(regressao1)
#mostrando os valores de AIC e SIC
suppressMessages(library(stargazer))

star.1 <- stargazer(regressao1,
                    title="Título: Resultado da Regressão",
                    align=TRUE,
                    type = "text", style = "all",
                    keep.stat=c("aic","bic","rsq", "adj.rsq","n"))

## 
## Título: Resultado da Regressão
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                QSOJA           
## -----------------------------------------------
## FERTILIZANTE                  -0.554           
##                               (1.059)          
##                             t = -0.523         
##                              p = 0.603         
## TRATOR                      -33.690***         
##                               (3.741)          
##                             t = -9.006         
##                              p = 0.000         
## MO                           -209.141*         
##                              (107.893)         
##                             t = -1.938         
##                              p = 0.056         
## Constant                    494.966***         
##                              (25.572)          
##                             t = 19.356         
##                              p = 0.000         
## -----------------------------------------------
## Observations                    117            
## R2                             0.465           
## Adjusted R2                    0.451           
## Akaike Inf. Crit.            1,209.807         
## Bayesian Inf. Crit.          1,223.617         
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

I(fitted(regressao1)^2)+
I(fitted(regressao1)^3)+
I(fitted(regressao1)^4)

reg_RESET_3<-lm(QSOJA~FERTILIZANTE+TRATOR+MO+I(fitted(regressao1)^2)+
                  I(fitted(regressao1)^3)+I(fitted(regressao1)^4),data=dados)

reg_RESET<-lm(QSOJA~FERTILIZANTE+TRATOR+MO+
                I(fitted(regressao1)^2)+I(fitted(regressao1)^3),data=dados)

results<-stargazer(list(regressao1,reg_RESET_3),
                   column.labels = c("regressao","RESET3"),
                   type="text",style="all",
                   keep.stat=c("aic","bic","rsq", "adj.rsq","n"))

## 
## ===================================================
##                            Dependent variable:     
##                        ----------------------------
##                                   QSOJA            
##                          regressao       RESET3    
##                             (1)           (2)      
## ---------------------------------------------------
## FERTILIZANTE              -0.554       304.130**   
##                           (1.059)      (135.047)   
##                         t = -0.523     t = 2.252   
##                          p = 0.603     p = 0.027   
## TRATOR                  -33.690***    18,591.290** 
##                           (3.741)     (8,231.767)  
##                         t = -9.006     t = 2.258   
##                          p = 0.000     p = 0.026   
## MO                       -209.141*   115,237.700** 
##                          (107.893)    (51,069.360) 
##                         t = -1.938     t = 2.256   
##                          p = 0.056     p = 0.027   
## I(fitted(regressao1)2)                  2.665**    
##                                         (1.165)    
##                                        t = 2.287   
##                                        p = 0.025   
## I(fitted(regressao1)3)                  -0.006**   
##                                         (0.002)    
##                                        t = -2.300  
##                                        p = 0.024   
## I(fitted(regressao1)4)                 0.00000**   
##                                        (0.00000)   
##                                        t = 2.303   
##                                        p = 0.024   
## Constant                494.966***   -230,604.700**
##                          (25.572)    (101,861.300) 
##                         t = 19.356     t = -2.264  
##                          p = 0.000     p = 0.026   
## ---------------------------------------------------
## Observations                117           117      
## R2                         0.465         0.532     
## Adjusted R2                0.451         0.507     
## Akaike Inf. Crit.        1,209.807                 
## Bayesian Inf. Crit.      1,223.617                 
## ===================================================
## Note:                   *p<0.1; **p<0.05; ***p<0.01

F<- ((0.532456 - 0.46511)/3)/((1 - 0.532456)/(117-7))
F

## [1] 5.281542

# probabilidade de F para H0: betas iguais a zero
probF<-1-pf(F,3,117-7)
probF                         

## [1] 0.001932346

RESETH0<-c("I(fitted(regressao1)^2)","I(fitted(regressao1)^3)",        
"I(fitted(regressao1)^4)") 

suppressMessages(library(car))
RESETH0<-c("I(fitted(regressao1)^2)","I(fitted(regressao1)^3)",
        "I(fitted(regressao1)^4)")
Tabela_RESET<-linearHypothesis(reg_RESET_3,RESETH0)
# outra alternativa é usar a linha abaixo com o matchCoefs
#Tabela_RESET<-linearHypothesis(reg_RESET_3, matchCoefs(reg_RESET,"fitted"))
Tabela_RESET

RESETH0<-c("I(fitted(regressao1)^2)","I(fitted(regressao1)^3)")

results<-stargazer(list(regressao1,reg_RESET),
                   type="text",style="all",
                   column.labels = c("regressao1","RESET"),
                   keep.stat=c("aic","bic","rsq", "adj.rsq","n") )

## 
## ===================================================
##                            Dependent variable:     
##                        ----------------------------
##                                   QSOJA            
##                          regressao1       RESET    
##                             (1)            (2)     
## ---------------------------------------------------
## FERTILIZANTE               -0.554        -5.647    
##                           (1.059)       (12.003)   
##                          t = -0.523    t = -0.470  
##                          p = 0.603      p = 0.639  
## TRATOR                   -33.690***     -291.821   
##                           (3.741)       (727.987)  
##                          t = -9.006    t = -0.401  
##                          p = 0.000      p = 0.690  
## MO                       -209.141*     -1,914.940  
##                          (107.893)     (4,501.175) 
##                          t = -1.938    t = -0.425  
##                          p = 0.056      p = 0.672  
## I(fitted(regressao1)2)                   -0.014    
##                                          (0.068)   
##                                        t = -0.205  
##                                         p = 0.839  
## I(fitted(regressao1)3)                   0.00000   
##                                         (0.0001)   
##                                         t = 0.057  
##                                         p = 0.955  
## Constant                 494.966***     3,168.566  
##                           (25.572)     (8,428.683) 
##                          t = 19.356     t = 0.376  
##                          p = 0.000      p = 0.708  
## ---------------------------------------------------
## Observations                117            117     
## R2                         0.465          0.510    
## Adjusted R2                0.451          0.488    
## Akaike Inf. Crit.        1,209.807                 
## Bayesian Inf. Crit.      1,223.617                 
## ===================================================
## Note:                   *p<0.1; **p<0.05; ***p<0.01

suppressMessages(library(car))
RESETH0<-c("I(fitted(regressao1)^2)","I(fitted(regressao1)^3)")
Tabela_RESET<-linearHypothesis(reg_RESET,RESETH0)
# outra alternativa é usar a linha abaixo com o matchCoefs
#Tabela_RESET<-linearHypothesis(reg_RESET, matchCoefs(reg_RESET,"fitted"))
Tabela_RESET

## Linear hypothesis test
## 
## Hypothesis:
## I(fitted(regressao1)^2) = 0
## I(fitted(regressao1)^3) = 0
## 
## Model 1: restricted model
## Model 2: QSOJA ~ FERTILIZANTE + TRATOR + MO + I(fitted(regressao1)^2) + 
##     I(fitted(regressao1)^3)
## 
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    113 194668                                
## 2    111 178360  2     16308 5.0746 0.007783 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

suppressMessages(library(lmtest))
# default é power = 2:3
TesteRESET<-resettest(regressao1, power = 2:3) 
TesteRESET

## 
##  RESET test
## 
## data:  regressao1
## RESET = 5.0746, df1 = 2, df2 = 111, p-value = 0.007783

#alterando as potencias
TesteRESET.power<-resettest(regressao1, power = 2:4)
TesteRESET.power

## 
##  RESET test
## 
## data:  regressao1
## RESET = 5.2816, df1 = 3, df2 = 110, p-value = 0.001932

4.2 Testes para omissão de variáveis

O teste para inclusão ou omissão de variáveis é um teste clássico de distribuição \(F\) em que se compara um modelo irrestrito com outro restrito, ou seja, cujas variáveis tenham coeficientes nulos. O modelo irrestrito será aquele com todas as variáveis que se deseja testar, enquanto o restrito terá menos variáveis.

Para exemplificar estes testes, faz-se uso do exercício 7.19 de Gujarati (2011, p.236) para os dados da tabela 7.9 com a Demanda por frangos nos Estados Unidos, de 1960 a 1982, com os dados do arquivo <gujarati 5ed p236 frangos tabela7_9.xlsx>. O leitor pode fazer uso também do arquivo complementar disponível em http://www.rpubs.com/amrofi/exercicio_gujarati_7_19. Observe que são colocadas quatro possíveis especificações, com as variáveis em logaritmos e a variável \(Y\) sendo dependente em todas. As alterações são para as variáveis \(X_2\) até \(X_5\). Neste exemplo, não se está preocupado com a interpretação econômica, para fins do exemplo e, portanto, apenas trata-se das variáveis como \(X\).

Figura 4.1 QR Code para o exemplo de Demanda de frangos. X.

• Exemplo gujarati_7_19_frangos:

7.19. Demanda por frangos nos Estados Unidos, 1960-1982. Para estudar o consumo per capita de frango nos Estados Unidos, use os dados da Tabela 7.9, em que Y = consumo per capita de frango em libras (peso) X2= renda real disponível per capita, em $ X3= preço real do frango no varejo, em centavos de dólar por libra (peso) ¢ X4= preço real da carne suína no varejo, em centavos de dólar por libra (peso) ¢ X5= preço real da carne bovina no varejo, em centavos de dólar por libra (peso) ¢ X6= preço real dos substitutos da carne de frango, em centavos de dólar por libra (peso), que é uma média ponderada dos preços reais das carnes suína e bovina, usando como pesos o consumo relativo de cada uma dessas carnes em relação ao consumo total delas. Agora, considere as seguintes funções de demanda:

1. \[ ln{Y_t} = \alpha _1 + \alpha_2 X_{2t} + \alpha _3 X_{3t} + u_t \]
2. \[ ln{Y_t} = \gamma _1 + \gamma _2 X_{2t} + \gamma _3 X_{3t} + \gamma _4 X_{4t} + u_t \]
3. \[ ln{Y_t} = \lambda _1 + \lambda _2 X_{2t} + \lambda _3 X_{3t} + \lambda _4 X_{5t} + u_t \]
4. \[ ln{Y_t} = \theta _1 + \theta _2 X_{2t} + \theta _3 X_{3t} + \theta _4 X_{4t} + \theta _5 X_{5t} + u_t \]
5. \[ ln{Y_t} = \beta _1 + \beta _2 X_{2t} + \beta _3 X_{3t} + \beta _4 X_{6t} + u_t \]

Chamando os dados de gujarati 5ed p236 frangos tabela7_9.xlsx:

suppressMessages(library(readxl))
frangos <- read_excel("gujarati 5ed p236 frangos tabela7_9.xlsx", 
                                 sheet = "dados")
#Y  Per Capita Consumption of Chickens, Pounds                      
#X2 Real Disposable Income Per Capita, $                        
#X3 Real Retail Price of Chicken Per Pound, Cents                       
#X4 Real Retail Price of Pork Per Pound, Cents                      
#X5 Real Retail Price of Beef Per Pound, Cents                      
#X6 Composite Real Price of Chicken Substitutes Per Pound, Cents                        
attach(dados)

Vamos ver como está a tabela importada:

library(knitr)
kable(frangos[,1:7],
      caption="Dados da tabela 7-9","latex", booktabs = T)

Estimando o modelo linear de regressao múltipla fazendo as regressões de \(Y\) contra as variáveis \(X_2, X_3, X_4, X_5\) e \(X_6\) com logs. Agora estimam-se as regressões com uso da função lm e com operadores de logaritmos. As equações tomam a forma:

EQ1<-lm(log(Y)~log(X2)+log(X3),data = frangos)
EQ2<-lm(log(Y)~log(X2)+log(X3)+log(X4),data = frangos)
EQ3<-lm(log(Y)~log(X2)+log(X3)+        log(X5),
        data = frangos)
EQ4<-lm(log(Y)~log(X2)+log(X3)+log(X4)+log(X5),
        data = frangos)
EQ5<-lm(log(Y)~log(X2)+log(X3)+                +log(X6),
        data = frangos)

As equações são como a seguir:

1. log(Y)~log(X2)+log(X3);

2. log(Y)~log(X2)+log(X3)+log(X4)

3. log(Y)~log(X2)+log(X3) +log(X5))

4. log(Y)~log(X2)+log(X3)+log(X4)+log(X5))

Os resultados de diferentes especificações pode ser visto:

Quadro 4.14. Resultados das estimações com diferentes especificações.

library(stargazer)
resultados<-stargazer(EQ1,EQ4,
                      type="text",style=("all"),
                      column.labels = c("EQ1","EQ4"),
                    keep.stat=c("aic","bic","rsq", "adj.rsq","n"))

## 
## =========================================
##                  Dependent variable:     
##              ----------------------------
##                         log(Y)           
##                   EQ1            EQ4     
##                   (1)            (2)     
## -----------------------------------------
## log(X2)         0.452***      0.343***   
##                 (0.025)        (0.083)   
##                t = 18.284     t = 4.114  
##                p = 0.000      p = 0.001  
## log(X3)        -0.372***      -0.505***  
##                 (0.063)        (0.111)   
##                t = -5.865    t = -4.550  
##               p = 0.00001    p = 0.0003  
## log(X4)                         0.149    
##                                (0.100)   
##                               t = 1.490  
##                               p = 0.154  
## log(X5)                         0.091    
##                                (0.101)   
##                               t = 0.905  
##                               p = 0.378  
## Constant        2.033***      2.190***   
##                 (0.116)        (0.156)   
##                t = 17.497    t = 14.063  
##                p = 0.000      p = 0.000  
## -----------------------------------------
## Observations       23            23      
## R2               0.980          0.982    
## Adjusted R2      0.978          0.978    
## =========================================
## Note:         *p<0.1; **p<0.05; ***p<0.01

Assim, formalmente, a equação irrestrita será a equação 4 (contendo todas as variáveis) e os testes de especificação serão para as hipóteses a saber:

\(H_0\): coeficientes iguais a zero para \(X_4\) e \(X_5\).

A rejeição de “\(H_0\): coeficientes iguais a zero para \(X_4\) e \(X_5\)” implica que a escolha pela equação 1 (portanto, restrita ao não incluir \(X_4\) e \(X_5\)) representará erro de especificação por omissão de variável e estará viesando o modelo.

O teste \(F\) é dado pela expressão a seguir em que \(SSR\) é a soma dos quadrados dos resíduos, respectivamente dos modelos restrito (\(SSR_r\)) e irrestrito (\(SSR_{ur}\)), para \(n\) observações, \(k\) variáveis explicativas e \(q\) restrições. A expressão pode ser calculada em termos das estatísticas \(R^2\) dos modelos irrestrito (\(R_{ur}^2\)) e restrito (\(R_r^2\)). O script está colocado na sequência das fórmulas.

\(F = \frac{{SSRr - SSRur}}{{SSRur}}.\frac{{n - k - 1}}{q} = \;\frac{{\frac{{SSRr - SSRur}}{q}}}{{\frac{{SSRur}}{{n - k - 1}}}}\)

ou

\(F = \frac{{R_{ur}^2 - R_r^2}}{{1 - R_{ur}^2}}.\frac{{n - k - 1}}{q} = \;\frac{{\frac{{R_{ur}^2 - R_r^2}}{q}}}{{\frac{{1 - R_{ur}^2}}{{n - k - 1}}}}\)

\(F ∼ F\left( {q,n - k - 1} \right)\) e nível de confiança.

#modelo irrestrito é o da EQ4
# Regressão OLS irrestrita
res.ur <- lm(log(Y)~log(X2)+log(X3)+log(X4)+log(X5),
             data = frangos)

#modelo restrito para X4 e X5 serem nulos, equivalente ao EQ1
# Regressão OLS restrita
res.r <- lm(log(Y)~log(X2)+log(X3),data = frangos)

# extração dos R2 de cada expressão:
( r2.ur <- summary(res.ur)$r.squared )

## [1] 0.9823133

( r2.r <- summary(res.r)$r.squared )

## [1] 0.9800744

# Cálculo da estatística F para (n-k-1) no numerador 
# e q restrições no denominador:
( F <- (r2.ur-r2.r) / (1-r2.ur) * (23-4-1)/2 )

## [1] 1.139245

# p value = 1-cdf da distribuição F:
1-pf(F, 2,18)

## [1] 0.3420835

res.ur<- lm(log(Y)~log(X2)+log(X3)+log(X4)+log(X5),
            data = frangos)
# Load package "car" (deve estar instalado em seu micro)
suppressMessages(library(car))

# F test: 
myH0 <- c("log(X4)","log(X5)")
linearHypothesis(res.ur, myH0)

## Linear hypothesis test
## 
## Hypothesis:
## log(X4) = 0
## log(X5) = 0
## 
## Model 1: restricted model
## Model 2: log(Y) ~ log(X2) + log(X3) + log(X4) + log(X5)
## 
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     20 0.015437                           
## 2     18 0.013703  2 0.0017345 1.1392 0.3421

4.3 Critérios de Informação de Akaike (\(AIC\)) e Schwarz (\(BIC\))

Outra forma de avaliar diferentes modelos é olhar os diferentes resultados e comparar o \(R^2ajustado\). Quanto mais próximo de 1 melhor será a estimação. Entretanto, este indicador é deficiente para o caso de variável omitida.

Outras opções são observar os valores dos chamados critérios de informação de Akaike (\(AIC\) de Akaike’s Information Criterion) e Schwarz (\(SIC\) de Schwarz’s Information Criterion ou \(BIC\) para a sigla de Bayesian Information criteria). Menores coeficientes \(AIC\) e \(SIC\) indicam melhores ajustamentos da regressão, mas só podem ser comparados se as unidades das variáveis das diferentes regressões forem as mesmas (por exemplo, não se aplica numa comparação entre \(Y\) e outra com \({\rm Log}Y\)). Deve-se olhar todos os critérios para melhor análise dos resultados.

Os Critérios de Informação de Akaike (ou \(AIC\)) (Akaike, 1974) ou de Schwarz (\(BIC\) - usaremos \(BIC\) para ficar aderente à função do R) (Schwarz, 1978) dependem da maximização da função de log verossimilhança do modelo. Desta forma, podem ser definidos como:

\(AIC = - \frac{{2l}}{n} + \frac{{2k}}{n}\)

\(SIC = - \frac{{2l}}{n} + \frac{{k.\log n}}{n}\)

\(l = - \frac{n}{2}\left( {1 + {\rm{log}}\left( {2\pi } \right) + \log \left( {\frac{{{{\hat \varepsilon }^{\rm{'}}}\hat \varepsilon }}{n}} \right)} \right)\)

em que \(k\) é o número de regressores incluindo-se o intercepto; \(n\) é o número de observações; \(l\) é o log Verossimilhança da regressão; e \(\hat \varepsilon\) são os resíduos estimados do modelo. No formato mais simplificado exposto por Greene (2002), tem-se:

\(AIC = \log \left( {\frac{{{{\hat \varepsilon }^{\rm{'}}}\hat \varepsilon }}{n}} \right) + \frac{{2k}}{n}\)

\(SIC = BIC = \log \left( {\frac{{{{\hat \varepsilon }^{\rm{'}}}\hat \varepsilon }}{n}} \right) + \frac{{k.\log n}}{n}\)

• critérios de informação pelo R

Já utilizamos as funções para o cálculo do critério de informação de Akaike e de Schwarz para os exemplos da soja e da demanda de frangos. As funções básicas chamam-se AIC e BIC. Assim, basta fazer, para obter os resultados de cada regressão:

AIC(regressao)
BIC(regressao)

Para o exemplo da soja, a regressão é reproduzida e o cálculo dos critérios e feito de forma a incluir no objeto de regressão os respectivos valores (a inclusão é feita pelo comando regressao1$AIC para o AIC e regressao1$BIC para o BIC). Esta inclusão no mesmo objeto da regressao1 facilitará a extração das tabelas contendo já o AIC e BIC na função stargazer. O leitor poderá observar que no summary da regressão não aparece nenhum critério de informação. Eles são calculados a partir do objeto resultado da função lm.

regressao1<-lm(QSOJA~FERTILIZANTE+TRATOR+MO)
summary(regressao1)

## 
## Call:
## lm(formula = QSOJA ~ FERTILIZANTE + TRATOR + MO)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.714 -33.171   1.768  24.894 149.637 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   494.9657    25.5723  19.356  < 2e-16 ***
## FERTILIZANTE   -0.5535     1.0589  -0.523   0.6022    
## TRATOR        -33.6899     3.7410  -9.006 6.09e-15 ***
## MO           -209.1407   107.8926  -1.938   0.0551 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41.51 on 113 degrees of freedom
## Multiple R-squared:  0.4651, Adjusted R-squared:  0.4509 
## F-statistic: 32.75 on 3 and 113 DF,  p-value: 2.608e-15

### fazendo os criterios de informacao de Akaike e Schwarz
regressao1$AIC <- AIC(regressao1)
regressao1$BIC <- BIC(regressao1)
regressao1$AIC

## [1] 1209.807

regressao1$BIC

## [1] 1223.617

Deve-se esclarecer que o objetivo Não é comparar o \(AIC\) com o \(BIC\). O objetivo é comparar o \(AIC\) de um modelo com o de outro modelo.

Assim, faremos uso do exemplo da demanda de frangos para comparar modelos alternativos. Sejam modelos com especificações de forma funcional distintas como log-log, lin-log, log-lin e translog. Para simplificar, considero apenas as variáveis \(X_2\) e \(X_3\) por questões didáticas. Os modelos serão do tipo:

# especificação
EQloglog<-lm(log(Y)~log(X2)+log(X3),data = frangos)
EQlinlog<-lm(Y~log(X2)+log(X3),data = frangos)
EQloglin<-lm(log(Y)~X2+X3,data = frangos)
EQtranslog<-lm(log(Y)~log(X2)+log(X3)+I(log(X2)*log(X3))
               +I(0.5*(log(X2)^2))+I(0.5*(log(X3)^2)),data = frangos)

Uma alternativa é extrair os resultados passo a passo para construir uma tabela resumo. Isto requer conhecimento das posições dos valores desejados em cada caso. Por exemplo, armazenei os quatro valores de AIC e os quatro de BIC dentro dos respectivos vetores AIC e BIC. Da mesma forma peguei os \(R^2\) e \(R^2ajustado\) de cada modelo e armazenei em vetores r2 e r2.adj. Veja, para o modelo 1, que o \(R^2\) estava dentro do summary(mod1)$r.squared; similarmente, para o \(R^2ajustado\) estava dentro de summary(mod1)$adj.r.squared. As colunas foram então combinadas (cbind) lado a lado na tab2 e posteriormente exibidas pelo knitr::kable.

# alternativa de extrair resultados
mod1<-EQloglog
mod2<-EQlinlog
mod3<-EQloglin
mod4<-EQtranslog
AIC<-rbind(AIC(mod1),AIC(mod2),AIC(mod3),AIC(mod4))
BIC<-rbind(BIC(mod1),BIC(mod2),BIC(mod3),BIC(mod4))
r2<-rbind(summary(mod1)$r.squared,
          summary(mod2)$r.squared,
          summary(mod3)$r.squared,
          summary(mod4)$r.squared)
r2.adj<-rbind(summary(mod1)$adj.r.squared,
              summary(mod2)$adj.r.squared,
              summary(mod3)$adj.r.squared,
              summary(mod4)$adj.r.squared)
tab2<-data.frame(cbind(r2,r2.adj,AIC,BIC))
row.names(tab2)<-c("log-log","lin-log","log-lin","translog")
library(knitr)
kable(tab2,
      caption = "Comparação de Modelos, 'Gujarati-frangos'", digits = 3,
      col.names=c("R²","R²Aj","AIC","BIC"),
      "latex", 
      booktabs = T)  

Primeira observação importante: Não se compara valores dos critérios AIC e BIC entre modelos com log e sem log. Assim, cabe comparação de AIC entre os modelos log-log e translog, assim como para o BIC. Já as estatísticas de \(R^2ajustado\) podem ser comparadas. Assim, o modelo de melhor ajuste considerando o \(R^2ajustado\) resultou empate entre o log-log e o translog com \(R^2ajustado=0.978\). Olhando o \(R^2\), o translog foi um pouco melhor que o log-log. Como ambos são log dos dois lados da expressão, podemos comparar seus AIC: o menor foi o do log-log (por causa do sinal negativo), \(AIC=-94.777\). Considerando o BIC, temos a mesma interpretação: log-log foi menor \(BIC=-90.235\). A escolha seria nesse caso para o log-log. Isto indica que a inclusão dos termos cruzados entre \(X_2\) e \(X_3\) assim como seus quadrados, não resultou em melhor especificação. Fazendo a emissão da tabela comparativa dos dois modelos temos a tabela abaixo. É possível constatar que nenhum dos termos cruzados e quadráticos não foram significativos estatisticamente, assim como das variáveis \(X_2\) e \(X_3\). Suspeita-se de uma multicolinearidade neste modelo translog.

mod1$AIC <- AIC(mod1) # incluir AIC no objeto da regressao
mod1$BIC <- BIC(mod1) # incluir BIC no objeto da regressao
mod4$AIC <- AIC(mod4) # incluir AIC no objeto da regressao
mod4$BIC <- BIC(mod4) # incluir BIC no objeto da regressao
resultados<-stargazer(mod1,mod4,type="text",style=("all"),
           keep.stat=c("aic","bic","rsq", "adj.rsq","n"),
           column.labels = c("log-log","translog"))

## 
## =================================================
##                          Dependent variable:     
##                      ----------------------------
##                                 log(Y)           
##                         log-log       translog   
##                           (1)            (2)     
## -------------------------------------------------
## log(X2)                 0.452***        1.050    
##                         (0.025)        (0.787)   
##                        t = 18.284     t = 1.333  
##                        p = 0.000      p = 0.200  
## log(X3)                -0.372***       -2.306    
##                         (0.063)        (3.161)   
##                        t = -5.865    t = -0.729  
##                       p = 0.00001     p = 0.476  
## I(log(X2) * log(X3))                    0.088    
##                                        (0.470)   
##                                       t = 0.187  
##                                       p = 0.855  
## I(0.5 * (log(X2)2))                    -0.142    
##                                        (0.182)   
##                                      t = -0.779  
##                                       p = 0.447  
## I(0.5 * (log(X3)2))                     0.364    
##                                        (1.514)   
##                                       t = 0.240  
##                                       p = 0.813  
## Constant                2.033***        3.698    
##                         (0.116)        (4.098)   
##                        t = 17.497     t = 0.902  
##                        p = 0.000      p = 0.380  
## -------------------------------------------------
## Observations               23            23      
## R2                       0.980          0.983    
## Adjusted R2              0.978          0.978    
## Akaike Inf. Crit.       -94.777        -92.306   
## Bayesian Inf. Crit.     -90.235        -84.357   
## =================================================
## Note:                 *p<0.1; **p<0.05; ***p<0.01

4.4 Variáveis especiais e outliers: dummies e tendência

4.4.1 Variáveis dummy, binárias ou categóricas

As variáveis dummy ou categóricas ou binárias, são variáveis atípicas que buscam diferenciar grupos ou estados ou características comuns a partes da amostra. São atípicas pois a ordem não indica superioridade ou grandeza do valor, por exemplo, uma variável categórica com valor igual 2 numa observação e 1 na outra não quer dizer que uma é duas vezes a outra, mas que uma é da categoria associada ao valor 2 e outra ao valor 1.

Explicando melhor, imagine um cenário em que a questão em uma entrevista pergunte a religião. Imagine que as possíveis respostas sejam: católico, presbiteriano, espírita ou muçulmano. Neste caso, o atributo religião possui quatro (4) possíveis respostas para este exemplo. Seria possível imaginar uma variável que assumiria valores 1, 2, 3 ou 4 conforme a religião respondida. Neste caso, seria uma variável multinomial, mas observe que o valor 4 (muçulmano) não quer dizer que é quatro vezes a resposta de uma observação com valor 1 (católico). Ainda, também não quer dizer que o 4 é maior que o 1. Apenas indica que um é do grupo 4 e outro do 1.

Neste caso, seria interessante trabalhar com variáveis binárias, normalmente associadas às respostas 0 e 1, no sentido de pertencer ou não a um grupo. Por exemplo, uma dummy (ou binária) para indicar católico (\(D_c\) = 1 se for e 0 se não for). Assim, similarmente para as demais respostas, têm-se outras três dummies binárias, uma para cada resposta (presbiteriano – \(D_p\), espírita – \(D_{es}\) ou muçulmano – \(D_{mu}\)).

Uma ideia seria fazer a regressão incluindo estas binárias na expressão:

\[ Y = {\alpha _0} + {\alpha _1}{X_1} + {\alpha _2}{X_2} + {\alpha _3}{D_1} + {\alpha _4}{D_1}{X_1} \\ + {\alpha _5}{D_2} + {\alpha _6}{D_2}{X_1} + {\alpha _7}{D_3} + {\alpha _8}{D_3}{X_1} \\ + {\alpha _9}{D_4} + {\alpha _{10}}{D_4}{X_1} + \varepsilon \]

ou

\[Y = {\alpha _0} + {\alpha _1}{X_1} + {\alpha _2}{X_2} + {\alpha _3}{D_c} + {\alpha _4}{D_c}{X_1} \\ + {\alpha _5}{D_p} + {\alpha _6}{D_p}{X_1} + {\alpha _7}{D_{es}} + {\alpha _8}{D_{es}}{X_1} \\ + {\alpha _9}{D_{mu}} + {\alpha _{10}}{D_{mu}}{X_1} + \varepsilon \]

O leitor poderia observar que em alguns termos a dummy aparece associada apenas ao coeficiente \(\alpha\) e noutros ao \(\alpha X_1\), ou seja, ao coeficiente e a variável \(X_1\). Imagine aqui que \(X_1\) seja uma variável quantitativa contínua. Em uma estimação cujo coeficiente \(\alpha\) fosse não significativo, seria o equivalente a ter um termo nulo, pois a dummy e a variável \(X\) multiplicada pelo coeficiente não significativo retornaria um valor não diferente de zero.

De outra forma, quando o \(\alpha\) for significativo para algum termo contendo a dummy, se está não estiver multiplicada por \(X_1\), seria o equivalente a uma alteração do intercepto \(\alpha_0\). Se estiver multiplicada por \(X_1\), seria o equivalente a uma alteração do coeficiente associado à inclinação de \(X_1\), somando-se o \(\alpha\) da dummy ao \(\alpha_1\). Neste caso, a presença da característica cujo coeficiente da dummy foi significativo, diferente de zero, indica uma curva estimada distinta das demais (no intercepto ou na inclinação ou em ambos).

Entretanto, uma estimação desta forma do exemplo, contendo quatro dummies para as quatro possíveis respostas, geraria um problema de multicolinearidade perfeita com o termo do intercepto. Recorde-se que na matriz de variáveis explicativas também está a coluna de 1 para o intercepto. Mas a soma das colunas das dummies também será 1. Ou seja, não será possível estimar o vetor de coeficientes, pois não será possível computar \((X’X)^{-1}\). Assim, o procedimento normal é excluir uma das dummies, ou seja, sempre que tiver \(m\) possíveis respostas, acrescentam-se \(m-1\) dummies. A dummy excluída do modelo será do grupo básico, e a estimação do modelo quando todos os coeficientes das dummies forem nulos retornará o resultado do grupo básico. Portanto, as observações do grupo básico terão valores zero para as dummies, por não fazer parte daqueles grupos.

Assim, para cada categoria, seria o equivalente a ter as equações como abaixo:

equação católico (\({D_c} = 1\)): \[Y = \left( {{\alpha _0} + {\alpha _3}} \right) + \left( {{\alpha _1} + {\alpha _4}} \right){X_1} + {\alpha _2}{X_2}\]

equação presbiteriano (\({D_p} = 1\)): \[Y = \left( {{\alpha _0} + {\alpha _5}} \right) + \left( {{\alpha _1} + {\alpha _6}} \right){X_1} + {\alpha _2}{X_2}\]

equação espírita (\({D_{es}} = 1\)): \[Y = \left( {{\alpha _0} + {\alpha _7}} \right) + \left( {{\alpha _1} + {\alpha _8}} \right){X_1} + {\alpha _2}{X_2}\]

equação muçulmano (\({D_c} = {D_p} = {D_{es}} = 0\))(BASE): \[Y = {\alpha _0} + {\alpha _1}{X_1} + {\alpha _2}{X_2}\]

Seja o exemplo abaixo, cujas observações respondem a uma pergunta indicando respostas de a até g. Imagine que os nomes sejam os setores de vendas em uma loja de varejo:

Nome	Resp
	a	b	c	d	e	f	g
Marcelo	x
Antônio		x
Tassiany	x
Elenir			x
Joice				x
Alexandre					x
Adriano						x
João							x
Pedro							x

Para representar estas respostas em variáveis binárias, no exemplo, adota-se o g como o grupo básico. Se forem feitas uma dummy para cada resposta, o quadro ficaria, colocando 1 para pertence e zero para não pertence:

Observação	\(D_a\)	\(D_b\)	\(D_c\)	\(D_d\)	\(D_e\)	\(D_f\)	\(D_g\)	Soma
Marcelo	1	0	0	0	0	0	0	1
Antônio	0	1	0	0	0	0	0	1
Tassiany	1	0	0	0	0	0	0	1
Elenir	0	0	1	0	0	0	0	1
Joice	0	0	0	1	0	0	0	1
Alexandre	0	0	0	0	1	0	0	1
Adriano	0	0	0	0	0	1	0	1
João	0	0	0	0	0	0	1	1
Pedro	0	0	0	0	0	0	1	1

Observe que a soma da linha é sempre igual a 1. Assim, uma coluna adicional para representar o intercepto seria exatamente igual a soma das colunas de dummies. A solução seria excluir a coluna de \(D_g\), e todas as observações do grupo g ficariam com zeros nas dummies de a até f (João e Pedro).

Observação	\(D_a\)	\(D_b\)	\(D_c\)	\(D_d\)	\(D_e\)	\(D_f\)
Marcelo	1	0	0	0	0	0
Antônio	0	1	0	0	0	0
Tassiany	1	0	0	0	0	0
Elenir	0	0	1	0	0	0
Joice	0	0	0	1	0	0
Alexandre	0	0	0	0	1	0
Adriano	0	0	0	0	0	1
João	0	0	0	0	0	0
Pedro	0	0	0	0	0	0

Estimando um modelo do tipo acima, alterando apenas intercepto, cada grupo teria a equação estimada como abaixo, ou seja, alterando o intercepto relativamente ao grupo básico g:

\[Y = {\alpha _0} + {\alpha _1}{X_1} + {\alpha _2}{X_2} + {\alpha _3}{D_a} + {\alpha _4}{D_b} \\ + {\alpha _5}{D_c} + {\alpha _6}{D_d} + {\alpha _7}{D_e} + {\alpha _8}{D_f} + \varepsilon \]

Quando as dummies forem todas com coeficientes nulos, \(D_a=D_b=D_c=D_d=D_e=D_f=0\), então o modelo se reduz a curva para o grupo g, quando as dummies são nulas:

\[Y = {\alpha _0} + {\alpha _1}{X_1} + {\alpha _2}{X_2} + \varepsilon \]

Assim, os coeficientes \({\alpha _0},{\alpha _1},{\alpha _2}\) representam exatamente o grupo básico g.

Suponha que os coeficientes estimados associados as dummies, sejam significativos, ou \({\alpha _d} \ne0\), d= 3, 4, 5, 6, 7 e 8. Então a expressão para os grupos serão:

para a: \({Y_a} = \left( {{{\hat \alpha }_0} + {{\hat \alpha }_3}} \right) + {{\hat \alpha }_1}{X_1} + {{\hat \alpha }_2}{X_2}\)

para b: \({Y_b} = \left( {{{\hat \alpha }_0} + {{\hat \alpha }_4}} \right) + {{\hat \alpha }_1}{X_1} + {{\hat \alpha }_2}{X_2}\)

para c: \({Y_c} = \left( {{{\hat \alpha }_0} + {{\hat \alpha }_5}} \right) + {{\hat \alpha }_1}{X_1} + {{\hat \alpha }_2}{X_2}\)

para d: \({Y_d} = \left( {{{\hat \alpha }_0} + {{\hat \alpha }_6}} \right) + {{\hat \alpha }_1}{X_1} + {{\hat \alpha }_2}{X_2}\)

para e: \({Y_e} = \left( {{{\hat \alpha }_0} + {{\hat \alpha }_7}} \right) + {{\hat \alpha }_1}{X_1} + {{\hat \alpha }_2}{X_2}\)

para f: \({Y_f} = \left( {{{\hat \alpha }_0} + {{\hat \alpha }_8}} \right) + {{\hat \alpha }_1}{X_1} + {{\hat \alpha }_2}{X_2}\)

Os valores de \({\alpha _d} \ne0\), d= 3, 4, 5, 6, 7 e 8, podem ser positivos ou negativos. Quando positivos, indicarão interceptos maiores para o grupo, comparativamente ao grupo básico g. Valores negativos indicarão interceptos menores do que o de g. Similarmente se teria para as inclinações, indicando reações maiores ou menores de \(Y\) com respeito a \(X\).

• Exemplo producao_de_algodão:

Figura 4.2 QR Code para o exemplo da produção de algodão nos estados. X.

Seja uma estimação da quantidade de algodão (\(Qalg\)) em função da área de algodão (\(Aalg\)), do preço do algodão (\(Palg\)), do preço da soja (\(Psoja\)), para municípios de Mato Grosso do Sul, Mato Grosso e Goiás. É possível dizer que a quantidade de algodão tem um comportamento distinto entre os estados MS, MT e GO?

A rotina para chamar os dados de uma arquivo Excel encontra-se a seguir, porém, optamos por colocá-los dentro do arquivo RMarkdown. A rotina do embeded data está logo na sequência.

# chamando os dados de producao de algodao dos estados do
# centro-oeste com dados da PAM IBGE 1990 a 2000, para 
# cada estado MS MT GO ou seja 33 observacoes
library(readxl); library(foreign)
dados <- read_excel("dados_centro_oeste.xlsx", 
                    sheet = "dados")
# Qalg = quantidade produzida de algodao
# Aalg = área total colhida de algodao 
# Palg = preco do algodao
# Psoja = preco da soja 
# DMS e DMT sao dummies para os estados em que GO é a base
# estado é uma variável indicadora contendo a sigla do estado

A opção gerada pelo dput() dos dados básicos está abaixo:

dados<-structure(list(ano = c(1990, 1991, 1992, 1993, 1994, 1995, 1996, 
1997, 1998, 1999, 2000, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 
1997, 1998, 1999, 2000, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 
1997, 1998, 1999, 2000), Qalg = c(59754, 83650, 83710, 94560, 
101368, 157031, 173796, 189699, 260452, 278363, 254476, 73559, 
90561, 85119, 64735, 77409, 105791, 87952, 56027, 93229, 114521, 
127839, 57634, 73458, 67862, 85641, 91828, 87458, 73553, 78376, 
271038, 630406, 1002836), Aalg = c(35459, 42990, 53772, 38172, 
53773, 69533, 81659, 83234, 186661, 117056, 96718, 44570, 51888, 
73333, 39643, 41135, 60011, 59637, 26604, 49151, 46231, 48450, 
43422, 68443, 53836, 69584, 66059, 69390, 55075, 42259, 106483, 
200182, 257762), Palg = c(0.356572484593691, 0.28548486348212, 
0.18460421844633, 0.206416737228048, 0.415961827791941, 0.306940214830232, 
0.355443864081562, 0.426290579971575, 0.345460690410503, 0.369425952314763, 
0.372826913527653, 0.370110836978266, 0.337634829251245, 0.19012070963983, 
0.18883532150463, 0.514823904914612, 0.340196812424667, 0.303950608442499, 
0.399129289604689, 0.340215816862215, 0.407992489922868, 0.362518549919648, 
0.268432114188467, 0.248378441722003, 0.171583759997412, 0.214788618467748, 
0.511497661369565, 0.304906963153987, 0.28341715278882, 0.359730772000928, 
0.355262834733223, 0.369540331161043, 0.326479031442301), 
Psoja = c(0.129725363498433, 
0.116121616325159, 0.112620111592226, 0.0873969828327483, 0.176972146441714, 
0.103769694520009, 0.149613452870967, 0.1663861049096, 0.143011987528068, 
0.14175993084962, 0.137533111256611, 0.150873051860415, 0.138236213975033, 
0.12167602921687, 0.129874022775396, 0.207511199152738, 0.108008706367883, 
0.154946820778127, 0.172922997949505, 0.138872367383235, 0.148515655529693, 
0.136549826243949, 0.113535605629564, 0.117219268828656, 0.0968498706796199, 
0.0719258443192821, 0.208422118749737, 0.104896138231911, 0.135136436052929, 
0.150488044032018, 0.123091268671333, 0.11957295449546, 0.130580697047997
), DMS = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), DMT = c(0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1), estado = c("go", "go", "go", "go", 
"go", "go", "go", "go", "go", "go", "go", "ms", "ms", "ms", "ms", 
"ms", "ms", "ms", "ms", "ms", "ms", "ms", "mt", "mt", "mt", "mt", 
"mt", "mt", "mt", "mt", "mt", "mt", "mt")), row.names = c(NA, 
-33L), class = c("tbl_df", "tbl", "data.frame"))
attach(dados)

Vamos ver como está a tabela importada:

library(knitr)
kable(dados,caption="Dados do Exercício")

Table 4.1: Dados do Exercício

ano	Qalg	Aalg	Palg	Psoja	DMS	DMT	estado
1990	59754	35459	0.3565725	0.1297254	0	0	go
1991	83650	42990	0.2854849	0.1161216	0	0	go
1992	83710	53772	0.1846042	0.1126201	0	0	go
1993	94560	38172	0.2064167	0.0873970	0	0	go
1994	101368	53773	0.4159618	0.1769721	0	0	go
1995	157031	69533	0.3069402	0.1037697	0	0	go
1996	173796	81659	0.3554439	0.1496135	0	0	go
1997	189699	83234	0.4262906	0.1663861	0	0	go
1998	260452	186661	0.3454607	0.1430120	0	0	go
1999	278363	117056	0.3694260	0.1417599	0	0	go
2000	254476	96718	0.3728269	0.1375331	0	0	go
1990	73559	44570	0.3701108	0.1508731	1	0	ms
1991	90561	51888	0.3376348	0.1382362	1	0	ms
1992	85119	73333	0.1901207	0.1216760	1	0	ms
1993	64735	39643	0.1888353	0.1298740	1	0	ms
1994	77409	41135	0.5148239	0.2075112	1	0	ms
1995	105791	60011	0.3401968	0.1080087	1	0	ms
1996	87952	59637	0.3039506	0.1549468	1	0	ms
1997	56027	26604	0.3991293	0.1729230	1	0	ms
1998	93229	49151	0.3402158	0.1388724	1	0	ms
1999	114521	46231	0.4079925	0.1485157	1	0	ms
2000	127839	48450	0.3625185	0.1365498	1	0	ms
1990	57634	43422	0.2684321	0.1135356	0	1	mt
1991	73458	68443	0.2483784	0.1172193	0	1	mt
1992	67862	53836	0.1715838	0.0968499	0	1	mt
1993	85641	69584	0.2147886	0.0719258	0	1	mt
1994	91828	66059	0.5114977	0.2084221	0	1	mt
1995	87458	69390	0.3049070	0.1048961	0	1	mt
1996	73553	55075	0.2834172	0.1351364	0	1	mt
1997	78376	42259	0.3597308	0.1504880	0	1	mt
1998	271038	106483	0.3552628	0.1230913	0	1	mt
1999	630406	200182	0.3695403	0.1195730	0	1	mt
2000	1002836	257762	0.3264790	0.1305807	0	1	mt

Neste exemplo, para \(m=3\) respostas, colocam-se \(m-1=2\) dummies: \(DMS\) e \(DMT\). Goiás ficou como grupo básico, ou seja, para \(DMS=DMT=0\). Conforme a estimação do Quadro 4.13, as dummies foram inseridas alterando o intercepto, com valores das probabilidades, respectivamente: 0,98 e 0,09, DMS e DMT. Assim, não foram coeficientes significativos diferentes de zero para \(log(Psoja)\) e \(DMS\), e pode-se dizer que não houve diferença entre os interceptos de MS e o grupo básico de GO. Mas pode-se dizer que existe diferença do intercepto de MT em relação a GO pois o valor da probabilidade de 0,09 indica que a \(DMT\) teve coeficiente significativo diferente de zero para MT.

attach(dados)
reg.empilhada<-lm(log(Qalg)~log(Aalg)+log(Palg)+log(Psoja))
summary(reg.empilhada)

## 
## Call:
## lm(formula = log(Qalg) ~ log(Aalg) + log(Palg) + log(Psoja))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46795 -0.14723  0.03913  0.14962  0.51140 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.08883    1.10940  -0.981   0.3345    
## log(Aalg)    1.15505    0.09816  11.767 1.46e-12 ***
## log(Palg)    0.63658    0.26442   2.407   0.0227 *  
## log(Psoja)  -0.35201    0.33259  -1.058   0.2986    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2711 on 29 degrees of freedom
## Multiple R-squared:  0.854,  Adjusted R-squared:  0.8389 
## F-statistic: 56.56 on 3 and 29 DF,  p-value: 3.118e-12

# incluindo Akaike e Schwarz no objeto reg.empilhada
reg.empilhada$AIC<- AIC(reg.empilhada)
reg.empilhada$BIC<- BIC(reg.empilhada)

reg.dint<-lm(log(Qalg)~log(Aalg)+log(Palg)+log(Psoja)+DMS+DMT, 
             data=dados)
summary(reg.dint)

## 
## Call:
## lm(formula = log(Qalg) ~ log(Aalg) + log(Palg) + log(Psoja) + 
##     DMS + DMT, data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58105 -0.13666 -0.00039  0.19076  0.57450 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.800698   1.251232  -1.439   0.1616    
## log(Aalg)    1.204696   0.103819  11.604 5.33e-12 ***
## log(Palg)    0.653454   0.258398   2.529   0.0176 *  
## log(Psoja)  -0.475624   0.336745  -1.412   0.1693    
## DMS         -0.002898   0.121463  -0.024   0.9811    
## DMT         -0.202506   0.114617  -1.767   0.0886 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2637 on 27 degrees of freedom
## Multiple R-squared:  0.8714, Adjusted R-squared:  0.8476 
## F-statistic: 36.58 on 5 and 27 DF,  p-value: 3.33e-11

reg.dint$AIC<- AIC(reg.dint)
reg.dint$BIC<- BIC(reg.dint)

stargazer::stargazer(list(reg.empilhada,reg.dint),
          title="Título: Resultados das Regressões",
          align=TRUE,
          type = "text", 
          style="commadefault",
          keep.stat=c("aic","bic","rsq", "adj.rsq","n")
          )

## 
## Título: Resultados das Regressões
## ================================================
##                         Dependent variable:     
##                     ----------------------------
##                              log(Qalg)          
##                          (1)            (2)     
## ------------------------------------------------
## log(Aalg)              1,155***      1,205***   
##                        (0,098)        (0,104)   
##                                                 
## log(Palg)              0,637**        0,653**   
##                        (0,264)        (0,258)   
##                                                 
## log(Psoja)              -0,352        -0,476    
##                        (0,333)        (0,337)   
##                                                 
## DMS                                   -0,003    
##                                       (0,121)   
##                                                 
## DMT                                   -0,203*   
##                                       (0,115)   
##                                                 
## Constant                -1,089        -1,801    
##                        (1,109)        (1,251)   
##                                                 
## ------------------------------------------------
## Observations              33            33      
## R2                      0,854          0,871    
## Adjusted R2             0,839          0,848    
## Akaike Inf. Crit.       13,233        13,061    
## Bayesian Inf. Crit.     20,716        23,536    
## ================================================
## Note:                *p<0,1; **p<0,05; ***p<0,01

A representação das equações para cada grupo serão, partindo da equação geral:

Log(Qalg) = 1.205* Log(Aalg) + 0.653* Log(Palg) - 0.476* Log(Psoja) - 0.003* DMS - 0.203*DMT - 1.801

Para GO, DMS=DMT=0, com PSOJA e intercepto C não significativos:

Log(Qalg) = 1.205* Log(Aalg) + 0.653*Log(Palg)

Para MS: DMT=0, mas DMS não teve coeficiente significativo:

Log(Qalg) = 1.205* Log(Aalg) + 0.653*Log(Palg)

Para MT: DMS=0, ou seja, esta curva está pouco abaixo das curvas de MS e GO:

Log(Qalg) = 1.205* Log(Aalg) + 0.653*Log(Palg) - 0.203

Para uma estimação contendo alterações nas inclinações de \(Aalg\) e no intercepto, faremos uso da variável estado e das opções com DMS e DMT. A variável estado é obtida pela função dummy do pacote dummies:

# o pacote cria estadogo, estadoms,estadomt
# 
library(dummies)
dummies<-dummy(estado, drop=FALSE)
colnames(dummies)<-c("estadogo", "estadoms","estadomt")
dados<-cbind(dados,dummies)
# o pacote cria dummies usando o nome da variável 
# justaposto com o codigo das observacoes: 
# estadogo, estadoms,estadomt
#print(dados)

reg.dum<-lm(log(Qalg)~log(Aalg)+log(Palg)+log(Psoja)+
              estado+I(log(Aalg)*DMS)+I(log(Aalg)*DMT), 
            data=dados)
reg.dum$AIC<- AIC(reg.dum)
reg.dum$BIC<- BIC(reg.dum)
#summary(reg.dum)
reg.dum2<-lm(log(Qalg)~log(Aalg)+log(Palg)+log(Psoja)+
               DMS+DMT+I(log(Aalg)*DMS)+I(log(Aalg)*DMT), 
             data=dados)
reg.dum2$AIC<- AIC(reg.dum2)
reg.dum2$BIC<- BIC(reg.dum2)
stargazer::stargazer(list(reg.dum,reg.dum2),
          title="Título: Resultados das Regressões",
          align=TRUE,
          type = "text", 
          style="commadefault",
          keep.stat=c("aic","bic","rsq", "adj.rsq","n")
          )

## 
## Título: Resultados das Regressões
## ================================================
##                         Dependent variable:     
##                     ----------------------------
##                              log(Qalg)          
##                          (1)            (2)     
## ------------------------------------------------
## log(Aalg)              0,926***      0,926***   
##                        (0,134)        (0,134)   
##                                                 
## log(Palg)              0,499**        0,499**   
##                        (0,207)        (0,207)   
##                                                 
## log(Psoja)              -0,339        -0,339    
##                        (0,272)        (0,272)   
##                                                 
## estadoms                2,855                   
##                        (3,272)                  
##                                                 
## estadomt              -7,201***                 
##                        (1,952)                  
##                                                 
## DMS                                    2,855    
##                                       (3,272)   
##                                                 
## DMT                                  -7,201***  
##                                       (1,952)   
##                                                 
## I(log(Aalg) * DMS)      -0,276        -0,276    
##                        (0,301)        (0,301)   
##                                                 
## I(log(Aalg) * DMT)     0,624***      0,624***   
##                        (0,174)        (0,174)   
##                                                 
## Constant                1,408          1,408    
##                        (1,625)        (1,625)   
##                                                 
## ------------------------------------------------
## Observations              33            33      
## R2                      0,926          0,926    
## Adjusted R2             0,906          0,906    
## Akaike Inf. Crit.       -1,353        -1,353    
## Bayesian Inf. Crit.     12,115        12,115    
## ================================================
## Note:                *p<0,1; **p<0,05; ***p<0,01

Observe que neste caso, não foram significativos os coeficientes de DMS*Log(Aalg), Log(Psoja), DMS e Constante. Assim, as expressões serão, omitindo os termos não significativos:

Log(Qalg) = 0.926* Log(Aalg) + 0.624* DMT* LOG(AALG) + 0.499* Log(Palg) - 7.201*DMT

Para GO: Log(Qalg) = 0.926* Log(Aalg) + 0.624* DMT* Log(Aalg) + 0.499* Log(Palg)
Para MS: Log(Qalg) = 0.926* Log(Aalg) + 0.624* DMT* Log(Aalg) + 0.499* Log(Palg)
Para MT: Log(Qalg) = (0.926+ 0.624)* Log(Aalg) + 0.499* Log(Palg) - 7.201

Assim, conclui-se que a reação de QALG a AALG para MT é maior que em GO e MS, mas a curva de MT está em nível mais baixo que os demais.

# sejam as variáveis v para a serie que se quer 
# a TGC e t para o ano

# para Qalg MT
vmt<-dados[23:33,2]
vmt

##  [1]   57634   73458   67862   85641   91828   87458   73553   78376  271038
## [10]  630406 1002836

t<-dados[23:33,1]   # anos 1990 a 2000
logv<-log(vmt)
eq=lm(logv~t)
summary(eq)

## 
## Call:
## lm(formula = logv ~ t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99807 -0.25869  0.05916  0.38977  0.82456 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -471.29928  113.40011  -4.156  0.00246 **
## t              0.24215    0.05684   4.260  0.00211 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5962 on 9 degrees of freedom
## Multiple R-squared:  0.6685, Adjusted R-squared:  0.6316 
## F-statistic: 18.15 on 1 and 9 DF,  p-value: 0.002111

# pego o segundo parâmetro - o primeiro é o intercepto
beta <- summary(eq)$coefficients[2]  
# Calculo da TGC
TGC=I((exp(beta)-1)*100)
TGC

## [1] 27.39809

tgc.dados<-structure(
  list(v = c(100, 112, 145, 178, 134, 104, 108, 112, 110), 
       t2 = c(91, 92, 93, 94, 95, 96, 97, 98, 99)), 
  row.names = c(NA, -9L), 
  class = c("tbl_df", "tbl", "data.frame"))
attach(tgc.dados)
logv<-log(v)
eq<- lm(logv~t2)
summary(eq)

## 
## Call:
## lm(formula = logv ~ t2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.23625 -0.11049 -0.04155  0.10611  0.37764 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  5.97187    2.43596   2.452    0.044 *
## t2          -0.01242    0.02563  -0.485    0.643  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1985 on 7 degrees of freedom
## Multiple R-squared:  0.03247,    Adjusted R-squared:  -0.1058 
## F-statistic: 0.2349 on 1 and 7 DF,  p-value: 0.6427

beta <- summary(eq)$coefficients[2]
(TGC=I((exp(beta)-1)*100))

## [1] -1.234577

#[1] -1.234577

4.4.3 Identificação de outliers do modelo

Os valores atípicos ou outliers são observações com valores muito diferentes da maioria, ou da média. Estes valores poderão alterar os valores estimados dos coeficientes e invalidar os resultados obtidos caso não sejam tratados adequadamente.

Em geral se utilizam gráficos tipo scatter (dispersão) ou box-plot associados aos valores do desvio-padrão e seus múltiplos como critérios para julgar se os valores estão muito distintos da média. Por exemplo, na Figura 4.4, observam-se alguns pontos muito destacados dos demais, ou valores atípicos. Tratam-se de observações dos setores censitários de Sinop-MT, Brasil, a partir do Censo Agropecuário de 2006. Maiores detalhes podem ser encontrados em http://rpubs.com/amrofi/exercicio_sinop.

• Exemplo sinop:

Figura 4.3 QR Code para o exemplo da produção de soja nos setores censitários de Sinop-MT em 2006. X.

O exercício traz dados dos 59 setores censitários de Sinop-MT, a partir do Censo Agropecuário de 2006. As variaveis foram logaritmizadas no excel. VP = valor total da produção; AREA = terra (medida através da área total colhida); QTMAQ = capital (medido através da quantidade total de maquinários), QMO = trabalho (medido através da mão de obra ocupada por pessoas com mais de 14 anos), PCALC = tecnologia (medida através da proporção do uso de calcário); PADUB = tecnologia (medida através da proporção do uso de adubo); PAGROT = tecnologia (medida através da proporção do uso de agrotóxico); PNIVEL = tecnologia (medida através da proporção do uso de curva de nível); PPRAGA = tecnologia (medida através da proporção do uso de controle geral de pragas); e, PROT = tecnologia (medida através da proporção do uso de rotação de culturas).

Gráfico de dispersão para identificação de outliers.

#Buscando os outliers
dispersao<-plot(PADUB~VP)

Os gráficos de box-plot mostram o centro e a dispersão dos dados. O box refere-se aos limites dados pelo primeiro e terceiro quartis. O ponto ao meio, dentro do box, indica a média. A mediana é o traço horizontal no box. Os traços horizontais denotam os inner fences e o staple: limites a partir dos quais são detectados outliers (ver Figura 4.4).

Figura 4 4. Detalhamento dos limites do box-plot.

Figura 4.5 Gráfico de box-plot para identificação de outliers.

Quando se tem dados com outliers, uma alternativa é distinguir grupos utilizando as dummies. Esta situação poderá auxiliar a melhor incorporar esta informação no modelo sem, no entanto excluir a observação. Esta opção de excluir a observação atípica poderá gerar problemas de viés de seleção da amostra.

No R, por meio do comando bp1_vp<-boxplot(VP) tem-se:

bp1_vp<-boxplot(VP, names=c("VP"))

# retorna os valores de VP outliers
(bp1_vp_out<-bp1_vp$out)  

## [1] 105352999  69309029  79804101 110434855  66886216  52673412  58868284
## [8]  73931129

Para as demais variáveis seria com o comando abaixo e gerando a Figura 4.7.

bp2_proportions<-boxplot(PCALC,PADUB,
      PAGROT,PNIVEL,PPRAGA,PROT,horizontal = TRUE,
      names = c("PCALC","PADUB","PAGROT","PNIVEL","PPRAGA","PROT"),
      main = "Figura 4 7. Box-plot das variáveis de proporções 
      de práticas tecnológicas PCALC, PADUB, PAGROT, 
      PNIVEL, PPRAGA, PROT.")

Um boxplot que julgo mais útil é por meio do pacote car, que reporta qual a observação atípica.

library(car)
Boxplot(VP,label=dados$obs)

## [1]  6  7 34 35 37 38 44 45

Vamos utilizar o pacote stargazer para organizar as saídas de resultados. Se a saída fosse apenas pelo comando summary, sairia da forma:

regressao1<-lm(log(VP)~log(AREA)+log(QMO)+log(PPRAGA)+log(PROT))
regressao1$AIC <- AIC(regressao1)
regressao1$BIC <- BIC(regressao1)
suppressMessages(library(stargazer))
stargazer(regressao1, 
          title = "Título: Resultado da Regressão OLS", 
          align = TRUE, type = "text", 
          style = "all", 
          keep.stat = c("AIC", "BIC", "rsq", "adj.rsq", "n"))

## 
## Título: Resultado da Regressão OLS
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               log(VP)          
## -----------------------------------------------
## log(AREA)                    1.041***          
##                               (0.025)          
##                             t = 41.000         
##                              p = 0.000         
## log(QMO)                      -0.020           
##                               (0.039)          
##                             t = -0.514         
##                              p = 0.610         
## log(PPRAGA)                   -0.006           
##                               (0.007)          
##                             t = -0.904         
##                              p = 0.370         
## log(PROT)                    -0.018**          
##                               (0.008)          
##                             t = -2.197         
##                              p = 0.033         
## Constant                     6.496***          
##                               (0.279)          
##                             t = 23.279         
##                              p = 0.000         
## -----------------------------------------------
## Observations                    59             
## R2                             0.971           
## Adjusted R2                    0.969           
## Akaike Inf. Crit.             36.063           
## Bayesian Inf. Crit.           48.528           
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

summary(regressao1)

## 
## Call:
## lm(formula = log(VP) ~ log(AREA) + log(QMO) + log(PPRAGA) + log(PROT))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35851 -0.15577 -0.03277  0.08441  1.80533 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.495877   0.279044  23.279   <2e-16 ***
## log(AREA)    1.041318   0.025398  41.000   <2e-16 ***
## log(QMO)    -0.020008   0.038909  -0.514   0.6092    
## log(PPRAGA) -0.006484   0.007170  -0.904   0.3698    
## log(PROT)   -0.018056   0.008217  -2.197   0.0323 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3101 on 54 degrees of freedom
## Multiple R-squared:  0.9713, Adjusted R-squared:  0.9691 
## F-statistic: 456.2 on 4 and 54 DF,  p-value: < 2.2e-16

Outra alternativa é observar os gráficos dos diagnósticos da regressão, em que se podem ver os outliers 51 e 54 na Figura 4.9.

split.screen(c(2,2))

## [1] 1 2 3 4

screen(1)
reg1.plot1<-plot(regressao1, which = 1)

screen(2)
reg1.plot2<-plot(regressao1, which = 2)

screen(3)
reg1.plot3<-plot(regressao1, which = 3) 

screen(4)
reg1.plot4<-plot(regressao1, which = 4)

close.screen(all = TRUE)    # exit split-screen mode

Estes outliers (observações 51 e 54) são outliers da regressão e não outliers da série de dados.

Outra opção para investigar outliers é o teste de Bonferroni do pacote ‘car’ do R, em que tem-se H0: a observação i não difere das demais observações, e portanto, i não é um outlier. Se o P-value foi menor que 0.10, rejeito H0 e tem-se um outlier! O ideal é retirar ou “controlar” a observação (por exemplo, por meio de uma dummy) e depois refazer o teste até não ter mais outliers (quando a probabilidade será detectada maior que 0.10).O teste acusa apenas a observação 54 como outlier.

library(car)
outlierTest(regressao1)

##    rstudent unadjusted p-value Bonferroni p
## 54 9.789123         1.7633e-13   1.0404e-11

Conforme o box-plot da variável VPROD, excluindo-se os outliers 1 e 2, tem-se a Figura 4.10. Para retirá-los, criaremos um novo dataset e farei a nova regressão.

# criando dataset excluindo as observacoes 51 e 54
dados2<-dados[dados$D51_54==0, ]
attach(dados2)
#View(dados2)
regressao2<-lm(log(VP)~log(AREA)+log(QMO)+
                 log(PPRAGA)+log(PROT),data=dados2)
results<-stargazer::stargazer(
  list(regressao1,regressao2),
  align = TRUE, type = "text", 
  style = "all", 
  keep.stat = c("AIC", "BIC", "rsq", "adj.rsq", "n"))

## 
## ================================================
##                         Dependent variable:     
##                     ----------------------------
##                               log(VP)           
##                          (1)            (2)     
## ------------------------------------------------
## log(AREA)              1.041***      1.052***   
##                        (0.025)        (0.013)   
##                       t = 41.000    t = 80.874  
##                       p = 0.000      p = 0.000  
## log(QMO)                -0.020        -0.038*   
##                        (0.039)        (0.020)   
##                       t = -0.514    t = -1.858  
##                       p = 0.610      p = 0.069  
## log(PPRAGA)             -0.006       -0.010**   
##                        (0.007)        (0.004)   
##                       t = -0.904    t = -2.655  
##                       p = 0.370      p = 0.011  
## log(PROT)              -0.018**      -0.019***  
##                        (0.008)        (0.004)   
##                       t = -2.197    t = -4.475  
##                       p = 0.033     p = 0.00005 
## Constant               6.496***      6.423***   
##                        (0.279)        (0.144)   
##                       t = 23.279    t = 44.591  
##                       p = 0.000      p = 0.000  
## ------------------------------------------------
## Observations              59            57      
## R2                      0.971          0.993    
## Adjusted R2             0.969          0.992    
## Akaike Inf. Crit.       36.063                  
## Bayesian Inf. Crit.     48.528                  
## ================================================
## Note:                *p<0.1; **p<0.05; ***p<0.01

# montando os vários gráficos de diagnosticos da regressao
# rodando novo teste de Bonferroni para outlier 
# na regressao sem obs 51 e 54
car::outlierTest(regressao2)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 8 -2.375175           0.021338           NA

u.hat2<-resid(regressao2)
t.test(u.hat2)

## 
##  One Sample t-test
## 
## data:  u.hat2
## t = 4.4084e-17, df = 56, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.04051635  0.04051635
## sample estimates:
##    mean of x 
## 8.916136e-19

hist.reg2<- hist(u.hat2, freq = FALSE)
curve(dnorm, add = TRUE,col="red")

split.screen(c(2,2))

## [1] 1 2 3 4

screen(1)
reg2.plot1<-plot(regressao2, which = 1)

screen(2)
reg2.plot2<-plot(regressao2, which = 2)

screen(3)
reg2.plot3<-plot(regressao2, which = 3) 

screen(4)
reg2.plot4<-plot(regressao2, which = 4)

close.screen(all = TRUE)    # exit split-screen mode

• Opção com dummies para os outliers

Neste caso, optaremos por manter o dataset com todas as observações e controlaremos os efeitos para os outliers 51 e 54. Observe que agora o modelo fica melhor que sem as dummies para outliers.

mod51_54<-lm(log(VP)~log(AREA)+log(QMO)+
               log(PPRAGA)+log(PROT)+D51_54, 
             data = dados)
summary(mod51_54)

## 
## Call:
## lm(formula = log(VP) ~ log(AREA) + log(QMO) + log(PPRAGA) + log(PROT) + 
##     D51_54, data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53609 -0.12275  0.02064  0.11632  0.53609 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.521749   0.170251  38.307  < 2e-16 ***
## log(AREA)    1.049127   0.015515  67.619  < 2e-16 ***
## log(QMO)    -0.052553   0.023977  -2.192 0.032811 *  
## log(PPRAGA) -0.009423   0.004385  -2.149 0.036211 *  
## log(PROT)   -0.018117   0.005013  -3.614 0.000671 ***
## D51_54       1.323510   0.137910   9.597 3.48e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1892 on 53 degrees of freedom
## Multiple R-squared:  0.9895, Adjusted R-squared:  0.9885 
## F-statistic: 999.2 on 5 and 53 DF,  p-value: < 2.2e-16

mod51_54$AIC <- AIC(mod51_54)
mod51_54$BIC <- BIC(mod51_54)
#suppressMessages(library(stargazer))
stargazer(list(regressao1,mod51_54), 
          title = "Título: Resultado da Regressão OLS", 
          align = TRUE, type = "text", 
          style = "all", 
          keep.stat = c("AIC", "BIC", "rsq", "adj.rsq", "n"))

## 
## Título: Resultado da Regressão OLS
## ================================================
##                         Dependent variable:     
##                     ----------------------------
##                               log(VP)           
##                          (1)            (2)     
## ------------------------------------------------
## log(AREA)              1.041***      1.049***   
##                        (0.025)        (0.016)   
##                       t = 41.000    t = 67.619  
##                       p = 0.000      p = 0.000  
## log(QMO)                -0.020       -0.053**   
##                        (0.039)        (0.024)   
##                       t = -0.514    t = -2.192  
##                       p = 0.610      p = 0.033  
## log(PPRAGA)             -0.006       -0.009**   
##                        (0.007)        (0.004)   
##                       t = -0.904    t = -2.149  
##                       p = 0.370      p = 0.037  
## log(PROT)              -0.018**      -0.018***  
##                        (0.008)        (0.005)   
##                       t = -2.197    t = -3.614  
##                       p = 0.033      p = 0.001  
## D51_54                               1.324***   
##                                       (0.138)   
##                                      t = 9.597  
##                                      p = 0.000  
## Constant               6.496***      6.522***   
##                        (0.279)        (0.170)   
##                       t = 23.279    t = 38.307  
##                       p = 0.000      p = 0.000  
## ------------------------------------------------
## Observations              59            59      
## R2                      0.971          0.990    
## Adjusted R2             0.969          0.989    
## Akaike Inf. Crit.       36.063        -21.358   
## Bayesian Inf. Crit.     48.528        -6.815    
## ================================================
## Note:                *p<0.1; **p<0.05; ***p<0.01

• Alterando interceptos e inclinações da AREA para 51 e 54

mod51_54incl<-lm(log(VP)~log(AREA)+log(QMO)+log(PPRAGA)+
                   log(PROT)+D51_54+I(log(AREA)*D51_54), 
                 data = dados)
summary(mod51_54incl)

## 
## Call:
## lm(formula = log(VP) ~ log(AREA) + log(QMO) + log(PPRAGA) + log(PROT) + 
##     D51_54 + I(log(AREA) * D51_54), data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30924 -0.11933  0.00000  0.09946  0.28607 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            6.422928   0.144040  44.591  < 2e-16 ***
## log(AREA)              1.052109   0.013009  80.874  < 2e-16 ***
## log(QMO)              -0.037739   0.020313  -1.858   0.0688 .  
## log(PPRAGA)           -0.009752   0.003673  -2.655   0.0105 *  
## log(PROT)             -0.018799   0.004201  -4.475 4.19e-05 ***
## D51_54                 9.003078   1.586565   5.675 6.26e-07 ***
## I(log(AREA) * D51_54) -0.862481   0.177712  -4.853 1.15e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1585 on 52 degrees of freedom
## Multiple R-squared:  0.9928, Adjusted R-squared:  0.9919 
## F-statistic:  1191 on 6 and 52 DF,  p-value: < 2.2e-16

mod51_54incl$AIC <- AIC(mod51_54incl)
mod51_54incl$BIC <- BIC(mod51_54incl)
#suppressMessages(library(stargazer))
stargazer(list(mod51_54incl,mod51_54), 
          title = "Título: Resultado da Regressão OLS", 
          align = TRUE, type = "text", 
          style = "all", 
          keep.stat = c("AIC", "BIC", "rsq", "adj.rsq", "n"))

## 
## Título: Resultado da Regressão OLS
## ==================================================
##                           Dependent variable:     
##                       ----------------------------
##                                 log(VP)           
##                            (1)            (2)     
## --------------------------------------------------
## log(AREA)                1.052***      1.049***   
##                          (0.013)        (0.016)   
##                         t = 80.874    t = 67.619  
##                         p = 0.000      p = 0.000  
## log(QMO)                 -0.038*       -0.053**   
##                          (0.020)        (0.024)   
##                         t = -1.858    t = -2.192  
##                         p = 0.069      p = 0.033  
## log(PPRAGA)              -0.010**      -0.009**   
##                          (0.004)        (0.004)   
##                         t = -2.655    t = -2.149  
##                         p = 0.011      p = 0.037  
## log(PROT)               -0.019***      -0.018***  
##                          (0.004)        (0.005)   
##                         t = -4.475    t = -3.614  
##                        p = 0.00005     p = 0.001  
## D51_54                   9.003***      1.324***   
##                          (1.587)        (0.138)   
##                         t = 5.675      t = 9.597  
##                        p = 0.00000     p = 0.000  
## I(log(AREA) * D51_54)   -0.862***                 
##                          (0.178)                  
##                         t = -4.853                
##                        p = 0.00002                
## Constant                 6.423***      6.522***   
##                          (0.144)        (0.170)   
##                         t = 44.591    t = 38.307  
##                         p = 0.000      p = 0.000  
## --------------------------------------------------
## Observations                59            59      
## R2                        0.993          0.990    
## Adjusted R2               0.992          0.989    
## Akaike Inf. Crit.        -41.400        -21.358   
## Bayesian Inf. Crit.      -24.780        -6.815    
## ==================================================
## Note:                  *p<0.1; **p<0.05; ***p<0.01

5.1 Média zero dos resíduos

A maior dificuldade é que não existe teste formal para essa pressuposição. É similar a um erro de especificação do modelo, como por exemplo, com variáveis relevantes omitidas do modelo. O modelo com uma correta especificação provavelmente não terá problemas com média dos resíduos não nula.

Normalmente se faz o teste simples de “H0: média igual a zero” para investigar a violação ou não da pressuposição. Valores elevados para a probabilidade indicarão a não-rejeição da hipótese nula e confirmação da pressuposição.

No R, assim como será visto também em outros softwares, faz-se a estimação da regressão, gera-se a série de resíduos e testa se a média desta série é zero. Para obter os resíduos, faz-se resid(regressao).

• Exemplo sinop (ver capítulo da violação da especificação):

#obtendo residuos do modelo
u.hat<-resid(regressao1)
u.hat

##            1            2            3            4            5            6 
## -0.174850609 -0.012678148 -0.020593163 -0.286070937  0.162888449 -0.096015009 
##            7            8            9           10           11           12 
## -0.120594387 -0.302767125 -0.223128321  0.045685710  0.246465497 -0.019311122 
##           13           14           15           16           17           18 
## -0.187851144 -0.115828370 -0.155105541 -0.230579713  0.096317031 -0.133384939 
##           19           20           21           22           23           24 
## -0.237323386  0.008075268  0.044066328  0.200870553  0.041872603 -0.061033219 
##           25           26           27           28           29           30 
## -0.184940788 -0.233421897 -0.289243771 -0.002670153 -0.156434505  0.072496831 
##           31           32           33           34           35           36 
## -0.113523308 -0.358511284  0.218644733  0.121759126 -0.153545584  0.172861769 
##           37           38           39           40           41           42 
## -0.033423018  0.004438535  0.006998865 -0.189545313 -0.064328543  0.148195994 
##           43           44           45           46           47           48 
## -0.141714385  0.017418670 -0.149733068  0.210476148  0.106280133 -0.022229403 
##           49           50           51           52           53           54 
##  0.211855004 -0.032769196  0.685646266 -0.184801630  0.018763768  1.805327879 
##           55           56           57           58           59 
##  0.206069915  0.195336677 -0.056312310 -0.091973193 -0.212575269

Uma primeira forma é ver os gráficos de diagnósticos já preparados pela função lm do R e fazendo graficamente os resíduos em função dos valores estimados da variável dependente (gráfico residual) ou em função dos valores de uma das variáveis independentes. Ou seja, faz-se o gráfico Residuals x fitted com: plot(regressao, which = 1).

Como exemplo, reproduz-se na Figura 5.1 o gráfico Residuals x fitted do exemplo da produção de soja em setores censitários de Sinop-MT:

reg1.plot1<-plot(regressao1, which = 1)

Figure 5.1: Figura 5.1 . Gráfico dos resíduos versus valores previstos da variável dependente.

Figura 5.1 . Gráfico dos resíduos versus valores previstos da variável dependente.

Na Figura 5.1, os resultados devem se distribuir aleatoriamente ao longo de uma ‘reta’ com ‘dispersão’ uniforme, ou seja, dentro de alguma ‘banda’ de variação, em torno da abscissa ‘zero’. Havendo algum padrão tem-se indicação de presença de autocorrelação residual. Havendo dispersão variável tem-se indicação de presença de heterocedasticidade residual. Fora da abscissa zero, a média é não nula. A hipótese de teste H0: a média dos resíduos é zero. O script para R está na Figura 5.2.

# teste dos residuos do modelo
u.hat_reg1<-resid(regressao1)
t.test(u.hat_reg1)

## 
##  One Sample t-test
## 
## data:  u.hat_reg1
## t = -2.2888e-17, df = 58, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.07798565  0.07798565
## sample estimates:
##     mean of x 
## -8.917103e-19

Figura 5.2. Script para o teste de média zero dos resíduos no R, exemplo dos setores censitários de Sinop-MS.

Olhando os resultados para a hipótese H0: a média dos resíduos é zero, tem-se logo abaixo a hipótese alternativa: alternative hypothesis: true mean is not equal to 0. Para o valor da estatística de teste \(t = -2.2888e-17\), a probabilidade p-value = 1 indica pela não rejeição de H0, e a média pode ser considerada como zero.

5.2 Normalidade dos resíduos

Nesta seção, analisa-se a Pressuposição: O resíduo estimado tem distribuição normal, com média zero e variância constante. Um parágrafo bastante útil é o de Fox (2016, p.297):

The assumption of normally distributed errors is almost always arbitrary. Nevertheless, the central limit theorem ensures that, under very broad conditions, inference based on the leastsquares estimator is approximately valid in all but small samples. Why, then, should we be concerned about non-normal errors?

Ou seja, em geral, a violação deste pressuposto pode ser justificada em muitos casos pela alegação do Teorema do Limite Central, ou no inglês, \(CLT\) (Central Limit Theorem). A ideia de que no limite, quando se aumenta o tamanho da amostra, a distribuição tende à normalidade.

A teoria estatística, de outro lado, assegura que o estimador de \(MQO\) é o melhor estimador linear não tendencioso somente em presença de resíduos normais. As distribuições de resíduos com “caudas pesadas” (heavy-tailed), portanto, não-normais, poderão significar presença de outliers (observações atípicas). Muitas vezes o argumento é de que os estimadores de mínimos quadrados necessitam as propriedades de Gauss-Markov para eficiência e não-tendenciosidade, mas que não requerem a normalidade dos resíduos (FOX, 2016).

É um fato comumente confirmado que os testes formais de normalidade rejeitam a hipótese de normalidade dos resíduos em conjuntos amostrais grandes, como os big data. Pode-se provar que quando o tamanho da amostra (\(n\)) fica muito grande, o menor desvio da perfeita normalidade levará a um resultado significativo. Considerando que todo conjunto de dados tem algum grau de aleatoriedade, será improvável que um conjunto de dados será uma amostra perfeita normalmente distribuída. Na estatística aplicada, a questão não é se os resíduos são perfeitamente normais, mas distribuídos normalmente o suficiente para as suposições serem válidas. De outro lado, alguns testes como os t-Student e a ANOVA (Análise de Variância) serão robustos à normalidade em grandes amostras.

Em geral, os testes questionam: há evidência convincente de qualquer desvio do ideal gaussiano? Com conjuntos de dados reais moderadamente grandes, a resposta é quase sempre sim. Usualmente o procedimento de duas etapas recomenda a inspeção gráfica e posteriormente o teste estatístico.

Um método gráfico usual é por meio da comparação quantílica, ou o chamado Q-Q plot, em que se compara a distribuição amostral dos resíduos padronizados com a distribuição normal dos quantis, \(N(0,1)\), ou alternativamente, com a distribuição \(t\) dos quantis para \(n-k-2\) graus de liberdade. As duas possibilidades gerarão resultados convergentes para amostras maiores. Os resíduos padronizados são preferidos por terem variâncias semelhantes. A comparação gráfica do Q-Q plot será efetiva na identificação de outliers, assimetria (skewness), e caudas leves/pesadas presentes na curtose (kurtosis). Em amostras grandes, o histograma poderá dar uma boa ideia da distribuição.

A assimetria é o grau de distorção relativamente à curva de distribuição normal, simétrica. De modo geral, avalia-se: \(-0,5 \le skewness \le 0,5\) diz-se que os dados são bastante simétricos; \(-1,0 \le skewness \le -0,5\) diz-se que tem inclinação negativa; ou entre \(0,5 \le skewness \le 1,0\) diz-se inclinação positiva, e os dados são inclinados moderadamente. Se a assimetria for menor que -1 (inclinada negativamente) ou maior que 1 (inclinada positivamente), os dados serão considerados assimétricos.

Já com respeito à curtose, pode-se dizer ser uma medida da presença de outliers na distribuição, que podem ser tanto negativos como positivos. Curtoses elevadas indicam distribuições de dados com caudas pesadas e presença de outliers. A curtose próxima a 3 indica uma distribuição próxima da normal. A curtose maior que 3 indica uma distribuição leptocúrtica enquanto valores de curtose abaixo de 3 indicam a distribuição platicúrtica.

A curtose está associada ao achatamento da distribuição, quanto mais chata menor o valor da curtose (\(K\)). Exemplo: \(K>3\) (distribuição mais “em pé”), \(K<3\) (distribuição mais “achatada”). A curtose (em inglês kurtosis) é definida pela expressão de \(K\):

para

\[ {\hat \sigma ^2} = Var\left( \varepsilon \right) = \hat s\sqrt {\frac{{N - 1}}{N}} \]

\[ K = \frac{1}{N}\sum\nolimits_{t = 1}^N {{{\left[ {\frac{{{Y_i} - \bar Y}}{{Var\left( \varepsilon \right)}}} \right]}^4}} \]

A simetria (do inglês skewness) é dada pela expressão \(S\):

\[ S = \frac{1}{N}\sum\nolimits_{t = 1}^N {{{\left[ {\frac{{{Y_i} - \bar Y}}{{Var\left( \varepsilon \right)}}} \right]}^3}} \]

A normalidade multivariada pode ser difícil de avaliar graficamente e a convergência para distribuições assintóticas pode ser demorada para estatísticas multivariadas. Portanto, os testes para normalidade serão mais úteis em uma configuração multivariada.

O teste para detecção mais usual é o Bera-Jarque, ou teste \(BJ\) (muitas vezes aparece como Jarque-Bera ou \(JB\), referentes aos autores Carlos M. Jarque e Anil K. Bera), o qual testa a simetria e a curtose da distribuição dos resíduos em relação à curva normal.

A estatística do teste de normalidade de Jarque-Bera (\(JB\)) será:

\[ JB = \frac{1}{6} \left( {n - p} \right) \left[ {S^2 + {\textstyle {\frac {1}{4}}}{{\left( {K - 3} \right)}^2}} \right] ∼ \chi _{gl = 2}^2 \]

Se rejeitar H0: erros normais, ou seja, p-value menor que 10%, 5% ou 1% de significância, tenho que descobrir qual é a distribuição real dos resíduos e fazer nova dedução do estimador dos parâmetros. Em geral, para amostras grandes, aplica-se o Teorema do Limite Central argumentando-se que no limite tem-se a normalidade da distribuição dos resíduos.

Outros testes de normalidade podem ser listados, do pacote nortest do R: Kolmogorov-Smirnov (KS), Lilliefors, Cramér-von Mises, Shapiro-Wilk, Shapiro-Francia, Anderson-Darling.

• Exemplo sinop:

Seja o projeto dummies.Rproj e o script dummies e outliers.R. Neste caso tem-se 59 observações dos setores censitários de Sinop-MT, o mesmo dataset utilizado para investigar os outliers. O script apresenta o teste de normalidade por meio do gráfico QQPlot, seguido do teste de Bera-Jarque e por testes alternativos variados. Vários pacotes podem ser utilizados em diferentes casos: car, tseries, moments, nortest, stats.

Existem dois tipos de gráficos de probabilidade normal: 1º tipo: representa a probabilidade acumulada que seria de esperar se a distribuição fosse normal, em função da probabilidade observada acumulada dos erros (Normal P-P Plot); 2º tipo: representa o quantil de probabilidade esperado se a distribuição fosse normal em função dos resíduos (Normal Q-Q Plot) (Figura 5.3.). Observe que pela função plot básica (à esquerda), o gráfico resultante apresenta a mesma informação que pelo pacote car e função qqPlot, mas esta última tem a vantagem de apresentar bandas de confiança. Neste caso específico, é possível suspeitar de que os pontos estão em torno da linha pontilhada, indicando pela normalidade, mas com outliers 51 e 54. A ilustração da direita, ao olhar as bandas tracejadas, verifica-se que a reta da normal sai das bandas, o que revela preocupação com a normalidade dos resíduos.

split.screen(c(1,2))

## [1] 1 2

screen(1)
reg1.plot2<-plot(regressao1, which = 2) 
screen(2)
car::qqPlot(regressao1)

## [1] 51 54

close.screen(all = TRUE)

Figura 5.3. Gráficos Q-Q Normal pela função plot da regressão (esquerda) e pelo pacote car, função qqPlot (direita), no R.

Pelo pacote tseries do R, o teste Jarque-Bera será como na Figura 5.4.

library(tseries)
JB.reg1<-jarque.bera.test(u.hat)
JB.reg1

## 
##  Jarque Bera Test
## 
## data:  u.hat
## X-squared = 1224.2, df = 2, p-value < 2.2e-16

Figura 5.4. Teste de normalidade Jarque-Bera pelo tseries em R.

A alternativa de teste pelo pacote moments usando o jarque.test(u.hat) para H0: resíduos normais, será como na Figura 5.5.

# criando dataset excluindo as observacoes 51 e 54
dados2<-dados[dados$D51_54==0, ]
attach(dados2)
#View(dados2)
regressao2<-lm(log(VP)~log(AREA)+log(QMO)+
                 log(PPRAGA)+log(PROT),data=dados2)
u.hat_reg2<-resid(regressao2)
library(moments)
(JB_reg1<-jarque.test(u.hat))

## 
##  Jarque-Bera Normality Test
## 
## data:  u.hat
## JB = 1224.2, p-value < 2.2e-16
## alternative hypothesis: greater

(JB_reg2<-jarque.test(u.hat_reg2))  # sem observacoes 51 e 54

## 
##  Jarque-Bera Normality Test
## 
## data:  u.hat_reg2
## JB = 2.1565, p-value = 0.3402
## alternative hypothesis: greater

Figura 5.5. Teste de normalidade Jarque-Bera pelo moments em R.

Outros testes são aqui apresentados: Kolmogorov-Smirnov (KS), Lilliefors, Cramér-von Mises, Shapiro-Wilk, Shapiro-Francia, Anderson-Darling, do pacote nortest do R (Figura 5.6.).

Figura 5.6. Script para diversos testes de normalidade em R.

# Pacote com alguns testes
library(nortest)
u.hat <- resid(regressao1)
# Testes para a regressao 1
t1.1 <- ks.test(u.hat,"pnorm",mean(u.hat),summary(regressao1)$sigma)# KS
t2.1 <- lillie.test(u.hat)       # Lilliefors
t3.1 <- cvm.test(u.hat)          # Cramér-von Mises
t4.1 <- shapiro.test(u.hat)      # Shapiro-Wilk
t5.1 <- sf.test(u.hat)           # Shapiro-Francia
t6.1 <- ad.test(u.hat)           # Anderson-Darling

O resultado do script será, para a regressão 1 do exemplo de Sinop-MT, contendo todas as observações, inclusive os outliers 51 e 54, a Tabela 4.2.

Tabela 4.2. Resultados de diversos testes de normalidade dos resíduos, em R, exemplo de Sinop-MT.

# Tabela de resultados
testes <- c(t1.1$method, t2.1$method, t3.1$method, 
             t4.1$method, t5.1$method,t6.1$method)
estt <- as.numeric(c(t1.1$statistic, t2.1$statistic, t3.1$statistic, 
                      t4.1$statistic, t5.1$statistic, t6.1$statistic))
valorp <- c(t1.1$p.value, t2.1$p.value, t3.1$p.value, 
             t4.1$p.value, t5.1$p.value,t6.1$p.value)
resultados <- cbind(estt, valorp)
rownames(resultados) <- testes
colnames(resultados) <- c("Estatística", "prob.")
knitr::kable(resultados, digits = 3)

Portanto, probabilidades baixas (menores que 0.10 = 10%) indicam a rejeição de H0 e conclui-se que a distribuição não é normal.

6.1 Teste de White de Heterocedasticidade dos resíduos

Por exemplo, para uma regressão da forma (6.3):

\[\begin{equation} Y_i =\beta_1 + \beta_2 X_{2i} + \beta_3 X_{3i}+ \varepsilon_i, \end{equation}\] \tag{6.3}

o teste de White é implementado manualmente da seguinte forma:

1. Estima-se a regressão inicial e obtêm-se os resíduos \(e_i\);

2. Faz-se uma regressão auxiliar do tipo (6.4):

\[\begin{equation} e_i^2=\alpha_1 + \alpha_2 X_{2i} + \alpha_3 X_{3i}+ \alpha_4 X_{2i}^2+ \alpha_5 X_{3i}^2+ \alpha_6 X_{2i} X_{3i}+ \nu_i, \end{equation}\] \tag{6.4}

ou seja, o quadrado dos resíduos estimados como função das variáveis explicativas (neste caso, \(X_2\) e \(X_3\)), dos quadrados das variáveis explicativas (\(X_2^2\) e \(X_3^2\)) e do produto cruzado das variáveis explicativas (\(X_2 X_3\)). Deve-se incluir o termo do intercepto (\(\alpha_1\)) mesmo que na regressão original não o tenha.

3. Analisa-se o \(R^2\) da regressão auxiliar multiplicado pelo tamanho da amostra (\(n\)) comparando com o valor da tabela qui-quadrado (\(\chi^2\)) para graus de liberdade iguais ao número total de regressores da equação auxiliar. No nosso exemplo,

\(nR^2 \sim \chi^2\) com gl = 5 (para \(X_{2i}, X_{3i},X_{2i}^2,X_{3i}^2,X_{2i}X_{3i}\))

• Se \(nR^2 > \chi^2\) tabelado, então existe heterocedasticidade;

• Se \(nR^2 < \chi^2\) tabelado, então \(\alpha_2 = \alpha_3 = \alpha_4 = \alpha_5 = \alpha_6 = 0\) , e não existe heterocedasticidade.

Deve-se tomar cuidado com este teste, pois ele prevê a inclusão de termos adicionais que, em presença de muitas variáveis \(X\), poderá comprometer o modelo com relação aos seus graus de liberdade.

library(readxl)
# library(foreign)
dados <- read_excel("soja_apostila.xlsx", 
                    sheet = "dados")
attach(dados)
# summary(dados) # omiti porque já consta do capítulo 4

# Teste de White para heterocedasticidade
regressao1<-lm(QSOJA~FERTILIZANTE+TRATOR+MO)
m <- regressao1
data <- dados
#rotina do teste com base em m e data
u2 <- m$residuals^2
#sem termos cruzados, no cross-terms
reg.auxiliar <- lm(u2 ~ FERTILIZANTE+TRATOR+MO+
                     I(FERTILIZANTE^2)+I(TRATOR^2)+I(MO^2))  
summary(reg.auxiliar)

## 
## Call:
## lm(formula = u2 ~ FERTILIZANTE + TRATOR + MO + I(FERTILIZANTE^2) + 
##     I(TRATOR^2) + I(MO^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3451.5 -1319.0  -345.4   365.3 18550.3 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          9516.00    6550.45   1.453   0.1491  
## FERTILIZANTE         -810.45     551.79  -1.469   0.1448  
## TRATOR              -2061.75    2106.30  -0.979   0.3298  
## MO                  76907.49   29850.12   2.576   0.0113 *
## I(FERTILIZANTE^2)      18.85      15.37   1.226   0.2227  
## I(TRATOR^2)           245.72     236.74   1.038   0.3016  
## I(MO^2)           -230942.82  103719.00  -2.227   0.0280 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2805 on 110 degrees of freedom
## Multiple R-squared:  0.1389, Adjusted R-squared:  0.09188 
## F-statistic: 2.956 on 6 and 110 DF,  p-value: 0.0103

Ru2<- summary(reg.auxiliar)$r.squared
LM <- nrow(data)*Ru2
#obtendo o numero de regressores menos o intercepto
k <- length(coefficients(reg.auxiliar))-1
# O Teste tem k termos regressores em 'reg.auxiliar'
p.value <- 1-pchisq(LM, k) 
#c("LM","p.value")
c(LM, p.value)

## [1] 16.24547461  0.01249533

library(lmtest)
bptest(regressao1,
       ~ FERTILIZANTE+TRATOR+MO+
         I(FERTILIZANTE^2)+I(TRATOR^2)+I(MO^2))

## 
##  studentized Breusch-Pagan test
## 
## data:  regressao1
## BP = 16.245, df = 6, p-value = 0.0125

reg.auxiliar2 <- lm(u2 ~ FERTILIZANTE+I(FERTILIZANTE*FERTILIZANTE)+
                     I(FERTILIZANTE*TRATOR)+I(FERTILIZANTE*MO)+
                     TRATOR+I(TRATOR*TRATOR)+I(TRATOR*MO)+
                     MO+I(MO*MO)+I(FERTILIZANTE*TRATOR*MO))
summary(reg.auxiliar2)

## 
## Call:
## lm(formula = u2 ~ FERTILIZANTE + I(FERTILIZANTE * FERTILIZANTE) + 
##     I(FERTILIZANTE * TRATOR) + I(FERTILIZANTE * MO) + TRATOR + 
##     I(TRATOR * TRATOR) + I(TRATOR * MO) + MO + I(MO * MO) + I(FERTILIZANTE * 
##     TRATOR * MO))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5430.3 -1048.5  -411.8   702.2 14557.9 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     -54297.38   17767.60  -3.056 0.002839 ** 
## FERTILIZANTE                      2199.87     972.89   2.261 0.025790 *  
## I(FERTILIZANTE * FERTILIZANTE)      31.90      15.78   2.021 0.045763 *  
## I(FERTILIZANTE * TRATOR)          -621.37     237.26  -2.619 0.010113 *  
## I(FERTILIZANTE * MO)            -40494.73   12480.66  -3.245 0.001574 ** 
## TRATOR                            9647.85    4935.84   1.955 0.053258 .  
## I(TRATOR * TRATOR)                 374.98     250.52   1.497 0.137419    
## I(TRATOR * MO)                 -152774.51   58702.12  -2.603 0.010578 *  
## MO                              933331.63  240560.44   3.880 0.000182 ***
## I(MO * MO)                     -354820.55  127882.92  -2.775 0.006534 ** 
## I(FERTILIZANTE * TRATOR * MO)     7015.27    3094.79   2.267 0.025432 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2551 on 106 degrees of freedom
## Multiple R-squared:  0.3139, Adjusted R-squared:  0.2492 
## F-statistic:  4.85 on 10 and 106 DF,  p-value: 9.664e-06

Ru2<- summary(reg.auxiliar2)$r.squared
LM <- nrow(data)*Ru2
#obtendo o numero de regressores menos o intercepto
k <- length(coefficients(reg.auxiliar2))-1
k

## [1] 10

p.value <- 1-pchisq(LM, k) 
c(LM, p.value)

## [1] 3.672653e+01 6.310037e-05

library(lmtest)
# reg.auxiliar <- lm(u2 ~ FERTILIZANTE+I(FERTILIZANTE*FERTILIZANTE)+
#                     I(FERTILIZANTE*TRATOR)+I(FERTILIZANTE*MO)+
#                     TRATOR+I(TRATOR*TRATOR)+I(TRATOR*MO)+
#                     MO+I(MO*MO)+I(FERTILIZANTE*TRATOR*MO))
bptest(regressao1,
       ~ FERTILIZANTE+I(FERTILIZANTE*FERTILIZANTE)+
         I(FERTILIZANTE*TRATOR)+I(FERTILIZANTE*MO)+
         TRATOR+I(TRATOR*TRATOR)+I(TRATOR*MO)+
         MO+I(MO*MO)+I(FERTILIZANTE*TRATOR*MO))

## 
##  studentized Breusch-Pagan test
## 
## data:  regressao1
## BP = 36.727, df = 10, p-value = 6.31e-05

#regressao1<-lm(QSOJA~FERTILIZANTE+TRATOR+MO)
#library(car) 
#possibilidades: hccm(regressao1,type=c("hc0","hc1","hc2","hc3","hc4"))
vcov.white<-car::hccm(regressao1,type=c("hc1"))
coeftest(regressao1,vcov.white)

## 
## t test of coefficients:
## 
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   494.96573   26.46010 18.7061 < 2.2e-16 ***
## FERTILIZANTE   -0.55354    1.27243 -0.4350   0.66437    
## TRATOR        -33.68994    3.79488 -8.8777 1.198e-14 ***
## MO           -209.14074  115.17696 -1.8158   0.07205 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(estimatr)
library(car)
# Erros padrões robustos HC1 – padrão do Stata
reg.robust <- lm_robust(QSOJA~FERTILIZANTE+TRATOR+MO, 
                   data = dados,
                   se_type = "stata")
summary(reg.robust)

## 
## Call:
## lm_robust(formula = QSOJA ~ FERTILIZANTE + TRATOR + MO, data = dados, 
##     se_type = "stata")
## 
## Standard error type:  HC1 
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper  DF
## (Intercept)   494.9657     26.460  18.706 2.147e-36  442.543  547.388 113
## FERTILIZANTE   -0.5535      1.272  -0.435 6.644e-01   -3.074    1.967 113
## TRATOR        -33.6899      3.795  -8.878 1.198e-14  -41.208  -26.172 113
## MO           -209.1407    115.177  -1.816 7.205e-02 -437.327   19.046 113
## 
## Multiple R-squared:  0.4651 ,    Adjusted R-squared:  0.4509 
## F-statistic: 33.94 on 3 and 113 DF,  p-value: 1.021e-15

# alternativa HC3 - padrão do R
options("scipen"=0, "digits"=7) 
model <- lm_robust(QSOJA~FERTILIZANTE+TRATOR+MO, 
                   data = dados, 
                   se_type = "HC3")
summary(model)

## 
## Call:
## lm_robust(formula = QSOJA ~ FERTILIZANTE + TRATOR + MO, data = dados, 
##     se_type = "HC3")
## 
## Standard error type:  HC3 
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper  DF
## (Intercept)   494.9657     27.746 17.8390 1.194e-34  439.995  549.936 113
## FERTILIZANTE   -0.5535      1.308 -0.4233 6.729e-01   -3.144    2.037 113
## TRATOR        -33.6899      3.896 -8.6463 4.069e-14  -41.410  -25.970 113
## MO           -209.1407    129.250 -1.6181 1.084e-01 -465.208   46.927 113
## 
## Multiple R-squared:  0.4651 ,    Adjusted R-squared:  0.4509 
## F-statistic: 32.44 on 3 and 113 DF,  p-value: 3.334e-15

regressao1<-lm(QSOJA~FERTILIZANTE+TRATOR+MO)
stargazer::stargazer(regressao1,
          se = starprep(regressao1,se_type = "HC3"),
          p = starprep(regressao1,se_type = "HC3", stat = "p.value"),
          type="text",style="all",
          omit.stat = "f")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                QSOJA           
## -----------------------------------------------
## FERTILIZANTE                  -0.554           
##                               (1.308)          
##                             t = -0.423         
##                              p = 0.673         
## TRATOR                      -33.690***         
##                               (3.896)          
##                             t = -8.646         
##                              p = 0.000         
## MO                           -209.141          
##                              (129.250)         
##                             t = -1.618         
##                              p = 0.109         
## Constant                    494.966***         
##                              (27.746)          
##                             t = 17.839         
##                              p = 0.000         
## -----------------------------------------------
## Observations                    117            
## R2                             0.465           
## Adjusted R2                    0.451           
## Residual Std. Error      41.506 (df = 113)     
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

dados <- structure(list(obs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 
    15, 16, 17, 18), VENDAS = c(6375.3, 11626.4, 14655.1, 21869.2, 26408.3, 
    32405.6, 35107.7, 40295.4, 70761.6, 80552.8, 95294, 101314.1, 116141.3, 
    122315.7, 141649.9, 175025.8, 230614.5, 293543), PD = c(62.5, 92.9, 178.3, 
    258.4, 494.7, 1083, 1620.6, 421.7, 509.2, 6620.1, 3918.6, 1595.3, 6107.5, 
    4454.1, 3163.8, 13210.7, 1703.8, 9528.2), LUCROS = c(185.1, 1569.5, 276.8, 
    2828.1, 225.9, 3751.9, 2884.1, 4645.7, 5036.4, 13869.9, 4487.8, 10278.9, 
    8787.29999999999, 16438.8, 9761.4, 19774.5, 22626.6, 18415.4)), 
    row.names = c(NA,-18L), class = c("tbl_df", "tbl", "data.frame"))

library(stargazer);library(car);library(lmtest);library(AER)
# VENDAS = vendas setoriais PD = pesquisa e desenvolvimento LUCRO = lucro
# setorial
knitr::kable(dados)

attach(dados)
options(scipen = 999, digits = 4)
regressao1 <- lm(PD ~ VENDAS)
summary(regressao1)

## 
## Call:
## lm(formula = PD ~ VENDAS)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -5846   -951   -402    346   7434 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 192.99311  990.98579    0.19   0.8480   
## VENDAS        0.03190    0.00833    3.83   0.0015 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2760 on 16 degrees of freedom
## Multiple R-squared:  0.478,  Adjusted R-squared:  0.446 
## F-statistic: 14.7 on 1 and 16 DF,  p-value: 0.00148

# gerar residuos
uhat <- regressao1$residuals
uhat 

##        1        2        3        4        5        6        7        8 
##  -333.87  -470.98  -482.20  -632.23  -540.73  -143.74   307.66 -1056.73 
##        9       10       11       12       13       14       15       16 
## -1941.11  3857.45   685.70 -1829.65  2209.56   359.20 -1547.87  7434.33 
##       17       18 
## -5845.87   -28.91

# residuos quadrados
uhat2 <- uhat^2
uhat2 

##          1          2          3          4          5          6          7 
##   111467.4   221821.3   232512.7   399712.1   292385.3    20661.9    94654.4 
##          8          9         10         11         12         13         14 
##  1116677.8  3767914.5 14879888.0   470179.8  3347606.6  4882159.0   129021.3 
##         15         16         17         18 
##  2395907.8 55269198.4 34174223.2      835.9

# residuos em valores absolutos
uhat.abs <- abs(uhat)
uhat.abs

##       1       2       3       4       5       6       7       8       9      10 
##  333.87  470.98  482.20  632.23  540.73  143.74  307.66 1056.73 1941.11 3857.45 
##      11      12      13      14      15      16      17      18 
##  685.70 1829.65 2209.56  359.20 1547.87 7434.33 5845.87   28.91

# Estimar “uhat.abs” em função de Vendas: ou seja, módulo dos resíduos em
# função de vendas:
reg.glejser <- lm(uhat.abs ~ VENDAS)
summary(reg.glejser)

## 
## Call:
## lm(formula = uhat.abs ~ VENDAS)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -4054   -714   -259    194   4766 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 578.5710   678.6950    0.85    0.407  
## VENDAS        0.0119     0.0057    2.09    0.053 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1890 on 16 degrees of freedom
## Multiple R-squared:  0.215,  Adjusted R-squared:  0.166 
## F-statistic: 4.38 on 1 and 16 DF,  p-value: 0.0526

options(scipen = 0,digits = 7)
library(lmtest)
bptest(regressao1, ~VENDAS + I(VENDAS^2))

## 
##  studentized Breusch-Pagan test
## 
## data:  regressao1
## BP = 5.2125, df = 2, p-value = 0.07381

m <- regressao1
data <- dados
# rotina do teste com base em m e data
u2 <- m$residuals^2
reg.auxiliar <- lm(u2 ~ VENDAS + I(VENDAS^2))  #Com termos cruzados, cross-terms
summary(reg.auxiliar)

## 
## Call:
## lm(formula = u2 ~ VENDAS + I(VENDAS^2))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -14823352  -8256873   -830623   3687265  37800060 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -6.220e+06  6.460e+06  -0.963   0.3509  
## VENDAS       2.294e+02  1.262e+02   1.817   0.0892 .
## I(VENDAS^2) -5.371e-04  4.495e-04  -1.195   0.2507  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13200000 on 15 degrees of freedom
## Multiple R-squared:  0.2896, Adjusted R-squared:  0.1949 
## F-statistic: 3.057 on 2 and 15 DF,  p-value: 0.07698

Ru2 <- summary(reg.auxiliar)$r.squared
LM <- nrow(data) * Ru2
# obtendo o numero de regressores menos o intercepto
k <- length(coefficients(reg.auxiliar)) - 1
p.value <- 1 - pchisq(LM, k)  # O TESTE TEM k TERMOS REGRESSORES EM reg.auxiliar
# c('LM','p.value')
out <- list(LM,p.value)  
names(out) <- c('LM','p.value')
out

## $LM
## [1] 5.212492
## 
## $p.value
## [1] 0.07381112

# regressao1<- lm(PD~VENDAS) library(car) possibilidades:
# hccm(regressao1,type=c('hc0','hc1','hc2','hc3','hc4'))
vcov.white0 <- car::hccm(regressao1, type = c("hc1"))
# 
coeftest(regressao1, vcov.white0)

## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 192.993110 533.931735  0.3615 0.722487   
## VENDAS        0.031900   0.010147  3.1438 0.006276 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

regressao1<- lm(PD~VENDAS)
cov <- car::hccm(regressao1, type = c("hc1"))
robust.se <- sqrt(diag(cov))

stargazer::stargazer(regressao1, regressao1,
            se=list(NULL, robust.se),
            column.labels=c("default","robusto HC1"), 
            align=TRUE, type="text",
            style="all",omit.stat = "f")

## 
## ==========================================================
##                                   Dependent variable:     
##                               ----------------------------
##                                            PD             
##                                  default      robusto HC1 
##                                    (1)            (2)     
## ----------------------------------------------------------
## VENDAS                           0.032***      0.032***   
##                                  (0.008)        (0.010)   
##                                 t = 3.830      t = 3.144  
##                                 p = 0.002      p = 0.002  
## Constant                         192.993        192.993   
##                                 (990.986)      (533.932)  
##                                 t = 0.195      t = 0.361  
##                                 p = 0.849      p = 0.718  
## ----------------------------------------------------------
## Observations                        18            18      
## R2                                0.478          0.478    
## Adjusted R2                       0.446          0.446    
## Residual Std. Error (df = 16)   2,759.153      2,759.153  
## ==========================================================
## Note:                          *p<0.1; **p<0.05; ***p<0.01

6.2 Teste de Breusch-Pagan-Godfrey

Neste caso, seguem-se os passos:

1. Estima por \(MQO\) e obtém resíduos;

2. Obtém a variância dos resíduos (\(Var(û)\));

3. Constrói variáveis \(p = û²/Var(û)\);

4. Faz regressão auxiliar de p contra variáveis \(Z\), em que podem ser alguns ou todos os X como parte de \(Z\);

5. Obtém a soma dos quadrados da regressão, faz a divisão pela metade e compara com distribuição qui-quadrado (\(\chi^2\));

1. Se a probabilidade de \(LM = SQReg/(2 \hat \sigma^4 )\) < 0.10, (10%), rejeita-se H0 e existirá heterocedasticidade;

2. H0: não tem heterocedasticidade dos resíduos.

6. Outra estatística é a LM (BP): \(LM (BP) = nR^2_{\hat u^2}\).
   No R, diferentes pacotes podem ser utilizados com a finalidade de fazer o teste de \(BPG\), ou às vezes chamado apenas de “BP Test”. São eles: library(AER) ou library(lmtest). Por padrão, o teste é feito com resíduos padronizados na forma de Student, ou ‘Studentizados’, assim como no EViews (seção 4.4.2.3) (Observar que no Stata o padrão do hettest (comando do Stata) é sem studentizar. Fazendo para o exemplo de Quantidade de soja em função de fertilizante, trator e Mão-de-obra, como feito para o teste RESET na seção 4.1.1.1.

regressao1<-lm(QSOJA~FERTILIZANTE+TRATOR+MO)
library(AER)
## teste de Breusch-Pagan - bptest simples
bp.het<-bptest(regressao1, studentize = TRUE)
bp.het

## 
##  studentized Breusch-Pagan test
## 
## data:  regressao1
## BP = 10.087, df = 3, p-value = 0.01784

Pelo pacote lmtest o resultado será o mesmo. Farei aqui com (studentized Breusch-Pagan test) e sem resíduos ‘Studentizados’(Breusch-Pagan test):

lmtest::bptest(regressao1, studentize = FALSE)

## 
##  Breusch-Pagan test
## 
## data:  regressao1
## BP = 15.652, df = 3, p-value = 0.001336

lmtest::bptest(regressao1)

## 
##  studentized Breusch-Pagan test
## 
## data:  regressao1
## BP = 10.087, df = 3, p-value = 0.01784

Para H0: erros são homocedásticos, o resultado de p-value inferior a 5% de significância leva, neste exemplo, à rejeição da hipótese nula. Deve-se então corrigir o modelo fazendo o MQP (mínimos quadrados ponderados).

7.1 Teste de Durbin-Watson e h de Durbin

O teste mais comum para detectar a presença de erros autocorrelacionados é o Teste de Durbin-Watson (\(DW\)). A hipótese nula a ser testada é que

\(H_0: \rho = 0\) , ou seja, não há autocorrelação residual de primeira ordem;

Contra a hipótese alternativa:

\(H_1: \rho \ne 0\) , ou seja, \(\rho>0\) e ocorre autocorrelação residual positiva de primeira ordem; ou \(\rho<0\) e ocorre autocorrelação negativa de primeira ordem.

A estatística de teste é o chamado DW, calculado como em (7.3):

\[\begin{equation} DW = \frac{{\sum\limits_{t = 2}^T {{{\left( {{{\hat \varepsilon }_t} - {{\hat \varepsilon }_{t - 1}}} \right)}^2}} }}{{\sum\limits_{t = 1}^T {{{\hat \varepsilon }_t}^2} }} = 2\left( {1 - \hat \rho } \right) \end{equation}\] \tag{7.3}:

em que
* \(\rho=0\) então \(DW=2\), indicando ausência de autocorrelação;
* \(\rho=+1\) então \(DW=0\), indicando autocorrelação positiva e perfeita; e,
* \(\rho=-1\) então \(DW=4\), indicando autocorrelação negativa e perfeita.
Portanto, deseja-se DW próximo de 2, ou seja, ausência de autocorrelação.
A análise requer a comparação dos valores de DW com valores tabelados, que prevêem duas distribuições de probabilidade entrelaçadas: uma distribuição inferior e outra superior. Elas determinam áreas de aceitação e rejeição da hipótese nula, como na figura a seguir:

em que
dL = limite inferior: vem da tabela para n observações e k variáveis explanatórias; e,

dU = limite superior: vem da tabela para n observações e k variáveis explanatórias.

Exemplo:

Para k= 3 (referente a um modelo com X1, X2 e X3), para n = 30 observações, a tabela de DW para 5% de significância nos fornece dL = 1,214 e dU=1,650, e portanto,

4-dL = 4 – 1,214 = 2,786

4 – dU = 4 – 1,650 = 2,350

Para 0<DW<1,214 = rejeição de Ho e autocorrelação positiva

Para 1,214<DW<1,650 = área inconclusiva

Para 1,650<DW<2,350 = aceitação de Ho e não-autocorrelação

Para 2,350<DW<2,786 = área inconclusiva

As deficiências neste método são:

• presença de áreas inconclusivas (esta deficiência tem sido reduzida pelas rotinas em R que permitem obter o valor da probabilidade associada à estatística de teste e assim facilitando a análise da rejeição de H0);    

• só testa autocorrelação de primeira ordem;    

• deve incluir intercepto na regressão; e,    

• o teste não é válido quando o modelo tem variável dependente defasada como variável explicativa.    

Uma alternativa para contornar essas deficiências é via teste h de Durbin, para casos de variável dependente defasada como variável explicativa. O teste é calculado por (7.4):

\[\begin{equation} h = \rho \sqrt {\frac{T}{{1 - T.Var(\beta )}}} ∼ N(0,1) \end{equation}\] \tag{7.4}:

em que \(T\) é o número de observações e \(\beta\) é o parâmetro da variável dependente defasada. A estatística \(h\) será comparada com os valores críticos de \(\pm1,96\) para 5% de significância.

Teste de Durbin-Watson: Exemplo: seja um modelo de PNB em função da oferta de moeda, para os dados da tabela 5.6, conforme Gujarati (2011, p.157), PNB e estoque de moeda segundo quatro critérios, Estados Unidos, 1970-1983. Pode-se proceder com o script como no quadro. Primeiro se define o modelo como série temporal, realiza a regressão e depois realiza o teste de DW.

dados <- 
  structure(list(ANO = c(1970, 1971, 1972, 1973, 1974, 1975, 1976, 
1977, 1978, 1979, 1980, 1981, 1982, 1983), PNB = c(992.7, 1077.6, 
1185.9, 1326.4, 1434.2, 1549.2, 1718, 1918.3, 2163.9, 2417.8, 
2631.7, 2957.8, 3069.3, 3304.8), M1 = c(216.6, 230.8, 252, 265.9, 
277.6, 291.2, 310.4, 335.4, 363.1, 389.1, 414.9, 441.9, 480.5, 
525.4), M2 = c(628.2, 712.8, 805.2, 861, 908.5, 1023.3, 1163.6, 
1286.7, 1389.1, 1498.5, 1632.6, 1796.6, 1965.4, 2196.3), M3 = c(677.5, 
776.2, 886, 985, 1070.5, 1174.2, 1311.9, 1472.9, 1647.1, 1804.8, 
1990, 2238.2, 2462.5, 2710.4), L = c(816.3, 903.1, 1023, 1141.7, 
1249.3, 1367.9, 1516.6, 1704.7, 1910.6, 2117.1, 2326.2, 2599.8, 
2870.8, 3183.1)), row.names = c(NA, -14L), class = c("tbl_df", 
"tbl", "data.frame"))
attach(dados)
library(knitr)
kable(dados,caption="Dados da tabela 5.6")

dados.st<-ts(dados, start=c(1970), end=c(1983))
PNB.ts<-ts(PNB,start=c(1970), end=c(1983),frequency = 1)
regressao1<-lm(PNB~M1+M2+M3+L, data=dados.st)
summary(regressao1)

## 
## Call:
## lm(formula = PNB ~ M1 + M2 + M3 + L, data = dados.st)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -79.60 -24.04 -13.39  25.84 108.84 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 913.48812  720.11266   1.269    0.236
## M1           -6.01416    6.40016  -0.940    0.372
## M2           -0.36474    0.87823  -0.415    0.688
## M3            2.33646    1.98902   1.175    0.270
## L             0.03084    1.68708   0.018    0.986
## 
## Residual standard error: 65.03 on 9 degrees of freedom
## Multiple R-squared:  0.9952, Adjusted R-squared:  0.9931 
## F-statistic: 467.7 on 4 and 9 DF,  p-value: 1.992e-10

O cálculo para DW será para a situação de autocorrelação de primeira ordem, de uma regressão auxiliar como em (7.3). A função dwtest do pacote lmtest permite acessar os resultados do teste de Durbin-Watson com uma distribuição de probabilidade.

lmtest::dwtest(regressao1)

## 
##  Durbin-Watson test
## 
## data:  regressao1
## DW = 0.70426, p-value = 8.024e-05
## alternative hypothesis: true autocorrelation is greater than 0

No exemplo, DW=0.7, com certeza muito afastado de 2, indicando que existe rejeição de \(H_0\) para primeira ordem. Da tabela dos valores críticos de DW (por exemplo o usuário pode consultar as tabelas estatísticas usuais em livros de estatística) tem-se, para:

K=4 (número de variáveis explanatórias excluindo termo constante);
N=14 (tamanho da amostra)
dL = 0.632, logo, 4-dL = 3.368
dU = 2.030, logo, 4-dU = 1.970
Para DW = 0.70, o leitor verificará que se trata de uma área inconclusiva – o leitor deve atentar que neste exemplo, a área inconclusiva se estende por sobre a área de não rejeição de H0!

O resultado foi pela rejeição de \(H_0\), ou seja, rejeita-se que não tenha autocorrelação até 1ª ordem. Neste caso, a deficiência original de que existem áreas inconclusivas é superada pela função dwtest ao associar a uma distribuição de probabilidade e permitir a rejeição ou não da hipótese nula por meio da comparação do p-value = 8.024e-05, portanto menor que o nível de significância de 10% (o leitor pode optar por 5% ou 1%, se desejar).

7.2 Teste de Breusch-Godfrey ou de LM de Correlação Serial

Outro procedimento que contorna as deficiências do teste DW e tem sido incorporado à maioria dos softwares econométricos é o teste de Breusch-Godfrey de autocorrelação superior, ou de Multiplicador de Lagrange para Correlação Serial (Serial Correlation LM Test).

Este teste é bastante útil para detectar autocorrelação de ordens maiores que a primeira, ou seja, para AR(P>1) em que P é a ordem da autocorrelação. A hipótese nula será:

• \(H_0\): não autocorrelação dos resíduos até ordem \(p\): todos os \(\rho_1,\rho_2,\rho_3,...,\rho_p\) são simultaneamente nulos;
  
• \(H_1\): existe autocorrelação dos resíduos até ordem \(p\): pelo menos um dos coeficientes \(\rho_1,\rho_2,\rho_3,...,\rho_p\) é diferente de zero.

Por exemplo, suponha o seguinte processo auto-regressivo (7.5):

\[\begin{equation} {\varepsilon _t} = {\rho _1}{\varepsilon _{t - 1}} + {\rho _2}{\varepsilon _{t - 2}} + {\rho _3}{\varepsilon _{t - 3}} + \ldots + {\rho _p}{\varepsilon _{t - p}} + {\upsilon _t} \end{equation}\] \tag{7.5}:

A hipótese nula será de que todos os coeficientes de autocorrelação são simultaneamente nulos, ou seja, todos os \(\rho_i = 0\), \(i=(1,2,3,...,p)\) e não há autocorrelação de qualquer ordem. A estatística de teste será um multiplicador de Lagrange (\(LM\)) do tipo (7.6):

\[\begin{equation} LM =(n - p)R_0^2 = (n - p)\left( {\frac{{e'{X_0}{{(X_0'X_0)}^{-1}}X_0'e}}{{e'e}}} \right)∼\chi _p^2 \end{equation}\] \tag{7.6}:

O procedimento será:

1. estimar o modelo de regressão (“a”) pelo método usual de MQO e obter resíduos \(\varepsilon_t\);
   
2. estimar o modelo de \(\varepsilon_t\) como função das demais variáveis \(X\) do modelo “a” e também de variáveis \(\varepsilon_t\) defasadas (\(\varepsilon_{t-1},\varepsilon_{t-2},...,etc\)) (formando assim a matriz de regressores \(X_0\)), utilizando para estas defasagens os resíduos obtidos em “a”;
   
3. obter o valor de \(R^2\) desta regressão “b” (\(R^2_0\));
   
4. A estatística de teste será \(LM_{BG} = (n-p).R^2 ∼ \chi _p^2\) graus de liberdade. \(P\) é o número de defasagens incluídas na regressão “b”.

Se a estatística de teste \(LM_{BG}\) for maior que o valor crítico de \(\chi _p^2\) então se rejeita a hipótese nula e existe autocorrelação serial de ordem \(P\), ou seja, pelo menos um \(\rho_i \ne 0\). Neste teste, podem-se ter variáveis \(X\) ou mesmo \(Y\) defasadas, o que representa uma vantagem sobre o teste \(DW\).

require(lmtest)
bgteste1<-bgtest(regressao1, order = 1, type = c("Chisq"), data = dados.st)
bgteste2<-bgtest(regressao1, order = 2, type = c("Chisq"), data = dados.st)
bgteste3<-bgtest(regressao1, order = 3, type = c("Chisq"), data = dados.st)
bgteste4<-bgtest(regressao1, order = 4, type = c("Chisq"), data = dados.st)
bgteste5<-bgtest(regressao1, order = 5, type = c("Chisq"), data = dados.st)
bgteste6<-bgtest(regressao1, order = 6, type = c("Chisq"), data = dados.st)
bgteste1

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  regressao1
## LM test = 9.4239, df = 1, p-value = 0.002142

bgteste2

## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  regressao1
## LM test = 9.95, df = 2, p-value = 0.006908

bgteste3

## 
##  Breusch-Godfrey test for serial correlation of order up to 3
## 
## data:  regressao1
## LM test = 12.242, df = 3, p-value = 0.006599

bgteste4

## 
##  Breusch-Godfrey test for serial correlation of order up to 4
## 
## data:  regressao1
## LM test = 12.352, df = 4, p-value = 0.01492

bgteste5

## 
##  Breusch-Godfrey test for serial correlation of order up to 5
## 
## data:  regressao1
## LM test = 12.556, df = 5, p-value = 0.02791

bgteste6

## 
##  Breusch-Godfrey test for serial correlation of order up to 6
## 
## data:  regressao1
## LM test = 12.707, df = 6, p-value = 0.04794

# padrao do teste de BG, com distribuição qui-quadrado
bgorder = 1:6  # definindo até a máxima ordem do bgtest
d=NULL
for (p in bgorder) {
   bgtest.chi<-bgtest(regressao1,
                      order = p,type=c("Chisq"), data = dados.st)
   print(bgtest.chi) 
   d = rbind(d,
             data.frame(bgtest.chi$statistic,bgtest.chi$p.value))
   }

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  regressao1
## LM test = 9.4239, df = 1, p-value = 0.002142
## 
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  regressao1
## LM test = 9.95, df = 2, p-value = 0.006908
## 
## 
##  Breusch-Godfrey test for serial correlation of order up to 3
## 
## data:  regressao1
## LM test = 12.242, df = 3, p-value = 0.006599
## 
## 
##  Breusch-Godfrey test for serial correlation of order up to 4
## 
## data:  regressao1
## LM test = 12.352, df = 4, p-value = 0.01492
## 
## 
##  Breusch-Godfrey test for serial correlation of order up to 5
## 
## data:  regressao1
## LM test = 12.556, df = 5, p-value = 0.02791
## 
## 
##  Breusch-Godfrey test for serial correlation of order up to 6
## 
## data:  regressao1
## LM test = 12.707, df = 6, p-value = 0.04794

d

##          bgtest.chi.statistic bgtest.chi.p.value
## LM test              9.423914        0.002141738
## LM test1             9.950036        0.006908395
## LM test2            12.241923        0.006598754
## LM test3            12.352221        0.014915519
## LM test4            12.556160        0.027912726
## LM test5            12.706543        0.047940113

u.hat<-residuals(regressao1)
library(dynlm)
reg.aux<-dynlm(u.hat~lag(u.hat,1), data=dados.st)
summary(reg.aux)

## 
## Time series regression with "ts" data:
## Start = 1, End = 14
## 
## Call:
## dynlm(formula = u.hat ~ lag(u.hat, 1), data = dados.st)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -8.337e-16 -3.118e-16 -2.053e-16 -8.430e-17  3.322e-15 
## 
## Coefficients:
##                Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)   0.000e+00  2.741e-16 0.000e+00        1    
## lag(u.hat, 1) 1.000e+00  5.257e-18 1.902e+17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.026e-15 on 12 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 3.619e+34 on 1 and 12 DF,  p-value: < 2.2e-16

-     chnimp: importações de cloreto de bário chinês    
-     bchlimp: importações totais de cloreto de bário dos EUA    
-     befile6: =1 para todos os 6 meses anteriores ao processo    
-     affile6: =1 para todos os 6 meses posteriores ao processo    
-     afdec6: =1 para todos os 6 meses posteriores à decisão    
-     befile12: =1 para todos os 12 meses anteriores ao processo    
-     affile12: =1 para todos os 12 meses posteriores ao processo    
-     afdec12: =1 para todos os 12 meses posteriores à decisão    
-     chempi: índice de produção de químicos     
-     gas: produção de gasolina    
-     rtwex: índice da taxa cambial    
-     spr: =1 para os meses da primavera dos EUA         
-     sum: =1 para os meses do verão dos EUA         
-     fall: =1 para os meses do outono dos EUA     
-     lchnimp: log(chnimp)    
-     lgas: log(gas)    
-     lrtwex: log(rtwex)    
-     lchempi: log(chempi)    
-     t: time trend    
-     feb: =1 se for mês de fevereiro    
-     mar: =1 se for mês de março    
-     apr: ...    
-     may:    
-     jun:    
-     jul:    
-     aug:    
-     sep:    
-     oct:    
-     nov: ...    
-     dec: =1 se for mês de dezembro    
-     percchn: porcentagem de importações oriundas da China

library(dynlm);library(stargazer);library(wooldridge)
data(barium, package='wooldridge')
dados=barium
attach(dados)
tsdata <- ts(dados, start=c(1978,2), frequency=12)
mod1 <- dynlm(log(chnimp)~log(chempi)+log(gas)+log(rtwex)+
                               befile6+affile6+afdec6,
              data=tsdata )
summary(mod1)

## 
## Time series regression with "ts" data:
## Start = 1978(2), End = 1988(12)
## 
## Call:
## dynlm(formula = log(chnimp) ~ log(chempi) + log(gas) + log(rtwex) + 
##     befile6 + affile6 + afdec6, data = tsdata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.03356 -0.39080  0.03048  0.40248  1.51720 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.80277   21.04537  -0.846   0.3992    
## log(chempi)   3.11719    0.47920   6.505 1.72e-09 ***
## log(gas)      0.19634    0.90662   0.217   0.8289    
## log(rtwex)    0.98302    0.40015   2.457   0.0154 *  
## befile6       0.05957    0.26097   0.228   0.8198    
## affile6      -0.03241    0.26430  -0.123   0.9026    
## afdec6       -0.56524    0.28584  -1.978   0.0502 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5974 on 124 degrees of freedom
## Multiple R-squared:  0.3049, Adjusted R-squared:  0.2712 
## F-statistic: 9.064 on 6 and 124 DF,  p-value: 3.255e-08

require(stargazer)
mod1$AIC <- AIC(mod1)
mod1$BIC <- BIC(mod1)
star.1 <- stargazer(mod1, 
            title="Título: Resultado da Regressão",
            align=TRUE,
            type = "text", style = "all",
            keep.stat=c("aic","bic","rsq", "adj.rsq","n"))

## 
## Título: Resultado da Regressão
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                             log(chnimp)        
## -----------------------------------------------
## log(chempi)                  3.117***          
##                               (0.479)          
##                              t = 6.505         
##                              p = 0.000         
## log(gas)                       0.196           
##                               (0.907)          
##                              t = 0.217         
##                              p = 0.829         
## log(rtwex)                    0.983**          
##                               (0.400)          
##                              t = 2.457         
##                              p = 0.016         
## befile6                        0.060           
##                               (0.261)          
##                              t = 0.228         
##                              p = 0.820         
## affile6                       -0.032           
##                               (0.264)          
##                             t = -0.123         
##                              p = 0.903         
## afdec6                        -0.565*          
##                               (0.286)          
##                             t = -1.978         
##                              p = 0.051         
## Constant                      -17.803          
##                              (21.045)          
##                             t = -0.846         
##                              p = 0.400         
## -----------------------------------------------
## Observations                    131            
## R2                             0.305           
## Adjusted R2                    0.271           
## Akaike Inf. Crit.             245.573          
## Bayesian Inf. Crit.           268.575          
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

library(car); library(lmtest);library(sandwich)
#armazenando o valor de n=numero de observacoes
n<-nrow(dados)  
#obtendo o numero de regressores menos o intercepto
k <- length(coefficients(mod1))-1   
c("núm. de observações: n=",n)

## [1] "núm. de observações: n=" "131"

c("núm. de regressores k=",k)

## [1] "núm. de regressores k=" "6"

# Breusch-Godfrey LM Serial Correlation test:
residual <- resid(mod1)
resreg <- dynlm(residual ~
                  L(residual)+L(residual,2)+L(residual,3)+
                   log(chempi)+log(gas)+log(rtwex)+befile6+
                   affile6+afdec6, data=tsdata )
summary(resreg)