Econometria: exercício preço do cobre nos EUA

Abstract

We analyse a linear multiple regression of copper price in the USA, according to Gujarati and Porter (2011, p.460), portuguese edition.

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Econometria: exercício preço do cobre nos EUA. Campo Grande-MS, Brasil: RStudio/Rpubs, 2020. Disponível em http://rpubs.com/amrofi/copper_exercise.

1 Introdução

Os primeiros passos são criar ou abrir um diretório de trabalho. Se optar por criar um novo projeto, haverá a possibilidade de criar em uma pasta vazia. Em seguida, sugere-se que coloque os dados nesta pasta, se possível em um arquivo MS Excel e chame a planilha de ‘dados’.

Seja o enunciado como em Gujarati e Porter (2011, p.460-461).

Enunciado do exercício. Fonte: Gujarati e Porter (2011, p.460-461).

O dataset pode ser obtido direto da editora do Grupo A em https://loja.grupoa.com.br/econometria-basica-ebook-p988374?tsid=34, ou em http://highered.mheducation.com/sites/dl/free/0070660050/37004/data_sets.zip. Neste caso, a tabela é a 12.7, e colocamos ‘embeded’ no code que pode ser baixado clicando no botão ao alto.

library(readxl)
dados <- read_excel("dados.xlsx")
View(dados)

dados <- structure(list(ANO = c(1951, 1952, 1953, 1954, 1955, 1956, 1957, 1958, 1959, 
    1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 
    1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980), C = c(21.89, 22.29, 19.63, 22.85, 
    33.77, 39.18, 30.58, 26.3, 30.7, 32.1, 30, 30.8, 30.8, 32.6, 35.4, 36.6, 38.6, 
    42.2, 47.9, 58.2, 52, 51.2, 59.5, 77.3, 64.2, 69.6, 66.8, 66.5, 98.3, 101.4), 
    PNB = c(330.2, 347.2, 366.1, 366.3, 399.3, 420.7, 442, 447, 483, 506, 523.3, 
        563.8, 594.7, 635.7, 688.1, 753, 796.3, 868.5, 935.5, 982.4, 1063.4, 1171.1, 
        1306.6, 1412.9, 1528.8, 1700.1, 1887.2, 2127.6, 2628.8, 2633.1), I = c(45.1, 
        50.9, 53.3, 53.6, 54.6, 61.1, 61.9, 57.9, 64.8, 66.2, 66.7, 72.2, 76.5, 81.7, 
        89.8, 97.8, 100, 106.3, 111.1, 107.8, 109.6, 119.7, 129.8, 129.3, 117.8, 
        129.8, 137.1, 145.2, 152.5, 147.1), L = c(220.4, 259.5, 256.3, 249.3, 352.3, 
        329.1, 219.6, 234.8, 237.4, 245.8, 229.2, 233.9, 234.2, 347, 468.1, 555, 
        418, 525.2, 620.7, 588.6, 444.4, 427.8, 727.1, 877.6, 556.6, 780.6, 750.7, 
        709.8, 935.7, 940.9), H = c(1491, 1504, 1438, 1551, 1646, 1349, 1224, 1382, 
        1553.7, 1296.1, 1365, 1492.5, 1634.9, 1561, 1509.7, 1195.8, 1321.9, 1545.4, 
        1499.5, 1469, 2084.5, 2378.5, 2057.5, 1352.5, 1171.4, 1547.6, 1989.8, 2023.3, 
        1749.2, 1298.5), A = c(19, 19.41, 20.93, 21.78, 23.68, 26.01, 27.52, 26.89, 
        26.85, 27.23, 25.46, 23.88, 22.62, 23.72, 24.5, 24.5, 24.98, 25.58, 27.18, 
        28.72, 29, 26.67, 25.33, 34.06, 39.79, 44.49, 51.23, 54.42, 61.01, 70.87)), 
    row.names = c(NA, -30L), class = c("tbl_df", "tbl", "data.frame"))

As variáveis são:

• Ano - ano de 1951 a 1980;

