"""Tests for Regression Diagnostics and Specification Tests

Created on Thu Feb 09 13:19:47 2012

Author: Josef Perktold
License: BSD-3

currently all tests are against R

"""
import json
import os

import numpy as np
from numpy.testing import (
    assert_,
    assert_allclose,
    assert_almost_equal,
    assert_array_equal,
    assert_equal,
)
import pandas as pd
from pandas.testing import assert_frame_equal
import pytest

from statsmodels.datasets import macrodata, sunspots
from statsmodels.regression.linear_model import OLS
import statsmodels.stats.diagnostic as smsdia
import statsmodels.stats.outliers_influence as oi
import statsmodels.stats.sandwich_covariance as sw
from statsmodels.tools.tools import Bunch, add_constant
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.tsa.arima.model import ARIMA

cur_dir = os.path.abspath(os.path.dirname(__file__))


@pytest.fixture(scope="module")
def diagnostic_data():
    rs = np.random.RandomState(93674328)
    e = rs.standard_normal(500)
    x = rs.standard_normal((500, 3))
    y = x.sum(1) + e
    c = np.ones_like(y)
    data = pd.DataFrame(np.c_[y, c, x], columns=["y", "c", "x1", "x2", "x3"])

    return data


def compare_to_reference(sp, sp_dict, decimal=(12, 12)):
    assert_allclose(
        sp[0],
        sp_dict["statistic"],
        atol=10 ** -decimal[0],
        rtol=10 ** -decimal[0],
    )
    assert_allclose(
        sp[1],
        sp_dict["pvalue"],
        atol=10 ** -decimal[1],
        rtol=10 ** -decimal[0],
    )


def test_gq():
    d = macrodata.load().data

    realinv = d["realinv"]
    realgdp = d["realgdp"]
    realint = d["realint"]
    endog = realinv
    exog = add_constant(np.c_[realgdp, realint])

    # goldfeld-quandt
    het_gq_greater = dict(
        statistic=13.20512768685082,
        df1=99,
        df2=98,
        pvalue=1.246141976112324e-30,
        distr="f",
    )

    gq = smsdia.het_goldfeldquandt(endog, exog, split=0.5)
    compare_to_reference(gq, het_gq_greater, decimal=(12, 12))
    assert_equal(gq[-1], "increasing")


