"""
Covariance and weight estimation for GMM IV estimators.
"""
from __future__ import annotations
from collections.abc import Sequence
from typing import cast
from numpy import array, empty, errstate, ndarray, repeat, sqrt, zeros_like
from linearmodels.asset_pricing.covariance import _HACMixin
from linearmodels.iv.covariance import kernel_optimal_bandwidth
from linearmodels.iv.gmm import _warn_default_bandwidth
from linearmodels.shared.utility import AttrDict
from linearmodels.system._utility import blocked_inner_prod
import linearmodels.typing.data
[docs]
class HomoskedasticWeightMatrix:
r"""
Homoskedastic (unadjusted) weight estimation.
Parameters
----------
center : bool, optional
Flag indicating whether to center the moment conditions by subtracting
the mean before computing the weight matrix. The default is False.
debiased : bool, optional
Flag indicating whether to use small-sample adjustments. The default is
False.
Notes
-----
The weight matrix estimator is
.. math::
Z'(\Sigma \otimes I_N)Z
where :math:`Z` is a block diagonal matrix containing both the exogenous
regressors and instruments and :math:`\Sigma` is the covariance of the
model residuals.
``center`` has no effect on this estimator since it is always centered.
"""
def __init__(self, center: bool = False, debiased: bool = False) -> None:
self._center = center
self._debiased = debiased
self._bandwidth: float | None = 0
self._name = "Homoskedastic (Unadjusted) Weighting"
self._config = AttrDict(center=center, debiased=debiased)
def __str__(self) -> str:
out = self._name
extra = [
": ".join([str(key), str(self._str_extra[key])]) for key in self._str_extra
]
if extra:
out += " (" + ", ".join(extra) + ")"
return out
def __repr__(self) -> str:
return self.__str__() + f", id: {hex(id(self))}"
@property
def _str_extra(self) -> AttrDict:
return AttrDict(Debiased=self._debiased, Center=self._center)
[docs]
def sigma(
self,
eps: linearmodels.typing.data.Float64Array,
x: Sequence[linearmodels.typing.data.Float64Array],
) -> linearmodels.typing.data.Float64Array:
"""
Estimate residual covariance.
Parameters
----------
eps : ndarray
The residuals from the system of equations, with shape (nobs, neq).
x : sequence of ndarray
The regressors of each equation in the system, with shape
(nobs, k_i) for equation i.
Returns
-------
ndarray
The estimated covariance matrix of the residuals.
"""
nobs = eps.shape[0]
eps = eps - eps.mean(0)
sigma = eps.T @ eps / nobs
scale: float | ndarray = 1.0
if self._debiased:
k: ndarray = array([a.shape[1] for a in x])[:, None]
k = sqrt(k)
scale = nobs / (nobs - k @ k.T)
sigma *= scale
return sigma
[docs]
def weight_matrix(
self,
x: Sequence[linearmodels.typing.data.Float64Array],
z: Sequence[linearmodels.typing.data.Float64Array],
eps: linearmodels.typing.data.Float64Array,
*,
sigma: ndarray,
) -> linearmodels.typing.data.Float64Array:
"""
Construct a GMM weight matrix for a model.
Parameters
----------
x : sequence of ndarray
The regressors of each equation in the system, with shape
(nobs, k_i) for equation i.
z : sequence of ndarray
The instruments of each equation in the system, with shape
(nobs, nz_i) for equation i.
eps : ndarray
Model errors, with shape (nobs, neqn). They are not used by this
estimator, which uses ``sigma``.
sigma : ndarray
Covariance of the model errors, with shape (neqn, neqn).
Returns
-------
ndarray
Weight matrix, which is the covariance of the GMM moment conditions
with shape (nmoments, nmoments).
"""
nobs = z[0].shape[0]
w = cast("ndarray", blocked_inner_prod(z, sigma) / nobs)
return w
@property
def config(self) -> AttrDict:
"""
Weight estimator configuration.
Returns
-------
AttrDict
Dictionary containing weight estimator configuration information.
"""
return self._config
[docs]
class HeteroskedasticWeightMatrix(HomoskedasticWeightMatrix):
r"""
Heteroskedasticity robust weight estimation.
Parameters
----------
center : bool, optional
Flag indicating whether to center the moment conditions by subtracting
the mean before computing the weight matrix. The default is False.
debiased : bool, optional
Flag indicating whether to use small-sample adjustments. The default is
False.
Notes
-----
The weight matrix estimator is
.. math::
W & = n^{-1}\sum_{i=1}^{n}g'_ig_i \\
g_i & = (z_{1i}\epsilon_{1i},z_{2i}\epsilon_{2i},\ldots,z_{ki}\epsilon_{ki})
where :math:`g_i` is the vector of scores across all equations for
observation i. :math:`z_{ji}` is the vector of instruments for equation
j and :math:`\epsilon_{ji}` is the error for equation j for observation
i. This form allows for heteroskedasticity and arbitrary cross-sectional
dependence between the moment conditions.
"""
def __init__(self, center: bool = False, debiased: bool = False) -> None:
super().__init__(center, debiased)
self._name = "Heteroskedastic (Robust) Weighting"
[docs]
def weight_matrix(
self,
x: Sequence[linearmodels.typing.data.Float64Array],
z: Sequence[linearmodels.typing.data.Float64Array],
eps: linearmodels.typing.data.Float64Array,
*,
sigma: ndarray | None = None,
) -> linearmodels.typing.data.Float64Array:
"""
Construct a GMM weight matrix for a model.
Parameters
----------
x : sequence of ndarray
The regressors (exog and endog) of each equation in the system, with
shape (nobs, k_i) for equation i.
z : sequence of ndarray
The instruments (exog and instruments) of each equation in the
system, with shape (nobs, nz_i) for equation i.
eps : ndarray
Model errors, with shape (nobs, neqn).
sigma : ndarray, optional
Not used by this estimator. It is accepted so that all the weight
matrix estimators have the same signature.
Returns
-------
ndarray
Weight matrix, which is the covariance of the GMM moment conditions
with shape (nmoments, nmoments).
"""
nobs = x[0].shape[0]
k = len(x)
k_total = sum(a.shape[1] for a in z)
ze = empty((nobs, k_total))
loc = 0
for i in range(k):
e = eps[:, [i]]
zk = z[i].shape[1]
ze[:, loc : loc + zk] = z[i] * e
loc += zk
mu = ze.mean(axis=0) if self._center else 0
ze -= mu
w = ze.T @ ze / nobs
scale = self._debias_scale(nobs, x, z)
w *= scale
return w
def _debias_scale(
self,
nobs: int,
x: Sequence[linearmodels.typing.data.Float64Array],
z: Sequence[linearmodels.typing.data.Float64Array],
) -> linearmodels.typing.data.Float64Array:
nvar: ndarray = array([a.shape[1] for a in x])
ninstr = array([a.shape[1] for a in z])
nvar = repeat(nvar, ninstr)
if not self._debiased:
nvar = zeros_like(nvar)
nvar = cast("linearmodels.typing.data.Float64Array", sqrt(nvar))[:, None]
scale = nobs / (nobs - nvar @ nvar.T)
return scale
[docs]
class KernelWeightMatrix(HeteroskedasticWeightMatrix, _HACMixin):
r"""
Heteroskedasticity and autocorrelation robust (HAC) weight estimation.
Parameters
----------
center : bool, optional
Flag indicating whether to center the moment conditions by subtracting
the mean before computing the weight matrix. The default is False.
debiased : bool, optional
Flag indicating whether to use small-sample adjustments. The default is
False.
kernel : str, optional
Name of the kernel. The default is "bartlett". The names are not case
sensitive. Supported kernels are:
* "bartlett", "newey-west" - Bartlett's kernel
* "parzen", "gallant" - Parzen's kernel
* "qs", "quadratic-spectral", "andrews" - Quadratic spectral kernel
bandwidth : float, optional
Non-negative bandwidth to use for the kernel. If not provided, the
bandwidth is selected by ``optimal_bw``.
optimal_bw : bool, optional
Flag indicating whether to estimate the bandwidth from the data when
``bandwidth`` is not provided. If False, nobs - 2 is used. If True, the
bandwidth is estimated using
:func:`~linearmodels.iv.covariance.kernel_optimal_bandwidth`, which
needs at least 4 observations with the Bartlett kernel, 5 with the
Parzen kernel and 6 with the quadratic spectral kernel. The default
also estimates the bandwidth, but it issues a
``DefaultBandwidthWarning`` since it changed in 7.1.
Notes
-----
The weight matrix estimator is
.. math::
W & = \hat{\Gamma}_0+\sum_{i=1}^{n-1} w_i (\hat{\Gamma}_i+\hat{\Gamma}_i^\prime) \\
\hat{\Gamma}_j & = n^{-1}\sum_{i=1}^{n-j} g'_ig_{i+j} \\
g_i & = (z_{1i}\epsilon_{1i},z_{2i}\epsilon_{2i},\ldots,z_{ki}\epsilon_{ki})
where :math:`g_i` is the vector of scores across all equations for
observation i and :math:`w_j` are the kernel weights which depend on the
selected kernel and bandwidth. :math:`z_{ji}` is the vector of instruments
for equation j and :math:`\epsilon_{ji}` is the error for equation j for
observation i. This form allows for heteroskedasticity and autocorrelation
between the moment conditions.
If neither ``bandwidth`` nor ``optimal_bw`` is given, the bandwidth is
estimated from the data every time that the weight matrix is computed.
In 7.0 and earlier the default was ``nobs - 2``, which is not a consistent
estimator of the covariance of the moment conditions. Using the default
issues a ``DefaultBandwidthWarning`` (a ``UserWarning``), once for each
instance. The previous behavior is restored by ``optimal_bw=False`` or by
``bandwidth=nobs - 2``. Passing ``bandwidth`` or ``optimal_bw`` silences the
warning, and ``optimal_bw=True`` selects the new default without it.
.. versionchanged:: 7.1
``weight_matrix`` raises a ``ValueError`` if the bandwidth is estimated
from the data and there are fewer observations than the kernel needs
(see ``optimal_bw``), or if the estimated bandwidth is not finite. Pass
``bandwidth`` or ``optimal_bw=False`` to use the weight matrix with
fewer observations. Earlier versions raised "cannot convert float
infinity to integer" (an ``OverflowError``) or "cannot convert float NaN
to integer" for some data when ``optimal_bw`` was True, and used a
bandwidth of nobs - 1, or another bandwidth, for other data.
"""
def __init__(
self,
center: bool = False,
debiased: bool = False,
kernel: str = "bartlett",
bandwidth: float | None = None,
optimal_bw: bool | None = None,
) -> None:
_HACMixin.__init__(self, kernel, bandwidth)
super().__init__(center, debiased)
self._name = "Kernel (HAC) Weighting"
self._check_kernel(kernel)
self._check_bandwidth(bandwidth)
self._predefined_bw = self._bandwidth
# Selecting the bandwidth from the data is the default since 7.1, and it
# warns since the default was nobs - 2
self._default_bandwidth = bandwidth is None and optimal_bw is None
self._optimal_bw = True if optimal_bw is None else optimal_bw
self._warned = False
[docs]
def weight_matrix(
self,
x: Sequence[linearmodels.typing.data.Float64Array],
z: Sequence[linearmodels.typing.data.Float64Array],
eps: linearmodels.typing.data.Float64Array,
*,
sigma: ndarray | None = None,
) -> linearmodels.typing.data.Float64Array:
"""
Construct a GMM weight matrix for a model.
Parameters
----------
x : sequence of ndarray
The regressors (exog and endog) of each equation in the system, with
shape (nobs, k_i) for equation i.
z : sequence of ndarray
The instruments (exog and instruments) of each equation in the
system, with shape (nobs, nz_i) for equation i.
eps : ndarray
Model errors, with shape (nobs, neqn).
sigma : ndarray, optional
Not used by this estimator. It is accepted so that all the weight
matrix estimators have the same signature.
Returns
-------
ndarray
Weight matrix, which is the covariance of the GMM moment conditions
with shape (nmoments, nmoments).
"""
nobs = x[0].shape[0]
k = len(x)
k_total = sum(a.shape[1] for a in z)
ze = empty((nobs, k_total))
loc = 0
for i in range(k):
e = eps[:, [i]]
zk = z[i].shape[1]
ze[:, loc : loc + zk] = z[i] * e
loc += zk
mu = ze.mean(axis=0) if self._center else 0
ze -= mu
self._optimal_bandwidth(ze)
w = self._kernel_cov(ze)
scale = self._debias_scale(nobs, x, z)
w *= scale
return w
def _optimal_bandwidth(
self, moments: linearmodels.typing.data.Float64Array
) -> float:
"""Compute optimal bandwidth used in estimation if needed."""
if self._predefined_bw is not None:
return self._predefined_bw
elif not self._optimal_bw:
self._bandwidth = moments.shape[0] - 2
else:
if self._default_bandwidth and not self._warned:
# Once for each instance, and not for each call
self._warned = True
_warn_default_bandwidth(
weights=True, covariance=False, stacklevel=4, model="IVSystemGMM"
)
# A moment that does not vary has a standard deviation of 0, and the
# series is NaN. kernel_optimal_bandwidth raises, so do not warn
with errstate(divide="ignore", invalid="ignore"):
m = moments / moments.std(0)[None, :]
m = m.sum(1)
self._bandwidth = kernel_optimal_bandwidth(m, kernel=self.kernel)
assert self._bandwidth is not None
return self._bandwidth
@property
def bandwidth(self) -> float:
"""Bandwidth used to estimate covariance of moment conditions."""
assert self._bandwidth is not None
return self._bandwidth
@property
def config(self) -> AttrDict:
"""
Weight estimator configuration.
Returns
-------
AttrDict
Dictionary containing weight estimator configuration information.
The bandwidth is None if the bandwidth was not given and the weight
matrix has not been computed, as in ``IVSystemGMM.fit`` with
``iter_limit=1``.
"""
out = AttrDict([(k, v) for k, v in self._config.items()])
out["bandwidth"] = self._bandwidth
return out