Econometria: exemplo lanchonetes de Hill et al (2008)

Abstract

This is an undergrad student level exercise for class use.

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: exemplo de lanchonetes de Hill et al (2008). Campo Grande-MS,Brasil: RStudio/Rpubs, 2020. Disponível em http://rpubs.com/amrofi/estimation_lanchonetes.

1 Introdução

Exercício baseado em rotina de Colonescu (2016) - livro do bookdown em https://bookdown.org/ccolonescu/RPoE4/multiplereg.html, baseado em exercício de Hill et al (2008).

1.1 Pacotes necessários

# limpando area de trabalho
rm(list = ls())
# chamando os pacotes necessários library(PoEdata) # não está instalando na
# versao R.4.0
library(knitr)
library(xtable)
library(printr)
library(effects)
library(car)
library(AER)
library(broom)
library(tibble)

2 Big Andy’s Hamburger Sales: 5.2 Example:

Para realização deste exercício, como estão ocorrendo erros para instalar a partir do git e não está no CRAN, recomendo fazer download do dataset a partir de https://github.com/ccolonescu/PoEdata/tree/master/data diretamente, invés de instalar e chamar o dataset. Colocarei embedded no arquivo do code (.Rmd). Se voce baixou o arquivo no diretorio de trabalho, pode chamar fazendo load("andy.rda").

Uma opção que funcionou para mim foi:

githubURL <- "https://github.com/ccolonescu/PoEdata/raw/master/data/andy.rda"
download.file(githubURL, "myfile")
load("myfile")

Coloquei a opção do dput abaixo, mas também pode-se usar o pacote mais atual do livro, que se chama POE5Rdata e pode ser instalado do github:

devtools::install_github("ccolonescu/POE5Rdata")
library(POE5Rdata)  #activates the package 'POE5Rdata'
data(andy)  #reads the dataset `andy` into the global environment

andy <- structure(list(sales = c(73.2, 71.8, 62.4, 67.4, 89.3, 70.3, 73.2, 86.1,
    81, 76.4, 76.6, 82.2, 82.1, 68.6, 76.5, 80.3, 70.7, 75, 73.7, 71.2, 84.7, 73.6,
    73.7, 78.1, 75.7, 74.4, 68.7, 83.9, 86.1, 73.7, 75.7, 78.8, 73.7, 80.2, 69.9,
    69.1, 83.8, 84.3, 66, 84.3, 79.5, 80.2, 67.6, 86.5, 87.6, 84.2, 75.2, 84.7, 73.7,
    81.2, 69, 69.7, 78.1, 88, 80.4, 79.7, 73.2, 85.9, 83.3, 73.6, 79.2, 88.1, 64.5,
    84.1, 91.2, 71.8, 80.6, 73.1, 81, 73.7, 82.2, 74.2, 75.4, 81.3, 75), price = c(5.69,
    6.49, 5.63, 6.22, 5.02, 6.41, 5.85, 5.41, 6.24, 6.2, 5.48, 6.14, 5.37, 6.45,
    5.35, 5.22, 5.89, 5.21, 6, 6.37, 5.33, 5.23, 5.88, 6.24, 5.59, 6.22, 6.41, 4.96,
    4.83, 6.35, 6.47, 5.69, 5.56, 6.41, 5.54, 6.47, 4.94, 6.16, 5.93, 5.2, 5.62,
    5.28, 5.46, 5.11, 5.04, 5.08, 5.86, 4.89, 5.68, 5.83, 6.33, 6.47, 5.7, 5.22,
    5.05, 5.76, 6.25, 5.34, 4.98, 6.39, 6.22, 5.1, 6.49, 4.86, 5.1, 5.98, 5.02, 5.08,
    5.23, 6.02, 5.73, 5.11, 5.71, 5.45, 6.05), advert = c(1.3, 2.9, 0.8, 0.7, 1.5,
    1.3, 1.8, 2.4, 0.7, 3, 2.8, 2.7, 2.8, 2.8, 2.3, 1.7, 1.5, 0.8, 2.9, 0.5, 2.1,
    0.8, 1.1, 1.9, 2.1, 1.3, 1.1, 1.1, 2.9, 1.4, 2.5, 3, 1, 3.1, 0.5, 2.7, 0.9, 1.5,
    2.8, 2.3, 1.2, 3.1, 1, 2.5, 2.1, 2.8, 3.1, 3.1, 0.9, 1.8, 3.1, 1.9, 0.7, 1.6,
    2.9, 2.3, 1.7, 1.8, 0.6, 3.1, 1.2, 2.1, 0.5, 2.9, 1.6, 1.5, 2, 1.3, 1.1, 2.2,
    1.7, 0.7, 0.7, 2, 2.2)), datalabel = "", time.stamp = "", formats = c("%10.0g",
    "%10.0g", "%10.0g"), types = c(255L, 255L, 255L), val.labels = c("", "", ""),
    var.labels = c("Monthly sales revenue ($1000s)", "A price index for all products sold in a given month.",
        "Expenditure on advertising ($1000s)"), row.names = c("1", "2", "3", "4",
        "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17",
        "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30",
        "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43",
        "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56",
        "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", "69",
        "70", "71", "72", "73", "74", "75"), version = 12L, class = "data.frame")

