1# -*- coding: utf-8 -*- 2"""Test for short_panel and panel sandwich 3 4Created on Fri May 18 13:05:47 2012 5 6Author: Josef Perktold 7 8moved example from main of random_panel 9""" 10 11import numpy as np 12from numpy.testing import assert_almost_equal 13import numpy.testing as npt 14import statsmodels.tools.eval_measures as em 15from statsmodels.stats.moment_helpers import cov2corr, se_cov 16from statsmodels.regression.linear_model import OLS 17 18from statsmodels.sandbox.panel.panel_short import ShortPanelGLS, ShortPanelGLS2 19from statsmodels.sandbox.panel.random_panel import PanelSample 20import statsmodels.sandbox.panel.correlation_structures as cs 21import statsmodels.stats.sandwich_covariance as sw 22 23def assert_maxabs(actual, expected, value): 24 npt.assert_array_less(em.maxabs(actual, expected, None), value) 25 26 27def test_short_panel(): 28 #this checks that some basic statistical properties are satisfied by the 29 #results, not verified results against other packages 30 #Note: the ranking of robust bse is different if within=True 31 #I added within keyword to PanelSample to be able to use old example 32 #if within is False, then there is no within group variation in exog. 33 nobs = 100 34 nobs_i = 5 35 n_groups = nobs // nobs_i 36 k_vars = 3 37 38 dgp = PanelSample(nobs, k_vars, n_groups, corr_structure=cs.corr_arma, 39 corr_args=([1], [1., -0.9],), seed=377769, within=False) 40 #print 'seed', dgp.seed 41 y = dgp.generate_panel() 42 noise = y - dgp.y_true 43 44 #test dgp 45 46 dgp_cov_e = np.array( 47 [[ 1. , 0.9 , 0.81 , 0.729 , 0.6561], 48 [ 0.9 , 1. , 0.9 , 0.81 , 0.729 ], 49 [ 0.81 , 0.9 , 1. , 0.9 , 0.81 ], 50 [ 0.729 , 0.81 , 0.9 , 1. , 0.9 ], 51 [ 0.6561, 0.729 , 0.81 , 0.9 , 1. ]]) 52 53 npt.assert_almost_equal(dgp.cov, dgp_cov_e, 13) 54 55 cov_noise = np.cov(noise.reshape(-1,n_groups, order='F')) 56 corr_noise = cov2corr(cov_noise) 57 npt.assert_almost_equal(corr_noise, dgp.cov, 1) 58 59 #estimate panel model 60 mod2 = ShortPanelGLS(y, dgp.exog, dgp.groups) 61 res2 = mod2.fit_iterative(2) 62 63 64 #whitened residual should be uncorrelated 65 corr_wresid = np.corrcoef(res2.wresid.reshape(-1,n_groups, order='F')) 66 assert_maxabs(corr_wresid, np.eye(5), 0.1) 67 68 #residual should have same correlation as dgp 69 corr_resid = np.corrcoef(res2.resid.reshape(-1,n_groups, order='F')) 70 assert_maxabs(corr_resid, dgp.cov, 0.1) 71 72 assert_almost_equal(res2.resid.std(),1, decimal=0) 73 74 y_pred = np.dot(mod2.exog, res2.params) 75 assert_almost_equal(res2.fittedvalues, y_pred, 13) 76 77 78 #compare with OLS 79 80 res2_ols = mod2._fit_ols() 81 npt.assert_(mod2.res_pooled is res2_ols) 82 83 res2_ols = mod2.res_pooled #TODO: BUG: requires call to _fit_ols 84 85 #fitting once is the same as OLS 86 #note: I need to create new instance, otherwise it continuous fitting 87 mod1 = ShortPanelGLS(y, dgp.exog, dgp.groups) 88 res1 = mod1.fit_iterative(1) 89 90 assert_almost_equal(res1.params, res2_ols.params, decimal=13) 91 assert_almost_equal(res1.bse, res2_ols.bse, decimal=13) 92 93 res_ols = OLS(y, dgp.exog).fit() 94 assert_almost_equal(res1.params, res_ols.params, decimal=13) 95 assert_almost_equal(res1.bse, res_ols.bse, decimal=13) 96 97 98 #compare with old version 99 mod_old = ShortPanelGLS2(y, dgp.exog, dgp.groups) 100 res_old = mod_old.fit() 101 102 assert_almost_equal(res2.params, res_old.params, decimal=13) 103 assert_almost_equal(res2.bse, res_old.bse, decimal=13) 104 105 106 mod5 = ShortPanelGLS(y, dgp.exog, dgp.groups) 107 res5 = mod5.fit_iterative(5) 108 109 #make sure it's different 110 #npt.assert_array_less(0.009, em.maxabs(res5.bse, res2.bse)) 111 112 cov_clu = sw.cov_cluster(mod2.res_pooled, dgp.groups.astype(int)) 113 clubse = se_cov(cov_clu) 114 pnwbse = se_cov(sw.cov_nw_panel(mod2.res_pooled, 4, mod2.group.groupidx)) 115 bser = np.vstack((res2.bse, res5.bse, clubse, pnwbse)) 116 bser_mean = np.mean(bser, axis=0) 117 118 #cov_cluster close to robust and PanelGLS 119 #is up to 24% larger than mean of bser 120 #npt.assert_array_less(0, clubse / bser_mean - 1) 121 npt.assert_array_less(clubse / bser_mean - 1, 0.25) 122 #cov_nw_panel close to robust and PanelGLS 123 npt.assert_array_less(pnwbse / bser_mean - 1, 0.1) 124 #OLS underestimates bse, robust at least 60% larger 125 npt.assert_array_less(0.6, bser_mean / res_ols.bse - 1) 126 127 #cov_hac_panel with uniform_kernel is the same as cov_cluster for balanced 128 #panel with full length kernel 129 #I fixe default correction to be equal 130 cov_uni = sw.cov_nw_panel(mod2.res_pooled, 4, mod2.group.groupidx, 131 weights_func=sw.weights_uniform, 132 use_correction='c') 133 assert_almost_equal(cov_uni, cov_clu, decimal=13) 134 135 #without correction 136 cov_clu2 = sw.cov_cluster(mod2.res_pooled, dgp.groups.astype(int), 137 use_correction=False) 138 cov_uni2 = sw.cov_nw_panel(mod2.res_pooled, 4, mod2.group.groupidx, 139 weights_func=sw.weights_uniform, 140 use_correction=False) 141 assert_almost_equal(cov_uni2, cov_clu2, decimal=13) 142 143 cov_white = sw.cov_white_simple(mod2.res_pooled) 144 cov_pnw0 = sw.cov_nw_panel(mod2.res_pooled, 0, mod2.group.groupidx, 145 use_correction='hac') 146 assert_almost_equal(cov_pnw0, cov_white, decimal=13) 147