from statsmodels.compat.numpy import lstsq
from statsmodels.compat.pandas import MONTH_END, YEAR_END, assert_index_equal
from statsmodels.compat.platform import PLATFORM_WIN
from statsmodels.compat.python import PYTHON_IMPL_WASM, lrange

import os
import warnings

import numpy as np
from numpy.testing import (
    assert_allclose,
    assert_almost_equal,
    assert_equal,
)
import pandas as pd
from pandas import DataFrame, Series, date_range
import pytest
from scipy import stats
from scipy.interpolate import interp1d

from statsmodels.datasets import macrodata, modechoice, nile, randhie, sunspots
from statsmodels.tools.sm_exceptions import (
    CollinearityWarning,
    InfeasibleTestError,
    InterpolationWarning,
    MissingDataError,
    SingularMatrixWarning,
)

# Remove imports when range unit root test gets an R implementation
from statsmodels.tools.validation import array_like, bool_like
from statsmodels.tsa.arima_process import arma_acovf
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import (
    acf,
    acovf,
    adfuller,
    arma_order_select_ic,
    breakvar_heteroskedasticity_test,
    ccf,
    ccovf,
    coint,
    grangercausalitytests,
    innovations_algo,
    innovations_filter,
    kpss,
    levinson_durbin,
    levinson_durbin_pacf,
    leybourne,
    pacf,
    pacf_burg,
    pacf_ols,
    pacf_yw,
    range_unit_root_test,
    zivot_andrews,
)

DECIMAL_8 = 8
DECIMAL_6 = 6
DECIMAL_5 = 5
DECIMAL_4 = 4
DECIMAL_3 = 3
DECIMAL_2 = 2
DECIMAL_1 = 1

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


@pytest.fixture(scope="module")
def acovf_data():
    rnd = np.random.RandomState(12345)
    return rnd.randn(250)


@pytest.fixture(scope="module")
def gc_data():
    mdata = macrodata.load_pandas().data
    mdata = mdata[["realgdp", "realcons"]].values
    data = mdata.astype(float)
    return np.diff(np.log(data), axis=0)


class CheckADF:
    """
    Test Augmented Dickey-Fuller

    Test values taken from Stata.
    """

    levels = ["1%", "5%", "10%"]
    data = macrodata.load_pandas()
    x = data.data["realgdp"].values
    y = data.data["infl"].values

    def test_teststat(self):
        assert_almost_equal(self.res1[0], self.teststat, DECIMAL_5)

    def test_pvalue(self):
        assert_almost_equal(self.res1[1], self.pvalue, DECIMAL_5)

    def test_critvalues(self):
        critvalues = [self.res1[4][lev] for lev in self.levels]
        assert_almost_equal(critvalues, self.critvalues, DECIMAL_2)


class TestADFConstant(CheckADF):
    """
    Dickey-Fuller test for unit root
    """

    @classmethod
    def setup_class(cls):
        cls.res1 = adfuller(cls.x, regression="c", autolag=None, maxlag=4)
        cls.teststat = 0.97505319
        cls.pvalue = 0.99399563
        cls.critvalues = [-3.476, -2.883, -2.573]


class TestADFConstantTrend(CheckADF):
    """"""

    @classmethod
    def setup_class(cls):
        cls.res1 = adfuller(cls.x, regression="ct", autolag=None, maxlag=4)
        cls.teststat = -1.8566374
        cls.pvalue = 0.67682968
        cls.critvalues = [-4.007, -3.437, -3.137]


# FIXME: do not leave commented-out
# class TestADFConstantTrendSquared(CheckADF):
#    """
#    """
#    pass
# TODO: get test values from R?


class TestADFNoConstant(CheckADF):
    """"""

    @classmethod
    def setup_class(cls):
        cls.res1 = adfuller(cls.x, regression="n", autolag=None, maxlag=4)
        cls.teststat = 3.5227498

        cls.pvalue = 0.99999
        # Stata does not return a p-value for noconstant.
        # Tau^max in MacKinnon (1994) is missing, so it is
        # assumed that its right-tail is well-behaved

        cls.critvalues = [-2.587, -1.950, -1.617]


# No Unit Root


class TestADFConstant2(CheckADF):
    @classmethod
    def setup_class(cls):
        cls.res1 = adfuller(cls.y, regression="c", autolag=None, maxlag=1)
        cls.teststat = -4.3346988
        cls.pvalue = 0.00038661
        cls.critvalues = [-3.476, -2.883, -2.573]


