ANN: QuantEcon.jl

[Spencer Lyon]

New package QuantEcon.jl.

This package collects code for quantitative economic modeling. It is currently comprised of two main parts:

1. A toolbox of routines useful when doing economics

2. Implementations of types and solution methods for common economic models.

This library has a python twin: QuantEcon.py. The same development team is working on both projects, so we hope to keep the two libraries in sync very closely as new functionality is added.

The library contains all the code necessary to do the computations found on http://quant-econ.net/, a website dedicated to providing lectures that each economics and programming. The website currently (as of 9/18/14) has only a python version, but the Julia version is in late stages of refinement and should be live very soon (hopefully within a week).

The initial version of the website will feature 6 lectures dedicated to helping a new user set up a working Julia environment and learn the basics of the language. In addition to this language specific section, the website will include 22 other lectures on topics including

• statistics: markov processes (continuous and discrete state), auto-regressive processes, the Kalman filter, covariance stationary proceses, ect.
• economic models: the income fluctuation problem, an asset pricing model, the classic optimal growth model, optimal (Ramsey) taxation , the McCall search model
• dynamic programming: shortest path, as well as recursive solutions to economic models

All the lectures have code examples in Julia and most of the 22 will display code from the QuantEcon.jl library.

​

[Gray Calhoun]

That's very exciting! I can't wait to see the new notes.

On Thursday, September 18, 2014 9:14:22 PM UTC-5, Spencer Lyon wrote:

New package QuantEcon.jl.

This package collects code for quantitative economic modeling. It is currently comprised of two main parts:

1. A toolbox of routines useful when doing economics

2. Implementations of types and solution methods for common economic models.

This library has a python twin: QuantEcon.py. The same development team is working on both projects, so we hope to keep the two libraries in sync very closely as new functionality is added.

[...cut...]

[David Anthoff]

This is fantastic!

 

Any experiences/opinions/pitfalls your group discovered on the Python vs Julia question? I guess you pretty much implemented the same algorithms in both languages. Did you find one is faster than the other? Were there major areas where Julia lagged behind Python?

 

Thanks,

David

 

From: julia...@googlegroups.com [mailto:julia...@googlegroups.com] On Behalf Of Spencer Lyon
Sent: Thursday, September 18, 2014 7:14 PM
To: julia...@googlegroups.com
Subject: [julia-users] ANN: QuantEcon.jl

 

New package QuantEcon.jl.

This package collects code for quantitative economic modeling. It is currently comprised of two main parts:

1.      A toolbox of routines useful when doing economics

2.      Implementations of types and solution methods for common economic models.

[Spencer Lyon]

Hi David,

Thanks for the questions, I’ll respond in chunks.

Any experiences/opinions/pitfalls your group discovered on the Python vs Julia question?

I personally love both languages and use both regularly. I find myself reaching for Julia for most things right now — I love the more advanced features like metaprogramming and parallel processing that are either non-existent or not as well integrated into the python language itself (I know there are many packages that provide similar features for python, but they don’t feel as natural as they do in Julia).

Did you find one is faster than the other?

As you might expect, we found that Julia was faster than Python for most things. This is probably an artifact of the type of algorithms we used. Often we would define an operator that loops over a grid, and then iterate on that operator until we find a fixed point. Because of this looping, Julia has a predictable speed advantage.

Were there major areas where Julia lagged behind Python?

I can think of two places where performance in Julia wasn’t as good as performance in python:

1. We often need to do quadrature in the innermost part of our loops (to approximate expected values) and we found that Julia’s quadgk routine was much slower scipy.integrate.fixed_quad. This is not a fair comparison because the algorithms employed by the two functions are very different. Specifically quadgk uses an adaptive algorithm while fixed_quad just uses standard Gaussian quadrature. To get around this we actually implemented a whole suite of quadrature routines in the file quad.jl. After switching the Julia code from using quadgk to our own internal quadrature methods, we got approximately the same solutions, but the code was between 35-115 faster.
2. Linear interpolation. The function numpy.interp does simple piecewise linear interpolation. We ended up using the CoordInterpGrid type from Grid.jl to accomplish this. As of right now the interpolation steps are the biggest bottleneck in most of the functions we have.

