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