Nonlinear prediction via Hermite transformation

ABSTRACT

General prediction formulas involving Hermite polynomials are developed for time series expressed as a transformation of a Gaussian process. The prediction gains over linear predictors are examined numerically, demonstrating the improvement of nonlinear prediction.

Keywords:

• Gaussian transformation
• Hermite polynomial
• nonlinear forecasting
• time series

1. Introduction

The general prediction problem is to compute the best predictor of a random variable Y given a data vector $\underset{\_}{X}$, where a joint distribution is presumed to exist for Y and $\underset{\_}{X}$. If we define ‘best’ according to mean squared error (MSE) loss,¹ the best predictor (when the random variables are square integrable) is the conditional expectation $\mathbb{E}[Y|\underset{\_}{X}]$, which in the case of Gaussian variables is a linear function of $\underset{\_}{X}$. This linear function is completely computable in terms of first and second moments of the joint vector $(Y,\underset{\_}{X})$, as discussed in Brockwell and Davis (2013, Chapter 2). The problem can also be generalised to projection on infinite data sets, which arise in forecasting and signal extraction problems.

The theory for linear predictors is very well understood and is commonly applied to non-Gaussian data because it is simple to compute. Nevertheless, there can be a substantial predictive loss when non-Gaussian features are present in the data, such as asymmetry and excess kurtosis (Brockett et al., 1988; Maravall, 1983). A common technique for handling such raw time series data is to apply a transformation that reduces asymmetry and kurtosis, thereby generating cumulants in the transformed space that more closely resemble Gaussian cumulants. Box-Cox transforms are an example of such functions, and are typically identified through exploratory analysis or via metadata; see discussion in McElroy (2016).

This paper provides exact formulas for non-linear prediction in scenarios where the non-Gaussian data process can be expressed as a univariate transformation of some Gaussian process. For some special cases, such as the log-normal distribution, exact formulas are already available for nonlinear predictors; here the general case is developed. The main result of the paper (Section 2) is an expansion of the conditional expectation in terms of Hermite coefficients of the transformation function – an idea that was utilised in Janicki and McElroy (2016) to model marginal quantiles. Here this technique is used to derive analytical expressions for predictors, along with the MSE; the solution is given as an explicit function of the Hermite coefficients of the transformation function, the various Hermite polynomials evaluated at the linear predictor, and further weights explicitly determined from the mean squared error of the linear predictor.

These results are general, in the sense that they can be applied to diverse contexts in statistics, such as linear models, spatial statistics, and multivariate analysis. But our applications are focussed on time series, and in particular on forecasting problems. In forecasting, the data vector $\underset{\_}{X}={[X_1,\dots ,X_T]}^{\mathrm{\prime }}$ is a sample of size T from a time series $\{X_t\}$, and $Y=X_{T+1}$ represents the next unobserved value of the process. Backcasting involves setting $Y=X_0$, and missing value problems can similarly be addressed by letting $Y=X_t$ and $\underset{\_}{X}={[X_1,\dots ,X_{t-1},X_{t+1},\dots ,X_T]}^{\mathrm{\prime }}$; see McElroy and McCracken (2017) for a recent treatment. These facets of the general methodology are developed in Section 2, and numerical comparisons are given in Section 3. The proofs are in an Appendix.

2. Nonlinear prediction

Our goal is to compute the minimal mean squared error (MSE) estimate of Y given data $\underset{\_}{X}={[X_1,\dots ,X_T]}^{\mathrm{\prime }}$ (finite-sample case) or $\underset{\_}{X}=\{X_s:s\le T\}$ (semi-infinite sample case). This estimate is the conditional expectation $\mathbb{E}[Y|\underset{\_}{X}]$, denoted $\hat{Y}$ for short. We presume that there exists an invertible function f such that $Z_t=f(X_t)$ yields a Gaussian process $\{Z_t\}$. Moreover, $f(Y)$ is a Gaussian random variable whose joint distribution with the process $\{Z_t\}$ is known. In the case of forecasting, backcasting, or missing value imputation, $f(Y)=Z_t$ for some t.

In the following development, it is important to impose that $f(Y)$ be standard normal, although this would rarely be the case if f is obtained by exploratory analysis. (For example, if $Y=X_{T+1}$ and $f(x)=\log (x)$ by exploratory analysis, it would rarely be the case that $\log (X_{T+1})$ would have unit variance.) Let $Z_{\star }$ denote the standardisation of $f(Y)$, i.e., $Z_{\star }=(f(Y)-\mathbb{E}[f(Y)])/\sqrt{\text{Var}[f(Y)]}$. We will define g as the inverse map from $Z_{\star }$ to Y, so that $\begin{array}{r}g(x)=f^{-1}(x\sqrt{\text{Var}[f(Y)]}+\mathbb{E}[f(Y)]).\end{array}$The mapping $Y=g(Z_{\star })$ allows us to obtain a Hermite expansion of g; if the marginal distribution of $Z_{\star }$ were non-Gaussian, we could instead have recourse to the Appell polynomials (Varma, 1951).

