Statsmodels: statistical modeling and econometrics in Python

@josef-pkt I am starting a PR . At the moment there is the gam file that contains the gam penalty class. Smooth_basis contains some functions to get bsplines and polynomial basis. This file will be removed when we will be able to use directly patsy.

There are also 2 files that contains examples or small scripts. We will remove them later.

Let me know what do you think about that.

Todo

• [ ] predict errors (stateful transform, patsy ?), Note fittedvalues are available only a problem in pirls example, predict after fit works, requires spline basis values
• [ ] get_prediction also errors, maybe consequence of predict error currently errors because weights is None
• [ ] check test coverage for offset and exposure

Interface

• [ ] partial_values plot_partial has inconvenient arguments (smoother and mask) instead of column index or term name or similar
• [ ] formula-like interface for predict (create spline basis values internally)
• [ ] adjust inherited methods like plot_partial_residuals (after GSOC?)
• [ ] default param names for splines (should be more informative than "xi" (i is range(k)), but less (?) verbose than patsy's), need them for test_terms

Here's a simple branch with the code added into statsmodels.tsa.statespace. A couple of thoughts

• At least in the dev process, I thought it might be nicer to keep it in its own module, rather than putting it with the kalmanf. I don't know what makes the most sense in the long run.
• I have unit tests that rely on the statespace model, but I'm rewriting them so the KF pull request can be done on its own without other dependencies, especially since the statespace model is likely to change.

Edit: Original line comments.

L72

Question: do we want to keep the single-precision version (and the complex single precision, below)? I don't really see a use case, and it appears from preliminary tests that the results tend to overflow. Maybe I'll post a unit test to demonstrate and we can go from there.

L332

Question: there are a bunch of ways to initialize the KF, depending on the theory underlying the data. This one is only valid for stationary processes. Probably best to move out to the Python wrapper?

L444 Question: I think we'll want to add the ability to specify whether or not to check for convergence and alter the tolerance.

L414

This inversion is using an LU decomposition, but I think in this case I can rely on f to be positive definite since it's the covariance matrix of the forecast error, in which case I could use the much faster Cholesky decomposition approach. This is something I'm looking into, but if you happen to know one way or the other, that would be great too.

Related to this is the idea that you should "never take inverses", and I guess I need to look into replacing this with a linear solver routine if possible, in the updating stage below.

This is a new version of the kernel density estimation. The purpose is to provide an implementation that is faster in the case of grid evaluation, and also works on bounded domains.

There is still some work to do, in particular, I need to add tests for the multi-dimensional and discrete densities.

This model include following models:

1. Generic Zero Inflated model
2. Zero-Inflated Poisson model
3. Zero-Inflated Generalized Poisson model (ZIGP-P)
4. Zero-Inflated Generalized Negative Binomial model (ZINB-P)

Each model include this parts:

• [x] LLF
• [x] Score
• [x] Hessian
• [x] Predict
• [x] Fit
• [x] Docs
• [x] Tests

Status: - reviewing, need to implement better way to generate start_params Last commit for end of GSoC17: Changed way to find start params

This PR include implementation of Generalized Poisson model This model include this parts:

• [x] Log-likelihood function
• [x] Score
• [x] Hessian
• [x] Fit
• [x] Result
• [x] Tests
• [x] Docs

Status - merged #3795

Hi, I have started implementing Kim Filter, outlined a basic functionality, as described in Kim-Nelson book (see diagram on p. 105). I didn't even run the code yet to check for errors. Coding style and class interface bother me more for the moment, as well as the possible ways to test it without implementing models.

xref #3773

I apologize for this taking longer than had hoped, but my busy period was exacerbated by some unforeseen events. But I'm somewhat back now.

This is the first step to solving the above issue.

Still a lot is missing... I've only really gotten loglike to work for all the families except Gamma, Binomial, Gaussian, and Tweedie (which has never had loglike). I also need to re-work the docstrings. I'm thinking its more logical to have a loglike_obs function called in each family and then have a loglike function in the Family class so that loglike will simply be inherited. I think you the doctrings could be elegantly written to handle this too.

I'm thinking I will work on deviance next and then ciricle back to loglike... I feel like R takes some computational shortcuts...

rebased version of #4575 which was rebased version of #2744 The original GSOC PR with most of the development discussion is #2435

rebase conflict in compat.python: This has now unneeded itertools.combinations import but I dropped the py 2.6 compat code.

Hello,

xref #3371

This (should) get var_weights to work for GLM. I added a Tweedie usecase against R.

I can't seem to get the Poisson with var_weights to HC0 test that was added (but disabled because of the lack of functionality) to work. The test is here:

https://github.com/thequackdaddy/statsmodels/blob/var_weight/statsmodels/genmod/tests/test_glm_weights.py#L142

It fails on the following assert

assert_allclose(res1.bse, corr_fact * res2.bse, atol= 1e-6, rtol=2e-6)

Its relatively close... consistently off by a factor of 0.98574... corr_fact brings it much closer... I'm wondering if another adjustment is necessary?

Honestly, I don't really understand much about sandwiches (except PBJ, meatball, and bánh mì).

Rest of the test seem to work pretty well.

I'd be happy to have this merged relatively soon, so thanks for the feedback and review!

PR as discussed in #3274. Tested with a SAS example and two R examples, produced the same results.

To-do:

• [X] Core stats computation
• [x] api
• [x] automatic create dummy variable and hypothesis testing for categorical type independent variables.
• [x] Input validation
• [x] More examples to be tested

references: [1] https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_introreg_sect012.htm) [2] GLM Algorithms ftp://public.dhe.ibm.com/software/analytics/spss/documentation/statistics/20.0/en/client/Manuals/IBM_SPSS_Statistics_Algorithms.pdf [3] http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.278.6976&rep=rep1&type=pdf

