"""
Test functions for models.regression
"""

from statsmodels.compat.python import lrange

# TODO: Test for LM
from statsmodels.compat.scipy import SP_LT_116

import warnings

import numpy as np
from numpy.testing import (
    assert_,
    assert_allclose,
    assert_almost_equal,
    assert_equal,
)
import pandas as pd
import pytest
from scipy.linalg import toeplitz
from scipy.stats import t as student_t

from statsmodels.datasets import longley
from statsmodels.formula._manager import FormulaManager
from statsmodels.regression.linear_model import (
    GLS,
    OLS,
    WLS,
    burg,
    yule_walker,
)
from statsmodels.tools.tools import add_constant

DECIMAL_4 = 4
DECIMAL_3 = 3
DECIMAL_2 = 2
DECIMAL_1 = 1
DECIMAL_7 = 7
DECIMAL_0 = 0

try:
    import cvxopt  # noqa:F401

    has_cvxopt = True
except ImportError:
    has_cvxopt = False


class CheckRegressionResults:
    """
    res2 contains results from Rmodelwrap or were obtained from a statistical
    packages such as R, Stata, or SAS and were written to model_results
    """

    decimal_params = DECIMAL_4

    def test_params(self):
        assert_almost_equal(self.res1.params, self.res2.params, self.decimal_params)

    decimal_standarderrors = DECIMAL_4

    def test_standarderrors(self):
        assert_allclose(self.res1.bse, self.res2.bse, self.decimal_standarderrors)

    decimal_confidenceintervals = DECIMAL_4

    def test_confidenceintervals(self):
        # NOTE: stata rounds residuals (at least) to sig digits so approx_equal
        conf1 = self.res1.conf_int()
        conf2 = self.res2.conf_int()
        for i in range(len(conf1)):
            assert_allclose(
                conf1[i][0],
                conf2[i][0],
                rtol=10**-self.decimal_confidenceintervals,
            )
            assert_allclose(
                conf1[i][1],
                conf2[i][1],
                rtol=10**-self.decimal_confidenceintervals,
            )

    decimal_conf_int_subset = DECIMAL_4

    def test_conf_int_subset(self):
        if len(self.res1.params) > 1:
            with pytest.warns(FutureWarning, match="cols is"):
                ci1 = self.res1.conf_int(cols=(1, 2))
            ci2 = self.res1.conf_int()[1:3]
            assert_almost_equal(ci1, ci2, self.decimal_conf_int_subset)
        else:
            pass

    decimal_scale = DECIMAL_4

    def test_scale(self):
        assert_almost_equal(self.res1.scale, self.res2.scale, self.decimal_scale)

    decimal_rsquared = DECIMAL_4

    def test_rsquared(self):
        assert_almost_equal(
            self.res1.rsquared, self.res2.rsquared, self.decimal_rsquared
        )

    decimal_rsquared_adj = DECIMAL_4

    def test_rsquared_adj(self):
        assert_almost_equal(
            self.res1.rsquared_adj,
            self.res2.rsquared_adj,
            self.decimal_rsquared_adj,
        )

    def test_degrees(self):
        assert_equal(self.res1.model.df_model, self.res2.df_model)
        assert_equal(self.res1.model.df_resid, self.res2.df_resid)

    decimal_ess = DECIMAL_4

    def test_ess(self):
        # Explained Sum of Squares
        assert_almost_equal(self.res1.ess, self.res2.ess, self.decimal_ess)

    decimal_ssr = DECIMAL_4

    def test_sumof_squaredresids(self):
        assert_almost_equal(self.res1.ssr, self.res2.ssr, self.decimal_ssr)

    decimal_mse_resid = DECIMAL_4

    def test_mse_resid(self):
        # Mean squared error of residuals
        assert_almost_equal(
            self.res1.mse_model, self.res2.mse_model, self.decimal_mse_resid
        )

    decimal_mse_model = DECIMAL_4

    def test_mse_model(self):
        assert_almost_equal(
            self.res1.mse_resid, self.res2.mse_resid, self.decimal_mse_model
        )

    decimal_mse_total = DECIMAL_4

    def test_mse_total(self):
        assert_almost_equal(
            self.res1.mse_total,
            self.res2.mse_total,
            self.decimal_mse_total,
            err_msg="Test class %s" % self,
        )

    decimal_fvalue = DECIMAL_4

    def test_fvalue(self):
        # did not change this, not sure it should complain -inf not equal -inf
        # if not (np.isinf(self.res1.fvalue) and np.isinf(self.res2.fvalue)):
        assert_almost_equal(self.res1.fvalue, self.res2.fvalue, self.decimal_fvalue)

    decimal_loglike = DECIMAL_4

    def test_loglike(self):
        assert_almost_equal(self.res1.llf, self.res2.llf, self.decimal_loglike)

    decimal_aic = DECIMAL_4

    def test_aic(self):
        assert_almost_equal(self.res1.aic, self.res2.aic, self.decimal_aic)
        # the following just checks the definition
        aicc1 = self.res1.info_criteria("aicc")
        k = self.res1.df_model + self.res1.model.k_constant
        nobs = self.res1.model.nobs
        aicc2 = self.res1.aic + 2 * (k**2 + k) / (nobs - k - 1)
        assert_allclose(aicc1, aicc2, rtol=1e-10)
        hqic1 = self.res1.info_criteria("hqic")
        hqic2 = (self.res1.aic - 2 * k) + 2 * np.log(np.log(nobs)) * k
        assert_allclose(hqic1, hqic2, rtol=1e-10)

    decimal_bic = DECIMAL_4

    def test_bic(self):
        assert_almost_equal(self.res1.bic, self.res2.bic, self.decimal_bic)

    decimal_pvalues = DECIMAL_4

    def test_pvalues(self):
        assert_almost_equal(self.res1.pvalues, self.res2.pvalues, self.decimal_pvalues)

    decimal_wresid = DECIMAL_4

    def test_wresid(self):
        assert_almost_equal(self.res1.wresid, self.res2.wresid, self.decimal_wresid)

    decimal_resids = DECIMAL_4

    def test_resids(self):
        assert_almost_equal(self.res1.resid, self.res2.resid, self.decimal_resids)

    decimal_norm_resids = DECIMAL_4

    def test_norm_resids(self):
        assert_almost_equal(
            self.res1.resid_pearson,
            self.res2.resid_pearson,
            self.decimal_norm_resids,
        )


