Econometria: Tabela 5.6 PNB e demanda de moeda, Gujarati e Porter (2011, p.157)

Abstract

We analyse different money supply concepts and its relation to the GDP. The example is draw following Gujarati and Porter (2011, p.157).

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

1 Primeiros passos

Esse post é o mesmo que uso no video do Youtube em: https://youtu.be/lpni-E1PAbU, ECONOMETRIA 2020.2 Slides 02a Instalação R e RStudio.

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’.Neste caso, a planilha chama-se Tabela 5.6aula_gujarati 5ed p157.xlsx.

A Tabela 5.6 apresenta dados sobre o PNB e quatro definições do estoque de moeda dos Estados Unidos no período 1970-1983. Fazendo as regressões do PNB contra as várias definições de moeda, obtemos os resultados apresentados na Tabela 5.7. Os monetaristas ou adeptos da teoria quantitativa afirmam que a renda nominal (o PNB nominal) é determinada, em grande parte, pela variação na quantidade ou estoque de moeda, embora não haja consenso quanto à definição “certa” de moeda. De acordo com os resultados, responda às seguintes perguntas:
a. Que definição de moeda parece apresentar relação mais estreita com o PNB nominal?
b. Como os \(R^2\) são todos elevados, isso significa que a escolha da definição de moeda não tem importância?
c. Se o FED (Banco Central dos Estados Unidos) quer controlar a oferta de moeda, qual desses indicadores de moeda seria o melhor objetivo para esse fim? Isso pode ser dito com base nos resultados da regressão?

Exercício 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.

Coloco aqui a opção de download do xlsx mas coloquei os dados embeded no code para quem quiser. É só executar o chunk a seguir que conseguem carregar os dados.

> 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"))

> library(readxl)
> library(foreign)
> dados <- read_excel("Tabela 5.6aula_gujarati 5ed p157.xlsx", sheet = "dados")
> attach(dados)

A opção do github pode ser obtida no chunk abaixo. Também pode ser feito abrindo o github em “https://github.com/amrofi/Gujarati-and-Porter-dataset/raw/master/table5-6.wf1.rds” e fazer o download manual (pelo botão de download, ou clique com o botao da direita do mouse sobre download e ‘save as’ no local do seu projeto de pesquisa). Feito o download, abra o arquivo fazendo uso do comando readRDS(“table5-6.wf1.rds”) como está no rpubs.

> library(curl)
> githubURL <- ("https://github.com/amrofi/Gujarati-and-Porter-dataset/raw/master/table5-6.wf1.rds")
> download.file(githubURL, "table5-6.wf1.rds", method = "auto")
> dados <- readRDS("table5-6.wf1.rds")

Vamos ver como está a tabela importada:

> library(knitr)
> kable(dados, caption = "Dados da tabela 5.6")

Dados da tabela 5.6

ANO	PNB	M1	M2	M3	L
1970	992.7	216.6	628.2	677.5	816.3
1971	1077.6	230.8	712.8	776.2	903.1
1972	1185.9	252.0	805.2	886.0	1023.0
1973	1326.4	265.9	861.0	985.0	1141.7
1974	1434.2	277.6	908.5	1070.5	1249.3
1975	1549.2	291.2	1023.3	1174.2	1367.9
1976	1718.0	310.4	1163.6	1311.9	1516.6
1977	1918.3	335.4	1286.7	1472.9	1704.7
1978	2163.9	363.1	1389.1	1647.1	1910.6
1979	2417.8	389.1	1498.5	1804.8	2117.1
1980	2631.7	414.9	1632.6	1990.0	2326.2
1981	2957.8	441.9	1796.6	2238.2	2599.8
1982	3069.3	480.5	1965.4	2462.5	2870.8
1983	3304.8	525.4	2196.3	2710.4	3183.1

2 Resultados

2.1 Estimação

Vamos estimar o modelo linear multipla fazendo as regressoes do PNB contra as varias definicoes de moeda.

> EQ1 <- lm(PNB ~ M1)
> EQ2 <- lm(PNB ~ M2)
> EQ3 <- lm(PNB ~ M3)
> EQ4 <- lm(PNB ~ L)

2.2 Saída gráfica

[1] 1 2 3 4