The main result of the paper is an expression for $\mathbb{E}[Y|\underset{\_}{X}]$ in terms of $\mathbb{E}[Z_{\star }|\underset{\_}{Z}]$, with $\underset{\_}{Z}={[Z_1,\dots ,Z_T]}^{\mathrm{\prime }}$. This is of interest because such a Gaussian conditional expectation has a well-known linear formula; see Chapter 4 of McElroy and Politis (2020). In applications, one might transform the data by applying f, then model the Gaussian process, and finally compute $\mathbb{E}[Z_{\star }|\underset{\_}{Z}]$ by plugging into the linear formulas. Then the main formula of this paper can be used to obtain $\mathbb{E}[Y|\underset{\_}{X}]$ by inserting $\mathbb{E}[Z_{\star }|\underset{\_}{Z}]$ and its MSE, denoted by V, which would also be available in applications. In particular, (1) $\begin{array}{r}V=\mathbb{E}\left[{(Z_{\star }-\mathbb{E}[Z_{\star }|\underset{\_}{Z}])}^2\right]\end{array}$(1) and is given by formulas in McElroy and Politis (2020); also, see (11(11) $\begin{array}{rl}V& ={\gamma }_Z(0)-[{\gamma }_Z(1),\dots ,{\gamma }_Z(T)]{\mathrm{\Gamma }}_Z^{-1}\\ & \times [{\gamma }_Z(1),\dots ,{\gamma }_Z(T){]}^{\mathrm{\prime }}\end{array}$(11) ) below. We define the space ${\mathbb{L}}_2(\mathrm{d}\mathrm{\Phi })$ (where Φ is the standard normal cumulative distribution function) as all functions that are square integrable with respect to the measure $\mathrm{d}\mathrm{\Phi }$. An inner product on this space is defined via $⟨f,h⟩={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(x)h(x)\varphi (x)\mathrm{d}x$, where ϕ is the standard normal probability density function. Hence, we can also write $⟨f,g⟩=\mathbb{E}[f(W)g(W)]$ for $W\sim \mathcal{N}(0,1)$. It follows that we can do a Hermite expansion on g (see Janicki & McElroy, 2016 for background): (2) $\begin{array}{r}g(x)=\sum _{k=0}^{\mathrm{\infty }}{J}_k{H}_k(x)\end{array}$(2) with $H_k$ the normalised Hermite polynomials (Roman, 1984), and $J_k=⟨g,H_k⟩$ the Hermite coefficients. The Hermite polynomials are defined via $\begin{array}{r}{H}_k(x)=\frac{1}{\sqrt{k!}}{(-1)}^k{\mathrm{e}}^{x^2/2}{\mathrm{\partial }}_x^k{\mathrm{e}}^{-x^2/2}\end{array}$for $k\ge 0$, and form a complete orthonormal system. (Hence, the coefficients $J_k$ tend to zero as $k\to \mathrm{\infty }$.) Plugging $Z_{\star }$ into (2(2) $\begin{array}{r}g(x)=\sum _{k=0}^{\mathrm{\infty }}{J}_k{H}_k(x)\end{array}$(2) ) yields $Y=\sum _{k=0}^{\mathrm{\infty }}{J}_k{H}_k(Z_{\star })$, and applying the conditional expectation operator (which is linear) yields (3) $\begin{array}{r}\hat{Y}=\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}{J}_k\mathbb{E}[H_k(Z_{\star })|\underset{\_}{X}].\end{array}$(3) Evidently, the nonlinear predictor can be computed in terms of conditional expectations of $H_k(Z_{\star })$, and the formula is summarised in our main theorem below. The proof relies upon the Hermite generating function (4) $\begin{array}{r}h(x,t)=\exp \{xt-t^2/2\}=\sum _{k=0}^{\mathrm{\infty }}\frac{t^k}{\sqrt{k!}}{H}_k(x).\end{array}$(4) (The second equality follows from the definition of the Hermite polynomials, and is discussed in Roman (1984).) From (4(4) $\begin{array}{r}h(x,t)=\exp \{xt-t^2/2\}=\sum _{k=0}^{\mathrm{\infty }}\frac{t^k}{\sqrt{k!}}{H}_k(x).\end{array}$(4) ) we see that (5) $\begin{array}{r}{\mathrm{\partial }}_t^{k}h(x,t){|}_{t=0}=\sqrt{k!}{H}_k(x).\end{array}$(5) Therefore, $J_k=\mathbb{E}[g(W)H_k(W)]=\frac{1}{\sqrt{k!}}{\mathrm{\partial }}_t^k\mathbb{E}[g(W)h(W,t)]{|}_{t=0}$. We next discuss the optimal predictor.

Theorem 2.1