Unit test in test_MultivariateOLS.py

compare_r_output_dogs_data()

It reproduces results from the following R code:

library(car)
Drug = c('Morphine', 'Morphine', 'Morphine', 'Morphine', 'Morphine', 'placebo', 'placebo', 'placebo', 'placebo', 'placebo', 'Trimethaphan', 'Trimethaphan', 'Trimethaphan', 'Trimethaphan', 'Trimethaphan')
Depleted = c('N', 'N', 'N', 'N', 'Y', 'Y', 'Y', 'N', 'N', 'N', 'N', 'Y', 'Y', 'Y', 'Y')
subject  = c(1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16)
Histamine0 = c(-3.218876, -3.912023, -2.65926, -1.771957, -2.302585, -2.65926, -2.995732, -3.506558, -3.506558, -2.65926, -2.407946, -2.302585, -2.525729, -2.040221, -2.813411)
Histamine1 = c(-1.609438, -2.813411, 0.336472, -0.562119, -2.407946, -2.65926, -2.65926, -0.478036, 0.04879, -0.18633, 1.141033, -2.407946, -2.407946, -2.302585, -2.995732)
Histamine3 = c(-2.302585, -3.912023, -0.733969, -1.049822, -2.040221, -2.813411, -2.813411, -1.171183, -0.314711, 0.067659, 0.722706, -2.407946, -2.407946, -2.120264, -2.995732)
Histamine5 = c(-2.525729, -3.912023, -1.427116, -1.427116, -1.966113, -2.65926, -2.65926, -1.514128, -0.510826, -0.223144, 0.207014, -2.525729, -2.302585, -2.120264, -2.995732)
data = data.frame(Histamine0, Histamine1, Histamine3, Histamine5)
hismat = as.matrix(data[,1:4])
result = lm(hismat ~ Drug * Depleted)
linearHypothesis(result, c(1, 0, 0, 0, 0, 0)) 
linearHypothesis(result, t(cbind(c(0, 1, 0, 0, 0, 0), c(0, 0, 1, 0, 0, 0))))
linearHypothesis(result, c(0, 0, 0, 1, 0, 0)) 
linearHypothesis(result, t(cbind(c(0, 0, 0, 0, 1, 0), c(0, 0, 0, 0, 0, 1))))
# Or ManRes <- Manova(result, type="III")

