I finally figured it out!
The problem was that unique(Yvec) does not yield a sorted vector. Once I changed the line lev = unique(Yvec) to lev = levels(Yvec) I got the same result as Econometrics.jl (and Stata).
Aside, I use abs to compute the increasing thresholds, while @nosferican you use exp. For values away from 0, it should be fine I think. I have used abs in other cases to ensure that some quantity remains positive, but is there any concern in using abs vs exp?