• C (média de 12 meses do preço doméstico de cobre nos EUA em cents por pound);

• PNB é o produto nacional bruto anual em bilhões de doláres;

• I (média de 12 meses do índice de produção industrial);

• L (média de 12 meses do preço do cobre na bolsa de Londres em libras esterlinas);

• H (número de prédios em construção por ano em milhares de unidades); e,

• A (média de 12 meses do preço do alumínio em cents por pound).

knitr::kable(dados)

ANO	C	PNB	I	L	H	A
1951	21.89	330.2	45.1	220.4	1491.0	19.00
1952	22.29	347.2	50.9	259.5	1504.0	19.41
1953	19.63	366.1	53.3	256.3	1438.0	20.93
1954	22.85	366.3	53.6	249.3	1551.0	21.78
1955	33.77	399.3	54.6	352.3	1646.0	23.68
1956	39.18	420.7	61.1	329.1	1349.0	26.01
1957	30.58	442.0	61.9	219.6	1224.0	27.52
1958	26.30	447.0	57.9	234.8	1382.0	26.89
1959	30.70	483.0	64.8	237.4	1553.7	26.85
1960	32.10	506.0	66.2	245.8	1296.1	27.23
1961	30.00	523.3	66.7	229.2	1365.0	25.46
1962	30.80	563.8	72.2	233.9	1492.5	23.88
1963	30.80	594.7	76.5	234.2	1634.9	22.62
1964	32.60	635.7	81.7	347.0	1561.0	23.72
1965	35.40	688.1	89.8	468.1	1509.7	24.50
1966	36.60	753.0	97.8	555.0	1195.8	24.50
1967	38.60	796.3	100.0	418.0	1321.9	24.98
1968	42.20	868.5	106.3	525.2	1545.4	25.58
1969	47.90	935.5	111.1	620.7	1499.5	27.18
1970	58.20	982.4	107.8	588.6	1469.0	28.72
1971	52.00	1063.4	109.6	444.4	2084.5	29.00
1972	51.20	1171.1	119.7	427.8	2378.5	26.67
1973	59.50	1306.6	129.8	727.1	2057.5	25.33
1974	77.30	1412.9	129.3	877.6	1352.5	34.06
1975	64.20	1528.8	117.8	556.6	1171.4	39.79
1976	69.60	1700.1	129.8	780.6	1547.6	44.49
1977	66.80	1887.2	137.1	750.7	1989.8	51.23
1978	66.50	2127.6	145.2	709.8	2023.3	54.42
1979	98.30	2628.8	152.5	935.7	1749.2	61.01
1980	101.40	2633.1	147.1	940.9	1298.5	70.87

2 Análise

A primeira questão pede para fazer a estimação e analisar os resultados, para a equação:

\(lnC_t=\beta_1+\beta_2lnI_t+\beta_3lnL_t+\beta_4lnH_t+\beta_5lnA_t+u_t\)

Portanto, estima-se a equação abaixo:

attach(dados)
reg1 <- lm(log(C) ~ log(I) + log(L) + log(H) + log(A))
summary(reg1)

Call:
lm(formula = log(C) ~ log(I) + log(L) + log(H) + log(A))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.23624 -0.06926  0.03014  0.07209  0.22490 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.500441   1.003020  -1.496 0.147192    
log(I)       0.467509   0.165987   2.817 0.009340 ** 
log(L)       0.279443   0.114726   2.436 0.022328 *  
log(H)      -0.005152   0.142947  -0.036 0.971538    
log(A)       0.441449   0.106508   4.145 0.000341 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1217 on 25 degrees of freedom
Multiple R-squared:  0.9361,    Adjusted R-squared:  0.9259 
F-statistic: 91.54 on 4 and 25 DF,  p-value: 1.491e-14