class TestADFConstantTrend2(CheckADF):
    @classmethod
    def setup_class(cls):
        cls.res1 = adfuller(cls.y, regression="ct", autolag=None, maxlag=1)
        cls.teststat = -4.425093
        cls.pvalue = 0.00199633
        cls.critvalues = [-4.006, -3.437, -3.137]


class TestADFNoConstant2(CheckADF):
    @classmethod
    def setup_class(cls):
        cls.res1 = adfuller(cls.y, regression="n", autolag=None, maxlag=1)
        cls.teststat = -2.4511596
        cls.pvalue = 0.013747
        # Stata does not return a p-value for noconstant
        # this value is just taken from our results
        cls.critvalues = [-2.587, -1.950, -1.617]
        _, _1, _2, cls.store = adfuller(
            cls.y, regression="n", autolag=None, maxlag=1, store=True
        )

    def test_store_str(self):
        assert_equal(self.store.__str__(), "Augmented Dickey-Fuller Test Results")


@pytest.mark.parametrize("x", [np.full(8, 5.0)])
def test_adfuller_resid_variance_zero(x):
    with pytest.raises(ValueError):
        adfuller(x)


class CheckCorrGram:
    """
    Set up for ACF, PACF tests.
    """

    data = macrodata.load_pandas()
    x = data.data["realgdp"]
    filename = os.path.join(CURR_DIR, "results", "results_corrgram.csv")
    results = pd.read_csv(filename, delimiter=",")


class TestACF(CheckCorrGram):
    """
    Test Autocorrelation Function
    """

    @classmethod
    def setup_class(cls):
        cls.acf = cls.results["acvar"]
        # cls.acf = np.concatenate(([1.], cls.acf))
        cls.qstat = cls.results["Q1"]
        cls.res1 = acf(cls.x, nlags=40, qstat=True, alpha=0.05, fft=False)
        cls.confint_res = cls.results[["acvar_lb", "acvar_ub"]].values

    def test_acf(self):
        assert_almost_equal(self.res1[0][1:41], self.acf, DECIMAL_8)

    def test_confint(self):
        centered = self.res1[1] - self.res1[1].mean(1)[:, None]
        assert_almost_equal(centered[1:41], self.confint_res, DECIMAL_8)

    def test_qstat(self):
        assert_almost_equal(self.res1[2][:40], self.qstat, DECIMAL_3)
        # 3 decimal places because of stata rounding

    # FIXME: enable/xfail/skip or delete
    # def pvalue(self):
    #    pass
    # NOTE: should not need testing if Q stat is correct


class TestACF_FFT(CheckCorrGram):
    # Test Autocorrelation Function using FFT
    @classmethod
    def setup_class(cls):
        cls.acf = cls.results["acvarfft"]
        cls.qstat = cls.results["Q1"]
        cls.res1 = acf(cls.x, nlags=40, qstat=True, fft=True)

    def test_acf(self):
        assert_almost_equal(self.res1[0][1:], self.acf, DECIMAL_8)

    def test_qstat(self):
        # todo why is res1/qstat 1 short
        assert_almost_equal(self.res1[1], self.qstat, DECIMAL_3)


class TestACFMissing(CheckCorrGram):
    # Test Autocorrelation Function using Missing
    @classmethod
    def setup_class(cls):
        cls.x = np.concatenate((np.array([np.nan]), cls.x))
        cls.acf = cls.results["acvar"]  # drop and conservative
        cls.qstat = cls.results["Q1"]
        cls.res_drop = acf(
            cls.x, nlags=40, qstat=True, alpha=0.05, missing="drop", fft=False
        )
        cls.res_conservative = acf(
            cls.x,
            nlags=40,
            qstat=True,
            alpha=0.05,
            fft=False,
            missing="conservative",
        )
        cls.acf_none = np.empty(40) * np.nan  # lags 1 to 40 inclusive
        cls.qstat_none = np.empty(40) * np.nan
        cls.res_none = acf(
            cls.x, nlags=40, qstat=True, alpha=0.05, missing="none", fft=False
        )

    def test_raise(self):
        with pytest.raises(MissingDataError):
            acf(
                self.x,
                nlags=40,
                qstat=True,
                fft=False,
                alpha=0.05,
                missing="raise",
            )

    def test_acf_none(self):
        assert_almost_equal(self.res_none[0][1:41], self.acf_none, DECIMAL_8)

    def test_acf_drop(self):
        assert_almost_equal(self.res_drop[0][1:41], self.acf, DECIMAL_8)

    def test_acf_conservative(self):
        assert_almost_equal(self.res_conservative[0][1:41], self.acf, DECIMAL_8)

    def test_qstat_none(self):
        # todo why is res1/qstat 1 short
        assert_almost_equal(self.res_none[2], self.qstat_none, DECIMAL_3)


