Choosing a numerical programming language for economic research: Julia,

mcreel · August 17, 2022, 3:01pm

A nice collection of tutorials is GitHub - PaulSoderlind/JuliaTutorial: Julia Tutorial for Finance and Econometrics Students

pdeffebach · August 17, 2022, 4:57pm

Fair enough. I agree with @jlperla that there is a massive benefit to a centralized econ organization which provides these useful tools with good documentation and tutorials (though SciML benefits from lots of pharma money, and I don’t think there is a demand for niche performance tools in the private sector of economics)

It’s still a nice package that hopefully can help out OP.

jlperla · August 17, 2022, 5:50pm

Not sure about that. But I will do one last list of issues for those who might care.

People primarily driven by research (which might not even use Julia if they change the research direction) probably aren’t the key here.

I think what would make a huge difference would be someone with professional coding experience and a big open-source software grant who can go through and fill in the ecosystem. The SciML crew is doing their best, but it is an enormous undertaking and they have papers to ship. On that note, it seems to me that the backbone of many open-source projects out there are people who are not pure researchers (e.g. researchers from organizations who appreciate this sort of service work, teaching faculty, research software engineers) who apply for big grants to do the work.

If anyone out there loves this kind of stuff, they should apply for grants (I suggest working with the SciML organization in doing so - and if you are relatively junior they could even be a PI). Unless one has a big tech company who requires your platform for their research and operations (e.g. pytorch or jax) somebody needs to pay for professional software development. Search up on the funding of jupyter+binder/jupyterbook, Stan, R, numpy/scipy, etc. sometime to see the resources they have dumped in.

Any problems in the julia ecosystem are about resources (and ensuring there is focus given limited resources), not the language itself.

First, a shout out to all of the people (e.g. @tbeason and @mcreel and many others) who have created amazing tutorials. I think those are all necessary but not sufficient for success.

As I said, my perspective here is that the source of issue is not training, discovery, or (directly) documentation. Instead, it is about inconsistent code quality and uniform code interfaces/testing/benchmarking/maintenance of key packages. Economists may have less training
in general, but even the ones who have plenty of experience (and are even told which package to try!) are having trouble. That certainly includes me.

Of course, people need to use different algorithms in different circumstances and there is no way to avoid teaching them those tradeoffs but with julia there is a signal extraction problem. Is it failing because you chose the wrong package, there is a bug in that package, a bug in an upstream package you didn’t know existed, because you used the package wrong, or because it wasn’t the ideal method for your problem? It is hard to triage that - even for people who know what they are doing.

In something like matlab it just amounts to changing flags in the optimizer and/or supplying gradients to choose between different algorithms. It won’t have bugs and you won’t be calling the interface wrong, so you can focus on trying out different algorithms. No signal extraction problem or triage process required. Similarly, in matlab or python if I want to try various types of interpolate I don’t need to search around between 5 different packages - each with their own missing features and quirks and possible performance regressions.

Not to say that these happen all the time and with every package, of course. But it is a lot more frequent than other environments.

