1# -*- coding: utf-8 -*-
2"""
3Created on Mon Aug 16 16:19:55 2021
4
5Author: Josef Perktod
6License: BSD-3
7"""
8import numpy as np
9from scipy import stats
10# import matplotlib.pyplot as plt
11from numpy.testing import assert_allclose
12
13from statsmodels.base.model import GenericLikelihoodModel
14from statsmodels.distributions.copula.api import (
15    ClaytonCopula, GaussianCopula, FrankCopula,
16    GumbelCopula, IndependenceCopula, CopulaDistribution)
17
18from statsmodels.distributions.copula import depfunc_ev as dep
19from statsmodels.distributions.copula.extreme_value import ExtremeValueCopula
20
21
22class CopulaModel(GenericLikelihoodModel):
23
24    def __init__(self, copula_distribution, endog, k_params=None):
25        self.copula_distribution = copula_distribution
26        self.endog = endog
27        self.exog = None
28        if k_params is None:
29            k_params = 1
30        self.nparams = k_params
31        self.k_copparams = 1
32        super().__init__(endog, self.exog)
33
34    def split_params(self, params):
35        pass
36
37    def loglike(self, params):
38        params = np.atleast_1d(params)
39        cd = self.copula_distribution
40        # ll = cd.logpdf(self.endog, args=(params[:2], params[2:]))
41        cop_args = params[:self.k_copparams]
42        if cop_args.size == 0:
43            cop_args = ()
44        if len(params) > self.k_copparams:
45            marg_args = np.split(params[self.k_copparams:], 2)
46        else:
47            marg_args = None
48
49        ll = cd.logpdf(self.endog,
50                       cop_args=cop_args, marg_args=marg_args
51                       ).sum()
52        return ll
53
54
55def get_data(nobs):
56    cop_f = FrankCopula(theta=2)
57    cd_f = CopulaDistribution(cop_f, [stats.norm, stats.norm])
58    # np.random.seed(98645713)
59    # at some seeds, parameters atol-differ from true
60    # TODO: setting seed doesn't work for copula,
61    # copula creates new randomly initialized random state, see #7650
62    rng = np.random.RandomState(98645713)
63    rvs = cd_f.rvs(nobs, random_state=rng)
64    assert_allclose(rvs.mean(0), [-0.02936002,  0.06658304], atol=1e-7)
65    return rvs
66
67
68data_ev = get_data(500)
69
70
71class CheckEVfit1(object):
72
73    def test(self):
74        cop = self.copula
75        args = self.cop_args
76
77        cev = CopulaDistribution(cop, [stats.norm, stats.norm], cop_args=None)
78        k_marg = 4
79        mod = CopulaModel(cev, data_ev + [0.5, -0.1],
80                          k_params=self.k_copparams + k_marg)
81
82        # TODO: patching for now
83        mod.k_copparams = self.k_copparams
84        mod.df_resid = len(mod.endog) - mod.nparams
85        mod.df_model = mod.nparams - 0
86        res = mod.fit(start_params=list(args) + [0.5, 1, -0.1, 1],
87                      method="bfgs")
88        res = mod.fit(method="newton", start_params=res.params)
89
90        assert mod.nparams == self.k_copparams + k_marg
91        assert res.nobs == len(mod.endog)
92        assert_allclose(res.params[[-4, -2]], [0.5, -0.1], atol=0.2)
93        res.summary()
94        assert res.mle_retvals["converged"]
95        assert not np.isnan(res.bse).any()
96
97    def test2m(self):
98        cop = self.copula
99        args = self.cop_args
100
101        cev = CopulaDistribution(cop, [stats.norm, stats.norm], cop_args=None)
102        k_marg = 2
103        mod = CopulaModel(cev, data_ev + [0.5, -0.1],
104                          k_params=self.k_copparams + k_marg)
105
106        # TODO: patching for now
107        mod.k_copparams = self.k_copparams
108        mod.df_resid = len(mod.endog) - mod.nparams
109        mod.df_model = mod.nparams - 0
110        res = mod.fit(start_params=list(args) + [0.5, -0.1],
111                      method="bfgs")
112        # the following fails in TestEVAsymLogistic with nan loglike
113        # res = mod.fit(method="newton", start_params=res.params)
114
115        assert mod.nparams == self.k_copparams + k_marg
116        assert res.nobs == len(mod.endog)
117        assert_allclose(res.params[[-2, -1]], [0.5, -0.1], atol=0.2)
118        res.summary()
119        assert res.mle_retvals["converged"]
120        assert not np.isnan(res.bse).any()
121
122
123# temporarily split for copulas that only have fixed cop_args
124class CheckEVfit0(object):
125
126    def test0(self):
127        # test with fixed copula params
128        cop = getattr(self, "copula_fixed", None)
129        if cop is None:
130            # skip test if not yet available
131            return
132        args = self.cop_args
133
134        cev = CopulaDistribution(cop, [stats.norm, stats.norm], cop_args=args)
135        k_marg = 2
136        mod = CopulaModel(cev, data_ev + [0.5, -0.1],
137                          k_params=0 + k_marg)
138
139        # TODO: patching for now
140        mod.k_copparams = 0
141        mod.df_resid = len(mod.endog) - mod.nparams
142        mod.df_model = mod.nparams - 0
143        res = mod.fit(start_params=[0.5, -0.1],
144                      method="bfgs")
145        # the following fails in TestEVAsymLogistic with nan loglike
146        # res = mod.fit(method="newton", start_params=res.params)
147
148        assert mod.nparams == 0 + k_marg
149        assert res.nobs == len(mod.endog)
150        assert_allclose(res.params, [0.5, -0.1], atol=0.2)
151        res.summary()
152        assert res.mle_retvals["converged"]
153        assert not np.isnan(res.bse).any()
154
155
156class CheckEVfit(CheckEVfit1, CheckEVfit0):
157    # unit test mainly for arg handling, not to verify parameter estimates
158    pass
159
160
161class TestEVHR(CheckEVfit):
162
163    @classmethod
164    def setup_class(cls):
165        cls.copula = ExtremeValueCopula(transform=dep.HR())
166        cls.cop_args = (1,)
167        cls.k_copparams = 1
168        cls.copula_fixed = ExtremeValueCopula(transform=dep.HR(),
169                                              args=cls.cop_args)
170
171
172class TestEVAsymLogistic(CheckEVfit):
173
174    @classmethod
175    def setup_class(cls):
176        cls.copula = ExtremeValueCopula(transform=dep.AsymLogistic())
177        cls.cop_args = (0.1, 0.7, 0.7)
178        cls.k_copparams = 3
179        cls.copula_fixed = ExtremeValueCopula(transform=dep.AsymLogistic(),
180                                              args=cls.cop_args)
181
182
183class TestEVAsymMixed(CheckEVfit):
184
185    @classmethod
186    def setup_class(cls):
187        cls.copula = ExtremeValueCopula(transform=dep.AsymMixed())
188        cls.cop_args = (0.5, 0.05)
189        cls.k_copparams = 2
190        cls.copula_fixed = ExtremeValueCopula(transform=dep.AsymMixed(),
191                                              args=cls.cop_args)
192
193
194class TestFrank(CheckEVfit):
195
196    @classmethod
197    def setup_class(cls):
198        cls.copula = FrankCopula()
199        cls.cop_args = (0.5,)
200        cls.k_copparams = 1
201        cls.copula_fixed = FrankCopula(*cls.cop_args)
202
203
204class TestGaussian(CheckEVfit):
205
206    @classmethod
207    def setup_class(cls):
208        cls.copula = GaussianCopula()
209        cls.cop_args = ()
210        cls.k_copparams = 0
211
212
213class TestClayton(CheckEVfit):
214
215    @classmethod
216    def setup_class(cls):
217        cls.copula = ClaytonCopula()
218        cls.cop_args = (1.01,)
219        cls.k_copparams = 1
220        cls.copula_fixed = ClaytonCopula(*cls.cop_args)
221
222
223class TestGumbel(CheckEVfit):
224
225    @classmethod
226    def setup_class(cls):
227        cls.copula = GumbelCopula()
228        cls.cop_args = (1.01,)
229        cls.k_copparams = 1
230        cls.copula_fixed = GumbelCopula(*cls.cop_args)
231
232
233class TestIndependence(CheckEVfit0):
234
235    @classmethod
236    def setup_class(cls):
237        cls.copula = IndependenceCopula()
238        cls.cop_args = ()
239        cls.k_copparams = 0
240        cls.copula_fixed = IndependenceCopula(*cls.cop_args)
241