# TODO: test fittedvalues and what else?


class TestOLS(CheckRegressionResults):
    @classmethod
    def setup_class(cls):
        from .results.results_regression import Longley

        data = longley.load()
        endog = np.asarray(data.endog)
        exog = np.asarray(data.exog)
        exog = add_constant(exog, prepend=False)
        res1 = OLS(endog, exog).fit()
        res2 = Longley()
        res2.wresid = res1.wresid  # workaround hack
        cls.res1 = res1
        cls.res2 = res2

        res_qr = OLS(endog, exog).fit(method="qr")

        model_qr = OLS(endog, exog)
        Q, R = np.linalg.qr(exog)
        model_qr.exog_Q, model_qr.exog_R = Q, R
        model_qr.normalized_cov_params = np.linalg.inv(np.dot(R.T, R))
        model_qr.rank = np.linalg.matrix_rank(R)
        res_qr2 = model_qr.fit(method="qr")

        cls.res_qr = res_qr
        cls.res_qr_manual = res_qr2

    def test_eigenvalues(self):
        eigenval_perc_diff = self.res_qr.eigenvals - self.res_qr_manual.eigenvals
        eigenval_perc_diff /= self.res_qr.eigenvals
        zeros = np.zeros_like(eigenval_perc_diff)
        assert_almost_equal(eigenval_perc_diff, zeros, DECIMAL_7)

    # Robust error tests.  Compare values computed with SAS
    def test_HC0_errors(self):
        # They are split up because the copied results do not have any
        # DECIMAL_4 places for the last place.
        assert_almost_equal(self.res1.HC0_se[:-1], self.res2.HC0_se[:-1], DECIMAL_4)
        assert_allclose(np.round(self.res1.HC0_se[-1]), self.res2.HC0_se[-1])

    def test_HC1_errors(self):
        assert_almost_equal(self.res1.HC1_se[:-1], self.res2.HC1_se[:-1], DECIMAL_4)
        # Note: tolerance is tight; rtol=3e-7 fails while 4e-7 passes
        assert_allclose(self.res1.HC1_se[-1], self.res2.HC1_se[-1], rtol=4e-7)

    def test_HC2_errors(self):
        assert_almost_equal(self.res1.HC2_se[:-1], self.res2.HC2_se[:-1], DECIMAL_4)
        # Note: tolerance is tight; rtol=4e-7 fails while 5e-7 passes
        assert_allclose(self.res1.HC2_se[-1], self.res2.HC2_se[-1], rtol=5e-7)

    def test_HC3_errors(self):
        assert_almost_equal(self.res1.HC3_se[:-1], self.res2.HC3_se[:-1], DECIMAL_4)
        # Note: tolerance is tight; rtol=1e-7 fails while 1.5e-7 passes
        assert_allclose(self.res1.HC3_se[-1], self.res2.HC3_se[-1], rtol=1.5e-7)

    def test_qr_params(self):
        assert_almost_equal(self.res1.params, self.res_qr.params, 6)

    def test_qr_normalized_cov_params(self):
        # todo: need assert_close
        assert_almost_equal(
            np.ones_like(self.res1.normalized_cov_params),
            self.res1.normalized_cov_params / self.res_qr.normalized_cov_params,
            5,
        )

    def test_missing(self):
        data = longley.load()
        data.exog = add_constant(data.exog, prepend=False)
        data.endog[[3, 7, 14]] = np.nan
        mod = OLS(data.endog, data.exog, missing="drop")
        assert_equal(mod.endog.shape[0], 13)
        assert_equal(mod.exog.shape[0], 13)

    def test_rsquared_adj_overfit(self):
        # Test that if df_resid = 0, rsquared_adj = 0.
        # This is a regression test for user issue:
        # https://github.com/statsmodels/statsmodels/issues/868
        with warnings.catch_warnings(record=True):
            x = np.random.randn(5)
            y = np.random.randn(5, 6)
            results = OLS(x, y).fit()
            rsquared_adj = results.rsquared_adj
            assert_equal(rsquared_adj, np.nan)

    def test_qr_alternatives(self):
        assert_allclose(self.res_qr.params, self.res_qr_manual.params, rtol=5e-12)

    def test_norm_resid(self):
        resid = self.res1.wresid
        norm_resid = resid / np.sqrt(np.sum(resid**2.0) / self.res1.df_resid)
        model_norm_resid = self.res1.resid_pearson
        assert_almost_equal(model_norm_resid, norm_resid, DECIMAL_7)

    def test_summary_slim(self):
        # check that slim summary is smaller, does not verify content
        with warnings.catch_warnings():
            msg = "kurtosistest only valid for n>=20"
            warnings.filterwarnings("ignore", message=msg, category=UserWarning)

            summ = self.res1.summary(slim=True)
        assert len(summ.tables) == 2
        assert len(str(summ)) < 6700

    def test_norm_resid_zero_variance(self):
        with warnings.catch_warnings(record=True):
            y = self.res1.model.endog
            res = OLS(y, y).fit()
            assert_allclose(res.scale, 0, atol=1e-20)
            assert_allclose(res.wresid, res.resid_pearson, atol=5e-11)


