The problem is the way you shift the transform of the rectangular window to compute the transform of the Hamming window. You can't just shift on the discrete frequency grid, but you have to evaluate the (continuously) shifted functions on the grid. If $W_r(\omega)$ is the transform of the rectangular window
$$W_r(\omega)=\frac{1}{N}e^{-j\omega (N-1)/2}\frac{\sin(N\omega /2)}{\sin(\omega/2)}\tag{1}$$
then the transform of the Hamming window can be written as
$$W_h(\omega)=\alpha W(\omega)-\frac{\beta}{2}\left[W_r(\omega-\omega_0)+W_r(\omega+\omega_0)\right]$$
with $\alpha=0.54$, $\beta=0.46$, and $\omega_0=2\pi/(N-1)$.
In Matlab/Octave you get:
w = -pi:0.01:pi; N = 20; w0 = 2*pi/(N-1); Wr = 1/N*exp(-j*w*(N-1)/2).*sin(N*w/2)./sin(w/2); W1 = 1/N*exp(-j*(w-w0)*(N-1)/2).*sin(N*(w-w0)/2)./sin((w-w0)/2); W2 = 1/N*exp(-j*(w+w0)*(N-1)/2).*sin(N*(w+w0)/2)./sin((w+w0)/2); Wh = .54*Wr -.23*(W1+W2); plot(w/pi,20*log10(abs(Wr)),w/pi,20*log10(abs(Wh)/max(abs(Wh)))) axis([-1,1,-60,0]), grid
EDIT:
Another error I've mentioned in the comments is the sign of the shifted windows in your formula. The positive sign that you use is valid for a zero-phase Hamming window, i.e. a window that is centered at $n=0$. However, your rectangular window starts at $n=0$, so it's not a zero-phase window. Consequently, combining the two windows as if they were both zero-phase doesn't work. That's where the negative sign in my formula comes from (both windows start at $n=0$). If you change the sign in your code the result is much closer to a Hamming window than before, but it's still off due to the shift on the grid I've mentioned above. In the figure below you can see the difference (blue: exact Hamming; red: your solution with corrected sign):
