'Score
Parameters
params : array_like, optional
Array of parameters at which to evaluate the loglikelihood
function.
**kwargs
Additional keyword arguments to pass to the Kalman filter. See
`KalmanFilter.filter` for more details.
Notes
This method is from Harvey (1989), section 3.4.5
References
Harvey, Andrew C. 1990.
Forecasting, Structural Time Series Models and the Kalman Filter.
Cambridge University Press.'

def _score_obs_harvey(self, params, approx_complex_step=True, approx_centered=False, **kwargs):

    params = np.array(params, ndmin=1)
    n = len(params)
    self.update(params, transformed=True, complex_step=approx_complex_step)
    if approx_complex_step:
        kwargs['inversion_method'] = (INVERT_UNIVARIATE | SOLVE_LU)
    res = self.ssm.filter(complex_step=approx_complex_step, **kwargs)
    (partials_forecasts_error, partials_forecasts_error_cov) = self._forecasts_error_partial_derivatives(params, transformed=True, approx_complex_step=approx_complex_step, approx_centered=approx_centered, res=res, **kwargs)
    partials = np.zeros((self.nobs, n))
    k_endog = self.k_endog
    for t in range(self.nobs):
        for i in range(n):
            inv_forecasts_error_cov = np.linalg.inv(res.forecasts_error_cov[:, :, t])
            partials[(t, i)] += np.trace(np.dot(np.dot(inv_forecasts_error_cov, partials_forecasts_error_cov[:, :, t, i]), (np.eye(k_endog) - np.dot(inv_forecasts_error_cov, np.outer(res.forecasts_error[:, t], res.forecasts_error[:, t])))))
            partials[(t, i)] += (2 * np.dot(partials_forecasts_error[:, t, i], np.dot(inv_forecasts_error_cov, res.forecasts_error[:, t])))
    return ((- partials) / 2.0)