class TestRTO(CheckRegressionResults):
    @classmethod
    def setup_class(cls):
        from .results.results_regression import LongleyRTO

        data = longley.load()
        endog = np.asarray(data.endog)
        exog = np.asarray(data.exog)
        res1 = OLS(endog, exog).fit()
        res2 = LongleyRTO()
        res2.wresid = res1.wresid  # workaround hack
        cls.res1 = res1
        cls.res2 = res2

        res_qr = OLS(endog, exog).fit(method="qr")
        cls.res_qr = res_qr


class TestFtest:
    """
    Tests f_test vs. RegressionResults
    """

    @classmethod
    def setup_class(cls):
        data = longley.load()
        data.exog = add_constant(data.exog, prepend=False)
        cls.res1 = OLS(data.endog, data.exog).fit()
        R = np.identity(7)[:-1, :]
        cls.Ftest = cls.res1.f_test(R)

    def test_F(self):
        assert_almost_equal(self.Ftest.fvalue, self.res1.fvalue, DECIMAL_4)

    def test_p(self):
        assert_almost_equal(self.Ftest.pvalue, self.res1.f_pvalue, DECIMAL_4)

    def test_Df_denom(self):
        assert_equal(self.Ftest.df_denom, self.res1.model.df_resid)

    def test_Df_num(self):
        assert_equal(self.Ftest.df_num, 6)


class TestFTest2:
    """
    A joint test that the coefficient on
    GNP = the coefficient on UNEMP  and that the coefficient on
    POP = the coefficient on YEAR for the Longley dataset.

    Ftest1 is from statsmodels.  Results are from Rpy using R's car library.
    """

    @classmethod
    def setup_class(cls):
        data = longley.load()
        columns = [f"x{i}" for i in range(1, data.exog.shape[1] + 1)]
        data.exog.columns = columns
        data.exog = add_constant(data.exog, prepend=False)
        res1 = OLS(data.endog, data.exog).fit()
        R2 = [[0, 1, -1, 0, 0, 0, 0], [0, 0, 0, 0, 1, -1, 0]]
        cls.Ftest1 = res1.f_test(R2)
        hyp = "x2 = x3, x5 = x6"
        cls.NewFtest1 = res1.f_test(hyp)

    def test_new_ftest(self):
        assert_equal(self.NewFtest1.fvalue, self.Ftest1.fvalue)

    def test_fvalue(self):
        assert_almost_equal(self.Ftest1.fvalue, 9.7404618732968196, DECIMAL_4)

    def test_pvalue(self):
        assert_almost_equal(self.Ftest1.pvalue, 0.0056052885317493459, DECIMAL_4)

    def test_df_denom(self):
        assert_equal(self.Ftest1.df_denom, 9)

    def test_df_num(self):
        assert_equal(self.Ftest1.df_num, 2)