Suppose that $Z_{\star }$ is standard normal, and $Y=g(Z_{\star })$ with g given by (2(2) $\begin{array}{r}g(x)=\sum _{k=0}^{\mathrm{\infty }}{J}_k{H}_k(x)\end{array}$(2) ), and $\widehat{Z_{\star }}=\mathbb{E}[Z_{\star }|\underset{\_}{Z}]$ is the linear prediction in the transformed space, with MSE given by V in (1(1) $\begin{array}{r}V=\mathbb{E}\left[{(Z_{\star }-\mathbb{E}[Z_{\star }|\underset{\_}{Z}])}^2\right]\end{array}$(1) ). Then the MSE optimal nonlinear predictor $\hat{Y}$ of Y given the data $\underset{\_}{X}$ is (6) $\begin{array}{r}\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}\sum _{\ell =0}^k(\genfrac{}{}{0ex}{}{k}{\ell })\sqrt{\ell !}{H}_{\ell }(\widehat{Z_{\star }}){\kappa }_{k-\ell },\end{array}$(6) where ${\kappa }_{\ell }$ is the ℓth moment of a Gaussian variable of variance V, i.e., ${\kappa }_{\ell }$ equals 1 if $\ell =0$, and equals zero if ℓ is odd and equals ${\sqrt{V}}^{\ell }(\ell -1)!!$ if ℓ is even. With $\epsilon =Y-\hat{Y}$, the optimal prediction MSE is (7) $\begin{array}{r}\mathbb{E}[{\epsilon }^2]=\sum _{k=1}^{\mathrm{\infty }}{J}_k^2(1-{(1-V)}^k).\end{array}$(7)

The optimal predictor (6(6) $\begin{array}{r}\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}\sum _{\ell =0}^k(\genfrac{}{}{0ex}{}{k}{\ell })\sqrt{\ell !}{H}_{\ell }(\widehat{Z_{\star }}){\kappa }_{k-\ell },\end{array}$(6) ) can be explicitly computed to any desired level of accuracy, truncating the first summation over k to a desirable level. Note that this result applies to finite-samples, semi-infinite samples, and bi-infinite samples, allowing us to apply (6(6) $\begin{array}{r}\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}\sum _{\ell =0}^k(\genfrac{}{}{0ex}{}{k}{\ell })\sqrt{\ell !}{H}_{\ell }(\widehat{Z_{\star }}){\kappa }_{k-\ell },\end{array}$(6) ) to forecasting (from an infinite past) and missing value problems.

Remark 2.1

Since $\sum _{k=0}^{\mathrm{\infty }}{J}_k^2=⟨g,g⟩$, we see that the Hermite coefficients are square summable if $g\in {\mathbb{L}}_2(\mathrm{d}\mathrm{\Phi })$, and the MSE is well-defined, being bounded above by $\sum _{k=1}^{\mathrm{\infty }}{J}_k^2{V}^k$. (Note that $V\le 1$ because $Z_{\star }$ has variance one.) Because $\text{Var}[{\hat{Y}}^2]=\text{Var}[Y^2]-\mathbb{E}[{\epsilon }^2]$, we see that $\hat{Y}$ is square integrable, and hence the formula (6(6) $\begin{array}{r}\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}\sum _{\ell =0}^k(\genfrac{}{}{0ex}{}{k}{\ell })\sqrt{\ell !}{H}_{\ell }(\widehat{Z_{\star }}){\kappa }_{k-\ell },\end{array}$(6) ) converges.

Remark 2.2

In applications, we typically will know $f(Y)|\underset{\_}{Z}$ rather than $Z_{\star }|\underset{\_}{Z}$, but the latter can be obtained from the former by applying the normalising transform described above, i.e., $\begin{array}{rl}\mathbb{E}[Z_{\star }|\underset{\_}{Z}]& =(\mathbb{E}[f(Y)|\underset{\_}{Z}]-\mathbb{E}[f(Y)])/\sqrt{\text{Var}[f(Y)]}\\ V& =\mathbb{E}[{(f(Y)-\mathbb{E}[f(Y)|\underset{\_}{Z}])}^2]/\text{Var}[f(Y)].\end{array}$

Remark 2.3

To apply Theorem 2.1 we must obtain the Hermite coefficients $J_k$. Either g is known (obtained by exploratory analysis) or estimated (see Janicki & McElroy, 2016), and then the Hermite coefficients can be obtained via Monte Carlo simulation: $\begin{array}{rl}{J}_k=⟨g,H_k⟩& =\mathbb{E}[g(W)H_k(W)]\\ & \approx M^{-1}\sum _{m=1}^{M}g(W_m)H_k(W_m)\end{array}$for $W_1,\dots ,W_m$ i.i.d. standard normal. Another approach to computation involves the generating function: $\begin{array}{rl}{J}_k& =\frac{1}{\sqrt{k!}}{\mathrm{\partial }}_t^k\mathbb{E}[g(W)h(W,t)]{|}_{t=0}\\ & =\frac{1}{\sqrt{2\pi k!}}{\mathrm{\partial }}_t^k\int g(w)\exp \{wt-t^2/2\}\\ & \times \exp \{-w^2/2\}\mathrm{d}w{|}_{t=0}\\ & =\frac{1}{\sqrt{2\pi k!}}{\mathrm{\partial }}_t^k\int g(w)\exp \{-{(w-t)}^2/2\}\mathrm{d}w{|}_{t=0}\\ & =\frac{1}{\sqrt{2\pi k!}}{\mathrm{\partial }}_t^k\int g(w+t)\exp \{-w^2/2\}\mathrm{d}w{|}_{t=0}\\ & =\frac{1}{\sqrt{k!}}{\mathrm{\partial }}_t^k\mathbb{E}[g(W+t)]{|}_{t=0}=\frac{1}{\sqrt{k!}}\mathbb{E}[g^{(k)}(W)].\end{array}$The last expression denotes the k-fold derivative of the function. In cases where g is explicitly known, this can be an easier approach to getting the Hermite coefficients.