And just to be clear: key packages here could just call out to existing C code like NLopt or Dierkx if it was found to be the most robust. Julia long ago fixed the binary dependency installation woes so there is no reason to fear binary dependencies and plenty to love about them. It means someone is maintaining a separate implementation which probably means more testing! Pure Julia implementations sometimes have benefits (e.g., it can achieve better inlining, work with static arrays, handle generic types better) but the baseline Float64 vector versions can all call out to C or Fortran and a user wouldn’t care. Even auto-differentiation often works fine with external packages because you want standardized custom AD rules (e.g. https://sensitivity.sciml.ai/dev/ calculations don’t require you to use the primal calculation in julia).

What matters, then, is having a consolidated interface to swap algorithms and good coverage/maintenance of essential algorithms. And to deprecate packages from that common interface when they fall into disrepair.

I think the first step is to focus. As far as I am concerned, what everyone is doing with DataFrames/FixedEffectModels/etc. is all awesome. I don’t see any reason to push hard on that because it isn’t Julia’s bread and butter. If all someone is doing is data work where specialized estiamtion pacakges already exist and they know R, they probably should just use R rather than switch to Julia.

If someone knows julia really well or the data stuff is a small part of a project, then I think those packages are all on the right track. Just make sure they don’t break (e.g. CSV woes earlier this year).

For bread-and-butter scientific computing there is no reason that Julia should be more frustrating than matlab/scipy. From my side, the high level list and coordination device for this is mostly there: https://docs.sciml.ai/dev/

The problem is that this is a lot more than just listing out packages and adding in documentation.

Some of the things listed (e.g. Interpolation) are deceiving because the interface doesn’t handle the diversity of the tasks.
- The DataInterpolations package does what it does well, but it doesn’t cover multi-dimensional interpolation. But it isn’t like you can just switch to Interpolations.jl because it doesn’t support irregular grid cubic spline - etc.
- So you need to go through and ensure all of the interface support the variations that are required without loss of structure (e.g. you need to distinguish between a regular and irregular grid) and then ensure all of the backend packages actually work, then benchmark and document the algorithm choice in different circumstances. Hopefully documentation + benchmarking + downstream tests will keep the ecosystem in check without monorepos.

Others might not implement the algorithms that economists really need and additional backend algorithms or features are required
- For example, in https://integrals.sciml.ai/stable/basics/IntegralProblem/ there are great options for adaptive quadrature schemes, but what many economists really need are fast, low-dimensional gaussian quadrature. The reason is that these things are often evaluated in an inner loop and often are implementing expectations of some functional form.
- Somehting like GitHub - QuantEcon/Expectations.jl: Expectation operators for Distributions.jl objects could be deprecated and rolled into the integrals setup, but it would take work.
The https://docs.sciml.ai/dev/modules/Optimization/ wrapper is essential, but it is just getting started. Finding a consistent call interface and solution return type is tricky in that domain.
- As one small examples, doesn’t have good ways to swap the sense for all of the optimizers consistently.
- It doesn’t have the interfaces to deal with constrained optimization especially well - especially if you want to use AD to get a sparse jacobian of the constraints and the objective separately. There is no support for Complementarity conditions - which are essential in a lot of economics applications.
- For people solving “real” problems, they often find the commerical optimizers significantly better. We need to show those examples and test things like accessing with Knitro (and add in software if the Optimizers.jl to MOI interface doesn’t work).
- This needs tonnes of tests and examples to ensure that it works seamlessly - which is a lot of work because every optimizer has different arguments and return types.
- It is necessary to reconcile it with other great ecosystems such as JuliaSmoothOptimizers and it isn’t obvious how to achieve that.
Showing consistent patterns of new things that are available by using Julia is a lot of work.
- For example, when and how can you use symbolic expressions and how does it generate jacobians/gradients/etc.? How to use “LabelledArrays” in each of them (if appropriate).
- How can we get gradients of solutions with respect to parameters, and what does that open up.
Also, let me point out the importance of gradients.
- The biggest benefit that economists can get is to start providing gradients to their algorithms, which lets them use new classes of algoritms. In principle, with things like Symbolics and upcoming projects like Enzyme it will be far easier to do that then anything in matlab or scipy.
- Furthermore, people may realize that if they can differentiate their components of the model with respect to the parameters themselves through smart use of AD rules, then Julia can enable whole new clasess of algorithms and methods.
On that note, cross your fingers that Enzyme.jl works out, and provide as many bug reports and use-cases to that project as you can when it is ready to be released into the wild. Zygote is coming towards its end-of-life and ForwardDiff is limited. Without good AD the best Julai could hope for is to match the stability of python/matlab packages. With good AD and the combination of symbolics+AD it can vastly surpass it.

So my list is basically that of Chris’s (and I am purposefully leaving out things like linear regressions which seem on track for now). Do things like pick up a copy of CompEcon – Paul L. Fackler and ensure that all of the key general-purpose algorithms are covered, etc… I think it starts with getting the downstream packages solid and maintained for all of the things people need, and then the next step is filling in that documentation/unit tests/benchmarking. In that process people will find a lot of bugs.

For sure. But SciML is an amazing organization that is willing to shoulder the burden of consolidating the scientific computing ecosystem here - and there are many complementarities and spillovers from pharma applications. It just isn’t a perfect fit, which is why it won’t happen without our involvement.

Just to repeat my statement at the top: package quality/documentation happened with Matlab because academics paid $$$ to them to write/test/document/maintain packages, and it happened with scipy/jupyter/etc. because someone went to ask grant agencies for a lot of money. My understanding is that these didn’t require huge numbers of people, but rather they required a few talented professional developers and enough time to focus rather than writing papers. If that is of interest to anyone who wants to take on components of the ecosystem and they can dedicate the time, then ask agencies for money (within the SciML umbrella)!

Sorry for the long list here. It is my last because I think it fully summarizes my feelings on where the source of the real issues are and the best approaches to fix them. Listen to other economists opinions as well (especially economists giving a “I tried Julia” report), since I may have a limited perspective and already know julia too well.

Palli · August 17, 2022, 7:10pm

What do you have in mind? At least both R and Julia are 1-based (and column-major), so seem like a good match. I’ve just tried to use RCall.jl but don’t use R regularly. I however also use Python (and through PyCall at the time, would now use PythonCall), and with Python (and e.g. C) 0-based, then that could be a source of bugs (as with OffsetArrays.jl). Python and C/C++ being row-major however, when using with or converting code to/from Julia, wouldn’t be a bug, only a performance pitfall to know about.

I would have thought RCall, with its R REPL mode being pretty sweet, for those knowing and sometimes preferring that language. [Is such also needed in PythonCall?]

I see (non-default, unlike for PythonCall that download Python for you by default):

When R_HOME is set to "*", RCall.jl will automatically install R for you using Conda.

@cjdoris PythonCall is pretty sweet (compared to PyCall), in the sense that it take care of Python package dependencies for you (or indirectly with Conda/Mamba). I suppose RCall could do similar, and even reuse the same infrastructure, to make installation transparent. I don’t know what most R users use for CRAN, Conda (or RStudio Package Manager, which also supports PyPi)?

I’m a big believer in using Julia and other languages, e.g. Python or R (for its libraries, but also Julia from other languages, if people can be persuaded to do so). Possibly we need more wrappers for their code?

@y.lin can you clarify:

The reason is that while Stata offers much better algorithms than any of the four languages, users engage with algorithms written by Stata rather than writing their own.

Which one or more algorithms to you have in mind (just curious what could be improved in e.g. Julia)? I’m a bit surprised that slower or somehow worse elsewhere (in other more established languages, e.g. R). Did you have something else in mind, than slower, something workflow-related and “algorithm” then not right or misleading term?

alfaromartino · August 18, 2022, 2:46am

Thanks for all the advice!!! I’ll definitely explore all the suggestions you made in the next few days, when I’ll retake the computation of the model.

Coming back to the topic here, I think there are two matters:

how to improve packages, which @jperla talks about and I can’t say anything meaningful about it.
how a new user could learn these tools more smoothly. I’ll refer exclusively to economists (rather than for general advice for any type of user), since others could have different needs.

I’ll refer only to 2). I cannot give a solution, but provide an example of a path I think would be worth exploring.

