PageRenderTime 38ms CodeModel.GetById 27ms app.highlight 8ms RepoModel.GetById 1ms app.codeStats 0ms

/statsmodels/sandbox/examples/ex_random_panel.py

http://github.com/statsmodels/statsmodels
Python | 151 lines | 93 code | 19 blank | 39 comment | 4 complexity | 37555fcfb34788e7eaa084da3ac3346b MD5 | raw file
  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'))