class TestFtestQ:
    """
    A joint hypothesis test that Rb = q.  Coefficient tests are essentially
    made up.  Test values taken from Stata.
    """

    @classmethod
    def setup_class(cls):
        data = longley.load()
        data.exog = add_constant(data.exog, prepend=False)
        res1 = OLS(data.endog, data.exog).fit()
        R = np.array(
            [
                [0, 1, 1, 0, 0, 0, 0],
                [0, 1, 0, 1, 0, 0, 0],
                [0, 1, 0, 0, 0, 0, 0],
                [0, 0, 0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0, 1, 0],
            ]
        )
        q = np.array([0, 0, 0, 1, 0])
        cls.Ftest1 = res1.f_test((R, q))

    def test_fvalue(self):
        assert_almost_equal(self.Ftest1.fvalue, 70.115557, 5)

    def test_pvalue(self):
        assert_almost_equal(self.Ftest1.pvalue, 6.229e-07, 10)

    def test_df_denom(self):
        assert_equal(self.Ftest1.df_denom, 9)

    def test_df_num(self):
        assert_equal(self.Ftest1.df_num, 5)


class TestTtest:
    """
    Test individual t-tests.  Ie., are the coefficients significantly
    different than zero.
    """

    @classmethod
    def setup_class(cls):
        data = longley.load()
        columns = [f"x{i}" for i in range(1, data.exog.shape[1] + 1)]
        data.exog.columns = columns
        data.exog = add_constant(data.exog, prepend=False)
        cls.res1 = OLS(data.endog, data.exog).fit()
        R = np.identity(7)
        cls.Ttest = cls.res1.t_test(R)
        hyp = "x1 = 0, x2 = 0, x3 = 0, x4 = 0, x5 = 0, x6 = 0, const = 0"
        cls.NewTTest = cls.res1.t_test(hyp)

    def test_new_tvalue(self):
        assert_equal(self.NewTTest.tvalue, self.Ttest.tvalue)

    def test_tvalue(self):
        assert_almost_equal(self.Ttest.tvalue, self.res1.tvalues, DECIMAL_4)

    def test_sd(self):
        assert_almost_equal(self.Ttest.sd, self.res1.bse, DECIMAL_4)

    def test_pvalue(self):
        assert_almost_equal(
            self.Ttest.pvalue,
            student_t.sf(np.abs(self.res1.tvalues), self.res1.model.df_resid) * 2,
            DECIMAL_4,
        )

    def test_df_denom(self):
        assert_equal(self.Ttest.df_denom, self.res1.model.df_resid)

    def test_effect(self):
        assert_almost_equal(self.Ttest.effect, self.res1.params)


class TestTtest2:
    """
    Tests the hypothesis that the coefficients on POP and YEAR
    are equal.

    Results from RPy using 'car' package.
    """

    @classmethod
    def setup_class(cls):
        R = np.zeros(7)
        R[4:6] = [1, -1]
        data = longley.load()
        data.exog = add_constant(data.exog, prepend=False)
        res1 = OLS(data.endog, data.exog).fit()
        cls.Ttest1 = res1.t_test(R)

    def test_tvalue(self):
        assert_almost_equal(self.Ttest1.tvalue, -4.0167754636397284, DECIMAL_4)

    def test_sd(self):
        assert_almost_equal(self.Ttest1.sd, 455.39079425195314, DECIMAL_4)

    def test_pvalue(self):
        assert_almost_equal(self.Ttest1.pvalue, 2 * 0.0015163772380932246, DECIMAL_4)

    def test_df_denom(self):
        assert_equal(self.Ttest1.df_denom, 9)

    def test_effect(self):
        assert_almost_equal(self.Ttest1.effect, -1829.2025687186533, DECIMAL_4)