First, let me emphasize some aspects. I imagine anyone who is reading this post is interested in the topic. By this, I mean that all of us will invest enough time until it learns the tools. However, it’s not that type of user that I’m concerned about. Even, I know I’ll spend weeks and months until I learn all the tricks.

Rather, I’m thinking of a new user, who doesn’t necessarily have time or will use Julia/R every day. And even for people who have to learn the tool, it’s about making these people’s lives easier and simultaneously not getting bad habits. Reaching a broad audience is what you eventually need to support Julia, since it creates the necessary scale to attract financing. So, in this sense, it’s also related to what @jlperla was saying.

Second, my perspective is not about defending Julia. As this post’s title, I’m thinking of whether people should use Julia/Python/R.

For the example, I imagine that someone asks me whether he should learn Latex, and how to do so. Given my experience with Latex, this is how I’d tackle the matter.

First, I’d ask him “what will be your use?” If you’re aiming at industry/use it too sporadically, don’t learn Latex. I don’t think you’ll need it, and quite likely everyone will be using Office/LibreOffice. So, learning a new tool is always good, but the time is finite, so I don’t think it’s optimal.

Second, assume you’ll work in academy or any sector that uses a lot of math. Then, learning Latex is a valuable tool. So, the question here is: how to spend your time learning, without wasting too much time on finding an answer? There’s no unique answer, but it’s really valuable to know different people’s experiences.