# Resumo estatístico
s = tidy(andy)[, c(1:5, 8, 9)]
s

column	n	mean	sd	median	min	max
sales	75	77.37467	6.4885365	76.50	62.40	91.20
price	75	5.68720	0.5184320	5.69	4.83	6.49
advert	75	1.84400	0.8316769	1.80	0.50	3.10

2.1 O modelo andy básico de regressão múltipla

mod1 <- lm(sales ~ price + advert, data = andy)
smod1 <- data.frame(xtable(summary(mod1)))
kable(smod1, caption = " O modelo básico de regressão múltipla", col.names = c("Coeficiente",
    "Erro padrão", "t-valor", "p-valor"), align = "c", digits = 3)

O modelo básico de regressão múltipla

	Coeficiente	Erro padrão	t-valor	p-valor
(Intercept)	118.914	6.352	18.722	0.000
price	-7.908	1.096	-7.215	0.000
advert	1.863	0.683	2.726	0.008

summary(mod1)

Call:
lm(formula = sales ~ price + advert, data = andy)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.4825  -3.1434  -0.3456   2.8754  11.3049 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 118.9136     6.3516  18.722  < 2e-16 ***
price        -7.9079     1.0960  -7.215 4.42e-10 ***
advert        1.8626     0.6832   2.726  0.00804 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.886 on 72 degrees of freedom
Multiple R-squared:  0.4483,    Adjusted R-squared:  0.4329 
F-statistic: 29.25 on 2 and 72 DF,  p-value: 5.041e-10

2.1.1 Efeitos

# usando o pacote effects
alleffandy <- allEffects(mod1)
plot(alleffandy)

plot(predictorEffects(mod1))

summary(alleffandy)

 model: sales ~ price + advert

 price effect
price
     4.8      5.2      5.7      6.1      6.5 
84.39052 81.22737 77.27345 74.11030 70.94716 

 Lower 95 Percent Confidence Limits
price
     4.8      5.2      5.7      6.1      6.5 
82.14947 79.67882 76.14838 72.66864 68.84513 

 Upper 95 Percent Confidence Limits
price
     4.8      5.2      5.7      6.1      6.5 
86.63156 82.77593 78.39851 75.55197 73.04919 

 advert effect
advert
     0.5      1.1      1.8      2.4      3.1 
74.87135 75.98890 77.29271 78.41026 79.71407 

 Lower 95 Percent Confidence Limits
advert
     0.5      1.1      1.8      2.4      3.1 
72.72299 74.47507 76.16640 77.05439 77.66686 

 Upper 95 Percent Confidence Limits
advert
     0.5      1.1      1.8      2.4      3.1 
77.01971 77.50274 78.41902 79.76613 81.76128 

2.2 Modelo quadrático de regressão múltipla

# Outro exemplo
mod2 <- lm(sales ~ price + advert + I(advert^2), data = andy)
summary(mod2)

