A Comparison of Prediction Algorithms and Nexting for Short Term Weather Forecasts

TL;DR

This paper compares various supervised learning algorithms, including Nexting, neural networks, and regression trees, for short-term weather forecasting, highlighting that the choice of algorithm depends on data characteristics and side information.

Contribution

It evaluates the performance of different regression algorithms on weather data, demonstrating that Nexting performs well with slowly varying signals and sufficient training data.

Findings

Nexting performs well with slowly varying signals.

Algorithm choice depends on available side information.

No single method is clearly superior across all cases.

Abstract

This report first provides a brief overview of a number of supervised learning algorithms for regression tasks. Among those are neural networks, regression trees, and the recently introduced Nexting. Nexting has been presented in the context of reinforcement learning where it was used to predict a large number of signals at different timescales. In the second half of this report, we apply the algorithms to historical weather data in order to evaluate their suitability to forecast a local weather trend. Our experiments did not identify one clearly preferable method, but rather show that choosing an appropriate algorithm depends on the available side information. For slowly varying signals and a proficient number of training samples, Nexting achieved good results in the studied cases.

Tables

Table 1. Table I: Wind speed, dry-bulb temperature, and direct normal irradiance approximated by different methods trained with 24 samples. The training root mean square error is given as well as the number of consecutive forecast samples that lie within a interval around the target signal.

Method	Wind Speed RMSE	Wind Speed $\pm 1,\pm 3\frac{m}{s}$	Temperature RMSE	Temperature $\pm 0.5,\pm {1.5}^{\circ }C$	Irradiance RMSE	Irradiance $\pm 100,\pm 300\frac{\mathrm{Wh}}{m^2}$
@a xhline Polynomial Regression	0.9337	2, 7	0.3467	2, 2	72.6259	1, 1
Ridge Regression	1.0306	2, 24	0.7389	6, 24	173.5266	5, 24
Radial Basis Function Network	0.9358	5, 8	0.4149	7, 9	44.6131	8, 9
Smoothing Splines	0.9286	3, 5	0.1441	3, 3	25.8115	9, 9
Seasonal ARIMA1	-	2, 15	-	3, 23	-	8, 9
Regression Tree	1.2096	2, 7	0.6020	0, 17	132.4193	8, 16
Nexting2	2.6362	1, 2	1.4651	2, 8	317.3383	5, 5
1Two training periods were used.
2The weight vector was fixed after one period. To compute the RMSE, the signal was rescaled and shifted in time.

Table 2. Table II: Settings for the Results in Table  I

Method	Parameters	Wind Speed	Temperature	Irradiance
Polynomial Regression	polynomial degree $L$	6	7	6
Ridge Regression	regularization parameter	0.1	0.1	5
	function $g_1$	$\cos \left(\frac{2\pi }{24}x-\frac{\pi }{4}\right)$	$\cos \left(\frac{2\pi }{24}x\right)$	$\cos \left(\frac{2\pi }{24}x\right)$
Radial Basis Function	number of basis functions $N$	4	5	4
	scaling parameter $\sigma$	6.3	6.6	3.7
Smoothing Splines	smoothing parameter $\lambda$	9	0.1	0.5
Seasonal ARIMA	AR orders $p,P$	0, 1	2, 0	0, 1
	MA order $q,Q$	3, 0	0, 0	1, 1
	seasonality $s$	24	24	24
	degrees of differencing $d,D$	0, 1	2, 1	0, 1
Regression Tree	number of nodes	5	5	5
Nexting	discount rate ${\gamma }^{\left(i\right)}$	0	0	0
	step-size parameter $\alpha$	0.3	0.3	0.2
	trace-decay parameter $\lambda$	0.9	0.9	0.9

Equations

$$f_{\mathrm{PR}}(x)=\sum _{l=0}^L{\theta }_l{x}^l={\theta }_0+{\theta }_{1}x+{\theta }_2{x}^2+\cdots +{\theta }_L{x}^L$$

$$\underset{\text{mathbit}y}{\underbrace{\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(D)}\end{array}\right]}}=\underset{\text{mathbit}X}{\underbrace{\left[\begin{array}{ccccc}1& {\xi }^{(1)}& {\xi }^{(1),2}& \dots & {\xi }^{(1),L}\\ 1& {\xi }^{(2)}& {\xi }^{(2),2}& \dots & {\xi }^{(2),L}\\ ⋮& ⋮& ⋮& & ⋮\\ 1& {\xi }^{(D)}& {\xi }^{(D),2}& \dots & {\xi }^{(D),L}\end{array}\right]}}\underset{\text{mathbit}\theta }{\underbrace{\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_L\end{array}\right]}}.$$

$${\widehat{\text{mathbit}\theta }}_{\mathrm{PR}}=\arg \underset{\text{mathbit}\theta }{\min }\Vert \text{mathbit}y-\text{mathbit}X\text{mathbit}\theta {\Vert }_2^2$$

$${\widehat{\text{mathbit}\theta }}_{\mathrm{PR}}=(\text{mathbit}{X}^T\text{mathbit}X{)}^{-1}\text{mathbit}{X}^T\text{mathbit}y$$

$${\hat{f}}_{\mathrm{PR}}(x)=\sum _{l=0}^L{\hat{\theta }}_{l,\mathrm{PR}}{x}^l.$$

