Python utils.spdot函数代码示例

 def gradient(self, par):
     beta = np.reshape(np.array(par), (self.k, 1))
     q = 2 * self.y - 1
     qxb = q * spdot(self.x, beta)
     lamb = q * norm.pdf(qxb) / norm.cdf(qxb)
     gradient = spdot(lamb.T, self.x)[0]
     return gradient

 def par_est(self):
     start = np.dot(la.inv(spdot(self.x.T, self.x)),
                    spdot(self.x.T, self.y))
     flogl = lambda par: -self.ll(par)
     if self.optim == 'newton':
         fgrad = lambda par: self.gradient(par)
         fhess = lambda par: self.hessian(par)
         par_hat = newton(flogl, start, fgrad, fhess, self.maxiter)
         warn = par_hat[2]
     else:
         fgrad = lambda par: -self.gradient(par)
         if self.optim == 'bfgs':
             par_hat = op.fmin_bfgs(
                 flogl, start, fgrad, full_output=1, disp=0)
             warn = par_hat[6]
         if self.optim == 'ncg':
             fhess = lambda par: -self.hessian(par)
             par_hat = op.fmin_ncg(
                 flogl, start, fgrad, fhess=fhess, full_output=1, disp=0)
             warn = par_hat[5]
     if warn > 0:
         warn = True
     else:
         warn = False
     return par_hat, warn

    def __init__(self, y, x, yend, q=None, h=None,
                 robust=None, gwk=None, sig2n_k=False):

        if issubclass(type(q), np.ndarray) and issubclass(type(h), np.ndarray):
            raise Exception, "Please do not provide 'q' and 'h' together"
        if q == None and h == None:
            raise Exception, "Please provide either 'q' or 'h'"

        self.y = y
        self.n = y.shape[0]
        self.x = x

        self.kstar = yend.shape[1]
        # including exogenous and endogenous variables
        z = sphstack(self.x, yend)
        if type(h).__name__ not in ['ndarray', 'csr_matrix']:
            # including exogenous variables and instrument
            h = sphstack(self.x, q)
        self.z = z
        self.h = h
        self.q = q
        self.yend = yend
        # k = number of exogenous variables and endogenous variables
        self.k = z.shape[1]
        hth = spdot(h.T, h)
        hthi = la.inv(hth)
        zth = spdot(z.T, h)
        hty = spdot(h.T, y)

        factor_1 = np.dot(zth, hthi)
        factor_2 = np.dot(factor_1, zth.T)
        # this one needs to be in cache to be used in AK
        varb = la.inv(factor_2)
        factor_3 = np.dot(varb, factor_1)
        betas = np.dot(factor_3, hty)
        self.betas = betas
        self.varb = varb
        self.zthhthi = factor_1

        # predicted values
        self.predy = spdot(z, betas)

        # residuals
        u = y - self.predy
        self.u = u

        # attributes used in property
        self.hth = hth     # Required for condition index
        self.hthi = hthi   # Used in error models
        self.htz = zth.T

        if robust:
            self.vm = ROBUST.robust_vm(reg=self, gwk=gwk, sig2n_k=sig2n_k)

        self._cache = {}
        if sig2n_k:
            self.sig2 = self.sig2n_k
        else:
            self.sig2 = self.sig2n

 def hessian(self, par):
     beta = np.reshape(np.array(par), (self.k, 1))
     q = 2 * self.y - 1
     xb = spdot(self.x, beta)
     qxb = q * xb
     lamb = q * norm.pdf(qxb) / norm.cdf(qxb)
     hessian = spdot(self.x.T, spbroadcast(self.x,-lamb * (lamb + xb)))
     return hessian

 def xb(self):
     try:
         return self._cache['xb']
     except AttributeError:
         self._cache = {}
         self._cache['xb'] = spdot(self.x, self.betas)
     except KeyError:
         self._cache['xb'] = spdot(self.x, self.betas)
     return self._cache['xb']

 def slopes_vm(self):
     if 'slopes_vm' not in self._cache:
         x = self.xmean
         b = self.betas
         dfdb = np.eye(self.k) - spdot(b.T, x) * spdot(b, x.T)
         slopes_vm = (self.scale ** 2) * \
             np.dot(np.dot(dfdb, self.vm), dfdb.T)
         self._cache['slopes_vm'] = slopes_vm[1:, 1:]
     return self._cache['slopes_vm']

 def _get_spat_diag_props(self,y, x, w, yend, q, w_lags, lag_q):
     self._cache = {}
     yend, q = set_endog(y, x, w, yend, q, w_lags, lag_q)
     x = USER.check_constant(x)
     x = REGI.regimeX_setup(x, self.regimes, [True] * x.shape[1], self.regimes_set)
     self.z = sphstack(x,REGI.regimeX_setup(yend, self.regimes, [True] * (yend.shape[1]-1)+[False], self.regimes_set))
     self.h = sphstack(x,REGI.regimeX_setup(q, self.regimes, [True] * q.shape[1], self.regimes_set))
     hthi = np.linalg.inv(spdot(self.h.T,self.h))
     zth = spdot(self.z.T,self.h)     
     self.varb = np.linalg.inv(spdot(spdot(zth,hthi),zth.T))

 def _get_spat_diag_props(self, results, regi_ids, x, yend, q):
     self._cache = {}
     x = USER.check_constant(x)
     x = REGI.regimeX_setup(
         x, self.regimes, [True] * x.shape[1], self.regimes_set)
     self.z = sphstack(x, REGI.regimeX_setup(
         yend, self.regimes, [True] * yend.shape[1], self.regimes_set))
     self.h = sphstack(
         x, REGI.regimeX_setup(q, self.regimes, [True] * q.shape[1], self.regimes_set))
     hthi = np.linalg.inv(spdot(self.h.T, self.h))
     zth = spdot(self.z.T, self.h)
     self.varb = np.linalg.inv(spdot(spdot(zth, hthi), zth.T))

 def AB(self):
     """
     Computes A and B matrices as in Cliff-Ord 1981, p. 203
     """
     if 'AB' not in self._cache:
         U = (self.w.sparse + self.w.sparse.T) / 2.
         z = spdot(U, self.reg.x, array_out=False)
         c1 = spdot(self.reg.x.T, z, array_out=False)
         c2 = spdot(z.T, z, array_out=False)
         G = self.reg.xtxi
         A = spdot(G, c1)
         B = spdot(G, c2)
         self._cache['AB'] = [A, B]
     return self._cache['AB']

    def __init__(self, y, x, w):

        #1a. OLS --> \tilde{betas}
        ols = OLS.BaseOLS(y=y, x=x)
        self.n, self.k = ols.x.shape
        self.x = ols.x
        self.y = ols.y

        #1b. GMM --> \tilde{\lambda1}
        moments = _momentsGM_Error(w, ols.u)
        lambda1 = optim_moments(moments)

        #2a. OLS -->\hat{betas}
        xs = get_spFilter(w, lambda1, self.x)
        ys = get_spFilter(w, lambda1, self.y)
        ols2 = OLS.BaseOLS(y=ys, x=xs)

        #Output
        self.predy = spdot(self.x, ols2.betas)
        self.u = y - self.predy
        self.betas = np.vstack((ols2.betas, np.array([[lambda1]])))
        self.sig2 = ols2.sig2n
        self.e_filtered = self.u - lambda1*w*self.u

        self.vm = self.sig2 * ols2.xtxi
        se_betas = np.sqrt(self.vm.diagonal())
        self._cache = {}

 def _get_spat_diag_props(self, x, sig2n_k):
     self.k = self.kr
     self._cache = {}
     x = np.hstack((np.ones((x.shape[0], 1)), x))
     self.x = REGI.regimeX_setup(x, self.regimes, [True] * x.shape[1], self.regimes_set)
     self.xtx = spdot(self.x.T, self.x)
     self.xtxi = np.linalg.inv(self.xtx)                

