Skip to content

Commit

Permalink
REF: Refactor common code out of asset pricing
Browse files Browse the repository at this point in the history 
Refactor covariance for HAC estimation to share code
Refactor formula code to share code

xref #89
  • Loading branch information
@bashtage
bashtage committed Nov 30, 2017
1 parent b9bc20f commit 25aec32
Show file tree
Hide file tree
Showing 2 changed files with 73 additions and 87 deletions.
117 changes: 53 additions & 64 deletions linearmodels/asset_pricing/covariance.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,49 @@
kernel_optimal_bandwidth)


class _HACMixin(object):

def __init__(self):
self._bandwidth = None

@property
def kernel(self):
"""Kernel used in estimation"""
return self._kernel

@property
def bandwidth(self):
"""Bandwidth used in estimation"""
if self._bandwidth is None:
moments = self._moments
m = moments / moments.std(0)[None, :]
m = m.sum(1)
bw = kernel_optimal_bandwidth(m, kernel=self.kernel)
self._bandwidth = int(bw)

return self._bandwidth

def _check_kernel(self, kernel):
self._kernel = kernel.lower()
if self._kernel not in KERNEL_LOOKUP:
raise ValueError('Unknown kernel')

def _check_bandwidth(self, bandwidth):
self._bandwidth = bandwidth
if bandwidth is not None:
if bandwidth < 0:
raise ValueError('bandwidth must be non-negative.')

def _kernel_cov(self, z):
nobs = z.shape[0]
bw = self.bandwidth
kernel = self._kernel
kernel = KERNEL_LOOKUP[kernel]
weights = kernel(bw, nobs - 1)
out = _cov_kernel(z, weights)
return (out + out.T) / 2


class HeteroskedasticCovariance(object):
"""
Heteroskedasticity robust covariance estimator
Expand All @@ -32,7 +75,7 @@ class HeteroskedasticCovariance(object):
def __init__(self, xe, *, jacobian=None, inv_jacobian=None,
center=True, debiased=False, df=0):

self._xe = xe
self._moments = self._xe = xe
self._jac = jacobian
self._inv_jac = inv_jacobian
self._center = center
Expand Down Expand Up @@ -116,7 +159,7 @@ def cov(self):
return out


class KernelCovariance(HeteroskedasticCovariance):
class KernelCovariance(HeteroskedasticCovariance, _HACMixin):
"""
Heteroskedasticity-autocorrelation (HAC) robust covariance estimator
Expand Down Expand Up @@ -149,13 +192,8 @@ def __init__(self, xe, *, jacobian=None, inv_jacobian=None,
inv_jacobian=inv_jacobian,
center=center,
debiased=debiased, df=df)
self._kernel = kernel.lower()
if self._kernel not in KERNEL_LOOKUP:
raise ValueError('Unknown kernel')
self._bandwidth = bandwidth
if bandwidth is not None:
if bandwidth < 0:
raise ValueError('bandwidth must be non-negative.')
self._check_kernel(kernel)
self._check_bandwidth(bandwidth)

