library("knitr") knitr::opts_knit$set(self.contained = FALSE) knitr::opts_chunk$set(tidy = TRUE, collapse=TRUE, comment = "#>", tidy.opts=list(blank=FALSE, width.cutoff=55))
The package GMMSensitivity
implements estimators and confidence interval for
sensitivity analysis in moment condition models considered in @ak18sensitivity.
In this vignette, we demonstrate the implementation of these confidence
intervals using the demand model for automobiles from @blp95.
The package includes the dataset blp
, which contains estimates of the @blp95
model, as implemented by @ags17. A description of the dataset can be obtained
via the documentation in R, using ?blp
. We give additional description here:
G
: Matrix with 31 rows and 17 columns, estimate of derivative of the moment
condition $G=\partial E[g(w_i,\theta)]/\partial \theta$, evaluated at initial
estimate $\hat{\theta}_{\text{initial}}$. The initial estimate corresponds to
a second-step GMM estimate, as computed by @ags17.
H
: Vector of length 17, estimate of derivative of average markup $h(\theta)$,
evaluated at $\hat{\theta}_{\text{initial}}$.
W
: Weight matrix used to obtain the estimate $\hat{\theta}_{\text{initial}}$. It
corresponds to an estimate of the inverse of the variance of the moment
condition $\Sigma$.
g_init
: Average moment condition $n^{-1}\sum_{i=1}^{n}g(w_i,\theta)$, evaluated at
estimate $\hat{\theta}_{\text{initial}}$ of $\theta$ from @blp95.
h_init
: Estimate of the average markup, $h(\hat{\theta}_{\text{initial}})$.
names
: Two lists, one for names of the moment conditions, and one for elements of theta
ZZ
: Gram matrix of the instruments, used to specify the set $\mathcal{C}$
Sig
: Variance of moment condition, estimate of $\Sigma$, given by the sample
variance of the moment condition evaluated at $\hat{\theta}_{\text{initial}}$.
sdZ
: vector of standard deviations of the instruments
perturb
: scaling parameter to give interpretable meaning to violations of moment
conditions. For demand-side moments, it corresponds to an estimate of
$\delta_{d}=0.01\overline{p}/(\overline{y}/\alpha)$, and for supply-side
moments, it corresponds to an estimate of
$\delta_{s}=-0.01\overline{p}/\overline{m c}$, where $\overline{p}$ is average
price, $\overline{y}$ is average income, $\overline{m c}$ is average marginal
cost, and $\alpha$ is a parameter in the utility function. With this scaling,
if a given demand-side instrument enters the utility function with coefficient
$\delta_{d}$, a consumer is willing to pay 1\% of the average car price for a
unit increase in the instrument. If a given supply-side instrument enters the
cost function with coefficient $\delta_s$, increasing the instrument by one
unit decreases the marginal cost by 1\% of the average car price.
n
: Sample size, number of model/years.
library("GMMSensitivity")
The package implements estimators, confidence intervals, and efficiency calculations for the model (in the notation of @ak18sensitivity) \begin{equation} g(\theta_{0})=c/\sqrt{n},\quad c\in \mathcal{C},\quad \mathcal{C}={B\gamma\colon \norm{\gamma}_{p}\leq M}. \end{equation}
Suppose that we want to allow all excluded instruments in the @blp95 application to be potentially invalid. Fix $B$ to a scaling matrix so that if the $j$th supply-side instrument is invalid with $\gamma_{s j}=1$, this means that changing the instrument by one standard deviation changes the marginal cost by $\gamma_{s j}$\% of the average car price, and if the $j$ the demand-side instrument is invalid with $\gamma_{d j}=1$, then the consumer willingness to pay for one standard deviation change in the instrument is $\gamma_{d j}\%$ of the average 1980 car price. Let $p=2$, and $M=\sqrt{# I}$, where $#I$ is the number of invalid instruments so that $\gamma=1$ is included in the set (this is the same scaling as described in @ak18sensitivity in the paper). Then the confidence interval can be constructed as follows:
## Construct estimate of initial sensitivity blp$k_init <- -drop(blp$H %*% solve(crossprod(blp$G, blp$W %*% blp$G), crossprod(blp$G, blp$W))) ## list collecting initial estimates of H, G, Sigma, n, g(thetahat), initial ## sensitivity k, and initial estimate of average markup h(thetahat) eo <- list(H=blp$H, G=blp$G, Sig=blp$Sig, n=blp$n, g_init=blp$g_init, k_init=blp$k_init, h_init= blp$h_init) ## Rows corresponding to invalid instruments I <- vector(mode="logical", length=nrow(eo$G)) I[c(6:13, 20:31)] <- TRUE ## Matrix B, scaled as described in the paper B0 <- blp$ZZ %*% diag(sqrt(blp$n)*abs(blp$perturb)/blp$sdZ) ## Value of M M0 <- sqrt(sum(I)) ## Select columns of B0 corresponding to invalid instruments OptEstimator(eo, B0[, I], M=M0, p=2, alpha=0.05, opt.criterion="FLCI")
The efficiency $\kappa_{*}$ for this confidence interval can be computed using the
EffBounds
function (which can also be used to compute efficiency of one-sided
confidence intervals):
EffBounds(eo, B0[, I], M=M0, p=2)$twosided
In contrast the CI based on the initial estimate is much wider:
OptEstimator(eo, B0[, I], M=M0, p=2, alpha=0.05, opt.criterion="Valid")
A specification test for whether the value $M=M0$ is too low, that is a test of
the hypothesis $H_0\colon M\leq M0$, can be conducted using the Jtest
function:
## Update eo so that Sig corresponds to the initial estimate of ## Sigma, so that thetahat_initial minimizes n*g(theta)Sig^{-1} g(theta), ## and the J statistic is given by the value of this minimum. eoJ <- eo eoJ$Sig <- solve(blp$W) jt <- Jtest(eoJ, B0[, I], M=M0, p=2, alpha=0.05)
Here J
is the $J$-statistic, p0
is the p-value of the usual $J$-test (that
assumes $M=0$ under the null), pC
is the $p$-value of the test, and Mmin
is
the smallest value of $M$ that is not rejected. Rescaling by $\sqrt{# I}$, we
obtain the result in Table 1 in @ak18sensitivity:
jt$Mmin/M0
If were only concerned about the validity of the demand-side instrument
"demand_firm_const"
(number of cars produced by the same firm), then since
this is the sixth instrument, an analogous analysis could be conducted as:
I <- vector(mode="logical", length=nrow(eo$G)) I[6] <- TRUE OptEstimator(eo, B0[, I, drop=FALSE], M=1, p=2, alpha=0.05, opt.criterion="FLCI") EffBounds(eo, B0[, I, drop=FALSE], M=1, p=2)$twosided OptEstimator(eo, B0[, I, drop=FALSE], M=1, p=2, alpha=0.05, opt.criterion="Valid") Jtest(eoJ, B0[, I, drop=FALSE], M=1, p=2, alpha=0.05)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.