Generating a time series given a transfer function
I think that procedure should be as follows:

1.Generate the white (Gaussian) noise  signal - randn function in MATLAB.

2.Having a transfer function of your filter divide it by $z^4$ (highest power in denominator) to get a proper representation for IIR filter:

$$H_d(z)=\dfrac{0.1584\cdot z^{-1} −0.3765\cdot z^{-2} + 0.2716\cdot z^{-3} −0.0534\cdot z^{-4}}{1 −2.9951\cdot z^{-1} + 3.0893\cdot z^{-2} −1.1930\cdot z^{-3} + 0.0988\cdot z^{-4}}$$

3.Now you can extract the coefficients according to equations that can be found here: IIR.

$$ H(z) = \dfrac{\sum_{i=0}^{P}b_i z^{-1}}{1+\sum_{j=1}^{Q}a_j z^{-j}}$$

This will give you following vectors:

a=[1 -2.9951 3.0893 -1.1930 0.0988] % denominator
b=[0 0.1584 -0.3765 0.2716 -0.0534] % numerator


4.These are coefficients of your filter that you can apply to the white noise signal. If you are using MATLAB/Octave, then use the filter function. Something like:

y=filter(b,a,x)


where x is your white noise signal generated previously. 

The amplitude response of your filter is shown below - because you are going to apply it to white noise (having all frequencies equally probable), then spectrum of your filtered signal should be corresponding to amplitude response of your filter. There is a high gain in lower frequencies which makes somewhat sense (wind gusts are mostly energetic in lower range):



Judging from your article, now you only have to normalise the filtered white noise and scale it accordingly.
Good luck!