# FIXME: enable/xfail/skip or delete
# how to do this test? the correct q_stat depends on whether nobs=len(x) is
# used when x contains NaNs or whether nobs<len(x) when x contains NaNs
#    def test_qstat_drop(self):
#        assert_almost_equal(self.res_drop[2][:40], self.qstat, DECIMAL_3)


class TestPACF(CheckCorrGram):
    @classmethod
    def setup_class(cls):
        cls.pacfols = cls.results["PACOLS"]
        cls.pacfyw = cls.results["PACYW"]

    def test_ols(self):
        pacfols, confint = pacf(self.x, nlags=40, alpha=0.05, method="ols")
        assert_almost_equal(pacfols[1:], self.pacfols, DECIMAL_6)
        centered = confint - confint.mean(1)[:, None]
        # from edited Stata ado file
        res = [[-0.1375625, 0.1375625]] * 40
        assert_almost_equal(centered[1:41], res, DECIMAL_6)
        # check lag 0
        assert_equal(centered[0], [0.0, 0.0])
        assert_equal(confint[0], [1, 1])
        assert_equal(pacfols[0], 1)

    def test_ols_inefficient(self):
        lag_len = 5
        pacfols = pacf_ols(self.x, nlags=lag_len, efficient=False)
        x = self.x.copy()
        x -= x.mean()
        n = x.shape[0]
        lags = np.zeros((n - 5, 5))
        lead = x[5:]
        direct = np.empty(lag_len + 1)
        direct[0] = 1.0
        for i in range(lag_len):
            lags[:, i] = x[5 - (i + 1) : -(i + 1)]
            direct[i + 1] = lstsq(lags[:, : (i + 1)], lead, rcond=None)[0][-1]
        assert_allclose(pacfols, direct, atol=1e-8)

    def test_yw(self):
        pacfyw = pacf_yw(self.x, nlags=40, method="mle")
        assert_almost_equal(pacfyw[1:], self.pacfyw, DECIMAL_8)

    @pytest.mark.skipif(PYTHON_IMPL_WASM, reason="No fp exception support in WASM")
    def test_yw_singular(self):
        with pytest.warns(SingularMatrixWarning):
            pacf(np.ones(30), nlags=6)

    def test_ld(self):
        pacfyw = pacf_yw(self.x, nlags=40, method="mle")
        pacfld = pacf(self.x, nlags=40, method="ldb")
        assert_almost_equal(pacfyw, pacfld, DECIMAL_8)

        pacfyw = pacf(self.x, nlags=40, method="yw")
        pacfld = pacf(self.x, nlags=40, method="lda")
        assert_almost_equal(pacfyw, pacfld, DECIMAL_8)

    def test_burg(self):
        pacfburg_, _ = pacf_burg(self.x, nlags=40)
        pacfburg = pacf(self.x, nlags=40, method="burg")
        assert_almost_equal(pacfburg_, pacfburg, DECIMAL_8)


class TestCCF:
    """
    Test cross-correlation function
    """

    data = macrodata.load_pandas()
    x = data.data["unemp"].diff().dropna()
    y = data.data["infl"].diff().dropna()
    filename = os.path.join(CURR_DIR, "results", "results_ccf.csv")
    results = pd.read_csv(filename, delimiter=",")
    nlags = 20

    @classmethod
    def setup_class(cls):
        cls.ccf = cls.results["ccf"]
        cls.res1 = ccf(cls.x, cls.y, nlags=cls.nlags, adjusted=False, fft=False)

    def test_ccf(self):
        assert_almost_equal(self.res1, self.ccf, DECIMAL_8)

    def test_confint(self):
        alpha = 0.05
        res2, confint = ccf(
            self.x, self.y, nlags=self.nlags, adjusted=False, fft=False, alpha=alpha
        )
        assert_equal(res2, self.res1)
        assert_almost_equal(res2 - confint[:, 0], confint[:, 1] - res2, DECIMAL_8)
        alpha1 = stats.norm.cdf(confint[:, 1] - res2, scale=1.0 / np.sqrt(len(self.x)))
        assert_almost_equal(alpha1, np.repeat(1 - alpha / 2.0, self.nlags), DECIMAL_8)