Call:
lm(formula = sales ~ price + advert + I(advert^2), data = andy)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.2553  -3.1430  -0.0117   2.8513  11.8050 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 109.7190     6.7990  16.137  < 2e-16 ***
price        -7.6400     1.0459  -7.304 3.24e-10 ***
advert       12.1512     3.5562   3.417  0.00105 ** 
I(advert^2)  -2.7680     0.9406  -2.943  0.00439 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.645 on 71 degrees of freedom
Multiple R-squared:  0.5082,    Adjusted R-squared:  0.4875 
F-statistic: 24.46 on 3 and 71 DF,  p-value: 5.6e-11

plot(effect("I(advert^2)", mod2))

2.2.1 Extraindo os estimadores e intervalo de confiança

Em geral, estas informações são obtidas diretamente do summary, mas aqui o leitor pode ver como extraí-las do objeto mod1:

# grau de liberdade do modelo 1
df <- mod1$df.residual
N <- nobs(mod1)
b1 <- coef(mod1)[[1]]  #or
b2 <- coef(mod1)[["price"]]
b3 <- coef(mod1)[["advert"]]
sighat2 <- smod1$sigma^2
anov <- anova(mod1)
SSE <- anov[3, 2]  # soma dos quadrados dos residuos SQRes
SST <- sum(anov[, 2])  #sum of column 2 in anova - SQTot
SSR <- SST - SSE  # soma dos quadrados da regressao - SQReg
kable(data.frame(vcov(mod1)), align = "c", digits = 3, caption = "A matriz de covariância dos coeficientes",
    col.names = c("(Intercept)", "price", "advert"))

A matriz de covariância dos coeficientes

	(Intercept)	price	advert
(Intercept)	40.343	-6.795	-0.748
price	-6.795	1.201	-0.020
advert	-0.748	-0.020	0.467

2.2.1.1 extraindo erros-padrões dos coeficientes

### extraindo erros-padrões dos coeficientes
varb1 <- vcov(mod1)[1, 1]
varb2 <- vcov(mod1)[2, 2]
varb3 <- vcov(mod1)[3, 3]
covb1b2 <- vcov(mod1)[1, 2]
covb1b3 <- vcov(mod1)[1, 3]
covb2b3 <- vcov(mod1)[2, 3]
seb2 <- sqrt(varb2)  #standard error of b2
seb3 <- sqrt(varb3)

2.2.1.2 intervalo de confiança para os parâmetros

# intervalo de confiança para os parâmetros
alpha <- 0.05  # escolhi 5% de significancia
tcr <- qt(1 - alpha/2, df)  # t critico
lowb2 <- b2 - tcr * seb2
upb2 <- b2 + tcr * seb2
lowb3 <- b3 - tcr * seb3
upb3 <- b3 + tcr * seb3
confints <- confint(mod1, parm = c("price", "advert"), level = 0.95)
kable(data.frame(confints), caption = "Confidence intervals for 'price' and 'advert'",
    align = "c", col.names = c("lowb", "upb"), digits = 4)

Confidence intervals for ‘price’ and ‘advert’

	lowb	upb
price	-10.0927	-5.7230
advert	0.5007	3.2245

2.2.1.3 Intervalos de confiança com confint

summary(mod1)

Call:
lm(formula = sales ~ price + advert, data = andy)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.4825  -3.1434  -0.3456   2.8754  11.3049 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 118.9136     6.3516  18.722  < 2e-16 ***
price        -7.9079     1.0960  -7.215 4.42e-10 ***
advert        1.8626     0.6832   2.726  0.00804 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.886 on 72 degrees of freedom
Multiple R-squared:  0.4483,    Adjusted R-squared:  0.4329 
F-statistic: 29.25 on 2 and 72 DF,  p-value: 5.041e-10

confint(mod1, level = 0.95)

	2.5 %	97.5 %
(Intercept)	106.251852	131.575368
price	-10.092676	-5.723032
advert	0.500659	3.224510

2.2.1.4 opção do pacote jtools

