1# -*- coding: utf-8 -*- 2""" 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 statsmodels.sandbox.panel.panel_short import ShortPanelGLS, ShortPanelGLS2 13from statsmodels.sandbox.panel.random_panel import PanelSample 14import statsmodels.sandbox.panel.correlation_structures as cs 15 16import statsmodels.stats.sandwich_covariance as sw 17#from statsmodels.stats.sandwich_covariance import ( 18# S_hac_groupsum, weights_bartlett, _HCCM2) 19from statsmodels.stats.moment_helpers import se_cov 20cov_nw_panel2 = sw.cov_nw_groupsum 21 22 23examples = ['ex1'] 24 25 26if 'ex1' in examples: 27 nobs = 100 28 nobs_i = 5 29 n_groups = nobs // nobs_i 30 k_vars = 3 31 32# dgp = PanelSample(nobs, k_vars, n_groups, corr_structure=cs.corr_equi, 33# corr_args=(0.6,)) 34# dgp = PanelSample(nobs, k_vars, n_groups, corr_structure=cs.corr_ar, 35# corr_args=([1, -0.95],)) 36 dgp = PanelSample(nobs, k_vars, n_groups, corr_structure=cs.corr_arma, 37 corr_args=([1], [1., -0.9],), seed=377769) 38 print('seed', dgp.seed) 39 y = dgp.generate_panel() 40 noise = y - dgp.y_true 41 print(np.corrcoef(y.reshape(-1,n_groups, order='F'))) 42 print(np.corrcoef(noise.reshape(-1,n_groups, order='F'))) 43 44 mod = ShortPanelGLS2(y, dgp.exog, dgp.groups) 45 res = mod.fit() 46 print(res.params) 47 print(res.bse) 48 #Now what? 49 #res.resid is of transformed model 50 #np.corrcoef(res.resid.reshape(-1,n_groups, order='F')) 51 y_pred = np.dot(mod.exog, res.params) 52 resid = y - y_pred 53 print(np.corrcoef(resid.reshape(-1,n_groups, order='F'))) 54 print(resid.std()) 55 err = y_pred - dgp.y_true 56 print(err.std()) 57 #OLS standard errors are too small 58 mod.res_pooled.params 59 mod.res_pooled.bse 60 #heteroscedasticity robust does not help 61 mod.res_pooled.HC1_se 62 #compare with cluster robust se 63 64 print(sw.se_cov(sw.cov_cluster(mod.res_pooled, dgp.groups.astype(int)))) 65 #not bad, pretty close to panel estimator 66 #and with Newey-West Hac 67 print(sw.se_cov(sw.cov_nw_panel(mod.res_pooled, 4, mod.group.groupidx))) 68 #too small, assuming no bugs, 69 #see Peterson assuming it refers to same kind of model 70 print(dgp.cov) 71 72 mod2 = ShortPanelGLS(y, dgp.exog, dgp.groups) 73 res2 = mod2.fit_iterative(2) 74 print(res2.params) 75 print(res2.bse) 76 #both implementations produce the same results: 77 from numpy.testing import assert_almost_equal 78 assert_almost_equal(res.params, res2.params, decimal=12) 79 assert_almost_equal(res.bse, res2.bse, decimal=13) 80 mod5 = ShortPanelGLS(y, dgp.exog, dgp.groups) 81 res5 = mod5.fit_iterative(5) 82 print(res5.params) 83 print(res5.bse) 84 #fitting once is the same as OLS 85 #note: I need to create new instance, otherwise it continuous fitting 86 mod1 = ShortPanelGLS(y, dgp.exog, dgp.groups) 87 res1 = mod1.fit_iterative(1) 88 res_ols = mod1._fit_ols() 89 assert_almost_equal(res1.params, res_ols.params, decimal=12) 90 assert_almost_equal(res1.bse, res_ols.bse, decimal=13) 91 92 #cov_hac_panel with uniform_kernel is the same as cov_cluster for balanced 93 #panel with full length kernel 94 #I fixe default correction to be equal 95 mod2._fit_ols() 96 cov_clu = sw.cov_cluster(mod2.res_pooled, dgp.groups.astype(int)) 97 clubse = se_cov(cov_clu) 98 cov_uni = sw.cov_nw_panel(mod2.res_pooled, 4, mod2.group.groupidx, 99 weights_func=sw.weights_uniform, 100 use_correction='cluster') 101 assert_almost_equal(cov_uni, cov_clu, decimal=7) 102 103 #without correction 104 cov_clu2 = sw.cov_cluster(mod2.res_pooled, dgp.groups.astype(int), 105 use_correction=False) 106 cov_uni2 = sw.cov_nw_panel(mod2.res_pooled, 4, mod2.group.groupidx, 107 weights_func=sw.weights_uniform, 108 use_correction=False) 109 assert_almost_equal(cov_uni2, cov_clu2, decimal=8) 110 111 cov_white = sw.cov_white_simple(mod2.res_pooled) 112 cov_pnw0 = sw.cov_nw_panel(mod2.res_pooled, 0, mod2.group.groupidx, 113 use_correction='hac') 114 assert_almost_equal(cov_pnw0, cov_white, decimal=13) 115 116 time = np.tile(np.arange(nobs_i), n_groups) 117 #time = mod2.group.group_int 118 cov_pnw1 = sw.cov_nw_panel(mod2.res_pooled, 4, mod2.group.groupidx) 119 cov_pnw2 = cov_nw_panel2(mod2.res_pooled, 4, time) 120 #s = sw.group_sums(x, time) 121 122 c2, ct, cg = sw.cov_cluster_2groups(mod2.res_pooled, time, dgp.groups.astype(int), use_correction=False) 123 ct_nw0 = cov_nw_panel2(mod2.res_pooled, 0, time, weights_func=sw.weights_uniform, use_correction=False) 124 cg_nw0 = cov_nw_panel2(mod2.res_pooled, 0, dgp.groups.astype(int), weights_func=sw.weights_uniform, use_correction=False) 125 assert_almost_equal(ct_nw0, ct, decimal=13) 126 assert_almost_equal(cg_nw0, cg, decimal=13) #pnw2 0 lags 127 assert_almost_equal(cov_clu2, cg, decimal=13) 128 assert_almost_equal(cov_uni2, cg, decimal=8) #pnw all lags 129 130 131 132 133 import pandas as pa 134 #pandas.DataFrame does not do inplace append 135 se = pa.DataFrame(res_ols.bse[None,:], index=['OLS']) 136 se = se.append(pa.DataFrame(res5.bse[None,:], index=['PGLSit5'])) 137 clbse = sw.se_cov(sw.cov_cluster(mod.res_pooled, dgp.groups.astype(int))) 138 se = se.append(pa.DataFrame(clbse[None,:], index=['OLSclu'])) 139 pnwse = sw.se_cov(sw.cov_nw_panel(mod.res_pooled, 4, mod.group.groupidx)) 140 se = se.append(pa.DataFrame(pnwse[None,:], index=['OLSpnw'])) 141 print(se) 142 #list(se.index) 143 from statsmodels.iolib.table import SimpleTable 144 headers = [str(i) for i in se.columns] 145 stubs=list(se.index) 146# print SimpleTable(np.round(np.asarray(se), 4), 147# headers=headers, 148# stubs=stubs) 149 print(SimpleTable(np.asarray(se), headers=headers, stubs=stubs, 150 txt_fmt=dict(data_fmts=['%10.4f']), 151 title='Standard Errors')) 152