É possível perceber um ajustamento geral de 93,61% das variações de log(C) explicadas pelas variáveis explicativas. Apenas a variável H (número de prédios em construção por ano em milhares de unidades) e o intercepto não tiveram parâmetros significativos. Podemos afirmar a princípio que I (média de 12 meses do índice de produção industrial), L (média de 12 meses do preço do cobre na bolsa de Londres em libras esterlinas) e A (média de 12 meses do preço do alumínio em cents por pound) afetam C (média de 12 meses do preço doméstico de cobre nos EUA em cents por pound). PNB é o produto nacional bruto anual em bilhões de doláres.

A segunda pergunta pede que se obtenha os resíduos padronizados da regressão e faça um gráfico para avaliar a presença de autocorrelação residual.

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

reg1.stdres = rstandard(reg1)
plot(fitted(reg1), reg1.stdres, ylab = "Resíduos padronizados", xlab = "log(C)", 
    main = "Resíduos padronizados x log(C) previsto")
abline(0, 0)  # a origem

Parece haver uma indicação de autocorrelação conforme o primeiro gráfico aponta na linha vermelha.

Calcularemos agora a estatística Durbin-Watson conforme solicitado na letra c.

car::durbinWatsonTest(reg1)

 lag Autocorrelation D-W Statistic p-value
   1       0.5196715     0.9549399       0
 Alternative hypothesis: rho != 0

Nesta opção do pacote car , ele fornece para quantos lags foi o teste, o valor do parâmetro de autocorrelação (0.5197), a estatística de teste DW=0.95 e o valor da probabilidade = 0, indicando a rejeição de H0 de ausência de autocorrelação de primeira ordem. Portanto, indica a presença de autocorrelação residual de primeira ordem. Outra opção seria com o pacote lmteste a função dwtest.

lmtest::dwtest(reg1)

    Durbin-Watson test

data:  reg1
DW = 0.95494, p-value = 6.388e-05
alternative hypothesis: true autocorrelation is greater than 0

Neste caso, a mesma estatística e o p-value são fornecidos.

O exercício também pede o teste das carreiras, também chamado de teste de Geary, menos usual e menos poderoso que o teste de Breusch-Godfrey (LM de correlação serial). Desta forma, passaremos para o teste de BG de correlação serial, função bgtest do pacote lmtest. Sugerimos começar com ordens altas e reduzir uma a uma para interpretar.