O pacote jtools traz algumas facilidades para o leitor que desejar ganhar tempo. Ver a vinheta em https://cran.r-project.org/web/packages/jtools/vignettes/summ.html

library(jtools)
summ(mod1, confint = TRUE, digits = 3)

Observations	75
Dependent variable	sales
Type	OLS linear regression

F(2,72)	29.248
R²	0.448
Adj. R²	0.433

	Est.	2.5%	97.5%	t val.	p
(Intercept)	118.914	106.252	131.575	18.722	0.000
price	-7.908	-10.093	-5.723	-7.215	0.000
advert	1.863	0.501	3.225	2.726	0.008
Standard errors: OLS

2.2.2 intervalo de confianca para uma combinacao de parametros

Para sales quando price decresce 40 cents e advert aumenta $800.

# na equacao andy basica
alpha <- 0.1  # neste caso, assumi 10% de significancia
tcr <- qt(1 - alpha/2, df)
a1 <- 0
a2 <- -0.4  # price decresce 40  cents
a3 <- 0.8  # advert aumenta $800
L <- a1 * b1 + a2 * b2 + a3 * b3  # L = 0*beta1-0.4*beta2+0.8*beta3
# variancia da combinacao linear
varL <- a1^2 * varb1 + a2^2 * varb2 + a3^2 * varb3 + 2 * a1 * a2 * covb1b2 + 2 *
    a1 * a3 * covb1b3 + 2 * a2 * a3 * covb2b3
# desvio-padrao da combinacao linear
seL <- sqrt(varL)
# intervalos inferior 'low' e superior 'upper'
lowbL <- L - tcr * seL
upbL <- L + tcr * seL
c(lowbL, upbL)

[1] 3.470785 5.835633

# (3.471,5.836)

Nesse cenario as vendas oscilarão entre 3.471 e 5.836.

2.2.3 outra forma mais rápida e direta

# vetor para price decresce 40 cents e advert aumenta $800
a <- c(0, -0.4, 0.8)
b <- as.numeric(coef(mod1))  # vetor de coeficientes
L <- L <- sum(a * b)  # soma dos produtos elementwise
V <- vcov(mod1)  # a matriz variancia-covariancia
A <- as.vector(a)  # (de fato, desnecessario)
varL <- as.numeric(t(A) %*% V %*% A)
seL <- sqrt(varL)
lowbL <- L - tcr * seL
upbL <- L + tcr * seL
c(lowbL, upbL)

[1] 3.470785 5.835633

2.3 Modelo log-log

modlog <- lm(log(sales) ~ log(price) + log(advert), data = andy)
summary(modlog)

Call:
lm(formula = log(sales) ~ log(price) + log(advert), data = andy)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.182691 -0.038026 -0.005745  0.040823  0.143404 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.31995    0.13795  38.564  < 2e-16 ***
log(price)  -0.57494    0.07935  -7.246 3.88e-10 ***
log(advert)  0.04544    0.01366   3.327  0.00139 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06237 on 72 degrees of freedom
Multiple R-squared:  0.4691,    Adjusted R-squared:  0.4544 
F-statistic: 31.81 on 2 and 72 DF,  p-value: 1.26e-10

2.4 Modelo quadrático em price e em advert

mod4 <- lm(sales ~ price + I(price^2) + advert + I(advert^2), data = andy)
summary(mod4)

Call:
lm(formula = sales ~ price + I(price^2) + advert + I(advert^2), 
    data = andy)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.1771  -2.9378   0.0949   2.8897  11.6856 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 248.8390    81.5712   3.051 0.003225 ** 
price       -57.0629    28.8988  -1.975 0.052262 .  
I(price^2)    4.3364     2.5340   1.711 0.091454 .  
advert       13.1624     3.5582   3.699 0.000427 ***
I(advert^2)  -3.0896     0.9469  -3.263 0.001708 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.583 on 70 degrees of freedom
Multiple R-squared:  0.528, Adjusted R-squared:  0.501 
F-statistic: 19.57 on 4 and 70 DF,  p-value: 7.553e-11

2.5 Simulacao do intervalo de previsão

