Statistics and Algorithms / Spectral Analysis /

Forecasting

Spectral methods can be used to forecast univariate time series. The @SPECFORE procedure implements the algorithm described here, which follows a suggestion made by John Geweke. The basis of the method is the moving average representation:

\begin{equation} X_t = c{\kern 1pt} {\kern 1pt} \left( L \right){\kern 1pt} {\kern 1pt} \varepsilon _t {\kern 1pt} {\kern 1pt} {\text{, where }}c{\kern 1pt} \left( 0 \right) = 1{\text{ and }}\varepsilon {\text{ is fundamental for }}X \end{equation}

Box-Jenkins techniques attempt to represent \(c\) as a rational function: the numerator is the moving average part, the denominator is the autoregressive. Spectral methods allow us to compute an estimate of the Fourier transform of \(c\), which can be used to compute forecasts. Their advantage over Box-Jenkins is that the same “model” applies to all data series, so the forecasting process is automatic. That, unfortunately, is one of the disadvantages as well—see the cautions below.

The following derivation can be found in Koopmans (1974), p. 235-237. Write the spectral density of \(X\) as the z-transform

\begin{equation} f_X \left( z \right) = c{\kern 1pt} \left( z \right){\kern 1pt} {\kern 1pt} c{\kern 1pt} \left( {z^{ - 1} } \right){\kern 1pt} {\kern 1pt} \sigma ^2 \label{eq:spectrum_fromMAR} \end{equation}

Under reasonable regularity conditions on \(c\), \(\log \left( {f_X } \right)\) has a Laurent expansion:

\begin{equation} \log {\kern 1pt} {\kern 1pt} {\kern 1pt} f_X \left( z \right) = d\left( z \right) + d\left( {z^{ - 1} } \right) + d_0 \label{eq:spectral_laurent} \end{equation}

where \(d\) is a one-sided polynomial in positive powers of \(z\). If we take the exp of both sides of \eqref{eq:spectral_laurent}, and combine it with \eqref{eq:spectrum_fromMAR}, we get:

\begin{equation} c{\kern 1pt} \left( z \right){\kern 1pt} {\kern 1pt} c{\kern 1pt} \left( {z^{ - 1} } \right){\kern 1pt} {\kern 1pt} \sigma ^2 = \exp \left( {d\left( z \right)} \right)\exp \left( {d\left( {z^{ - 1} } \right)} \right)\exp \left( {d_0 } \right) \label{eq:spectral_expLaurent} \end{equation}

To get the Fourier transform of \(c\) we need to

1.compute the log spectral density of \(X\)

2.mask the negative and zero frequencies to get the Fourier transform of \(d\)

3.take its exp, frequency by frequency (again, omitting the negative and zero frequencies)

We now have a “non-parametric” estimate of the moving average representation of the series. Once we have an estimate of the FT of \(c\), we divide the FT of \(X\) by the FT of \(c\) to get the FT of \(e\) process, set \(e\) to zero outside of the range of actual data, then filter \(e\) by \(c\) (in the frequency domain) to get the forecasts.

Note that the formulas \eqref{eq:spectral_laurent} and \eqref{eq:spectral_expLaurent} apply only to univariate processes—there’s no straightforward extension of this to multivariate processes.

Cautions

As with other frequency domain filtering techniques, the series must be well–padded. Also, please note the following:

•Data series which can be adequately represented by a Box-Jenkins model with few parameters will usually be better forecast using that technique, as the model is more sharply estimated.

•Spectral methods require more data (100+ data points to be safe) than other techniques, since they do not use it as efficiently.

•The series you run through the model must be stationary: proper differencing is vital.

Code (Out of SPECFORE.SRC)

Compute the log spectral density. Scale factors will wash out.

fft 1

cmult 1 1 / 2

window 2 / 3

cln 3

Back to sequence-space to shut out the negative powers in the Laurent expansion. FT back and exponentiate to get FT of exp(D(Z)).

ift 3

cset 3 = %z(t,3)*(t<=ordinate/2)

fft 3

cexp 3

IFT to get the sequence c. Normalize so that C(0)=1

ift 3

cset(scratch) 3 = %z(t,3)/%z(1,3)

Filter the input series by 1/C(L) to get residuals. Mask it outside the data range.

fft 3

cset 1 = %z(t,1)/%z(t,3)

ift 1 / 2

cset 2 = %z(t,2)*(t>=startl.and.t<=endl)

Filter the residuals E by C(L) in the frequency domain

fft 2 / 1

cset 1 = %z(t,1)*%conjg(%z(t,3))

ift 1