class TestGLS:
    """
    These test results were obtained by replication with R.
    """

    @classmethod
    def setup_class(cls):
        from .results.results_regression import LongleyGls

        data = longley.load()
        exog = add_constant(
            np.column_stack((data.exog.iloc[:, 1], data.exog.iloc[:, 4])),
            prepend=False,
        )
        tmp_results = OLS(data.endog, exog).fit()
        rho = np.corrcoef(tmp_results.resid[1:], tmp_results.resid[:-1])[0][
            1
        ]  # by assumption
        order = toeplitz(np.arange(16))
        sigma = rho**order
        GLS_results = GLS(data.endog, exog, sigma=sigma).fit()
        cls.res1 = GLS_results
        cls.res2 = LongleyGls()
        # attach for test_missing
        cls.sigma = sigma
        cls.exog = exog
        cls.endog = data.endog

    def test_aic(self):
        # Note: tolerance is tight; rtol=3e-3 fails while 4e-3 passes
        assert_allclose(self.res1.aic + 2, self.res2.aic, rtol=4e-3)

    def test_bic(self):
        # Note: tolerance is tight; rtol=1e-2 fails while 1.5e-2 passes
        assert_allclose(self.res1.bic, self.res2.bic, rtol=1.5e-2)

    def test_loglike(self):
        assert_almost_equal(self.res1.llf, self.res2.llf, DECIMAL_0)

    def test_params(self):
        assert_almost_equal(self.res1.params, self.res2.params, DECIMAL_1)

    def test_resid(self):
        assert_almost_equal(self.res1.resid, self.res2.resid, DECIMAL_4)

    def test_scale(self):
        assert_almost_equal(self.res1.scale, self.res2.scale, DECIMAL_4)

    def test_tvalues(self):
        assert_almost_equal(self.res1.tvalues, self.res2.tvalues, DECIMAL_4)

    def test_standarderrors(self):
        assert_almost_equal(self.res1.bse, self.res2.bse, DECIMAL_4)

    def test_fittedvalues(self):
        assert_almost_equal(self.res1.fittedvalues, self.res2.fittedvalues, DECIMAL_4)

    def test_pvalues(self):
        assert_almost_equal(self.res1.pvalues, self.res2.pvalues, DECIMAL_4)

    def test_missing(self):
        endog = self.endog.copy()  # copy or changes endog for other methods
        endog[[4, 7, 14]] = np.nan
        mod = GLS(endog, self.exog, sigma=self.sigma, missing="drop")
        assert_equal(mod.endog.shape[0], 13)
        assert_equal(mod.exog.shape[0], 13)
        assert_equal(mod.sigma.shape, (13, 13))


class TestGLS_alt_sigma(CheckRegressionResults):
    """
    Test that GLS with no argument is equivalent to OLS.
    """

    @classmethod
    def setup_class(cls):
        data = longley.load()
        endog = np.asarray(data.endog)
        exog = np.asarray(data.exog)
        exog = add_constant(exog, prepend=False)
        ols_res = OLS(endog, exog).fit()
        gls_res = GLS(endog, exog).fit()
        gls_res_scalar = GLS(endog, exog, sigma=1)
        cls.endog = endog
        cls.exog = exog
        cls.res1 = gls_res
        cls.res2 = ols_res
        cls.res3 = gls_res_scalar

    #        self.res2.conf_int = self.res2.conf_int()

    def test_wrong_size_sigma_1d(self):
        n = len(self.endog)
        with pytest.raises(ValueError):
            GLS(self.endog, self.exog, sigma=np.ones(n - 1))

    def test_wrong_size_sigma_2d(self):
        n = len(self.endog)
        with pytest.raises(ValueError):
            GLS(
                self.endog,
                self.exog,
                sigma=np.ones((n - 1, n - 1)),
            )

    @pytest.mark.skip("Test does not raise but should")
    def test_singular_sigma(self):
        n = len(self.endog)
        sigma = np.ones((n, n)) + np.diag(np.ones(n))
        sigma[0, 1] = sigma[1, 0] = 2
        assert np.linalg.matrix_rank(sigma) == n - 1
        with pytest.raises(np.linalg.LinAlgError):
            GLS(self.endog, self.exog, sigma=sigma)


# FIXME: do not leave commented-out, use or move/remove
#    def check_confidenceintervals(self, conf1, conf2):
#        assert_almost_equal(conf1, conf2, DECIMAL_4)