In the following examples we suppose that ${\hat{Z}}_{\star }=\mathbb{E}[Z_{\star }|\underset{\_}{Z}]$ and $V=\mathbb{E}[{({\hat{Z}}_{\star }-Z_{\star })}^2]$ are known and available to the practitioner; we present various cases of transforms f, and apply the results of Theorem 2.1.

Example 2.1

Gaussian

A simple affine transformation $g(x)=\sigma x+\mu$ ensures that Y is still Gaussian, with mean μ and variance ${\sigma }^2$. In this case $J_0=\mu$ and $J_1=\sigma$, and we more simply have $\mathbb{E}[Y|\underset{\_}{X}]=\sigma \widehat{Z_{\star }}+\mu$.

Example 2.2

Lognormal

Suppose $g(x)={\mathrm{e}}^x$; applying the generating function method of Remark 2.3, we obtain $\begin{array}{r}{J}_k=\frac{1}{\sqrt{k!}}\mathbb{E}[\exp \{W\}]=\frac{{\mathrm{e}}^{1/2}}{\sqrt{k!}}\end{array}$for all $k\ge 0$. This can be utilised in (3(3) $\begin{array}{r}\hat{Y}=\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}{J}_k\mathbb{E}[H_k(Z_{\star })|\underset{\_}{X}].\end{array}$(3) ), together with (5(5) $\begin{array}{r}{\mathrm{\partial }}_t^{k}h(x,t){|}_{t=0}=\sqrt{k!}{H}_k(x).\end{array}$(5) ), to yield $\begin{array}{rl}\mathbb{E}[Y|\underset{\_}{X}]& ={\mathrm{e}}^{1/2}\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}{\mathrm{\partial }}_t^k(h(\widehat{Z_{\star }},t){\mathrm{e}}^{V{t}^2/2}){|}_{t=0}\\ & ={\mathrm{e}}^{1/2}(h(\widehat{Z_{\star }},t){\mathrm{e}}^{V{t}^2/2}){|}_{t=1}\\ & ={\mathrm{e}}^{1/2}\exp \{\widehat{Z_{\star }}+(V-1)/2\}\\ & =\exp \{\widehat{Z_{\star }}+V/2\},\end{array}$which corresponds to the result of McElroy (2010). Applying Theorem 2.1 to compute the optimal MSE, we see that $\mathbb{E}[{(\hat{Y}-Y)}^2]={\mathrm{e}}^2(1-{\mathrm{e}}^{-V}).$

Example 2.3

Uniform

For $\{Z_t\}$ with standard normal marginal, set $g=\mathrm{\Phi }$, so that $\{X_t\}$ has a marginal distribution that is uniform on $(0,1)$. Then by Remark 2.3, we find that $g^{(k)}(x)={\varphi }^{(k-1)}(x)$, and hence $\begin{array}{rl}{J}_k& =\frac{1}{\sqrt{k!}},\mathbb{E}[{\varphi }^{(k-1)}(W)]\\ & =k^{-1/2}{(-1)}^{k-1}\mathbb{E}[H_{k-1}(W)\varphi (W)],\end{array}$where $W\sim \mathcal{N}(0,1)$.

Example 2.4

Logistic

Consider a logistic transform given by $g(x)={\mathrm{e}}^x/(1+{\mathrm{e}}^x)$. The first few derivatives of g are $\begin{array}{rl}{g}^{(1)}(x)& ={\mathrm{e}}^x(1+{\mathrm{e}}^x{)}^{-2}\\ g^{(2)}(x)& =({\mathrm{e}}^x-{\mathrm{e}}^{2x})(1+{\mathrm{e}}^x{)}^{-3}\\ g^{(3)}(x)& =({\mathrm{e}}^x-4{\mathrm{e}}^{2x}+{\mathrm{e}}^{3x})(1+{\mathrm{e}}^x{)}^{-4}.\end{array}$By Monte Carlo, we obtain $J_0=0.500$, $J_1=0.207$, $J_2=0$, and $J_3=-.025$.

Example 2.5

Square