test_affine_hypothesis()

It reproduces results from the following R code:

result = lm(hismat ~ Drug*Depleted)
fml = t(cbind(c(0, 1.1, 1.2, 1.3, 1.4, 1.5), c(0, 2.1, 3.2, 3.3, 4.4, 5.5)))
linearHypothesis(result, fml, t(cbind(c(1, 2, 3, 4), c(5, 6, 7, 8))), verbose=TRUE)

Ok, I wanted to make a PR for this, there are still a couple of things that I need to fix but things are pretty close to done and I'm happy with the current approach.

The key part of the current approach is a function distributed_estimation. This function works by taking a generator for endog and exog, endog_generator and exog_generator, as well as a series of functions and key word arguments to be run on each machine and then used to recombine the results. The generator approach allows for a variety of use cases and can handle a lot of the ideas discussed in the initial proposal. For each data set yielded by the generators a model is initialized using model_class and init_kwds and then the function estimation_method is applied to the model along with the key words fit_kwds and estimation_kwds. Finally, the results are recombined from each data set using join_method.

Currently, this defaults to the distributed regularized approach discussed here:

http://arxiv.org/pdf/1503.04337v3.pdf

but the way I've set things up means that the user should be able to apply any number of procedures here.

The current todo list:

• [x] Fix hess_obs
• [x] Fix joblib fit
• [x] Add dask fit
• [x] Add WLS/GLS for debiased regularized fit
• [x] Add likelihood result
• [x] Add data tests

Let me know if there are any comments/questions/criticisms, as I've said, it certainly isn't complete but I wanted to get this out there so I can start integrating any changes as I finish it up.

Is your feature request related to a problem? Please describe

Whenever I need to perform a non-parametric test using Chebyshev's Inequality, I have to write the equation myself. It would be good to have it as a part of statsmodel

Describe the solution you'd like

A function which takes in a distribution(array), a scalar value to be used for the test, the type of test(one tailed/two tailed), confidence and returns whether the scalar value is within the limits

Describe alternatives you have considered

Currently I have to write the code myself so I don't have any other alternative, would be better to have it as a part of statsmodel

Additional context

If anyone thinks this will be helpful, do let me know, I can raise a pull request.

I am currently looking for how to implement robust poisson regression in Python. I found one article explaining the way to implement the robust poisson regression in R language using the similar method implemented in statsmodels.

• POISSON REGRESSION | R DATA ANALYSIS EXAMPLES

The following code creates the output in R,

library(sandwich)
p <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
p <- within(p, {
  prog <- factor(prog, levels=1:3, labels=c("General", "Academic", 
                                                     "Vocational"))
  id <- factor(id)
})
summary(p)
summary(m1 <- glm(num_awards ~ prog + math, family="poisson", data=p))
cov.m1 <- vcovHC(m1, type="HC0")
std.err <- sqrt(diag(cov.m1))
r.est <- cbind(Estimate= coef(m1), "Robust SE" = std.err,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err,
UL = coef(m1) + 1.96 * std.err)

r.est

and here is the corresponding code in statsmodels.

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
data = pd.read_csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
res = smf.poisson("num_awards ~ C(prog) + math", data=data).fit(cov_type="HC0")
res.summary()

It seems the results from both codes are matched quite well. However, exploring the past threads about heteroscedasticity for poisson regression, for example

• Standard Error differences Binomial Regression vs Poisson Regression
• ENH: handling cluster_weights, freq_weights in models, cov_type, moments, descriptive statistics

I am afraid even cov_type HC0 is correctly implemented or not. If you have time, could you explain the current implementation situation around heteroscedasticity of the poisson regression? What is the potential problem not covered by the current package, and what is currently firmly implemented?

see https://github.com/statsmodels/statsmodels/issues/8538#issuecomment-1366169937

errors in pre-release testing of scipy extradoc has been removed from keyword arguments

this comments out extradoc keyword in sandbox distributions no alternative How to modify docs to add extra info is open.

Hi all!

This PR closes #8461.

It adds support for CRV3 cluster robust inference for OLS and WLS via a cluster jackknife.

In short, the PR adds a new cov_type, cluster-jackknife, which can be run in the following way:

crv3 = model.fit(cov_type='cluster-jackknife', cov_kwds={'groups': id})

I have updated the docstring in linear_model.py.

I have not yet added tests, as I'd first like to get feedback on where and how to test the new cov_type (and I also have to get up to speed with pytest) . There are multiple strategies how to test the new feature:

• test against implementations in other languages (e.g. Stata's summclust, R's summclust, or R's sandwich::vcovJK).
• internal tests: under appropriate small sample corrections and when clusters are singletons, the cluster jackknife is equivalent to the HC3 estimator.

Below is a small example of the ladder test, using the new cluster-jackknife covtype and HC3 inference:

from statsmodels.regression.linear_model import OLS
from statsmodels.datasets import grunfeld
from statsmodels.tools.tools import add_constant
import pandas as pd
import numpy as np

dtapa = grunfeld.data.load_pandas()
dtapa_endog = dtapa.endog[:200]
dtapa_exog = dtapa.exog[:200]
exog = add_constant(dtapa_exog[["value", "capital"]], prepend=False)
N = exog.shape[0]
id = pd.Series(range(0, N))

res_hc3 = OLS(dtapa_endog, exog).fit(cov_type = "HC3")
res_crv3 = OLS(dtapa_endog, exog).fit(cov_type = "cluster-jackknife", cov_kwds={'groups': id})

res = pd.DataFrame({'tstat hc3:': res_hc3.tvalues * np.sqrt((N-1) / (N)), 'tstat crv3:' : res_crv3.tvalues})
res
#          tstat hc3:  tstat crv3:
# value     16.093573    16.093573
# capital    3.932688     3.932688
# const     -3.040458    -3.040458

Let me know which types of tests you'd like to see =)

Is your feature request related to a problem? Please describe

I want to use "X-13ARIMA-SEATS Seasonal Adjustment Program" via Python. I know statsmodels.tsa.x13.x13_arima_analysis is the Python wrapper.

Using statsmodels, is it possible to provide the inputs (specification) files into the Seasonal adjustment procedure? Eg, I want to add/drop certain outliers or regressors. How do I go about doing it in statsmodels?

Thank you

Describe the solution you'd like

A clear and concise description of what you want to happen.

Describe alternatives you have considered

A clear and concise description of any alternative solutions or features you have considered.

Additional context

Add any other context about the feature request here.

(semi-random idea while browsing SAS docs for logistic and firth penalization)

related to Firth penalization #3561

reference Heinze, G., and Schemper, M. (2002). “A Solution to the Problem of Separation in Logistic Regression.” Statistics in Medicine 21:2409–2419.

penalized loglikelihood is log L*(beta) = log L*(beta) + 1/2 log | I(beta)| where I(beta) is the information matrix evaluated at beta

That means we need the penalization class to have access to the model so we can reuse the hessian (OIM or EIM) This might currently be tricky or impossible because we define the penalization instance before creating a penalized model instance.

One possibility may be to assign a _model attribute to the penalty instance after creating the model. This makes it a bit circulate The penalized model instance holds the penalty instance which holds the penalized model instance as attributes. I guess the circularity cannot be avoided unless we use extra args in each method of the penalty class.

Note we need the hessian/EIM of the original model and not the penalized hessian of the penalized model.