bccorr | R Documentation |
With a model like y = X\beta + D\theta + F\psi + \epsilon
, where D
and F
are matrices with dummy encoded factors, one application of lfe is to
study the correlation cor(D\theta, F\psi)
. However, if we use
estimates for \theta
and \psi
, the resulting correlation is biased.
The function bccorr
computes a bias corrected correlation
as described in Gaure (2014).
bccorr(
est,
alpha = getfe(est),
corrfactors = 1L:2L,
nocovar = (length(est$X) == 0) && length(est$fe) == 2,
tol = 0.01,
maxsamples = Inf,
lhs = NULL
)
est |
an object of class '"felm"', the result of a call to
|
alpha |
a data frame, the result of a call to |
corrfactors |
integer or character vector of length 2. The factors to correlate. The default is fine if there are only two factors in the model. |
nocovar |
logical. Assume no other covariates than the two factors are present, or that they are uncorrelated with them. |
tol |
The absolute tolerance for the bias-corrected correlation. |
maxsamples |
Maximum number of samples for the trace sample means estimates |
lhs |
character. Name of left hand side if multiple left hand sides. |
The bias expressions from Andrews et al. are of the form tr(AB^{-1}C)
where A
, B
, and C
are matrices too large to be handled
directly. bccorr
estimates the trace by using the formula tr(M) = E(x^t M x)
where x is a vector with coordinates drawn uniformly from the set \{-1,1\}
.
More specifically, the expectation is estimated by
sample means, i.e. in each sample a vector x is drawn, the
equation Bv = Cx
is solved by a conjugate gradient method, and the
real number x^t Av
is computed.
There are three bias corrections, for the variances of D\theta
(vD
) and
F\psi
(vF
), and their covariance (vDF
).The correlation is computed as
rho <- vDF/sqrt(vD*vF)
. The variances are estimated to a
relative tolerance specified by the argument tol
. The covariance
bias is estimated to an absolute tolerance in the correlation rho
(conditional on the already bias corrected vD
and vF
) specified by
tol
. The CG algorithm does not need to be exceedingly precise,
it is terminated when the solution reaches a precision which is
sufficient for the chosen precision in vD, vF, vDF
.
If est
is the result of a weighted felm()
estimation,
the variances and correlations are weighted too.
bccorr
returns a named integer vector with the following fields:
corr |
the bias corrected correlation. |
v1 |
the bias corrected variance for the first factor specified
by |
v2 |
the bias corrected variance for the second factor. |
cov |
the bias corrected covariance between the two factors. |
d1 |
the bias correction for the first factor. |
d2 |
the bias correction for the second factor. |
d12 |
the bias correction for covariance. |
The bias corrections have been subtracted from the bias estimates. E.g. v2 = v2' - d2, where v2' is the biased variance.
Bias correction for IV-estimates are not supported as of now.
Note that if est
is the result of a call to felm()
with keepX=FALSE
(the default), the correlation will be computed
as if the covariates X are independent of the two factors. This will be
faster (typically by a factor of approx. 4), and possibly wronger.
Note also that the computations performed by this function are non-trivial, they may take quite some time. It would be wise to start out with quite liberal tolerances, e.g. tol=0.1, to get an idea of the time requirements.
The algorithm used is not very well suited for small datasets with only a few thousand levels in the factors.
Gaure, S. (2014), Correlation bias correction in two-way fixed-effects linear regression, Stat 3(1):379:390, 2014.
fevcov()
x <- rnorm(500)
x2 <- rnorm(length(x))
## create individual and firm
id <- factor(sample(40, length(x), replace = TRUE))
firm <- factor(sample(30, length(x), replace = TRUE, prob = c(2, rep(1, 29))))
foo <- factor(sample(20, length(x), replace = TRUE))
## effects
id.eff <- rnorm(nlevels(id))
firm.eff <- rnorm(nlevels(firm))
foo.eff <- rnorm(nlevels(foo))
## left hand side
y <- x + 0.25 * x2 + id.eff[id] + firm.eff[firm] + foo.eff[foo] + rnorm(length(x))
# make a data frame
fr <- data.frame(y, x, x2, id, firm, foo)
## estimate and print result
est <- felm(y ~ x + x2 | id + firm + foo, data = fr, keepX = TRUE)
# find bias corrections
bccorr(est)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.