'Test for heteroskedasticity of standardized residuals
Tests whether the sum-of-squares in the first third of the sample is
significantly different than the sum-of-squares in the last third
of the sample. Analogous to a Goldfeld-Quandt test. The null hypothesis
is of no heteroskedasticity.
Parameters
method : string {\'breakvar\'} or None
The statistical test for heteroskedasticity. Must be \'breakvar\'
for test of a break in the variance. If None, an attempt is
made to select an appropriate test.
alternative : string, \'increasing\', \'decreasing\' or \'two-sided\'
This specifies the alternative for the p-value calculation. Default
is two-sided.
use_f : boolean, optional
Whether or not to compare against the asymptotic distribution
(chi-squared) or the approximate small-sample distribution (F).
Default is True (i.e. default is to compare against an F
distribution).
Returns
output : array
An array with `(test_statistic, pvalue)` for each endogenous
variable. The array is then sized `(k_endog, 2)`. If the method is
called as `het = res.test_heteroskedasticity()`, then `het[0]` is
an array of size 2 corresponding to the first endogenous variable,
where `het[0][0]` is the test statistic, and `het[0][1]` is the
p-value.
Notes
The null hypothesis is of no heteroskedasticity. That means different
things depending on which alternative is selected:
- Increasing: Null hypothesis is that the variance is not increasing
throughout the sample; that the sum-of-squares in the later
subsample is *not* greater than the sum-of-squares in the earlier
subsample.
- Decreasing: Null hypothesis is that the variance is not decreasing
throughout the sample; that the sum-of-squares in the earlier
subsample is *not* greater than the sum-of-squares in the later
subsample.
- Two-sided: Null hypothesis is that the variance is not changing
throughout the sample. Both that the sum-of-squares in the earlier
subsample is not greater than the sum-of-squares in the later
subsample *and* that the sum-of-squares in the later subsample is
not greater than the sum-of-squares in the earlier subsample.
For :math:`h = [T/3]`, the test statistic is:
.. math::
H(h) = \sum_{t=T-h+1}^T  \tilde v_t^2
\Bigg / \sum_{t=d+1}^{d+1+h} \tilde v_t^2
where :math:`d` is the number of periods in which the loglikelihood was
burned in the parent model (usually corresponding to diffuse
initialization).
This statistic can be tested against an :math:`F(h,h)` distribution.
Alternatively, :math:`h H(h)` is asymptotically distributed according
to :math:`\chi_h^2`; this second test can be applied by passing
`asymptotic=True` as an argument.
See section 5.4 of [1]_ for the above formula and discussion, as well
as additional details.
TODO
- Allow specification of :math:`h`
References
.. [1] Harvey, Andrew C. 1990.
Forecasting, Structural Time Series Models and the Kalman Filter.
Cambridge University Press.'

def test_heteroskedasticity(self, method, alternative='two-sided', use_f=True):

    if (method is None):
        method = 'breakvar'
    if (method == 'breakvar'):
        squared_resid = (self.filter_results.standardized_forecasts_error ** 2)
        d = self.loglikelihood_burn
        test_statistics = []
        p_values = []
        for i in range(self.model.k_endog):
            h = int(np.round((self.nobs_effective / 3)))
            numer_resid = squared_resid[i, (- h):]
            numer_resid = numer_resid[(~ np.isnan(numer_resid))]
            numer_dof = len(numer_resid)
            denom_resid = squared_resid[i, d:(d + h)]
            denom_resid = denom_resid[(~ np.isnan(denom_resid))]
            denom_dof = len(denom_resid)
            if (numer_dof < 2):
                warnings.warn(('Early    subset    of    data    for    variable    %d        has    too    few    non-missing    observations    to    calculate    test    statistic.' % i))
                numer_resid = np.nan
            if (denom_dof < 2):
                warnings.warn(('Later    subset    of    data    for    variable    %d        has    too    few    non-missing    observations    to    calculate    test    statistic.' % i))
                denom_resid = np.nan
            test_statistic = (np.sum(numer_resid) / np.sum(denom_resid))
            if use_f:
                from scipy.stats import f
                pval_lower = (lambda test_statistics: f.cdf(test_statistics, numer_dof, denom_dof))
                pval_upper = (lambda test_statistics: f.sf(test_statistics, numer_dof, denom_dof))
            else:
                from scipy.stats import chi2
                pval_lower = (lambda test_statistics: chi2.cdf((numer_dof * test_statistics), denom_dof))
                pval_upper = (lambda test_statistics: chi2.sf((numer_dof * test_statistics), denom_dof))
            alternative = alternative.lower()
            if (alternative in ['i', 'inc', 'increasing']):
                p_value = pval_upper(test_statistic)
            elif (alternative in ['d', 'dec', 'decreasing']):
                test_statistic = (1.0 / test_statistic)
                p_value = pval_upper(test_statistic)
            elif (alternative in ['2', '2-sided', 'two-sided']):
                p_value = (2 * np.minimum(pval_lower(test_statistic), pval_upper(test_statistic)))
            else:
                raise ValueError('Invalid    alternative.')
            test_statistics.append(test_statistic)
            p_values.append(p_value)
        output = np.c_[(test_statistics, p_values)]
    else:
        raise NotImplementedError('Invalid    heteroskedasticity    test    method.')
    return output