We also didn’t really find that Julia was lagging behind Python in terms of libraries that we needed. All the functionality we needed was already available to us. We are using PyPlot.jl to generate graphics, so I guess the some of the Julia code is dependent on Python. Come to think of it the one thing we would like to have that we haven’t been able to find is arbitrary precision linear algebra. In Python this can be achieved through SymPy’s wrapping of mpmath, but as far as I know we don’t yet have arbitrary precision linear algebra in Julia.

[Jameson Nash]

> arbitrary precision linear algebra

i think most linear algebra operations are defined in julia and thus will work for bignums. sticking a call to big() somewhere in the code is usually enough to promote everything

[Spencer Lyon]

Hey Jameson,

Thanks for the tip, but I don’t think we have “most linear algebra operations” defined on BigFloat types. See the example below (on Julia master, one day old):

julia> P1
3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.2  0.5  0.3
 0.0  0.0  1.0

julia> big_P1 = big(P1)
3x3 Array{BigFloat,2}:
 1e+00                                                        …  0e+00
 2.00000000000000011102230246251565404236316680908203125e-01     2.99999999999999988897769753748434595763683319091796875e-01
 0e+00                                                           1e+00

julia> eig(P1)
([0.5,1.0,1.0],
3x3 Array{Float64,2}:
 0.0  0.928477  0.0
 1.0  0.371391  0.514496
 0.0  0.0       0.857493)

julia> eig(big_P1)
ERROR: `eigfact!` has no method matching eigfact!(::Array{BigFloat,2})
 in eigfact at linalg/factorization.jl:451

[Andreas Noack]

...bug you can do

julia> A = big(randn(3,3));