2.3 Saída tabular

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

> summary(EQ1)

Call:
lm(formula = PNB ~ M1)

Residuals:
     Min       1Q   Median       3Q      Max 
-156.272  -27.294   -5.445   25.318  171.934 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -787.4724    77.9665   -10.1 3.22e-07 ***
M1             8.0863     0.2197    36.8 1.04e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 76.28 on 12 degrees of freedom
Multiple R-squared:  0.9912,    Adjusted R-squared:  0.9905 
F-statistic:  1354 on 1 and 12 DF,  p-value: 1.039e-13

Agora, criando uma tabela com as várias saídas de modelos, com o pacote stargazer tem-se:

==========================================================================
                                          Dependent variable:             
                              --------------------------------------------
                                                  PNB                     
                                  (1)        (2)        (3)        (4)    
--------------------------------------------------------------------------
M1                             8.086***                                   
                                (0.220)                                   
                              t = 36.801                                  
                               p = 0.000                                  
M2                                         1.587***                       
                                           (0.045)                        
                                          t = 35.409                      
                                          p = 0.000                       
M3                                                    1.204***            
                                                      (0.026)             
                                                     t = 45.817           
                                                     p = 0.000            
L                                                                1.029*** 
                                                                 (0.023)  
                                                                t = 43.937
                                                                p = 0.000 
Constant                      -787.472***  -44.063   158.788*** 164.207***
                               (77.966)    (61.013)   (43.041)   (44.766) 
                              t = -10.100 t = -0.722 t = 3.689  t = 3.668 
                              p = 0.00000 p = 0.485  p = 0.004  p = 0.004 
--------------------------------------------------------------------------
Observations                      14          14         14         14    
R2                               0.991      0.991      0.994      0.994   
Adjusted R2                      0.990      0.990      0.994      0.993   
Residual Std. Error (df = 12)   76.282      79.255     61.368     63.977  
==========================================================================
Note:                                          *p<0.1; **p<0.05; ***p<0.01

Vou usar o pacote effects para ver os efeitos de cada uma:

> library(effects)
> # usando o pacote effects
> alleff1 <- allEffects(EQ1)
> plot(alleff1)

> summary(alleff1)

 model: PNB ~ M1

 M1 effect
M1
      220       290       370       450       530 
 991.5146 1557.5560 2204.4603 2851.3647 3498.2690 

 Lower 95 Percent Confidence Limits
M1
      220       290       370       450       530 
 917.9501 1506.5214 2158.1283 2783.3756 3398.1086 

 Upper 95 Percent Confidence Limits
M1
     220      290      370      450      530 
1065.079 1608.591 2250.792 2919.354 3598.430 

> alleff2 <- allEffects(EQ2)
> plot(alleff2)

> summary(alleff2)

 model: PNB ~ M2

 M2 effect
M2
      630      1000      1400      1800      2200 
 956.0393 1543.4007 2178.3860 2813.3713 3448.3566 

 Lower 95 Percent Confidence Limits
M2
      630      1000      1400      1800      2200 
 877.8396 1489.9386 2130.6787 2744.4718 3347.0075 

 Upper 95 Percent Confidence Limits
M2
     630     1000     1400     1800     2200 
1034.239 1596.863 2226.093 2882.271 3549.706 

> alleff3 <- allEffects(EQ3)
> plot(alleff3)

> summary(alleff3)

 model: PNB ~ M3

 M3 effect
M3
      680      1200      1700      2200      2700 
 977.2223 1603.0843 2204.8746 2806.6649 3408.4552 

 Lower 95 Percent Confidence Limits
M3
      680      1200      1700      2200      2700 
 917.5564 1563.0635 2167.6002 2753.6075 3331.7819 

 Upper 95 Percent Confidence Limits
M3
     680     1200     1700     2200     2700 
1036.888 1643.105 2242.149 2859.722 3485.128 

> alleff4 <- allEffects(EQ4)
> plot(alleff4)

> summary(alleff4)

 model: PNB ~ L

 L effect
L
     820     1400     2000     2600     3200 
1008.031 1604.882 2222.314 2839.746 3457.179 

 Lower 95 Percent Confidence Limits