class TestLM:
    @classmethod
    def setup_class(cls):
        # TODO: Test HAC method
        rs = np.random.RandomState(1234)
        x = rs.randn(100, 3)
        b = np.ones((3, 1))
        e = rs.randn(100, 1)
        y = np.dot(x, b) + e
        # Cases?
        # Homoskedastic
        # HC0
        cls.res1_full = OLS(y, x).fit()
        cls.res1_restricted = OLS(y, x[:, 0]).fit()

        cls.res2_full = cls.res1_full.get_robustcov_results("HC0")
        cls.res2_restricted = cls.res1_restricted.get_robustcov_results("HC0")

        cls.x = x
        cls.Y = y

    def test_LM_homoskedastic(self):
        resid = self.res1_restricted.wresid
        n = resid.shape[0]
        x = self.x
        S = np.dot(resid, resid) / n * np.dot(x.T, x) / n
        Sinv = np.linalg.inv(S)
        s = np.mean(x * resid[:, None], 0)
        LMstat = n * np.dot(np.dot(s, Sinv), s.T)
        LMstat_OLS = self.res1_full.compare_lm_test(self.res1_restricted)
        LMstat2 = LMstat_OLS[0]
        assert_almost_equal(LMstat, LMstat2, DECIMAL_7)

    def test_LM_heteroskedastic_nodemean(self):
        resid = self.res1_restricted.wresid
        n = resid.shape[0]
        x = self.x
        scores = x * resid[:, None]
        S = np.dot(scores.T, scores) / n
        Sinv = np.linalg.inv(S)
        s = np.mean(scores, 0)
        LMstat = n * np.dot(np.dot(s, Sinv), s.T)
        LMstat_OLS = self.res2_full.compare_lm_test(self.res2_restricted, demean=False)
        LMstat2 = LMstat_OLS[0]
        assert_almost_equal(LMstat, LMstat2, DECIMAL_7)

    def test_LM_heteroskedastic_demean(self):
        resid = self.res1_restricted.wresid
        n = resid.shape[0]
        x = self.x
        scores = x * resid[:, None]
        scores_demean = scores - scores.mean(0)
        S = np.dot(scores_demean.T, scores_demean) / n
        Sinv = np.linalg.inv(S)
        s = np.mean(scores, 0)
        LMstat = n * np.dot(np.dot(s, Sinv), s.T)
        LMstat_OLS = self.res2_full.compare_lm_test(self.res2_restricted)
        LMstat2 = LMstat_OLS[0]
        assert_almost_equal(LMstat, LMstat2, DECIMAL_7)

    def test_LM_heteroskedastic_LRversion(self):
        resid = self.res1_restricted.wresid
        resid_full = self.res1_full.wresid
        n = resid.shape[0]
        x = self.x
        scores = x * resid[:, None]
        s = np.mean(scores, 0)
        scores = x * resid_full[:, None]
        S = np.dot(scores.T, scores) / n
        Sinv = np.linalg.inv(S)
        LMstat = n * np.dot(np.dot(s, Sinv), s.T)
        LMstat_OLS = self.res2_full.compare_lm_test(self.res2_restricted, use_lr=True)
        LMstat2 = LMstat_OLS[0]
        assert_almost_equal(LMstat, LMstat2, DECIMAL_7)

    def test_LM_nonnested(self):
        with pytest.raises(ValueError):
            self.res2_restricted.compare_lm_test(self.res2_full)


class TestOLS_GLS_WLS_equivalence:
    @classmethod
    def setup_class(cls):
        data = longley.load()
        data.exog = add_constant(data.exog, prepend=False)
        y = data.endog
        x = data.exog
        n = y.shape[0]
        w = np.ones(n)
        cls.results = []
        cls.results.append(OLS(y, x).fit())
        cls.results.append(WLS(y, x, w).fit())
        # scaling weights does not change main results (except scale)
        cls.results.append(GLS(y, x, 100 * w).fit())
        cls.results.append(GLS(y, x, np.diag(0.1 * w)).fit())

    def test_ll(self):
        llf = np.array([r.llf for r in self.results])
        llf_1 = np.ones_like(llf) * self.results[0].llf
        assert_almost_equal(llf, llf_1, DECIMAL_7)

        ic = np.array([r.aic for r in self.results])
        ic_1 = np.ones_like(ic) * self.results[0].aic
        assert_almost_equal(ic, ic_1, DECIMAL_7)

        ic = np.array([r.bic for r in self.results])
        ic_1 = np.ones_like(ic) * self.results[0].bic
        assert_almost_equal(ic, ic_1, DECIMAL_7)

    def test_params(self):
        params = np.array([r.params for r in self.results])
        params_1 = np.array([self.results[0].params] * len(self.results))
        assert_allclose(params, params_1)

    def test_ss(self):
        bse = np.array([r.bse for r in self.results])
        bse_1 = np.array([self.results[0].bse] * len(self.results))
        assert_allclose(bse, bse_1)

    def test_rsquared(self):
        rsquared = np.array([r.rsquared for r in self.results])
        rsquared_1 = np.array([self.results[0].rsquared] * len(self.results))
        assert_almost_equal(rsquared, rsquared_1, DECIMAL_7)