With $g(x)=x^2$, we have $J_0=1$, $J_1=0$, $J_2=\sqrt{2}$, and $J_k=0$ for k>2. Then the optimal nonlinear predictor is $\mathbb{E}[Y|\underset{\_}{X}]={\widehat{Z_{\star }}}^2+V$. It is simple to check that the error is $\begin{array}{r}ϵ=(Z_{\star }-\widehat{Z_{\star }})(Z_{\star }+\widehat{Z_{\star }})-V,\end{array}$which has mean zero and is independent of all functions of the data. The prediction MSE is $4V-2V^2$.

3. Comparing linear and nonlinear prediction

It is of interest to understand how much benefit nonlinear prediction provides. Clearly, if g is affine (Example 2.1) then the minimal MSE is equal to ${\sigma }^{2}V$, the same as the linear prediction, but we can expect gains to the degree that g differs from the affine case.

Remark 3.1

A related nonlinear predictor that is sometimes used in applications is defined via $\tilde{Y}=g(\hat{Z})$, but unfortunately this estimator can be biased. Following the same arguments used in the proof of Theorem 2.1, $\begin{array}{rl}Y-\tilde{Y}& =\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}{\mathrm{\partial }}_t^k(h(Z,t)-h(\widehat{Z_{\star }},t)){|}_{t=0}\\ & =\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}(\exp \{\widehat{Z_{\star }}t-t^2/2\}\\ & \times [\exp \{(Z-\widehat{Z_{\star }}t\}-1]){|}_{t=0},\end{array}$so that the expectation of the quantity in parentheses is $\begin{array}{r}\mathbb{E}[\exp \{\widehat{Z_{\star }}t-t^2/2\}][\exp \{V{t}^2/2\}-1]>0.\end{array}$Hence there is no guarantee that the bias is zero.

More properly, a comparison can be made to the best linear predictor. When $\{Z_t\}$ is a strictly stationary time series, then $\{X_t\}$ is as well, and we can in some cases determine the best linear estimator's MSE for comparison. Note that in this special case, the mean and variance of each variable $Z_t$ is the same, and hence without any loss of generality we may assume that $\{Z_t\}$ is standardised, i.e., each $Z_t$ is standard normal. Because $X_t=g(Z_t)$, it follows from Taniguchi and Kakizawa (2000, p. 319) that $\mathbb{E}[X_t]=J_0$ and (8) $\begin{array}{r}{\gamma }_X(h)=\sum _{k=0}^{\mathrm{\infty }}{J}_k^2{{\gamma }_Z(h)}^k-J_0^2=\sum _{k=1}^{\mathrm{\infty }}{J}_k^2{{\gamma }_Z(h)}^k,\end{array}$(8) where $\{{\gamma }_X(h)\}$ and $\{{\gamma }_Z(h)\}$ are the autocovariance sequences of $\{X_t\}$ and $\{Z_t\}$, respectively. So in principle we can understand the second order structure of $\{X_t\}$ in terms of the Hermite coefficients and the original autocovariances.

Suppose further that we are interested in one-step ahead forecasting from a sample of size T. The best linear predictor is obtained by solving the Yule-Walker equations in $\{{\gamma }_X(h)\}$, and the MSE of such is given by (9) $\begin{array}{r}{\gamma }_X(0)-[{\gamma }_X(1),\dots ,{\gamma }_X(T)]{\mathrm{\Gamma }}_X^{-1}[{\gamma }_X(1),\dots ,{\gamma }_X(T){]}^{\mathrm{\prime }},\end{array}$(9) where ${\mathrm{\Gamma }}_X$ is the T-dimensional Toeplitz covariance matrix with jkth entry ${\gamma }_X(j-k)$. We know that such an MSE must be greater than the minimal MSE provided in Theorem 2.1, with equality occuring only in the case that a linear estimator is globally optimal (e.g., the time series is linear, or is Gaussian). If instead we are forecasting from an infinite past, then the MSE of the linear predictor is the innovation variance ${\sigma }^2$ given by Kolmogorov's formula: (10) $\begin{array}{r}{\sigma }^2=\exp \{{(2\pi )}^{-1}{\int }_{-\pi }^{\pi }\log f_X(\lambda )\mathrm{d}\lambda \}.\end{array}$(10) Here, $f_X(\lambda )=\sum _{h=-\mathrm{\infty }}^{\mathrm{\infty }}{\gamma }_X(h){\mathrm{e}}^{-i\lambda h}$ is the spectral density of $\{X_t\}$. Therefore, for either a finite sample or for an infinite past, we can determine the linear predictor MSE for a stationary process $\{X_t\}$ by first computing ${\gamma }_X(h)$ from ${\gamma }_Z(h)$ via (8(8) $\begin{array}{r}{\gamma }_X(h)=\sum _{k=0}^{\mathrm{\infty }}{J}_k^2{{\gamma }_Z(h)}^k-J_0^2=\sum _{k=1}^{\mathrm{\infty }}{J}_k^2{{\gamma }_Z(h)}^k,\end{array}$(8) ), followed by application of (9(9) $\begin{array}{r}{\gamma }_X(0)-[{\gamma }_X(1),\dots ,{\gamma }_X(T)]{\mathrm{\Gamma }}_X^{-1}[{\gamma }_X(1),\dots ,{\gamma }_X(T){]}^{\mathrm{\prime }},\end{array}$(9) ) or (10(10) $\begin{array}{r}{\sigma }^2=\exp \{{(2\pi )}^{-1}{\int }_{-\pi }^{\pi }\log f_X(\lambda )\mathrm{d}\lambda \}.\end{array}$(10) ) as each case requires. As for the best (nonlinear) predictor, its MSE is given by (7(7) $\begin{array}{r}\mathbb{E}[{\epsilon }^2]=\sum _{k=1}^{\mathrm{\infty }}{J}_k^2(1-{(1-V)}^k).\end{array}$(7) ) of Theorem 2.1, where (11) $\begin{array}{rl}V& ={\gamma }_Z(0)-[{\gamma }_Z(1),\dots ,{\gamma }_Z(T)]{\mathrm{\Gamma }}_Z^{-1}\\ & \times [{\gamma }_Z(1),\dots ,{\gamma }_Z(T){]}^{\mathrm{\prime }}\end{array}$(11) is the analogue of (9(9) $\begin{array}{r}{\gamma }_X(0)-[{\gamma }_X(1),\dots ,{\gamma }_X(T)]{\mathrm{\Gamma }}_X^{-1}[{\gamma }_X(1),\dots ,{\gamma }_X(T){]}^{\mathrm{\prime }},\end{array}$(9) ) for the $\{Z_t\}$ process.

We provide an illustration in the case of an MA(1) process with various values of θ, and sample size T = 100. The innovation variance is set equal to ${(1+{\theta }^2)}^{-1}$ so that ${\gamma }_Z(0)=1$, as required by the above discussion. For the MA(1) process with T = 100, the value of V given by (11(11) $\begin{array}{rl}V& ={\gamma }_Z(0)-[{\gamma }_Z(1),\dots ,{\gamma }_Z(T)]{\mathrm{\Gamma }}_Z^{-1}\\ & \times [{\gamma }_Z(1),\dots ,{\gamma }_Z(T){]}^{\mathrm{\prime }}\end{array}$(11) ) is the same (up to the fourth decimal place) as the innovation variance ${(1+{\theta }^2)}^{-1}$. For transformations, we study $g(x)=x^2$, $g(x)={\mathrm{e}}^x$, and the logistic. Observe that from (8(8) $\begin{array}{r}{\gamma }_X(h)=\sum _{k=0}^{\mathrm{\infty }}{J}_k^2{{\gamma }_Z(h)}^k-J_0^2=\sum _{k=1}^{\mathrm{\infty }}{J}_k^2{{\gamma }_Z(h)}^k,\end{array}$(8) ) the process $\{X_t\}$ will be m-dependent if $\{Z_t\}$ is (although the converse need not be true). Hence, if we obtained a sample from $\{X_t\}$ it would likely be identified with an MA(q) model, and the parameter estimates (e.g., obtained using a Whittle likelihood, which is valid for non-Gaussian processes so long as the cumulants are summable) would likely converge to those corresponding to the spectral factorisation of $f_X$. Thus, our illustration provides an accurate rendition of the prediction MSE one would obtain in the case of linear or nonlinear predictors, only with the impact of parameter estimation error completely removed (Tables 1–3).

Table 1. MSE for linear and non-linear predictors applied to a squared MA(1) process of parameter θ.

θ	0	$.2$	$.4$	$.6$	$.8$
V	$1.000$	0.9615	$0.8621$	0.7353	0.6098
Linear	$2.000$	1.9973	$1.9713$	$1.9211$	$1.8795$
Non-linear	$2.000$	1.9970	$1.9620$	$1.8599$	$1.6954$

Table 2. MSE for linear and non-linear predictors applied to an exponential MA(1) process of parameter θ.

θ	0	$.2$	$.4$	$.6$	$.8$
V	$1.000$	0.9615	$0.8621$	0.7353	0.6098
Linear	$4.6708$	4.5985	$4.3851$	$4.1192$	$3.9269$
Non-linear	$4.6708$	4.5642	$4.2688$	$3.8470$	$3.3732$

Table 3. MSE for linear and non-linear predictors applied to a logistic MA(1) process of parameter θ.

θ	0	$.2$	$.4$	$.6$	$.8$
V	$1.000$	0.9615	$0.8621$	0.7353	0.6098
Linear	$0.0433$	0.0417	$0.0375$	0.0323	$0.0274$
Non-linear	$0.0433$	0.0417	$0.0374$	0.0320	$0.0266$

In each case, the degree of benefit to nonlinear prediction increases with θ, as to be expected; however, there are large discrepancies between the three functions. When $\theta =.8$, the logistic transformation offers only a $3\mathrm{\%}$ improvement with nonlinear prediction, indicating that linear prediction is almost just as good as the conditional expectation. With $g(x)=x^2$ the analogous improvement is $11\mathrm{\%}$, and is $16\mathrm{\%}$ for $g(x)={\mathrm{e}}^x$, indicating some real benefit to nonlinear prediction.

We end with a remark on how a confidence interval can be constructed using simulations. Let the formula given by (6(6) $\begin{array}{r}\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}\sum _{\ell =0}^k(\genfrac{}{}{0ex}{}{k}{\ell })\sqrt{\ell !}{H}_{\ell }(\widehat{Z_{\star }}){\kappa }_{k-\ell },\end{array}$(6) ) be denoted $h({\hat{Z}}_{\star })$. Then the optimal prediction error can be written $\begin{array}{r}\epsilon =g({\hat{Z}}_{\star }+\delta )-h({\hat{Z}}_{\star }),\end{array}$where $\delta =Z_{\star }-{\hat{Z}}_{\star }$ is the (linear) prediction error for the Gaussian variable, and is uncorrelated with (and hence independent of) ${\hat{Z}}_{\star }$. Also $\delta \sim \mathcal{N}(0,V)$. Therefore, in cases where it is easy to simulate ${\hat{Z}}_{\star }$ (e.g., suppose we are forecasting from a Gaussian ARMA process) we can independently draw δ and compute ε for repeated Monte Carlo draws, thereby obtaining a confidence interval for Y.

