The idea is quite simple: whenever we call compute_fixed_point(T, ...) we should update the definition of T so that any temporary arrays can be optionally supplied by the caller. See the new compute_lt_price function in the LucasTree class of quantecon/models/lucastree.py for an example (there we avoid allocating Tf on each iteration).

This test could use a more economically meaningful test.

Hi,
in the study group session today, we found some unexpected behavior in the script given in the solution to Exercise 3 in "A First Look at the Kalman Filter".
I think I identified the problem: kn.current_x_hat has to be "flattened":

    # Wrong
    e1[t] = np.sum((x - kn.current_x_hat)**2)

    # Correct
    e1[t] = np.sum((x - kn.current_x_hat.flatten())**2)

kn.current_x_hat is a column vector (2x1), while x is a row vector (1x2), so x - kn.current_x_hat becomes 2x2, and hence np.sum((x - kn.current_x_hat)**2) gets doubled.

Here's how we got the "unexpected behavior":
http://nbviewer.ipython.org/gist/oyamad/1f9b66fda3679f85f794
(Please see the second half part starting with "Exercise 3".)

One may say that the problem is caused by the implementation design of Kalman, where current_xhat is implemented as a column vector. Mathematically it is legitimate, but in python we do not have to care in many cases as numpy guesses correctly to compute e.g. np.dot(A, x) even if x is a row vector. It might be better to consistently stick to row vectors.

We have talked about writing tests, so I thought we should open up an issue to organize this effort. Below is a todo list of modules that are currently (as of 7/9/14) being imported into the main init.py file for the package. My thoughts for this issue are to organize test writing efforts so that we can track which tests have been written, which are being worked on right now, and who is working on them. This way we can write all the tests, but not duplicate work (duplicating work is already unpleasant, but when it is duplication of writing tests, the pain is magnified 😡)

• asset_pricing
• career
• compute_fp
• discrete_rv
• ecdf
• estspec
• ifp
• jv
• kalman
• lae
• linproc
• lqcontrol
• lss
• lucastree
• mc_tools [Matt]
• odu
• optgrowth
• quadsums
• rank_nullspace
• riccati
• robustlq
• tauchen
• quant-econ/examples/ [Matt]

I suggest that as someone starts working on tests for a particular module that they say so in a comment on this issue so that we don't end up duplicating work.

Any other suggestions for how we should proceed are encouraged.

On my laptop, I got the following output when running the tests this morning:

======================================================================
ERROR: quantecon.tests.test_ricatti.dare_test_tjm_2
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/pablo/.local/opt/anaconda/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/home/pablo/Programming/quantecon/quant-econ/quantecon/tests/test_ricatti.py", line 50, in dare_test_tjm_2
    X = solve_discrete_riccati(A, B, Q, R)
  File "/home/pablo/Programming/quantecon/quant-econ/quantecon/matrix_eqn.py", line 162, in solve_discrete_riccati
    Q_tilde = - Q + dot(C.T, solve(Z, C + gamma * BTA)) + gamma * I
  File "/home/pablo/.local/opt/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 381, in solve
    r = gufunc(a, b, signature=signature, extobj=extobj)
  File "/home/pablo/.local/opt/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 90, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
LinAlgError: Singular matrix

======================================================================
ERROR: quantecon.tests.test_ricatti.dare_test_tjm_3
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/pablo/.local/opt/anaconda/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/home/pablo/Programming/quantecon/quant-econ/quantecon/tests/test_ricatti.py", line 66, in dare_test_tjm_3
    X = solve_discrete_riccati(A, B, Q, R)
  File "/home/pablo/Programming/quantecon/quant-econ/quantecon/matrix_eqn.py", line 162, in solve_discrete_riccati
    Q_tilde = - Q + dot(C.T, solve(Z, C + gamma * BTA)) + gamma * I
  File "/home/pablo/.local/opt/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 381, in solve
    r = gufunc(a, b, signature=signature, extobj=extobj)
  File "/home/pablo/.local/opt/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 90, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
