# Choosing a numerical programming language for economic research: Julia,

**URL:** <https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697>\
**Category:** Community\
**Tags:** blog, blog-post\
**Created:** [August 13, 2022, 7:37am UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697 "2022-08-13T07:37:53Z")\
**Posts on this page:** 20\
**Page:** 4

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 15, 2022, 4:37am UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/64 "2022-08-15T04:37:46Z")

</div>

> [@samerb](#):
>
> I’m wondering how far away Julia is from those for panel data and causal inference econometrics or whatever other uses draws these people to R/Stata.

I’d almost guarantee there’s a perfectly easy way to do almost all of it already in Julia, just people don’t know what it is (possibly me included). I mostly don’t do regressions with regression formulae unless it’s really simple. The next step up for me from very basic `lm(@formula(a ~ b + c + d),...)` is to immediately go to Turing and specify a full Bayesian model and either sample or optimize it. I really think it’s a mistake to spend a lot of time in that intermediate space with simplified formula syntax, but I know this is a matter of opinion.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 15, 2022, 4:46am UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/65 "2022-08-15T04:46:27Z")

</div>

> [@Albert\_Zevelev](#):
>
> Some data cleaning tasks (generating variables) are much easier in STATA

Sorry I wasn’t specific enough, I was wondering specifically can you give some data cleaning examples that are easy in STATA but you struggle with in Julia.

We could maybe split this out into its own topic, like “Data manipulation in Julia vs other standard tools. Examples” or something.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 15, 2022, 5:35am UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/66 "2022-08-15T05:35:02Z")

</div>

@Jorge_Vieyra, on Slack, also managed to optimize the likelihood calculation. @Jorge_Vieyra , I think it would be great if you posted your code here.

---

<div class="post-metadata">