class TestGLS_WLS_equivalence(TestOLS_GLS_WLS_equivalence):
    # reuse test methods

    @classmethod
    def setup_class(cls):
        data = longley.load()
        data.exog = add_constant(data.exog, prepend=False)
        y = data.endog
        x = data.exog
        n = y.shape[0]
        np.random.seed(5)
        w = np.random.uniform(0.5, 1, n)
        w_inv = 1.0 / w
        cls.results = []
        cls.results.append(WLS(y, x, w).fit())
        # scaling weights does not change main results (except scale)
        cls.results.append(WLS(y, x, 0.01 * w).fit())
        cls.results.append(GLS(y, x, 100 * w_inv).fit())
        cls.results.append(GLS(y, x, np.diag(0.1 * w_inv)).fit())


class TestNonFit:
    @classmethod
    def setup_class(cls):
        data = longley.load()
        data.exog = add_constant(data.exog, prepend=False)
        cls.endog = data.endog
        cls.exog = data.exog
        cls.ols_model = OLS(data.endog, data.exog)

    def test_df_resid(self):
        df_resid = self.endog.shape[0] - self.exog.shape[1]
        assert_equal(self.ols_model.df_resid, df_resid)


class TestWLS_CornerCases:
    @classmethod
    def setup_class(cls):
        cls.exog = np.ones((1,))
        cls.endog = np.ones((1,))
        weights = 1
        cls.wls_res = WLS(cls.endog, cls.exog, weights=weights).fit()

    def test_wrong_size_weights(self):
        weights = np.ones((10, 10))
        with pytest.raises(ValueError):
            WLS(self.endog, self.exog, weights=weights)


class TestWLSExogWeights(CheckRegressionResults):
    # Test WLS with Greene's credit card data
    # reg avgexp age income incomesq ownrent [aw=1/incomesq]
    @classmethod
    def setup_class(cls):
        from statsmodels.datasets.ccard import load

        from .results.results_regression import CCardWLS

        dta = load()
        endog = np.asarray(dta.endog)
        exog = np.asarray(dta.exog)
        exog = add_constant(exog, prepend=False)
        nobs = 72.0

        weights = 1 / exog[:, 2]
        # for comparison with stata analytic weights
        scaled_weights = (weights * nobs) / weights.sum()

        cls.res1 = WLS(endog, exog, weights=scaled_weights).fit()
        cls.res2 = CCardWLS()
        cls.res2.wresid = scaled_weights**0.5 * cls.res2.resid

        # correction because we use different definition for loglike/llf
        corr_ic = 2 * (cls.res1.llf - cls.res2.llf)
        cls.res2.aic -= corr_ic
        cls.res2.bic -= corr_ic
        cls.res2.llf += 0.5 * np.sum(np.log(cls.res1.model.weights))


def test_wls_example():
    # example from the docstring, there was a note about a bug, should
    # be fixed now
    Y = [1, 3, 4, 5, 2, 3, 4]
    x = lrange(1, 8)
    x = add_constant(x, prepend=False)
    wls_model = WLS(Y, x, weights=lrange(1, 8)).fit()
    # taken from R lm.summary
    assert_almost_equal(wls_model.fvalue, 0.127337843215, 6)
    assert_almost_equal(wls_model.scale, 2.44608530786**2, 6)


def test_wls_tss():
    y = np.array([22, 22, 22, 23, 23, 23])
    x = [[1, 0], [1, 0], [1, 1], [0, 1], [0, 1], [0, 1]]

    ols_mod = OLS(y, add_constant(x, prepend=False)).fit()

    yw = np.array([22, 22, 23.0])
    Xw = [[1, 0], [1, 1], [0, 1]]
    w = np.array([2, 1, 3.0])

    wls_mod = WLS(yw, add_constant(Xw, prepend=False), weights=w).fit()
    assert_equal(ols_mod.centered_tss, wls_mod.centered_tss)


class TestWLSScalarVsArray(CheckRegressionResults):
    @classmethod
    def setup_class(cls):
        from statsmodels.datasets.longley import load

        dta = load()
        endog = np.asarray(dta.endog)
        exog = np.asarray(dta.exog)
        exog = add_constant(exog, prepend=True)
        wls_scalar = WLS(endog, exog, weights=1.0 / 3).fit()
        weights = [1 / 3.0] * len(endog)
        wls_array = WLS(endog, exog, weights=weights).fit()
        cls.res1 = wls_scalar
        cls.res2 = wls_array