class TestBreakvarHeteroskedasticityTest:
    from scipy.stats import chi2, f

    def test_1d_input(self):

        input_residuals = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
        expected_statistic = (4.0**2 + 5.0**2) / (0.0**2 + 1.0**2)
        # ~ F(2, 2), two-sided test
        expected_pvalue = 2 * min(
            self.f.cdf(expected_statistic, 2, 2),
            self.f.sf(expected_statistic, 2, 2),
        )
        actual_statistic, actual_pvalue = breakvar_heteroskedasticity_test(
            input_residuals
        )

        assert actual_statistic == expected_statistic
        assert actual_pvalue == expected_pvalue

    def test_2d_input_with_missing_values(self):

        input_residuals = np.array(
            [
                [0.0, 0.0, np.nan],
                [1.0, np.nan, 1.0],
                [2.0, 2.0, np.nan],
                [3.0, 3.0, 3.0],
                [4.0, 4.0, 4.0],
                [5.0, 5.0, 5.0],
                [6.0, 6.0, 6.0],
                [7.0, 7.0, 7.0],
                [8.0, 8.0, 8.0],
            ]
        )
        expected_statistic = np.array(
            [
                (8.0**2 + 7.0**2 + 6.0**2) / (0.0**2 + 1.0**2 + 2.0**2),
                (8.0**2 + 7.0**2 + 6.0**2) / (0.0**2 + 2.0**2),
                np.nan,
            ]
        )
        expected_pvalue = np.array(
            [
                2
                * min(
                    self.f.cdf(expected_statistic[0], 3, 3),
                    self.f.sf(expected_statistic[0], 3, 3),
                ),
                2
                * min(
                    self.f.cdf(expected_statistic[1], 3, 2),
                    self.f.sf(expected_statistic[1], 3, 2),
                ),
                np.nan,
            ]
        )
        actual_statistic, actual_pvalue = breakvar_heteroskedasticity_test(
            input_residuals
        )

        assert_equal(actual_statistic, expected_statistic)
        assert_equal(actual_pvalue, expected_pvalue)

    @pytest.mark.parametrize(
        "subset_length,expected_statistic,expected_pvalue",
        [
            (2, 41, 2 * min(f.cdf(41, 2, 2), f.sf(41, 2, 2))),
            (0.5, 10, 2 * min(f.cdf(10, 3, 3), f.sf(10, 3, 3))),
        ],
    )
    def test_subset_length(self, subset_length, expected_statistic, expected_pvalue):

        input_residuals = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
        actual_statistic, actual_pvalue = breakvar_heteroskedasticity_test(
            input_residuals,
            subset_length=subset_length,
        )

        assert actual_statistic == expected_statistic
        assert actual_pvalue == expected_pvalue

    @pytest.mark.parametrize(
        "alternative,expected_statistic,expected_pvalue",
        [
            ("two-sided", 41, 2 * min(f.cdf(41, 2, 2), f.sf(41, 2, 2))),
            ("decreasing", 1 / 41, f.sf(1 / 41, 2, 2)),
            ("increasing", 41, f.sf(41, 2, 2)),
        ],
    )
    def test_alternative(self, alternative, expected_statistic, expected_pvalue):

        input_residuals = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
        actual_statistic, actual_pvalue = breakvar_heteroskedasticity_test(
            input_residuals,
            alternative=alternative,
        )
        assert actual_statistic == expected_statistic
        assert actual_pvalue == expected_pvalue

    def test_use_chi2(self):

        input_residuals = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
        expected_statistic = (4.0**2 + 5.0**2) / (0.0**2 + 1.0**2)
        expected_pvalue = 2 * min(
            self.chi2.cdf(2 * expected_statistic, 2),
            self.chi2.sf(2 * expected_statistic, 2),
        )
        actual_statistic, actual_pvalue = breakvar_heteroskedasticity_test(
            input_residuals,
            use_f=False,
        )
        assert actual_statistic == expected_statistic
        assert actual_pvalue == expected_pvalue


class CheckCoint:
    """
    Test Cointegration Test Results for 2-variable system

    Test values taken from Stata
    """

    levels = ["1%", "5%", "10%"]
    data = macrodata.load_pandas()
    y1 = data.data["realcons"].values
    y2 = data.data["realgdp"].values

    def test_tstat(self):
        assert_almost_equal(self.coint_t, self.teststat, DECIMAL_4)


# this does not produce the old results anymore
class TestCoint_t(CheckCoint):
    """
    Get AR(1) parameter on residuals
    """

    @classmethod
    def setup_class(cls):
        # cls.coint_t = coint(cls.y1, cls.y2, trend="c")[0]
        cls.coint_t = coint(cls.y1, cls.y2, trend="c", maxlag=0, autolag=None)[0]
        cls.teststat = -1.8208817
        cls.teststat = -1.830170986148