class TestDiagnosticG:
    @classmethod
    def setup_class(cls):
        d = macrodata.load_pandas().data
        # growth rates
        gs_l_realinv = 400 * np.diff(np.log(d["realinv"].values))
        gs_l_realgdp = 400 * np.diff(np.log(d["realgdp"].values))
        lint = d["realint"][:-1].values
        tbilrate = d["tbilrate"][:-1].values

        endogg = gs_l_realinv
        exogg = add_constant(np.c_[gs_l_realgdp, lint])
        exogg2 = add_constant(np.c_[gs_l_realgdp, tbilrate])
        exogg3 = add_constant(np.c_[gs_l_realgdp])

        res_ols = OLS(endogg, exogg).fit()
        res_ols2 = OLS(endogg, exogg2).fit()

        res_ols3 = OLS(endogg, exogg3).fit()

        cls.res = res_ols
        cls.res2 = res_ols2
        cls.res3 = res_ols3
        cls.endog = cls.res.model.endog
        cls.exog = cls.res.model.exog

    def test_basic(self):
        # mainly to check I got the right regression
        # > mkarray(fm$coefficients, "params")
        params = np.array(
            [-9.48167277465485, 4.3742216647032, -0.613996969478989]
        )

        assert_almost_equal(self.res.params, params, decimal=12)

    def test_hac(self):
        res = self.res
        # > nw = NeweyWest(fm, lag = 4, prewhite = FALSE, verbose=TRUE)
        # > nw2 = NeweyWest(fm, lag=10, prewhite = FALSE, verbose=TRUE)

        # > mkarray(nw, "cov_hac_4")
        cov_hac_4 = np.array(
            [
                1.385551290884014,
                -0.3133096102522685,
                -0.0597207976835705,
                -0.3133096102522685,
                0.1081011690351306,
                0.000389440793564336,
                -0.0597207976835705,
                0.000389440793564339,
                0.0862118527405036,
            ]
        ).reshape((3, 3), order="F")

        # > mkarray(nw2, "cov_hac_10")
        cov_hac_10 = np.array(
            [
                1.257386180080192,
                -0.2871560199899846,
                -0.03958300024627573,
                -0.2871560199899845,
                0.1049107028987101,
                0.0003896205316866944,
                -0.03958300024627578,
                0.0003896205316866961,
                0.0985539340694839,
            ]
        ).reshape((3, 3), order="F")

        cov = sw.cov_hac_simple(res, nlags=4, use_correction=False)
        bse_hac = sw.se_cov(cov)
        assert_almost_equal(cov, cov_hac_4, decimal=12)
        assert_almost_equal(bse_hac, np.sqrt(np.diag(cov)), decimal=12)

        cov = sw.cov_hac_simple(res, nlags=10, use_correction=False)
        bse_hac = sw.se_cov(cov)
        assert_almost_equal(cov, cov_hac_10, decimal=12)
        assert_almost_equal(bse_hac, np.sqrt(np.diag(cov)), decimal=12)

    def test_het_goldfeldquandt(self):
        # TODO: test options missing

        # > gq = gqtest(fm, alternative='greater')
        # > mkhtest_f(gq, 'het_gq_greater', 'f')
        het_gq_greater = dict(
            statistic=0.5313259064778423,
            pvalue=0.9990217851193723,
            parameters=(98, 98),
            distr="f",
        )

        # > gq = gqtest(fm, alternative='less')
        # > mkhtest_f(gq, 'het_gq_less', 'f')
        het_gq_less = dict(
            statistic=0.5313259064778423,
            pvalue=0.000978214880627621,
            parameters=(98, 98),
            distr="f",
        )

        # > gq = gqtest(fm, alternative='two.sided')
        # > mkhtest_f(gq, 'het_gq_two_sided', 'f')
        het_gq_two_sided = dict(
            statistic=0.5313259064778423,
            pvalue=0.001956429761255241,
            parameters=(98, 98),
            distr="f",
        )

        # > gq = gqtest(fm, fraction=0.1, alternative='two.sided')
        # > mkhtest_f(gq, 'het_gq_two_sided_01', 'f')
        het_gq_two_sided_01 = dict(
            statistic=0.5006976835928314,
            pvalue=0.001387126702579789,
            parameters=(88, 87),
            distr="f",
        )

        endogg, exogg = self.endog, self.exog
        # tests
        gq = smsdia.het_goldfeldquandt(endogg, exogg, split=0.5)
        compare_to_reference(gq, het_gq_greater, decimal=(12, 12))
        assert_equal(gq[-1], "increasing")

        gq = smsdia.het_goldfeldquandt(
            endogg, exogg, split=0.5, alternative="decreasing"
        )
        compare_to_reference(gq, het_gq_less, decimal=(12, 12))
        assert_equal(gq[-1], "decreasing")

        gq = smsdia.het_goldfeldquandt(
            endogg, exogg, split=0.5, alternative="two-sided"
        )
        compare_to_reference(gq, het_gq_two_sided, decimal=(12, 12))
        assert_equal(gq[-1], "two-sided")

        # TODO: forcing the same split as R 202-90-90-1=21
        gq = smsdia.het_goldfeldquandt(
            endogg, exogg, split=90, drop=21, alternative="two-sided"
        )
        compare_to_reference(gq, het_gq_two_sided_01, decimal=(12, 12))
        assert_equal(gq[-1], "two-sided")
        # TODO other options ???

    def test_het_breusch_pagan(self):
        res = self.res

        bptest = dict(
            statistic=0.709924388395087,
            pvalue=0.701199952134347,
            parameters=(2,),
            distr="f",
        )

        bp = smsdia.het_breuschpagan(res.resid, res.model.exog)
        compare_to_reference(bp, bptest, decimal=(12, 12))

    def test_het_breusch_pagan_1d_err(self):
        res = self.res
        x = np.asarray(res.model.exog)[:, -1]
        with pytest.raises(ValueError, match="The Breusch-Pagan"):
            smsdia.het_breuschpagan(res.resid, x)
        x = np.ones_like(x)
        with pytest.raises(ValueError, match="The Breusch-Pagan"):
            smsdia.het_breuschpagan(res.resid, x)
        x = np.asarray(res.model.exog).copy()
        x[:, 0] = 0
        with pytest.raises(ValueError, match="The Breusch-Pagan"):
            smsdia.het_breuschpagan(res.resid, x)

    def test_het_breusch_pagan_nonrobust(self):
        res = self.res

        bptest = dict(
            statistic=1.302014063483341,
            pvalue=0.5215203247110649,
            parameters=(2,),
            distr="f",
        )

        bp = smsdia.het_breuschpagan(res.resid, res.model.exog, robust=False)
        compare_to_reference(bp, bptest, decimal=(12, 12))

    def test_het_white(self):
        res = self.res

        # TODO: regressiontest, compare with Greene or Gretl or Stata
        hw = smsdia.het_white(res.resid, res.model.exog)
        hw_values = (
            33.503722896538441,
            2.9887960597830259e-06,
            7.7945101228430946,
            1.0354575277704231e-06,
        )
        assert_almost_equal(hw, hw_values)

    def test_het_white_error(self):
        res = self.res

        with pytest.raises(ValueError, match="White's heteroskedasticity"):
            smsdia.het_white(res.resid, res.model.exog[:, :1])

    def test_het_arch(self):
        # test het_arch and indirectly het_lm against R
        # > library(FinTS)
        # > at = ArchTest(residuals(fm), lags=4)
        # > mkhtest(at, 'archtest_4', 'chi2')
        archtest_4 = dict(
            statistic=3.43473400836259,
            pvalue=0.487871315392619,
            parameters=(4,),
            distr="chi2",
        )

        # > at = ArchTest(residuals(fm), lags=12)
        # > mkhtest(at, 'archtest_12', 'chi2')
        archtest_12 = dict(
            statistic=8.648320999014171,
            pvalue=0.732638635007718,
            parameters=(12,),
            distr="chi2",
        )

        at4 = smsdia.het_arch(self.res.resid, nlags=4)
        at12 = smsdia.het_arch(self.res.resid, nlags=12)
        compare_to_reference(at4[:2], archtest_4, decimal=(12, 13))
        compare_to_reference(at12[:2], archtest_12, decimal=(12, 13))

    def test_het_arch2(self):
        # test autolag options, this also test het_lm
        # unfortunately optimal lag=1 for this data
        resid = self.res.resid

        res1 = smsdia.het_arch(resid, nlags=5, store=True)
        rs1 = res1[-1]
        res2 = smsdia.het_arch(resid, nlags=5, store=True)
        rs2 = res2[-1]

        assert_almost_equal(rs2.resols.params, rs1.resols.params, decimal=12)
        assert_almost_equal(res2[:4], res1[:4], decimal=12)

        # test that smallest lag, nlags=1 works
        res3 = smsdia.het_arch(resid, nlags=5)
        assert_almost_equal(res3[:4], res1[:4], decimal=12)

    def test_acorr_breusch_godfrey(self):
        res = self.res

        # bgf = bgtest(fm, order = 4, type="F")
        breuschgodfrey_f = dict(
            statistic=1.179280833676792,
            pvalue=0.321197487261203,
            parameters=(
                4,
                195,
            ),
            distr="f",
        )

        # > bgc = bgtest(fm, order = 4, type="Chisq")
        # > mkhtest(bgc, "breuschpagan_c", "chi2")
        breuschgodfrey_c = dict(
            statistic=4.771042651230007,
            pvalue=0.3116067133066697,
            parameters=(4,),
            distr="chi2",
        )

        bg = smsdia.acorr_breusch_godfrey(res, nlags=4)
        bg_r = [
            breuschgodfrey_c["statistic"],
            breuschgodfrey_c["pvalue"],
            breuschgodfrey_f["statistic"],
            breuschgodfrey_f["pvalue"],
        ]
        assert_almost_equal(bg, bg_r, decimal=11)

        # check that lag choice works
        bg2 = smsdia.acorr_breusch_godfrey(res, nlags=None)
        bg3 = smsdia.acorr_breusch_godfrey(res, nlags=10)
        assert_almost_equal(bg2, bg3, decimal=12)

    def test_acorr_breusch_godfrey_multidim(self):
        res = Bunch(resid=np.empty((100, 2)))
        with pytest.raises(ValueError, match="Model resid must be a 1d array"):
            smsdia.acorr_breusch_godfrey(res)

    def test_acorr_breusch_godfrey_exogs(self):
        data = sunspots.load_pandas().data["SUNACTIVITY"]
        res = ARIMA(data, order=(1, 0, 0), trend="n").fit()
        smsdia.acorr_breusch_godfrey(res, nlags=1)

    def test_acorr_ljung_box(self):

        # unit-test which may be useful later
        # ddof correction for fitted parameters in ARMA(p,q) fitdf=p+q
        # bt = Box.test(residuals(fm), lag=4, type = "Ljung-Box", fitdf=2)
        # mkhtest(bt, "ljung_box_4df2", "chi2")
        # ljung_box_4df2 = dict(statistic=5.23587172795227,
        #                       pvalue=0.0729532930400377,
        #                       parameters=(2,), distr='chi2')

        # bt = Box.test(residuals(fm), lag=4, type = "Box-Pierce", fitdf=2)
        # mkhtest(bt, "ljung_box_bp_4df2", "chi2")
        # ljung_box_bp_4df2 = dict(statistic=5.12462932741681,
        #                          pvalue=0.0771260128929921,
        #                          parameters=(2,), distr='chi2')

        res = self.res
        # general test
        # bt = Box.test(residuals(fm), lag=4, type = "Ljung-Box")
        # mkhtest(bt, "ljung_box_4", "chi2")
        ljung_box_4 = dict(
            statistic=5.23587172795227,
            pvalue=0.263940335284713,
            parameters=(4,),
            distr="chi2",
        )

        # bt = Box.test(residuals(fm), lag=4, type = "Box-Pierce")
        # mkhtest(bt, "ljung_box_bp_4", "chi2")
        ljung_box_bp_4 = dict(
            statistic=5.12462932741681,
            pvalue=0.2747471266820692,
            parameters=(4,),
            distr="chi2",
        )

        df = smsdia.acorr_ljungbox(res.resid, 4, boxpierce=True)
        compare_to_reference(
            [df.loc[4, "lb_stat"], df.loc[4, "lb_pvalue"]],
            ljung_box_4,
            decimal=(12, 12),
        )
        compare_to_reference(
            [df.loc[4, "bp_stat"], df.loc[4, "bp_pvalue"]],
            ljung_box_bp_4,
            decimal=(12, 12),
        )

    def test_acorr_ljung_box_big_default(self):
        res = self.res
        # R test with big dataset and default lag
        # bt = Box.test(residuals(fm), type = "Ljung-Box")
        # mkhtest(bt, "ljung_box_none", "chi2")
        ljung_box_none = dict(
            statistic=51.03724531797195, pvalue=0.11334744923390, distr="chi2"
        )

        # bt = Box.test(residuals(fm), type = "Box-Pierce")
        # mkhtest(bt, "ljung_box_bp_none", "chi2")
        ljung_box_bp_none = dict(
            statistic=45.12238537034000, pvalue=0.26638168491464, distr="chi2"
        )
        lags = min(40, res.resid.shape[0] // 2 - 2)
        df = smsdia.acorr_ljungbox(res.resid, boxpierce=True, lags=lags)
        idx = df.index.max()
        compare_to_reference(
            [df.loc[idx, "lb_stat"], df.loc[idx, "lb_pvalue"]],
            ljung_box_none,
            decimal=(12, 12),
        )
        compare_to_reference(
            [df.loc[idx, "bp_stat"], df.loc[idx, "bp_pvalue"]],
            ljung_box_bp_none,
            decimal=(12, 12),
        )

    def test_acorr_ljung_box_small_default(self):
        res = self.res
        # R test with small dataset and default lag
        # bt = Box.test(residuals(fm), type = "Ljung-Box")
        # mkhtest(bt, "ljung_box_small", "chi2")
        ljung_box_small = dict(
            statistic=9.61503968281915,
            pvalue=0.72507000996945,
            parameters=(0,),
            distr="chi2",
        )

        # bt = Box.test(residuals(fm), type = "Box-Pierce")
        # mkhtest(bt, "ljung_box_bp_small", "chi2")
        ljung_box_bp_small = dict(
            statistic=7.41692150864936,
            pvalue=0.87940785887006,
            parameters=(0,),
            distr="chi2",
        )
        if isinstance(res.resid, np.ndarray):
            resid = res.resid[:30]
        else:
            resid = res.resid.iloc[:30]
        df = smsdia.acorr_ljungbox(
            resid, boxpierce=True, lags=13
        )
        idx = df.index.max()
        compare_to_reference(
            [df.loc[idx, "lb_stat"], df.loc[idx, "lb_pvalue"]],
            ljung_box_small,
            decimal=(12, 12),
        )
        compare_to_reference(
            [df.loc[idx, "bp_stat"], df.loc[idx, "bp_pvalue"]],
            ljung_box_bp_small,
            decimal=(12, 12),
        )

    def test_acorr_ljung_box_against_r(self, reset_randomstate):
        rs = np.random.RandomState(9876543)
        y1 = rs.standard_normal(100)
        e = rs.standard_normal(201)
        y2 = np.zeros_like(e)
        y2[0] = e[0]
        for i in range(1, 201):
            y2[i] = 0.5 * y2[i - 1] - 0.4 * e[i - 1] + e[i]
        y2 = y2[-100:]

        r_results_y1_lb = [
            [0.15685, 1, 0.6921],
            [5.4737, 5, 0.3608],
            [10.508, 10, 0.3971],
        ]

        r_results_y2_lb = [
            [2.8764, 1, 0.08989],
            [3.8104, 5, 0.577],
            [8.4779, 10, 0.5823],
        ]
        res_y1 = smsdia.acorr_ljungbox(y1, 10)
        res_y2 = smsdia.acorr_ljungbox(y2, 10)
        for i, loc in enumerate((1, 5, 10)):
            row = res_y1.loc[loc]
            assert_allclose(
                r_results_y1_lb[i][0], row.loc["lb_stat"], rtol=1e-3
            )
            assert_allclose(
                r_results_y1_lb[i][2], row.loc["lb_pvalue"], rtol=1e-3
            )

            row = res_y2.loc[loc]
            assert_allclose(
                r_results_y2_lb[i][0], row.loc["lb_stat"], rtol=1e-3
            )
            assert_allclose(
                r_results_y2_lb[i][2], row.loc["lb_pvalue"], rtol=1e-3
            )

        res = smsdia.acorr_ljungbox(y2, 10, boxpierce=True)
        assert_allclose(res.loc[10, "bp_stat"], 7.8935, rtol=1e-3)
        assert_allclose(res.loc[10, "bp_pvalue"], 0.639, rtol=1e-3)

        res = smsdia.acorr_ljungbox(y2, 10, boxpierce=True, model_df=1)
        assert_allclose(res.loc[10, "bp_pvalue"], 0.5449, rtol=1e-3)

    def test_harvey_collier(self):
        # > hc = harvtest(fm, order.by = NULL, data = list())
        # > mkhtest_f(hc, 'harvey_collier', 't')
        harvey_collier = dict(
            statistic=0.494432160939874,
            pvalue=0.6215491310408242,
            parameters=198,
            distr="t",
        )

        hc = smsdia.linear_harvey_collier(self.res)
        compare_to_reference(hc, harvey_collier, decimal=(12, 12))
        hc_skip = smsdia.linear_harvey_collier(self.res, skip=20)
        assert not np.allclose(hc[0], hc_skip[0])

    def test_rainbow(self):
        # rainbow test
        # > rt = raintest(fm)
        # > mkhtest_f(rt, 'raintest', 'f')
        raintest = dict(
            statistic=0.6809600116739604,
            pvalue=0.971832843583418,
            parameters=(101, 98),
            distr="f",
        )

        # > rt = raintest(fm, fraction=0.4)
        # > mkhtest_f(rt, 'raintest_fraction_04', 'f')
        raintest_fraction_04 = dict(
            statistic=0.565551237772662,
            pvalue=0.997592305968473,
            parameters=(122, 77),
            distr="f",
        )

        rb = smsdia.linear_rainbow(self.res)
        compare_to_reference(rb, raintest, decimal=(12, 12))
        rb = smsdia.linear_rainbow(self.res, frac=0.4)
        compare_to_reference(rb, raintest_fraction_04, decimal=(12, 12))

    def test_compare_lr(self):
        res = self.res
        res3 = self.res3  # nested within res
        # lrtest
        # lrt = lrtest(fm, fm2)
        # Model 1: ginv ~ ggdp + lint
        # Model 2: ginv ~ ggdp

        lrtest = dict(
            loglike1=-763.9752181602237,
            loglike2=-766.3091902020184,
            chi2value=4.66794408358942,
            pvalue=0.03073069384028677,
            df=(4, 3, 1),
        )
        lrt = res.compare_lr_test(res3)
        assert_almost_equal(lrt[0], lrtest["chi2value"], decimal=11)
        assert_almost_equal(lrt[1], lrtest["pvalue"], decimal=11)

        waldtest = dict(
            fvalue=4.65216373312492,
            pvalue=0.03221346195239025,
            df=(199, 200, 1),
        )

        wt = res.compare_f_test(res3)
        assert_almost_equal(wt[0], waldtest["fvalue"], decimal=11)
        assert_almost_equal(wt[1], waldtest["pvalue"], decimal=11)

    def test_compare_j(self):
        res = self.res
        res2 = self.res2
        # jt = jtest(fm, lm(ginv ~ ggdp + tbilrate))
        # Estimate         Std. Error  t value Pr(>|t|)
        jtest = [
            (
                "M1 + fitted(M2)",
                1.591505670785873,
                0.7384552861695823,
                2.155182176352370,
                0.032354572525314450,
                "*",
            ),
            (
                "M2 + fitted(M1)",
                1.305687653016899,
                0.4808385176653064,
                2.715438978051544,
                0.007203854534057954,
                "**",
            ),
        ]

        jt1 = smsdia.compare_j(res2, res)
        assert_almost_equal(jt1, jtest[0][3:5], decimal=12)

        jt2 = smsdia.compare_j(res, res2)
        assert_almost_equal(jt2, jtest[1][3:5], decimal=12)

    @pytest.mark.parametrize("comp", [smsdia.compare_cox, smsdia.compare_j])
    def test_compare_error(self, comp, diagnostic_data):
        data = diagnostic_data
        res1 = OLS(data.y, data[["c", "x1"]]).fit()
        res2 = OLS(data.x3, data[["c", "x2"]]).fit()
        with pytest.raises(ValueError, match="endogenous variables"):
            comp(res1, res2)

    @pytest.mark.parametrize("comp", [smsdia.compare_cox, smsdia.compare_j])
    def test_compare_nested(self, comp, diagnostic_data):
        data = diagnostic_data
        res1 = OLS(data.y, data[["c", "x1"]]).fit()
        res2 = OLS(data.y, data[["c", "x1", "x2"]]).fit()
        with pytest.raises(ValueError, match="The exog in results_x"):
            comp(res1, res2)
        with pytest.raises(ValueError, match="The exog in results_x"):
            comp(res2, res1)

    def test_compare_cox(self):
        res = self.res
        res2 = self.res2
        # Estimate        Std. Error  z value   Pr(>|z|)
        coxtest = [
            (
                "fitted(M1) ~ M2",
                -0.782030488930356,
                0.599696502782265,
                -1.304043770977755,
                1.922186587840554e-01,
                " ",
            ),
            (
                "fitted(M2) ~ M1",
                -2.248817107408537,
                0.392656854330139,
                -5.727181590258883,
                1.021128495098556e-08,
                "***",
            ),
        ]

        ct1 = smsdia.compare_cox(res, res2)
        assert_almost_equal(ct1, coxtest[0][3:5], decimal=12)

        ct2 = smsdia.compare_cox(res2, res)
        assert_almost_equal(ct2, coxtest[1][3:5], decimal=12)

        _, _, store = smsdia.compare_cox(res, res2, store=True)
        assert isinstance(store, smsdia.ResultsStore)

    def test_cusum_ols(self):
        # R library(strucchange)
        # > sc = sctest(ginv ~ ggdp + lint, type="OLS-CUSUM")
        # > mkhtest(sc, 'cusum_ols', 'BB')
        cusum_ols = dict(
            statistic=1.055750610401214,
            pvalue=0.2149567397376543,
            parameters=(),
            distr="BB",
        )  # Brownian Bridge

        k_vars = 3
        cs_ols = smsdia.breaks_cusumolsresid(self.res.resid, ddof=k_vars)
        compare_to_reference(cs_ols, cusum_ols, decimal=(12, 12))

    def test_breaks_hansen(self):
        # > sc = sctest(ginv ~ ggdp + lint, type="Nyblom-Hansen")
        # > mkhtest(sc, 'breaks_nyblom_hansen', 'BB')
        breaks_nyblom_hansen = dict(
            statistic=1.0300792740544484,
            pvalue=0.1136087530212015,
            parameters=(),
            distr="BB",
        )

        bh = smsdia.breaks_hansen(self.res)
        assert_almost_equal(
            bh[0], breaks_nyblom_hansen["statistic"], decimal=12
        )
        # TODO: breaks_hansen does not return pvalues

    def test_recursive_residuals(self):

        reccumres_standardize = np.array(
            [
                -2.151,
                -3.748,
                -3.114,
                -3.096,
                -1.865,
                -2.230,
                -1.194,
                -3.500,
                -3.638,
                -4.447,
                -4.602,
                -4.631,
                -3.999,
                -4.830,
                -5.429,
                -5.435,
                -6.554,
                -8.093,
                -8.567,
                -7.532,
                -7.079,
                -8.468,
                -9.320,
                -12.256,
                -11.932,
                -11.454,
                -11.690,
                -11.318,
                -12.665,
                -12.842,
                -11.693,
                -10.803,
                -12.113,
                -12.109,
                -13.002,
                -11.897,
                -10.787,
                -10.159,
                -9.038,
                -9.007,
                -8.634,
                -7.552,
                -7.153,
                -6.447,
                -5.183,
                -3.794,
                -3.511,
                -3.979,
                -3.236,
                -3.793,
                -3.699,
                -5.056,
                -5.724,
                -4.888,
                -4.309,
                -3.688,
                -3.918,
                -3.735,
                -3.452,
                -2.086,
                -6.520,
                -7.959,
                -6.760,
                -6.855,
                -6.032,
                -4.405,
                -4.123,
                -4.075,
                -3.235,
                -3.115,
                -3.131,
                -2.986,
                -1.813,
                -4.824,
                -4.424,
                -4.796,
                -4.000,
                -3.390,
                -4.485,
                -4.669,
                -4.560,
                -3.834,
                -5.507,
                -3.792,
                -2.427,
                -1.756,
                -0.354,
                1.150,
                0.586,
                0.643,
                1.773,
                -0.830,
                -0.388,
                0.517,
                0.819,
                2.240,
                3.791,
                3.187,
                3.409,
                2.431,
                0.668,
                0.957,
                -0.928,
                0.327,
                -0.285,
                -0.625,
                -2.316,
                -1.986,
                -0.744,
                -1.396,
                -1.728,
                -0.646,
                -2.602,
                -2.741,
                -2.289,
                -2.897,
                -1.934,
                -2.532,
                -3.175,
                -2.806,
                -3.099,
                -2.658,
                -2.487,
                -2.515,
                -2.224,
                -2.416,
                -1.141,
                0.650,
                -0.947,
                0.725,
                0.439,
                0.885,
                2.419,
                2.642,
                2.745,
                3.506,
                4.491,
                5.377,
                4.624,
                5.523,
                6.488,
                6.097,
                5.390,
                6.299,
                6.656,
                6.735,
                8.151,
                7.260,
                7.846,
                8.771,
                8.400,
                8.717,
                9.916,
                9.008,
                8.910,
                8.294,
                8.982,
                8.540,
                8.395,
                7.782,
                7.794,
                8.142,
                8.362,
                8.400,
                7.850,
                7.643,
                8.228,
                6.408,
                7.218,
                7.699,
                7.895,
                8.725,
                8.938,
                8.781,
                8.350,
                9.136,
                9.056,
                10.365,
                10.495,
                10.704,
                10.784,
                10.275,
                10.389,
                11.586,
                11.033,
                11.335,
                11.661,
                10.522,
                10.392,
                10.521,
                10.126,
                9.428,
                9.734,
                8.954,
                9.949,
                10.595,
                8.016,
                6.636,
                6.975,
            ]
        )

        rr = smsdia.recursive_olsresiduals(self.res, skip=3, alpha=0.95)
        assert_equal(
            np.round(rr[5][1:], 3), reccumres_standardize
        )  # extra zero in front
        # assert_equal(np.round(rr[3][4:], 3), np.diff(reccumres_standardize))
        assert_almost_equal(rr[3][4:], np.diff(reccumres_standardize), 3)
        assert_almost_equal(rr[4][3:].std(ddof=1), 10.7242, decimal=4)
        order = np.arange(self.res.model.endog.shape[0])
        rr_order = smsdia.recursive_olsresiduals(
            self.res, skip=3, alpha=0.95, order_by=order
        )
        assert_allclose(rr[0], rr_order[0], atol=1e-6)

        # regression number, visually checked with graph from gretl
        ub0 = np.array(
            [13.37318571, 13.50758959, 13.64199346, 13.77639734, 13.91080121]
        )
        ub1 = np.array(
            [39.44753774, 39.58194162, 39.7163455, 39.85074937, 39.98515325]
        )
        lb, ub = rr[6]
        assert_almost_equal(ub[:5], ub0, decimal=7)
        assert_almost_equal(lb[:5], -ub0, decimal=7)
        assert_almost_equal(ub[-5:], ub1, decimal=7)
        assert_almost_equal(lb[-5:], -ub1, decimal=7)

        # test a few values with explicit OLS
        endog = self.res.model.endog
        exog = self.res.model.exog
        params = []
        ypred = []
        for i in range(3, 10):
            resi = OLS(endog[:i], exog[:i]).fit()
            ypred.append(resi.model.predict(resi.params, exog[i]))
            params.append(resi.params)
        assert_almost_equal(rr[2][3:10], ypred, decimal=12)
        assert_almost_equal(rr[0][3:10], endog[3:10] - ypred, decimal=12)
        assert_almost_equal(rr[1][2:9], params, decimal=12)

    def test_normality(self):
        res = self.res

        # > library(nortest) # Lilliefors (Kolmogorov-Smirnov) normality test
        # > lt = lillie.test(residuals(fm))
        # > mkhtest(lt, "lilliefors", "-")
        lilliefors1 = dict(
            statistic=0.0723390908786589,
            pvalue=0.01204113540102896,
            parameters=(),
            distr="-",
        )

        # > lt = lillie.test(residuals(fm)**2)
        # > mkhtest(lt, "lilliefors", "-")
        lilliefors2 = dict(
            statistic=0.301311621898024,
            pvalue=1.004305736618051e-51,
            parameters=(),
            distr="-",
        )

        # > lt = lillie.test(residuals(fm)[1:20])
        # > mkhtest(lt, "lilliefors", "-")
        lilliefors3 = dict(
            statistic=0.1333956004203103,
            pvalue=0.455683,
            parameters=(),
            distr="-",
        )

        lf1 = smsdia.lilliefors(res.resid, pvalmethod="approx")
        lf2 = smsdia.lilliefors(res.resid ** 2, pvalmethod="approx")
        if isinstance(res.resid, np.ndarray):
            resid = res.resid[:20]
        else:
            resid = res.resid.iloc[:20]
        lf3 = smsdia.lilliefors(resid, pvalmethod="approx")

        compare_to_reference(lf1, lilliefors1, decimal=(12, 12))
        compare_to_reference(
            lf2, lilliefors2, decimal=(12, 12)
        )  # pvalue very small
        assert_allclose(lf2[1], lilliefors2["pvalue"], rtol=1e-10)
        compare_to_reference(lf3, lilliefors3, decimal=(12, 1))
        # R uses different approximation for pvalue in last case

        # > ad = ad.test(residuals(fm))
        # > mkhtest(ad, "ad3", "-")
        adr1 = dict(
            statistic=1.602209621518313,
            pvalue=0.0003937979149362316,
            parameters=(),
            distr="-",
        )

        # > ad = ad.test(residuals(fm)**2)
        # > mkhtest(ad, "ad3", "-")

        # > ad = ad.test(residuals(fm)[1:20])
        # > mkhtest(ad, "ad3", "-")
        adr3 = dict(
            statistic=0.3017073732210775,
            pvalue=0.5443499281265933,
            parameters=(),
            distr="-",
        )

        ad1 = smsdia.normal_ad(res.resid)
        compare_to_reference(ad1, adr1, decimal=(11, 13))
        ad2 = smsdia.normal_ad(res.resid ** 2)
        assert_(np.isinf(ad2[0]))
        ad3 = smsdia.normal_ad(resid)