LinAlgError: Singular matrix

The problem disappeared when I upgraded Anaconda, and removed the MKL and accelerate packages (since I don't have a license anymore):

The following packages will be REMOVED:

    accelerate:   1.5.1-np18py27_p0 
    mkl:          11.1-np18py27_p2  

The following packages will be UPDATED:

    numexpr:      2.3.1-np19py27_p0  [mkl] --> 2.3.1-np19py27_0 
    numpy:        1.9.0-py27_p0      [mkl] --> 1.9.0-py27_0     
    scikit-learn: 0.15.2-np19py27_p0 [mkl] --> 0.15.2-np19py27_0
    scipy:        0.14.0-np19py27_p0 [mkl] --> 0.14.0-np19py27_0

Has anybody run the tests with the accelerate/mkl packages ?

Step 1:

1. Build py writer in setup.py that writes a version.py file upon built or install
2. Add version to class as an attribute

Step 2:

1. Write some documentation & lectures about checking versions, and upgrading etc.

Currently these equations live in different modules and some of them are not easily discoverable.

For all of them there are also several possible implementations:

• in scipy (calling linpack)
• in quantecon (i.e. doublej)
• in Slycot, which is based on the last opensource version of SLICOT

Slycot, used to be complicated to build, but this has changed, and it is now pip-installable. It is supposed to be efficient and AFAIK is the only option for Generalized Sylvester equations.
The python-control library has a mateqn.py file (https://github.com/python-control/python-control/blob/master/src/mateqn.py) that is also quite complete and that we could partly reuse or depend on. It seems to do extensive error handling and input checks.

My point is that we should deal with all of these matrix equations for optimal control together in a systematic way. For instance we could

• move them to a common module (like matrix_equations.py)
• provide wrapper functions (solve_lyapunov, solve_generalized_sylvester, solve_riccati, etc...) that can switch between different implementations
• investigate which implementation is the best.

Line 38 is repeated here:

self.assertTrue(abs(sssigx - self.ss.C/(1 - self.ss.A**2)))

This doesn't actually do anything. In the preceeding 3 lines we have statements look similar and produce a bool by taking abs( expression) < 2e-8. The problem is that we have sssigx = array([[ 0.02564102]]) and self.ss.C/(1 - self.ss.A**2) = array([[ 0.51282036]]) -- not even close.

@cc7768, do you know what we were supposed to be doing here?

The readthedocs URL is currently

http://quant-econ.readthedocs.org/en/latest/

Need to get rid of the hyphen or supply a custom one.

We should write some tests that check the accuracy of the returns in the future. They currently check to make sure shapes and some other basic information is in line with what is expected.

This is testing several of the functions, but in the future it might be useful to have more thorough tests. Also only handles a scalar case, would be helpful to have matrix cases.

This is something we do a lot of. At present we use SciPy's solve_discrete_lyapunov. It might be worth comparing this against Tom's doublej routine.

All,

I have a good bit of code for solving, simulating, and calibrating continuous time versions of the Solow and Ramsey/Cass/Koopmans growth models. I use the code primarily to teach grad students how to solve IVPs and BVPs.

Would there be interest in incorporating these models into quant-econ? I would need to refactor the code a bit (and write some tests) before submitting a pull request.

Scipy's fixed_quad routine is (as far as I have been able to tell) equivalent to our legendre routine with 5 points.

When I used our own routines in Julia instead of the quadrature that comes in Julia (quadgk) I got huge speedups. I'm hoping we can see a similar performance boost in python, but we would have to do some benchmarks. One place we could really get speedups is when we integrate over the same region inside a for loop, or nested for loops. In this case we could just compute the integration nodes and weights one time and to the quadrature by simply evaluating the integrand at the nodes and having numpy compute a dot product between that result and the weights. Not recomputing nodes/weights each iteration can only help.

Disclaimer: the quadgk routine in Julia is an adaptive routine, so my little performance comparison wasn't exactly fair. However, for all the tests I threw at the two routines they produced effectively equivalent results for the integral .

I think we should strive for 100% python 3 compatibility. Guido has officially announced (more than once) that the python 2 branch will only be maintained and that all new language features will only be implemented in python 3. Doing a quick 2to3 pass over the library shows that we will have to make very minor changes to the existing code for it to be python 3 compliant (mostly just swapping out print statements for a print function call).

One of the more exciting features forthcoming in the next python release (3.5) is the new syntax @ that will mean matrix multiplication. So, instead of writing np.dot(a, b) or a.dot(b), we can write the much more compact a @ b.

@jstac I believe that the quantecon/tests directory needs and __init__.py script.

After a discussion with @jstac, we decided to move the plotting functions in arma.py and ecdf.py out of the main library.

The ones from arma should go into a file examples/arma_plots.py. We can probably do away with the ecdf plot unless anyone objects.

Currently the jv.py file is having issues calling the scipy.optimize minimization routine SLSQP in python 3. Currently, we call the SLSQP in python 2.7 and COBYLA for python 3. Unforunately, the COBYLA routine creates convergence problems for compute_fixed_point and never comes within 1e-2 in either python 2 or python 3 (I terminated it at ~250 iterations when it was oscillating between 0.03 and 0.011). This causes the test test_jv.py to fail because we expect a higher degree of accuracy. We need to figure out why the call to SLSQP isn't working in python 3 even though it works in python 2.

Here is some example code.

import numpy as np
from quantecon.models.jv import JvWorker
from quantecon.compute_fp import compute_fixed_point


jv = JvWorker()
sze = jv.x_grid.size

v_init = jv.x_grid * .65
sol = compute_fixed_point(jv.bellman_operator, v_init,
                          max_iter=3000,
                          error_tol=1e-5)

jv.bellman_operator(sol)

Since the quad branch is closed, I'm opening a new issue with the comments I made about it, so that we can start a fresh discussion.
There are currently a few issues with the quadrature routines. The two items can have side-effects on future code and should be fixed quickly.

• the gridmake function does not enumerate points in a way that is consistent with default Python conventions. If you compute gridmake(array([0,1]), ([-1,2])) it produces:

array([[ 0, -1],
       [ 1, -1],
       [ 0, -2],
       [ 1, -2],

meaning that the first order varies faster. Now, if you want to represent values on this grid by a 2d array vals such that vals[i,j] contains values, then when you do vals.ravel(), you don't enumerate points in the same order because last index is supposed to vary faster. This one is quite annoying.

• Multidimensional functions don't always return an array with the same number of dimensions. That is not a problem in Matlab, but in Python, if I do quad.qnwnorm([2,2]) I get 2-dimensional arrays and with qnwnorm([2]) I get 1-dimensional vectors. This is problematic, since in generic code, it will force one to always distinguish dimension 1 from higher dimensions. My opinion here is that multidimensional routines, should always return multidimensional objects in a predictible fashion. (maybe the 1-d routines could be exported too)
• (minor) It looks like qnwnorm will fail if one column of the covariance matrix (and the corresponding line) is full of zeros. That was already the case in the compecon toolbox, but it is still a common usecase.
• (minor) in the future, we may want to have more engaging names than qnwnorm and co, don't we ?

A possible way to deal with these issues would be to rename quad.py into ce_quad.py so that the latter can be left untouched and remain as close as possible from the original version (including Fortran order). A quadrature.py could then contain Python compliant versions, possibly with more explicit names.

As for the gridmake replacement, the function cartesian does the required thing. It is also faster. (cf http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays)

Actually all these issues concern only the multidimensional functions.

I am wondering how much of the functionality of the compecon toolkit already exists in the scientific Python eco-system? Particularly for optimization: PyOpt is just one of many packages available.

When I was at ZICE, I took notes on various Python implementations of the tools and techniques that I was taught. I have a repo with all of the lecture notes.

As part of the python 3 compatibility push I have been running through all the solutions notebooks. I noticed that the following notebooks took an extra long time to run:

• ifp_solutions
• odu_solutions

This issue is merely a placeholder for us to document where we have noticed opportunities for performance improvements.

Incorporate the linear quadratic nash solution.

Hello

In the NumPy chapter of quant-econ, Part 2: The Scientific Libraries, you wrote
"(Python allocates 3 contiguous 46 bit pieces of memory, and the existing contents of those memory slots are interpreted as float64 values)".

I presume that you meant "... 3 contiguous 64 bit pieces...".

Rudi Farkas

Right now none of our classes print helpful messages when queried at the repl. Let's remedy this by implementing __str__, __repr__, and __unicode__ methods. The easiest thing would probably be to just implement __unicode__ and have the others just call and return its value:

def __str__(self):
    return unicode(self).encode('utf-8')

def __repr__(self):
    return unicode(self).encode('utf-8')

Here is a list of classes we need to update:

• ARMA (arma.py)
• AssetPrices (asset_pricing.py)
• CareerWorkerProblem (career.py)
• DiscreteRV (discrete_rv.py)
• ECDF (ecdf.py)
• ConsumerProblem (ifp.py)
• JvWorker (jv.py)
• Kalman (kalman.py)
• LAE (lae.py)
• LQ (lqcontrol.py)
• LSS (lss.py)
• SearchProblem (odu.py)
• GrowthModel (optgrowth.py)
• RBLQ (robustlq.py)

At this stage we have no discussion of profiling code to find bottlenecks. Here's the IPython solution:

http://pynash.org/2013/03/06/timing-and-profiling.html

I also hear good things about this visual profiling library:

http://pycallgraph.slowchop.com/en/master/

This is a placeholder issue for discussing the possibility of including code using various scientific computing libraries developing by U.S. national labs. Examples include...

• PETSc
• Trilinos
• TAO
• COOPR
• PATH

...amongst others. Python bindings for the first four of these projects already exist:

• petsc4py
• tao4py
• PyTrilinos
• Pyomo

Most of these packages provide state of the art routines for solving large systems of linear, and non-linear equations, non-linear optimization, etc. These routines are significantly better than those available via SciPy. In case you are wondering, I was exposed to most of these technologies as ZICE this past February.

Another interesting alternative is to make use of the NEOS. NEOS has a Python API which one can use to push an optimization/root finding problem to the server and retrieve the results for post processing.

I am not suggesting that we should make any of these tools dependencies of QuantEcon: there are too many architecture specific build requirements (I am not even sure that all of the above tools can be built on a Windows PC). We should consider, however, showing how they can be leveraged from within Python.

I think this one is the current one:

https://scipy-lectures.github.io/intro/matplotlib/matplotlib.html

I have started incorporating a generic initial value problem (IVP) solver that I developed during my PhD for my reseach and teaching. The solver uses finite-difference methods as implemented in the scipy.ode module to do the integration and then uses parametric B-spline interpolation to interpolate value of the solution between the grid points. My IVP solver is a base class from which classes for solving two-point boundary value problems (2PBVP) using shooting methods inherit.

Before I spend time polishing the code and writing tests I want to make sure that similar functionality does not already exists within quantecon and that no one else is actively developing anything is this area.

Shift a number of files now in the top level of quantecon into a new subpackage called models. To be included:

• asset_pricing.py
• career.py
• ifp.py
• jv.py
• lucastree.py
• odu.py
• optgrowth.py

Note: The corresponding lectures in quant-econ.net might need to be edited
slightly to reflect this change, and the notebooks in the solutions directory will certainly need to be changed so as not to break.

We currently only test it with univariate functions. We should add some multivariate cases in the future.

This could also use a more economically interesting test case and a solution to test against. It currently does mostly checks to make sure the robust and simple rule agree, K2F and F2K agree, etc...

@jstac

Have you considering including other repos as part of QuantEcon using git submodules? I think this would allow authors of submodules to contribute whilst still connecting a DOI to their own repo for citation tracking purposes.

I am far from convinced that this would be the ideal way to structure things. I am more interested in whether or not it has been considered already.

I have reviewed the quadrature branch and I think that code is ready to merge.

There are a couple of things that we should still clean up soon, but that shouldn't prevent the branch from being merged. Here are the residual clean-up items:

• In the docstring for qnwequi we don't know what the equidist_pp argument does.
• We need to fill in a description of the docstring for ce_util.gridmake
• The example notebook has a lot of code and pretty pictures, but not much exposition. It could use some love to make it a better teaching tool.

The following is a message from Daisuke Oyama at Tokyo University. It shows instances where the output of mc_tools/mc_compute_stationary.py produces negative values. (The code also contains a nice sets of tests for that function --- I think we can adopt them.)

**

I reproduced instances of outputs of mc_compute_stationary with
negative elements:

https://github.com/oyamad/test_mc_compute_stationary

It may be inevitable given the limited precision in floating point numbers
(especially in executing np.linalg.solve),
but it seems relevant in some contexts such as the Kandori-Mailath-Rob model
with a large number of players and/or a very small probability of mutations.

Should I add a method to the quantecon.ivp.IVP class that computes a linear approximation of the solution to the ODE around a specified point in state/phase space? Under the hood the method would use routines from scipy.linalg to compute eigenvalues and eigenvectors of the Jacobian and use them to construct the linearized solution. I think that this could be done in a fairly (if not completely) general way, but I will need to do a bit more reading to know for sure.

I suppose the real utility if such a method would be largely pedagogical: many, many papers use linearization around a steady state as the sole technique for analyzing dynamics of a model; accuracy of this approach degrades quickly away from steady state, other numerical solution methods are more accurate globally; with method for linearizing we could then demonstrate exactly how much it matters in practice.

Find some simple examples that we can compute answer by hand to test against.

Some basic information needs to be added to this page:

https://pypi.python.org/pypi/quantecon

As a model, see, for example, this page:

https://pypi.python.org/pypi/gunicorn/

Does Github markdown allow the inclusion of Latex?
The Wiki on Writing Tests doesn't parse the matrices correctly (using $$$) in the simple test example.

Hi,
I (re)started a study group with students on quant-econ.
Today we found the following:

Regarding illustrates_lln.py, lines 44-46, if I am not mistaken,

    sample_mean = np.empty(n)
    for i in range(1,n):
        sample_mean[i] = np.mean(data[:i])

should be

    sample_mean = np.empty(n)
    for i in range(n):  # range(1,n) -> range(n)
        sample_mean[i] = np.mean(data[:i+1])  # data[:i] -> data[:i+1]

http://nbviewer.ipython.org/gist/oyamad/4789a8bcaf5aa867e92d

(Should I submit a pull request myself?)

In the code shown in Simulation 1, the following lines are missing:

import numpy as np
from scipy.stats import expon, norm, poisson

Compare the original code.

Issue with 1 x nosetest failure in python3.4:

The issue is that the anaconda installer for python3.4 doesn't package statsmodels or scikits. There are 141 packages in python 3.4 version of anaconda vs. 197 packages in the python 2.7 version.

Also running conda install statsmodels renders an “Unsatisfiable package specifications” error.

http://stackoverflow.com/questions/23343484/python-3-statsmodels

It does look like the majority of statsmodels will however run in python 3.
http://statsmodels.sourceforge.net/devel/introduction.html#python-3

It also looks like statsmodels is available in the continuum repository of packages for py3.3 – but there must be some conflict with python3.4. [If it’s in python 3.3 I can’t imagine it will be long] I have emailed continuum to see if there is a timeline for inclusion. Then perhaps we can compile a list of dependencies and check them against 3.4 standard inclusions.

In the meantime - should we put a warning box on page http://quant-econ.net/getting_started.html that the majority of the package will work with python 3.4 but from time to time a few errors may crop up due to the evolving package support in python3.4?

Update folder name so it isn't picked up by nosetests

Current: test-script
Suggestion: scripts

To unify syntax. Make change in career.py and update the lecture and solution file appropriately.

There can be multiple stationary distributions (different absorbing states for example) in a discrete markov transition matrix. The current code in mc_tools/mc_compute_stationary.py cannot handle this example. We can solve this using eigenvalues instead and this possibly could solve #18 as well. Creating a separate branch and will update on progress.

import numpy as np
import quantecon as qe

P = np.array([[1, 0], [0, 1]])

qe.mc_tools.mc_compute_stationary(P)

>>> LinAlgError: Singular matrix

I am building a conda environment as a precursor to contributing to the project and need a complete listing of the dependencies for quantecon. I poked about a bit on the wiki but couldn't find anything.

I suspect that dependencies are at least:

• NumPy
• SciPy
• Pandas
• IPython
• SymPy

Am I missing anything else?

This Issue is to collect any Data Requirements we may think when writing tests, that will demand infrastructure for importing a data file etc.

1. , , Brief reason, Brief description of data requirement

Current behavior raises an AttributeError. This is undesirable. I think fix, just requires changing the setter for the attributes as follows...

def f_args(self, new_f_args):
        """Set new values for the parameters."""
        if not isinstance(new_f_args, (tuple, type(None))):
            mesg = "IVP.args must be a tuple, not a {}."
            raise AttributeError(mesg.format(new_f_args.__class__))
        elif new_f_args is not None:
            self.set_f_params(*new_f_args)
            self._f_args = new_f_args
        else:
            pass

to

def f_args(self, new_f_args):
        """Set new values for the parameters."""
        if not isinstance(new_f_args, (tuple, type(None))):
            mesg = "IVP.args must be a tuple, not a {}."
            raise AttributeError(mesg.format(new_f_args.__class__))
        else:
            self.set_f_params(*new_f_args)
            self._f_args = new_f_args

Similar fix for the jac_args setter. PR for this will follow shortly.

The discrete Riccati solve in QuantEcon.py/quantecon/matrix_eqn.py can handle cross product terms in the Bellman equation (terms of the form x'Mu). We can exploit this in a future version so the LQ class can solve such problems.

There are some issues on testing by following the Wiki page Creating a Conda development environment.
I may be missing some important points, but let me report them anyway.

1. Should Step 4 come before Step 3?

2. In order to execute unit testing, I had to do conda install also for

   • pytables
   • numba
   • statsmodels

   It would be helpful if this is mentioned in the Wiki page.

   (Actually, I first tried to do testing by using virtualenv on the machine with python installed via a distribution different from anaconda, but it was quite messy to install some numba related stuff, and I gave up this time.)

3. I first installed quantecon into the site-packages directory (without using -e in pip install /PATH/TO/LOCAL/CLONE/quant-econ), but the file "matlab_quad.mat" in tests/data was not installed. Is it as supposed?

This is an issue to track the progress of updating the quantecon documentation to numpydoc conventions. Submit pull requests (PRs) as documentation is updated and mark them out here as you begin working on them so that we don't duplicate work.

Feel free to work on the docs branch. We can submit PRs from there.

• asset_pricing.py
• career.py
• compute_fp.py
• discrete_rv.py
• ecdf.py
• estspec.py
• ifp.py
• jv.py
• kalman.py
• lae.py
• linproc.py
• lqcontrol.py
• lss.py
• lucastree.py
• mc_tools.py
• odu.py
• optgrowth.py
• quadsums.py
• rank_nullspace.py
• riccati.py
• robustlq.py
• tauchen.py

Set up a proper website for QuantEcon. Something along the lines of

http://www.astropy.org/index.html

One option would be to use Pelican to build the site, or some other static website generator.

Sloan's logo should be prominent. A better URL needs to be chosen and set up.

As a side note, it would be nice to have a logo for QuantEcon too --- in the same way that SciPy has a logo.