Power spectrum of downsampled image
I have solved the problem.

Since I was taking the sum of 4 pixels in downsampling, what I was doing in effect was convolving with a sum of delta functions (that is, one delta function for each pixel). In Fourier space this translates to multiplying with a cosine, which explains the increasingly strong (relative) discrepancy between input and observed power spectrum.

To mitigate this, Wikipedia helped me out. If we first implement a lowpass filter in Fourier space (filtering out high frequencies, to avoid aliasing later on), and then reduce the sampling rate (to 1 in 4 pixels), we get the desired result, and indeed: the observed power spectrum looks nice up until the cutoff: