Fixed. Forgot about the z^0 case. Again, this is just proof of concept. The parsing can be done better and in a much more robust way. Any other bugs, please let me know.

As I explained in my previous comments, assuming a behaviour is not good enough for me. But even with this assumption, in

Expand[(y + z)^4] 

and

Expand[(y + z)^4] /. y -> 1/z

the order is changing with the replacement. The binomial expansion order is not kept. Mathematica recognizes that the result is of a special form, and does undocumented (or at least I did not find the relevant documentation) reordering for the output after the replacement. Again, this is not a problem for computational use of Mathematica, rather it is a strength that it recognizes so many forms of expressions, but it would be nice to know a way of forcing it not to do the reordering in this case after the replacement.

Yes, of course I can write my own algorithms, which is what I ended up doing in this case and also in several other cases (e.g. I am drawing circular arcs using ParametricPlot). I now use

CoefficientList[(z + 1/z)^4*z^4, z]

because it is well documented, and I can get all the information I need from the output.

Polynomials are always output in the increasing order of powers. If you wish to see the highest power to the left (as we might do when handwriting), one might try the following.

ClearAttributes[Plus, Orderless](*Remove the "Orderless" attribute of Plus*);
(*Apply the orderless*)Plus @@(List @@ (*Change head from Plus to List*)Expand[(z + 1/z)^5])

By removing the builtin "Orderless", we tell the Plus sign that it must keep the user provided order.

I don't want the order changed for display purposes. I want to be sure, that what I get as an output now, it will be the same in future versions of Mathematica. Apparently Mathematica uses some algorithms to order expressions, but these are hidden from the user.

Thanks, but this link did not answer my question. The "OrderedForm" operation there did not put the terms of (z+1/z)^4 in increasing powers of z. However, my main question is: Can I rely on the behaviour I see now to stay like this in future versions of Mathematica? In the documentation I do not see any reference as to how Expand[] orders the output, and in this case it is not the normal order of the binomial expansion. I am writing demonstrations to use by others, so I cannot rely on experimenting with the output. I need a guarantee that it will always be the same. Mathematica is very powerful, but because of this, calculations happen behind the scene, and although the output is nice, I don't see how I can control it.