def _optimal_weight(reg, sig2n_k, warn=True):
    try:
        Hu = reg.h.toarray() * reg.u ** 2
    except:
        Hu = reg.h * reg.u ** 2
    if sig2n_k:
        S = spdot(reg.h.T, Hu, array_out=True) / (reg.n - reg.k)
    else:
        S = spdot(reg.h.T, Hu, array_out=True) / reg.n
    Si = np.linalg.inv(S)
    ZtH = spdot(reg.z.T, reg.h)
    ZtHSi = spdot(ZtH, Si)
    fac2 = np.linalg.inv(spdot(ZtHSi, ZtH.T, array_out=True))
    fac3 = spdot(ZtHSi, spdot(reg.h.T, reg.y), array_out=True)
    betas = np.dot(fac2, fac3)
    if sig2n_k:
        vm = fac2 * (reg.n - reg.k)
    else:
        vm = fac2 * reg.n
    RegressionProps_basic(reg, betas=betas, vm=vm, sig2=False)
    reg.title += " (Optimal-Weighted GMM)"
    if warn:
        set_warn(
            reg, "Residuals treated as homoskedastic for the purpose of diagnostics.")
    return

 def _optimal_weight(self,sig2n_k):
     ''' Computes optimal weights GMM. '''
     H = spdot(spdot(self.h,self.hthi),self.htz)
     fac2, ZtHSi = self._get_fac2_het(H,self.u,sig2n_k)
     fac3 = spdot(ZtHSi,spdot(H.T,self.y),array_out=True)
     betas = np.dot(fac2,fac3)
     """
     # Updating S_hat to S_tilde
     predy = spdot(self.z, betas)
     u = self.y - predy
     fac2, ZtHSi = self._get_fac2_het(u,sig2n_k)
     """
     if sig2n_k:
         vm = fac2*(self.n-self.k)
     else:
         vm = fac2*self.n
     #return betas, predy, u, vm #Use this line instead of next if updating S_hat to S_tilde.
     return betas, vm

 def slopes_vm(self):
     try:
         return self._cache['slopes_vm']
     except AttributeError:
         self._cache = {}
         x = self.xmean
         b = self.betas
         dfdb = np.eye(self.k) - spdot(b.T, x) * spdot(b, x.T)
         slopes_vm = (self.scale ** 2) * \
             np.dot(np.dot(dfdb, self.vm), dfdb.T)
         self._cache['slopes_vm'] = slopes_vm[1:, 1:]
     except KeyError:
         x = self.xmean
         b = self.betas
         dfdb = np.eye(self.k) - spdot(b.T, x) * spdot(b, x.T)
         slopes_vm = (self.scale ** 2) * \
             np.dot(np.dot(dfdb, self.vm), dfdb.T)
         self._cache['slopes_vm'] = slopes_vm[1:, 1:]
     return self._cache['slopes_vm']

 def j(self):
     if 'j' not in self._cache:
         wxb = self.w.sparse * self.reg.predy
         wxb2 = np.dot(wxb.T, wxb)
         xwxb = spdot(self.reg.x.T, wxb)
         num1 = wxb2 - np.dot(xwxb.T, np.dot(self.reg.xtxi, xwxb))
         num = num1 + (self.t * self.reg.sig2n)
         den = self.reg.n * self.reg.sig2n
         self._cache['j'] = num / den
     return self._cache['j']