**Author:** ![JackStrauss](https://avatars.discourse-cdn.com/v4/letter/j/f14d63/32.png) [@JackStrauss](https://discourse.julialang.org/u/JackStrauss)\
**Post date:** [August 15, 2022, 5:53am UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/67 "2022-08-15T05:53:26Z")

</div>

> [@dlakelan](#):
>
> Can you give an example of something that is very easy in STATA that you struggle with in Julia?

Panel models with aggregate fixed effects. i.e. where the fixed effect is spans a number of values (in our case all observations for each variable in 6 hour blocks.) This is not possible either in R’s plm() which otherwise is quite good. One needs STATA.

There is a huge number of very useful, and even more importantly trusted, functions in STATA that exist nowhere else. Will take a long time, and a lot of determination and resources to make Julia statsmodels compete.

And last time I checked, there was noting like R’s plm() in Julia (and google just now did not suggest that had changed)

---

<div class="post-metadata">

**Author:** ![y.lin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/y.lin/32/38569_2.png) [@y.lin](https://discourse.julialang.org/u/y.lin)\
**Post date:** [August 15, 2022, 8:31am UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/68 "2022-08-15T08:31:46Z")

</div>

> [@jlperla](#):
>
> My point is that we should all be careful making unconditional statements ranking languages. They should always be conditional on a particular class of tasks, as you do in your conclusions talking about Julia but not R.

We agree, and we tried to be as careful as possible in our piece. That said, it was not directed at an expert user like yourself who is quite able to evaluate these languages without outside advise.

We surmise that the vast majority of people who use one of these programming languages, once they have learned one, never switch to another language. I think this is the reason for the enduring popularity of Matlab. Students learn it in graduate school and then continue doing so as professors and make their PhD students also learn Matlab. We know of professors who produce world leading research with statistical software that has not been updated in over 30 years.

Most users want a language to solve their problems, they don’t care about languages. And in that case, all you need is one. And then any of the four, Matlab, Python, Julia, and R will usually get the job done.

Those we are writing for, and why we make such a recommendation are those who are deciding between languages for important projects. Research departments in central banks, PhD students starting out, etc. Once they commit to a language, it is very costly to switch to another. And what they care about is very different from what people on the Julia discourse channel tend to care about.

And in such cases, the best language is generally that which can do most things, not a language that is really good in some things and bad in others. This is where R wins out over Julia. And if all you need to statistics, STATA is usually the best.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 15, 2022, 2:11pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/69 "2022-08-15T14:11:43Z")

</div>

> [@y.lin](#):
>
> And in such cases, the best language is generally that which can do most things, not a language that is really good in some things and bad in others. This is where R wins out over Julia

I strongly disagree, but I’m thinking you probably just have a really narrow idea of “most things” for example try coding an agent based model in R, that will be slower than molasses. Or try solving PDEs for each time step in a time series model or write a model with a delay differential equation. Or do an optimization for each week in a 20 year time series… Etc

I worked on a model where after each week of collecting new data we sampled a Bayesian model for a few thousand samples and then ran an optimization the choose a portfolio of bets to maximize an expected outcome at the end of the week. It would have been ridiculously painful to do in R.

---

<div class="post-metadata">

**Author:** ![Jorge\_Vieyra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorge_vieyra/32/6527_2.png) [@Jorge\_Vieyra](https://discourse.julialang.org/u/Jorge_Vieyra)\
**Post date:** [August 15, 2022, 2:37pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/70 "2022-08-15T14:37:00Z")

</div>

Sure, I was thinking of making a github repo and polish all the codes (to be fair) and speed them all up to the best of my ability. I also made a `Project.toml` instead of having all the libraries dangling.  
For the Julia code this is a faster version:

```julia
using StatsBase
using DelimitedFiles
using BenchmarkTools
using LoopVectorization
using Base: Iterators
using IterTools

y = readdlm("data.dat")
y = y[:, 1]
N = length(y)
S = 1e3
o, a, b = [0.001, 0.85, 0.01]
y2 = y .^ 2;
v = var(y);

function likelihood(o, a, b, h, y2, N)
    local lik = 0.0
    for i in 2:N
        h = o + a * y2[i-1] + b * h
        lik += log(h) + y2[i] / h
    end
    return lik
end

function likelihood2(o, a, b, h₀, y2)
    N = length(y2)
    h = similar(y2)
    lik = 0.0
    h[begin] = h₀
    @inbounds for i in 2:N
        h[i] = o + a * y2[i-1] + b * h[i-1]
    end
    @turbo for i in 2:N
        lik += log(h[i]) + y2[i] / h[i]
    end
    return lik
end

lik = likelihood(o, a, b, v, y2, N)
lik2 = likelihood2(o, a, b, v, y2)
lik2 ≈ lik && @info "Original and modified version results match"

bench_o = @benchmark likelihood($o, $a, $b, $v, $y2, $N)
println("output,Julia-inbounds ", VERSION, ",", lik, ",", minimum(bench_o.times))
original_time = bench_o.times |> minimum

bench2 = @benchmark likelihood2($o, $a, $b, $v, $y2)
println("output,Julia-double-loopvec ", VERSION, ",", lik2, ",", minimum(bench2.times))

speedup = original_time / (bench2.times |> minimum) 
@info "Speedup $(round(speedup, digits=1))X"

```

This is the output while running it in my system:

```julia
el_oso@galen:~/Downloads/GARCH$ julia --project garch_modified.jl 
[ Info: Original and modified version results match
output,Julia-inbounds 1.8.0-rc4,206810.2364004015,91227.0
output,Julia-double-loopvec 1.8.0-rc4,206810.23640040183,34652.0
[ Info: Speedup 2.6X

```

---

<div class="post-metadata">

**Author:** ![Jorge\_Vieyra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorge_vieyra/32/6527_2.png) [@Jorge\_Vieyra](https://discourse.julialang.org/u/Jorge_Vieyra)\
**Post date:** [August 15, 2022, 2:49pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/71 "2022-08-15T14:49:42Z")

</div>

This is something that has been criticized at my company. When users want to learn they hit a “wall” because they cannot learn by example for most cases and are confronted with a lot of reference manuals without many examples. Sometimes the reference is the code itself, which is not very useful for a starter.

Also, people want to google and find an already worked out example in StackOverflow and copy/paste it. Finally, people don’t know where to look for documentation and/or help. It is generally available, but not on the usual places, e.g. Discourse instead of StackOverflow, examples in Discourse instead of the package documentation, and available on “chats” like Slack, Zulip, Discord … instead of readily posted and digested for general consumption.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 15, 2022, 3:40pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/72 "2022-08-15T15:40:16Z")

</div>

> [@JackStrauss](#):
>
> Panel models with aggregate fixed effects. i.e. where the fixed effect is spans a number of values (in our case all observations for each variable in 6 hour blocks.) This is not possible either in R’s plm() which otherwise is quite good. One needs STATA.

This is exactly the kind of thing where I think the emphasis on “special syntax” is totally misplaced. Turing will let you write general models in Julia itself and I strongly recommend people with complicated models jump directly to doing that. Even if you don’t want to MCMC sample from it, you can have Turing optimize the model and give you a point estimate.

Most of these comparisons discuss things from the perspective of a beginner. But I honestly think the more important perspective is of an intermediate or advanced user. Of course you have to be able to get there, so beginner stuff is important, but the walls you hit trying to shoehorn everything into R’s lmer syntax or stuff like that are real. Having the full Julia language available to write a model in is what makes Julia so amazing, so if we never consider that particularly perspective we get a very biased view of Julia’s advantages.

---

<div class="post-metadata">

**Author:** ![y.lin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/y.lin/32/38569_2.png) [@y.lin](https://discourse.julialang.org/u/y.lin)\
**Post date:** [August 15, 2022, 5:24pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/73 "2022-08-15T17:24:32Z")

</div>

> [@Jorge\_Vieyra](#):
>
> Sure, I was thinking of making a github repo and polish all the codes (to be fair) and speed them all up to the best of my ability. I also made a `Project.toml` instead of having all the libraries dangling.  
> For the Julia code this is a faster version:

An excellent initiative. Do you mind if we include this on our web appendix?

One observation though, doesn’t @turbo use multithreading? We should have been more clear, but all the other comparisons are single thread. It is not a surprise you get a 2.6 speed up.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 15, 2022, 5:29pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/74 "2022-08-15T17:29:21Z")

</div>

`@turbo` is single-threaded, while `@tturbo` is multi-threaded. The extra `t` stands for threading.

---

<div class="post-metadata">

**Author:** ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)\
**Post date:** [August 15, 2022, 5:43pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/76 "2022-08-15T17:43:21Z")

</div>

Inspired by @Sukera 's comment a bit up the chain about `test` documentation and @johnmyleswhite 's comments on identifying and trying to tackle issues, I opened a PR to the Julia docs that gives an example workflow on how to use `pkg> test` within one’s own package going through how to create tests, run the tests, and add the test suite to one’s own package.

> <https://github.com/JuliaLang/julia/pull/46357>
>
> Hi folks,
> 
> Here is a possible docs PR to the Unit Tests section that explains …an example workflow for creating tests for one's own package. This was inspired by my own difficulties sometimes with creating tests as well as discussion on \[this Discourse post about lack of information on how to create tests\](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/8?u=thecedarprince) . This was a possible idea and it may be better suited for Pkg.jl's docs but I was torn as I thought for discoverability and beginners, this may be a better location.
> 
> Either way, feel free to reject, modify, accept, or request revisions.
> 
> Thanks and have an awesome day!

Feel free to add some review/revisions to it. Found myself on a long bus ride today and figured I would try to contribute another actionable item from this interesting discussion.

---

<div class="post-metadata">

**Author:** ![y.lin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/y.lin/32/38569_2.png) [@y.lin](https://discourse.julialang.org/u/y.lin)\
**Post date:** [August 15, 2022, 5:45pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/77 "2022-08-15T17:45:05Z")

</div>

> [@DNF](#):
>
> `@turbo` is single-threaded, while `@tturbo` is multi-threaded. The extra `t` stands for threading.

Yes, of course. Still it uses SIMD and we could have used that for (some of) the other languages, so doesn’t quite make for a fair comparison. But certainly worth mentioning.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 15, 2022, 5:52pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/78 "2022-08-15T17:52:15Z")

</div>

I disagree. Automatic simd is definitely fair, since it doesn’t use extra resources, just exploits the allocated ones better. In fact, this is not at all limited to `@turbo`, but is something Julia uses automatically in many different places. Would you block simd-vectorization somehow deliberately?

Anyway, numpy and Matlab also surely uses simd internally, probably R too.

Limiting multi-threading makes sense for comparisons, but preventing simd does not make sense.

---

<div class="post-metadata">

**Author:** ![JackStrauss](https://avatars.discourse-cdn.com/v4/letter/j/f14d63/32.png) [@JackStrauss](https://discourse.julialang.org/u/JackStrauss)\
**Post date:** [August 15, 2022, 5:52pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/79 "2022-08-15T17:52:35Z")

</div>

> [@TheCedarPrince](#):
>
> Inspired by @Sukera 's comment a bit up the chain about `test` documentation and @johnmyleswhite 's comments on identifying and trying to tackle issues, I opened a PR to the Julia docs that gives an example workflow on how to use `pkg> test` within one’s own package going through how to create tests, run the tests, and add the test suite to one’s own package.

A fantastic idea.

some time my misunderstanding of Pkg let me to initiate a discussion of here, where I showed myself to be rather ignorant of Julia, alas, I have no shame and here it is, if of use. [Best practice for packages on shared drives](https://discourse.julialang.org/t/best-practice-for-packages-on-shared-drives/83374)

---

<div class="post-metadata">

**Author:** ![TheCedarPrince](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thecedarprince/32/17323_2.png) [@TheCedarPrince](https://discourse.julialang.org/u/TheCedarPrince)\
**Post date:** [August 15, 2022, 5:55pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/80 "2022-08-15T17:55:36Z")

</div>

Thanks and no worries - it’s discussions like these that can lead us to solutions. 😃 I figure even if that isn’t a perfect PR there, it is at least something now we can all critique and make better. Feel free to comment on that PR with additional ideas/thoughts.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [August 15, 2022, 6:06pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/81 "2022-08-15T18:06:15Z")

</div>

I would agree on most cases. However, I believe the authors specifically used an algorithm that can’t be vectorised. The improved Julia code uses a modified algorithm which can be vectorised in part. It’s not really a fair comparison in the spirit of the article.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 15, 2022, 6:12pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/82 "2022-08-15T18:12:37Z")

</div>

Changing the algorithm might not be fair game, I agree, but any automatic simd that might happen should be fine.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [August 15, 2022, 6:14pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/83 "2022-08-15T18:14:43Z")

</div>

Agreed, and @turbo is no less fair than @jit.

---

<div class="post-metadata">

**Author:** ![Jorge\_Vieyra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorge_vieyra/32/6527_2.png) [@Jorge\_Vieyra](https://discourse.julialang.org/u/Jorge_Vieyra)\
**Post date:** [August 15, 2022, 6:21pm UTC](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697/84 "2022-08-15T18:21:02Z")

</div>

In case you are interested I also made a functional version of the code that maybe can be parallelized. It is not faster, but I like it because it has zero allocations and the GC is not called.

```julia
using StatsBase
using DelimitedFiles
using BenchmarkTools
using Base: Iterators
using IterTools

y = readdlm("data.dat")
y = y[:, 1]
N = length(y)
S = 1e3
o, a, b = [0.001, 0.85, 0.01]
y2 = y .^ 2;
v = var(y);

function likelihood_functional(o, a, b, h0, y2)
    N = length(y2)
    function updateh(t)
        newh, i = t
        @inbounds o + a * y2[i] + b * newh, i + 1
    end

    function loginv(t)
        hh, i = t
        @inbounds (i>1) * (log(hh) + y2[i] / hh)
    end

    return mapreduce(loginv, +, Iterators.take(iterated(updateh, (h0, 1, false)), N))
end
lik = likelihood_functional(o, a, b, v, y2)

bench = @benchmark likelihood_functional($o, $a, $b, $v, $y2)
println("output,Julia-functional ", VERSION, ",", lik, ",", minimum(bench.times))

```

It is not really faster, but I like about this is that the timings are more consistent

```julia
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 92.477 μs … 126.451 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 94.428 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 94.960 μs ± 2.134 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

          ▅██▄▇▇▅▂▂▁▂▂ ▂▁▁▂▂▃▁▁▁ ▁ ▂
  ▆▄▃▇▅▄▇▆█████████████▇██████████▇▅▆▇▆▆▆▆▇▆▇▆▇▅▄▅▅▄▄▅▄▄▅▅▄▇██ █
  92.5 μs Histogram: log(frequency) by time 102 μs <

```

[Previous page](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697.md?page=3)

[Next page](https://discourse.julialang.org/t/choosing-a-numerical-programming-language-for-economic-research-julia/85697.md?page=5)