julia> lu(A);
julia> qr(A);
julia> chol(A'A);
julia> big(randn(10,2))\ones(10)
2-element Array{BigFloat,1}:
 -3.159225931277445120604659524255654976434759649748354505550941362505085384003288e-01
  4.105614943861208993931724574474979425594621393874386147118234751387185525595618e-01

I have done some of the work on eigenvalues, but it is not done yet.

Med venlig hilsen

Andreas Noack

[Spencer Lyon]

Ok  that is exciting -- I'm looking forward to seeing some eigenvalue routines on BigFloat Matrices!

<div title="MDH:SGV5IEphbWVzb24sPGRpdj48YnI+PC9kaXY+PGRpdj5UaGFua3MgZm9yIHRoZSB0aXAsIGJ1dCBJ IGRvbid0IHRoaW5rIHdlIGhhdmUgIm1vc3QgbGluZWFyIGFsZ2VicmEgb3BlcmF0aW9ucyIgZGVm aW5lZCBvbiBgQmlnRmxvYXRgIHR5cGVzLiBTZWUgdGhlIGV4YW1wbGUgYmVsb3cgKG9uIEp1bGlh IG1hc3Rlciwgb25lIGRheSBvbGQpOjwvZGl2PjxkaXY+PGJyPjwvZGl2PjxkaXY+YGBgPC9kaXY+ PGRpdj48ZGl2Pmp1bGlhJmd0OyBQMTwvZGl2PjxkaXY+M3gzIEFycmF5e0Zsb2F0NjQsMn06PC9k aXY+PGRpdj4mbmJzcDsxLjAgJm5ic3A7MC4wICZuYnNwOzAuMDwvZGl2PjxkaXY+Jm5ic3A7MC4y ICZuYnNwOzAuNSAmbmJzcDswLjM8L2Rpdj48ZGl2PiZuYnNwOzAuMCAmbmJzcDswLjAgJm5ic3A7 MS4wPC9kaXY+PGRpdj48YnI+PC9kaXY+PGRpdj5qdWxpYSZndDsgYmlnX1AxID0gYmlnKFAxKTwv ZGl2PjxkaXY+M3gzIEFycmF5e0JpZ0Zsb2F0LDJ9OjwvZGl2PjxkaXY+Jm5ic3A7MWUrMDAgJm5i c3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJz cDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNw OyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7 ICZuYnNwOyAmbmJzcDsgJm5ic3A74oCmICZuYnNwOzBlKzAwPC9kaXY+PGRpdj4mbmJzcDsyLjAw MDAwMDAwMDAwMDAwMDExMTAyMjMwMjQ2MjUxNTY1NDA0MjM2MzE2NjgwOTA4MjAzMTI1ZS0wMSAm bmJzcDsgJm5ic3A7IDIuOTk5OTk5OTk5OTk5OTk5ODg4OTc3Njk3NTM3NDg0MzQ1OTU3NjM2ODMz MTkwOTE3OTY4NzVlLTAxPC9kaXY+PGRpdj4mbmJzcDswZSswMCAmbmJzcDsgJm5ic3A7ICZuYnNw OyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7 ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsg Jm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAmbmJzcDsgJm5ic3A7ICZuYnNwOyAm bmJzcDsgJm5ic3A7IDFlKzAwPC9kaXY+PGRpdj48YnI+PC9kaXY+PGRpdj5qdWxpYSZndDsgZWln KFAxKTwvZGl2PjxkaXY+KFswLjUsMS4wLDEuMF0sPC9kaXY+PGRpdj4zeDMgQXJyYXl7RmxvYXQ2 NCwyfTo8L2Rpdj48ZGl2PiZuYnNwOzAuMCAmbmJzcDswLjkyODQ3NyAmbmJzcDswLjA8L2Rpdj48 ZGl2PiZuYnNwOzEuMCAmbmJzcDswLjM3MTM5MSAmbmJzcDswLjUxNDQ5NjwvZGl2PjxkaXY+Jm5i c3A7MC4wICZuYnNwOzAuMCAmbmJzcDsgJm5ic3A7ICZuYnNwOyAwLjg1NzQ5Myk8L2Rpdj48ZGl2 Pjxicj48L2Rpdj48ZGl2Pmp1bGlhJmd0OyBlaWcoYmlnX1AxKTwvZGl2PjxkaXY+RVJST1I6IGBl aWdmYWN0IWAgaGFzIG5vIG1ldGhvZCBtYXRjaGluZyBlaWdmYWN0ISg6OkFycmF5e0JpZ0Zsb2F0 LDJ9KTwvZGl2PjxkaXY+Jm5ic3A7aW4gZWlnZmFjdCBhdCBsaW5hbGcvZmFjdG9yaXphdGlvbi5q bDo0NTE8L2Rpdj48ZGl2PmBgYDwvZGl2Pjxicj5PbiBGcmlkYXksIFNlcHRlbWJlciAxOSwgMjAx NCA4OjI3OjUwIFBNIFVUQy00LCBKYW1lc29uIHdyb3RlOjxibG9ja3F1b3RlIGNsYXNzPSJnbWFp bF9xdW90ZSIgc3R5bGU9Im1hcmdpbjogMDttYXJnaW4tbGVmdDogMC44ZXg7Ym9yZGVyLWxlZnQ6 IDFweCAjY2NjIHNvbGlkO3BhZGRpbmctbGVmdDogMWV4OyI+PGRpdiBkaXI9Imx0ciI+PHNwYW4g c3R5bGU9ImZvbnQtZmFtaWx5OmFyaWFsLHNhbnMtc2VyaWY7Zm9udC1zaXplOjEzcHgiPiZndDsg YXJiaXRyYXJ5IHByZWNpc2lvbiBsaW5lYXIgYWxnZWJyYTwvc3Bhbj48YnI+PGRpdj48c3BhbiBz dHlsZT0iZm9udC1mYW1pbHk6YXJpYWwsc2Fucy1zZXJpZjtmb250LXNpemU6MTNweCI+PGJyPjwv c3Bhbj48L2Rpdj48ZGl2PjxzcGFuIHN0eWxlPSJmb250LWZhbWlseTphcmlhbCxzYW5zLXNlcmlm O2ZvbnQtc2l6ZToxM3B4Ij5pIHRoaW5rIG1vc3QgbGluZWFyIGFsZ2VicmEgb3BlcmF0aW9ucyBh cmUgZGVmaW5lZCBpbiBqdWxpYSBhbmQgdGh1cyB3aWxsIHdvcmsgZm9yIGJpZ251bXMuIHN0aWNr aW5nIGEgY2FsbCB0byBiaWcoKSBzb21ld2hlcmUgaW4gdGhlIGNvZGUgaXMgdXN1YWxseSBlbm91 Z2ggdG8gcHJvbW90ZSBldmVyeXRoaW5nPC9zcGFuPjwvZGl2PjxkaXY+PGJyPjxkaXYgY2xhc3M9 ImdtYWlsX3F1b3RlIj5PbiBGcmksIFNlcCAxOSwgMjAxNCBhdCA4OjA2IFBNLCBTcGVuY2VyIEx5 b24gPHNwYW4gZGlyPSJsdHIiPiZsdDs8YSBocmVmPSJqYXZhc2NyaXB0OiIgdGFyZ2V0PSJfYmxh bmsiIGdkZi1vYmZ1c2NhdGVkLW1haWx0bz0idGEwSm96Z1F2UHNKIiBvbm1vdXNlZG93bj0idGhp cy5ocmVmPSdqYXZhc2NyaXB0Oic7cmV0dXJuIHRydWU7IiBvbmNsaWNrPSJ0aGlzLmhyZWY9J2ph dmFzY3JpcHQ6JztyZXR1cm4gdHJ1ZTsiPnNwZW5jZS4uLkBnbWFpbC5jb208L2E+Jmd0Ozwvc3Bh bj4gd3JvdGU6PGJyPjxibG9ja3F1b3RlIGNsYXNzPSJnbWFpbF9xdW90ZSIgc3R5bGU9Im1hcmdp bjowIDAgMCAuOGV4O2JvcmRlci1sZWZ0OjFweCAjY2NjIHNvbGlkO3BhZGRpbmctbGVmdDoxZXgi PjxkaXYgZGlyPSJsdHIiPjxkaXY+PHAgc3R5bGU9Im1hcmdpbjoxLjJlbSAwcHghaW1wb3J0YW50 Ij5IaSBEYXZpZCw8L3A+CjxwIHN0eWxlPSJtYXJnaW46MS4yZW0gMHB4IWltcG9ydGFudCI+VGhh bmtzIGZvciB0aGUgcXVlc3Rpb25zLCBJ4oCZbGwgcmVzcG9uZCBpbiBjaHVua3MuPC9wPjxzcGFu Pgo8YmxvY2txdW90ZSBzdHlsZT0ibWFyZ2luOjEuMmVtIDBweDtib3JkZXItbGVmdC13aWR0aDo0 cHg7Ym9yZGVyLWxlZnQtc3R5bGU6c29saWQ7Ym9yZGVyLWxlZnQtY29sb3I6cmdiKDIyMSwyMjEs MjIxKTtwYWRkaW5nOjBweCAxZW07Y29sb3I6cmdiKDExOSwxMTksMTE5KTtxdW90ZXM6bm9uZSI+ CjxwIHN0eWxlPSJtYXJnaW46MS4yZW0gMHB4IWltcG9ydGFudCI+QW55IGV4cGVyaWVuY2VzL29w aW5pb25zL3BpdGZhbGxzIHlvdXIgZ3JvdXAgZGlzY292ZXJlZCBvbiB0aGUgUHl0aG9uIHZzIEp1 bGlhIHF1ZXN0aW9uPyA8L3A+CjwvYmxvY2txdW90ZT4KPC9zcGFuPjxwIHN0eWxlPSJtYXJnaW46 MS4yZW0gMHB4IWltcG9ydGFudCI+SSBwZXJzb25hbGx5IGxvdmUgYm90aCBsYW5ndWFnZXMgYW5k IHVzZSBib3RoIHJlZ3VsYXJseS4gSSBmaW5kIG15c2VsZiByZWFjaGluZyBmb3IgSnVsaWEgZm9y IG1vc3QgdGhpbmdzIHJpZ2h0IG5vdyDigJQgSSBsb3ZlIHRoZSBtb3JlIGFkdmFuY2VkIGZlYXR1 cmVzIGxpa2UgbWV0YXByb2dyYW1taW5nIGFuZCBwYXJhbGxlbCBwcm9jZXNzaW5nIHRoYXQgYXJl IGVpdGhlciBub24tZXhpc3RlbnQgb3Igbm90IGFzIHdlbGwgaW50ZWdyYXRlZCBpbnRvIHRoZSBw eXRob24gbGFuZ3VhZ2UgaXRzZWxmIChJIGtub3cgdGhlcmUgYXJlIG1hbnkgcGFja2FnZXMgdGhh dCBwcm92aWRlIHNpbWlsYXIgZmVhdHVyZXMgZm9yIHB5dGhvbiwgYnV0IHRoZXkgZG9u4oCZdCBm ZWVsIGFzIG5hdHVyYWwgYXMgdGhleSBkbyBpbiBKdWxpYSkuPC9wPjxzcGFuPgo8YmxvY2txdW90 ZSBzdHlsZT0ibWFyZ2luOjEuMmVtIDBweDtib3JkZXItbGVmdC13aWR0aDo0cHg7Ym9yZGVyLWxl ZnQtc3R5bGU6c29saWQ7Ym9yZGVyLWxlZnQtY29sb3I6cmdiKDIyMSwyMjEsMjIxKTtwYWRkaW5n OjBweCAxZW07Y29sb3I6cmdiKDExOSwxMTksMTE5KTtxdW90ZXM6bm9uZSI+CjxwIHN0eWxlPSJt YXJnaW46MS4yZW0gMHB4IWltcG9ydGFudCI+RGlkIHlvdSBmaW5kIG9uZSBpcyBmYXN0ZXIgdGhh biB0aGUgb3RoZXI/PC9wPgo8L2Jsb2NrcXVvdGU+Cjwvc3Bhbj48cCBzdHlsZT0ibWFyZ2luOjEu MmVtIDBweCFpbXBvcnRhbnQiPkFzIHlvdSBtaWdodCBleHBlY3QsIHdlIGZvdW5kIHRoYXQgSnVs aWEgd2FzIGZhc3RlciB0aGFuIFB5dGhvbiBmb3IgbW9zdCB0aGluZ3MuIFRoaXMgaXMgcHJvYmFi bHkgYW4gYXJ0aWZhY3Qgb2YgdGhlIHR5cGUgb2YgYWxnb3JpdGhtcyB3ZSB1c2VkLiBPZnRlbiB3 ZSB3b3VsZCBkZWZpbmUgYW4gb3BlcmF0b3IgdGhhdCBsb29wcyBvdmVyIGEgZ3JpZCwgYW5kIHRo ZW4gaXRlcmF0ZSBvbiB0aGF0IG9wZXJhdG9yIHVudGlsIHdlIGZpbmQgYSBmaXhlZCBwb2ludC4g QmVjYXVzZSBvZiB0aGlzIGxvb3BpbmcsIEp1bGlhIGhhcyBhICBwcmVkaWN0YWJsZSBzcGVlZCBh ZHZhbnRhZ2UuPC9wPjxzcGFuPgo8YmxvY2txdW90ZSBzdHlsZT0ibWFyZ2luOjEuMmVtIDBweDti b3JkZXItbGVmdC13aWR0aDo0cHg7Ym9yZGVyLWxlZnQtc3R5bGU6c29saWQ7Ym9yZGVyLWxlZnQt Y29sb3I6cmdiKDIyMSwyMjEsMjIxKTtwYWRkaW5nOjBweCAxZW07Y29sb3I6cmdiKDExOSwxMTks MTE5KTtxdW90ZXM6bm9uZSI+CjxwIHN0eWxlPSJtYXJnaW46MS4yZW0gMHB4IWltcG9ydGFudCI+ V2VyZSB0aGVyZSBtYWpvciBhcmVhcyB3aGVyZSBKdWxpYSBsYWdnZWQgYmVoaW5kIFB5dGhvbj88 L3A+CjwvYmxvY2txdW90ZT4KPC9zcGFuPjxwIHN0eWxlPSJtYXJnaW46MS4yZW0gMHB4IWltcG9y dGFudCI+SSBjYW4gdGhpbmsgb2YgdHdvIHBsYWNlcyB3aGVyZSBwZXJmb3JtYW5jZSBpbiBKdWxp YSB3YXNu4oCZdCBhcyBnb29kIGFzIHBlcmZvcm1hbmNlIGluIHB5dGhvbjo8L3A+CjxvbCBzdHls ZT0ibWFyZ2luOjEuMmVtIDBweDtwYWRkaW5nLWxlZnQ6MmVtIj4KPGxpIHN0eWxlPSJtYXJnaW46 MC41ZW0gMHB4Ij5XZSBvZnRlbiBuZWVkIHRvIGRvIHF1YWRyYXR1cmUgaW4gdGhlIGlubmVybW9z dCBwYXJ0IG9mIG91ciBsb29wcyAodG8gYXBwcm94aW1hdGUgZXhwZWN0ZWQgdmFsdWVzKSBhbmQg d2UgZm91bmQgdGhhdCBKdWxpYeKAmXMgPGNvZGUgc3R5bGU9ImZvbnQtc2l6ZTowLjg1ZW07Zm9u dC1mYW1pbHk6Q29uc29sYXMsSW5jb25zb2xhdGEsQ291cmllcixtb25vc3BhY2U7bWFyZ2luOjBw eCAwLjE1ZW07cGFkZGluZzowcHggMC4zZW07d2hpdGUtc3BhY2U6cHJlLXdyYXA7Ym9yZGVyOjFw eCBzb2xpZCByZ2IoMjM0LDIzNCwyMzQpO2JvcmRlci10b3AtbGVmdC1yYWRpdXM6M3B4O2JvcmRl ci10b3AtcmlnaHQtcmFkaXVzOjNweDtib3JkZXItYm90dG9tLXJpZ2h0LXJhZGl1czozcHg7Ym9y ZGVyLWJvdHRvbS1sZWZ0LXJhZGl1czozcHg7ZGlzcGxheTppbmxpbmU7YmFja2dyb3VuZC1jb2xv cjpyZ2IoMjQ4LDI0OCwyNDgpIj5xdWFkZ2s8L2NvZGU+IHJvdXRpbmUgd2FzIG11Y2ggc2xvd2Vy IDxjb2RlIHN0eWxlPSJmb250LXNpemU6MC44NWVtO2ZvbnQtZmFtaWx5OkNvbnNvbGFzLEluY29u c29sYXRhLENvdXJpZXIsbW9ub3NwYWNlO21hcmdpbjowcHggMC4xNWVtO3BhZGRpbmc6MHB4IDAu M2VtO3doaXRlLXNwYWNlOnByZS13cmFwO2JvcmRlcjoxcHggc29saWQgcmdiKDIzNCwyMzQsMjM0 KTtib3JkZXItdG9wLWxlZnQtcmFkaXVzOjNweDtib3JkZXItdG9wLXJpZ2h0LXJhZGl1czozcHg7 Ym9yZGVyLWJvdHRvbS1yaWdodC1yYWRpdXM6M3B4O2JvcmRlci1ib3R0b20tbGVmdC1yYWRpdXM6 M3B4O2Rpc3BsYXk6aW5saW5lO2JhY2tncm91bmQtY29sb3I6cmdiKDI0OCwyNDgsMjQ4KSI+c2Np cHkuaW50ZWdyYXRlLmZpeGVkX3F1YWQ8L2NvZGU+LiBUaGlzIGlzIDxlbT5ub3Q8L2VtPiBhIGZh aXIgY29tcGFyaXNvbiBiZWNhdXNlIHRoZSBhbGdvcml0aG1zIGVtcGxveWVkIGJ5IHRoZSB0d28g ZnVuY3Rpb25zIGFyZSB2ZXJ5IGRpZmZlcmVudC4gU3BlY2lmaWNhbGx5IDxjb2RlIHN0eWxlPSJm b250LXNpemU6MC44NWVtO2ZvbnQtZmFtaWx5OkNvbnNvbGFzLEluY29uc29sYXRhLENvdXJpZXIs bW9ub3NwYWNlO21hcmdpbjowcHggMC4xNWVtO3BhZGRpbmc6MHB4IDAuM

...

[Tim Holy]

On Friday, September 19, 2014 05:06:41 PM Spencer Lyon wrote:
> 2. Linear interpolation. The function numpy.interp does simple piecewise

> linear interpolation. We ended up using the CoordInterpGrid type from
> Grid.jl to accomplish this. As of right now the interpolation steps are the
> biggest bottleneck in most of the functions we have.

Performance of Grid is not terrific. There's an effort to change that (current
efforts are focused at https://github.com/tlycken/Interpolations.jl). But I
think it's held up now by the need to clarify exactly how to handle
https://github.com/timholy/Grid.jl/issues/37.

Of course, there are quite a few algorithms where, even if optimized,
interpolation is the main bottleneck. But it would be less of a bottleneck if
it were faster.

--Tim

[Chase Coleman]

Just wanted to bump this because the site is now live.

[Florian Oswald]

If I may, I've got an implementation of linear interpolation that I'm quite happy with. It uses a GSL kind of "accelerator" that remembers the index position of the last evaluation. the binary search for the correct bracket of x in the grid is the main bottleneck in my experience. So taking advantage of this (particularly with monotonic functions) speeds up things considerably. 

the documentation is sparse and it's not quite ready to be advertised, but the tests are quite complete, so check it out at 

https://github.com/floswald/ApproXD.jl

[Nils Gudat]

Hi Spencer,

Can I just hop in with a question on your integration/quadrature routine? I'm working on a value function iteration in a fairly large state space (think of the order of 500,000 to 1m grid points per period) and hence have to do calculate a lot of expected values. Naturally, I tried out the QuantEcon package right away after seeing this thread here and I'm a bit puzzled by its performance. When timing the inner part of my loop, in which the quadrature is used (and around 400 integrations are performes), using `quadgk` takes between 3 and 4 seconds to do the job, while using `qnwlege` and `do_quad` clocks in at 18 to 20 seconds. However, for some iterations, the time used with `quadgk` explodes to 350 to 450 seconds.

I checked the examples on the quant-econ website and tried to reproduce the integration of the cos function found here. While I can reproduce the result on the website, the relative advantage of `do_quad` seems to vanish when I narrow the interval, e.g. when integrating from pi/2 to pi, `quadgk` seems to be an order of magnitude faster than `do_quad`, especially when timing both steps for `do_quad` (i.e. including the calculation of nodes and weights).

Do you have a general idea about when `do_quad` outperforms `quadgk` and vice versa? Ideally, I'd like to write my code in a way that chooses the better option in each iteration, but at the moment I'm stumped as for how to figure out which situation suits which routine.

Thanks!

[Spencer Lyon]

Hey Nils,

That is very interesting. I hadn't noticed any cases where `quadgk` is faster. 

I'm not sure I can give a great answer here without seeing your code and really digging into the `quadgk` implementation. I can say, however, that [`do_quad` is a very simple function](https://github.com/QuantEcon/QuantEcon.jl/blob/master/src/quad.jl#L602): it simply evaluates the the function at the nodes and computes the dot product between that result and the weights. You should get consistent timing results from `do_quad` for an interval of any length, given a fixed number of nodes. 

One more comment. The example with integrating over cos used 65 nodes. This is much more than needed for "reasonable" accuracy:

```

julia> nodes, weights = qnwlege(65, -2pi, 2pi);

julia> integral = do_quad(x -> cos(x), nodes, weights)

-2.9819896552041314e-15

julia> nodes, weights = qnwlege(21, -2pi, 2pi);

julia> do_quad(x -> cos(x), nodes, weights)

-1.887379141862766e-14

```

As the computational complexity of the function being integrated goes up, using less nodes will likely be a good tradeoff between accuracy and efficiency. In [this example](https://github.com/QuantEcon/QuantEcon.jl/commit/f89abfed0bb4f478cb131c2055916ea6c527c29c) I changed from using `quadgk` to computing the nodes/weights once and then calling `do_quad` and achieved a speed up of 115x for total time for the `compute_lt_price` function to converge.

If you would be willing to share your code I might be able to give better advice. 

[Hans W Borchers]

The integration function appears to have difficulties when the value of the
integral is (very near to) zero. Perhaps it attempts to keep the error term
smaller than the value and therefore needs many more iterations.

I consider this a kind of bug or oversight in the implementation of quadgk, and
not of the algorithm itself. In R the same algorithm (adaptive Gauss-Kronrod
with 7/15 nodes) is implemented as a function with the same name, but gives
results very quickly:

    R> library(pracma)
    R> system.time(q <- quadgk(cos, -2*pi, 2*pi, tol=eps()))
     user  system elapsed
    0.000   0.000   0.001
    > q
    [1] -5.601367e-16

(Not as fast as Julia's quadgk() when the value is farther away from zero)

[Nils Gudat]

Hi Spencer, hi Hans,

Thanks for both of your replies.

Spencer, you make a good point about the nodes - of course 65 are way too many here, and stuffing that many nodes in a small interval is bound to make your `do_quad` function unnecessarily slow. In fact, in my own implementation, I get reasonable accuracy (this being defined as the distance between quadgk and do_quad <1e-10, or the distance in the resulting saving choice in my decision problem being smaller than 0.0001) with just 9 nodes. Then, the time required to complete one round of the loop goes down to 2.9 seconds, narrowly beating quadgk (and obviously beating it by quite a lot when you take into account those iterations at which quadgk implodes...)

Hans, that's an interesting observation as well, but I'm not sure it explains my issue. I'll try to investigate this further though, as this might potentially be an issue with `quadgk` that would merit being raised on Git.

[Spencer Lyon]

That’s great!

One more tip. If you know that the integration bounds are going to be constant over certain iterations (which is likely given that you are doing value function iteration), it helps performance significantly if you precompute the nodes and weights before the iterations begin, and then simply call do_quad using these same nodes and weights each iteration.

If you do have to have changing bounds, there is a helper function called quadrect(f::Function, n, a, b, kind="lege") that takes a function, number of nodes, bounds (a, b), and an optional argument for the type of nodes and computes nodes/weights for you then callsdo_quad`.

[Nils Gudat]

I do have changing bounds unfortunately, so I'm stuck with computing the nodes and weights each time. Using the helper function now, but this doesn't affect performance - is this just meant to reduce the number of lines of code I have to write?

[Spencer Lyon]

Oh ok. Yeas, the helper function is just that: a helper to reduce repetitive code. It should have no impact on performance. 

[Nils Gudat]

In any case, very cool project, I read through the QuantEcon book a while ago when learning Python to do some data analysis and visualization, but now that I'm back to heavier computational stuff Julia is the language of choice.
Do you have any sort of documentation anywhere that includes all the functions you've written? I could obviously go through your Git repository and check every single file, but I was wondering whether there's something more accessible.

[Spencer Lyon]

Right now we don’t have any documentation specific to QuantEcon.jl. We are waiting for this issue to be closed in the main Julia repo so that we can have a more stable backend for generating documentation.

All the functionality from QuantEcon.py is implemented in QuantEcon.jl. The only differences are cases where we write the code in a more pythonic or Julian style. As such, the documentation for the python library should provide a pretty good sense of what QuantEcon.jl offers.

In addition to that I would say that there are three other sources for seeing how the code works:

1. The solutions notebooks for the exercises on the website
2. The tests for the package. We have about 95% coverage, so you should be able to see how just about every function works by looking at the associated tests
3. The lectures on quant-econ.net

In the time between now and when we get dedicated Julia documentation up feel free to post any questions/issues you have at the new QuantEcon google group or on the QuantEcon.jl issue list.

[Hans W Borchers]

If calculating the nodes and weights is the bottleneck, another suggestion might
be on spot:

Sometimes Clenshaw-Curtis is considered to be faster than Gauss-Legendre as
it computes nodes and weights using the Fast Fourier Transform (FFT).

In my own tests some time ago (not in Julia at that time) Clenshaw-Curtis was
up to 5 times faster for the same number of nodes and comparable accuracy.

[Spencer Lyon]

Florian, this sounds very cool.

I'll keep my eye on it and check it out when you feel it is more ready to go.

Thanks for letting me know about it.