$${\widehat{\text{mathbit}\theta }}_{\mathrm{RR}}=\arg \underset{\text{mathbit}\theta }{\min }\Vert \text{mathbit}y-\text{mathbit}X\text{mathbit}\theta {\Vert }_2^2+{\lambda }^2\Vert \text{mathbit}L\text{mathbit}\theta {\Vert }_2^2$$

$${\widehat{\text{mathbit}\theta }}_{\mathrm{RR}}=(\text{mathbit}{X}^T\text{mathbit}X+{\lambda }^2\text{mathbit}{L}^T\text{mathbit}L{)}^{-1}\text{mathbit}{X}^T\text{mathbit}y$$

$$\text{mathbit}X=\left[\begin{array}{cccc}{g}_0({\xi }^{(1)})& g_1({\xi }^{(1)})& \dots & g_K({\xi }^{(1)})\\ g_0({\xi }^{(2)})& g_1({\xi }^{(2)})& \dots & g_K({\xi }^{(2)})\\ ⋮& ⋮& & ⋮\\ g_0({\xi }^{(D)})& g_1({\xi }^{(D)})& \dots & g_K({\xi }^{(D)})\end{array}\right].$$

$${\hat{y}}^{(i)}={\hat{f}}_{\mathrm{RR}}({\xi }^{(i)})=\sum _{k=0}^K{\hat{\theta }}_{k,\mathrm{RR}}{g}_k({\xi }^{(i)}).$$

$$f_{\mathrm{RB}}(\text{mathbit}x)=b_0+\sum _{n=1}^N{\theta }_{n,\mathrm{RB}}{g}_n(\Vert \text{mathbit}x-\text{mathbit}{c}_i{\Vert }_2^2)$$

$$\sum _{d=1}^D(y^{(i)}-f_{\mathrm{SS}}({\xi }^{(i)}){)}^2+\lambda \int (f_{\mathrm{SS}}^{\prime \prime }(t){)}^2\mathrm{d}t$$

$$f_{\mathrm{SS}}(x)=\sum _{d=1}^D{\theta }_d{N}_d(x).$$

$${\Delta }_d(x)=\frac{(x-{\xi }^{(d)}{)}_{+}^3-(x-{\xi }^{(D)}{)}_{+}^3}{{\xi }^{(D)}-{\xi }^{(d)}}$$

$${\widehat{\text{mathbit}\theta }}_{\mathrm{SS}}=\arg \underset{\text{mathbit}\theta }{\min }(\text{mathbit}y-\text{mathbit}X\text{mathbit}\theta {)}^T(\text{mathbit}y-\text{mathbit}X\text{mathbit}\theta )+\lambda \text{mathbit}{\theta }^T\text{mathbit}\Omega \text{mathbit}\theta .$$

$${\hat{y}}^{(i)}={\hat{f}}_{\mathrm{SS}}({\xi }^{(i)})=\sum _{d=1}^D{\hat{\theta }}_{d,\mathrm{SS}}{N}_d({\xi }^{(i)}).$$

$${\hat{f}}_{\mathrm{KR}}(\text{mathbit}x)=\frac{\sum _{i=1}^D{y}^{(i)}\left(\prod _{j=1}^P{K}_{\lambda }(x_j,{\xi }_j^{(i)})\right)}{\sum _{j=1}^D\left(\prod _{j=1}^P{K}_{\lambda }(x_j,{\xi }_j^{(i)})\right)}$$

$$K_{\lambda }(x,{\xi }_i)\ge 0$$

$$\int x{K}_{\lambda }(0,x)\mathrm{d}x=0$$

$$0<\int x^2{K}_{\lambda }(0,x)\mathrm{d}x<\infty .$$

$$K_{\lambda }(x,{\xi }_i)=\frac{1}{\sqrt{2\pi }}\exp \left(\frac{-(x-{\xi }_i{)}^2}{2{\lambda }^2}\right).$$

$$\tilde{y}[t]={\varphi }_1\tilde{y}[t-1]+{\varphi }_2\tilde{y}[t-2]+\cdots +{\varphi }_p\tilde{y}[t-p]+a[t]$$

$$\tilde{y}[t]=a[t]-{\theta }_{1}a[t-1]-{\theta }_{2}a[t-2]-\cdots -{\theta }_{q}a[t-q].$$

$$\phi (B)y[t]={\varphi }_p(B){\nabla }^{d}y[t]={\theta }_0+{\theta }_q(B)a[t]$$

$$\begin{gathered} y[t]={\phi }_{1}y[t-1]+\cdots +{\phi }_{p+d}y[t-p-d] \\ -{\theta }_{1}a[t-1]-\cdots -{\theta }_{q}a[t-q]+a[t]. \end{gathered}$$

$$\hat{y}[t+l]=E[y[t+l]|y[t],y[t-1],\dots ]=E_t[y[t+l]]$$

$$\begin{gathered} E_t[y[t+l]]={\phi }_1{E}_t[y[t+l-1]]+\cdots +{\phi }_{p+d}{E}_t[y[t+l-p-d]] \\ -{\theta }_1{E}_t[a[t+l-1]]-\cdots -{\theta }_q{E}_t[a[t+l-q]]+E_t[a[t+l]]. \end{gathered}$$