Disclaimer

This report is released to inform interested parties of research and to encourage discussion. The views expressed on statistical issues are those of the authors and not those of the U.S. Census Bureau.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Notes on contributors

Tucker McElroy

Tucker McElroy is Senior Time Series Mathematical Statistician at the U.S. Census Bureau.

Srinjoy Das

Srinjoy Das is a Postdoctoral Scholar in the Mathematics department at the University of California, San Diego.

Notes

1 Other loss functions could of course be envisioned (mean absolute loss yields the conditional median, for example).

Appendix

Proof of Theorem 2.1

From (3(3) $\begin{array}{r}\hat{Y}=\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}{J}_k\mathbb{E}[H_k(Z_{\star })|\underset{\_}{X}].\end{array}$(3) ) and (5(5) $\begin{array}{r}{\mathrm{\partial }}_t^{k}h(x,t){|}_{t=0}=\sqrt{k!}{H}_k(x).\end{array}$(5) ) we obtain $\begin{array}{r}\mathbb{E}[H_k(Z_{\star })|\underset{\_}{X}]=\frac{1}{\sqrt{k!}}{\mathrm{\partial }}_t^k\mathbb{E}[h(Z_{\star },t)|\underset{\_}{X}]{|}_{t=0}.\end{array}$We can write $Z_{\star }|\underset{\_}{X}\sim \mathcal{N}(\widehat{Z_{\star }},V)$, and using the property that $Z_{\star }-\widehat{Z_{\star }}$ is independent of all functions of the data, we obtain $\begin{array}{rl}\mathbb{E}[h(Z_{\star },t)|\underset{\_}{X}]& =\mathbb{E}[\exp \{Z_{\star }t-\widehat{Z_{\star }}t\}\exp \{\widehat{Z_{\star }}t-t^2/2\}|\underset{\_}{X}]\\ & =\mathbb{E}[\exp \{Z_{\star }t-\widehat{Z_{\star }}t\}|\underset{\_}{X}]\exp \{\widehat{Z_{\star }}t-t^2/2\}\\ & =\exp \{V{t}^2/2\}h(\widehat{Z_{\star }},t).\end{array}$Hence $\begin{array}{rl}\mathbb{E}[H_k(Z_{\star })|\underset{\_}{X}]& =\frac{1}{\sqrt{k!}}{\mathrm{\partial }}_t^k(h(\widehat{Z_{\star }},t){\mathrm{e}}^{V{t}^2/2}){|}_{t=0}\\ & =\frac{1}{\sqrt{k!}}\sum _{\ell =0}^k(\genfrac{}{}{0ex}{}{k}{\ell }){\mathrm{\partial }}_t^{\ell }h(\widehat{Z_{\star }},t){\mathrm{\partial }}_t^{k-\ell }{\mathrm{e}}^{V{t}^2/2}{|}_{t=0}\\ & =\frac{1}{\sqrt{k!}}\sum _{\ell =0}^k(\genfrac{}{}{0ex}{}{k}{\ell })\sqrt{\ell !}{H}_{\ell }(\widehat{Z_{\star }}){\kappa }_{k-\ell }.\end{array}$The prediction error is $\begin{array}{rl}\epsilon & =Y-\mathbb{E}[Y|\underset{\_}{X}]=\sum _{k=0}^{\mathrm{\infty }}{J}_k(H_k(Z_{\star })-\mathbb{E}[H_k(Z_{\star })|\underset{\_}{X}])\\ & =\sum _{k=1}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}{\mathrm{\partial }}_t^k(h(Z_{\star },t)-\mathbb{E}[h(Z_{\star },t)|\underset{\_}{X}]){|}_{t=0}\\ & =\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}{\mathrm{\partial }}_t^k(h(\widehat{Z_{\star }},t)\\ & \cdot [\exp \{(Z_{\star }-\widehat{Z_{\star }})t\}-\exp \{V{t}^2/2\}]){|}_{t=0}.\end{array}$Note that $Z_{\star }-\widehat{Z_{\star }}$ is orthogonal to all linear functions of the data; because this error is Gaussian, it is moreover independent of all functions of the data $\underset{\_}{X}$. It follows that $\mathbb{E}[\epsilon ]=0$, because $\mathbb{E}[\exp \{(Z_{\star }-\widehat{Z_{\star }})t\}]=\exp \{V{t}^2/2\}.$ Moreover, for any function $\ell (\underset{\_}{X})$ of the data, $\begin{array}{rl}\mathbb{E}[ϵ\ell (\underset{\_}{X})]& =\sum _{k=0}^{\mathrm{\infty }}\frac{J_k}{\sqrt{k!}}{\mathrm{\partial }}_t^k(h(\widehat{Z_{\star }},t)\ell (\underset{\_}{X})\\ & \cdot \mathbb{E}[\exp \{(Z_{\star }-\widehat{Z_{\star }})t\}\\ & -\exp \{V{t}^2/2\}]){|}_{t=0}=0.\end{array}$This verifies optimality. To compute the MSE, first observe that $\begin{array}{rl}{\epsilon }^2& =\sum _{j,k=1}^{\mathrm{\infty }}\frac{J_j{J}_k}{\sqrt{j!}\sqrt{k!}}{\mathrm{\partial }}_s^j{\mathrm{\partial }}_t^k(\exp \{\widehat{Z_{\star }}(s+t)-(s^2+t^2)/2\}\\ & \cdot [\exp \{(Z_{\star }-\widehat{Z_{\star }})(s+t)\}\\ & -\exp \{V{t}^2/2+(Z_{\star }-\widehat{Z_{\star }})s\}\\ & -\exp \{V{s}^2/2+(Z_{\star }-\widehat{Z_{\star }})t\}\\ & +\exp \{V(s^2+t^2)/2\}]){|}_{s,t=0}.\end{array}$Note that $\mathbb{E}[\widehat{Z_{\star }}]=\mathbb{E}[Z_{\star }]=0$, because $Z_{\star }$ is standard normal. Moreover, due to orthogonality, $Z_{\star }=(Z_{\star }-\widehat{Z_{\star }})+\widehat{Z_{\star }}$ with the two summands orthogonal, and hence $1=\mathbb{E}[Z_{\star }^2]=V+\mathbb{E}[{\widehat{Z_{\star }}}^2]$. Using these facts and again using the independence property, we take the expectation of ${\epsilon }^2$ and obtain $\begin{array}{rl}\mathbb{E}[{\epsilon }^2]& =\sum _{j,k=1}^{\mathrm{\infty }}\frac{J_j{J}_k}{\sqrt{j!}\sqrt{k!}}{\mathrm{\partial }}_s^j{\mathrm{\partial }}_t^k(\exp \{\mathbb{E}[\widehat{Z_{\star }}](s+t)\\ & +\mathbb{E}[{\widehat{Z_{\star }}}^2]{(s+t)}^2/2-(s^2+t^2)/2\}\\ & \cdot (\exp \{V(s+t{)}^2/2\}-\exp \{V(s^2+t^2)/2\})){|}_{s,t=0}\\ & =\sum _{j,k=1}^{\mathrm{\infty }}\frac{J_j{J}_k}{\sqrt{j!}\sqrt{k!}}{\mathrm{\partial }}_s^j{\mathrm{\partial }}_t^k\\ & \times (\exp \{(1-V)(s+t{)}^2/2-(s^2+t^2)/2\}\\ & \cdot (\exp \{V(s+t{)}^2/2\}-\exp \{V(s^2+t^2)/2\})){|}_{s,t=0}\\ & =\sum _{j,k=1}^{\mathrm{\infty }}\frac{J_j{J}_k}{\sqrt{j!}\sqrt{k!}}{\mathrm{\partial }}_s^j{\mathrm{\partial }}_t^k(\exp \{st\}\cdot (1-\exp \{-Vst\})){|}_{s,t=0}.\end{array}$Now, it is straight-forward to show that for any constant A, $\begin{array}{r}{\mathrm{\partial }}_s^j{\mathrm{\partial }}_t^k\exp \{Ast\}{|}_{s,t=0}=\{\begin{array}{ll}{A}^{k}k!& \text{if}j=k\\ 0& \text{else}.\end{array}\end{array}$Applying this with A = 1 and A = −V yields $\begin{array}{rl}\mathbb{E}[{\epsilon }^2]& =\sum _{k=1}^{\mathrm{\infty }}\frac{J_k^2}{k!}(1-{(1-V)}^k)k!\\ & =\sum _{k=1}^{\mathrm{\infty }}{J}_k^2(1-{(1-V)}^k).\end{array}$