Source code for linearmodels.system.gmm

"""
Covariance and weight estimation for GMM IV estimators
"""
from __future__ import annotations

from typing import Optional, Sequence, cast

from numpy import array, empty, ndarray, repeat, sqrt, zeros_like

from linearmodels.asset_pricing.covariance import _HACMixin
from linearmodels.iv.covariance import kernel_optimal_bandwidth
from linearmodels.shared.utility import AttrDict
from linearmodels.system._utility import blocked_inner_prod
from linearmodels.typing import Float64Array


[docs]class HomoskedasticWeightMatrix(object): r""" Homoskedastic (unadjusted) weight estimation Parameters ---------- center : bool Flag indicating whether to center the moment conditions by subtracting the mean before computing the weight matrix. debiased : bool Flag indicating whether to use small-sample adjustments 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: Optional[float] = 0 self._name = "Homoskedastic (Unadjusted) Weighting" self._config = AttrDict(center=center, debiased=debiased) def __str__(self) -> str: out = self._name extra = [] for key in self._str_extra: extra.append(": ".join([str(key), str(self._str_extra[key])])) if extra: out += " (" + ", ".join(extra) + ")" return out def __repr__(self) -> str: return self.__str__() + ", id: {0}".format(hex(id(self))) @property def _str_extra(self) -> AttrDict: return AttrDict(Debiased=self._debiased, Center=self._center)
[docs] def sigma(self, eps: Float64Array, x: Sequence[Float64Array]) -> Float64Array: """ Estimate residual covariance. Parameters ---------- eps : ndarray The residuals from the system of equations. x : list[ndarray] A list of the regressor matrices for each equation in the system. Returns ------- ndarray The estimated covariance matrix of the residuals. """ nobs = eps.shape[0] eps = eps - eps.mean(0) sigma = eps.T @ eps / nobs scale = 1.0 if self._debiased: k = 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[Float64Array], z: Sequence[Float64Array], eps: Float64Array, *, sigma: ndarray, ) -> Float64Array: """ Construct a GMM weight matrix for a model. Parameters ---------- x : list[ndarray] List of containing model regressors for each equation in the system z : list[ndarray] List of containing instruments for each equation in the system eps : ndarray Model errors (nobs by neqn) sigma : ndarray Fixed covariance of model errors. If None, estimated from eps. Returns ------- ndarray Covariance of GMM moment conditions. """ 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 Flag indicating whether to center the moment conditions by subtracting the mean before computing the weight matrix. debiased : bool Flag indicating whether to use small-sample adjustments 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(HeteroskedasticWeightMatrix, self).__init__(center, debiased) self._name = "Heteroskedastic (Robust) Weighting"
[docs] def weight_matrix( self, x: Sequence[Float64Array], z: Sequence[Float64Array], eps: Float64Array, *, sigma: Optional[ndarray] = None, ) -> Float64Array: """ Construct a GMM weight matrix for a model. Parameters ---------- x : list[ndarray] Model regressors (exog and endog), (nobs by nvar) z : list[ndarray] Model instruments (exog and instruments), (nobs by ninstr) eps : ndarray Model errors (nobs by 1) sigma : ndarray Fixed covariance of model errors. If None, estimated from eps. Returns ------- ndarray Covariance of GMM moment conditions. """ nobs = x[0].shape[0] k = len(x) k_total = sum(map(lambda a: a.shape[1], 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[Float64Array], z: Sequence[Float64Array] ) -> Float64Array: nvar = 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(Float64Array, sqrt(nvar))[:, None] scale = nobs / (nobs - nvar @ nvar.T) return scale
[docs]class KernelWeightMatrix(HeteroskedasticWeightMatrix, _HACMixin): r""" Heteroskedasticity robust weight estimation Parameters ---------- center : bool Flag indicating whether to center the moment conditions by subtracting the mean before computing the weight matrix. debiased : bool Flag indicating whether to use small-sample adjustments kernel : str Name of kernel to use. Supported kernels include: * "bartlett", "newey-west" : Bartlett's kernel * "parzen", "gallant" : Parzen's kernel * "qs", "quadratic-spectral", "andrews" : Quadratic spectral kernel bandwidth : float Bandwidth to use for the kernel. If not provided the optimal bandwidth will be estimated. optimal_bw : bool Flag indicating whether to estimate the optimal bandwidth, when bandwidth is None. If False, nobs - 2 is used 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. """ def __init__( self, center: bool = False, debiased: bool = False, kernel: str = "bartlett", bandwidth: Optional[float] = None, optimal_bw: bool = False, ) -> None: _HACMixin.__init__(self, kernel, bandwidth) super(KernelWeightMatrix, self).__init__(center, debiased) self._name = "Kernel (HAC) Weighting" self._check_kernel(kernel) self._check_bandwidth(bandwidth) self._predefined_bw = self._bandwidth self._optimal_bw = optimal_bw
[docs] def weight_matrix( self, x: Sequence[Float64Array], z: Sequence[Float64Array], eps: Float64Array, *, sigma: Optional[ndarray] = None, ) -> Float64Array: """ Construct a GMM weight matrix for a model. Parameters ---------- x : List[ndarray] Model regressors (exog and endog) z : List[ndarray] Model instruments (exog and instruments) eps : ndarray Model errors (nobs by nequation) sigma : ndarray Fixed covariance of model errors. If None, estimated from eps. Returns ------- ndarray Covariance of GMM moment conditions. """ nobs = x[0].shape[0] k = len(x) k_total = sum(map(lambda a: a.shape[1], 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: 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: 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 """ out = AttrDict([(k, v) for k, v in self._config.items()]) out["bandwidth"] = self.bandwidth return out