class TestWLS_GLS(CheckRegressionResults):
    @classmethod
    def setup_class(cls):
        from statsmodels.datasets.ccard import load

        data = load()
        endog = np.asarray(data.endog)
        exog = np.asarray(data.exog)
        sigma = exog[:, 2]
        cls.res1 = WLS(endog, exog, weights=1 / sigma).fit()
        cls.res2 = GLS(endog, exog, sigma=sigma).fit()

    def check_confidenceintervals(self, conf1, conf2):  # FIXME: never called
        assert_almost_equal(conf1, conf2(), DECIMAL_4)


def test_wls_missing():
    from statsmodels.datasets.ccard import load

    data = load()
    endog = data.endog
    endog[[10, 25]] = np.nan
    mod = WLS(data.endog, data.exog, weights=1 / data.exog.iloc[:, 2], missing="drop")
    assert_equal(mod.endog.shape[0], 70)
    assert_equal(mod.exog.shape[0], 70)
    assert_equal(mod.weights.shape[0], 70)


class TestWLS_OLS(CheckRegressionResults):
    @classmethod
    def setup_class(cls):
        data = longley.load()
        endog = np.asarray(data.endog)
        exog = np.asarray(data.exog)
        exog = add_constant(exog, prepend=False)
        cls.res1 = OLS(endog, exog).fit()
        cls.res2 = WLS(endog, exog).fit()

    def check_confidenceintervals(self, conf1, conf2):  # FIXME: never called
        assert_almost_equal(conf1, conf2(), DECIMAL_4)


class TestGLS_OLS(CheckRegressionResults):
    @classmethod
    def setup_class(cls):
        data = longley.load()
        endog = np.asarray(data.endog)
        exog = np.asarray(data.exog)
        exog = add_constant(exog, prepend=False)
        cls.res1 = GLS(endog, exog).fit()
        cls.res2 = OLS(endog, exog).fit()

    def check_confidenceintervals(self, conf1, conf2):  # FIXME: never called
        assert_almost_equal(conf1, conf2(), DECIMAL_4)


# FIXME: do not leave this commented-out sitting here
# TODO: test AR
# why the two-stage in AR?
# class TestAR:
#     from statsmodels.datasets.sunspots import load
#     data = load()
#     model = AR(data.endog, rho=4).fit()
#     R_res = RModel(data.endog, aic="FALSE", order_max=4)#

#     def test_params(self):
#         assert_almost_equal(self.model.rho,
#         pass

#     def test_order(self):
# In R this can be defined or chosen by minimizing the AIC if aic=True
#        pass


class TestYuleWalker:
    @classmethod
    def setup_class(cls):
        from statsmodels.datasets.sunspots import load

        data = load()
        cls.rho, cls.sigma = yule_walker(data.endog, order=4, method="mle")
        cls.R_params = [
            1.2831003105694765,
            -0.45240924374091945,
            -0.20770298557575195,
            0.047943648089542337,
        ]

    def test_params(self):
        assert_almost_equal(self.rho, self.R_params, DECIMAL_4)


class TestDataDimensions(CheckRegressionResults):
    @classmethod
    def setup_class(cls):
        np.random.seed(54321)
        cls.endog_n_ = np.random.uniform(0, 20, size=30)
        cls.endog_n_one = cls.endog_n_[:, None]
        cls.exog_n_ = np.random.uniform(0, 20, size=30)
        cls.exog_n_one = cls.exog_n_[:, None]
        cls.degen_exog = cls.exog_n_one[:-1]
        cls.mod1 = OLS(cls.endog_n_one, cls.exog_n_one)
        cls.mod1.df_model += 1
        cls.res1 = cls.mod1.fit()
        # Note that these are created for every subclass..
        # A little extra overhead probably
        cls.mod2 = OLS(cls.endog_n_one, cls.exog_n_one)
        cls.mod2.df_model += 1
        cls.res2 = cls.mod2.fit()

    def check_confidenceintervals(self, conf1, conf2):  # FIXME: never called
        assert_almost_equal(conf1, conf2(), DECIMAL_4)


class TestGLS_large_data(TestDataDimensions):
    @classmethod
    def setup_class(cls):
        super().setup_class()
        nobs = 1000
        y = np.random.randn(nobs, 1)
        x = np.random.randn(nobs, 20)
        sigma = np.ones_like(y)
        cls.gls_res = GLS(y, x, sigma=sigma).fit()
        cls.gls_res_scalar = GLS(y, x, sigma=1).fit()
        cls.gls_res_none = GLS(y, x).fit()
        cls.ols_res = OLS(y, x).fit()

    def test_large_equal_params(self):
        assert_almost_equal(self.ols_res.params, self.gls_res.params, DECIMAL_7)

    def test_large_equal_loglike(self):
        assert_almost_equal(self.ols_res.llf, self.gls_res.llf, DECIMAL_7)

    def test_large_equal_params_none(self):
