Polyphase decimation filter plot in Octave
I think I might have finally solved it myself.

Afaik, while it is true that the actual allpass filters are 1st order, the trick is to see that the original formula actually just says that the two chains of allpasses run interleaved. That is the reason for the z^-2 and the reason for the extra multiply by z^-1.
Hence it is perfectly reasonable to simply use z^-1 filters in the implementation, and feed the allpasses the data interleaved. 

However, for plotting the response, I had to first compute the cascaded filter coefficients, and then finally apply the z^-1 delay to the A1 chain.

This is the result:


Which is exactly what I wanted.
This is the final code:

c = [0.036681502163648017, 0.13654762463195794, 0.27463175937945444, 0.42313861743656711, 0.56109869787919531, 0.67754004997416184, 0.76974183386322703, 0.83988962484963892, 0.89226081800387902, 0.9315419599631839, 0.96209454837808417, 0.98781637073289585 ];
K = length(c);
i=0; 

b0=[1];
a0=[1];
b1=[1];
a1=[1];

while (i&lt;K)

  beta0 = c(i+1);
  b0=conv(b0,[beta0,0,1]);
  a0=conv(a0,[1,0,beta0]);

  beta1= c(i+2);
  b1=conv(b1,[beta1,0,1]); 
  a1=conv(a1,[1,0,beta1]);

  b1z1 = [0,b1];

  [A0,w] = freqz(b0,a0,512,2); 
  [A1,w] = freqz(b1z1,a1,512,2);

  A=(A0+A1)/2;

  plot(w, 20*log10(abs(A))); % amplitude plot in decibel  
  hold on;
  ylabel("Gain (dB)"), axis([0 1 -120 20])
  xlabel("Normalized frequency");
  grid on

  i+=2;
end


I think this should be the accepted answer, but if anyone has anything to add, please do.