$${\varphi }_p(B){\Phi }_P(B^s){\nabla }^d{\nabla }_s^{D}y[t]={\theta }_q(B){\Theta }_Q(B^s)a_t$$

$${\mathbb{I}}_{{\mathcal{R}}_m}(\text{mathbit}x)=\{\begin{array}{l}1\text{if }\text{mathbit}x\in {\mathcal{R}}_m\\ 0\text{otherwise}\end{array}$$

$$f_{\mathrm{RT}}(\text{mathbit}x)=\sum _{m=1}^M{c}_m{\mathbb{I}}_{{\mathcal{R}}_m}(\text{mathbit}x)$$

$$\underset{j,s}{\min }(\underset{c_1}{\min }\sum _{i:\text{mathbit}{\xi }^{(i)}\in {\mathcal{R}}_1(j,s)}(y^{(i)}-c_1{)}^2+\underset{c_2}{\min }\sum _{i:\text{mathbit}{\xi }^{(i)}\in {\mathcal{R}}_2(j,s)}(y^{(i)}-c_2{)}^2).$$

$${\hat{c}}_m$$

Full text

I Introduction

Typically, regression algorithms follow a two stage process. In a first step, the algorithms are trained on historic data and optimized. Then, they are used to predict outcomes based on the learned rules. This approach is known as supervised learning.

In this report, we compare a number of such algorithms. First, we briefly introduce them in Section II and illustrate their main idea by applying them to a sinusoidal signal. One of the chosen algorithms is the recently introduced Nexting [1]. It was designed to predict numerous target signals at multiple timescales. We compare Nexting to well-known methods like polynomial regression.

In Section III, we apply the most promising candidates to historical weather measurements. The results are compared to assess the suitability of the algorithms as a predictor of future data. Being able to predict a local trend might be interesting for applications that depend on weather conditions. For example, this could help to improve forecasts of the energy production of solar systems or wind power plants. The last section then concludes our findings.

II Algorithms

We write vectors like $\mathbit{x}\in\mathbb{R}^{P}$ in bold face and denote their elements $x_{i}$ by subscripts.

Consider a physical system whose output $y\in\mathbb{R}$ can be described by a function $f(\mathbit{x})=y$ where $\mathbit{x}\in\mathbb{R}^{P}$ collects the $P$ system inputs. Assuming that the function $f(\mathbit{x})$ is unknown, we try to approximate it by some $\hat{f}(\mathbit{x})$. To determine this approximation, we first take $D+F$ measurements $(\mathbit{\xi}^{(i)},y^{(i)})\in\mathbb{R}^{P}\times\mathbb{R}$ with $y^{(i)}=f(\mathbit{\xi}^{(i)})$, $i=1,2,\dots,D,D+1,\dots,D+F$. Then, we apply the regression methods described below to the data set $\{(\mathbit{\xi}^{(i)},y^{(i)})\}_{i=1}^{D}$, which yields $\hat{f}(\mathbit{x})$. Having done so, we are interested in the behavior of $\hat{f}(\mathbit{x})$ for input $\mathbit{\xi}^{(i)}$ with $i>D$. We thus aim at assessing the suitability of $\hat{f}(\mathbit{x})$ as a predictor of future data $y^{(D+1)},y^{(D+2)},\dots,y^{(D+F)}$.

In the special case that $P=1$ and $f(x)=f(t)$ is a function of time $t$ which we sample with period $T$, we write $y[t+j]=f(t+jT)$ for integers $j$. This is a convenient notation for some of the following algorithms and in case of periodic signals.

II-A Polynomial Regression

In polynomial regression, we model the unknown function $f(x)$ by a polynomial

[TABLE]

of degree $L$. Note that $x\in\mathbb{R}$ in this case. Given our training data $\{(\xi^{(i)},y^{(i)})\}_{i=1}^{D}$, we have the system of equations

[TABLE]

To fit the polynomial coefficients $\mathbit{\theta}$, we use the least squares approach, i.e., we solve the optimization problem

[TABLE]

where $\|.\|_{2}$ denotes the Euclidean norm. Writing the norm as an inner product, setting the derivative with respect to $\mathbit{\theta}$ to zero, and solving for $\mathbit{\theta}$ leads to the unique solution

[TABLE]

provided that the inverse $(\mathbit{X}^{\operatorname{T}}\mathbit{X})^{-1}$ exists (cf. [2]).

Summing up the above, our approximating function is

[TABLE]

Predicting the $y^{(i)}$ is therefore done by computing $\hat{y}^{(i)}=\hat{f}_{\mathrm{PR}}(\xi^{(i)})$. More information on polynomial regression can be found in [2] and [3].

As an example, Fig. 1 shows the approximation of the test function $f(x)=\sin(\frac{2\pi}{100}x)$ (solid red) by a polynomial of degree $L=3$ (blue dots). To compute the polynomial coefficients, $D=100$ data points $\{(i,\sin(i))\}_{i=1}^{D}$ were used, i.e., in this case, we have $\xi^{(i)}=i$ and $y^{(i)}=f(\xi^{(i)})=\sin(\frac{2\pi}{100}i)$. While the approximation seems to be acceptable for $x\in[1,100]$, it is poor for $x>100$. The author of [4] addresses further problems of polynomial regression.

II-B Ridge Regression

In statistics, ridge regression is a synonym for Tikhonov regularization [5]. Often, this method is applied if solving a system of equations like $\mathbit{X}\mathbit{\theta}=\mathbit{y}$ (cf. (2)) is not well posed. In the general case, a regularization term is introduced to form the optimization problem

[TABLE]

where $\lambda$ is the regularization parameter (cf. [5]). The matrix $\mathbit{L}$ allows to represent prior assumptions on the size of the solution by the norm $\|\mathbit{L}\mathbit{\theta}\|_{2}$. In most cases, $\mathbit{L}=\mathbf{I}$ is a good choice.

An alternative formulation of the problem is to minimize $(\mathbit{X}^{\operatorname{T}}\mathbit{X}+\lambda^{2}\mathbit{L}^{\operatorname{T}}\mathbit{L})\mathbit{\theta}=\mathbit{X}^{\operatorname{T}}\mathbit{y}$ (cf. [5]). If the intersection of the null-spaces of $\mathbit{X}$ and $\mathbit{L}$ contains $\boldsymbol{0}$ only, the unique solution is given by

[TABLE]

which is derived in [5].

For our purposes, ridge regression fits the weights $\theta_{k}$ of the function $f_{\mathrm{RR}}(x)=\sum_{k=0}^{K}\theta_{k}g_{k}(x)$ by means of the minimization in (6). The $g_{k}(x)$ denote arbitrary functions. So, e.g., they can be sinusoidal or exponential functions. The structure of $\mathbit{X}$ therefore looks as follows:

[TABLE]

Note that for $\lambda=0$, $K=L$, and $g_{k}(x)=x^{k}$ equation (7) corresponds to (4), the least squares polynomial regression.

We predict data with

[TABLE]

More information on ridge regression can be found in [5] and [6].

As the functions $g_{k}$ can be chosen arbitrarily, ridge regression allows to incorporate some knowledge of the input signal into the approximation. Considering that, the test function $f(x)=\sin(\frac{2\pi}{100}x)$ is an interesting task. If we know the input function can be modeled by a sinusoidal signal, we can first compute the frequency from the given data and then use an according sine function as $g_{0}(x)$ in (8). In this special case, the choice $K=1,$ $\lambda=0$ leads to no deviation between the approximation and $f(x)$.

II-C Radial Basis Function Network

According to [7], multilayer feedforward networks with only one hidden layer can approximate any measurable function with arbitrary accuracy. This makes feedforward networks an interesting tool for approximation.

The radial basis function (RBF) network is a single hidden layer feedforward network (cf. [8]). The output transfer functions are linear, the hidden-layer transfer functions $g_{n}(z)$ are nonlinear. Typically, we have Gaussian $g_{n}(z)=\exp(-\frac{z}{2\sigma^{2}})$. The network output is then given by

[TABLE]

where $b_{0}$ is a bias term, the $\mathbit{c}_{i}$ are called centers, and $N$ denotes the number of basis functions (hidden neurons) (cf. [8]). The bias is sometimes incorporated into the summation by introducing an extra constant basis function $g_{0}(z)=1$.

There are different training schemes for RBF networks. A simple way is to fix the $\mathbit{c}_{i}$ in advance. For example, this can be done by evenly spacing the centers over the input space or choosing the centers equal to the input vectors. If the scaling parameters $\sigma$ are also fixed, only the weights $\theta_{n,\mathrm{RB}}$ have to be determined. This can be done by minimizing the squared error between the network output $f_{\mathrm{RB}}(\mathbit{\xi}^{(i)})$ and the desired $y^{(i)}$ (cf. [9]). More information on RBF networks can be found in [8, 9, 10].

Fig. 1 illustrates the principle of RBF networks (green triangles). We chose $N=2$ Gaussian basis functions with $\sigma=10$, which can be recognized in the graph. By increasing $N$, the approximation can be improved. As argued in [11], RBF networks are capable of universal approximation.

II-D Smoothing Splines

The spline basis approach aims at finding a twice continuously differentiable function $\hat{f}_{\mathrm{SS}}(x)$ that minimizes

[TABLE]

(cf. [2]). Usually, the smoothing parameter $\lambda$ is positive. With $\lambda=0$, the minimizer can be any interpolating function.

The unique minimizer of (11) is known to be a natural cubic spline (cf. [2]). A natural cubic spline $f_{\mathrm{SS}}(x)$ with $D$ knots can be represented by $D$ basis functions $N_{d}(x)$ as

[TABLE]

The basis functions are given by $N_{1}(x)=1$, $N_{2}(x)=x$, $N_{d+2}(x)=\Delta_{d}(x)-\Delta_{D-1}(x)$ with

[TABLE]

where the data points $\xi^{(d)}$ are knots and $(.)_{+}$ denotes the positive part (cf. [2]).

To find the optimal weights in (12), we define the matrices $[\mathbit{X}]_{ij}=N_{j}(\xi^{(i)})$ and $[\mathbit{\Omega}]_{jk}=\int N_{j}^{\prime\prime}(t)N_{k}^{\prime\prime}(t)\mathrm{d}t$ for $i,j,k=1,2,\dots,D$. This allows us to formulate the optimization problem

[TABLE]

The solution is then given by $\hat{\mathbit{\theta}}_{\mathrm{SS}}=(\mathbit{X}^{\operatorname{T}}\mathbit{X}+\lambda\mathbit{\Omega})^{-1}\mathbit{X}^{\operatorname{T}}\mathbit{y}$ (cf. [2]). As for large $D$ the matrices get big, the authors of [2] suggest to apply thinning strategies to simplify the computations. That is, not all internal knots are used. The effect on the fit is negligible.

Having found $\hat{\mathbit{\theta}}_{\mathrm{SS}}$, we predict data with

[TABLE]

Fig. 2 shows the result (green triangles) of this method applied to the test function $f(x)=\sin(\frac{2\pi}{100}x)$. Increasing the smoothing parameter $\lambda$ puts more emphasis on the second derivative penalty in (11), i.e., smooths the fit.

II-E Kernel Regression

Kernel regression is also known as kernel smoothing. The basic idea was introduced by Nadaraya [12] and Watson [13]. The approximating function can be written as

[TABLE]

where $K_{\lambda}(x,\xi_{i})$ is a kernel parameterized by the smoothing parameter $\lambda$ (cf. [14]). Each kernel needs to meet three properties:

[TABLE]

A common choice is the Gaussian kernel

[TABLE]

More kernel functions can be found in [14]. However, choosing an appropriate kernel is not as important as properly selecting the bandwidth $\lambda>0$. With increasing $\lambda$ the approximation smooths. Setting the bandwidth equal to the variance of all $y^{(i)}$ is a simple approach (see [14]).

Typically, $\hat{f}_{\mathrm{KR}}(x)$ shows bad behavior on the boundaries. This can also be observed in Fig. 2 (blue dots) for $x$ near 0 and 100. To overcome this issue, local polynomial fits at the borders have been suggested (see, e.g., [2]). For the signal in Fig. 2, we selected a Gaussian kernel and bandwidth $\lambda=4$.

II-F Autoregressive Integrated Moving Average

The autoregressive integrated moving average (ARIMA) model is a linear non-stationary model. In this context, non-stationarity means that the data does not move around a fixed mean. The idea of an integrated model is to reduce this non-stationary behavior by a differentiation step. After this, an autoregressive moving average (ARMA) model is applied to the differentiated data (cf. [15]).

ARIMA is often used in time series analysis. For this reason, we treat our data as a time series $y[t]$ with time index $t$. We define two operators to simplify the notation (similar to [15]). The backward shift or lag operator $B$ shifts the time index: $B^{m}y[t]=y[t-m]$. And the backward difference operator $\nabla$ enables a short notation for the difference between two time series values: $\nabla y[t]=y[t]-y[t-1]=(1-B)y[t]$. Before we discuss the general ARIMA model, we first introduce its two components, namely autoregressive and moving average processes, separately.

The value of an autoregressive process of order $p$ at time $t$ depends on a linear combination of $p$ previous values [15]:

[TABLE]

where $a[t]$ is a random shock from a white noise process with variance $\sigma_{a}^{2}$. For convenience, the series $\tilde{y}[t]=y[t]-\mu$ of deviation from the mean $\mu$ is introduced. Using the autoregressive operator $\phi_{p}(B)=1-\phi_{1}B-\phi_{2}B^{2}-\dots-\phi_{p}B^{p}$ of order $p$, the model can be written as $\phi_{p}(B)\tilde{y}[t]=a[t]$ (cf. [15]). This model has $p+2$ unknown parameters: $\mu$, $\sigma_{a}^{2}$, and the $\phi_{i}$.

In contrast, the value of a moving average process of order $q$ at time $t$ depends on a linear combination of $q$ previous random shocks [15]:

[TABLE]

Using the moving average operator $\theta_{q}(B)=1-\theta_{1}B-\theta_{2}B^{2}-\dots-\theta_{q}B^{q}$, the model is given by $\tilde{y}[t]=\theta_{q}(B)a[t]$. It has $q+2$ unknown parameters: $\mu$, $\sigma_{a}^{2}$, and the $\theta_{i}$.

The general ARIMA($p,d,q$) model is formulated as

[TABLE]

where we use the non-stationarity operator $\varphi(B)=\phi_{p}(B)\nabla^{d}=1-\varphi_{1}B-\varphi_{2}B^{2}-\dots-\varphi_{p+d}B^{p+d}$ (cf. [15]). With $\theta_{0}\neq 0$ a deterministic polynomial trend of order $d$ can be realized. We can also express the model by means of the difference equation:

[TABLE]

To determine a forecast $\hat{y}[t+l]$ for lead time $l>0$ at origin $t$, the mean squared error $\operatorname{E}[(y[t+l]-\hat{y}[t+l])^{2}]$ is minimized. If the $a[t]$ form a sequence of independent random variables, it can be shown that the MMSE forecast is given by the conditional expectation

[TABLE]

where $\operatorname{E}_{t}[y[t+l]]$ is a short notation for the conditional expectation. We can then use the difference equation (24) to compute a forecast:

[TABLE]

The conditional expectations in (26) are computed by inserting $y$’s and $a$’s when they are given. Future $a$’s are set to zero, and intermediate values $\operatorname{E}_{t}[y[t+1]]$, $\operatorname{E}_{t}[y[t+2]]$, …, $\operatorname{E}_{t}[y[t+l-1]]$ are calculated iteratively. Using the difference equation (24) is only one way to obtain forecasts. Two more possibilities are given in [15] which may be advantageous in certain situations.

A general method for obtaining initial parameter estimates of an ARIMA process is described in [15]. In a preparation step, a new series $w[t]=\nabla^{d}y[t]$ of differences is generated, then, a three stage procedure is applied. However, an outline of the whole method is beyond the scope of this report as much context and many definitions are needed. The interested reader is referred to the chapters on model identification in [15].

In some applications, a seasonal behavior of the data can be observed. For example, if we model hourly changes of the temperature, we expect similarities between measurements that are 24 hours apart. To incorporate such a periodicity in an ARIMA model, the authors of [15] introduce the so-called multiplicative model

[TABLE]

as an extension of (23). Here, $s$ describes the periodicity. So, e.g., in the above temperature example we have $s=24$. The operator $\nabla_{s}^{D}$ is defined by $\nabla_{s}=1-B^{s}$, and $\Phi_{P}(B^{s})$ and $\Theta_{Q}(B^{s})$ are polynomials in $B^{s}$ of orders $P$ and $Q$, respectively. Forecasting can again be done by means of the according difference equation. For more information, please refer to the chapter on seasonal models in [15].

In the example with the test function $f(x)=\sin(\frac{2\pi}{100}x)$, it turns out to be sufficient to use an ARIMA(2,0,0) model to predict one period of the sine with a mean squared error of less than $10^{-10}$.

II-G Regression Tree

The following overview is based on the considerations in [2]. The idea of binary regression trees is to partition the input space into $M$ rectangular regions $\mathcal{R}_{m}$ and apply a simple fit, e.g., a constant, to each region. Using the indicator function

[TABLE]

the output of a regression tree is given by

[TABLE]

i.e., the constant $c_{m}$ of the region $\mathcal{R}_{m}$ to which the element $\mathbit{x}$ belongs. By minimizing the sum of squares $\sum_{i=1}^{D}(y^{(i)}-f_{\mathrm{RT}}(\mathbit{\xi}^{(i)}))^{2}$, the constant $c_{m}$ associated with region $\mathcal{R}_{m}$ is determined as the average of all $y^{(i)}$ in that region: $c_{m}=\mathrm{average}(y^{(i)}|\mathbit{\xi}^{(i)}\in\mathcal{R}_{m})$ (cf. [2]). Thus, computing the constants is simple, and the main challenge lies in determining the regions.

We describe a greedy algorithm from [2] to achieve binary partitions in the following. Let $x_{j}$ be a splitting variable and $s$ a split joint. We define the two half spaces $\mathcal{R}_{1}(j,s)=\{\mathbit{x}|x_{j}\leq s\}$ and $\mathcal{R}_{2}(j,s)=\{\mathbit{x}|x_{j}>s\}$. Then, $j$ and $s$ are found by solving

[TABLE]

For any $j$ and $s$, the solution to the inner problem is given by $\hat{c}_{1}=\mathrm{average}(y^{(i)}|\mathbit{\xi}^{(i)}\in R_{1}(j,s))$ and $\hat{c}_{2}=\mathrm{average}(y^{(i)}|\mathbit{\xi}^{(i)}\in R_{2}(j,s))$. Hence, computing the split joint can be done quickly if $j$ is given. By browsing through all inputs, the best pair of $j$ and $s$ is found. The result of this partitioning are two regions on which the same process is repeated.

Basically, this is done until a predetermined tree size is reached. However, it might be difficult to choose a proper tree size. One strategy is to first compute a tree $T_{0}$ with a minimum node-size and then prune it by means of cost-complexity pruning which we illustrate briefly. By pruning $T_{0}$, i.e., collapsing internal nodes, we get a subtree $T\subset T_{0}$. Let $|T|$ be the number of terminal nodes (regions $\mathcal{R}_{m}$) of $T$. Further, denote by $N_{m}$ the number of elements $\mathbit{\xi}^{(i)}$ in region $R_{m}$. Additionally, we define the following two quantities:

[TABLE]

Finally, our cost complexity criterion is defined as

[TABLE]

For each $\alpha$, we aim at finding the subtree $T_{\alpha}\subseteq T_{0}$ which minimizes the cost criterion. The resulting subtree is unique. With the parameter $\alpha\geq 0$, we control the tree size and fit quality. Large values lead to small trees, whereas $\alpha=0$ is the original tree $T_{0}$. The parameter $\alpha$ can be chosen in accordance to preferences or it can be estimated by cross-validation (cf. [2]).

One technique to reduce the variance of the prediction is called bagging (cf. [2]) Basically, the number of given training data is artificially increased by random sampling with replacement. Given a set $\mathcal{S}$ of $D$ training samples, we generate $B$ sets $\mathcal{S}_{1},\mathcal{S}_{2},\dots,\mathcal{S}_{B}$ with $D$ data points each. Then, a tree model $f_{\mathrm{RT}}^{(b)}(x)$ is fit to every set $\mathcal{S}_{b}$ and the bagging estimate $f_{\mathrm{bag}}(x)$ is defined by

[TABLE]

In Fig. 2, the regression tree approximation (orange squares) of $f(x)=\sin(\frac{2\pi}{100}x)$ is given, where $D=100$ training samples were used. Note that the regions $\mathcal{R}_{m}$ and the constant fits can be recognized. As can be seen, evaluating the tree for $x>100$ gives the constant of the rightmost region. Thus, in the given one-dimensional case a regression tree without further modifications is not a good choice to forecast data.

II-H Nexting

Basically, Nexting, as introduced in [1], treats the given data as a time series. Assume, we have $P$ signals $y^{(i)}[t]$, $i=1,2,\dots,P$, where $t$ denotes the time. The idea of Nexting is to approximate the return

[TABLE]

of signal $i$ by a function $f_{\mathrm{nexting}}^{(i)}[t]$. The discount-rate $\gamma^{(i)}$ controls the timescale. For $\gamma^{(i)}=0$, the return $G^{(i)}[t]=y^{(i)}[t+1]$ reduces to the value in the immediate next time step, which is very short-sighted. With $\gamma^{(i)}>0$ all (discounted) future signal values are taken into account.

The approximating function has the form

[TABLE]

This is the inner product of the feature vector $\mathbit{\phi}[t]$ and the weight vector $\mathbit{\theta}^{(i)}[t]$. Note that there is only one feature vector for all $P$ signals. This vector is constructed by tile coding after normalizing all signals to lie in the interval $[0,1]$ (cf. [1]). As a result, $\mathbit{\phi}[t]\in\{0,1\}^{N}$ has a constant number of 1 features.

The weight vectors $\mathbit{\theta}^{(i)}[t]$ are determined by a linear temporal-difference approach [1]:

[TABLE]

The so-called step-size parameter $\alpha$ is positive and the eligibility vector $\mathbit{e}^{(i)}[t]$ follows the update rule $\mathbit{e}^{(i)}[t]=\gamma^{(i)}\lambda\mathbit{e}^{(i)}[t-1]+\mathbit{\phi}[t]$. Here, $\lambda\in[0,1]$ is the trace-decay parameter. A special case is $\lambda=1$. For this choice, $\mathbit{\theta}^{(i)}[t]$ converges asymptotically to the vector which minimizes the mean squared error between the prediction and $G^{(i)}[t]$.

For the choice $\alpha=0.1$, $\lambda=0.9$ and the two $\gamma$-values $\gamma=0$ and $\gamma=0.9375$, Fig. 3 depicts the result of Nexting for the sinusoid signal $f(x)=\sin(\frac{2\pi}{100}x)$. As stated above, choosing $\gamma=0$ should lead to a prediction of the function value in the immediate next time step. That is, $\hat{f}_{\mathrm{nexting},\gamma=0}(x)$ is a forecast of $f(x+1)$. However, in the situation shown in Fig. 3, this behavior is hardly recognizable. Similarly, larger values of $\gamma$ should indicate changes in $f(x)$ earlier. Again, this can hardly be recognized in our example. We remark, however, that in Section III, where we apply Nexting to weather data, the expected forecasting behavior can be observed.

Note that the principle of Nexting differs from the algorithms we introduced so far. According to (36), in order to compute $\hat{f}_{\mathrm{nexting}}(x)$ for $x>100$, we use $\mathbit{\phi}[x]$ and $\mathbit{\theta}[x]$, i.e., the prediction $\hat{f}_{\mathrm{nexting}}(x)$ is calculated online. In contrast, polynomial regression, for example, gives a function $\hat{f}_{\mathrm{PR}}(x)$ which we evaluate for $x>100$ (cf. (5)). Thus, Nexting forecasts are based on more input data than, e.g., those of polynomial regression which makes comparisons difficult. Further, in the example with the sinusoid input signal of Fig. 3, Nexting requires several training periods to reach a steady state. For the plot, we used 10 training periods and thus a set of $D=1000$ data points.

III Applying the Algorithms to Weather Data

In the previous section, we had $f(x)=\sin(\frac{2\pi}{100}x)$ as an example function to illustrate the idea of the algorithms we use. In this section, we take hourly measurements of wind speed, dry-bulb temperature, and direct normal irradiance as input. The data points stem from a typical meteorological year 3 (TMY3) set provided by the national renewable energy laboratory [16] where we chose the weather station of Los Angeles.

In a first step, we apply the methods from Section II to the three sets of weather measurements in order to compare the different algorithms. Following that, we take a closer look at Nexting and regression trees.

In Fig. 2 (blue dots), we already saw the critical behavior of kernel smoothing at the boundaries. This method is thus not suited for data forecasting and therefore excluded from the following considerations.

III-A Comparison

For a qualitative statement of the forecast accuracy, we define intervals of a constant deviation around the input target signals. As an illustration, the shaded areas in Fig. 4 represent a deviation of $\pm 1\frac{\mathrm{m}}{\mathrm{s}}$ and $\pm 3\frac{\mathrm{m}}{\mathrm{s}}$ around the wind speed. We try to train all algorithms such that as many predictions as possible lie within the smaller interval. The idea is to gain an understanding of the algorithms’ performance.

We use $D=24$ data samples (one day) to train the algorithms and then forecast one day. Tab. I lists our findings. The first column of the three data sets represents the root mean square error (RMSE) between the actual values and the approximations for the 24 training samples. Column two then gives the number of consecutive forecast samples (the number of hours) which lie within the two intervals defined above. The settings which lead to the table entries can be found in Tab. II.

As an example, Fig. 4 depicts the polynomial regression and radial basis function (RBF) network approximations of the wind speed (solid red). From (10) we know that RBF networks approach a constant bias term for large $x$-values. In the depicted case, the wind speed curve shape coincides with the decaying behavior of the approximation which is the reason for the seemingly good forecast. In general, RBF networks are more fit to interpolate than to forecast data.

There are many ways to select the functions $g_{k}$ in (9) for ridge regression. We chose a constant $g_{0}$ and sinusoid $g_{1}$ for the three data sets. Fig. 5 shows the result for this simple choice on the example of the dry-bulb temperature. While the training RMSE is comparatively large, it is interesting to see that a sinusoidal function indicates the trend fairly well.

As we saw in Fig. 2, regression trees basically split up the abscissa into a number of intervals and fit a constant to these. Forecasting then gives the value of the rightmost interval which is in general not very accurate. Instead, one idea might be to view the tree as a prototype representation for one signal period. To make predictions, we then evaluate this prototype period. So, e.g., to forecast the temperature for $t=26$, we evaluate the prototype at $t=2$, which was done in Tab. I for all data sets. In other words, we perform a modulo 24 operation on $t$ prior to evaluating the tree.

Following the idea of tree bagging (see (34)), the described prototype approach can be used to incorporate multiple training periods in the predictions. We investigate the benefit of this in the next subsection.

Tab. I also presents the results of a seasonal ARIMA model applied to the three data sets wind speed, dry-bulb temperature, and direct normal irradiance. However, to train a model with seasonality $s=24$, we had to use two days of training data. Further, it is difficult to determine the orders $p,q,P,Q$ of the associated polynomials. The sample autocorrelation and partial autocorrelation functions give a first idea (cf. [15]). For our data, we combined the insights delivered by these functions with the practical guideline that orders of $P,Q>1$ are rarely necessary (cf. [15]) and compared the result of different parameter values.

We already mentioned the difficulty that arises in comparing Nexting to the other algorithms in Sec. II. In our studies, several training periods were necessary to reach a steady state. Thus, the bad performance of Nexting in Tab. I where only one period was used to train is not surprising. For the results in Tab. I, we fixed the weight vector $\mathbit{\theta}$ in (36) after 24 training samples, i.e., the algorithm’s learning process was stopped after 24 training samples. This makes a comparison to the other algorithms more fair.

III-B Nexting and Regression Trees with Multiple Training Periods

As stated above, Nexting achieves better results if applied to multiple training periods. Fig. 6 shows Nexting approximations for the three cases wind speed, dry-bulb temperature, and direct normal irradiance. In contrast to the previous experiments, we did not fix the weight vector. Further, 20 training periods as opposed to only one were provided. The parameter values were again those of Tab. II — except for direct normal irradiance where $\alpha=0.1$ was chosen. As by default the Nexting output lies in between zero and one, we scaled it to match the target signal for a better illustration.

The discount rate is $\gamma=0$ in all cases. Hence, we expect the value of $\hat{f}_{\mathrm{nexting}}(t)$ to be a forecast of the target signal at $t+1$ (see (35)). In Fig. 6, this behavior can generally be observed. However, there are some points where the approximation lags behind. An example is given by the left chart of Fig. 6 around $t=15$. The wind speed shows fast variations which seem to be a challenge for Nexting.

In the previous subsection, we suggested to compute regression trees for multiple signal training periods and average the predictions of those as in (34). In combination with the modulo 24 operation which we described above, this might improve the root mean square error.

The best choice for the wind speed as target turns out to be 6 training periods which improves the RMSE to 1.1893 as apposed to 1.2096 in Tab. I. Further, using two training periods in the case of direct normal irradiance results in a RMSE of 120.4159 as apposed to 132.4193. However, there is no improvement for the temperature as target.

IV Conclusion

The results of the last paragraph show that regression trees can benefit from multiple training periods. However, too many periods can worsen the performance and it is difficult to decide for the right number. Nevertheless, regression trees are suited to forecast a trend as Tab. I suggests.

We saw that ridge regression requires knowledge of the target signal model to properly select the functions $g_{i}$ in (8). The possibility to incorporate such knowledge can be an advantage. However, in contrast to smoothing splines and polynomial regression where basically only one parameter needs to be tuned, it might be difficult to decide for suitable functions. While the computation time for these three methods is about the same, ridge regression requires more preliminary work.

In our studies, choosing appropriate parameters for the seasonal ARIMA models was difficult. What is more, the computation time was high in comparison to the algorithms mentioned so far. Further, in contrast to all methods of Sec. II, the ARIMA model does not provide predictions for the training interval. Nevertheless, it is worthwhile investigating such models in future research as the forecasts were quite accurate.

We already mentioned that radial basis function networks can approximate functions with arbitrary accuracy. However, as argued above, we deem those not suited for forecasting tasks. Additionally, the training phase was time-consuming in our experiments.

Lastly, Nexting’s online learning strategy is advantageous in our weather data context as new incoming measurements can easily be incorporated to update the forecasts. However, the output needs to be scaled in order to provide useful numerical values. Also, in contrast to the other algorithms, Nexting predictions cannot be calculated for arbitrary points in time without further effort.