Faremos agora uma simulacao do intervalo de previsão de outra forma, com resultados da estimação por lm(). Seja mod1 <- lm(sales~price+advert, data=andy).

Estimamos mod1, em seguida envolvemos os parâmetros dentro de um novo data frame da variável ‘newdata’ com valores desejados das variáveis X. Exemplo, para price de 4.5 e advert de 2.0, quanto seria sales, e limites para 95% de confiança?

Vamos agora aplicar a predict function e definir a variável preditora no argumento ‘newdata’. Também definimos o tipo de intervalo como “confidence” e usamos o nível de confiança de 0.95 (alpha=0.05):

mod1 <- lm(sales ~ price + advert, data = andy)
newdata = data.frame(price = 4.5, advert = 2)
predict(mod1, newdata, interval = "confidence", level = 0.99)

fit	lwr	upr
87.05343	83.28366	90.82321

O resultado seria sales de 87.05, variando entre 84.21 e 89.89.

3 Testando Hipóteses Simultâneas

Agora, seguindo a seção 6.2 do Colonescu (2016), testaremos a hipótese de dois ou mais parâmetros simultaneamente não significativos. Por exemplo, para a equação do modelo quadrático acima estimado,

\[ sales=\beta_{1}+\beta_{2}price+\beta_{3}advert+\beta_{4}advert^2+\varepsilon \]

uma hipótese nula do tipo:

\[ H_{0}:\beta_{3}=0\;\; E \;\;\beta_{4}=0;\;\;\; \\H_{A}:\beta_{3}\neq 0\;\; OU\;\; \beta_{4}\neq 0 \]

Primeiro, observaremos a saída do modelo estimado (mod2):

# mod2 <- lm(sales~price+advert+I(advert^2), data=andy)
summary(mod2)

Call:
lm(formula = sales ~ price + advert + I(advert^2), data = andy)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.2553  -3.1430  -0.0117   2.8513  11.8050 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 109.7190     6.7990  16.137  < 2e-16 ***
price        -7.6400     1.0459  -7.304 3.24e-10 ***
advert       12.1512     3.5562   3.417  0.00105 ** 
I(advert^2)  -2.7680     0.9406  -2.943  0.00439 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.645 on 71 degrees of freedom
Multiple R-squared:  0.5082,    Adjusted R-squared:  0.4875 
F-statistic: 24.46 on 3 and 71 DF,  p-value: 5.6e-11

# plot(effect('I(advert^2)', mod2))

Olhando a saída acima, percebe-se a identificação da significância, mas a questão é saber se realmente \(\beta_3\) e \(\beta_4\) são simultaneamente nulos. Usaremos a função LinearHypothesis para tal.

library(car)
myH0 <- c("advert=0", "I(advert^2)=0")
linearHypothesis(mod2, myH0)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
73	1896.391	NA	NA	NA	NA
71	1532.084	2	364.3064	8.44136	0.0005142

Interpretação do teste: é um teste F cuja expressão irrestrita é \(sales=\beta_{1}+\beta_{2}price+\beta_{3}advert+\beta_{4}advert^2+\varepsilon\) e a restrita é \(sales=\beta_{1}+\beta_{2}price+\varepsilon\). Os resultados indicaram uma estatística \(F=8.44\) e um p-valor \(Pr(>F) = 0.0005\), portanto, muito menor que 1% de significância, indicando a rejeição de H0 com 99% de confiança. Portanto, ou \(\beta_{3}\neq 0\), ou \(\beta_{4}\neq 0\). Como vimos nos resultados, ambos foram estatisticamente diferentes de zero, e portanto, significativos. A omissão de \(advert\) ou \(advert^2\) levarão a viés de omissão de variáveis relevantes.

4 Referências

COLONESCU, C. Principles of Econometrics with R. 2016. https://bookdown.org/ccolonescu/RPoE4/
HILL, R. Carter; GRIFFITHS, William E.; JUDGE, George G. Econometria. 2. ed. São Paulo, SP: Saraiva, 2008. 471 p. ISBN 8502039040. Número de chamada: 330.015195 H647e.2 (CBC)