library(AER)
data("CigarettesSW", package = "AER")
# turn off scientific notation except for big numbers
options(scipen = 9)CigarettesSW: é a base dados utilizada. Cada observação é um estado dos EUA no ano de 1995 (48 no total).Taxar o uso de cigarros é importante pq tem externalidade negativa. Mas determinar o valor do tributo é algo que depende da elasticidade da demanda pelo fumo. Uma demanda mais elástica (-Inf) precisa de menos imposto para reduzir o consumo da diamba.
Vamos começar organizando a base de dados a base de dados:
library(tidyverse)
cigs <- CigarettesSW %>%
filter(year == 1995) %>%
mutate(rprice = price/cpi, # do preço do bagulho (deflacionado)
rincome = (income/population)/cpi, # a renda do usuário
lprice = log(rprice), # log(preço)
lquant = log(packs), # quantidade
lincome = log(rincome), # renda da população
tdiff = (taxs - tax)/cpi) %>% # imposto cobrado por qtde (taxa real)
glimpse()## Observations: 48
## Variables: 15
## $ state <fctr> AL, AR, AZ, CA, CO, CT, DE, FL, GA, IA, ID, IL, IN...
## $ year <fctr> 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 19...
## $ cpi <dbl> 1.524, 1.524, 1.524, 1.524, 1.524, 1.524, 1.524, 1....
## $ population <dbl> 4262731, 2480121, 4306908, 31493524, 3738061, 32652...
## $ packs <dbl> 101.08543, 111.04297, 71.95417, 56.85931, 82.58292,...
## $ income <dbl> 83903280, 45995496, 88870496, 771470144, 92946544, ...
## $ tax <dbl> 40.50000, 55.50000, 65.33333, 61.00000, 44.00000, 7...
## $ price <dbl> 158.3713, 175.5425, 198.6075, 210.5047, 167.3500, 2...
## $ taxs <dbl> 41.90467, 63.85917, 74.79082, 74.77133, 44.00000, 8...
## $ rprice <dbl> 103.91821, 115.18538, 130.31989, 138.12643, 109.809...
## $ rincome <dbl> 12.91535, 12.16907, 13.53964, 16.07359, 16.31556, 2...
## $ lprice <dbl> 4.643604, 4.746543, 4.869992, 4.928169, 4.698749, 4...
## $ lquant <dbl> 4.615966, 4.709917, 4.276029, 4.040580, 4.413803, 4...
## $ lincome <dbl> 2.558416, 2.498898, 2.605622, 2.777178, 2.792119, 3...
## $ tdiff <dbl> 0.9216975, 5.4850193, 6.2057067, 9.0363074, 0.00000...De um modo geral, a elasticidade é calculada como:
\[elasticity = \epsilon = \frac{dQ_d/Q_d}{dP/P} = \frac{dln(Q_d)}{dln(P)}\] Mas as combinações de preços e quantidades refletem tanto os efeitos da demanda como da oferta. E normalmente a oferta e a demanda por um bem são modeladas como um sistema de equações simultaneas que relacionam separadamente a quanditade e o preço de equilíbzZZzZZZzz…. e isso resulta em uma equação reduzida e ‘logaritmizada’ que pode ser estimado através de um modelo de regressão linear, onde \(\widehat{\beta}\) é o parâmetro de estimação da eslasticidade da demanda.
\[lnQ_{d} = \alpha + \beta lnP + u\] A algebrera consegue nos dar um parâmetro razoável para estimar a elasticidade da demanda, mas ela não consegue dissociar a correlação entre \(\widehat{\beta}\) e \(u\), pois este último contém em sua estrutura os efeitos da oferta que podem refletir no preço. Vejamos a regressão por Mínimos Quadrados Ordinários, onde é esperado que o coeficiente \(\widehat{\beta}_{OLS}\) seja viesado:
ols <- lm(lquant ~ lprice, data = cigs)
# function to calculate corrected SEs for OLS regression
cse <- function(reg) {
rob <- sqrt(diag(vcovHC(reg, type = "HC1")))
return(rob)
}
library("stargazer")
stargazer(ols,
se = list(cse(ols)), # os erros são robustos. A função cse() faz esse trabalho.
title = "OLS Regression", type = "text", df = FALSE, digits = 5)##
## OLS Regression
## ===============================================
## Dependent variable:
## ---------------------------
## lquant
## -----------------------------------------------
## lprice -1.21306***
## (0.19459)
##
## Constant 10.33892***
## (0.93482)
##
## -----------------------------------------------
## Observations 48
## R2 0.40575
## Adjusted R2 0.39283
## Residual Std. Error 0.18962
## F Statistic 31.40859***
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01(a.k.a Two Stage Least Squares (TSLS))
A solução para o modelo requer o emprego de variáveis instrumentais, pois o preço (\(P\) ou lprice ou, genericamente \(X\)) é correlacionado com o o termo de erro (\(u\)). Logo, o \(\widehat{\beta}_{OLS}\) estimado é inconsistente ou viesasdo e o que precisamos é de uma variável \(Z\) que seja altamente correlacionada com \(X\) mas não correlacionada com \(u\).
Formalmente, conforme Stock & Watson (2010), o instrumento utilizado deve atender a condição de:
Uma variável instrumental que é relevante e endógena pode capturar movimentos em \(X_{i}\) que são exógenos e essa variação exógena pode ser utilizada para estimar \(\beta_{i}\).
É por isso que no primeiro estágio da estimação, partimos da relação entre \(X\) e \(Z\):
\[X_{i} = \pi_{0} + \pi_{1}Z_{i} + v_{i}\]
A variável dependente (no modelo original) \(lnP\) é regredida contra a variável instrumental \(Z_{i}\) (= tdiff). O preço é desmembrado em três componentes, um da própria demanda (\(\pi_{0}\)), outro dos impostos cobrados (\(Z_{i}\)) e mais um com os fatores não observados (\(v_{i}\)). Assumimos, portanto, que os impostos não tem relação com a demanda, já que o governo é quem o faz arbitráriamente (exogeneidade!). Adicionalmente, é razoável assumir que os estados que cobram impostos mais altos sobre o preço do paiêro tem os maiores preços (\(X_{i}\) = \(lnP_{i}\)) registrados para esse bem (relevãncia!).
### IV regression
# Let's take a look at the first-stage regression (12.9)
# Regress X on Z (and any other exogenous variables W if we had them)
# As an instrument we use tdiff
first1 = lm(lprice ~ tdiff, data = cigs)
stargazer(first1,
se=list(cse(first1)),
title="First-stage Regression", type="text",
df=FALSE, digits=5)##
## First-stage Regression
## ===============================================
## Dependent variable:
## ---------------------------
## lprice
## -----------------------------------------------
## tdiff 0.03073***
## (0.00484)
##
## Constant 4.61655***
## (0.02892)
##
## -----------------------------------------------
## Observations 48
## R2 0.47100
## Adjusted R2 0.45950
## Residual Std. Error 0.09394
## F Statistic 40.95588***
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Desta forma, temos a primeira equação:
\[\widehat{ln(P^{_{i}^{cigaretets}})} = \underset{(0.030)}{4.62} - \underset{(0.005)}{0.031}SalesTax_{i}\]
Estocaremos os valores estimados \(\widehat{ln(P^{_{i}^{cigaretets}})}\) no data.frame:
cigs$p_hat <- first1$fitted.valuesE utilizamos ele mesmo para fazer o cálculo do segundo estágio:
# Reduced form: Regress Y on Z (and any other exogenous variables W if we had them)
reduced = lm(lquant ~ p_hat, data = cigs)
stargazer(reduced,
se = list(cse(reduced)),
title = "Reduced-form Regression", type = "text",
df = FALSE, digits = 5)##
## Reduced-form Regression
## ===============================================
## Dependent variable:
## ---------------------------
## lquant
## -----------------------------------------------
## p_hat -1.08359***
## (0.33370)
##
## Constant 9.71988***
## (1.59712)
##
## -----------------------------------------------
## Observations 48
## R2 0.15249
## Adjusted R2 0.13407
## Residual Std. Error 0.22645
## F Statistic 8.27665***
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01O que nos dá uma demanda bastante inelástica (-1,08): um aumento de 1% no preço da paranga reduz a o consumo em 1,08%
\[\widehat{ln(Q^{_{i}^{cigaretets}})} = \underset{(1.597)}{9.72} - \underset{(0.333)}{1.08}ln(P_{i}^{cigarette})\]
AER::ivreg()Função AER::ivreg(): O pacote AER tem uma infraestrutura consistente para estimação em 2 estágios. Como exemplo de implementação para variáveis instrumental, é trabalhado o modelo desenvolvido no capítulo 12 do livro Introduction to Econometrics (3th ed.) (Stock & Watson, 2010).
ivpack::robust.se: fornece infraestrutura para inclusão de erros padrões rocustos.
O mesmo resultado pode ser obtido com o comando ivreg:
iv1 = ivreg(lquant ~ lprice | tdiff , data = cigs)
library(ivpack)
# corrected SEs for IV regressions... slight difference from S&W method
ivse = function(reg) {
rob = robust.se(reg)[,2]
return(rob)
}
# note: use the function ivse for corrected SEs in IV
stargazer(iv1,
se=list(ivse(iv1)),
title="IV Regression", type="text",
df=FALSE, digits=5)## [1] "Robust Standard Errors"
##
## IV Regression
## ===============================================
## Dependent variable:
## ---------------------------
## lquant
## -----------------------------------------------
## lprice -1.08359***
## (0.31220)
##
## Constant 9.71988***
## (1.49614)
##
## -----------------------------------------------
## Observations 48
## R2 0.40113
## Adjusted R2 0.38811
## Residual Std. Error 0.19035
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01\[y = X\beta + Z\delta + \varepsilon\] O modelo precisa ser estimado por 2SLS. Escreva a função que estime o modelo com quatro argumentos: \(y\),\(X\),\(Z\); e a matriz de instrumentos \(H\):
Resposta:
Criamos a função abaixo para estimação em 2 estágios:
TSLS <- function(y, X, Z, H){
y <- as.matrix(y) # Variável dependente
X <- as.matrix(X) # Variável explicativa exógena
Z <- as.matrix(Z) # Variável explicativa endógena
H <- as.matrix(H) # Matriz de Instrumentos (inclui X e os demais)
EXP <- cbind(X, Z)
betas <- solve(t(EXP) %*% H %*% solve(t(H) %*% H) %*% t(H) %*% EXP) %*%
t(EXP) %*% H %*% solve(t(H) %*% H) %*% t(H) %*% y
output <- data.frame(betas)
colnames(output) <- "Coeficientes"
rownames(output) <- NULL
print(output)
}Um teste:
# Teste: # Dados
set.seed(44)
y <- rnorm(500)
Xa <- rnorm(500)
Int <- rep(1,500)
X <- cbind(Int, Xa)
Z <- rnorm(500)
X1 <- rnorm(500) # Instrumento 1
X2 <- rnorm(500) # Instrumento 2
H <- cbind(Xa, X1, X2)
head(H) ## Xa X1 X2
## [1,] 1.7495518 0.74350693 -0.7160694
## [2,] -0.4433087 -0.55804034 0.3438573
## [3,] -0.3431100 -0.69004880 1.2326432
## [4,] 1.2491928 -1.17804117 -0.5013019
## [5,] -0.2989443 -2.28623172 -2.6151014
## [6,] -0.6738205 -0.07962112 0.6173945TSLS(y,X,Z,H)## Coeficientes
## 1 0.16282615
## 2 0.02228528
## 3 -0.34957993