def test_coint():
    nobs = 200
    scale_e = 1
    const = [1, 0, 0.5, 0]
    np.random.seed(123)
    unit = np.random.randn(nobs).cumsum()
    y = scale_e * np.random.randn(nobs, 4)
    y[:, :2] += unit[:, None]
    y += const
    y = np.round(y, 4)

    # FIXME: enable/xfail/skip or delete
    for trend in []:  # ['c', 'ct', 'ctt', 'n']:
        print("\n", trend)
        print(coint(y[:, 0], y[:, 1], trend=trend, maxlag=4, autolag=None))
        print(coint(y[:, 0], y[:, 1:3], trend=trend, maxlag=4, autolag=None))
        print(coint(y[:, 0], y[:, 2:], trend=trend, maxlag=4, autolag=None))
        print(coint(y[:, 0], y[:, 1:], trend=trend, maxlag=4, autolag=None))

    # results from Stata egranger
    res_egranger = {}
    # trend = 'ct'
    res = res_egranger["ct"] = {}
    res[0] = [
        -5.615251442239,
        -4.406102369132,
        -3.82866685109,
        -3.532082997903,
    ]
    res[1] = [
        -5.63591313706,
        -4.758609717199,
        -4.179130554708,
        -3.880909696863,
    ]
    res[2] = [
        -2.892029275027,
        -4.758609717199,
        -4.179130554708,
        -3.880909696863,
    ]
    res[3] = [-5.626932544079, -5.08363327039, -4.502469783057, -4.2031051091]

    # trend = 'c'
    res = res_egranger["c"] = {}
    # first critical value res[0][1] has a discrepancy starting at 4th decimal
    res[0] = [
        -5.760696844656,
        -3.952043522638,
        -3.367006313729,
        -3.065831247948,
    ]
    # manually adjusted to have higher precision as in other cases
    res[0][1] = -3.952321293401682
    res[1] = [
        -5.781087068772,
        -4.367111915942,
        -3.783961136005,
        -3.483501524709,
    ]
    res[2] = [
        -2.477444137366,
        -4.367111915942,
        -3.783961136005,
        -3.483501524709,
    ]
    res[3] = [
        -5.778205811661,
        -4.735249216434,
        -4.152738973763,
        -3.852480848968,
    ]

    # trend = 'ctt'
    res = res_egranger["ctt"] = {}
    res[0] = [
        -5.644431269946,
        -4.796038299708,
        -4.221469431008,
        -3.926472577178,
    ]
    res[1] = [-5.665691609506, -5.111158174219, -4.53317278104, -4.23601008516]
    res[2] = [-3.161462374828, -5.111158174219, -4.53317278104, -4.23601008516]
    res[3] = [
        -5.657904558563,
        -5.406880189412,
        -4.826111619543,
        -4.527090164875,
    ]

    # The following for 'n' are only regression test numbers
    # trend = 'n' not allowed in egranger
    # trend = 'n'
    res = res_egranger["n"] = {}
    nan = np.nan  # shortcut for table
    res[0] = [-3.7146175989071137, nan, nan, nan]
    res[1] = [-3.8199323012888384, nan, nan, nan]
    res[2] = [-1.6865000791270679, nan, nan, nan]
    res[3] = [-3.7991270451873675, nan, nan, nan]

    for trend in ["c", "ct", "ctt", "n"]:
        res1 = {}
        res1[0] = coint(y[:, 0], y[:, 1], trend=trend, maxlag=4, autolag=None)
        res1[1] = coint(y[:, 0], y[:, 1:3], trend=trend, maxlag=4, autolag=None)
        res1[2] = coint(y[:, 0], y[:, 2:], trend=trend, maxlag=4, autolag=None)
        res1[3] = coint(y[:, 0], y[:, 1:], trend=trend, maxlag=4, autolag=None)

        for i in range(4):
            res = res_egranger[trend]

            assert_allclose(res1[i][0], res[i][0], rtol=1e-11)
            r2 = res[i][1:]
            r1 = res1[i][2]
            assert_allclose(r1, r2, rtol=0, atol=6e-7)

    # use default autolag #4490
    res1_0 = coint(y[:, 0], y[:, 1], trend="ct", maxlag=4)
    assert_allclose(res1_0[2], res_egranger["ct"][0][1:], rtol=0, atol=6e-7)
    # the following is just a regression test
    assert_allclose(
        res1_0[:2],
        [-13.992946638547112, 2.270898990540678e-27],
        rtol=1e-10,
        atol=1e-27,
    )