L
      820      1400      2000      2600      3200 
 947.0352 1563.1979 2183.1996 2783.2023 3375.0845 

 Upper 95 Percent Confidence Limits
L
     820     1400     2000     2600     3200 
1069.027 1646.566 2261.429 2896.291 3539.273 

2.4 Teste de autocorrelação

Primeiro farei a definição da série temporal.

> library(datasets)
> library(graphics)
> library(lmtest)
> library(stats)
> # como a base é anual, faz-se
> dados.st <- ts(dados, start = c(1970), end = c(1983))
> PNB.ts <- ts(PNB, start = c(1970), end = c(1983), frequency = 1)
> 
> summary(dados.st)

      ANO            PNB               M1              M2        
 Min.   :1970   Min.   : 992.7   Min.   :216.6   Min.   : 628.2  
 1st Qu.:1973   1st Qu.:1353.3   1st Qu.:268.8   1st Qu.: 872.9  
 Median :1976   Median :1818.2   Median :322.9   Median :1225.2  
 Mean   :1976   Mean   :1982.0   Mean   :342.5   Mean   :1276.3  
 3rd Qu.:1980   3rd Qu.:2578.2   3rd Qu.:408.4   3rd Qu.:1599.1  
 Max.   :1983   Max.   :3304.8   Max.   :525.4   Max.   :2196.3  
       M3               L         
 Min.   : 677.5   Min.   : 816.3  
 1st Qu.:1006.4   1st Qu.:1168.6  
 Median :1392.4   Median :1610.7  
 Mean   :1514.8   Mean   :1766.4  
 3rd Qu.:1943.7   3rd Qu.:2273.9  
 Max.   :2710.4   Max.   :3183.1  

> plot(dados.st)

> 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

> # analise grafica da autocorrelacao
> plot(regressao1)

2.4.1 Teste de Durbin-Watson

> dwtest(regressao1)

    Durbin-Watson test

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

2.4.2 Teste de Breusch-Godfrey para ordens superiores

Também chamado de LM de correlação serial

> 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", "F"), data = dados.st)
> bgteste4 <- bgtest(regressao1, order = 4, type = c("Chisq", "F"), data = dados.st)
> bgteste5 <- bgtest(regressao1, order = 5, type = c("Chisq", "F"), data = dados.st)
> bgteste6 <- bgtest(regressao1, order = 6, type = c("Chisq", "F"), 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

Fazendo a regressão auxiliar do LM (BG) de correlação serial para 1 lag.

> 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

2.4.3 correcao por cochrane-orcutt

> library(dynlm)
> library(car)
> library(orcutt)
> library(sandwich)
> olsres <- dynlm(PNB ~ M1 + M2 + M3 + L, data = dados.st)
> # Cochrane-Orcutt estimation (para primeira ordem)
> cochrane.orcutt(olsres)

Cochrane-orcutt estimation for first order autocorrelation 
 
Call:
dynlm(formula = PNB ~ M1 + M2 + M3 + L, data = dados.st)

 number of interaction: 100
 rho 0.828336

Durbin-Watson statistic 
(original):    0.70426 , p-value: 8.024e-05
(transformed): NA , p-value: NA
 
 coefficients: 
[1] NA

2.4.4 Correcao por Newey-West - Inferência Robusta

> # resultados com erros-padrões usuais
> coeftest(olsres)

t test of coefficients:

              Estimate Std. Error t value Pr(>|t|)
(Intercept) 913.488117 720.112663  1.2685   0.2364
M1           -6.014155   6.400161 -0.9397   0.3719
M2           -0.364738   0.878230 -0.4153   0.6876
M3            2.336460   1.989021  1.1747   0.2703
L             0.030841   1.687079  0.0183   0.9858

> # resultados com Erros-padrões HAC - robusto a heterocedasticidade e
> # autocorrelação
> coeftest(olsres, vcovHAC)

t test of coefficients:

              Estimate Std. Error t value Pr(>|t|)  
(Intercept) 913.488117 320.349897  2.8515  0.01904 *
M1           -6.014155   2.609987 -2.3043  0.04667 *
M2           -0.364738   0.152061 -2.3986  0.03999 *
M3            2.336460   1.042247  2.2418  0.05170 .
L             0.030841   0.657567  0.0469  0.96362  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1