What would I recommend? I use a hybrid method between Lyx and pure Latex. I use pure Latex for the preamble and some specific functions. Instead, Lyx allows me to use keybord shortcuts. Right now, I press ctrl+o for all the operators I commonly use, ctrl+g for all Greek letters, etc. You can even modify its menu and calling it with keyboard shortcuts. For instance, if I press alt+1, I get the menu “Decor” I created for my version of Lyx:

or you can press ctrl+d and get:

Why I’m recommending this? Because you’ll increase your productivity. For instance, I do all the algebra in Latex. Yesterday, I literally wrote 60 pages of pure algebra.

And if you need pure Latex, e.g. to create an environment, you can do it. Moreover, Lyx allows you to create then a keyboard shortcut for that environment (menu “Remark” in my case).

Point of the example: This way of using Latex worked for me—it’s not the truth, but just an approach that suited me and could work for others. Having notes with different opinions is more valuable than discussions about whether Julia is faster than Python. And I think laying out the motives behind the recommendation is key: it can let you know whether the approach suits your needs. For example, people used to tell me that I should learn pure Latex. Why? just because. And, until I found Lyx, I was finding that was answering my concerns about productivity.

Right now, regarding the choice of Julia/R/Python, I feel like I’m spending more time finding those opinions. You only tend to find discussions oriented to package creators/experienced programmers, which are not key for new users from economics. I think I started using Julia for the wrong reasons, but kept using it because of good reasons (eg, I chose it because of its speed, but stayed because it has a more natural syntax for scientific matters).

jmboehm · August 18, 2022, 9:43am

Maybe I’m stating the obvious, but it’s generally a good idea for new researchers (=PhD students) to ask their advisors about choice of software, as the advisors will have a better understanding of the problems that the student will be facing (compared to, say, people on an internet forum, or the authors of a VoxEU column).

nilshg · August 18, 2022, 10:20am

Although there’s also a good chance you then end up progamming in FORTRAN or GAUSS

jmboehm · August 18, 2022, 11:16am

Yes, that thought crossed my mind as well, but I wouldn’t feel that it’s right to recommend a student to disregard the instructions of their advisors even in that case. I’ve seen advisors (in structural labor, where Fortran is still being used a lot) go over their advisees Fortran codes to an extent that probably would not have happened if they had chosen a more modern language.

Bottom line (and then I’m out), when you decide which language to learn to do economics research, I’d advise speaking to people who understand what you want to do (in the medium term, and in the long term), what your background is, and what your situation is (including the support you can get from advisors, colleagues, friends, documentation, and the internet).

nilshg · August 18, 2022, 11:36am

Actually I always wondered why there weren’t more structural people around in the Julia community, it seems their usescases would be well addressed by a language that allows coding high performance models from scratch (not unlike quant macro)? I’m aware of the pyBLP project, something like that feels like it would have had a more natural home in Julia?

tbeason · August 18, 2022, 12:30pm

Macro, macrofinance, and corporate finance structural models are beginning to pop up in Julia. I think you are right that micro models like demand estimation that you mentioned have not come around as much. I would consider it only a matter of time.

Palli · August 18, 2022, 1:19pm

I get 2.85x speedup on top of your 2.6x for 7.4x over @y.lin’s, and yet 2x more should be possible. I was trying out different things (e.g. log/exp to get rid of division, scroll to right to see), that didn’t pan out, and also just lower precision (did work but uncovered a minor bug in LoopVectorization):

julia> @time ( function likelihood4(o, a, b, h₀, y2)
           h = zeros(Float32, length(y2))  # h = similar(y2)
           lik = Float16(0.0)
           h[begin] = h₀; my_range = eachindex(y2)[begin+1:end]
           @inbounds for i in my_range
               h[i] = o + a * y2[i-1] + b * h[i-1]
           end
           @turbo for i in my_range
               lik += log(convert(Float16, h[i])) + convert(Float16, y2[i]) / convert(Float16, h[i]) # log_h = log(h[i]); lik += log_h + ℯ^(log(y2[i]) - log_h) 
           end
           return lik
       end )
  0.000485 seconds (54 allocations: 9.308 KiB)  # that timing is bogus!

julia> @btime lik2 = likelihood4($o, $a, $b, $v, $y2)
  22.606 μs (1 allocation: 11.88 KiB)
