Power spectral density of sigma-delta modulator I am simulating a 3rd order SDM and I want to plot the PSD of the output.

As you may know, the output is a sequence of pulses ranging from -3 to +4.
The number of samples is $N$.

The way I approached the problem is the following:


Use fft on the output sequence 
normalize (divide by $N$)
Square the above signal. 
multiply by 2 ( because I need the one-sided band -- not sure about this ).


The plot should then be $10\log_{10}(signal)$ if I am not mistaken.

Problem is that I get a big deviation from the ideal curve of the expected noise.

Am I doing something wrong?
Everything you have written sounds correct so I assume you forgot to consider the BW for each FFT bin and are comparing to an expected spectral density that is per Hz.  Without windowing or zero padding the BW in each bin will be 1/T where T is the time length of your signal (so in other words 1 bin wide).  Therefore you want to reduce the total power measured in each bin by 10Log(1/T) to get the equivalent power per Hz spectral density.

If you are using windows you can refer to fred harris' classis paper on windowing for the Equivalent Noise BW of all common bins as well as windowing implementation loss for an accurate computation.

https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf