tag:blogger.com,1999:blog-1282744976876626082024-03-20T02:33:22.374-04:00Scipy Stats ProjectThe statsmodels project started as part of the Google Summer of Code 2009. Now that the GSoC is officially over, this blog will be a place to learn about updates to the project. Any comments and questions are welcome. Anyone who wishes to help with development is very welcome! Discussion of the project will take place on the scipy-dev mailing list.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.comBlogger23125tag:blogger.com,1999:blog-128274497687662608.post-53637112493759051812011-08-22T23:51:00.000-04:002011-08-22T23:51:16.284-04:00SVAR EstimationFor the next phase of my GSoC project, I integrated SVAR estimation with restrictions on the within period effects and shock identification. This follows the method outlined in section 11.6 of Hamilton (1994) [1]. In order to show how the system works, I will work through an example.<br />
<br />
Structural VARs of the type I will be positing here are best suited to situations where impulse response function shocks need to be identified, but there is little to no theoretical justification for the within period dynamics required for orthogonalization. While this allows greater flexibility to the researcher in specifying the model, it also leaves greater room for error, especially when specifying what parameters can be assumed zero and what should be estimated, and whether or not this specification meets the criteria for shock identification.<br />
<br />
Given time series data, SVAR class, initiation requires (at a minimum):<br />
<br />
1) svar_type<br />
2) A, B matrices marking with an 'E' or 'e' where parameters are unknown<br />
<br />
This must be A,B, or AB. An 'A' system assumed that the matrix premulitplying the error matrix is the identity matrix and vice-versa for the 'B' system. In an 'AB' system, elements in both the A and B matrix need to be estimated. The system can be summarized as:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=Ay_t = c @plus; A_1y_{t-1} @plus; A_2y_{t-2} @plus; \cdots @plus; A_py_{t-p} @plus; B\varepsilon_t" target="_blank"><img src="http://latex.codecogs.com/gif.latex?Ay_t = c + A_1y_{t-1} + A_2y_{t-2} + \cdots + A_py_{t-p} + B\varepsilon_t" title="Ay_t = c + A_1y_{t-1} + A_2y_{t-2} + \cdots + A_py_{t-p} + B\varepsilon_t" /></a><br />
<br />
We estimate A, B in a second stage after A_1-A_n, /Sigma_u are estimated via OLS, by maximizing the log-likelihood:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=\dpi{80} \log{L} = - \frac{Tn}{2}\log(2π) @plus; \frac{T}{2}\log|A|^2 - \frac{T}{2}\ln|B|^2 - \frac{T}{2}trace(A'B'^{-1}B^{-1}A\Sigma_u)" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\dpi{80} \log{L} = - \frac{Tn}{2}\log(2π) + \frac{T}{2}\log|A|^2 - \frac{T}{2}\ln|B|^2 - \frac{T}{2}trace(A'B'^{-1}B^{-1}A\Sigma_u)" title="\dpi{80} \log{L} = - \frac{Tn}{2}\log(2π) + \frac{T}{2}\log|A|^2 - \frac{T}{2}\ln|B|^2 - \frac{T}{2}trace(A'B'^{-1}B^{-1}A\Sigma_u)" /></a><br />
<br />
We could set up an example system as:<br />
<br />
<pre><div class="mycode">In [1]: import numpy as np
In [2]: A = np.array([[1, 'E', 0], [1, 'E', 'E'], [0, 0, 1]])
In [3]: A
Out[3]:
array([['1', 'E', '0'],
['1', 'E', 'E'],
['0', '0', '1']],
dtype='|S1')
In [4]: B = np.array([['E', 0, 0], [0, 'E', 0], [0, 0, 'E']])
In [5]: B
Out[5]:
array([['E', '0', '0'],
['0', 'E', '0'],
['0', '0', 'E']],
dtype='|S1')
</div></pre><br />
In order to aid numerical maximum likelihood estimation, the SVAR class fit can also be passed guess matrices for both A and B parameters.<br />
<br />
Building upon the previously used three variable macro system, a simple example of an AB system estimation is as follows:<br />
<br />
<pre><div class="mycode">import numpy as np
import scikits.statsmodels.api as sm
from scikits.statsmodels.tsa.api import VAR, SVAR
import matplotlib.pyplot as plt
mdata = sm.datasets.macrodata.load().data
mdata = mdata[['realgdp','realcons','realinv']]
names = mdata.dtype.names
data = mdata.view((float,3))
data = np.diff(np.log(data), axis=0)
#define structural inputs
A = np.asarray([[1, 0, 0],['E', 1, 0],['E', 'E', 1]])
B = np.asarray([['E', 0, 0], [0, 'E', 0], [0, 0, 'E']])
A_guess = np.asarray([0.5, 0.25, -0.38])
B_guess = np.asarray([0.5, 0.1, 0.05])
mymodel = SVAR(data, svar_type='AB', A=A, B=B, names=names)
res = mymodel.fit(maxlags=3, maxiter=10000, maxfun=10000, solver='bfgs')
res.irf(periods=30).plot(impulse='realgdp', plot_stderr=False)
</div></pre><br />
From here, we can access the estimates of both A and B:<br />
<br />
<pre><div class="mycode">In [2]: res.A
Out[2]:
array([[ 1. , 0. , 0. ],
[-0.50680204, 1. , 0. ],
[-5.53605672, 3.04117688, 1. ]])
In [3]: res.B
Out[3]:
array([[ 0.00757566, 0. , 0. ],
[ 0. , 0.00512051, 0. ],
[ 0. , 0. , 0.02070894]])
</div></pre><br />
This also produces a SVAR irf, resulting from a impulse shock to realgdp:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgcavkSB9xfUyL50GKKqnAEmY3xUCPM4RB9GVoZTBknRho2qIRVQOecwVdFpP8EFd0uForM8XlWHOLo-b4fW3OtsdR1IJBxp2U9NBZovBRUm4aiRx6O4wobe9k5D96up_aBFNSr5CdBM7U/s1600/svar1.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgcavkSB9xfUyL50GKKqnAEmY3xUCPM4RB9GVoZTBknRho2qIRVQOecwVdFpP8EFd0uForM8XlWHOLo-b4fW3OtsdR1IJBxp2U9NBZovBRUm4aiRx6O4wobe9k5D96up_aBFNSr5CdBM7U/s320/svar1.png" /></a></div><br />
We can check our work by comparing our work to comparable packages (see below) and also performing some simple calculations. If we calculated A and B correctly, then:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=\Sigma_u = A^{-1}BB(A^{-1})'" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\Sigma_u = A^{-1}BB(A^{-1})'" title="\Sigma_u = A^{-1}BB(A^{-1})'" /></a><br />
<br />
We find:<br />
<br />
<pre><div class="mycode">In [10]: P = np.dot(npl.inv(res.A), res.B)
In [11]: np.dot(P,P.T)
Out[11]:
array([[ 5.73906806e-05, 2.90857141e-05, 2.29263262e-04],
[ 2.90857141e-05, 4.09603703e-05, 3.64524316e-05],
[ 2.29263262e-04, 3.64524316e-05, 1.58721639e-03]])
In [12]: res.sigma_u
Out[12]:
array([[ 5.73907395e-05, 2.90857557e-05, 2.29263451e-04],
[ 2.90857557e-05, 4.09604397e-05, 3.64524456e-05],
[ 2.29263451e-04, 3.64524456e-05, 1.58721761e-03]])
</div></pre><br />
Through much trial, it seems like the 'bfgs' estimator is the best for solving the maximum likelihood problem. I would like to investigate in the future why this may be. Estimating the likelihood was aided greatly by bringing in the likelihood methods already present in the base component of the statsmodels package.<br />
<br />
This equivalent system can be estimated and plotted in the R [2] script:<br />
<br />
<pre><div class="mycode">library("vars")
#data <- read.csv("/home/skipper/statsmodels/statsmodels-skipper/scikits/statsmodels/datasets/macrodata/macrodata.csv")
data <- read.csv("/home/bart/statsmodels/scikits/statsmodels/datasets/macrodata/macrodata.csv")
names <- colnames(data)
data <- log(data[c('realgdp','realcons','realinv')])
data <- sapply(data, diff)
data = ts(data, start=c(1959,2), frequency=4)
amat <- matrix(0,3,3)
amat[1,1] <- 1
amat[2,1] <- NA
amat[3,1] <- NA
amat[2,2] <- 1
amat[3,2] <- NA
amat[3,3] <- 1
bmat <- diag(3)
diag(bmat) <- NA
svar <- SVAR(var, estmethod = 'scoring', Bmat=bmat, Amat=amat)
plot(irf(svar, n.ahead=30, impulse = 'realgdp'))
#myirf <- plot(irf(myvar, impulse = "realgdp", response = c("realgdp", "realcons", "realinv"), boot=TRUE, n.ahead=30, ci=0.95))
#plot.irf()
</div>
</pre><br />
This is the end of my Google Summer of Code project. In the future, I hope to continue work on SVAR, bring in long-run restrictions ala Blanchard and Quah, and further test solvers and their performance. I have benefited a lot from this project and would like to sincerely thank my mentors: Skipper Seabold, Josef Perktold, and Alan Isaac. They have been a great support and have answered my questions swiftly and completely. All blame for failure to complete the goals I set for myself at the beginning of the summer rests on me alone. Not only have I learned a lot about time series econometrics, but also quite a bit about how community software development works and especially realistic timelines. This has been an invaluable experience and I plan to further improve my contributions to the project in the coming year.<br />
<br />
Bart Baker<br />
<br />
[1] Hamilton, James. 1994. Time Series Analysis. Princeton University Press: Princeton. <br />
[2] R Development Core Team. 2011. R: A Language and Environment for Statistical Computing. R FOundation for Statistical Copmuting. <http://www.R-project.org/><br />
Barthttp://www.blogger.com/profile/12636757026304914432noreply@blogger.com2tag:blogger.com,1999:blog-128274497687662608.post-85997665856174632942011-07-19T00:29:00.001-04:002011-07-19T12:32:16.426-04:00Sims and Zha IRF Error BandsAfter completing the Monte Carlo error bands, I moved onto integrating Sims and Zha error bands into the statsmodels package.These are based on pages 1127 to 1129 of Chris Sims and Tao Zha's 1999 Econometrica (Vol. 67 No. 5) article, "Error Bands for Impulse Responses."<br />
<br />
This method took a long time just to get my head around and a lot of trial and error. While Sims and Zha focus on Bayesian sampling methods from the joint distribution of the coefficients and covariance matrix to generate the draws of the MA(n) representations, the method used here to generate these draws are simpler Monte Carlo simulations. In order for these error bands to truly follow the prescription of Sims and Zha (SZ), the Bayesian sampling methods would need to be employed.<br />
<br />
Here's a quick overview of the theory.<br />
<br />
Given a covariance matrix Sigma, we can perform eigenvalue decomposition as such:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=W{\Lambda}W'=\Sigma" target="_blank"><img src="http://latex.codecogs.com/gif.latex?W{\Lambda}W'=\Sigma" title="W{\Lambda}W'=\Sigma" /></a><br />
<br />
<br />
where Lambda is diagonal and each diagonal element of Lambda corresponds to an eigenvalue of Sigma. Column 'k' of W is the eigenvector corresponding to the 'k'th diagonal element of Lambda or the 'k'th eigenvector.<br />
<br />
We will also define our moving average representation or impulse response function as:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=c_{ij}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?c_{ij}" title="c_{ij}" /></a><br />
<br />
where c is the time series vector ('t' to 't+h') response of variable 'i' to a shock in period 't' to variable 'j'.<br />
<br />
SZ make the case when the time series model (VAR) is fit to data that is not smooth (differenced, etc.), most of the variation will be contained in the principal components of W. With this information, SZ propose three methods for characterizing the uncertainty around the data.<br />
<br />
The following three arrays of graphs can be produced by the following code:<br />
<br />
<pre><div class="mycode">import numpy as np
import scikits.statsmodels.api as sm
from scikits.statsmodels.tsa.api import VAR
import matplotlib.pyplot as plt
mdata = sm.datasets.macrodata.load().data
mdata = mdata[['realgdp','realcons','realinv']]
names = mdata.dtype.names
data = mdata.view((float,3))
data = np.diff(np.log(data), axis=0)
mymodel = VAR(data,names=names)
res = mymodel.fit(maxlags=3,ic=None)
res.irf(periods=20).plot(impulse='realgdp', stderr_type='sz1', repl=1000, seed=30)
res.irf(periods=20).plot(impulse='realgdp', stderr_type='sz2',repl=1000, seed=30)
res.irf(periods=30).plot(impulse='realgdp', stderr_type='sz3', repl=1000, seed=30)
plt.show()
</div></pre><br />
At this point, I have only implemented this for the non-orthogalized impulse responses, but Sims and Zha explicitly address this case in their paper and it is analogous to the methods defined below.<br />
<br />
1) Symmetric, assumes Gaussian uncertainty. These error bands add and subtract the estimated impulse response functions with error completely defined by the principal component<sup><a href="#fn1" id="ref1">1</a></sup>:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=c_{ij}\pm W_{\cdot{k}}(t)\sqrt{\lambda_k}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?c_{ij}\pm W_{\cdot{k}}(t)\sqrt{\lambda_k}" title="c_{ij}\pm W_{\cdot{k}}(t)\sqrt{\lambda_k}" /></a><br />
<br />
where W_{.k} is the column of W corresponding to the 'k'th eigenvalue of Sigma. The above equation would be the 68% probability bands. 95% would simply be mulitplied by a scalar of 1.96 on both sides.<br />
<br />
In our three variable model, this method produces the following error bands a response to a gdp shock:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhDqBHBYRrRbUtT8qsWDCtbCFYMVMCxaKaxCHjZO_glZ8EEBSsy0ad_lbEqHbWdL2dkIKuNtBqykZmArbzWUWEEZUQ44Vi2hnjjt62Eo5fiRHSQldptzZs0gSxgJ6EQ7V23_escK6M2ZEY/s1600/SZ1.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhDqBHBYRrRbUtT8qsWDCtbCFYMVMCxaKaxCHjZO_glZ8EEBSsy0ad_lbEqHbWdL2dkIKuNtBqykZmArbzWUWEEZUQ44Vi2hnjjt62Eo5fiRHSQldptzZs0gSxgJ6EQ7V23_escK6M2ZEY/s320/SZ1.png" /></a></div><br />
Looks nice enough. These are defined as implemented to be symmetric, so I wasn't too worried about these.<br />
<br />
2) Non-symmetric error bands generated by Monte Carlo draws where only covariance between time but not across variables exists.<br />
<br />
In this case, instead of assuming Gaussian uncertainty, we retain the draws used to estimate the initial covariance matrix, Sigma, and for each of the these draws we calculate the vector gamma_k:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=\gamma_k=W_{k\cdot }\times c_{ij}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\gamma_k=W_{k\cdot }\times c_{ij}" title="\gamma_k=W_{k\cdot }\times c_{ij}" /></a><br />
<br />
where W_{k.} is the 'k'th row of W, and k refers to the largest eigenvalue of Sigma. Using these quantiles of the individual elements of gamma_k across the MC draws, we can generate 68% probability bands as follows:<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=\hat{c}_{ij}@plus;\gamma_{k,.16}, \hat{c}_{ij}@plus;\gamma_{k,.84}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\hat{c}_{ij}+\gamma_{k,.16}, \hat{c}_{ij}+\gamma_{k,.84}" title="\hat{c}_{ij}+\gamma_{k,.16}, \hat{c}_{ij}+\gamma_{k,.84}" /></a><br />
<br />
where the subscripts on gamma_k refer to the 16th and 84th percentile of the gamma draws. In our three variable case, this produces the following graphical representation:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj9vO060OqDEzY2Q2AvYIZJjXsXFpkyqnr9o6nsFXYnaTc-KO7gUbr5CelDdNbXxSAI7B23Ej4Jiodi5P7-KtU05D_9o8jFVhdobhBxknqZeN2s8bqWLrlrOcrvWSB68HXTIDc-pVhS6Zc/s1600/SZ2.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj9vO060OqDEzY2Q2AvYIZJjXsXFpkyqnr9o6nsFXYnaTc-KO7gUbr5CelDdNbXxSAI7B23Ej4Jiodi5P7-KtU05D_9o8jFVhdobhBxknqZeN2s8bqWLrlrOcrvWSB68HXTIDc-pVhS6Zc/s320/SZ2.png" /></a></div><br />
We can see that in this case the uncertainty clearly drops precipitously once we hit a certain t in the time series representation.<br />
<br />
Just for completeness, if we look at the response of the same variables to a shock to real consupmtion, we notice that this method does not even guarantee that the probability bands contain the estimated impulse response function:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUJ5pMWJNCHljJkWhsaIaGKs_x1DaA5IT0gYk4yuUPoc65XOMQk92lAX0YHC4NYN70YeCOlWZYbeRUW7w_3xpqFBYsUOtmcH1OVMeuZEDnHH2mt1fEzLfTEI-2ExzrFvjwgH3nJgoIvSY/s1600/SZ2a.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUJ5pMWJNCHljJkWhsaIaGKs_x1DaA5IT0gYk4yuUPoc65XOMQk92lAX0YHC4NYN70YeCOlWZYbeRUW7w_3xpqFBYsUOtmcH1OVMeuZEDnHH2mt1fEzLfTEI-2ExzrFvjwgH3nJgoIvSY/s320/SZ2a.png" /></a></div><br />
While this was extremely worrisome at first, Sims and Zha do give some examples where the probability bands do not contain the estimated impulse response function for certain 't'. These error bands are supposed to give the researcher an idea as to the symmetry (or lackthereof) of the posterior distribution of the impulse response functions.<br />
<br />
3) Non-symmetric error bands generated by Monte Carlo draws where covariance between time and between variables is identified.<br />
<br />
Here, instead of treating each vector of responses individually, we consider each set of impulse response functions as a single system. Sims and Zha note that while in most cases the majority of the covariance will be between intertemporal observations of a single variable, considering inter-variable time series covariance may be valuable in certain situations.<br />
<br />
In order to investigate how different c_ij related to each other, we stack each set of impulse response functions that respond to a single shock j and then compute the eigenvector decomposition of this stacked vector for each MC draw. This allows us to compute eigenvectors that contain the variation both across time periods and across variables. We calculate our gamma_k in an analagous manner to the second method and add the appropriate gamma_k quantiles to the estimated response.<br />
<br />
It happens that in our system, the cross-variable covariance does not reveal much addition information about a shock to GDP:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjbnf4qm2XyLVYguwHtMuDQeiz2T33ql4qki756xwBIqbIAt3F4sYc2XtY7sLOaZfN_xIjbuWxzh5Ai-qxciqijZKWS_pIRTK9spvW7rHQ6ZGwGoQhncTnn6YvwQdaoSsIzQtudVf0hPlQ/s1600/SZ3.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjbnf4qm2XyLVYguwHtMuDQeiz2T33ql4qki756xwBIqbIAt3F4sYc2XtY7sLOaZfN_xIjbuWxzh5Ai-qxciqijZKWS_pIRTK9spvW7rHQ6ZGwGoQhncTnn6YvwQdaoSsIzQtudVf0hPlQ/s320/SZ3.png" /></a></div><br />
It seems as though the variation in real GDP dominates the variation in the other variables.<br />
<br />
This can also be seen by examining the response of real GDP to all three shocks:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhMVJhKY4QL17pwT2el3lzci_u2xWopedQlffwOdvmJ8qyEgJlnWhy6yRqmtK9IZFSQf6xi2CY8weffkcSwYeRqZpB53WSIZxwKMhnxLUFyIZOl4bIQC2-3bsEYm1q7WKENW_kiBjD7uxg/s1600/SZ3a.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhMVJhKY4QL17pwT2el3lzci_u2xWopedQlffwOdvmJ8qyEgJlnWhy6yRqmtK9IZFSQf6xi2CY8weffkcSwYeRqZpB53WSIZxwKMhnxLUFyIZOl4bIQC2-3bsEYm1q7WKENW_kiBjD7uxg/s320/SZ3a.png" /></a></div><br />
Altogether these methods are meant to be used to examine the characteristics of the time series representation of the data. The first Sims and Zha method would most likely be the one to publish in a paper, but the other two methods give the researcher a more complete picture of the nature of the posterior distribution.<br />
<br />
A few small things still need to be worked out for the Sims and Zha methods to be a complete part of the VAR package. As I mentioned before, they still need to be implemented for orthogonalized IRFs, but this will not be difficult (they are very clear in their paper in moving towards this implementation). Also, it will be important to bring in user choice when it comes to which principal component to use when characterizing the data. For Sims and Zha error bands 1 and 2, the user can pass in a matrix of integers that correspond to the chosen principal components of the variance-covariance matrix.<br />
<br />
Scheduling here on out (let me know what you think, updated):<br />
<br />
1) I'd like to completely integrate Sims and Zha error bands. Specifically, this means:<br />
<br />
a) Component choice for SZ3<br />
<br />
b) Orthogonalized error bands<br />
<br />
c) Clean up code<br />
<br />
I will aim to finish the above tasks by this Saturday (7/23).<br />
<br />
2) After this, I would like to move on to structural VAR implementation. This in itself will feed back on the error band methods that I have been working on. SVAR will draw from the ML methods already present in statsmodels. I'd like to finish a SVAR estimation method by the August 3rd.<br />
<br />
3) Once SVAR is completely integrated into the package, per Skipper's suggestion, I will be using pure Python to generalize the Kalman filter. I'll have more questions once I reach that point.<br />
<br />
I'd really like to move much quicker than I did in the first half of GSOC and hit these goals.<br />
<br />
Bart<br />
<br />
<br />
<br />
<br />
<br />
<hr></hr><br />
<sup id="fn1">1. SZ suggest using the largest eigenvalue(s) of Lambda, as they will most likely identify the majority of the variation in this type of data.<a href="#ref1" title="Jump back to footnote 1 in the text.">↩</a></sup>Barthttp://www.blogger.com/profile/12636757026304914432noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-50228116820853625132011-06-11T14:01:00.001-04:002011-06-27T12:06:39.310-04:00Monte Carlo Standard Errors for Impulse Response FunctionsI've come to a major checkpoint in integrating Monte Carlo error bands<br />
for impulse response functions (this is only non-orthogonal right now).<br />
<br />
Here is some quick code to get VAR IRF MC standard errors:<br />
<br />
<pre><div class="mycode">import numpy as np
import scikits.statsmodels.api as sm
from scikits.statsmodels.tsa.api import VAR
import matplotlib.pyplot as plt
mdata = sm.datasets.macrodata.load().data
mdata = mdata[['realgdp','realcons','realinv']]
#mdata = mdata[['realgdp','realcons','realinv', 'pop', 'realdpi', 'tbilrate']]
names = mdata.dtype.names
data = mdata.view((float,3))
data = np.diff(np.log(data), axis=0)
mymodel = VAR(data,names=names)
res = mymodel.fit(maxlags=3,ic=None)
#first generate asymptotic standard errors (to compare)
res.irf(periods=30).plot(orth=False, stderr_type='asym')
#then generate monte carlo standard errors
res.irf(periods=30).plot(orth=False, stderr_type='mc', seed=29, repl=1000)
plt.show()
</div></pre><br />
Produces the following plots of a shock to (delta) realgdp.<br />
<br />
Asymptotic:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhi6Lnvv6JGBNHxTV7i-bz4mJEcaqGwMdaXh18D7cSq8XDt5iKjDxro4LGgGXs4j0G_ZIHNaX7y9jAKtBmC6p93tqij_UZ7sjyQ8FHwLfxM-56_UFT11qjYgM31OspMG032ev8jKv_rsWk/s1600/asym.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhi6Lnvv6JGBNHxTV7i-bz4mJEcaqGwMdaXh18D7cSq8XDt5iKjDxro4LGgGXs4j0G_ZIHNaX7y9jAKtBmC6p93tqij_UZ7sjyQ8FHwLfxM-56_UFT11qjYgM31OspMG032ev8jKv_rsWk/s320/asym.png" /></a></div><br />
Monte Carlo:<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj274CsqjYXiN3EbWQYOwpWY2CtiVrA1X7eR372YHxjsOTh2mAewLx7NCBFC6LcBFI9aNK7yf25oeyG6nqixUMYdCM7LTxe3nhX4YKkjlecghhbalabDYCFgSWGie2QkUzUSyehHYbKxDU/s1600/mc.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj274CsqjYXiN3EbWQYOwpWY2CtiVrA1X7eR372YHxjsOTh2mAewLx7NCBFC6LcBFI9aNK7yf25oeyG6nqixUMYdCM7LTxe3nhX4YKkjlecghhbalabDYCFgSWGie2QkUzUSyehHYbKxDU/s320/mc.png" /></a></div><br />
I added functions to the VARResults, IRAnalysis, and also added to<br />
some of the pre-existing functions in these classes and also to<br />
plotting.py. Because Monte Carlo standard errors are in general not<br />
symmetric, I had to alter the plot_with_error function from the<br />
tsa.vector_ar.plotting.py file and number of other functions.<br />
<br />
Functions added:<br />
<br />
VARResults.stderr_MC_irf(self, orth=False, repl=1000, T=25,<br />
signif=0.05, seed=None)<br />
This function generates a tuple that holds the lower and upper error<br />
bands generated by Monte Carlo simulations.<br />
<br />
IRAnalysis.cov_mc(self, orth=False, repl=1000, signif=0.05, seed=None)<br />
This just calls the stderr_MC_irf function from the original model<br />
using the number of periods specified when the irf() class is called<br />
from the VAR class.<br />
<br />
Modified functions:<br />
<br />
BaseIRAnalysis.plot<br />
Added specification of error type stderr_type='asym' or 'mc' (see<br />
example). Also repl (replications) and seed options added if 'mc'<br />
errors are specified.<br />
<br />
in tsa.vector_ar.plotting.py<br />
plot_with_error()<br />
irf_grid_plot()<br />
<br />
These functions now take the error type specified in the plot()<br />
function above and treats the errors accordingly. Because of how all<br />
of the VAR work (esp. plotting) assumed asymptotic standard errors,<br />
all of these irf plot functions assumed that errors were passed in as<br />
a single matrix with each standard error depending on the MA lag<br />
length and shock variable. Now, depending on which error is specified,<br />
the function will take a tuple of arrays as the standard error if<br />
(stderr_type='mc') rather than a single array.<br />
<br />
A serious issue right now is speed. While the asymptotic standard<br />
errors take about a half second on my home laptop to run, the Monte<br />
Carlo standard errors with 3 variables and 1000 replications takes<br />
about 13 seconds to run. Each replication discards the first 100<br />
observations. The most taxing aspect of generating the errors is<br />
re-simulating the data assuming normally distributed standard errors<br />
1000 times (using the util.varsim() function).<br />
<br />
BartBarthttp://www.blogger.com/profile/12636757026304914432noreply@blogger.com3tag:blogger.com,1999:blog-128274497687662608.post-66197143009561312432011-05-23T09:55:00.001-04:002011-05-24T10:54:31.357-04:005/6 - 5/20: Getting Acclimated and VARRresults.reorder() functionJust finished up the official first two weeks of GSoC. I spent most of the first week attempting to get a feel for how all of the time series methods for statsmodels are organized. Once I felt comfortable within the VAR package, I began work on adding a reorder() function to the VARResults class, which allows the user to specify the order of the endogenous variables in a vector auto-regression system. While the order of the variables plays no role in estimating the system, if the shocks are to be identified, variable order is used to specify the within period impact of shocks to individual variables in the system.<br />
<br />
For example, let us say that we have a 3-variable VAR system, originally ordered realGDP, realcons, and realinv, contained within the VARResults class 'res'. Reordering the variables in the system is as simple as follows:<br />
<br />
<pre><div class="mycode">In [1]: import scikits.statsmodels.api as sm
In [2]: import numpy as np
In [3]: mdata = sm.datasets.macrodata.load().data
In [4]: mdata = mdata[['realgdp','realcons','realinv']]
In [5]: names = mdata.dtype.names
In [6]: data = mdata.view((float,3))
In [7]: from scikits.statsmodels.tsa.api import VAR
In [8]: res = VAR(data, names=names).fit(maxlags=3,ic=None)
In [9]: res.names
Out[9]: ['realgdp', 'realcons', 'realinv']
In [10]: res_re = res.reorder(['realinv','realcons','realgdp'])
In [11]: res_re.names
Out[11]: ['realinv', 'realcons', 'realgdp']
</div></pre><br />
The reorder function reuses all of the results from the original VAR class, but rearranges them to be in line with the new system. If working with large a # of observations, the computational advantage becomes useful pretty quickly. For example, with a 100,000 observation system with three variables, re-estimating the system after changing the variable order took 3.37 seconds, while using the reorder function took 0.57 seconds.<br />
<br />
In the next few weeks I am planning on adding more impulse response function error band estimation methods. The current package only includes analytical error bands.Barthttp://www.blogger.com/profile/12636757026304914432noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-15845399984868704422011-04-23T15:18:00.003-04:002011-04-25T09:11:35.318-04:00Statsmodels: GSoC PrelimThis is my first entry in my Statsmodels Project Summer 2011 blog. I will update this blog weekly as the summer goes by with information regarding progress on my work for the scikits.statsmodels Python library.<br />
<br />
In order to complete the preparation process for the statsmodels Google Summer of Code sponsorship, I wrote a quick patch that included a cointegration test. As of now, the test can only be run on a bivariate system with a simple Dickey-Fuller test on the residuals, using the Mackinnon [1] critical values. I would like to expand the functionality of this test to allow for Augmented Dickey-Fuller tests of the residuals and also tests of multivariate cointegrated systems. This patch served as a nice way to dive into what time series methods scikits.statsmodels currently includes in its toolbox and where to go from here.Barthttp://www.blogger.com/profile/12636757026304914432noreply@blogger.com2tag:blogger.com,1999:blog-128274497687662608.post-54353991881658281832010-06-22T10:42:00.000-04:002010-06-22T10:42:19.252-04:00Statsmodels: GSoC Week 4 UpdateI spent the last week finishing up the paper that I submitted to accompany my talk at the <a href="http://conference.scipy.org/scipy2010/">SciPy conference</a>. I am really looking forward to going to Austin and hearing all the great talks (plus I hear the beer is cheap and the food and music are good, which doesn't hurt). In addition to finishing up the paper, I have started to clean up our time series code.<br />
<br />
So far this has included finishing the augmented Dickey-Fuller (ADF) test for unit roots. The big time sink here is that the ADF test-statistic has a non-standard distribution in most cases. The ADF test statistic is obtained by running the following regression<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=%5CDelta%20y_%7Bt%7D%20=%20%5Calpha@plus;%5Cbeta%20t@plus;%5Cgamma%20y_%7Bt-1%7D@plus;%5Cdelta_%7B1%7D%5CDelta%20y_%7Bt-1%7D%20@plus;%20%5Ccdots%20@plus;%5Cdelta_%7Bp%7D%5CDelta%20y_%7Bt-p%7D" target="_blank"><img src="http://latex.codecogs.com/gif.latex?%5CDelta%20y_%7Bt%7D%20=%20%5Calpha+%5Cbeta%20t+%5Cgamma%20y_%7Bt-1%7D+%5Cdelta_%7B1%7D%5CDelta%20y_%7Bt-1%7D%20+%20%5Ccdots%20+%5Cdelta_%7Bp%7D%5CDelta%20y_%7Bt-p%7D" title="\Delta y_{t} = \alpha+\beta t+\gamma y_{t-1}+\delta_{1}\Delta y_{t-1} + \cdots +\delta_{p}\Delta y_{t-p}" /></a><br />
<br />
One approach to testing for a unit root means testing the t-stat on the coefficient on the lagged level of <i>y</i>. The actual distribution for this statistic, however, is not Student's t. Many software packages use the tables in Fuller (1976, updated to 1996 version or not) in order to get the critical values for the test statistic depending on the sample size. They use linear interpolation for sample sizes not included in the table. The p-values for the obtained test statistic are usually obtained using MacKinnon's (1994) study that estimated regression surfaces of these distributions via Monte Carlo simulation.<br />
<br />
While we do use MacKinnon's approximate p-values from the 1994 paper, MacKinnon wrote a note updating this paper in early 2010, which gives new regression surface results for obtaining the critical values. We use these new results for the critical values. Therefore, when using our ADF test, it is advised that if the p-value is close to the reject/accept region then the critical values should be used in place of the p-value to make the ultimate decision.<br />
<br />
We can illlustrate the use of ADF. Note that this version is only in my branch and that it is still in the sandbox, even though it has now been tested, because the API and returned results may change. We will demonstrate on a series that we can easily guess is non-stationary, real GDP.<br />
<br />
<br />
<pre><div class="mycode">In [1]: import scikits.statsmodels as sm
In [2]: from scikits.statsmodels.sandbox.tsa.stattools import adfuller
In [3]: data = sm.datasets.macrodata.load()
In [4]: realgdp = data.data['realgdp']
In [5]: adf = adfuller(realgdp, maxlag=4, autolag=None, regression="ct")
In [6]: adf
Out[6]:
(-1.8566384063254346,
0.67682917510440099,
4,
198,
{'1%': -4.0052351400496136,
'10%': -3.1402115863254525,
'5%': -3.4329000694218998})
</div></pre><br />
<br />
The return values are the test statistic, its p-value (the null-hypothesis here is that the series <i>does</i> contain a unit root), the number of lags of the differences used, the number of observations for the regression, and a dictionary containing the critical values at the respective confidence levels. The regression option controls the type of regression (ie., whether to include a constant or a linear or quadratic time trend), and the autolag option has three options for choosing the lag length to help correct for serial correlation in the regression. There are 'AIC', 'BIC', and 't-stat'. The former two choose the lag length that maximizes the infofrmation criterion, the latter chooses the lag length based on the significance of the lag. This starts with maxlag and works its way down. The docstring has more detailed information.<br />
<br />
<br />
<br />
<br />
Beyond this, I have been working on an autocorrelation function (acf), a partial autocorrelation function (pacf), and Q-Statistics (Box-Ljung test). Next up for this week is finishing my VAR class with identification schemes. After this, I will work to integrate post-estimation tests into our results classes, most likely using some sort of mix-in classes and attach test containers to the results objects for test results. Then it's off to the SciPy conference. There I will hopefully be participating in the stats sprint, helping out with the docs marathon and discussing what we need for the future of statistics and Python.<br />
<br />
<br />
<pre><div class="mycode">Fuller, W.A. 1996. <i>Introduction to Statistical Time Series.</i> 2nd ed. Wiley.
MacKinnon, J.G. 1994. "Approximate asymptotic distribution functions for
unit-root and cointegration tests. <i>Journal of Business and Economic
Statistics</i> 12, 167-76.
MacKinnon, J.G. 2010. "Critical Values for Cointegration Tests."
Queen's University, Dept of Economics, Working Papers. Available at
http://ideas.repec.org/p/qed/wpaper/1227.html
</div></pre>jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com2tag:blogger.com,1999:blog-128274497687662608.post-37748952510873894002010-06-11T19:11:00.004-04:002010-06-14T11:54:36.211-04:00Statsmodels: GSoC Week 3 Update[<b>Edit</b>: Formatting should be fixed now. I will not be reformatting old posts though, so that they don't get reposted at <a href="http://planet.scipy.org">Planet SciPy</a>]<br /><br /><br />
<br />
Last week was spent mainly ensuring that I pass my comps and remain a PhD student. This week was much more productive for coding. For now, all changes are in my branch and have not been merged to trunk, but I will describe the two big changes.<br /><br /><br />
<br />
The first concerns the datasets package. This one is not all that exciting, but suffice it to say that the datasets are now streamlined and use the <a href="http://code.activestate.com/recipes/52308/">Bunch pattern</a> to load the data. Thanks, Gaël, for pointing this out. I also rewrote a bit of David's datasets proposal from scikits-learn to reflect the current design of our datasets and thoughts. You can see it <a href="http://bazaar.launchpad.net/%7Ejsseabold/statsmodels/statsmodels-skipper/annotate/head:/scikits/statsmodels/datasets/DATASET_PROPOSAL.rst">here</a> (soon to be on the docs page). We are making an effort to ensure that our datasets are going to be similar to those of scikits-learn.<br /><br /><br />
<br />
The second change was an improvement of the fitting of maximum likelihood models and the start of a GenericLikelihoodModel class. Maximum likelihood based models (mainly discrete choice models in the main code base right now) can now be fit using any of the unconstrained solvers from scipy.optimize (Nelder-Mead, BFGS, CG, Newton-CG, Powell) plus Newton-Raphson. To take a simple example to see how it works, we can fit a Probit model.<br /><br /><br />
<br />
<div class = "mycode">In [1]: import scikits.statsmodels as sm<br />
<br />
In [2]: data = sm.datasets.spector.load()<br />
<br />
In [3]: data.exog = sm.add_constant(data.exog)<br />
<br />
In [4]: res_newton = sm.Probit(data.endog, data.exog).fit(method="newton")<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations 6<br />
<br />
In [5]: res_nm = sm.Probit(data.endog, data.exog).fit(method="nm", maxiter=500)<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations: 439<br />
Function evaluations: 735<br />
<br />
In [6]: res_bfgs = sm.Probit(data.endog, data.exog).fit(method="bfgs")<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations: 15<br />
Function evaluations: 21<br />
Gradient evaluations: 21<br />
<br />
In [7]: res_cg = sm.Probit(data.endog, data.exog).fit(method="cg", maxiter=250)<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations: 188<br />
Function evaluations: 428<br />
Gradient evaluations: 428<br />
<br />
In [8]: res_ncg = sm.Probit(data.endog, data.exog).fit(method="ncg", avextol=1e-8)<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations: 12<br />
Function evaluations: 14<br />
Gradient evaluations: 12<br />
Hessian evaluations: 12<br />
<br />
In [9]: res_powell = sm.Probit(data.endog, data.exog).fit(method="powell", ftol=1e-8)<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations: 12<br />
Function evaluations: 568</div><br /><br /><br />
<br />
All of the options for the solvers are available and are documented in the fit method. As you can see, some of the default values need to be changed to ensure (accurate) convergence. The Results objects that are returned have two new attributes.<br /><br /><br />
<br />
<div class = "mycode"><br />
In [10]: res_powell.mle_retvals<br />
Out[10]: <br />
{'converged': True,<br />
'direc': array([[ 7.06629660e-02, -3.07499922e-03, 5.38418734e-01,<br />
-4.19910465e-01],<br />
[ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00,<br />
0.00000000e+00],<br />
[ 1.49194876e+00, -6.64992809e-02, -6.96792443e-03,<br />
-3.22306873e+00],<br />
[ -5.36227277e-02, 1.18544093e-01, -8.75205765e-02,<br />
-2.42149981e+00]]),<br />
'fcalls': 568,<br />
'fopt': 12.818804069990534,<br />
'iterations': 12,<br />
'warnflag': 0}<br />
<br />
In [11]: res_powell.mle_settings<br />
Out[11]: <br />
{'callback': None,<br />
'disp': 1,<br />
'fargs': (),<br />
'ftol': 1e-08,<br />
'full_output': 1,<br />
'maxfun': None,<br />
'maxiter': 35,<br />
'optimizer': 'powell',<br />
'retall': 0,<br />
'start_direc': None,<br />
'start_params': [0, 0, 0, 0],<br />
'xtol': 0.0001}</div><br />
The dict mle_retvals contains all of the values that are returned from the solver if the full_output keyword is True. The dict mle_settings contains all of the arguments passed to the solver, including the defaults so that these can be checked after the fit. Again, all settings and returned values are documented in the fit method and in the results class, respectively.<br /><br /><br />
<br />
Lastly, I started a GenericLikelihoodModel class. This is currently unfinished, though the basic idea is laid out. Take again the Probit example above using Lee Spector's <a href="http://www.jstor.org/pss/1182446">educational program data</a>. And assume we didn't have the Probit model from statsmodels. We could use the new GenericLikelihoodModel class. There are two ways (probably more) to proceed. For those comfortable with object oriented programming and inheritance in Python, we could subclass the GenericLikelihoodModel, defining our log-likelihood method.<br /><br /><br />
<br />
<div class="mycode">from scikits.statsmodels import GenericLikelihoodModel as LLM<br />
from scipy import stats<br />
import numpy as np<br />
<br />
class MyModel(LLM):<br />
def loglike(self, params):<br />
"""<br />
Probit log-likelihood<br />
"""<br />
q = 2*self.endog - 1<br />
X = self.exog<br />
return np.add.reduce(stats.norm.logcdf(q*np.dot(X,params)))</div><br /><br /><br />
<br />
Now this model could be fit, using any of the methods that only require an objective function, i.e., Nelder-Mead or Powell.<br /><br /><br />
<br />
<div class="mycode">In [43]: mod = MyModel(data.endog, data.exog)<br />
<br />
In [44]: res = mod.fit(method="nm", maxiter=500)<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations: 439<br />
Function evaluations: 735<br />
<br />
In [45]: res_nm.params<br />
Out[45]: array([ 1.62580058, 0.05172931, 1.42632242, -7.45229725])<br />
<br />
In [46]: res.params<br />
Out[46]: array([ 1.62580058, 0.05172931, 1.42632242, -7.45229725])</div><br /><br /><br />
<br />
The main drawback right now is that all statistics that rely on the covariance of the parameters, etc. will use numeric gradients and Hessians, which can lessen that accuracy of those statistics. This can be overcome by providing score and hessian methods as loglike was provided above. Of course, for more complicated likelihood functions this can soon become cumbersome. We are working towards more accurate numerical differentiation and discussing options for automatic or symbolic differentiation.<br /><br /><br />
<br />
The main advantage as opposed to just writing your likelihood and passing it to a solver is that you have all of the (growing number of) statistics and tests available to statsmodels right in the generic model.<br /><br /><br />
<br />
I would also like to accommodate those who are less familiar with OOP and inheritance in Python. I haven't quite worked out the final design for how this would go yet. Right now, you could do the following, though I don't think it quite meets the less complicated goal.<br /><br />
<br />
<div class="mycode">In [4]: from scikits.statsmodels.model import GenericLikelihoodModel as LLM<br />
<br />
In [5]: import scikits.statsmodels as sm<br />
<br />
In [6]: from scipy import stats<br />
<br />
In [7]: import numpy as np<br />
<br />
In [8]: <br />
<br />
In [9]: data = sm.datasets.spector.load()<br />
<br />
In [10]: data.exog = sm.add_constant(data.exog)<br />
<br />
In [11]: <br />
<br />
In [12]: def probitloglike(params, endog, exog):<br />
....: """<br />
....: Log likelihood for the probit<br />
....: """<br />
....: q = 2*endog - 1<br />
....: X = exog<br />
....: return np.add.reduce(stats.norm.logcdf(q*np.dot(X,params)))<br />
....: <br />
<br />
In [13]: mod = LLM(data.endog, data.exog, loglike=probitloglike)<br />
<br />
In [14]: res = mod.fit(method="nm", fargs=(data.endog,data.exog), maxiter=500)<br />
Optimization terminated successfully.<br />
Current function value: 12.818804<br />
Iterations: 439<br />
Function evaluations: 735<br />
<br />
In [15]: res.params<br />
Out[15]: array([ 1.62580058, 0.05172931, 1.42632242, -7.45229725])</div><br /><br /><br />
<br />
There are still a few design issues and bugs that need to be worked out with the last example, but the basic idea is there. That's all for now.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com4tag:blogger.com,1999:blog-128274497687662608.post-5931209765639074392010-05-31T10:33:00.000-04:002010-05-31T10:33:34.820-04:00Week 1 GSoC UpdateLast week was the first of the Google Summer of Code. I spent most of the week in a Bayesian econometrics class led by <a href="http://www.amazon.com/gp/product/0471679321?ie=UTF8&ref_=ntt_at_ep_dpt_1">John Geweke</a> and studying for a comprehensive exam that I take this week, so progress on statsmodels was rather slow. That said, I have been able to take care of some low hanging fruit.<br />
<br />
There are a few name changes:<br />
<br />
statsmodels/family -> statsmodels/families<br />
statsmodels/lib/io.py -> statsmodels/iolib/foreign.py<br />
<br />
Also Vincent has done a good bit of work on improving our output using the SimpleTable class from <a href="http://code.google.com/p/econpy/">econpy</a>. I will post some examples over the coming weeks, but SimpleTable provides an easy way to make tables in ASCII text, HTML, or LaTeX. The SimpleTable class has been moved<br />
<br />
statsmodels/sandbox/ouput.py -> statsmodels/iolib/table.py<br />
<br />
Beyond the renames, I have removed the soft dependency on RPy for running our tests in favor of hard-coded results, refactored our tests, and added a few additional ones along the way.<br />
<br />
We are also making an effort to keep our <a href="http://statsmodels.sourceforge.net/">online documentation</a> synced with the current trunk. The biggest change to our documentation is the addition of a <a href="http://statsmodels.sourceforge.net/developernotes.html">developer's page</a> for those who might like to get involved. As always, please report problems with the docs on either the <a href="http://www.scipy.org/Mailing_Lists">scipy-user list</a> or join in the discussions of statsmodels, pandas, larry, and other topics on statistics and Python at the <a href="http://groups.google.ca/group/pystatsmodels">pystatsmodels Google group</a>.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com1tag:blogger.com,1999:blog-128274497687662608.post-64329066440641825842010-05-01T12:59:00.000-04:002010-05-01T12:59:10.053-04:00Plans for the SummerA quick update on the plans for statsmodels over the next few months.<br />
<br />
I have been accepted for my second Google Summer of Code, which means that we will have a chance to make a big push to get a lot of our work out of the sandbox, tested, and included in the main code base.<br />
<br />
You can see the roadmap on Google's GSoC site <a href="http://socghop.appspot.com/gsoc/student_proposal/show/google/gsoc2010/jseabold/t127082935700">here</a>. You might have to log in to view it. <br />
<br />
The quick version follows. As far as general issues, I will be getting the code ready for Python 3 and focusing on some design issues including an improved generic maximum likelihood framework, post-estimation testing, variable name handling, and output in text tables, LaTeX, and html. I will then be working to get a lot of our code out of the sandbox. This includes timeseries convenience functions and models such as GARCH, VARMA, Hodrick-Prescott filter, and a state space model that uses the Kalman filter. I will be polishing the systems of equation framework and panel (longitudinal) data estimators. We have also been working on some nonparametric estimators including univariate kernel density estimators and kernel regression estimators. Finally, as part of my coursework I have been working toward (generalized) maximum entropy models that I hope to include as well as some work on the scipy.maxentropy module.<br />
<br />
I will give a quick talk on the project for the SciPy Conference in Austin.<br />
<br />
It looks like we are set to make a good deal of progress on the code this summer.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-87061377107681992912010-03-08T19:37:00.003-05:002010-03-08T22:06:55.731-05:00Sparse Least Squares SolverI have a homework doing some monte carlo experiments of an autoregressive process of order 1, and I thought I would use it as an example to demonstrate the sparse least squares solver that Stefan committed to scipy revision 6251.<br />
<br />
All mistakes are mine...<br />
<br />
Given an AR(1) process<br />
<br />
<a href="http://www.codecogs.com/eqnedit.php?latex=y_{t}=.9y_{t-1} @plus; .01\epsilon_{t} \newline \newline \text{where} \newline \newline \epsilon_{t} \sim \boldsymbol{\text{NID} }\left(0,1\right)" target="_blank"><img src="http://latex.codecogs.com/gif.latex?y_{t}=.9y_{t-1} + .01\epsilon_{t} \newline \newline \text{where} \newline \newline \epsilon_{t} \sim \boldsymbol{\text{NID} }\left(0,1\right)" title="y_{t}=.9y_{t-1} + .01\epsilon_{t} \newline \newline \text{where} \newline \newline \epsilon_{t} \sim \boldsymbol{\text{NID} }\left(0,1\right)" /></a><br />
<br />
<br />
We can estimate the following autoregressive coefficients.<br />
<br />
<div = "mycode">In [1]: import numpy as np<br />
<br />
In [2]: np.set_printoptions(threshold=25)<br />
<br />
In [3]: np.random.seed(1)<br />
<br />
In [4]: # make 1000 autoregressive series of length 100<br />
<br />
In [5]: # y_0 = 0 by assumption<br />
<br />
In [6]: samples = np.zeros((100,1000))<br />
<br />
In [7]: for i in range(1,100):<br />
...: error = np.random.randn(1,1000)<br />
...: samples[i] = .9 * samples[i-1] + .01 * error<br />
...: <br />
<br />
In [8]: from scipy import sparse<br />
<br />
In [9]: from scipy.sparse import linalg as spla<br />
<br />
In [10]: # make block diagonal sparse matrix of y_t-1<br />
<br />
In [11]: # recommended to build as a linked list<br />
<br />
In [12]: spX = sparse.lil_matrix((99*1000,1000))<br />
<br />
In [13]: for i in range(1000):<br />
....: spX[i*99:(i+1)*99,i] = samples[:-1,i][:,None]<br />
....: <br />
<br />
In [14]: spX = spX.tocsr() # convert it to csr for performance<br />
<br />
In [15]: # do the least squares<br />
<br />
In [16]: retval = spla.isolve.lsqr(spX, samples[1:,:].ravel('F'), calc_var=True) <br />
In [17]: retval[0]<br />
Out[17]: <br />
array([ 0.88347438, 0.8474124 , 0.85282674, ..., 0.91019165,<br />
0.89698465, 0.76895806])<br />
</div><br />
I'm curious if there's any downside to using sparse least squares whenever the RHS of a least squares can be written in block diagonal form.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com2tag:blogger.com,1999:blog-128274497687662608.post-49505206528934787212010-02-16T14:23:00.001-05:002010-02-16T14:24:38.969-05:00scikits.statsmodels 0.2.0 releaseWhile I find no time to blog, I thought I'd post our newest release announcement here.<br />
<br />
We are happy to announce the 0.2.0 (beta) release of scikits.statsmodels. This is both a bug-fix and new feature release.<br />
<br />
Download<br />
------------<br />
<br />
You can easy_install (or PyPI URL:<br />
<a href="http://pypi.python.org/pypi/scikits.statsmodels/">http://pypi.python.org/pypi/scikits.statsmodels/</a>)<br />
<br />
Source downloads: <a href="http://sourceforge.net/projects/statsmodels/">http://sourceforge.net/projects/statsmodels/</a><br />
<br />
Development branches: <a href="http://code.launchpad.net/statsmodels">http://code.launchpad.net/statsmodels</a><br />
<br />
Note that the trunk branch on launchpad is almost always stable and has the most up to date changes since our releases are so few and far between.<br />
<br />
Documentation<br />
------------------<br />
<br />
<a href="http://statsmodels.sourceforge.net/">http://statsmodels.sourceforge.net/</a><br />
<br />
We invite you to install, kick the tires, and make bug reports and feature requests.<br />
<br />
Feedback can either be on scipy-user or the mailing list at<br />
<a href="http://groups.google.com/group/pystatsmodels?hl=en">http://groups.google.com/group/pystatsmodels?hl=en</a><br />
Bug tracker: <a href="https://bugs.launchpad.net/statsmodels">https://bugs.launchpad.net/statsmodels</a><br />
<br />
Main Changes in 0.2.0<br />
---------------------------------<br />
<br />
* Improved documentation and expanded and more examples<br />
* Added four discrete choice models: Poisson, Probit, Logit, and Multinomial Logit.<br />
* Added PyDTA. Tools for reading Stata binary datasets (*.dta) and putting them into numpy arrays.<br />
* Added four new datasets for examples and tests.<br />
* Results classes have been refactored to use lazy evaluation.<br />
* Improved support for maximum likelihood estimation.<br />
* bugfixes<br />
* renames for more consistency<br />
RLM.fitted_values -> RLM.fittedvalues<br />
GLMResults.resid_dev -> GLMResults.resid_deviance<br />
<br />
Sandbox<br />
-------------<br />
<br />
We are continuing to work on support for systems of equations models, panel data models, time series analysis, and information and entropy econometrics in the sandbox. This code is often merged into trunk as it becomes more robust.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-4390660609869315152009-08-31T21:58:00.003-04:002009-08-31T22:01:23.914-04:00scikits.statsmodels Release AnnouncementWe have been working hard to get a release ready for general consumption for the statsmodels code. Well, we're happy to announce that a (very) beta release is ready.<br />
<br />
Background<br />
==========<br />
<br />
The statsmodels code was started by Jonathan Taylor and was formerly included as part of scipy. It was taken up to be tested, corrected, and extended as part of the Google Summer of Code 2009.<br />
<br />
What it is<br />
==========<br />
<br />
We are now releasing the efforts of the last few months under the scikits namespace as scikits.statsmodels. Statsmodels is a pure python package that requires numpy and scipy. It offers a convenient interface for fitting parameterized statistical models with growing support for displaying univariate and multivariate summary statistics, regression summaries, and (postestimation) statistical tests.<br />
<br />
Main Feautures<br />
==============<br />
<br />
* regression: Generalized least squares (including weighted least squares and least squares with autoregressive errors), ordinary least squares.<br />
* glm: Generalized linear models with support for all of the one-parameter exponential family distributions.<br />
* rlm: Robust linear models with support for several M-estimators.<br />
* datasets: Datasets to be distributed and used for examples and in testing.<br />
<br />
There is also a sandbox which contains code for generalized additive models (untested), mixed effects models, cox proportional hazards model (both are untested and still dependent on the nipy formula framework), generating descriptive statistics, and printing table output to ascii, latex, and html. None of this code is considered "production ready".<br />
<br />
Where to get it<br />
===============<br />
<br />
Development branches will be on LaunchPad. This is where to go to get the most up to date code in the trunk branch. Experimental code will also be hosted here in different branches.<br />
<br />
<a href="https://code.launchpad.net/statsmodels">https://code.launchpad.net/statsmodels</a><br />
<br />
Source download of stable tags will be on SourceForge.<br />
<br />
<a href="https://sourceforge.net/projects/statsmodels/">https://sourceforge.net/projects/statsmodels/</a><br />
<br />
or<br />
<br />
PyPi: <a href="http://pypi.python.org/pypi/scikits.statsmodels/">http://pypi.python.org/pypi/scikits.statsmodels/</a><br />
<br />
License<br />
=======<br />
<br />
Simplified BSD<br />
<br />
Documentation<br />
=============<br />
<br />
The official documentation is hosted on SourceForge.<br />
<br />
<a href="http://statsmodels.sourceforge.net/">http://statsmodels.sourceforge.net/</a><br />
<br />
The sphinx docs are currently undergoing a lot of work. They are not yet comprehensive, but should get you started.<br />
<br />
This blog will continue to be updated as we make progress on the code.<br />
<br />
Discussion and Development<br />
==========================<br />
<br />
All chatter will take place on the or <a href="http://www.scipy.org/ Mailing_Lists">scipy-user mailing list</a>. We are very interested in receiving feedback about usability, suggestions for improvements, and bug reports via the mailing list or the bug tracker at <a href="https://bugs.launchpad.net/statsmodels">https://bugs.launchpad.net/statsmodels</a>.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com3tag:blogger.com,1999:blog-128274497687662608.post-79176515449739929862009-08-27T00:15:00.004-04:002009-08-27T19:59:32.678-04:00GSoC Is OverWhoa, where did the last month go? The Google Summer of Code 2009 officially ended this Monday. Though I haven't taken a breath to update the blog, we (Josef and I) have been hard at work on the models code. <br /><br />We have working and tested versions of Generalized Least Squares, Weighted Least Squares, Ordinary Least Squares, Robust Linear Models with several M-estimators, and Generalized Linear Models with support for all (almost all?) one parameter exponential family distributions. We have also provided some more convenience functions, created a standalone python package for the models code, and obtained permissions to distribute a few more datasets. Due to a lack of time, there is only experimental (read untested) support for autoregressive models, mixed effects models, generalized additive models, and convenience functions for returning strings (possibly html and latex output as well) with regression results and descriptive statistics. I will continue to work on these as I find time.<br /><br />I will soon post a note on the progress that was made in the robust linear models code. Also, look out for a (semi-) official release of the code in the next few days. We have decided to name the project statsmodels and distribute it as a <a href="http://scikits.appspot.com/">scikit</a>. We need to finalize the documentation (should be ready to go in the next day or so...I am back taking courses) and clean up some of the usage examples, so people can jump right in and use the code, give feedback, and hopefully contribute extensions and new models. <br /><br />As for the future of statsmodels, we are discussing over the next few weeks the immediate extensions that we know would like to make. It's looking like I will be wearing my microeconometrician hat this semester in my own coursework. More specifically, I will probably be working with <a href="http://en.wikipedia.org/wiki/Cross-sectional_data">cross-sectional</a> and <a href="http://en.wikipedia.org/wiki/Panel_data">panel data</a> models for household survey data in my own research and finding some time for time series models as part of my teaching assistantship. Josef has also mentioned wanting to work more with time series models. <br /><br />If anyone (especially those from other disciplines) would like to contribute or see some extensions (my apologies to those who have made requests that I haven't yet been able to accomodate) feel free to post to the <a href="http://scipy.org/Mailing_Lists">scipy-dev mailing list</a>. I'm more than happy to discuss/debate with users and potential developers the design decisions that have been made, as I think the code is still in an unsettled enough state to merit some discussion.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com4tag:blogger.com,1999:blog-128274497687662608.post-37254581089014995972009-07-25T23:22:00.007-04:002009-07-26T15:10:17.075-04:00Iterated Reweighted Least SquaresI have spent the last two weeks putting the "finishing" touches on the generalized linear models and starting to go over the robust linear models (RLM). The test suite for GLM is not complete yet, but all of the exponential families are covered with at least their default link functions tested and are looking good. So in an effort to make a first pass over all of the existing code, I moved on to RLM. After the time spent with GLM theory, the RLMs theory and code was much more manageable. <br /><br />Before discussing the RLMs, their implementation, and the extensions I have made. I will describe the iterated reweighted least squares (IRLS) algorithm for the GLMs to demonstrate the theory and the solution method in the models code. A very similar iteration is done for the RLMs as well.<br /><br />The main idea of GLM, as noted, is to relate a response variable to a linear model via a link function, which allows us to use least squares regression. Let us take as an example, the binomial family (which is written to handle Bernoulli and binomial data). In this case, the response variable is Bernoulli, 1 indicates a "success" and 0 a "failure".<br /><br />The default link for the binomial family is the logit link. So we have<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=\eta=g(u)=\ln\frac{\mu}{1-\mu}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\eta=g(u)=\ln\frac{\mu}{1-\mu}" title="\eta=g(u)=\ln\frac{\mu}{1-\mu}" /></a><br /><br /><em/>η</em> is our linear predictor and <em>μ</em> is our actual mean response. The first thing that we need for the algorithm is to compute a first guess on <em>μ</em> (IRLS as opposed to Newton-Raphson makes a guess on the mean response rather than the parameter vector <em>β</em>). The algorithm is fairly robust to this first guess; however, a common choice is<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=\mu=\frac{y@plus;\bar{y}}{2}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\mu=\frac{y+\bar{y}}{2}" title="\mu=\frac{y+\bar{y}}{2}" /></a><br /><br />For the binomial family, we specifically use<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=\mu=\frac{y@plus;.5}{2}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\mu=\frac{y+.5}{2}" title="\mu=\frac{y+.5}{2}" /></a><br /><br />where <em>y</em> is our given response variable. We then use this first guess to initialize our linear predictor <em>η</em> via the link function given above. With these estimates we are able to start the iteration. Our convergence criteria is based on the deviance function, which is simply twice the log-likelihood ratio of our current guess on the fitted model versus the saturated log-likelihood. <br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=D=2\phi\left[ \mathcal{L}\left(y;y \right ) - \mathcal{L}\left(y;\mu \right )\right ]" target="_blank"><img src="http://latex.codecogs.com/gif.latex?D=2\phi\left[ \mathcal{L}\left(y;y \right ) - \mathcal{L}\left(y;\mu \right )\right ]" title="D=2\phi\left[ \mathcal{L}\left(y;y \right ) - \mathcal{L}\left(y;\mu \right )\right ]" /></a><br /><br />Where <em>Φ</em> is a dispersion (scale) parameter. Note that the saturated log-likelihood is simply the likelihood of the perfectly fitted model where <em>y</em> = <em>μ</em>. For the binomial family the deviance function is<br /> <br /><a href="http://www.codecogs.com/eqnedit.php?latex=\sum_{i}D_{i}=\begin{cases} -2\ln\left(1-\mu_{i}\right) & \text{, if }y_{i}=0\\ -2\ln\mu_{i} & \text{, if }y_{i}=1\end{cases}" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\sum_{i}D_{i}=\begin{cases} -2\ln\left(1-\mu_{i}\right) & \text{, if }y_{i}=0\\ -2\ln\mu_{i} & \text{, if }y_{i}=1\end{cases}" title="\sum_{i}D_{i}=\begin{cases} -2\ln\left(1-\mu_{i}\right) & \text{, if }y_{i}=0\\ -2\ln\mu_{i} & \text{, if }y_{i}=1\end{cases}" /></a><br /><br />The iteration continues while the deviance function evaluated at the updated <em>μ</em> differs from the previous by the given convergence tolerance (default is 1e-08) and the number of iterations is less than the given maximum (default is 100).<br /><br />The actual iterations, as the name of the algorithm suggests, run a weighted least squares fit of the actual regressors on the adjusted linear predictor (our transformed guess on the response variable). The adjustment is given by<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=\eta=\frac{\partial\eta}{\partial\mu}\left(y-\mu \right ) \\ \text{Recalling that }\eta=g(\mu)" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\eta=\frac{\partial\eta}{\partial\mu}\left(y-\mu \right ) \\ \text{Recalling that }\eta=g(\mu)" title="\eta=\frac{\partial\eta}{\partial\mu}\left(y-\mu \right ) \\ \text{Recalling that }\eta=g(\mu)" /></a><br /><br />which moves us closer to the root of the estimating equation (see the previously mentioned Gill or Hardin and Hilbe for the details of root-finding using the Newton-Raphson method. IRLS is simply Newton-Raphson using some simplifications.). The only remaining ingredient is the choice of weights for the weighted least squares. Similar to other methods such as general least squares (GLS) that correct for heteroskedasticity in the error terms, a diagonal matrix containing estimates of the variance of the error terms is used. However, in our case, the exact form of this variance is unknown and difficult to estimate, but the error each estimate is assumed to vary with the mean response variable. Thus, we improve the weights at each step given our best guess for the mean response variable <em>μ</em> and the known variance function of each family. For the binomial family, as is well known, the variance function is<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=V_{i}=\mu(1-\mu)" target="_blank"><img src="http://latex.codecogs.com/gif.latex?V_{i}=\mu(1-\mu)" title="V_{i}=\mu(1-\mu)" /></a><br /><br />Thus we can obtain an estimate of the parameters using weighted least squares. Using these estimates we update <em>η</em><br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=\eta=X\beta @plus; \left(\text{offset}\right)" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\eta=X\beta + \left(\text{offset}\right)" title="\eta=X\beta + \left(\text{offset}\right)" /></a><br /><br />Correcting by an offset if one is specified (the code currently does not support the use of offsets, but this would be a simple extension). Using this linear predictor (remember it was originally given by the linear transformation <em>η</em> = <em>g</em>(<em>μ</em>)) we update our guess on the mean response variable <em>μ</em> and use this to update our estimate of the deviance function. This is continued until convergence is reached.<br /><br />Once the fit is obtained, the covariance matrix is obtained as always though in the case of GLM it is weighted by an estimate of the scale chosen when the fit method is called. The default is Pearson's Chi-Squared. The standard errors of the estimate are obtained from the diagonal elements of this matrix.<br /><br />And that's basically it. Much of my time with the GLM code was spent getting my head around the theory and then this algorithm. Then I had to obtain data (<a href="http://artsci.wustl.edu/~jgill/Site/Homepage.html">Jeff Gill</a> was kind enough to give permission for SciPy to use and distribute the datasets from his nice monograph) and write tests to ensure that all of the intermediate and final results were correct for each family. This was no small feat considering there are 6 families after I added the negative binomial family and around 25 possible combinations of families and link functions. Figuring out the correct use of the estimated scale (or dispersion) parameters for each family was particularly challenging. As I mentioned, in the interest of time, I haven't written tests for the noncanonical links for each and every family, but the initial results look good and these tests will come.<br /><br />GLM provided a good base for understanding the remaining code and allowed me to more or less plow through the robust linear models estimator. I had some mathematical difficulties in extending the models to include the correct covariance matrices, since there is no strong theoretical consensus on what they should actually be! More on RLM in my next post, but before then I'll just give a quick look at how the GLM estimator is used.<br /><br />First, the algorithm described above has been made flexible enough to estimate truly binomial data. That is, we can have a vector with each row containing (# of successes, # of failures) as is the case in the star98 data from Dr. Gill, described and included in the models.datasets. It will be useful to have a look at the syntax for this type of data as it's slightly different than the other families.<br /><br /><div = "mycode"><br />In [1]: import models<br /><br />In [2]: import numpy as np<br /><br />In [3]: from models.datasets.star98.data import load<br /><br />In [4]: data = load()<br /><br />In [5]: data.exog = models.functions.add_constant(data.exog)<br /><br />In [6]: trials = data.endog[:,:2].sum(axis=1)<br /><br />In [7]: data.endog[0] # (successes, failures)<br />Out[7]: array([ 452., 355.])<br /><br />In [8]: trials[0] # total trials<br />Out[8]: 807.0 <br /><br />In [9]: from models.glm import GLMtwo as GLM # the name will change<br /><br />In [10]: bin_model = GLM(data.endog, data.exog, family=models.family.Binomial())<br /><br />In [11]: bin_results = bin_model.fit(data_weights = trials)<br /><br />In [12]: bin_results.params<br />Out[12]:<br />array([ -1.68150366e-02, 9.92547661e-03, -1.87242148e-02,<br /> -1.42385609e-02, 2.54487173e-01, 2.40693664e-01,<br /> 8.04086739e-02, -1.95216050e+00, -3.34086475e-01,<br /> -1.69022168e-01, 4.91670212e-03, -3.57996435e-03,<br /> -1.40765648e-02, -4.00499176e-03, -3.90639579e-03,<br /> 9.17143006e-02, 4.89898381e-02, 8.04073890e-03,<br /> 2.22009503e-04, -2.24924861e-03, 2.95887793e+00])<br /><br /></div><br /><br />The main difference between the above and the rest of the families is that you must specify the total number of "trials" (which as you can see is just the sum of success and failures for each observation) as the data_weights argument to fit. This was done so that the current implementation of the Bernouilli family could be extended without rewriting its other derived functions. The code could easily (and might be) extended to calculate these trials so that this argument doesn't need to be specified, but it's sometimes better to be explicit.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-40215589931924568372009-07-11T18:10:00.002-04:002009-08-10T09:25:34.578-04:00GLM Residuals and The Beauty of Stats with Python + SciPyI just finished including the Anscombe residuals for the families in the generalized linear models. The Anscombe residuals for the Binomial family were particularly tricky. It took me a while to work through the math and then figure out the SciPy syntax for what I need (some docs clarification coming...), but if I had had to implement these functions myself (presumably not in Python), it would have taken me more than a week! I thought it provided a good opportunity introduce the residuals and to show off how easy things are in Python with NumPy/SciPy.<br /><br />Note that there is not really a unified terminology for residual analysis in GLMs. I will try to use the most common names for these residuals in introducing the basics and point out any deviations from the norm both here and in the SciPy documentation. <br /><br />In general, residuals should be defined such that their distribution is approximately normal. This is achieved through the use of linear equations, transformed linear equations, and deviance contributions. Below you will find the five most common types of residuals that rely mainly on transformations and one that relies on deviance contributions, though there are as many as 9+ different types of residuals in the literature.<br /><h2>Response Residuals</h2><br />These are simply the observed response values minus the predicted values.<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=r_%7Bi%7D%5E%7BR%7D=y_%7Bi%7D-%5Chat%7B%5Cmu%7D_%7Bi%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?r_%7Bi%7D%5E%7BR%7D=y_%7Bi%7D-%5Chat%7B%5Cmu%7D_%7Bi%7D" title="r_{i}^{R}=y_{i}-\hat{\mu}_{i}" /></a><br /><br />In a classic linear model these are just the expected residuals<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=Y-X%5Cbeta" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?Y-X%5Cbeta" title="Y-X\beta" /></a><br /><br />However, for GLM, these become<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=Y-g%5E%7B-1%7D%28X%5Cbeta%29" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?Y-g%5E%7B-1%7D%28X%5Cbeta%29" title="Y-g^{-1}(X\beta)" /></a><br /><br />where <a href="http://www.codecogs.com/eqnedit.php?latex=g%28%5Ccdot%29" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?g%28%5Ccdot%29" title="g(\cdot)" /></a> is the link function that makes our model's linear predictor comparable to response vector.<br /><br />It is, however, common in GLMs to produce residuals that deviate greatly from the classical assumptions needed in residual analysis, so in addition we have these other residuals which attempt to mitigate deviations from the needed assumptions.<br /><h2>Pearson Residuals</h2><br />The Pearson residuals are a version of the response residuals, scaled by the standard deviation of the prediction. The name comes from the fact that the sum of the Pearson residuals for a Poisson GLM is equal to Pearson's <a href="http://www.codecogs.com/eqnedit.php?latex=%5Cchi%5E%7B2%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?%5Cchi%5E%7B2%7D" title="\chi^{2}" /></a> statistic, a goodness of fit measure.<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=r_%7Bi%7D%5E%7BP%7D=%5Cfrac%7B%28y_%7Bi%7D-%5Chat%7B%5Cmu_%7Bi%7D%7D%29%7D%7B%5Csqrt%7B%5Ctext%7BVAR%7D%5B%5Chat%7B%5Cmu%7D_%7Bi%7D%5D%7D%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?r_%7Bi%7D%5E%7BP%7D=%5Cfrac%7B%28y_%7Bi%7D-%5Chat%7B%5Cmu_%7Bi%7D%7D%29%7D%7B%5Csqrt%7B%5Ctext%7BVAR%7D%5B%5Chat%7B%5Cmu%7D_%7Bi%7D%5D%7D%7D" title="r_{i}^{P}=\frac{(y_{i}-\hat{\mu_{i}})}{\sqrt{\text{VAR}[\hat{\mu}_{i}]}}" /></a><br /><br />The scaling allows plotting of these residuals versus an individual predictor or the outcome to identify outliers. The problem with these residuals though is that asymptotically they are normally distributed, but in practice they can be quite skewed leading to a misinterpretation of the actual dispersion.<br /><h2>Working Residuals</h2><br />These are the difference between the working response and the linear predictor at convergence (ie., the last step in the iterative process).<br /><a href="http://www.codecogs.com/eqnedit.php?latex=r_{i}^{W}=\left(y_{i}-\hat{\mu}_{i}\right)\left.\left(\frac{\partial\eta}{\partial\mu}\right)_{i}\right|_{\mu=\hat{\mu}}" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?r_{i}^{W}=\left(y_{i}-\hat{\mu}_{i}\right)\left.\left(\frac{\partial\eta}{\partial\mu}\right)_{i}\right|_{\mu=\hat{\mu}}" title="r_{i}^{W}=\left(y_{i}-\hat{\mu}_{i}\right)\left.\left(\frac{\partial\eta}{\partial\mu}\right)_{i}\right|_{\mu=\hat{\mu}}" /></a><br /><h2>Anscombe Residuals</h2><br />Anscombe (1960, 1961) describes a general transformation <a href="http://www.codecogs.com/eqnedit.php?latex=A%28y%29" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?A%28y%29" title="A(y)" /></a> in place of <a href="http://www.codecogs.com/eqnedit.php?latex=y" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?y" title="y" /></a> so that they are as close to a normal distribution as possible. (Note: "There is a <strong>maddeningly</strong> great diversity of the forms that the Anscombe residuals take in the literature." I have included the simplest as described below. (Gill 54, emphasis added)) The function <a href="http://www.codecogs.com/eqnedit.php?latex=A%28%5Ccdot%29" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?A%28%5Ccdot%29" title="A(\cdot)" /></a> is given by<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=A%28%5Ccdot%29=%5Cint%5Ctext%7BVAR%7D%5B%5Cmu%5D%5E%7B-%5Cfrac%7B1%7D%7B3%7D%7Dd%5Cmu" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?A%28%5Ccdot%29=%5Cint%5Ctext%7BVAR%7D%5B%5Cmu%5D%5E%7B-%5Cfrac%7B1%7D%7B3%7D%7Dd%5Cmu" title="A(\cdot)=\int\text{VAR}[\mu]^{-\frac{1}{3}}d\mu" /></a><br /><br />This is done for both the response and the predictions. This difference is then scaled by dividing by<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=A%5E%7B%5Cprime%7D%5Cleft%28%5Chat%7B%5Cmu%7D_%7Bi%7D%20%5Cright%20%29%5Csqrt%7B%5Ctext%7BVAR%7D%5Cleft%28%5Chat%7B%5Cmu%7D_i%20%5Cright%20%29%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?A%5E%7B%5Cprime%7D%5Cleft%28%5Chat%7B%5Cmu%7D_%7Bi%7D%20%5Cright%20%29%5Csqrt%7B%5Ctext%7BVAR%7D%5Cleft%28%5Chat%7B%5Cmu%7D_i%20%5Cright%20%29%7D" title="A^{\prime}\left(\hat{\mu}_{i} \right )\sqrt{\text{VAR}\left(\hat{\mu}_i \right )}" /></a><br /><br />so that the Anscombe Residuals are<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=r_{i}^{A}=\frac{A\left(y_{i} \right )-A\left(\hat{\mu}_{i} \right )}{A^{\prime}\left(\hat{\mu}_{i}\right)\sqrt{\text{VAR}\left(\hat{\mu}_{i} \right )} \right )}" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?r_{i}^{A}=\frac{A\left(y_{i} \right )-A\left(\hat{\mu}_{i} \right )}{A^{\prime}\left(\hat{\mu}_{i}\right)\sqrt{\text{VAR}\left(\hat{\mu}_{i} \right )} \right )}" title="r_{i}^{A}=\frac{A\left(y_{i} \right )-A\left(\hat{\mu}_{i} \right )}{A^{\prime}\left(\hat{\mu}_{i}\right)\sqrt{\text{VAR}\left(\hat{\mu}_{i} \right )} \right )}" /></a><br /><br />The Poisson distribution has constant variance <a href="http://www.codecogs.com/eqnedit.php?latex=%5Cmu" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?%5Cmu" title="\mu" /></a> so that it's Anscombe Residuals are simply<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=%5Cfrac%7B3%5Cleft%28y_%7Bi%7D%5E%7B%5Cfrac%7B2%7D%7B3%7D%7D-%5Chat%7B%5Cmu%7D_i%5E%7B%5Cfrac%7B2%7D%7B3%7D%7D%20%5Cright%20%29%7D%7B2%5Chat%7B%5Cmu%7D_%7Bi%7D%5E%7B%5Cfrac%7B1%7D%7B6%7D%7D%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?%5Cfrac%7B3%5Cleft%28y_%7Bi%7D%5E%7B%5Cfrac%7B2%7D%7B3%7D%7D-%5Chat%7B%5Cmu%7D_i%5E%7B%5Cfrac%7B2%7D%7B3%7D%7D%20%5Cright%20%29%7D%7B2%5Chat%7B%5Cmu%7D_%7Bi%7D%5E%7B%5Cfrac%7B1%7D%7B6%7D%7D%7D" title="\frac{3\left(y_{i}^{\frac{2}{3}}-\hat{\mu}_i^{\frac{2}{3}} \right )}{2\hat{\mu}_{i}^{\frac{1}{6}}}" /></a><br /><br />Easy right? Sure was until I ran into the binomial distribution. The Anscombe residuals are built up in a different way for the binomial family. Indeed, the McCullagh and Nelder text does not even provide the general formula for the binomial Anscombe residuals and refers the reader to Anscombe (1953) and Cox & Snell (1968). The problem is that following this transformation for the binomial distribution leads to a rather nasty solution involving the hypergeometric 2F1 function or equivalently the incomplete beta function multiplied by the beta function as shown by Cox and Snell (1968).<br /><br />Gill writes "A partial reason that Anscombe residuals are less popular than other forms is the difficulty in obtaining these tabular values." The tabular values to which he is referring are found in Cox and Snell (1968) p. 260. It is a table of the incomplete beta function that was tabulated numerically for an easy reference. How difficult would it be to get this table with NumPy and SciPy?<br /><br /><div class="mycode"><br />import numpy as np<br />from scipy import special<br />betainc = lambda x: special.betainc(2/3.,2/3.,x)<br /><br />table = np.arange(0,.5,.001).reshape(50,10)<br />results = []<br />for i in table:<br /> results.append(betainc(i))<br /><br />results = np.asarray(results).reshape(50,10)<br /></div><br /><br />That's it!<br /><br />Now the Anscombe residuals for the binomial distribution are<br /><a href="http://www.codecogs.com/eqnedit.php?latex=r_%7Bi%7D%5E%7BA%7D=%5Cfrac%7B%5Cleft%28%5Cphi%5Cleft%28%5Cfrac%7BY_i%7D%7Bn_i%7D%20%5Cright%20%29-%5Cphi%7B%5Cleft%28%20%5Chat%7B%5Cmu%7D_%7Bi%7D%5Cright%20%29%7D%20%5Cright%20%29%7D%7B%5Chat%7B%5Cmu%7D_%7Bi%7D%5E%7B%5Cfrac%7B1%7D%7B6%7D%7D%5Cleft%281-%5Chat%7B%5Cmu%7D_%7Bi%7D%20%5Cright%20%29%5E%7B%5Cfrac%7B1%7D%7B6%7D%7D%7D%5C%20%5C%5C%20%5Ctext%7B%20where%20%7D%20%5Cphi%5Cleft%28%5Cmu%5Cright%29=%5Cint_%7B0%7D%5E%7B%5Cmu%7Dt%5E%7B-%5Cfrac%7B1%7D%7B3%7D%7D%5Cleft%281-t%20%5Cright%20%29%5E%7B-%5Cfrac%7B1%7D%7B3%7D%7D=I_%7B%5Cmu%7D%5Cleft%28%5Cfrac%7B2%7D%7B3%7D,%5Cfrac%7B2%7D%7B3%7D%20%5Cright%20%29B%5Cleft%28%5Cfrac%7B2%7D%7B3%7D,%5Cfrac%7B2%7D%7B3%7D%20%5Cright%20%29" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?r_%7Bi%7D%5E%7BA%7D=%5Cfrac%7B%5Cleft%28%5Cphi%5Cleft%28%5Cfrac%7BY_i%7D%7Bn_i%7D%20%5Cright%20%29-%5Cphi%7B%5Cleft%28%20%5Chat%7B%5Cmu%7D_%7Bi%7D%5Cright%20%29%7D%20%5Cright%20%29%7D%7B%5Chat%7B%5Cmu%7D_%7Bi%7D%5E%7B%5Cfrac%7B1%7D%7B6%7D%7D%5Cleft%281-%5Chat%7B%5Cmu%7D_%7Bi%7D%20%5Cright%20%29%5E%7B%5Cfrac%7B1%7D%7B6%7D%7D%7D%5C%20%5C%5C%20%5Ctext%7B%20where%20%7D%20%5Cphi%5Cleft%28%5Cmu%5Cright%29=%5Cint_%7B0%7D%5E%7B%5Cmu%7Dt%5E%7B-%5Cfrac%7B1%7D%7B3%7D%7D%5Cleft%281-t%20%5Cright%20%29%5E%7B-%5Cfrac%7B1%7D%7B3%7D%7D=I_%7B%5Cmu%7D%5Cleft%28%5Cfrac%7B2%7D%7B3%7D,%5Cfrac%7B2%7D%7B3%7D%20%5Cright%20%29B%5Cleft%28%5Cfrac%7B2%7D%7B3%7D,%5Cfrac%7B2%7D%7B3%7D%20%5Cright%20%29" title="r_{i}^{A}=\frac{\left(\phi\left(\frac{Y_i}{n_i} \right )-\phi{\left( \hat{\mu}_{i}\right )} \right )}{\hat{\mu}_{i}^{\frac{1}{6}}\left(1-\hat{\mu}_{i} \right )^{\frac{1}{6}}}\ \\ \text{ where } \phi\left(\mu\right)=\int_{0}^{\mu}t^{-\frac{1}{3}}\left(1-t \right )^{-\frac{1}{3}}=I_{\mu}\left(\frac{2}{3},\frac{2}{3} \right )B\left(\frac{2}{3},\frac{2}{3} \right )" /></a><br /><br />Where n is the number of trials for each observation (ie., 1, <a href="http://www.codecogs.com/eqnedit.php?latex=n" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?n" title="n" /></a>, or <a href="http://www.codecogs.com/eqnedit.php?latex=n_%7Bi%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?n_%7Bi%7D" title="n_{i}" /></a>)<br /><br />To implement this in the GLM binomial family class, I defined an intermediate function cox_snell similar to the above. It now looks like<br /><br /><div class="mycode"><br />from scipy import special<br />import numpy as np<br /><br />def resid_anscombe(self, Y, mu):<br /> cox_snell = lambda x: special.betainc(2/3.,2/3.,x)*special.beta(2/3.,2/3.)<br /> return np.sqrt(self.n)*(cox_snell(Y)-cox_snell(mu))/(mu**(1/6.)*(1-mu)**(1/6.))<br /></div><br /><br />We multiply the above formula by np.sqrt(self.n) the square root of the number of trials in order to correct for possible use of proportional outcomes Y and mu.<br /><br />Also, see the ticket <a href="http://projects.scipy.org/scipy/ticket/823">here</a> for a note about the incomplete beta function.<br /><h2>Deviance Residuals</h2><br />One other important type of residual in GLMs is the deviance residual. The deviance residual is the most general and also the most useful of the GLM residuals. The IRLS algorithm (as will be shown in a future post) depends on the convergence of the deviance function. The deviance residual then is just the increment to the overall deviance of each observation.<br /><br /><a href="http://www.codecogs.com/eqnedit.php?latex=r_%7Bi%7D%5E%7BD%7D=%5Ctext%7Bsign%7D%5Cleft%28y_i-%5Chat%7B%5Cmu%7D_%7Bi%7D%20%5Cright%20%29%5Csqrt%7B%5Chat%7Bd%7D_i%5E%7B2%7D%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?r_%7Bi%7D%5E%7BD%7D=%5Ctext%7Bsign%7D%5Cleft%28y_i-%5Chat%7B%5Cmu%7D_%7Bi%7D%20%5Cright%20%29%5Csqrt%7B%5Chat%7Bd%7D_i%5E%7B2%7D%7D" title="r_{i}^{D}=\text{sign}\left(y_i-\hat{\mu}_{i} \right )\sqrt{\hat{d}_i^{2}}" /></a><br /><br />where <a href="http://www.codecogs.com/eqnedit.php?latex=%5Chat%7Bd%7D_%7Bi%7D%5E%7B2%7D" target="_blank"><img src="http://codecogs.izyba.com/gif.latex?%5Chat%7Bd%7D_%7Bi%7D%5E%7B2%7D" title="\hat{d}_{i}^{2}" /></a> are defined for each family.<br /><br />Note that most of these residuals also come in variations such as <em>modified</em>, <em>standardized</em>, <em>studentized</em>, and <em>adjusted</em>.<br /><br /><h2>Selected References</h2><br />Anscombe, FJ (1953) "Contribution to the discussion of H. Hotelling's paper."<br /> <em>Journal of the Royal Statistical Society</em>, B, 15, 229-30.<br /><br />Anscombe, FJ (1960) "Rejection of outliers." <em>Technometrics</em>, 2,<br /> 123-47.<br /><br />Anscombe, FJ (1961) "Examination of residuals." <em>Proceedings of<br /> the Fourth Berkeley Symposium on Mathematical Statistics and<br /> Probability.</em> Berkeley: University of California Press.<br /><br />Cox, DR and Snell, EJ (1968). "A generalized definition of<br /> residuals." <em>Journal of the Royal Statistical Society</em> B, 30, 248-65.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-78177896374263849022009-07-11T13:48:00.009-04:002009-07-11T15:28:34.187-04:00Test PostI am testing out my new ASCIIMathML script, so I can type LaTeX in blogger.<br /><br />Note that to correctly view this page requires Internet Explorer 6 + MathPlayer or Firefox/Mozilla/Netscape (with MathML fonts). So the equations won't work in a feed reader.<br /><br />$A(\cdot)=\int\,$VAR$\left[\mu\right]^{-\frac{1}{3}}\, d\mu$<br /><br />Well it appears to work pretty well, except the \text{} tag doesn't work. Hmm, though it should.<br /><br />This works.<br />amath <br />` text(TEXT) `<br />endamath<br /><br /><br />This doesn't work.<br />$\text{TEXT}$<br /><br />You can get the ASCIIMathML script <a href="http://www1.chapman.edu/~jipsen/mathml/asciimath.html">here</a> and follow <a href="http://pleasemakeanote.blogspot.com/2008/09/how-to-post-math-equations-in-blogger.html">these directions</a> to add a javascript widget.<br /><br />Another possible solution for typing equations in Blogger relies on typing LaTeX equations into a remote site. This didn't really appeal to me, as I don't foresee needing any overly complicated mathematical formulae, though it seems like a perfectly workable solution. You can read about it <a href="http://pleasemakeanote.blogspot.com/2008/09/how-to-post-math-equations-in-blogger_06.html">here</a><br />[Edit: I might switch to this approach, since the equations don't display if you don't have the above mentioned requirements or are in a feed reader.]<br /><br />Hat tip to <a href="http://pleasemakeanote.blogspot.com/">Please Make a Note</a> for the pointers.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-34265307225958044432009-07-02T12:04:00.002-04:002009-07-02T12:31:08.499-04:00Generalized Linear ModelsAs I have mentioned, I have spent the last few weeks both in stats books, finding my way around <a href="http://http://www.r-project.org/">R</a>, and cleaning up and refactoring the code for the generalized linear models in the NiPy models code. I have recently hit a wall in this code, so I am trying to clear out some unposted blog drafts. I intended for this post to introduce the generalized linear models approach to estimation; however the full post will have to wait. For now, I will give an introduction to the theory and then explain where I am with the code. <br /><br />Generalized linear models was a topic that was completely foreign to me a few weeks ago, but after a little (okay a lot of) reading the approach seems almost natural. I have found the following references useful:<br /><ul class="disc"><br /><li>Jeff Gill's <em>Generalized Linear Models: A Unified Approach</em>.</li><br /><li>James Hardin and Joseph Hilbe's <em>Generalized Linear Models and Extensions</em>, 2nd edition.</li><br /><li>P. McCullagh and John Nelder's <em>Generalized Linear Models</em>, 2nd edition.</li><br /></ul><br />The basic point of the generalized linear model is to extend the approach taken in classical linear regression to models that have more complex outcomes but ultimately share the linearity property. In this respect, GLM subsumes classical linear regression, probit and logit analysis, loglinear and multinomial response models, and some models that deal with survival data to name a few.<br /><br />In my experience, I have found that econometrics is taught in a compartmentalized manner. This makes sense to a certain extent, as different estimators are tailored to particular problems and data. GLM on the other hand allows the use of a common a technique for obtaining parameter estimates so that it can be studied as a single technique rather than as a collection of distinct approaches.<br /><br />If interested in my ramblings, you can find a draft of my notes as an introduction to GLM <a href="http://eagle1.american.edu/~js2796a/GLMTheoryFirstSection.pdf">here</a>, as Blogger does not support LaTeX... Please note that this is a preliminary and incomplete draft (corrections and clarifications are very welcome). One thing it could definitely use is some clarification by example. However, as I noted, I have run into a bit of a wall trying to extend the binomial family to accept a vector of proportional data, and this is my intended example to walk through the theory and algorithm, so... a subsequent post will have to lay this out once I've got it sorted myself.<br /><br />Generally speaking, there are two basic algorithms for GLM estimation: one is a maximum likelihood optimization based on Newton's method the other is commonly refered to as iteratively (re)weighted least squares (IRLS or IWLS). Our implementation now only covers IRLS. As will be shown, the algorithm itself is pretty simple. It boils down to regressing the transformed (and updated) outcome variable on the untransformed design matrix weighted by the variance of the transformed observations. This is done until we have convergence of the deviance function (twice the log-likelihood ratio of the current and previous estimates). The problem that I am running into with updating the binomial family to accept proportional data (ie., a vector of pairs (successes, total trials) instead of a vector of 1s and 0s for success or failure) is more mathematical than computational. I have either calculated the variance (and therefore the weights) incorrectly, or I am updating the outcome variable incorrectly. Of course, there's always the remote possibility that my data is not well behaved, but I don't think this is the case here.<br /><br />More to come...jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-62013802692437511312009-07-02T01:28:00.006-04:002009-07-02T10:45:23.872-04:00Project StatusI have been making slow but steady progress on the NiPy models code. Right now for the midterm review, we have been focusing on design issues including the user interface and refactoring, test coverage/bug fixing, and some extensions for postestimation statistics. Other than this, I have spent the last month or so with anywhere from ten to fifteen stats, econometrics, or numerical linear algebra and optimization texts open on my desk.<br /><br />The main estimators currently included in the code are generalized least squares, ordinary least squares, weighted least squares, autoregressive AR(p), generalized linear models (with several available distribution families and corresponding link functions), robust linear models, general additive models, and mixed effects models. The test coverage is starting to look pretty good, then there is just squashing the few remaining bugs and improving the postestimation statistics.<br /><br />Some enhancements have also been made to the code. I have started to include some public domain or appropriately copyrighted datasets for testing purposes that could also be useful for examples and tutorials, so that every usage example doesn't have to start with generating your own random data. I have followed pretty closely to the datasets proposal in the <a href="http://www.scipy.org/scipy/scikits/browser/trunk/learn/scikits/learn/datasets">Scikits Learn</a> package. <br /><br />We have also decided to break from the formula framework that is used in NiPy. It was in flux (being changed to take advantage of <a href="http://code.google.com/p/sympy/">SymPy</a> the last I heard) and is intended to be somewhat similar to the formula framework in R. In its place for now, I have written some convenience functions to append a constant to a design matrix or to handle categorical variables for estimation. For the moment, a typical model/estimator is used as<br /><br /><div class="mycode"><br />In [1]: from models.regression import OLS<br /><br />In [2]: from models.datasets.longley.data import load<br /><br />In [3]: from models.functions import add_constant<br /><br />In [4]: data = load()<br /><br />In [5]: data.exog = add_constant(data.exog)<br /><br />In [6]: model = OLS(data.endog, data.exog)<br /><br />In [7]: results = model.fit()<br /><br />In [8]: results.params<br />Out[8]:<br />array([ 1.50618723e+01, -3.58191793e-02, -2.02022980e+00,<br /> -1.03322687e+00, -5.11041057e-02, 1.82915146e+03,<br /> -3.48225863e+06])<br /></div><br /><br />Barring any unforeseen difficulties, the models code <em>should</em> be available as a standalone package shortly after the midterm evaluation rapidly approaching in ten days. The second half of the summer will then be focused on optimizing the code, finalizing design issues, extending the models, and writing good documentation and tutorials so that the code can be included in SciPy!jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-82548484016847364852009-07-01T09:34:00.005-04:002009-07-01T21:03:41.760-04:00Econometrics with Python<blockquote cite="http://uninsubria.academia.edu/documents/0023/6072/Choirat-Seri-JAE-2009.pdf">There is as yet no equivalent of R in applied econometrics. Therefore, the econometric community can still decide to go along the Python path.<br /></blockquote><br /><br />That is Drs. Christine Choirat and Raffello Seri writing in the April issue of the <em>Journal of Applied Econometrics</em>. They have been kind enough to provide me with an <a href="http://uninsubria.academia.edu/documents/0023/6072/Choirat-Seri-JAE-2009.pdf">ungated copy</a> of their review, "Econometrics with Python." Mentioning the, quite frankly, redundant general programming functions and tools that had to be implemented for <a href="http://www.r-project.org/">R</a>, the authors make a nice case for Python as the programming language of choice for applied econometrics. The article provides a quick overview of some of the advantages of using Python and its many built-in libraries, extensions, and tools, gives some speed comparisons, and also mentions a few of the many tools out there in Python community for econometrics including <a href="http://rpy.sourceforge.net/">RPy</a> (RPy2 is now available), and of course NumPy and SciPy. Having spent the last week or more trying to master the basic syntax and usage of R, I very much sympathize with this position. The one complaint I hear most often from my fellow students is that Python is not an industry standard. I hope this can change and is changing, because it's much more of a pleasure to work with Python than the alternatives and that makes for increased productivity.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com8tag:blogger.com,1999:blog-128274497687662608.post-41100906581713717092009-06-15T21:09:00.002-04:002009-06-15T21:30:46.235-04:00Legendre on Least SquaresI found the epigraph to Åke Björk's <span style="font-style:italic;">Numerical Methods for Least Squares Problems</span> to be a nice intersection of my interests.<br /><blockquote cite="http://books.google.com/books?hl=en&lr=&id=PRcOAAAAQAAJ&oi=fnd&pg=PA2&ots=cIayVws04m&sig=Qx53POQcYhjCjwWKwCsP20vRJSc#PPR1,M1"><br />De tous les principes qu'on peut proposer pour cet object, je pense qu'il n'en est pas de plus general, de plus exact, ni d'une application plus facile que celui qui consiste à rendre <span style="font-style:italic;">minimum</span> la somme de carrés des erreurs.</blockquote><br /><blockquote cite="">Of all the principles that can be proposed, I think there is none more general, more exact, or of an easier application than that which consists of rendering the sum of squared errors a minimum.</blockquote><br /> <span style="font-style:italic;">Adrien Marie Legendre, Nouvelles méthodes pour la détermination des orbites des comètes. Appendice. Paris, 1805.</span>jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-37928393992570558022009-06-12T12:33:00.013-04:002009-06-15T21:34:03.980-04:00Design Issues: Understanding Python's superThe current statistical models package is housed in the <a href="http://neuroimaging.scipy.org/site/index.html">NiPy, Neuroimaging in Python,</a> project. Right now, it is designed to rely on Python's built-in super to handle class inheritance. This post will dig a little more into the super function and what it means for the design of the project and future extensions. Note that there are plenty of good places to learn about super and that this post is to help me as much as anyone else. [*Edit: With this in mind, I direct you to <a href="http://www.artima.com/weblogs/viewpost.jsp?thread=236275">Things to Know about Python Super</a> if you really want a deeper and <span style="font-style:italic;">correct</span> understanding of super. This post is mainly a heuristic approach that has helped me in understanding basic usage of super.] You can find the documentation for super <a href="http://docs.python.org/library/functions.html#super">here</a>. If this is a bit confusing, it will, I hope, become clearer after I demonstrate the usage.<br /><br />First, let's take a look at how super actually works for the simple case of single inheritance (right now, we are not planning on using multiple inheritance in the project) and an __init__ chain (note that super can call any of its parent class's methods, but using __init__ is my current use case).<br /><br />The following examples were adapted from some code provided by mentors (thank you!).<br /><br /><div class="mycode"><br />class A(object):<br /> def __init__(self, a):<br /> self.a = a<br /> print 'executing A().__init__'<br /><br />class B(A):<br /> def __init__(self, a):<br /> self.ab = a*2<br /> print 'executing B().__init__'<br /> super(B,self).__init__(a)<br /><br /><br />class C(B):<br /> def __init__(self, a):<br /> self.ac = a*3<br /> print 'executing C().__init__'<br /> super(C,self).__init__(a)<br /></div><br /><br />Now let's have a look at creating an instance of C.<br /><br /><div class="mycode"><br />In [2]: cinst = C(10)<br />executing C().__init__<br />executing B().__init__<br />executing A().__init__<br /><br />In [3]: vars(cinst)<br />Out[3]: {'a': 10, 'ab': 20, 'ac': 30}<br /><br /></div><br /><br />That seems simple enough. Creating an instance of C with a = 10 will also give C the attributes of B(10) and A(10). This means our one instance of C has three attributes: cinst.ac, cinst.ab, cinst.a. The latter two were created by its parent classes (or superclasses) __init__ method. Note that A is also a new-style class. It subclasses the 'object' type.<br /><br />The actual calls to super pass the generic class 'C' and a handle to that class 'self', which is 'cinst' in our case. Super returns the literal parent of the class instance C since we passed it 'self'. It should be noted that A and B were created when we initialized cinst and are, therefore, 'bound' class objects (bound to cinst in memory through the actual instance of class C) and not referring to the class A and class B instructions defined at the interpreter (assuming you are typing along at the interpreter). <br /><br />Okay now let's define a few more classes to look briefly at multiple inheritance.<br /><br /><div class="mycode"><br />class D(A):<br /> def __init__(self, a):<br /> self.ad = a*4<br /> print 'executing D().__init__'<br /> # if super is commented out then __init__ chain ends<br /> #super(D,self).__init__(a)<br /><br /><br />class E(D):<br /> def __init__(self, a):<br /> self.ae = a*5<br /> print 'executing E().__init__'<br /> super(E,self).__init__(a)<br /></div><br /><br />Note that the call to super in D is commented out. This breaks the __init__ chain.<br /><br /><div class="mycode"><br />In [4]: einst = E(10)<br />executing E().__init__<br />executing D().__init__<br /><br />In [5]: vars(einst)<br />Out[5]: {'ad': 40, 'ae': 50}<br /></div><br /><br />If we uncomment the super in D, we get as we would expect<br /><br /><div class="mycode"><br />In [6]: einst = E(10)<br />executing E().__init__<br />executing D().__init__<br />executing A().__init__<br /><br />In [7]: vars(einst)<br />Out[7]: {'a': 10, 'ad': 40, 'ae': 50}<br /></div><br /><br />Ok that's pretty straightforward. In this way super is used to pass off something to its parent class. For instance, say we have a little more realistic example and the instance of C takes some timeseries data that exhibits serial correlation. Then we can have C correct for the covariance structure of the data and "pass it up" to B where B can then perform OLS on our data now that it meets the assumptions of OLS. Further B can pass this data to A and return some descriptive statistics for our data. But remember these are 'bound' class objects, so they're all attributes to our original instance of C. Neat huh? Okay, now let's look at a pretty simple example of multiple inheritance.<br /><br /><div class="mycode"><br />class F(C,E):<br /> def __init__(self, a):<br /> self.af = a*6<br /> print 'executing F().__init__'<br /> super(F,self).__init__(a)<br /></div><br /><br />For this example we are using the class of D, that has super commented out.<br /><br /><div class="mycode"><br />In [8]: finst = F(10)<br />executing F().__init__<br />executing C().__init__<br />executing B().__init__<br />executing E().__init__<br />executing D().__init__<br /><br />In [8]: vars(finst)<br />Out[8]: {'ab': 20, 'ac': 30, 'ad': 40, 'ae': 50, 'af': 60}<br /></div><br /><br />The first time I saw this gave me pause. Why isn't there an finst.a? I was expecting the MRO to go F -> C -> B -> A -> E -> D -> A. Let's take a closer look. The F class has multiple inheritance. It inherits from both C and E. We can see F's method resolution order by doing<br /><br /><div class="mycode"><br />In [9]: F.__mro__<br />Out[9]: <br />(<class '__main__.F'>,<br /> <class '__main__.C'>,<br /> <class '__main__.B'>,<br /> <class '__main__.E'>,<br /> <class '__main__.D'>,<br /> <class '__main__.A'>,<br /> <type 'object'>) <br /></div><br /><br />Okay, so we can see that for F A is a subclass of D but not B. But why?<br /><br /><div class="mycode"><br />In [10]: A.__subclasses__()<br />Out[10]: [<class '__main__.B'>, <class '__main__.D'>]<br /></div> <br /><br />The reason is that A does not have a call to super, so the chain doesn't exist here. When you instantiate F, the hierarchy goes F -> C -> B -> E -> D -> A. The reason that it goes from B -> E is because A does not have a call to super, so it can't pass anything to E (It couldn't pass anything to E because the object.__init__ doesn't take a parameter "a" and because you cannot have a MRO F -> C -> B -> A -> E -> D -> A as this is inconsistent and will give an error!), so A does not cause a problem and the chain ends after D (remember that D's super is commented out, but if it were not then there would be finst.a = 10 as expected). Whew. <br /><br />I'm sure you're thinking "Oh that's (relatively) easy. I'm ready to go crazy with super." But there are a number of things must keep in mind when using super, which makes it necessary for the users of super to proceed carefully.<br /><br />1. super() only works with <a href="http://docs.python.org/glossary.html#term-new-style-class">new-style classes</a>. You can read more about classic/old-style vs new-style classes <a href="http://docs.python.org/reference/datamodel.html#newstyle">here</a>. From there you can click through or just go <a href="http://www.python.org/doc/newstyle/">here for more information on new-style classes</a>. Therefore, you must know that the base classes are new-style. This isn't a problem for our project right now, because I have access to all of the base classes.<br /><br />2. Subclasses must use super if their superclasses do. This is why the user of super must be well-documented. If we have to classes A and B that both use super and a class C that inherits from them, but does not know about super then we will have a problem. Consider the slightly different case<br /><br /><div class="mycode"><br />class A(object):<br /> def __init__(self):<br /> print "executing A().__init__"<br /> super(A, self).__init__()<br /> <br />class B(object):<br /> def __init__(self):<br /> print "executing B().__init__"<br /> super(B, self).__init__()<br /><br />class C(A,B):<br /> def __init__(self):<br /> print "executing C().__init__"<br /> A.__init__(self)<br /> B.__init__(self)<br /># super(C, self).__init__()<br /></div><br /><br />Say class C was defined by someone who couldn't see class A and B, then they wouldn't know about super. Now if we do<br /><br /><div class="mycode"><br />In [11]: C.__mro__<br />Out[11]:<br />(<class '__main__.C'>,<br /> <class '__main__.A'>,<br /> <class '__main__.B'>,<br /> <type 'object'>)<br /><br />In [12]: c = C()<br />executing C().__init__<br />executing A().__init__<br />executing B().__init__<br />executing B().__init__<br /></div><br /><br />B got called twice, but by now this should be expected. A's super calls __init__ on the next object in the MRO which is B (it works this time unlike above because there is no parameter passed with __init__), then C explicitly calls B again.<br /><br />If we uncomment super and comment out the calls to the parent __init__ methods in C then this works as expected.<br /><br />3. Superclasses probably should use super if their subclasses do.<br /><br />We saw this earlier with class D's super call commented out. Note also that A does not have a call to super. The last class in the MRO does not need super *if* there is only one such class at the end.<br /><br />4. Classes must have the exact same call signature.<br /><br />This should be obvious but is important for people to be able to subclass. It is possible however for subclasses to add additional arguments so *args and **kwargs should be probably always be included in the methods that are accessible to subclasses.<br /><br />5. Because of these last three points, the use of super must be explicitly documented, as it has become a part of the interface to our classes.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com0tag:blogger.com,1999:blog-128274497687662608.post-90781270966408870352009-05-09T22:26:00.015-04:002009-05-17T17:15:00.306-04:00Working with Data: Record ArraysA few posts are going to be directed at people who are not as familiar with Python and/or SciPy. I will try to assume as little as possible about what the user knows. Over the summer these types of posts might be put together as a larger tutorial.<br /><br />One of the goals of the SciPy stats project is to provide a transparent way to work with the different arrays types. Towards this end, we are going to work with some data as record arrays (note: right now this is only supported for files containing ASCII characters. If this is not the case, you must do some data cleaning beforehand with Python).<br /><br />We have four data files: a comma-delimited .csv file with variable labels (headers) with all numeric data (<a href="http://eagle1.american.edu/%7Ejs2796a/data/educ_data.csv">here</a>, <a href="http://eagle1.american.edu/%7Ejs2796a/data/CPS_info.txt">info here</a>), a comma-delimited file with headers with a mix of numeric and string variables (<a href="http://eagle1.american.edu/%7Ejs2796a/data/handguns_data.csv">here</a>, <a href="http://eagle1.american.edu/%7Ejs2796a/data/handguns_info.txt">info here</a>), this same file without headers (<a href="http://eagle1.american.edu/%7Ejs2796a/data/handguns_data_noheaders.csv">here</a>), and this same file without headers that is tab delimited (<a href="http://eagle1.american.edu/%7Ejs2796a/data/handguns_data_tab_noheaders.csv">here</a>).<br /><br />First, we are going to put the data into record arrays. Record arrays are simply structured arrays that allow access to the data through attribute access. If you are interested, you can have a look <a href="http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html">here</a> for a bit more about the basic data types in NumPy. Then see <a href="http://www.scipy.org/Cookbook/Recarray">here</a> for more on structured and record arrays (and why access to record arrays is slower) and <a href="http://docs.scipy.org/doc/numpy/user/basics.rec.html">here</a> for more details on creating and working with record arrays.<br /><br />First download the data files and put them in a directory. I am using Linux, so my directory paths will look a little different than those using Windows, but all of the other details should be the same.<br /><br /><div class="mycode"><br />>>>import numpy as np<br />>>>educ_dta=np.recfromcsv('/home/skipper/scipystats/scipystats/data/educ_data.csv')<br /></div><br /><br />We now have our data in a record array. We can take a closer look at what's going on by typing<br /><br /><div class="mycode"><br />>>>educ_dta.dtype<br />dtype([('ahe', '<f8'), ('female', '<i8'), ('ne', '<i8'), ('midwest', '<i8'), ('south', '<i8'), ('west', '<i8'), ('race', '<i8'), ('yrseduc', '<i8'), ('ba', '<i8'), ('hsdipl', '<i8'), ('age', '<i8')])<br />>>>educ_dta.dtype.names<br />('ahe', 'female', 'ne', 'midwest', 'south', 'west', 'race', 'yrseduc', 'ba', 'hsdipl', 'age')<br />>>>educ_dta.ahe<br />array([ 12.5 , 6.25, 17.31, ..., 9.13, 11.11, 14.9 ])<br />>>>educ_dta['ahe']<br />array([ 12.5 , 6.25, 17.31, ..., 9.13, 11.11, 14.9 ])<br /></div><br /><br />We can do the same for each of our other data files. First we have another dataset that contains a string variable for the state name. We proceed just as above, though the string variables will have to be handled differently when used in a statistical model. Then we have the same file, but without any headers information. Last we have a file without headers that is tab-delimited.<br /><br /><div class="mycode"><br />>>>gun_dta=np.recfromcsv('/home/skipper/scipystats/scipystats/data/handguns_data.csv')<br />>>>gun_dta_nh=np.recfromcsv('/home/skipper/scipystats/scipystats/data/handguns_data_noheaders.csv', names=None)<br />>>>gun_dta_tnh=np.recfromtxt('home/skipper/scipystats/scipystats/data/handguns_data_tab_noheaders', names=None, delimiter='\t')<br /></div><br /><br />You can learn more about the possibilities for loading data into arrays <a href="http://docs.scipy.org/numpy/docs/numpy.lib.io.genfromtxt/">here</a>, or by having a look at the doc string for np.genfromtxt in the Python interpreter. All of the functions loadtxt, recfromcsv, recfromtxt, etc. use genfromtxt they just have different default values for the keyword arguments.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com2tag:blogger.com,1999:blog-128274497687662608.post-13920299944846770362009-04-24T09:35:00.005-04:002009-04-26T21:30:12.208-04:00Getting Started with GSoC and SciPyI'm trying to resist making the obligatory "hello world" post here. The best I can do is only mentioning the urge.<br /><br />First, a little bit about myself. My name is Skipper Seabold. I am finishing my first year as a PhD student in economics at American University in Washington, DC, and I have recently been accepted to the <a href="http://code.google.com/soc/">Google Summer of Code 2009</a> to work on the SciPy project. I have been a computer hardware and programming hobbyist since my middle school days. I have built my computers my whole life and back in high school tinkered around with Visual Basic (Apps for AOL 3.0 on Windows 3.x and Windows 95 anyone?), Turbo Pascal, C++, and Java mostly in the context of coursework. Two years ago I was introduced to the <a href="http://www.python.org/">Python programming language</a>, and I haven't looked back. Needless to say I'm very happy to have two of my interests, economics and programming, overlap.<br /><br />This is where SciPy comes in. For those who are unfamiliar with SciPy, I direct you to the homepage <a href="http://www.scipy.org/">here</a>. In short, SciPy is an open source library of algorithms for numerical analysis for those working in engineering or the sciences more broadly defined. The SciPy library depends on NumPy. <a href="http://www.scipy.org/Tentative_NumPy_Tutorial">The Tentative NumPy Tutorial</a> is a good place to start learning about the capabilities of NumPy. And likewise, the <a href="http://www.scipy.org/Getting_Started">Getting Started page</a> has plenty of resources to introduce you to the power of SciPy. In particular the tutorials, documentation, and cookbook are good to look at.<br /><br />What I will be working on this summer is providing a consistent user interface for statistical models and appropriate statistical tests in SciPy similar to those found in other statistics/econometric software packages. I will also provide a unified development framework for those who would like to add to this effort in the future. Updates may be less regular over the next few weeks, but check here for at least weekly updates on the work over the summer.jseaboldhttp://www.blogger.com/profile/16419160807917835179noreply@blogger.com3