# -*- coding: utf-8 -*-
"""Unit tests for generic score/LM tests and conditional moment tests

Created on Mon Nov 17 08:44:06 2014

Author: Josef Perktold
License: BSD-3

"""

import numpy as np
from numpy.testing import assert_allclose

from statsmodels.regression.linear_model import OLS
from statsmodels.stats._diagnostic_other import CMTNewey, CMTTauchen
import statsmodels.stats._diagnostic_other as diao


class CheckCMT(object):

    def test_score(self):
        expected = self.results_score
        for msg, actual in self.res_score():
            # not all cases provide all 3 elements,
            # TODO: fix API, returns of functions
            assert_allclose(actual, expected[:np.size(actual)], rtol=1e-13,
                            err_msg=msg)

    def test_scorehc0(self):
        expected = self.results_hc0
        for msg, actual in self.res_hc0():
            assert_allclose(actual, expected[:np.size(actual)], rtol=1e-13,
                            err_msg=msg)

    def test_scoreopg(self):
        expected = self.results_opg
        for msg, actual in self.res_opg():
            assert_allclose(actual, expected[:np.size(actual)], rtol=1e-13,
                            err_msg=msg)


class TestCMTOLS(CheckCMT):

    @classmethod
    def setup_class(cls):

        np.random.seed(864259)

        nobs, k_vars = 100, 4
        sig_e = 1
        x0 = np.random.randn(nobs, k_vars)
        x0[:,0] = 1
        y_true = x0.sum(1)
        y = y_true + sig_e * np.random.randn(nobs)

        x1 = np.random.randn(nobs, 2)
        x = np.column_stack((x0, x1))

        cls.exog_full = x
        cls.exog_add = x1
        cls.res_ols = OLS(y, x0).fit()

        cls.attach_moment_conditions()

        # results from initial run, not reference package (drooped 2 digits)
        cls.results_score = (1.6857659627548, 0.43046770240535, 2)
        cls.results_hc0 = (1.6385932313952, 0.4407415561953, 2)
        cls.results_opg = (1.72226002418488, 0.422684174119544, 2.0)

    @classmethod    # TODO: a better structure ?
    def attach_moment_conditions(self):
        res_ols = self.res_ols
        # assumes x = column_stack([x0, x1])
        x = self.exog_full
        #x0 = self.res_ols.model.exog  # not used here
        x1 = self.exog_add

        nobs, k_constraints = x1.shape

        # TODO: cleanup after initial copy past
        moms_obs = res_ols.resid[:, None] * x
        moms = moms_obs.sum(0)
        cov_moms =  res_ols.mse_resid * x.T.dot(x) #np.linalg.inv(x.T.dot(x))
        cov_moms *=  res_ols.df_resid / nobs

        # weights used for GMM to replicate OLS
        weights = np.linalg.inv(cov_moms)
        # we don't use last two variables
        weights[:, -k_constraints:] = 0
        weights[-k_constraints:, :] = 0

        k_moms = moms.shape[0]
        # TODO: Newey has different versions that all produce the same result
        #       in example
        L = np.eye(k_moms)[-k_constraints:] #.dot(np.linalg.inv(cov_moms))

        moms_deriv = cov_moms[:, :-k_constraints]

        covm = moms_obs.T.dot(moms_obs)

        #attach
        self.nobs = nobs
        self.moms = moms
        self.moms_obs = moms_obs
        self.cov_moms = cov_moms
        self.covm = covm
        self.moms_deriv = moms_deriv
        self.weights = weights
        self.L = L

    def res_score(self):
        res_ols = self.res_ols
        nobs = self.nobs
        moms = self.moms
        moms_obs = self.moms_obs
        cov_moms = self.cov_moms
        covm = self.covm
        moms_deriv = self.moms_deriv
        weights = self.weights
        L = self.L
        x = self.exog_full  # for auxiliary regression only

        res_all = []

        # auxiliary regression
        stat = nobs * OLS(res_ols.resid, x).fit().rsquared
        res_all.append(('ols R2', stat))

        stat = moms.dot(np.linalg.solve(cov_moms, moms))
        res_all.append(('score simple', stat))

        tres = diao.lm_robust(moms, np.eye(moms.shape[0])[-2:],
                              np.linalg.inv(cov_moms), cov_moms)
        res_all.append(('score mle', tres))

        tres = CMTNewey(moms, cov_moms, cov_moms[:,:-2], weights, L).chisquare
        res_all.append(('Newey', tres))

        tres = CMTTauchen(moms[:-2], cov_moms[:-2, :-2], moms[-2:],
                          cov_moms[-2:, :-2], cov_moms).chisquare
        res_all.append(('Tauchen', tres))

        return res_all

    def res_opg(self):
        res_ols = self.res_ols
        nobs = self.nobs
        moms = self.moms
        moms_obs = self.moms_obs
        covm = self.covm
        moms_deriv = self.moms_deriv
        weights = self.weights
        L = self.L
        x = self.exog_full

        res_ols2_hc0 = OLS(res_ols.model.endog, x).fit(cov_type='HC0')

        res_all = []

        # auxiliary regression
        ones = np.ones(nobs)
        stat = nobs * OLS(ones, moms_obs).fit().rsquared
        res_all.append(('ols R2', stat))

        tres = res_ols2_hc0.compare_lm_test(res_ols, demean=False)
        res_all.append(('comp_lm uc', tres))

        tres = CMTNewey(moms, covm, covm[:,:-2], weights, L).chisquare
        res_all.append(('Newey', tres))

        tres = CMTTauchen(moms[:-2], covm[:-2, :-2], moms[-2:], covm[-2:, :-2],
                          covm).chisquare
        res_all.append(('Tauchen', tres))

        tres = diao.lm_robust_subset(moms[-2:], 2, covm, covm)
        res_all.append(('score subset QMLE', tres))

        tres = diao.lm_robust(moms, np.eye(moms.shape[0])[-2:],
                              np.linalg.inv(covm), covm, cov_params=None)
        res_all.append(('scoreB QMLE', tres))

        tres = diao.lm_robust(moms, np.eye(moms.shape[0])[-2:],
                              np.linalg.inv(covm), None,
                              cov_params=np.linalg.inv(covm))
        res_all.append(('scoreV QMLE', tres))

        return res_all

    def res_hc0(self):
        res_ols = self.res_ols
        nobs = self.nobs
        moms = self.moms
        moms_obs = self.moms_obs
        cov_moms = self.cov_moms   # Hessian with scale
        covm = self.covm
        moms_deriv = self.moms_deriv
        weights = self.weights
        L = self.L

        x0 = res_ols.model.exog
        x1 = self.exog_add

        res_all = []

        tres = diao.cm_test_robust(resid=res_ols.resid, resid_deriv=x0,
                                   instruments=x1, weights=1)
        # TODO: extra return and no df in cm_test_robust Wooldridge
        res_all.append(('Wooldridge', tres[:2]))

        tres = CMTNewey(moms, covm, moms_deriv, weights, L).chisquare
        res_all.append(('Newey', tres))

        tres = CMTTauchen(moms[:-2], cov_moms[:-2, :-2], moms[-2:],
                          cov_moms[-2:, :-2], covm).chisquare
        res_all.append(('Tauchen', tres))

        tres = diao.lm_robust_subset(moms[-2:], 2, cov_moms, covm)
        res_all.append(('score subset QMLE', tres))

        tres = diao.lm_robust(moms, np.eye(moms.shape[0])[-2:],
                              np.linalg.inv(cov_moms), covm)
        res_all.append(('scoreB QMLE', tres))

        # need sandwich cov_params V
        Ainv = np.linalg.inv(cov_moms)
        vv = Ainv.dot(covm).dot(Ainv)
        tres = diao.lm_robust(moms, np.eye(moms.shape[0])[-2:],
                              np.linalg.inv(cov_moms), None,
                              cov_params=vv)
        res_all.append(('scoreV QMLE', tres))

        tres = diao.conditional_moment_test_generic(moms_obs[:, -2:],
                                                    cov_moms[-2:, :-2],
                                                    moms_obs[:,:-2],
                                                    cov_moms[:-2, :-2])
        tres_ = (tres.stat_cmt, tres.pval_cmt)
        res_all.append(('cmt', tres_))

        # using unscaled hessian instead of scaled
        x = self.exog_full
        hess_unscaled = x.T.dot(x)
        tres = diao.conditional_moment_test_generic(moms_obs[:, -2:],
                    hess_unscaled[-2:, :-2], moms_obs[:,:-2],
                    hess_unscaled[:-2, :-2])#, covm)
        tres_ = (tres.stat_cmt, tres.pval_cmt)
        res_all.append(('cmt', tres_))

        score_deriv_uu = cov_moms[:-2, :-2]
        score_deriv_cu = cov_moms[-2:, :-2]
        cov_score_cc = covm[-2:, -2:]
        cov_score_cu = covm[-2:, :-2]
        cov_score_uu = covm[:-2, :-2]
        moms[-2:], 2, cov_moms, covm
        tres = diao.lm_robust_subset_parts(moms[-2:], 2, score_deriv_uu,
                                     score_deriv_cu, cov_score_cc,
                                     cov_score_cu, cov_score_uu)

        res_all.append(('score subset_parts QMLE', tres))

        params_deriv = np.eye(x.shape[1], x.shape[1] - 2)
        #params_deriv[[-2, -1], [-2, -1]] = 0
        score = moms
        score_deriv = cov_moms
        cov_score = covm

        tres = diao.lm_robust_reparameterized(score, params_deriv,
                           score_deriv, cov_score)

        res_all.append(('score reparam QMLE', tres))

        return res_all