26360.732f0

julia> @btime lik2 = likelihood2($o, $a, $b, $v, $y2)
  64.436 μs (1 allocation: 6.00 KiB)
26360.52849330753

@elrod, For some reason, I get Float32 number back, but should get Float16 back, and do if @tturbo or @turbo aren’t used, so that’s the minor bug (could lead to type-instability elsewhere, not a bug, just inconvenience).

I see you updated the appendix:

Post publication p.s.

After our article came out, Jorge Vieyra suggested here some ways to speed up the Julia calculation. It makes use of the @turbo macro […] While we think this proposal is interesting, it is not a fair comparison with the calculations above because they depend on specialist knowledge […]

Note how @turbo requires rewriting the likelihood function by summing the likelihood separately

I did first:

y2 = rand(Float16, 3000); y = y2; # the rest is the same, I didn't have "data.dat" or sensible data, not sure it matters, or the size)

@time using LoopVectorization
  3.918259 seconds (5.64 M allocations: 372.844 MiB, 36.23% compilation time: 80% of which was recompilation)

To be fair, that time alone swamps the runtime of all the actual likelihood functions, but you should be able to eliminate that fixed small startup cost with a sysimage. Actually compiling the code took also almost 10 sec, for combined 12.44 sec. while compiling the original likelihood function was immediate. That compilation overhead should also be possible to reduce to zero. It’s unclear to me if

Note, for the original @time function likelihood ... is immediate, and LoopVectorization is missing some precompilation statements, at least if the code were in a package code would be precompiled (not fully for packages, a TODO for Julia, but I think already if that package were also in the sysimage).

Is all this “specialist knowledge”, I suppose so for the sysimage, how to reduce latency, and reducing precision. But you only need to do it once in your code (for the h array). For the y2 input array, it’s best it’s already in the reduced precision or you’re just shifting work around. For some code, that not O(n), it might still be worthwhile to make a lower-precision copy (e.g. in the function), for other stuff, it would be best to handle with something like:

readdlm(Float16, "data.dat") # Not (yet), possible, Float64 is a sensible default, and I'm not sure this is possible with some other syntax, likely in some other package, CSV.jl?

“Half” (aka Float16) isn’t defined in C standard (so unfair to use? I don’t think so, just Julia better). It is available for, and at least in GCC extension.

y.lin · August 18, 2022, 2:02pm

Thank you all for the thoughtful discussion. If we can insert a link into the published VoxEU piece with a curated list of computational resources for economists, what would you recommend? If no such list exists, we would be happy to host it on our server if you send us recommendations, but it would be better if there was a GitHub page we all could contribute to.

jacobadenbaum · August 18, 2022, 2:40pm

We are here! I think that a lot of this is inertia. More senior people tend to use fortran still, but there’s a solid contingent of younger researchers who mostly tend to use julia. Almost all of the structural people in my phd cohort (who just went on the market this past year) use julia for instance…

TheCedarPrince · August 19, 2022, 9:26am

That is very fair. I will try to draft some language for that, although, after thinking about it, I think it would be better suited for the documentation in Pkg.jl as a separate PR. Will post back later.

Also, to refresh my memory, if you are using a Project.toml to manage the test environment of some package environment (I’ll call mypkg), do you not have to specify dependent packages in the test Project.toml even if they are used in the mypkg environment per you point in 1? My understanding was that, yes, you can have conditional test dependencies that are only available in a test environment but that the activated test environment could not “pull in” packages from the mypkg environment without explicitly adding them (except for mypkg itself). I may be misunderstanding your first point but I do understand your second.

Thanks!

liuyxpp · August 19, 2022, 9:38am

I am very confused about the test environment, especially after I have read the Pkg.jl’s documentation which says if I remember correctly test environment is experimental. Therefore, the tutorial you have drafted shall help me a lot. I never actually have faith in using test functionality correctly.

I tend to run include("any_test_scripts.jl") in mypkg environment instead of using test command.

nilshg · May 30, 2023, 8:19am

Hopefully relevant enough to necropost: just come across this

https://nittanylion.github.io/Grumps.jl/dev/

which implements a cutting-edge estimator for BLP-style models described in this paper:

so it seems that yes, structural IO people are using Julia after all!