def test_coint_identical_series():
    nobs = 200
    scale_e = 1
    np.random.seed(123)
    y = scale_e * np.random.randn(nobs)
    warnings.simplefilter("always", CollinearityWarning)
    with pytest.warns(CollinearityWarning):
        c = coint(y, y, trend="c", maxlag=0, autolag=None)
    assert_equal(c[1], 0.0)
    assert (np.isneginf(c[0]))


def test_coint_perfect_collinearity():
    # test uses nearly perfect collinearity
    nobs = 200
    scale_e = 1
    np.random.seed(123)
    x = scale_e * np.random.randn(nobs, 2)
    y = 1 + x.sum(axis=1) + 1e-7 * np.random.randn(nobs)
    warnings.simplefilter("always", CollinearityWarning)
    with warnings.catch_warnings(record=True):
        c = coint(y, x, trend="c", maxlag=0, autolag=None)
    assert_equal(c[1], 0.0)
    assert (np.isneginf(c[0]))


class TestGrangerCausality:
    def test_grangercausality(self):
        # some example data
        mdata = macrodata.load_pandas().data
        mdata = mdata[["realgdp", "realcons"]].values
        data = mdata.astype(float)
        data = np.diff(np.log(data), axis=0)

        # R: lmtest:grangertest
        r_result = [0.243097, 0.7844328, 195, 2]  # f_test
        with pytest.warns(FutureWarning, match="verbose is"):
            gr = grangercausalitytests(data[:, 1::-1], 2, verbose=False)
        assert_almost_equal(r_result, gr[2][0]["ssr_ftest"], decimal=7)
        assert_almost_equal(gr[2][0]["params_ftest"], gr[2][0]["ssr_ftest"], decimal=7)

    def test_grangercausality_single(self):
        mdata = macrodata.load_pandas().data
        mdata = mdata[["realgdp", "realcons"]].values
        data = mdata.astype(float)
        data = np.diff(np.log(data), axis=0)
        with pytest.warns(FutureWarning, match="verbose is"):
            gr = grangercausalitytests(data[:, 1::-1], 2, verbose=False)
        with pytest.warns(FutureWarning, match="verbose is"):
            gr2 = grangercausalitytests(data[:, 1::-1], [2], verbose=False)
        assert 1 in gr
        assert 1 not in gr2
        assert_almost_equal(gr[2][0]["ssr_ftest"], gr2[2][0]["ssr_ftest"], decimal=7)
        assert_almost_equal(gr[2][0]["params_ftest"], gr2[2][0]["ssr_ftest"], decimal=7)

    def test_granger_fails_on_nobs_check(self, reset_randomstate):
        # Test that if maxlag is too large, Granger Test raises a clear error.
        x = np.random.rand(10, 2)
        with pytest.warns(FutureWarning, match="verbose is"):
            grangercausalitytests(x, 2, verbose=False)  # This should pass.
        with pytest.raises(ValueError):
            with pytest.warns(FutureWarning, match="verbose is"):
                grangercausalitytests(x, 3, verbose=False)

    def test_granger_fails_on_finite_check(self, reset_randomstate):
        x = np.random.rand(1000, 2)
        x[500, 0] = np.nan
        x[750, 1] = np.inf
        with pytest.raises(ValueError, match="x contains NaN"):
            grangercausalitytests(x, 2)

    def test_granger_fails_on_zero_lag(self, reset_randomstate):
        x = np.random.rand(1000, 2)
        with pytest.raises(
            ValueError,
            match="maxlag must be a non-empty list containing only positive integers",
        ):
            grangercausalitytests(x, [0, 1, 2])


