1""" 2Tests for VARMAX models 3 4Author: Chad Fulton 5License: Simplified-BSD 6""" 7import os 8import re 9import warnings 10 11import numpy as np 12from numpy.testing import assert_equal, assert_raises, assert_allclose 13import pandas as pd 14import pytest 15 16from statsmodels.tsa.statespace import dynamic_factor 17from .results import results_varmax, results_dynamic_factor 18from statsmodels.iolib.summary import forg 19 20current_path = os.path.dirname(os.path.abspath(__file__)) 21 22output_path = os.path.join('results', 'results_dynamic_factor_stata.csv') 23output_results = pd.read_csv(os.path.join(current_path, output_path)) 24 25 26class CheckDynamicFactor(object): 27 @classmethod 28 def setup_class(cls, true, k_factors, factor_order, cov_type='approx', 29 included_vars=['dln_inv', 'dln_inc', 'dln_consump'], 30 demean=False, filter=True, **kwargs): 31 cls.true = true 32 # 1960:Q1 - 1982:Q4 33 dta = pd.DataFrame( 34 results_varmax.lutkepohl_data, columns=['inv', 'inc', 'consump'], 35 index=pd.date_range('1960-01-01', '1982-10-01', freq='QS')) 36 37 dta['dln_inv'] = np.log(dta['inv']).diff() 38 dta['dln_inc'] = np.log(dta['inc']).diff() 39 dta['dln_consump'] = np.log(dta['consump']).diff() 40 41 endog = dta.loc['1960-04-01':'1978-10-01', included_vars] 42 43 if demean: 44 endog -= dta.iloc[1:][included_vars].mean() 45 46 cls.model = dynamic_factor.DynamicFactor(endog, k_factors=k_factors, 47 factor_order=factor_order, 48 **kwargs) 49 50 if filter: 51 cls.results = cls.model.smooth(true['params'], cov_type=cov_type) 52 53 def test_params(self): 54 # Smoke test to make sure the start_params are well-defined and 55 # lead to a well-defined model 56 self.model.filter(self.model.start_params) 57 # Similarly a smoke test for param_names 58 assert_equal(len(self.model.start_params), len(self.model.param_names)) 59 # Finally make sure the transform and untransform do their job 60 actual = self.model.transform_params( 61 self.model.untransform_params(self.model.start_params)) 62 assert_allclose(actual, self.model.start_params) 63 # Also in the case of enforce stationarity = False 64 self.model.enforce_stationarity = False 65 actual = self.model.transform_params( 66 self.model.untransform_params(self.model.start_params)) 67 self.model.enforce_stationarity = True 68 assert_allclose(actual, self.model.start_params) 69 70 def test_results(self, close_figures): 71 # Smoke test for creating the summary 72 with warnings.catch_warnings(): 73 warnings.simplefilter("ignore") 74 self.results.summary() 75 76 # Test cofficient matrix creation 77 # (via a different, more direct, method) 78 if self.model.factor_order > 0: 79 model = self.model 80 k_factors = model.k_factors 81 pft_params = self.results.params[model._params_factor_transition] 82 coefficients = np.array(pft_params).reshape( 83 k_factors, k_factors * model.factor_order) 84 coefficient_matrices = np.array([ 85 coefficients[:self.model.k_factors, 86 i*self.model.k_factors:(i+1)*self.model.k_factors] 87 for i in range(self.model.factor_order) 88 ]) 89 assert_equal( 90 self.results.coefficient_matrices_var, 91 coefficient_matrices) 92 else: 93 assert_equal(self.results.coefficient_matrices_var, None) 94 95 @pytest.mark.matplotlib 96 def test_plot_coefficients_of_determination(self, close_figures): 97 # Smoke test for plot_coefficients_of_determination 98 self.results.plot_coefficients_of_determination() 99 100 def test_no_enforce(self): 101 return 102 # Test that nothing goes wrong when we do not enforce stationarity 103 params = self.model.untransform_params(self.true['params']) 104 params[self.model._params_transition] = ( 105 self.true['params'][self.model._params_transition]) 106 self.model.enforce_stationarity = False 107 results = self.model.filter(params, transformed=False) 108 self.model.enforce_stationarity = True 109 assert_allclose(results.llf, self.results.llf, rtol=1e-5) 110 111 def test_mle(self, init_powell=True): 112 with warnings.catch_warnings(record=True): 113 warnings.simplefilter('always') 114 start_params = self.model.start_params 115 if init_powell: 116 results = self.model.fit(method='powell', 117 maxiter=100, disp=False) 118 start_params = results.params 119 results = self.model.fit(start_params, maxiter=1000, disp=False) 120 results = self.model.fit(results.params, method='nm', maxiter=1000, 121 disp=False) 122 if not results.llf > self.results.llf: 123 assert_allclose(results.llf, self.results.llf, rtol=1e-5) 124 125 def test_loglike(self): 126 assert_allclose(self.results.llf, self.true['loglike'], rtol=1e-6) 127 128 def test_aic(self): 129 # We only get 3 digits from Stata 130 assert_allclose(self.results.aic, self.true['aic'], atol=3) 131 132 def test_bic(self): 133 # We only get 3 digits from Stata 134 assert_allclose(self.results.bic, self.true['bic'], atol=3) 135 136 def test_predict(self, **kwargs): 137 # Tests predict + forecast 138 self.results.predict(end='1982-10-01', **kwargs) 139 assert_allclose( 140 self.results.predict(end='1982-10-01', **kwargs), 141 self.true['predict'], 142 atol=1e-6) 143 144 def test_dynamic_predict(self, **kwargs): 145 # Tests predict + dynamic predict + forecast 146 assert_allclose( 147 self.results.predict(end='1982-10-01', dynamic='1961-01-01', 148 **kwargs), 149 self.true['dynamic_predict'], 150 atol=1e-6) 151 152 153class TestDynamicFactor(CheckDynamicFactor): 154 """ 155 Test for a dynamic factor model with 1 AR(2) factor 156 """ 157 @classmethod 158 def setup_class(cls): 159 true = results_dynamic_factor.lutkepohl_dfm.copy() 160 true['predict'] = output_results.iloc[1:][[ 161 'predict_dfm_1', 'predict_dfm_2', 'predict_dfm_3']] 162 true['dynamic_predict'] = output_results.iloc[1:][[ 163 'dyn_predict_dfm_1', 'dyn_predict_dfm_2', 'dyn_predict_dfm_3']] 164 super(TestDynamicFactor, cls).setup_class( 165 true, k_factors=1, factor_order=2) 166 167 def test_bse_approx(self): 168 bse = self.results._cov_params_approx().diagonal()**0.5 169 assert_allclose(bse, self.true['bse_oim'], atol=1e-5) 170 171 172class TestDynamicFactor2(CheckDynamicFactor): 173 """ 174 Test for a dynamic factor model with two VAR(1) factors 175 """ 176 @classmethod 177 def setup_class(cls): 178 true = results_dynamic_factor.lutkepohl_dfm2.copy() 179 true['predict'] = output_results.iloc[1:][[ 180 'predict_dfm2_1', 'predict_dfm2_2', 'predict_dfm2_3']] 181 true['dynamic_predict'] = output_results.iloc[1:][[ 182 'dyn_predict_dfm2_1', 'dyn_predict_dfm2_2', 'dyn_predict_dfm2_3']] 183 super(TestDynamicFactor2, cls).setup_class( 184 true, k_factors=2, factor_order=1) 185 186 def test_mle(self): 187 # Stata's MLE on this model does not converge, so no reason to check 188 pass 189 190 def test_bse(self): 191 # Stata's MLE on this model does not converge, and four of their 192 # params do not even have bse (possibly they are still at starting 193 # values?), so no reason to check this 194 pass 195 196 def test_aic(self): 197 # Stata uses 9 df (i.e. 9 params) here instead of 13, because since the 198 # model did not coverge, 4 of the parameters are not fully estimated 199 # (possibly they are still at starting values?) so the AIC is off 200 pass 201 202 def test_bic(self): 203 # Stata uses 9 df (i.e. 9 params) here instead of 13, because since the 204 # model did not coverge, 4 of the parameters are not fully estimated 205 # (possibly they are still at starting values?) so the BIC is off 206 pass 207 208 def test_summary(self): 209 with warnings.catch_warnings(): 210 warnings.simplefilter("ignore") 211 summary = self.results.summary() 212 tables = [str(table) for table in summary.tables] 213 params = self.true['params'] 214 215 # Make sure we have the right number of tables 216 assert_equal( 217 len(tables), 218 2 + self.model.k_endog + self.model.k_factors + 1) 219 220 # Check the model overview table 221 assert re.search( 222 r'Model:.*DynamicFactor\(factors=2, order=1\)', 223 tables[0]) 224 225 # For each endogenous variable, check the output 226 for i in range(self.model.k_endog): 227 offset_loading = self.model.k_factors * i 228 table = tables[i + 2] 229 230 # -> Make sure we have the right table / table name 231 name = self.model.endog_names[i] 232 assert re.search('Results for equation %s' % name, table) 233 234 # -> Make sure it's the right size 235 assert_equal(len(table.split('\n')), 7) 236 237 # -> Check that we have the right coefficients 238 assert re.search( 239 'loading.f1 +' + forg(params[offset_loading + 0], prec=4), 240 table) 241 assert re.search( 242 'loading.f2 +' + forg(params[offset_loading + 1], prec=4), 243 table) 244 245 # For each factor, check the output 246 for i in range(self.model.k_factors): 247 offset = (self.model.k_endog * (self.model.k_factors + 1) + 248 i * self.model.k_factors) 249 table = tables[self.model.k_endog + i + 2] 250 251 # -> Make sure we have the right table / table name 252 name = self.model.endog_names[i] 253 assert re.search('Results for factor equation f%d' % (i+1), table) 254 255 # -> Make sure it's the right size 256 assert_equal(len(table.split('\n')), 7) 257 258 # -> Check that we have the right coefficients 259 assert re.search('L1.f1 +' + forg(params[offset + 0], prec=4), 260 table) 261 assert re.search('L1.f2 +' + forg(params[offset + 1], prec=4), 262 table) 263 264 # Check the Error covariance matrix output 265 table = tables[2 + self.model.k_endog + self.model.k_factors] 266 267 # -> Make sure we have the right table / table name 268 name = self.model.endog_names[i] 269 assert re.search('Error covariance matrix', table) 270 271 # -> Make sure it's the right size 272 assert_equal(len(table.split('\n')), 8) 273 274 # -> Check that we have the right coefficients 275 offset = self.model.k_endog * self.model.k_factors 276 for i in range(self.model.k_endog): 277 iname = self.model.endog_names[i] 278 iparam = forg(params[offset + i], prec=4) 279 assert re.search('sigma2.%s +%s' % (iname, iparam), table) 280 281 282class TestDynamicFactor_exog1(CheckDynamicFactor): 283 """ 284 Test for a dynamic factor model with 1 exogenous regressor: a constant 285 """ 286 @classmethod 287 def setup_class(cls): 288 true = results_dynamic_factor.lutkepohl_dfm_exog1.copy() 289 true['predict'] = output_results.iloc[1:][[ 290 'predict_dfm_exog1_1', 291 'predict_dfm_exog1_2', 292 'predict_dfm_exog1_3']] 293 true['dynamic_predict'] = output_results.iloc[1:][[ 294 'dyn_predict_dfm_exog1_1', 295 'dyn_predict_dfm_exog1_2', 296 'dyn_predict_dfm_exog1_3']] 297 exog = np.ones((75, 1)) 298 super(TestDynamicFactor_exog1, cls).setup_class( 299 true, k_factors=1, factor_order=1, exog=exog) 300 301 def test_predict(self): 302 exog = np.ones((16, 1)) 303 super(TestDynamicFactor_exog1, self).test_predict(exog=exog) 304 305 def test_dynamic_predict(self): 306 exog = np.ones((16, 1)) 307 super(TestDynamicFactor_exog1, self).test_dynamic_predict(exog=exog) 308 309 def test_bse_approx(self): 310 bse = self.results._cov_params_approx().diagonal()**0.5 311 assert_allclose(bse**2, self.true['var_oim'], atol=1e-5) 312 313 314class TestDynamicFactor_exog2(CheckDynamicFactor): 315 """ 316 Test for a dynamic factor model with 2 exogenous regressors: a constant 317 and a time-trend 318 """ 319 @classmethod 320 def setup_class(cls): 321 true = results_dynamic_factor.lutkepohl_dfm_exog2.copy() 322 true['predict'] = output_results.iloc[1:][[ 323 'predict_dfm_exog2_1', 324 'predict_dfm_exog2_2', 325 'predict_dfm_exog2_3']] 326 true['dynamic_predict'] = output_results.iloc[1:][[ 327 'dyn_predict_dfm_exog2_1', 328 'dyn_predict_dfm_exog2_2', 329 'dyn_predict_dfm_exog2_3']] 330 exog = np.c_[np.ones((75, 1)), (np.arange(75) + 2)[:, np.newaxis]] 331 super(TestDynamicFactor_exog2, cls).setup_class( 332 true, k_factors=1, factor_order=1, exog=exog) 333 334 def test_bse_approx(self): 335 bse = self.results._cov_params_approx().diagonal()**0.5 336 assert_allclose(bse**2, self.true['var_oim'], atol=1e-5) 337 338 def test_predict(self): 339 exog = np.c_[np.ones((16, 1)), 340 (np.arange(75, 75+16) + 2)[:, np.newaxis]] 341 super(TestDynamicFactor_exog2, self).test_predict(exog=exog) 342 343 def test_dynamic_predict(self): 344 exog = np.c_[np.ones((16, 1)), 345 (np.arange(75, 75+16) + 2)[:, np.newaxis]] 346 super(TestDynamicFactor_exog2, self).test_dynamic_predict(exog=exog) 347 348 def test_summary(self): 349 with warnings.catch_warnings(): 350 warnings.simplefilter("ignore") 351 summary = self.results.summary() 352 tables = [str(table) for table in summary.tables] 353 params = self.true['params'] 354 355 # Make sure we have the right number of tables 356 assert_equal( 357 len(tables), 358 2 + self.model.k_endog + self.model.k_factors + 1) 359 360 # Check the model overview table 361 assert re.search(r'Model:.*DynamicFactor\(factors=1, order=1\)', 362 tables[0]) 363 assert_equal(re.search(r'.*2 regressors', tables[0]) is None, False) 364 365 # For each endogenous variable, check the output 366 for i in range(self.model.k_endog): 367 offset_loading = self.model.k_factors * i 368 offset_exog = self.model.k_factors * self.model.k_endog 369 table = tables[i + 2] 370 371 # -> Make sure we have the right table / table name 372 name = self.model.endog_names[i] 373 assert re.search('Results for equation %s' % name, table) 374 375 # -> Make sure it's the right size 376 assert_equal(len(table.split('\n')), 8) 377 378 # -> Check that we have the right coefficients 379 assert re.search( 380 'loading.f1 +' + forg(params[offset_loading + 0], prec=4), 381 table) 382 assert re.search( 383 'beta.const +' + forg(params[offset_exog + i*2 + 0], prec=4), 384 table) 385 assert re.search( 386 'beta.x1 +' + forg(params[offset_exog + i*2 + 1], prec=4), 387 table) 388 389 # For each factor, check the output 390 for i in range(self.model.k_factors): 391 offset = (self.model.k_endog * (self.model.k_factors + 3) + 392 i * self.model.k_factors) 393 table = tables[self.model.k_endog + i + 2] 394 395 # -> Make sure we have the right table / table name 396 name = self.model.endog_names[i] 397 assert re.search('Results for factor equation f%d' % (i+1), table) 398 399 # -> Make sure it's the right size 400 assert_equal(len(table.split('\n')), 6) 401 402 # -> Check that we have the right coefficients 403 assert re.search('L1.f1 +' + forg(params[offset + 0], prec=4), 404 table) 405 406 # Check the Error covariance matrix output 407 table = tables[2 + self.model.k_endog + self.model.k_factors] 408 409 # -> Make sure we have the right table / table name 410 name = self.model.endog_names[i] 411 assert re.search('Error covariance matrix', table) 412 413 # -> Make sure it's the right size 414 assert_equal(len(table.split('\n')), 8) 415 416 # -> Check that we have the right coefficients 417 offset = self.model.k_endog * (self.model.k_factors + 2) 418 for i in range(self.model.k_endog): 419 iname = self.model.endog_names[i] 420 iparam = forg(params[offset + i], prec=4) 421 assert re.search('sigma2.%s +%s' % (iname, iparam), table) 422 423 424class TestDynamicFactor_general_errors(CheckDynamicFactor): 425 """ 426 Test for a dynamic factor model where errors are as general as possible, 427 meaning: 428 429 - Errors are vector autocorrelated, VAR(1) 430 - Innovations are correlated 431 """ 432 @classmethod 433 def setup_class(cls): 434 true = results_dynamic_factor.lutkepohl_dfm_gen.copy() 435 true['predict'] = output_results.iloc[1:][[ 436 'predict_dfm_gen_1', 'predict_dfm_gen_2', 'predict_dfm_gen_3']] 437 true['dynamic_predict'] = output_results.iloc[1:][[ 438 'dyn_predict_dfm_gen_1', 439 'dyn_predict_dfm_gen_2', 440 'dyn_predict_dfm_gen_3']] 441 super(TestDynamicFactor_general_errors, cls).setup_class( 442 true, k_factors=1, factor_order=1, error_var=True, 443 error_order=1, error_cov_type='unstructured') 444 445 def test_bse_approx(self): 446 bse = self.results._cov_params_approx().diagonal() 447 assert_allclose(bse[:3], self.true['var_oim'][:3], atol=1e-5) 448 assert_allclose(bse[-10:], self.true['var_oim'][-10:], atol=3e-4) 449 450 @pytest.mark.skip("Known failure, no sequence of optimizers has been " 451 "found which can achieve the maximum.") 452 def test_mle(self): 453 # The following gets us to llf=546.53, which is still not good enough 454 # llf = 300.842477412 455 # res = mod.fit(method='lbfgs', maxiter=10000) 456 # llf = 460.26576722 457 # res = mod.fit(res.params, method='nm', maxiter=10000, maxfev=10000) 458 # llf = 542.245718508 459 # res = mod.fit(res.params, method='lbfgs', maxiter=10000) 460 # llf = 544.035160955 461 # res = mod.fit(res.params, method='nm', maxiter=10000, maxfev=10000) 462 # llf = 557.442240083 463 # res = mod.fit(res.params, method='lbfgs', maxiter=10000) 464 # llf = 558.199513262 465 # res = mod.fit(res.params, method='nm', maxiter=10000, maxfev=10000) 466 # llf = 559.049076604 467 # res = mod.fit(res.params, method='nm', maxiter=10000, maxfev=10000) 468 # llf = 559.049076604 469 # ... 470 pass 471 472 def test_summary(self): 473 with warnings.catch_warnings(): 474 warnings.simplefilter("ignore") 475 summary = self.results.summary() 476 tables = [str(table) for table in summary.tables] 477 params = self.true['params'] 478 479 # Make sure we have the right number of tables 480 assert_equal( 481 len(tables), 482 2 + self.model.k_endog + self.model.k_factors + 483 self.model.k_endog + 1) 484 485 # Check the model overview table 486 assert re.search(r'Model:.*DynamicFactor\(factors=1, order=1\)', 487 tables[0]) 488 assert re.search(r'.*VAR\(1\) errors', tables[0]) 489 490 # For each endogenous variable, check the output 491 for i in range(self.model.k_endog): 492 offset_loading = self.model.k_factors * i 493 table = tables[i + 2] 494 495 # -> Make sure we have the right table / table name 496 name = self.model.endog_names[i] 497 assert re.search('Results for equation %s' % name, table) 498 499 # -> Make sure it's the right size 500 assert_equal(len(table.split('\n')), 6) 501 502 # -> Check that we have the right coefficients 503 pattern = 'loading.f1 +' + forg(params[offset_loading + 0], prec=4) 504 assert re.search(pattern, table) 505 506 # For each factor, check the output 507 for i in range(self.model.k_factors): 508 offset = (self.model.k_endog * self.model.k_factors + 509 6 + i * self.model.k_factors) 510 table = tables[2 + self.model.k_endog + i] 511 512 # -> Make sure we have the right table / table name 513 name = self.model.endog_names[i] 514 assert re.search('Results for factor equation f%d' % (i+1), table) 515 516 # -> Make sure it's the right size 517 assert_equal(len(table.split('\n')), 6) 518 519 # -> Check that we have the right coefficients 520 assert re.search('L1.f1 +' + forg(params[offset + 0], prec=4), 521 table) 522 523 # For each error equation, check the output 524 for i in range(self.model.k_endog): 525 offset = (self.model.k_endog * (self.model.k_factors + i) + 526 6 + self.model.k_factors) 527 table = tables[2 + self.model.k_endog + self.model.k_factors + i] 528 529 # -> Make sure we have the right table / table name 530 name = self.model.endog_names[i] 531 assert re.search(r'Results for error equation e\(%s\)' % name, 532 table) 533 534 # -> Make sure it's the right size 535 assert_equal(len(table.split('\n')), 8) 536 537 # -> Check that we have the right coefficients 538 for j in range(self.model.k_endog): 539 name = self.model.endog_names[j] 540 pattern = r'L1.e\(%s\) +%s' % (name, forg(params[offset + j], 541 prec=4)) 542 assert re.search(pattern, table) 543 544 # Check the Error covariance matrix output 545 table = tables[2 + self.model.k_endog + 546 self.model.k_factors + self.model.k_endog] 547 548 # -> Make sure we have the right table / table name 549 name = self.model.endog_names[i] 550 assert re.search('Error covariance matrix', table) 551 552 # -> Make sure it's the right size 553 assert_equal(len(table.split('\n')), 11) 554 555 # -> Check that we have the right coefficients 556 offset = self.model.k_endog * self.model.k_factors 557 assert re.search( 558 r'cov.chol\[1,1\] +' + forg(params[offset + 0], prec=4), 559 table) 560 assert re.search( 561 r'cov.chol\[2,1\] +' + forg(params[offset + 1], prec=4), 562 table) 563 assert re.search( 564 r'cov.chol\[2,2\] +' + forg(params[offset + 2], prec=4), 565 table) 566 assert re.search( 567 r'cov.chol\[3,1\] +' + forg(params[offset+3], prec=4), 568 table) 569 assert re.search( 570 r'cov.chol\[3,2\] +' + forg(params[offset+4], prec=4), 571 table) 572 assert re.search( 573 r'cov.chol\[3,3\] +' + forg(params[offset + 5], prec=4), 574 table) 575 576 577class TestDynamicFactor_ar2_errors(CheckDynamicFactor): 578 """ 579 Test for a dynamic factor model where errors are as general as possible, 580 meaning: 581 582 - Errors are vector autocorrelated, VAR(1) 583 - Innovations are correlated 584 """ 585 @classmethod 586 def setup_class(cls): 587 true = results_dynamic_factor.lutkepohl_dfm_ar2.copy() 588 true['predict'] = output_results.iloc[1:][[ 589 'predict_dfm_ar2_1', 'predict_dfm_ar2_2', 'predict_dfm_ar2_3']] 590 true['dynamic_predict'] = output_results.iloc[1:][[ 591 'dyn_predict_dfm_ar2_1', 592 'dyn_predict_dfm_ar2_2', 593 'dyn_predict_dfm_ar2_3']] 594 super(TestDynamicFactor_ar2_errors, cls).setup_class( 595 true, k_factors=1, factor_order=1, error_order=2) 596 597 def test_bse_approx(self): 598 bse = self.results._cov_params_approx().diagonal() 599 assert_allclose(bse, self.true['var_oim'], atol=1e-5) 600 601 def test_mle(self): 602 with warnings.catch_warnings(record=True): 603 # Depending on the system, this test can reach a greater precision, 604 # but for cross-platform results keep it at 1e-2 605 mod = self.model 606 res1 = mod.fit(maxiter=100, optim_score='approx', disp=False) 607 res = mod.fit( 608 res1.params, method='nm', maxiter=10000, 609 optim_score='approx', disp=False) 610 # Added rtol to catch spurious failures on some platforms 611 assert_allclose(res.llf, self.results.llf, atol=1e-2, rtol=1e-4) 612 613 614class TestDynamicFactor_scalar_error(CheckDynamicFactor): 615 """ 616 Test for a dynamic factor model where innovations are uncorrelated and 617 are forced to have the same variance. 618 """ 619 @classmethod 620 def setup_class(cls): 621 true = results_dynamic_factor.lutkepohl_dfm_scalar.copy() 622 true['predict'] = output_results.iloc[1:][[ 623 'predict_dfm_scalar_1', 'predict_dfm_scalar_2', 624 'predict_dfm_scalar_3']] 625 true['dynamic_predict'] = output_results.iloc[1:][[ 626 'dyn_predict_dfm_scalar_1', 'dyn_predict_dfm_scalar_2', 627 'dyn_predict_dfm_scalar_3']] 628 exog = np.ones((75, 1)) 629 super(TestDynamicFactor_scalar_error, cls).setup_class( 630 true, k_factors=1, factor_order=1, 631 exog=exog, error_cov_type='scalar') 632 633 def test_bse_approx(self): 634 bse = self.results._cov_params_approx().diagonal() 635 assert_allclose(bse, self.true['var_oim'], atol=1e-5) 636 637 def test_predict(self): 638 exog = np.ones((16, 1)) 639 super(TestDynamicFactor_scalar_error, self).test_predict(exog=exog) 640 641 def test_dynamic_predict(self): 642 exog = np.ones((16, 1)) 643 super(TestDynamicFactor_scalar_error, 644 self).test_dynamic_predict(exog=exog) 645 646 647class TestStaticFactor(CheckDynamicFactor): 648 """ 649 Test for a static factor model (i.e. factors are not autocorrelated). 650 """ 651 @classmethod 652 def setup_class(cls): 653 true = results_dynamic_factor.lutkepohl_sfm.copy() 654 true['predict'] = output_results.iloc[1:][[ 655 'predict_sfm_1', 'predict_sfm_2', 'predict_sfm_3']] 656 true['dynamic_predict'] = output_results.iloc[1:][[ 657 'dyn_predict_sfm_1', 'dyn_predict_sfm_2', 'dyn_predict_sfm_3']] 658 super(TestStaticFactor, cls).setup_class( 659 true, k_factors=1, factor_order=0) 660 661 def test_bse_approx(self): 662 bse = self.results._cov_params_approx().diagonal() 663 assert_allclose(bse, self.true['var_oim'], atol=1e-5) 664 665 def test_bic(self): 666 # Stata uses 5 df (i.e. 5 params) here instead of 6, because one param 667 # is basically zero. 668 pass 669 670 671class TestSUR(CheckDynamicFactor): 672 """ 673 Test for a seemingly unrelated regression model (i.e. no factors) with 674 errors cross-sectionally, but not auto-, correlated 675 """ 676 @classmethod 677 def setup_class(cls): 678 true = results_dynamic_factor.lutkepohl_sur.copy() 679 true['predict'] = output_results.iloc[1:][[ 680 'predict_sur_1', 'predict_sur_2', 'predict_sur_3']] 681 true['dynamic_predict'] = output_results.iloc[1:][[ 682 'dyn_predict_sur_1', 'dyn_predict_sur_2', 'dyn_predict_sur_3']] 683 exog = np.c_[np.ones((75, 1)), (np.arange(75) + 2)[:, np.newaxis]] 684 super(TestSUR, cls).setup_class( 685 true, k_factors=0, factor_order=0, 686 exog=exog, error_cov_type='unstructured') 687 688 def test_bse_approx(self): 689 bse = self.results._cov_params_approx().diagonal() 690 assert_allclose(bse[:6], self.true['var_oim'][:6], atol=1e-5) 691 692 def test_predict(self): 693 exog = np.c_[np.ones((16, 1)), 694 (np.arange(75, 75+16) + 2)[:, np.newaxis]] 695 super(TestSUR, self).test_predict(exog=exog) 696 697 def test_dynamic_predict(self): 698 exog = np.c_[np.ones((16, 1)), 699 (np.arange(75, 75+16) + 2)[:, np.newaxis]] 700 super(TestSUR, self).test_dynamic_predict(exog=exog) 701 702 703class TestSUR_autocorrelated_errors(CheckDynamicFactor): 704 """ 705 Test for a seemingly unrelated regression model (i.e. no factors) where 706 the errors are vector autocorrelated, but innovations are uncorrelated. 707 """ 708 @classmethod 709 def setup_class(cls): 710 true = results_dynamic_factor.lutkepohl_sur_auto.copy() 711 true['predict'] = output_results.iloc[1:][[ 712 'predict_sur_auto_1', 'predict_sur_auto_2']] 713 true['dynamic_predict'] = output_results.iloc[1:][[ 714 'dyn_predict_sur_auto_1', 'dyn_predict_sur_auto_2']] 715 exog = np.c_[np.ones((75, 1)), (np.arange(75) + 2)[:, np.newaxis]] 716 super(TestSUR_autocorrelated_errors, cls).setup_class( 717 true, k_factors=0, factor_order=0, exog=exog, 718 error_order=1, error_var=True, 719 error_cov_type='diagonal', 720 included_vars=['dln_inv', 'dln_inc']) 721 722 def test_bse_approx(self): 723 bse = self.results._cov_params_approx().diagonal() 724 assert_allclose(bse, self.true['var_oim'], atol=1e-5) 725 726 def test_predict(self): 727 exog = np.c_[np.ones((16, 1)), 728 (np.arange(75, 75+16) + 2)[:, np.newaxis]] 729 super(TestSUR_autocorrelated_errors, self).test_predict(exog=exog) 730 731 def test_dynamic_predict(self): 732 exog = np.c_[np.ones((16, 1)), 733 (np.arange(75, 75+16) + 2)[:, np.newaxis]] 734 super(TestSUR_autocorrelated_errors, 735 self).test_dynamic_predict(exog=exog) 736 737 def test_mle(self): 738 super(TestSUR_autocorrelated_errors, self).test_mle(init_powell=False) 739 740 741def test_misspecification(): 742 # Tests for model specification and misspecification exceptions 743 endog = np.arange(20).reshape(10, 2) 744 745 # Too few endog 746 assert_raises( 747 ValueError, 748 dynamic_factor.DynamicFactor, endog[:, 0], k_factors=0, factor_order=0) 749 750 # Too many factors 751 assert_raises( 752 ValueError, 753 dynamic_factor.DynamicFactor, endog, k_factors=2, factor_order=1) 754 755 # Bad error_cov_type specification 756 assert_raises( 757 ValueError, 758 dynamic_factor.DynamicFactor, 759 endog, 760 k_factors=1, factor_order=1, order=(1, 0), error_cov_type='') 761 762 763def test_miscellaneous(): 764 # Initialization with 1-dimensional exog array 765 exog = np.arange(75) 766 mod = CheckDynamicFactor() 767 mod.setup_class(true=None, k_factors=1, factor_order=1, 768 exog=exog, filter=False) 769 exog = pd.Series(np.arange(75), 770 index=pd.date_range(start='1960-04-01', 771 end='1978-10-01', freq='QS')) 772 mod = CheckDynamicFactor() 773 mod.setup_class( 774 true=None, k_factors=1, factor_order=1, exog=exog, filter=False) 775 776 777def test_predict_custom_index(): 778 np.random.seed(328423) 779 endog = pd.DataFrame(np.random.normal(size=(50, 2))) 780 mod = dynamic_factor.DynamicFactor(endog, k_factors=1, factor_order=1) 781 res = mod.smooth(mod.start_params) 782 out = res.predict(start=1, end=1, index=['a']) 783 assert_equal(out.index.equals(pd.Index(['a'])), True) 784 785 786def test_forecast_exog(): 787 # Test forecasting with various shapes of `exog` 788 nobs = 100 789 endog = np.ones((nobs, 2)) * 2.0 790 exog = np.ones(nobs) 791 792 mod = dynamic_factor.DynamicFactor(endog, exog=exog, k_factors=1, 793 factor_order=1) 794 res = mod.smooth(np.r_[[0] * 2, 2.0, 2.0, 1, 1., 0.]) 795 796 # 1-step-ahead, valid 797 exog_fcast_scalar = 1. 798 exog_fcast_1dim = np.ones(1) 799 exog_fcast_2dim = np.ones((1, 1)) 800 801 assert_allclose(res.forecast(1, exog=exog_fcast_scalar), 2.) 802 assert_allclose(res.forecast(1, exog=exog_fcast_1dim), 2.) 803 assert_allclose(res.forecast(1, exog=exog_fcast_2dim), 2.) 804 805 # h-steps-ahead, valid 806 h = 10 807 exog_fcast_1dim = np.ones(h) 808 exog_fcast_2dim = np.ones((h, 1)) 809 810 assert_allclose(res.forecast(h, exog=exog_fcast_1dim), 2.) 811 assert_allclose(res.forecast(h, exog=exog_fcast_2dim), 2.) 812 813 # h-steps-ahead, invalid 814 assert_raises(ValueError, res.forecast, h, exog=1.) 815 assert_raises(ValueError, res.forecast, h, exog=[1, 2]) 816 assert_raises(ValueError, res.forecast, h, exog=np.ones((h, 2))) 817 818 819def check_equivalent_models(mod, mod2): 820 attrs = [ 821 'k_factors', 'factor_order', 'error_order', 'error_var', 822 'error_cov_type', 'enforce_stationarity', 'mle_regression', 'k_params'] 823 824 ssm_attrs = [ 825 'nobs', 'k_endog', 'k_states', 'k_posdef', 'obs_intercept', 'design', 826 'obs_cov', 'state_intercept', 'transition', 'selection', 'state_cov'] 827 828 for attr in attrs: 829 assert_equal(getattr(mod2, attr), getattr(mod, attr)) 830 831 for attr in ssm_attrs: 832 assert_equal(getattr(mod2.ssm, attr), getattr(mod.ssm, attr)) 833 834 assert_equal(mod2._get_init_kwds(), mod._get_init_kwds()) 835 836 837def test_recreate_model(): 838 nobs = 100 839 endog = np.ones((nobs, 3)) * 2.0 840 exog = np.ones(nobs) 841 842 k_factors = [0, 1, 2] 843 factor_orders = [0, 1, 2] 844 error_orders = [0, 1] 845 error_vars = [False, True] 846 error_cov_types = ['diagonal', 'scalar'] 847 848 import itertools 849 names = ['k_factors', 'factor_order', 'error_order', 'error_var', 850 'error_cov_type'] 851 for element in itertools.product(k_factors, factor_orders, error_orders, 852 error_vars, error_cov_types): 853 kwargs = dict(zip(names, element)) 854 855 mod = dynamic_factor.DynamicFactor(endog, exog=exog, **kwargs) 856 mod2 = dynamic_factor.DynamicFactor(endog, exog=exog, 857 **mod._get_init_kwds()) 858 check_equivalent_models(mod, mod2) 859 860 861def test_append_results(): 862 endog = np.arange(200).reshape(100, 2) 863 exog = np.ones(100) 864 params = [0.1, -0.2, 1., 2., 1., 1., 0.5, 0.1] 865 866 mod1 = dynamic_factor.DynamicFactor( 867 endog, k_factors=1, factor_order=2, exog=exog) 868 res1 = mod1.smooth(params) 869 870 mod2 = dynamic_factor.DynamicFactor( 871 endog[:50], k_factors=1, factor_order=2, exog=exog[:50]) 872 res2 = mod2.smooth(params) 873 res3 = res2.append(endog[50:], exog=exog[50:]) 874 875 assert_equal(res1.specification, res3.specification) 876 877 assert_allclose(res3.cov_params_default, res2.cov_params_default) 878 for attr in ['nobs', 'llf', 'llf_obs', 'loglikelihood_burn']: 879 assert_equal(getattr(res3, attr), getattr(res1, attr)) 880 881 for attr in [ 882 'filtered_state', 'filtered_state_cov', 'predicted_state', 883 'predicted_state_cov', 'forecasts', 'forecasts_error', 884 'forecasts_error_cov', 'standardized_forecasts_error', 885 'forecasts_error_diffuse_cov', 'predicted_diffuse_state_cov', 886 'scaled_smoothed_estimator', 887 'scaled_smoothed_estimator_cov', 'smoothing_error', 888 'smoothed_state', 889 'smoothed_state_cov', 'smoothed_state_autocov', 890 'smoothed_measurement_disturbance', 891 'smoothed_state_disturbance', 892 'smoothed_measurement_disturbance_cov', 893 'smoothed_state_disturbance_cov']: 894 assert_equal(getattr(res3, attr), getattr(res1, attr)) 895 896 assert_allclose(res3.forecast(10, exog=np.ones(10)), 897 res1.forecast(10, exog=np.ones(10))) 898 899 900def test_extend_results(): 901 endog = np.arange(200).reshape(100, 2) 902 exog = np.ones(100) 903 params = [0.1, -0.2, 1., 2., 1., 1., 0.5, 0.1] 904 905 mod1 = dynamic_factor.DynamicFactor( 906 endog, k_factors=1, factor_order=2, exog=exog) 907 res1 = mod1.smooth(params) 908 909 mod2 = dynamic_factor.DynamicFactor( 910 endog[:50], k_factors=1, factor_order=2, exog=exog[:50]) 911 res2 = mod2.smooth(params) 912 res3 = res2.extend(endog[50:], exog=exog[50:]) 913 914 assert_allclose(res3.llf_obs, res1.llf_obs[50:]) 915 916 for attr in [ 917 'filtered_state', 'filtered_state_cov', 'predicted_state', 918 'predicted_state_cov', 'forecasts', 'forecasts_error', 919 'forecasts_error_cov', 'standardized_forecasts_error', 920 'forecasts_error_diffuse_cov', 'predicted_diffuse_state_cov', 921 'scaled_smoothed_estimator', 922 'scaled_smoothed_estimator_cov', 'smoothing_error', 923 'smoothed_state', 924 'smoothed_state_cov', 'smoothed_state_autocov', 925 'smoothed_measurement_disturbance', 926 'smoothed_state_disturbance', 927 'smoothed_measurement_disturbance_cov', 928 'smoothed_state_disturbance_cov']: 929 desired = getattr(res1, attr) 930 if desired is not None: 931 desired = desired[..., 50:] 932 assert_equal(getattr(res3, attr), desired) 933 934 assert_allclose(res3.forecast(10, exog=np.ones(10)), 935 res1.forecast(10, exog=np.ones(10))) 936 937 938def test_apply_results(): 939 endog = np.arange(200).reshape(100, 2) 940 exog = np.ones(100) 941 params = [0.1, -0.2, 1., 2., 1., 1., 0.5, 0.1] 942 943 mod1 = dynamic_factor.DynamicFactor( 944 endog[:50], k_factors=1, factor_order=2, exog=exog[:50]) 945 res1 = mod1.smooth(params) 946 947 mod2 = dynamic_factor.DynamicFactor( 948 endog[50:], k_factors=1, factor_order=2, exog=exog[50:]) 949 res2 = mod2.smooth(params) 950 951 res3 = res2.apply(endog[:50], exog=exog[:50]) 952 953 assert_equal(res1.specification, res3.specification) 954 955 assert_allclose(res3.cov_params_default, res2.cov_params_default) 956 for attr in ['nobs', 'llf', 'llf_obs', 'loglikelihood_burn']: 957 assert_equal(getattr(res3, attr), getattr(res1, attr)) 958 959 for attr in [ 960 'filtered_state', 'filtered_state_cov', 'predicted_state', 961 'predicted_state_cov', 'forecasts', 'forecasts_error', 962 'forecasts_error_cov', 'standardized_forecasts_error', 963 'forecasts_error_diffuse_cov', 'predicted_diffuse_state_cov', 964 'scaled_smoothed_estimator', 965 'scaled_smoothed_estimator_cov', 'smoothing_error', 966 'smoothed_state', 967 'smoothed_state_cov', 'smoothed_state_autocov', 968 'smoothed_measurement_disturbance', 969 'smoothed_state_disturbance', 970 'smoothed_measurement_disturbance_cov', 971 'smoothed_state_disturbance_cov']: 972 assert_equal(getattr(res3, attr), getattr(res1, attr)) 973 974 assert_allclose(res3.forecast(10, exog=np.ones(10)), 975 res1.forecast(10, exog=np.ones(10))) 976 977 978def test_start_params_nans(): 979 ix = pd.date_range('1960-01-01', '1982-10-01', freq='QS') 980 dta = np.log(pd.DataFrame( 981 results_varmax.lutkepohl_data, columns=['inv', 'inc', 'consump'], 982 index=ix)).diff().iloc[1:] 983 984 endog1 = dta.iloc[:-1] 985 mod1 = dynamic_factor.DynamicFactor(endog1, k_factors=1, factor_order=1) 986 endog2 = dta.copy() 987 endog2.iloc[-1:] = np.nan 988 mod2 = dynamic_factor.DynamicFactor(endog2, k_factors=1, factor_order=1) 989 990 assert_allclose(mod2.start_params, mod1.start_params) 991