# padrao do teste de BG, com distribuição qui-quadrado definindo até a máxima
# ordem do bgtest
library(lmtest)
bgorder = 1:12
d = NULL
for (p in bgorder) {
    bgtest.chi <- bgtest(reg1, order = p, type = c("Chisq"), data = dados)
    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:  reg1
LM test = 8.8406, df = 1, p-value = 0.002946


    Breusch-Godfrey test for serial correlation of order up to 2

data:  reg1
LM test = 13.507, df = 2, p-value = 0.001167


    Breusch-Godfrey test for serial correlation of order up to 3

data:  reg1
LM test = 13.707, df = 3, p-value = 0.003332


    Breusch-Godfrey test for serial correlation of order up to 4

data:  reg1
LM test = 13.754, df = 4, p-value = 0.008123


    Breusch-Godfrey test for serial correlation of order up to 5

data:  reg1
LM test = 14.273, df = 5, p-value = 0.01396


    Breusch-Godfrey test for serial correlation of order up to 6

data:  reg1
LM test = 15.356, df = 6, p-value = 0.01766


    Breusch-Godfrey test for serial correlation of order up to 7

data:  reg1
LM test = 15.393, df = 7, p-value = 0.03128


    Breusch-Godfrey test for serial correlation of order up to 8

data:  reg1
LM test = 18.336, df = 8, p-value = 0.01884


    Breusch-Godfrey test for serial correlation of order up to 9

data:  reg1
LM test = 18.609, df = 9, p-value = 0.02873


    Breusch-Godfrey test for serial correlation of order up to 10

data:  reg1
LM test = 19.959, df = 10, p-value = 0.02964


    Breusch-Godfrey test for serial correlation of order up to 11

data:  reg1
LM test = 19.96, df = 11, p-value = 0.04589


    Breusch-Godfrey test for serial correlation of order up to 12

data:  reg1
LM test = 21.774, df = 12, p-value = 0.04013

knitr::kable(d)

	bgtest.chi.statistic	bgtest.chi.p.value
LM test	8.840634	0.0029460
LM test1	13.506831	0.0011669
LM test2	13.706948	0.0033324
LM test3	13.754136	0.0081225
LM test4	14.273477	0.0139625
LM test5	15.356109	0.0176607
LM test6	15.392859	0.0312806
LM test7	18.336310	0.0188417
LM test8	18.609052	0.0287300
LM test9	19.959052	0.0296424
LM test10	19.959820	0.0458948
LM test11	21.773878	0.0401336

É possível concluir pela rejeição de H0 em todos os lags testados, até 12.

Podemos testar maiores detalhes fazendo regressões auxiliares com os lags para comportar autocorrelações de ordens superiores. Infelizmente os testes não oferecem as saídas das regressões auxiliares diretamente, e portanto, precisaremos rodar os vários modelos. Existe uma opção de usar coeftest para extrair os coeficientes das regressões auxiliares.

Seja o teste para ordem 12:

bgtest.chi <- bgtest(reg1, order = 12, type = c("Chisq"), data = dados)
coeftest(bgtest.chi)

z test of coefficients:

               Estimate Std. Error z value Pr(>|z|)   
(Intercept)    0.283118   0.914993  0.3094 0.757001   
log(I)         0.078832   0.163492  0.4822 0.629680   
log(L)        -0.076886   0.116984 -0.6572 0.511031   
log(H)        -0.053007   0.129545 -0.4092 0.682409   
log(A)         0.063121   0.097293  0.6488 0.516486   
lag(resid)_1   0.543726   0.262826  2.0688 0.038568 * 
lag(resid)_2  -0.934762   0.354281 -2.6385 0.008328 **
lag(resid)_3  -0.043661   0.365592 -0.1194 0.904939   
lag(resid)_4  -0.384382   0.366618 -1.0485 0.294430   
lag(resid)_5  -0.096566   0.399222 -0.2419 0.808869   
lag(resid)_6  -0.708319   0.363597 -1.9481 0.051404 . 
lag(resid)_7   0.154905   0.388500  0.3987 0.690096   
lag(resid)_8  -0.902708   0.425155 -2.1232 0.033733 * 
lag(resid)_9   0.201572   0.428629  0.4703 0.638160   
lag(resid)_10 -0.912078   0.445163 -2.0489 0.040475 * 
lag(resid)_11  0.263836   0.431326  0.6117 0.540746   
lag(resid)_12 -0.630576   0.372424 -1.6932 0.090424 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

É possível ver significância dos termos lag(resid)-1 (que indica ser em t-1, ou seja, \(resid_{t-1}\), também \(resid_{t-2}\), \(resid_{t-6}\) , \(resid_{t-8}\) , \(resid_{t-10}\) e \(resid_{t-12}\). Ou seja, existem evidências para autocorrelações em ordens superiores.

O ideal é rodar cada regressão auxiliar e ver os resultados. Podemos colocar dentro do loop para ver as saídas de coeficientes.

library(lmtest)
bgorder = 1:12
d = NULL
for (p in bgorder) {
    bgtest.chi <- bgtest(reg1, order = p, type = c("Chisq"), data = dados)
    print(bgtest.chi)
    d = rbind(d, data.frame(bgtest.chi$statistic, bgtest.chi$p.value))
    print(coeftest(bgtest.chi))
}

    Breusch-Godfrey test for serial correlation of order up to 1

data:  reg1
LM test = 8.8406, df = 1, p-value = 0.002946


z test of coefficients:

               Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.3683855  0.8675702 -0.4246 0.671116   
log(I)       -0.1048768  0.1460790 -0.7179 0.472791   
log(L)        0.0892880  0.1022994  0.8728 0.382766   
log(H)        0.0453949  0.1233623  0.3680 0.712888   
log(A)       -0.0097854  0.0913455 -0.1071 0.914690   
lag(resid)_1  0.5670655  0.1790759  3.1666 0.001542 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 2

data:  reg1
LM test = 13.507, df = 2, p-value = 0.001167


z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  -0.552588   0.785758 -0.7033   0.48190    
log(I)       -0.036780   0.134421 -0.2736   0.78438    
log(L)        0.028677   0.095271  0.3010   0.76341    
log(H)        0.082020   0.112179  0.7311   0.46469    
log(A)       -0.017166   0.082432 -0.2082   0.83504    
lag(resid)_1  0.795135   0.184599  4.3074 1.652e-05 ***
lag(resid)_2 -0.489015   0.191703 -2.5509   0.01074 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 3

data:  reg1
LM test = 13.707, df = 3, p-value = 0.003332


z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  -0.512105   0.802317 -0.6383    0.5233    
log(I)       -0.029746   0.137274 -0.2167    0.8284    
log(L)        0.025964   0.096960  0.2678    0.7889    
log(H)        0.078075   0.114254  0.6833    0.4944    
log(A)       -0.025278   0.085213 -0.2966    0.7667    
lag(resid)_1  0.854409   0.219535  3.8919 9.946e-05 ***
lag(resid)_2 -0.592317   0.278293 -2.1284    0.0333 *  
lag(resid)_3  0.126950   0.244219  0.5198    0.6032    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 4

data:  reg1
LM test = 13.754, df = 4, p-value = 0.008123


z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  -0.490677   0.824585 -0.5951 0.5518043    
log(I)       -0.028277   0.140427 -0.2014 0.8404162    
log(L)        0.023640   0.099543  0.2375 0.8122785    
log(H)        0.073984   0.117942  0.6273 0.5304674    
log(A)       -0.020451   0.089258 -0.2291 0.8187755    
lag(resid)_1  0.858280   0.224923  3.8159 0.0001357 ***
lag(resid)_2 -0.625742   0.314986 -1.9866 0.0469701 *  
lag(resid)_3  0.180194   0.329817  0.5463 0.5848287    
lag(resid)_4 -0.063512   0.257158 -0.2470 0.8049276    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 5

data:  reg1
LM test = 14.273, df = 5, p-value = 0.01396


z test of coefficients:

                Estimate  Std. Error z value  Pr(>|z|)    
(Intercept)  -0.35148550  0.84879301 -0.4141 0.6788006    
log(I)       -0.01933142  0.14200382 -0.1361 0.8917161    
log(L)        0.01315185  0.10118451  0.1300 0.8965832    
log(H)        0.04913760  0.12277508  0.4002 0.6889911    
log(A)       -0.00042325  0.09330224 -0.0045 0.9963806    
lag(resid)_1  0.82345285  0.23077755  3.5682 0.0003595 ***
lag(resid)_2 -0.57683867  0.32321491 -1.7847 0.0743115 .  
lag(resid)_3  0.04484388  0.37189353  0.1206 0.9040217    
lag(resid)_4  0.10795368  0.33426332  0.3230 0.7467255    
lag(resid)_5 -0.22046667  0.27128024 -0.8127 0.4163960    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 6

data:  reg1
LM test = 15.356, df = 6, p-value = 0.01766


z test of coefficients:

              Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.314436   0.840915 -0.3739 0.708463   
log(I)       -0.033678   0.141109 -0.2387 0.811361   
log(L)        0.015800   0.100201  0.1577 0.874704   
log(H)        0.036635   0.122008  0.3003 0.763974   
log(A)        0.030550   0.095998  0.3182 0.750308   
lag(resid)_1  0.753322   0.236016  3.1918 0.001414 **
lag(resid)_2 -0.549915   0.320799 -1.7142 0.086491 . 
lag(resid)_3  0.073362   0.368973  0.1988 0.842397   
lag(resid)_4 -0.093621   0.372078 -0.2516 0.801338   
lag(resid)_5  0.024416   0.338858  0.0721 0.942560   
lag(resid)_6 -0.317351   0.267763 -1.1852 0.235941   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 7

data:  reg1
LM test = 15.393, df = 7, p-value = 0.03128


z test of coefficients:

               Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.3148894  0.8628760 -0.3649  0.71516   
log(I)       -0.0277728  0.1474289 -0.1884  0.85058   
log(L)        0.0081459  0.1089275  0.0748  0.94039   
log(H)        0.0364991  0.1251959  0.2915  0.77064   
log(A)        0.0368479  0.1028547  0.3583  0.72015   
lag(resid)_1  0.7388075  0.2516003  2.9364  0.00332 **
lag(resid)_2 -0.5595752  0.3322910 -1.6840  0.09218 . 
lag(resid)_3  0.0747707  0.3786652  0.1975  0.84347   
lag(resid)_4 -0.0889694  0.3824195 -0.2326  0.81603   
lag(resid)_5 -0.0160528  0.3963108 -0.0405  0.96769   
lag(resid)_6 -0.2686314  0.3576347 -0.7511  0.45257   
lag(resid)_7 -0.0647599  0.3043143 -0.2128  0.83148   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 8

data:  reg1
LM test = 18.336, df = 8, p-value = 0.01884


z test of coefficients:

              Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.517747   0.799428 -0.6476  0.51721   
log(I)       -0.023717   0.135574 -0.1749  0.86113   
log(L)       -0.015177   0.100789 -0.1506  0.88030   
log(H)        0.061137   0.115729  0.5283  0.59731   
log(A)        0.079924   0.096834  0.8254  0.40916   
lag(resid)_1  0.693964   0.232355  2.9867  0.00282 **
lag(resid)_2 -0.636719   0.307800 -2.0686  0.03858 * 
lag(resid)_3 -0.015259   0.350882 -0.0435  0.96531   
lag(resid)_4 -0.077449   0.351675 -0.2202  0.82569   
lag(resid)_5  0.009218   0.364608  0.0253  0.97983   
lag(resid)_6 -0.576040   0.360783 -1.5966  0.11035   
lag(resid)_7  0.328776   0.338224  0.9721  0.33102   
lag(resid)_8 -0.550453   0.265757 -2.0713  0.03833 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 9

data:  reg1
LM test = 18.609, df = 9, p-value = 0.02873


z test of coefficients:

               Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.5451522  0.8155430 -0.6685  0.50384  
log(I)       -0.0498137  0.1443954 -0.3450  0.73011  
log(L)        0.0047342  0.1075907  0.0440  0.96490  
log(H)        0.0616000  0.1178905  0.5225  0.60131  
log(A)        0.0860577  0.0991364  0.8681  0.38535  
lag(resid)_1  0.6321378  0.2569036  2.4606  0.01387 *
lag(resid)_2 -0.5884319  0.3231015 -1.8212  0.06858 .
lag(resid)_3 -0.0669995  0.3670720 -0.1825  0.85517  
lag(resid)_4 -0.1115638  0.3624502 -0.3078  0.75823  
lag(resid)_5  0.0250742  0.3722918  0.0674  0.94630  
lag(resid)_6 -0.5713833  0.3675895 -1.5544  0.12009  
lag(resid)_7  0.2008712  0.4017535  0.5000  0.61708  
lag(resid)_8 -0.3965798  0.3675454 -1.0790  0.28059  
lag(resid)_9 -0.1930864  0.3119581 -0.6189  0.53595  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 10

data:  reg1
LM test = 19.959, df = 10, p-value = 0.02964


z test of coefficients:

                Estimate Std. Error z value Pr(>|z|)  
(Intercept)   -0.1550335  0.8371588 -0.1852  0.85308  
log(I)         0.0014658  0.1445965  0.0101  0.99191  
log(L)        -0.0332259  0.1076969 -0.3085  0.75769  
log(H)         0.0044209  0.1211978  0.0365  0.97090  
log(A)         0.0941486  0.0962978  0.9777  0.32823  
lag(resid)_1   0.5776427  0.2520487  2.2918  0.02192 *
lag(resid)_2  -0.7632528  0.3366175 -2.2674  0.02336 *
lag(resid)_3   0.0357213  0.3632121  0.0983  0.92166  
lag(resid)_4  -0.2293107  0.3611031 -0.6350  0.52541  
lag(resid)_5  -0.0909405  0.3701265 -0.2457  0.80591  
lag(resid)_6  -0.5300646  0.3576243 -1.4822  0.13829  
lag(resid)_7   0.1740704  0.3900232  0.4463  0.65537  
lag(resid)_8  -0.7054087  0.4175043 -1.6896  0.09111 .
lag(resid)_9   0.1494190  0.3868736  0.3862  0.69933  
lag(resid)_10 -0.4540255  0.3197091 -1.4201  0.15557  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 11

data:  reg1
LM test = 19.96, df = 11, p-value = 0.04589


z test of coefficients:

                Estimate Std. Error z value Pr(>|z|)  
(Intercept)   -0.1434176  0.9364427 -0.1532  0.87828  
log(I)         0.0039303  0.1675589  0.0235  0.98129  
log(L)        -0.0348233  0.1216991 -0.2861  0.77477  
log(H)         0.0029415  0.1333501  0.0221  0.98240  
log(A)         0.0934710  0.1018034  0.9182  0.35854  
lag(resid)_1   0.5743952  0.2791354  2.0578  0.03961 *
lag(resid)_2  -0.7664678  0.3620142 -2.1172  0.03424 *
lag(resid)_3   0.0328263  0.3862216  0.0850  0.93227  
lag(resid)_4  -0.2275466  0.3776334 -0.6026  0.54680  
lag(resid)_5  -0.0969602  0.4250064 -0.2281  0.81954  
lag(resid)_6  -0.5296179  0.3704140 -1.4298  0.15277  
lag(resid)_7   0.1711476  0.4134661  0.4139  0.67892  
lag(resid)_8  -0.7072013  0.4356020 -1.6235  0.10448  
lag(resid)_9   0.1423660  0.4547913  0.3130  0.75425  
lag(resid)_10 -0.4483524  0.3736064 -1.2001  0.23011  
lag(resid)_11 -0.0138925  0.4246841 -0.0327  0.97390  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


    Breusch-Godfrey test for serial correlation of order up to 12

data:  reg1
LM test = 21.774, df = 12, p-value = 0.04013


z test of coefficients:

               Estimate Std. Error z value Pr(>|z|)   
(Intercept)    0.283118   0.914993  0.3094 0.757001   
log(I)         0.078832   0.163492  0.4822 0.629680   
log(L)        -0.076886   0.116984 -0.6572 0.511031   
log(H)        -0.053007   0.129545 -0.4092 0.682409   
log(A)         0.063121   0.097293  0.6488 0.516486   
lag(resid)_1   0.543726   0.262826  2.0688 0.038568 * 
lag(resid)_2  -0.934762   0.354281 -2.6385 0.008328 **
lag(resid)_3  -0.043661   0.365592 -0.1194 0.904939   
lag(resid)_4  -0.384382   0.366618 -1.0485 0.294430   
lag(resid)_5  -0.096566   0.399222 -0.2419 0.808869   
lag(resid)_6  -0.708319   0.363597 -1.9481 0.051404 . 
lag(resid)_7   0.154905   0.388500  0.3987 0.690096   
lag(resid)_8  -0.902708   0.425155 -2.1232 0.033733 * 
lag(resid)_9   0.201572   0.428629  0.4703 0.638160   
lag(resid)_10 -0.912078   0.445163 -2.0489 0.040475 * 
lag(resid)_11  0.263836   0.431326  0.6117 0.540746   
lag(resid)_12 -0.630576   0.372424 -1.6932 0.090424 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(d)

	bgtest.chi.statistic	bgtest.chi.p.value
LM test	8.840634	0.0029460
LM test1	13.506831	0.0011669
LM test2	13.706948	0.0033324
LM test3	13.754136	0.0081225
LM test4	14.273477	0.0139625
LM test5	15.356109	0.0176607
LM test6	15.392859	0.0312806
LM test7	18.336310	0.0188417
LM test8	18.609052	0.0287300
LM test9	19.959052	0.0296424
LM test10	19.959820	0.0458948
LM test11	21.773878	0.0401336

Olhando as significâncias dos termos defasados dos resíduos em cada teste, é possível concluir por uma estabilidade na relação entre os resíduos de lags 1 e 2 em todos os casos. Em alguns casos o lag 8 também aparece significativo. Ou seja, existem evidências razoáveis para autocorrelação de 2 lags ou AR(2).

Faremos a correção para dois lags pelo método de Newey-West, uma vez que os métodos de Cochrane-Orcutt e Prais-Winsten não corrigem para segunda ordem, e utilizaremos o pacote sandwich com a função NeweyWest.

# library(sandwich) matriz de var-cov robusta
library(sandwich)
# Using NeweyWest():
NW_VCOV <- NeweyWest(reg1, lag = 2, prewhite = F, adjust = T)
NW_VCOV

            (Intercept)       log(I)       log(L)       log(H)       log(A)
(Intercept)  1.68584254  0.032811653  0.028316546 -0.225314059 -0.104834655
log(I)       0.03281165  0.036798278 -0.014441301 -0.009310409 -0.012782349
log(L)       0.02831655 -0.014441301  0.011322445 -0.003016851 -0.002440085
log(H)      -0.22531406 -0.009310409 -0.003016851  0.031781861  0.015524529
log(A)      -0.10483466 -0.012782349 -0.002440085  0.015524529  0.018305041

Podemos comparar esta saída robusta com a regressão inicial:

# fazer stargazer
library(sandwich)
library(stargazer)
robust <- coeftest(reg1, vcov = NW_VCOV)
stargazer(list(reg1, reg1), se = list(NULL, robust[, 2]), p = list(NULL, robust[, 
    4]), p.auto = F, column.labels = c("MQO-reg1", "reg1 robusto HAC"), 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(C)           
              MQO-reg1   reg1 robusto HAC
                 (1)           (2)       
-----------------------------------------
log(I)        0.468***       0.468**     
               (0.166)       (0.192)     
              t = 2.817     t = 2.437    
              p = 0.010     p = 0.023    
log(L)         0.279**       0.279**     
               (0.115)       (0.106)     
              t = 2.436     t = 2.626    
              p = 0.023     p = 0.015    
log(H)         -0.005         -0.005     
               (0.143)       (0.178)     
             t = -0.036     t = -0.029   
              p = 0.972     p = 0.978    
log(A)        0.441***       0.441***    
               (0.107)       (0.135)     
              t = 4.145     t = 3.263    
             p = 0.0004     p = 0.004    
Constant       -1.500         -1.500     
               (1.003)       (1.298)     
             t = -1.496     t = -1.156   
              p = 0.148     p = 0.259    
-----------------------------------------
Observations     30             30       
R2              0.936         0.936      
Adjusted R2     0.926         0.926      
=========================================
Note:         *p<0.1; **p<0.05; ***p<0.01

Ou fazendo pelo coeftest para a saída robusta de Newey-West (atentar para a diferença entre o que significa cada asterisco - ver notas dos códigos de significância):

coeftest(reg1, NW_VCOV)

t test of coefficients:

              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -1.5004409  1.2984000 -1.1556 0.258767   
log(I)       0.4675085  0.1918288  2.4371 0.022260 * 
log(L)       0.2794425  0.1064070  2.6262 0.014530 * 
log(H)      -0.0051516  0.1782747 -0.0289 0.977176   
log(A)       0.4414489  0.1352961  3.2628 0.003184 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3 Referências

GUJARATI, Damodar N.; PORTER, Dawn C. Econometria básica. 5.ed. Porto Alegre: MGH/Bookman/McGraw-Hill do Brasil, 2011.

GUJARATI, Damodar N. Basic Econometrics. 4th edition. The McGraw−Hill Companies, 2004.