Sunday, February 12, 2012

Using rst2blogger and Breush-Godfrey test

This is mainly a test to see if rst2blogger works. The description sounded like just what I needed to make going from rst file to blogger less work. Some setup inconveniences: Why do python programs need to install distribute and destroy my setuptools? Lxml didn't have windows binaries, and easy_install fails because the source doesn't build. In contrast to sphinx, the plain docutils do not highlight python code, at least not right away. If there is an exception in rendering the html, then the message is completely uninformative: except Exception

(Reminds me of the sound I'm hearing when the kids are just guessing in a numbers game: "Try again", "Try again", ...)

try a rst list: some categories of tests that I have been working on recently.

  • diagnostic tests
    • autocorrelation
    • heteroscedasticity
    • multicollinearity
    • functional form
    • normality
    • outliers
    • influence
  • tests for coefficients/parameters
    • t test for linear restrictions
    • F test for linear restrictions
    • Wald test for linear restrictions missing
    • Likelihood Ratio test to compare nested models
    • Wald/F test test to compare nested models
    • comparing non-nested models: J, Cox

Writing Breush-Godfrey autocorrelation test in 5 minutes

The background story: I have written several diagnostics tests a while ago, maybe one and a half years ago. Many of the diagnostic test have a similar structure and it's easy to just copy the pattern. For one test, acorr_lm, I took the code and basic idea of the augmented dickey fuller test for unit roots, and wanted to make a Lagrange Multiplier test for autocorrelation out of it. Some time later, I added some comments that this could be used as Engle's ARCH test but that there is a difference to Breush-Godfrey's serial correlation test.

While I was writing tests comparing my diagnostic tests with R, I saw that Breush-Godfrey was still missing. Since this time I had unit tests, it was very quick work, copy and paste, and checking Wikipedia .

The test passed after a few changes, and I barely had to think.

Build the matrix of lagged residuals, column_stack them with the original regressors, get the matrix for the linear restrictions that the joint effect of the lagged residuals is zero. Most of the work is done by existing functions.

I haven't changed the docstring yet, and there is still cleanup work to do.

def acorr_breush_godfrey(results, nlags=None, store=False):
    '''Lagrange Multiplier tests for autocorrelation

    not checked yet, copied from unitrood_adf with adjustments
    check array shapes because of the addition of the constant.
    written/copied without reference
    This is not Breush-Godfrey. BG adds lags of residual to exog in the
    design matrix for the auxiliary regression with residuals as endog,
    see Greene 12.7.1.

    Notes
    -----
    If x is calculated as y^2 for a time series y, then this test corresponds
    to the Engel test for autoregressive conditional heteroscedasticity (ARCH).
    TODO: get details and verify

    '''

    x = np.concatenate((np.zeros(nlags), results.resid))
    exog_old = results.model.exog
    x = np.asarray(x)
    nobs = x.shape[0]
    if nlags is None:
        #for adf from Greene referencing Schwert 1989
        nlags = 12. * np.power(nobs/100., 1/4.)#nobs//4  #TODO: check default, or do AIC/BIC

    xdall = lagmat(x[:,None], nlags, trim='both')
    nobs = xdall.shape[0]
    xdall = np.c_[np.ones((nobs,1)), xdall]
    xshort = x[-nobs:]
    exog = np.column_stack((exog_old, xdall))
    k_vars = exog.shape[1]

    if store: resstore = ResultsStore()

    resols = sm.OLS(xshort, exog).fit()
    ft = resols.f_test(np.eye(nlags, k_vars, k_vars - nlags))
    fval = ft.fvalue
    fpval = ft.pvalue
    fval = np.squeeze(fval)[()]   #TODO: fix this in ContrastResults
    fpval = np.squeeze(fpval)[()]
    lm = nobs * resols.rsquared
    lmpval = stats.chi2.sf(lm, nlags)
    # Note: degrees of freedom for LM test is nvars minus constant = usedlags
    #return fval, fpval, lm, lmpval

    if store:
        resstore.resols = resols
        resstore.usedlag = nlags
        return fval, fpval, lm, lmpval, resstore
    else:
        return fval, fpval, lm, lmpval

Just a brief explanation to the main tools:

lagmat just creates a matrix of lagged values without having to worry about how to shift and cut the arrays.

>>> from scikits.statsmodels.tsa.tsatools import lagmat
>>> lagmat(np.arange(8), 3, trim='both')
array([[ 2.,  1.,  0.],
       [ 3.,  2.,  1.],
       [ 4.,  3.,  2.],
       [ 5.,  4.,  3.],
       [ 6.,  5.,  4.]])

numpy.eye is very flexible and allows shifting of the diagonal

>>> np.eye(3, 5, 5-3)
array([[ 0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  1.]])

The rest are mainly calls to OLS and f_test. One part I needed to change was that R takes the initial errors to be zero instead of truncating the regression, so I also had to concatenate some zeros in front of the regression residuals array.

Friday, February 10, 2012

Starting a list of diagnostic and specification tests

I started already several times an inventory of statistical tests in python, scipy.stats and statsmodels, compared to R.

Here is another try.

It is mainly the comparison with the R package lmtest. I just spend most of two days, writing tests against R, and before that some days of writing tests against Gretl, and before that outlier measures against SAS, as described in the previous posts.

I don't know yet how easy it will be to maintain a list like this, the current version is mainly based on parsing the lmtest index page. lmtest is not very complete, and there are tests covered in other packages and additional tests covered in Gretl.

For now I just keep it in a boring python module, with a string that's easy to manipulate and convert.

#cols: category, name, statsmodels, r_lmtest, gretl
s4 = '''\
acorr | Breusch-Godfrey Test | acorr_breush_godfrey | bgtest
acorr | Durbin-Watson Test | location_no_pvalues | dwtest
het | Breusch-Pagan Test | het_breush_pagan | bptest
het | Goldfeld-Quandt Test | het_goldfeld_quandt | gqtest
het | Harrison-McCabe test | missing | hmctest
het | White test | het_white | - |
causality | Test for Granger Causality | grangercausalitytest | grangertest
linear | Harvey-Collier Test | missing | harvtest
linear | PE Test for Linear vs. Log-Linear Specifications | missing | petest
linear | Rainbow Test | missing | raintest
func form | RESET Test | with outliers | resettest
compare nonnested | Cox Test | compare_cox | coxtest
compare nonnested | J Test | compare_j | jtest
compare nonnested | Encompassing Test | missing | encomptest
compare nested | Likelihood Ratio Test nested | compare_lr | lrtest
compare nested | Wald Test nested | compare_ftest | waldtest
coef | Testing Estimated Coefficients | t_test | coeftest
coef | Testing Estimated Coefficients | missing | coeftest.breakpointsfull
'''

add another separator
print '\n'.join(line + '|' for line in s4.split('\n'))
convert to list of list 
def str2list(ss, sep='|', keep_empty=4):
Unfortunately copying into blogger doesn't preserve intend, so skip this. And convert separator to tabs, so that google spreadsheet separates the cells:
print '\n'.join(line + '|' for line in s4.split('\n'))

Tomorrow, I will start looking for the diagnostic tests that are not yet on the list. R stats has some, for example (fm is my test case linear model result)

> names(ls.diag(fm))
 [1] "std.dev"      "hat"          "std.res"      "stud.res"     "cooks"        "dfits"      
 [7] "correlation"  "std.err"      "cov.scaled"   "cov.unscaled"
I figured out that json works pretty well transferring data from some R animals to python.

In other news:
The basic R doesn't save automatically the sessionlog or history. I was playing with Rcommander last night to see what default diagnostics they are proposing, and it crashed R after a while. Unfortunately, I hadn't saved my sessionlog and script file, so a day or two of work died with it. I didn't think about safeguarding against crashes anymore, Windows never crashes, not even the kids manage to turn it off anymore, spyder and firefox always recover after a crash with an existing history or session log. R never crashed before.


Tuesday, January 31, 2012

The nice thing about seeing zeros

Now, why am I happy to see a result like this?
>>> np.round(table_sm - table_sas, 3)
array([[ 0.,  0.,  0., -0.,  0., -0., -0.,  0., -0.,  0., -0.,  0.],
       [ 0.,  0., -0., -0.,  0., -0., -0.,  0.,  0.,  0., -0.,  0.],
       [ 0.,  0., -0., -0., -0.,  0., -0., -0.,  0., -0.,  0., -0.],
       [ 0.,  0., -0.,  0., -0.,  0., -0.,  0.,  0.,  0., -0.,  0.],
       [ 0.,  0., -0., -0.,  0., -0., -0.,  0.,  0., -0.,  0.,  0.],
       [ 0.,  0., -0.,  0.,  0.,  0., -0., -0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -0.,  0., -0., -0., -0., -0.,  0.],
       [ 0.,  0., -0., -0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -0., -0.,  0., -0., -0.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -0., -0., -0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -0.,  0., -0.,  0., -0.,  0.],
       [ 0.,  0.,  0., -0., -0.,  0., -0., -0., -0.,  0., -0., -0.],
       [ 0.,  0., -0.,  0.,  0.,  0.,  0., -0.,  0., -0., -0., -0.],
       [ 0.,  0.,  0.,  0., -0., -0.,  0., -0., -0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -0., -0.,  0.,  0., -0., -0., -0.,  0.,  0.],
       [ 0.,  0., -0.,  0.,  0., -0.,  0., -0.,  0.,  0., -0.,  0.],
       [ 0.,  0.,  0., -0.,  0., -0.,  0., -0., -0., -0., -0., -0.],
       [ 0.,  0.,  0., -0.,  0.,  0.,  0.,  0., -0.,  0.,  0., -0.],
       [ 0.,  0., -0., -0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0.]])
It means I get the same results, up to print precision, as SAS in this table http://www.sfu.ca/sasdoc/sashtml/stat/chap55/sect4.htm#regg2f
table_sm = np.column_stack([
                              np.arange(res.nobs) + 1,
                              res.model.endog,
                              res.fittedvalues,
                              predict_mean_se,
                              predict_mean_ci[:,0],
                              predict_mean_ci[:,1],
                              predict_ci[:,0],
                              predict_ci[:,1],
                              res.resid,
                              resid_se,
                              infl.resid_studentized_internal,
                              infl.cooks_distance()[0]
                              ])
from numpy.testing import assert_almost_equal
assert_almost_equal(table_sm, table_sas, decimal=4)
>>> print st
=====================================================================================================================================
       Obs    Dep Var  Predicted    Std Error    Mean ci    Mean ci Predict ci Predict ci   Residual  Std Error    Student     Cook's
           Population      Value Mean Predict    95% low    95% upp    95% low    95% upp              Residual   Residual          D
-------------------------------------------------------------------------------------------------------------------------------------
         1      3.929      5.038        1.729      1.373      8.704     -1.903     11.980     -1.109      2.178     -0.509      0.055
         2      5.308      5.039        1.391      2.090      7.987     -1.553     11.631      0.269      2.408      0.112      0.001
         3      7.239      6.309        1.130      3.912      8.705     -0.055     12.672      0.930      2.541      0.366      0.009
         4      9.638      8.847        0.957      6.818     10.876      2.612     15.082      0.791      2.611      0.303      0.004
         5     12.866     12.655        0.872     10.806     14.504      6.476     18.834      0.211      2.641      0.080      0.000
         6     17.069     17.732        0.858     15.913     19.550     11.562     23.901     -0.663      2.645     -0.251      0.002
         7     23.191     24.078        0.883     22.205     25.951     17.892     30.264     -0.887      2.637     -0.336      0.004
         8     31.443     31.693        0.920     29.742     33.644     25.483     37.903     -0.250      2.624     -0.095      0.000
         9     39.818     40.577        0.949     38.566     42.589     34.348     46.807     -0.759      2.614     -0.290      0.004
        10     50.155     50.731        0.959     48.697     52.764     44.494     56.967     -0.576      2.610     -0.221      0.002
        11     62.947     62.153        0.949     60.142     64.164     55.924     68.382      0.794      2.614      0.304      0.004
        12     75.994     74.845        0.920     72.894     76.796     68.635     81.055      1.149      2.624      0.438      0.008
        13     91.972     88.806        0.883     86.933     90.679     82.620     94.992      3.166      2.637      1.201      0.054
        14    105.710    104.035        0.858    102.217    105.854     97.866    110.205      1.675      2.645      0.633      0.014
        15    122.775    120.534        0.872    118.686    122.383    114.356    126.713      2.241      2.641      0.849      0.026
        16    131.669    138.303        0.957    136.274    140.332    132.068    144.537     -6.633      2.611     -2.541      0.289
        17    151.325    157.340        1.130    154.943    159.736    150.976    163.704     -6.015      2.541     -2.367      0.370
        18    179.323    177.646        1.391    174.697    180.595    171.054    184.238      1.677      2.408      0.696      0.054
        19    203.211    199.221        1.729    195.556    202.887    192.280    206.163      3.990      2.178      1.831      0.704
===================================================================================================================================== 
 
This is again just produced with SimpleTable, and I didn't spend any time on formatting yet.

(aside: I'm still trying to figure out how to create an article without a large administrative overhead. This time I used sphinx but it doesn't quite work yet.)