class TestKPSS:
    """
    R-code
    ------
    library(tseries)
    kpss.stat(x, "Level")
    kpss.stat(x, "Trend")

    In this context, x is the vector containing the
    macrodata['realgdp'] series.
    """

    def setup_method(self):
        self.data = macrodata.load_pandas()
        self.x = self.data.data["realgdp"].values

    def test_fail_nonvector_input(self, reset_randomstate):
        # should be fine
        with pytest.warns(InterpolationWarning):
            kpss(self.x, nlags="legacy")

        x = np.random.rand(20, 2)
        with pytest.raises(ValueError):
            kpss(x)

    def test_fail_unclear_hypothesis(self):
        # these should be fine,
        with pytest.warns(InterpolationWarning):
            kpss(self.x, "c", nlags="legacy")
        with pytest.warns(InterpolationWarning):
            kpss(self.x, "C", nlags="legacy")
        with pytest.warns(InterpolationWarning):
            kpss(self.x, "ct", nlags="legacy")
        with pytest.warns(InterpolationWarning):
            kpss(self.x, "CT", nlags="legacy")

        with pytest.raises(ValueError):
            kpss(self.x, "unclear hypothesis", nlags="legacy")

    def test_teststat(self):
        with pytest.warns(InterpolationWarning):
            kpss_stat, _, _, _ = kpss(self.x, "c", 3)
        assert_almost_equal(kpss_stat, 5.0169, DECIMAL_3)

        with pytest.warns(InterpolationWarning):
            kpss_stat, _, _, _ = kpss(self.x, "ct", 3)
        assert_almost_equal(kpss_stat, 1.1828, DECIMAL_3)

    def test_pval(self):
        with pytest.warns(InterpolationWarning):
            _, pval, _, _ = kpss(self.x, "c", 3)
        assert_equal(pval, 0.01)

        with pytest.warns(InterpolationWarning):
            _, pval, _, _ = kpss(self.x, "ct", 3)
        assert_equal(pval, 0.01)

    def test_store(self):
        with pytest.warns(InterpolationWarning):
            _, _, _, store = kpss(self.x, "c", 3, True)

        # assert attributes, and make sure they're correct
        assert_equal(store.nobs, len(self.x))
        assert_equal(store.lags, 3)

    # test autolag function _kpss_autolag against SAS 9.3
    def test_lags(self):
        # real GDP from macrodata data set
        with pytest.warns(InterpolationWarning):
            res = kpss(self.x, "c", nlags="auto")
        assert_equal(res[2], 9)
        # real interest rates from macrodata data set
        res = kpss(sunspots.load().data["SUNACTIVITY"], "c", nlags="auto")
        assert_equal(res[2], 7)
        # volumes from nile data set
        with pytest.warns(InterpolationWarning):
            res = kpss(nile.load().data["volume"], "c", nlags="auto")
        assert_equal(res[2], 5)
        # log-coinsurance from randhie data set
        with pytest.warns(InterpolationWarning):
            res = kpss(randhie.load().data["lncoins"], "ct", nlags="auto")
        assert_equal(res[2], 75)
        # in-vehicle time from modechoice data set
        with pytest.warns(InterpolationWarning):
            res = kpss(modechoice.load().data["invt"], "ct", nlags="auto")
        assert_equal(res[2], 18)

    def test_kpss_fails_on_nobs_check(self):
        # Test that if lags exceeds number of observations KPSS raises a
        # clear error
        # GH5925
        nobs = len(self.x)
        msg = r"lags \({}\) must be < number of observations \({}\)".format(nobs, nobs)
        with pytest.raises(ValueError, match=msg):
            kpss(self.x, "c", nlags=nobs)

    def test_kpss_autolags_does_not_assign_lags_equal_to_nobs(self):
        # Test that if *autolags* exceeds number of observations, we set
        # suitable lags
        # GH5925
        base = np.array([0, 0, 0, 0, 0, 1, 1.0])
        data_which_breaks_autolag = np.r_[np.tile(base, 297 // 7), [0, 0, 0]]
        kpss(data_which_breaks_autolag, nlags="auto")

    def test_legacy_lags(self):
        # Test legacy lags are the same
        with pytest.warns(InterpolationWarning):
            res = kpss(self.x, "c", nlags="legacy")
        assert_equal(res[2], 15)

    def test_unknown_lags(self):
        # Test legacy lags are the same
        with pytest.raises(ValueError):
            kpss(self.x, "c", nlags="unknown")

    def test_none(self):
        with pytest.warns(FutureWarning):
            kpss(self.x, nlags=None)


class TestRUR:
    """
    Simple implementation
    ------
    Since an R implementation of the test cannot be found, the method is tested against
    a simple implementation using a for loop.
    In this context, x is the vector containing the
    macrodata['realgdp'] series.
    """

    def setup_method(self):
        self.data = macrodata.load_pandas()
        self.x = self.data.data["realgdp"].values

    # To be removed when range unit test gets an R implementation
    def simple_rur(self, x, store=False):
        x = array_like(x, "x")
        store = bool_like(store, "store")

        nobs = x.shape[0]

        # if m is not one, n != m * n
        if nobs != x.size:
            raise ValueError(f"x of shape {x.shape} not understood")

        # Table from [1] has been replicated using 200,000 samples
        # Critical values for new n_obs values have been identified
        pvals = [0.01, 0.025, 0.05, 0.10, 0.90, 0.95]
        n = np.array([25, 50, 100, 150, 200, 250, 500, 1000, 2000, 3000, 4000, 5000])
        crit = np.array(
            [
                [0.6626, 0.8126, 0.9192, 1.0712, 2.4863, 2.7312],
                [0.7977, 0.9274, 1.0478, 1.1964, 2.6821, 2.9613],
                [0.907, 1.0243, 1.1412, 1.2888, 2.8317, 3.1393],
                [0.9543, 1.0768, 1.1869, 1.3294, 2.8915, 3.2049],
                [0.9833, 1.0984, 1.2101, 1.3494, 2.9308, 3.2482],
                [0.9982, 1.1137, 1.2242, 1.3632, 2.9571, 3.2482],
                [1.0494, 1.1643, 1.2712, 1.4076, 3.0207, 3.3584],
                [1.0846, 1.1959, 1.2988, 1.4344, 3.0653, 3.4073],
                [1.1121, 1.2200, 1.3230, 1.4556, 3.0948, 3.4439],
                [1.1204, 1.2295, 1.3318, 1.4656, 3.1054, 3.4632],
                [1.1309, 1.2347, 1.3318, 1.4693, 3.1165, 3.4717],
                [1.1377, 1.2402, 1.3408, 1.4729, 3.1252, 3.4807],
            ]
        )

        # Interpolation for nobs
        inter_crit = np.zeros((1, crit.shape[1]))
        for i in range(crit.shape[1]):
            f = interp1d(n, crit[:, i])
            inter_crit[0, i] = f(nobs)

        # Calculate RUR stat
        count = 0

        max_p = x[0]
        min_p = x[0]

        for v in x[1:]:
            if v > max_p:
                max_p = v
                count = count + 1
            if v < min_p:
                min_p = v
                count = count + 1

        rur_stat = count / np.sqrt(len(x))

        k = len(pvals) - 1
        for i in range(len(pvals) - 1, -1, -1):
            if rur_stat < inter_crit[0, i]:
                k = i
            else:
                break

        p_value = pvals[k]

        warn_msg = """\
        The test statistic is outside of the range of p-values available in the
        look-up table. The actual p-value is {direction} than the p-value returned.
        """
        direction = ""
        if p_value == pvals[-1]:
            direction = "smaller"
        elif p_value == pvals[0]:
            direction = "larger"

        if direction:
            warnings.warn(
                warn_msg.format(direction=direction), InterpolationWarning, stacklevel=2
            )

        crit_dict = {
            "10%": inter_crit[0, 3],
            "5%": inter_crit[0, 2],
            "2.5%": inter_crit[0, 1],
            "1%": inter_crit[0, 0],
        }

        if store:
            from statsmodels.stats.diagnostic import ResultsStore

            rstore = ResultsStore()
            rstore.nobs = nobs

            rstore.H0 = "The series is not stationary"
            rstore.HA = "The series is stationary"

            return rur_stat, p_value, crit_dict, rstore
        else:
            return rur_stat, p_value, crit_dict

    def test_fail_nonvector_input(self, reset_randomstate):
        with pytest.warns(InterpolationWarning):
            range_unit_root_test(self.x)

        x = np.random.rand(20, 2)
        with pytest.raises(ValueError):
            range_unit_root_test(x)

    def test_teststat(self):
        with pytest.warns(InterpolationWarning):
            rur_stat, _, _ = range_unit_root_test(self.x)
        simple_rur_stat, _, _ = self.simple_rur(self.x)
        assert_almost_equal(rur_stat, simple_rur_stat, DECIMAL_3)

    def test_pval(self):
        with pytest.warns(InterpolationWarning):
            _, pval, _ = range_unit_root_test(self.x)
        _, simple_pval, _ = self.simple_rur(self.x)
        assert_equal(pval, simple_pval)

    def test_store(self):
        with pytest.warns(InterpolationWarning):
            _, _, _, store = range_unit_root_test(self.x, True)

        # assert attributes, and make sure they're correct
        assert_equal(store.nobs, len(self.x))


def test_pandasacovf():
    s = Series(lrange(1, 11))
    assert_almost_equal(acovf(s, fft=False), acovf(s.values, fft=False))


def test_acovf2d(reset_randomstate):
    dta = sunspots.load_pandas().data
    dta.index = date_range(start="1700", end="2009", freq=YEAR_END)[:309]
    del dta["YEAR"]
    res = acovf(dta, fft=False)
