Bring Sims' function to estimate linear expectations models to Python
cbrunos
My name is Bruno Rodrigues and at the moment I'm writing my Masters
Thesis at Strasbourg University which deals with estimating DSGE
models (I'm discussing the Great Moderation). For now, I'm using my
thesis' advisor matlab code, which consists of Sims' functions to
estimate linear expectations models (gensys, qzswitch, qzdiv, etc) and
some custom written code to do Bayesian inference and compute irfs.
I wish to translate all this to Python. I came here to discuss this,
following Skipper Seabold advice (we already had a discussion on Scipy-
dev mailing list: http://mail.scipy.org/pipermail/scipy-dev/2012-March/017222.html)
For now, I only have qzswitch and qzdiv, which I didn't even
translate, I found them on some Dynare-Python archives. These
functions are useless by themselves but are needed in gensys, which is
the function that actually solves the linear expectations models. I am
translating this right now, and think I might finish soon. The thing
is, I am a beginner with Python and am simply translating Matlab code
to Python code, and I am sure the resulting code is neither efficient,
nor Pythonic.
It would be awesome if interested people could help, with the
objective of having these functions then ported to Cython for maximum
speed.
What do you think?
Skipper Seabold
> Hi everybody!
>
> My name is Bruno Rodrigues and at the moment I'm writing my Masters
> Thesis at Strasbourg University which deals with estimating DSGE
> models (I'm discussing the Great Moderation). For now, I'm using my
> thesis' advisor matlab code, which consists of Sims' functions to
> estimate linear expectations models (gensys, qzswitch, qzdiv, etc) and
> some custom written code to do Bayesian inference and compute irfs.
>
Are you doing Bayesian inference on the full DSGE model using the kalman filter?
> I wish to translate all this to Python. I came here to discuss this,
> following Skipper Seabold advice (we already had a discussion on Scipy-
> dev mailing list: http://mail.scipy.org/pipermail/scipy-dev/2012-March/017222.html)
>
> For now, I only have qzswitch and qzdiv, which I didn't even
> translate, I found them on some Dynare-Python archives. These
> functions are useless by themselves but are needed in gensys, which is
> the function that actually solves the linear expectations models. I am
> translating this right now, and think I might finish soon. The thing
> is, I am a beginner with Python and am simply translating Matlab code
> to Python code, and I am sure the resulting code is neither efficient,
> nor Pythonic.
>
> It would be awesome if interested people could help, with the
> objective of having these functions then ported to Cython for maximum
> speed.
>
> What do you think?
I think *some* of this could definitely find a home in statsmodels.
Things that are too DSGE specific might best be left to pymaclab
(parsing of model files, etc.) but things that are pure estimation -
ie., Bayesian kalman filtering, structural VAR, impulse response
functions, maybe even a full gensys solver implementation could find a
home here given that we already have some code of interest. Then
statsmodels might just be a dependency for your code and/or pymaclab.
One thing to be careful about is the license for Chris Sims' code. I
don't recall ever seeing anything explicit in his work. We are only
able to accept open source code that is BSD compatible. Loosely, this
is code not based on GPL code. You might ask Chris if he has released
his code under any particular license.
In any event, you'll find some willing eyes to look over your python
code here (given that people have the time to do so) and give you tips
for making it Pythonic and speeding it up, etc.
You might also be interested in perusing our GSoC ideas page to see if
there's anything there that overlaps with your university work that
you might like to participate on for the google summer of code. Google
pays students a stipend to contribute to open source software projects
over the summer.
https://github.com/statsmodels/statsmodels/wiki/GSoC-Ideas
http://code.google.com/soc/
Skipper
cbrunos
Are you doing Bayesian inference on the full DSGE model using the kalman filter?
Yes, the DSGE model's likelihood is evaluated using the kalman filter. My thesis advisor's
assistant translated Frank Scohrfheide's GAUSS code to matlab to do that. Then, the model
is estimated using a random-walk metropolis hastings algorithm, and again, this is code
translated for Frank Schorfheide's GAUSS code. For now, I am still only translating Sims'
code to solve the linear expectations model.
I think *some* of this could definitely find a home in statsmodels.
Things that are too DSGE specific might best be left to pymaclab
(parsing of model files, etc.) but things that are pure estimation -
ie., Bayesian kalman filtering, structural VAR, impulse response
functions, maybe even a full gensys solver implementation could find a
home here given that we already have some code of interest. Then
statsmodels might just be a dependency for your code and/or pymaclab.
One thing to be careful about is the license for Chris Sims' code. I
don't recall ever seeing anything explicit in his work. We are only
able to accept open source code that is BSD compatible. Loosely, this
is code not based on GPL code. You might ask Chris if he has released
his code under any particular license.
my code by LGPL and accepted by statsmodels? This way, the code can be
added to proprietary software, but modifications would have to be released to
the public. As for Frank Schorfheide's code, AFAIK, it doesn't have any license
neither.
The thing is, in my gensys function I use qzdiv and qzswitch from dynare-python
and these functions are GPL-licensed. This could be problematic as gensys
should then be GPL too.
In any event, you'll find some willing eyes to look over your python
code here (given that people have the time to do so) and give you tips
for making it Pythonic and speeding it up, etc.
You might also be interested in perusing our GSoC ideas page to see if
there's anything there that overlaps with your university work that
you might like to participate on for the google summer of code. Google
pays students a stipend to contribute to open source software projects
over the summer.
cbrunos
here is a first draft of gensys' port to python. The file contains:
- gensys.m: Sims' original code
- qzmanip.py: Contains qzswitch and qzdiv, needed for gensys, originally from python-dynare
- qz.py: Provides Generalized Schur decomposition that works just like in matlab by Sven Schreiber (http://econ.schreiberlin.de/schreibersoftware.html)
- gensys.py: Tentative translation of Sims' gensys.m to Python
- example.py: This file contains the necessary arrays to test gensys.py. The matrices in this file is what you get if you solve the model presented in An & Schorfheide - Bayesian analysis of DSGE models (2006)
The code, however, still doesn't work. I have a bug in the svd decomposition, line 65 in gensys.py: ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)
It seems like the matrix I try to feed to svd is of dimension 0 and I don't know why. I will try to look for the bug, but thought that people might be interested in taking a look even if it's still not completed and maybe help
me find the bug (if they have time of course).
Cheers!
Skipper Seabold
My inclination is not to accept LGPL code. Mainly because I am not a
licensing expert, and it sounds only like it will increase the burden
on us, the package distributors, in the future. Would you be willing
to relicense your code?
> The thing is, in my gensys function I use qzdiv and qzswitch from
> dynare-python
> and these functions are GPL-licensed. This could be problematic as gensys
> should then be GPL too.
>
This could be problematic. They didn't look too hard to code from
scratch given that you know what they're supposed to do.
Skipper Seabold
Ah, cool. I will try to have a look soon-ish. In the meantime, you
might consider getting familiar with git and posting your code to
github. So that it's version-controlled and easy for others to grab
and browse to stay updated.
I'm a bit concerned about the translation aspect, if the intention is
to include this in statsmodels, since the license of Sims' code is not
clear. If not, probably no big deal.
Skipper
cbrunos
My inclination is not to accept LGPL code. Mainly because I am not a
licensing expert, and it sounds only like it will increase the burden
on us, the package distributors, in the future. Would you be willing
to relicense your code?
I believe firmly in Free Software as defined by the FSF, but believe more in availability of
code.
> The thing is, in my gensys function I use qzdiv and qzswitch from
> dynare-python
> and these functions are GPL-licensed. This could be problematic as gensys
> should then be GPL too.
>This could be problematic. They didn't look too hard to code from
scratch given that you know what they're supposed to do.
use their implementation of qzdiv and qzswitch, and license with
a BSD-type license?
Ah, cool. I will try to have a look soon-ish. In the meantime, you
might consider getting familiar with git and posting your code to
github. So that it's version-controlled and easy for others to grab
and browse to stay updated.
Yep, I'm already somewhat familiar with github, I'll try to read the
guidelines to use it efficiently.
I'm a bit concerned about the translation aspect, if the intention is
to include this in statsmodels, since the license of Sims' code is not
clear. If not, probably no big deal.
Sims is probably very busy, but I could try to ask him if he is ok with this.
cbrunos
I'm going to ask them what the license is, and if we can use it.
Skipper Seabold
>> My inclination is not to accept LGPL code. Mainly because I am not a
>> licensing expert, and it sounds only like it will increase the burden
>> on us, the package distributors, in the future. Would you be willing
>> to relicense your code?
>
> statsmodels.
Awesome, thanks.
> I believe firmly in Free Software as defined by the FSF, but believe more in
> availability of
> code.
>
>
>>
>> > The thing is, in my gensys function I use qzdiv and qzswitch from
>> > dynare-python
>> > and these functions are GPL-licensed. This could be problematic as
>> > gensys
>> > should then be GPL too.
>> >
>>
>> This could be problematic. They didn't look too hard to code from
>> scratch given that you know what they're supposed to do.
>
> Should I ask the Dynare team if they're willing to accept that we
> use their implementation of qzdiv and qzswitch, and license with
> a BSD-type license?
>
It might be worth a shot. Though the original function was written by
Chris Sims again as far as I can tell. I wouldn't be surprised if all
his code were in the public domain. Nobel laureates can be generous.
Though honestly, it looks like these are just reordering functions,
which I *think* can be achieved by using the SELECT argument to dgges,
or using my ordqz wrapper. I'm pretty sure qzdiv and qzswitch were
written before MATLAB had an ordqz function.
Pablo W.
Pablo W.
cbrunos
In the meantime, youmight consider getting familiar with git and posting your code to
github. So that it's version-controlled and easy for others to grab
and browse to stay updated.
created a branched called DSGE, added a folder named DSGE with the functions and then did a
git add .
git commit -am 'Added DSGE branch'
git push origin dsge
Now the code appears on my fork. Is that how it works?
cbrunos
thank you for your reply. I will take a look at dolo, I already did this past days, but never really took the time
to analyze the code in detail and see how it works.
Is Sven Schreiber active in developing dynare-python?
What is the status of Bayesian estimation of DSGE models with Python? My thesis advisor's assistant translated
a bunch of code originally written by Frank Schorfheide for this (from GAUSS to Matlab) and I plan on bringing all
this to Python/Cython to speed it up (sampling is really slooow for now).
Cheers!
Pablo W.
Le mercredi 14 mars 2012 18:09:52 UTC+1, cbrunos a écrit :
Hello,
thank you for your reply. I will take a look at dolo, I already did this past days, but never really took the time
to analyze the code in detail and see how it works.
Is Sven Schreiber active in developing dynare-python?
What is the status of Bayesian estimation of DSGE models with Python? My thesis advisor's assistant translated
a bunch of code originally written by Frank Schorfheide for this (from GAUSS to Matlab) and I plan on bringing all
this to Python/Cython to speed it up (sampling is really slooow for now).
cbrunos
[...]
- gensys.py: Tentative translation of Sims' gensys.m to Python
The code, however, still doesn't work. I have a bug in the svd decomposition, line 65 in gensys.py: ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)
It seems like the matrix I try to feed to svd is of dimension 0 and I don't know why.[...]
Ok, so I had some spare time this evening and looked further into the issue. Turns out the matrix I'm giving to svd must be an empty array of shape (0,9) but the svd routine in matlab knows how to handle such a matrix: by returning an empty matrix of size 0 by 1 which scipy.linalg.svd does not (an error is returned). Any ideas on how to solve this?
Cheers!
Skipper Seabold
Can you post a self-contained example of what you mean?
Skipper
cbrunos
Yes for example:
>>> a = numpy.zeros((0,5))
and then
>>> scipy.linalg.svd(a)
returns an error: ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)
while Matlab returns an empty array of size 0x1. I've been busy with my masters thesis and haven't had time to look further into the issue.
Skipper Seabold
Yes, for example:
a=numpy.zeros((0,5))
and then
scipy.linalg.svd(a)
returns an error: ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)
while Matlab returns an empty array of size 0x1. I'm busy at the moment with my masters thesis, and haven't had time to further look into the issue.
Something like this would mimic the matlab behavior though I suspect there's something else going on if you're passing around empty arrays.
def svd_matlab(a):
a = np.asarray(a)
if not a.size:
return np.zeros((0,1))
else:
return linalg.svd(a)
Skipper
On Wednesday, April 4, 2012 5:01:02 PM UTC+2, jseabold wrote:
cbrunos
Something like this would mimic the matlab behavior though I suspect there's something else going on if you're passing around empty arrays.
def svd_matlab(a):
a = np.asarray(a)
if not a.size:
return np.zeros((0,1))
else:
return linalg.svd(a)
Skipper
I did that (had to change it a bit), and it seems to work; at least, I don't get that error message anymore. Now there are other bugs, line 84. I'm going to bed for tonight, I'll check it out soon. Thanks for your help!
Cheers!