I see what the upsampling has done but I fail to understand how the amplitude distribution is periodic?
The reason they can say this is because the symbols are assumed to be equally likely. Think about an upsampled BPSK stream, you'll have approximately the same amount of $+1$ pulses as $-1$ pulses so if you take a large enough block of samples it is periodic with period equal to the symbol period. In fact, this periodicity can be exploited to get an estimate of the symbol period (check out the paper Blind symbol rate estimation using autocorrelation and zero crossing detection). For enough symbols, the method does quite well for a blind estimator.
The author goes on to say "yet within each of the four time offsets the distribution is stationary, leading to the name cyclostationary" : is it stationary within the time offsets?
Yes, as the author says: "within each of the four time offsets the distribution is stationary..." Cyclostationary means that the statistics change for different time offsets but has cycles. Say the peak of your pulse has a value equal to $1.0$. In noise, you sample the peak every symbol (cycle $=$ symbol period here) and you get a bunch of values: $1.1, 0.95, 1.04, 0.98, 0.89$ and so on. You can then collect statistics on these samples like the mean and variance. Cyclostationary means that changing the sample time by a half a symbol period will give you different statistics, but changing the sample time by a multiple of the symbol period (remember cyclo/cyclical) will give the same statistics. This is different from stationary process, where the time shift does not at all affect the statistics.
"In the frequency domain four perfect images of the original spectrum are created due to the upsampling" how has upsampling caused this to happen?
You can see from the figure you posted, going from the plot labeled "Stationary Symbols" to "Upsampling" we are putting $3$ zeros between every symbol, or every symbol now occupies $4 \times$ as much samples (actual symbol + $3$ zeros = $4$ total samples). This is upsampling by a factor of $4$. When we do this operation, we get a compressed frequency domain axis and $4$ copies of the original spectrum (since we upsampled by $4$). I'm not deriving anything since this is a pretty standard DSP thing, check out
DSP Lecture 14: Continuous-time filtering with digital systems; upsampling and downsampling on Youtube starting at about the 56 minute mark for more in-depth upsampling explanation.