def __str__(self):
descr = ', Kernel: {0}, Bandwidth: {1}'.format(self._kernel,
Expand All @@ -169,23 +207,6 @@ def config(self):
out['bandwidth'] = self.bandwidth
return out

@property
def kernel(self):
"""Kernel used in estimation"""
return self._kernel

@property
def bandwidth(self):
"""Bandwidth used in estimation"""
if self._bandwidth is None:
xe = self._xe
x = xe / xe.std(0)[None, :]
x = x.sum(1)
bw = kernel_optimal_bandwidth(x, kernel=self.kernel)
self._bandwidth = int(bw)

return self._bandwidth

@property
def s(self):
"""
Expand All @@ -197,12 +218,7 @@ def s(self):
Covariance of the scores or moment conditions
"""
xe = self._xe
nobs = xe.shape[0]
bw = self.bandwidth
kernel = self._kernel
kernel = KERNEL_LOOKUP[kernel]
weights = kernel(bw, nobs - 1)
out = _cov_kernel(xe, weights)
out = self._kernel_cov(xe)

return (out + out.T) / 2

Expand Down Expand Up @@ -245,7 +261,7 @@ def w(self, moments):
return inv((out + out.T) / 2.0)


class KernelWeight(HeteroskedasticWeight):
class KernelWeight(HeteroskedasticWeight, _HACMixin):
"""
HAC GMM weighing matrix estimation
Expand All @@ -264,30 +280,8 @@ class KernelWeight(HeteroskedasticWeight):

def __init__(self, moments, center=True, kernel='bartlett', bandwidth=None):
super(KernelWeight, self).__init__(moments, center=center)
self._kernel = kernel.lower()
if self._kernel not in KERNEL_LOOKUP:
raise ValueError('Unknown kernel')
self._bandwidth = bandwidth
if bandwidth is not None:
if bandwidth < 0:
raise ValueError('bandwidth must be non-negative.')

@property
def kernel(self):
"""Kernel used in estimation"""
return self._kernel

@property
def bandwidth(self):
"""Bandwidth used in estimation"""
if self._bandwidth is None:
moments = self._moments
m = moments / moments.std(0)[None, :]
m = m.sum(1)
bw = kernel_optimal_bandwidth(m, kernel=self.kernel)
self._bandwidth = int(bw)

return self._bandwidth
self._check_kernel(kernel)
self._check_bandwidth(bandwidth)

def w(self, moments):
"""
Expand All @@ -305,11 +299,6 @@ def w(self, moments):
"""
if self._center:
moments = moments - moments.mean(0)[None, :]
nobs = moments.shape[0]
bw = self.bandwidth
kernel = self._kernel
kernel = KERNEL_LOOKUP[kernel]
weights = kernel(bw, nobs - 1)
out = _cov_kernel(moments, weights)
out = self._kernel_cov(moments)

return inv((out + out.T) / 2.0)
return inv(out)
43 changes: 20 additions & 23 deletions linearmodels/asset_pricing/model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,21 @@ def formula(self):
def formula(self, value):
self._formula = value

@staticmethod
def _prepare_data_from_formula(formula, data, portfolios):
na_action = NAAction(on_NA='raise', NA_types=[])
orig_formula = formula
if portfolios is not None:
factors = dmatrix(formula + ' + 0', data, return_type='dataframe', NA_action=na_action)
else:
formula = formula.split('~')
portfolios = dmatrix(formula[0].strip() + ' + 0', data,
return_type='dataframe', NA_action=na_action)
factors = dmatrix(formula[1].strip() + ' + 0', data,
return_type='dataframe', NA_action=na_action)

return factors, portfolios, orig_formula

@classmethod
def from_formula(cls, formula, data, *, portfolios=None):
"""
Expand Down Expand Up @@ -153,18 +168,9 @@ def from_formula(cls, formula, data, *, portfolios=None):
>>> formula = 'MktRF + SMB + HML'
>>> mod = TradedFactorModel.from_formula(formula, data, portfolios=portfolios)
"""
na_action = NAAction(on_NA='raise', NA_types=[])
orig_formula = formula
if portfolios is not None:
factors = dmatrix(formula + ' + 0', data, return_type='dataframe', NA_action=na_action)
else:
formula = formula.split('~')
portfolios = dmatrix(formula[0].strip() + ' + 0', data,
return_type='dataframe', NA_action=na_action)
factors = dmatrix(formula[1].strip() + ' + 0', data,
return_type='dataframe', NA_action=na_action)
factors, portfolios, formula = cls._prepare_data_from_formula(formula, data, portfolios)
mod = cls(portfolios, factors)
mod.formula = orig_formula
mod.formula = formula
return mod

def fit(self, cov_type='robust', debiased=True, **cov_config):
Expand Down Expand Up @@ -395,21 +401,12 @@ def from_formula(cls, formula, data, *, portfolios=None, risk_free=False,
>>> formula = 'MktRF + SMB + HML'
>>> mod = LinearFactorModel.from_formula(formula, data, portfolios=portfolios)
"""
na_action = NAAction(on_NA='raise', NA_types=[])
orig_formula = formula
if portfolios is not None:
factors = dmatrix(formula + ' + 0', data, return_type='dataframe', NA_action=na_action)
else:
formula = formula.split('~')
portfolios = dmatrix(formula[0].strip() + ' + 0', data,
return_type='dataframe', NA_action=na_action)
factors = dmatrix(formula[1].strip() + ' + 0', data,
return_type='dataframe', NA_action=na_action)
factors, portfolios, formula = cls._prepare_data_from_formula(formula, data, portfolios)
if sigma is not None:
mod = cls(portfolios, factors, risk_free=risk_free, sigma=sigma)
else:
mod = cls(portfolios, factors, risk_free=risk_free)
mod.formula = orig_formula
mod.formula = formula
return mod

def fit(self, cov_type='robust', debiased=True, **cov_config):
Expand Down Expand Up @@ -560,11 +557,11 @@ def _moments(self, eps, betas, lam, alphas, pricing_errors):
f = self.factors.ndarray
nobs, nf, nport, nrf, s1, s2, s3 = self._boundaries()
fc = np.c_[np.ones((nobs, 1)), f]
# Moments
f_rep = np.tile(fc, (1, nport))
eps_rep = np.tile(eps, (nf + 1, 1))
eps_rep = np.reshape(eps_rep.T, (nport * (nf + 1), nobs)).T

# Moments
g1 = f_rep * eps_rep
g2 = pricing_errors @ sigma_inv @ betas
g3 = pricing_errors - alphas.T
Expand Down

0 comments on commit 25aec32

Please sign in to comment.