Feedforward Networks for Time Series

"They slide a window over the past, hand me the last twenty numbers, and ask what comes next. I never see number twenty-one, only its corpse in the loss. I have no memory between windows; every batch, I am born again, briefly omniscient about exactly twenty lags, and blind to everything before them."

A Multilayer Perceptron Staring at a Fixed Window

Part II taught you to specify a model and then estimate its parameters: an AR($p$) was a linear recurrence over lags (Chapter 5), the Kalman filter was a recurrence over a hidden state (Chapter 7), and in every case you wrote down the functional form by hand and read its coefficients off the data. Part III inverts that contract. Beginning here, the functional form itself is learned: we fix only a flexible architecture and let optimization discover what mapping from past to future it should implement. The gentlest entry point is the feedforward network, because it requires no new temporal machinery at all. We borrow the sliding window of Section 2.5, which already turned a series into a matrix of overlapping lag vectors, and we feed each lag vector to a plain multilayer perceptron. Everything you know about feedforward networks from Appendix C applies unchanged; the only temporal content is in how the data is cut and how the train/test split is drawn. We will use $L$ for the window length, $\mathbf{x}_t = (y_{t-1}, \dots, y_{t-L})$ for the input window, and $\hat{y}_t$ for the one-step forecast, notation collected in Appendix A and shared with every forecaster in the book.

By the end of this section you should be able to do five concrete things: turn any univariate series into a supervised tensor $(\mathbf{X} \in \mathbb{R}^{N \times L}, \mathbf{y} \in \mathbb{R}^{N})$ by windowing without leaking the future; write the multilayer-perceptron forecaster as an explicit nonlinear autoregression and say what it can represent that a linear AR cannot; state what the universal-approximation theorem does and does not promise for time series; set up a training loop with a leakage-safe split, train-only normalization, and early stopping; and implement the whole forecaster twice, once from scratch in numpy with a hand-derived gradient and once in PyTorch, and measure the line-count reduction. Those five skills are the foundation every later chapter in Part III builds on.

1. From Linear AR to Nonlinear Function Approximation Beginner

Start from the model you already trust. The autoregressive model of order $L$ from Chapter 5 predicts the next value as an affine function of the last $L$ observations:

This is linear regression on a lag vector, nothing more. Its strength is that the coefficients $\boldsymbol{\phi}$ are interpretable and estimable in closed form; its weakness is that the response to the past is forced to be a fixed weighted sum, the same blend of lags no matter what regime the series is in. A multilayer perceptron keeps the input $\mathbf{x}_t$ identical, the very same lag window, but replaces the single affine map with a composition of affine maps and nonlinearities:

where each $\mathbf{W}^{(k)}$ is a weight matrix, each $\mathbf{b}^{(k)}$ a bias, and $\sigma$ a pointwise nonlinearity such as $\tanh$ or the rectified linear unit $\mathrm{ReLU}(z) = \max(0, z)$. Set every $\sigma$ to the identity and collapse the matrices and you recover the AR model exactly, so the linear autoregression is the degenerate special case of the network with no nonlinearity. The whole content of "neural forecasting" at this first level is that the activation $\sigma$ is not the identity. This is why we call the multilayer perceptron over a lag window a nonlinear autoregression, often written NAR: it predicts the next value from past values, like AR, but through a learned nonlinear function rather than a fixed linear one.

The windowing step that produces $\mathbf{x}_t$ is the only piece of temporal data engineering involved, and it is exactly the operation of Section 2.5: slide a window of width $L$ along the series, and for each landing position record the $L$ values inside the window as one input row and the value immediately after the window as the target. A series of length $T$ yields $N = T - L$ supervised pairs stacked into a design matrix $\mathbf{X} \in \mathbb{R}^{N \times L}$ and a target vector $\mathbf{y} \in \mathbb{R}^{N}$. From the network's point of view there is no time axis at all anymore; there is a tabular regression problem with $L$ features, and the fact that those features happen to be consecutive lags is invisible to the multilayer perceptron. That invisibility is both the convenience that makes this section easy and the limitation that subsection three will make precise. Figure 9.1.3 shows the windowing operation that manufactures the tensor.

Why does the nonlinearity matter at all for a series? Because real temporal responses are rarely a fixed linear blend of lags. Consider a sensor whose reading saturates: below some level it rises linearly, above it the response flattens. A linear AR must pick a single slope and will overshoot in the saturated regime; a network can bend the response, holding the high-lag region flat while keeping the low-lag region steep. Or consider a series whose short-term momentum reverses sign depending on its recent average: the effect of $y_{t-1}$ on $\hat{y}_t$ is positive in one regime and negative in another. A linear model assigns $y_{t-1}$ one coefficient for all time and cannot represent this; a network can, because the hidden units let the influence of $y_{t-1}$ depend on the values of the other lags. That dependence between lags, the influence of one lag conditioned on another, is exactly the interaction structure the next subsection argues is the central thing nonlinearity buys.

It is worth tabulating the two models against each other, because reading the columns side by side makes both the shared skeleton and the single point of divergence unmistakable. The input, the target, and the data engineering are identical; the only row that differs is the map applied inside the window, and every downstream consequence flows from that one row.

2. A Minimal Neural Forecaster Beginner

The architecture is deliberately plain: an input of $L$ lags, one or two hidden layers of modest width with a nonlinear activation, and a single output unit that emits the one-step forecast $\hat{y}_t$. With a window of $L = 20$ lags, a hidden layer of $H = 32$ units, and a $\tanh$ or ReLU activation, the network has $L \cdot H + H$ weights into the hidden layer and $H + 1$ into the output, a few hundred parameters, which is enough to capture meaningful nonlinearity on a single series without overfitting wildly. The activation choice is not cosmetic. The $\tanh$ saturates smoothly and was the historical default; the ReLU is piecewise linear, trains faster, and makes the network a continuous piecewise-linear function of its input, which is often a good match for series with regime breaks. We will use $\tanh$ in the from-scratch derivation because its derivative is clean, and note where ReLU would change the picture.

What does a hidden unit actually do for a time series? Each hidden unit computes one affine combination of the lags, $z_j = \mathbf{w}_j^{\top}\mathbf{x}_t + b_j$, and passes it through $\sigma$. You can read $\mathbf{w}_j$ as a learned feature detector over the window: one unit might respond to an upward trend (positive weights on recent lags, negative on older ones), another to a local spike, another to the window's average level. The output layer then combines these learned features. This is the sense in which the network "learns its own features": where the AR model was handed the raw lags and could only weight them, the network synthesizes intermediate features from the lags and weights those. The interactions the previous subsection promised live here. Because a feature like "recent average level" enters the output through a nonlinear unit, the marginal effect of any single lag depends on the others, which is precisely what a linear model forbids.

So far the network emits a single number, the one-step-ahead forecast $\hat{y}_t$. Multi-step forecasting, predicting $h > 1$ steps ahead, splits into the same two strategies you met for classical models in Chapter 5. The recursive strategy trains a one-step network and then feeds its own prediction back in as the newest lag to predict step two, and so on, reusing one model $h$ times. It is simple and uses all the data for one task, but it compounds errors: a small mistake at step one corrupts the input for step two, and the error can snowball over the horizon. The direct strategy trains a separate network (or a single network with $h$ output units) to predict each horizon directly from the same window, so the $h$-step forecaster never consumes its own output. Direct forecasting avoids error accumulation but ignores the dependency structure across horizons and needs more models or more output capacity. The trade is identical to the classical one, which is the point: the windowing and multi-step machinery of Section 2.5 and Chapter 5 carries over to neural forecasters unchanged; only the function inside the window became nonlinear.

3. Universal Approximation and Its Temporal Catch Intermediate

There is a famous theorem that seems to promise the multilayer perceptron can do anything. The universal-approximation theorem (Cybenko 1989, Hornik 1991) states that a feedforward network with a single hidden layer and a suitable nonlinear activation can approximate any continuous function on a compact set to arbitrary accuracy, given enough hidden units. Applied to our setting, it says: for any continuous target mapping $g^\star: \mathbb{R}^{L} \to \mathbb{R}$ from a window of $L$ lags to the next value, there exists a multilayer perceptron $f_\theta$ that matches $g^\star$ as closely as you like on any bounded region of window space. The forecaster is, in this precise sense, a universal nonlinear autoregression: whatever the true one-step map from the last $L$ values to the next is, a wide enough network can represent it.

Read the theorem carefully and three caveats jump out, and for time series the third is decisive. First, it is an existence result: it guarantees a good network exists but says nothing about whether gradient descent will find it, or how much data you need. Second, "enough hidden units" can mean exponentially many; the theorem promises representability, not efficient representability. Third, and this is the catch that defines the rest of Part III, the function being approximated has a fixed, finite input of exactly $L$ lags. The universal-approximation theorem says nothing about temporal structure because there is no time in its statement; it is a theorem about functions on $\mathbb{R}^{L}$, and the multilayer perceptron inherits exactly that limitation.

Make the limitation concrete. The forecaster has a fixed, finite receptive field: it can only see $L$ lags, so any dependence on something that happened $L + 1$ steps ago is invisible, no matter how wide or deep the network. A yearly seasonal pattern needs $L \geq 365$ daily lags even to be representable, and a window that long makes the input huge and the windowing wasteful. The forecaster has no notion of order beyond the window: as the fun note observed, permuting the lag columns consistently gives an equivalent model, so the network does not know "yesterday" comes after "the day before" except insofar as the loss happens to reward weights that act as if it does. And the forecaster has no memory between windows: each prediction is computed from its own window in isolation, with nothing carried forward, so a state that built up over a long stretch of history (a slowly accumulating trend, a regime that has persisted for months) cannot be summarized and reused; it must be re-read from the window every time, and if it does not fit in $L$ lags it is simply lost.

Figure 9.1.5 draws the three gaps and pairs each with the architecture that closes it, so the structure of the entire part is visible at a glance before you meet any of the models individually.

These three gaps are not flaws to be patched within the feedforward family; they are the design requirements that motivate the rest of Part III. Recurrence (Chapter 10) adds a carried hidden state so the model has memory between steps and an unbounded effective receptive field, the exact recurrence the Kalman filter foreshadowed in Chapter 7. Convolution (Chapter 11) keeps the feedforward, parallel nature of the multilayer perceptron but shares weights across time and stacks dilations to grow the receptive field cheaply while respecting order. Attention (Chapter 12) abandons the fixed window entirely, letting each prediction look at any past position directly. Seen this way, the humble multilayer perceptron is the perfect first model of Part III: it is the baseline whose specific failures name the three great architectural ideas that follow. Route forward through the table of contents and you will find each gap closed in turn.

4. Training Essentials Without Leaking the Future Intermediate

A neural forecaster is trained by minimizing a loss over the supervised pairs, but the temporal origin of the data imposes constraints that ordinary tabular training ignores at its peril. Start with the loss. For point forecasting the default is mean squared error, $\mathcal{L}_{\text{MSE}} = \frac{1}{N}\sum_t (\hat{y}_t - y_t)^2$, which corresponds to predicting the conditional mean and penalizes large errors heavily. When the series has outliers or heavy tails (financial returns, sensor glitches) the mean absolute error $\mathcal{L}_{\text{MAE}} = \frac{1}{N}\sum_t |\hat{y}_t - y_t|$ is more robust and targets the conditional median. When you care about a specific quantile, say the 90th percentile for a capacity-planning forecast, the pinball or quantile loss $\rho_\tau(u) = \max(\tau u, (\tau - 1)u)$ with $u = y_t - \hat{y}_t$ trains the network to emit that quantile directly, the seed of the probabilistic forecasting developed fully in Chapter 19. The architecture does not change; only the loss does.

Now the constraint that makes temporal training different from tabular training: the split must respect time. In ordinary supervised learning you shuffle the rows and split randomly. Do that here and you leak the future into the past. Because windows overlap, a randomly chosen test window can contain values that also appear in a training window, so the model has effectively seen its test answers, and a random split also lets the model train on data from after the test period, which is impossible at deployment. The fix, established in Chapter 2 and routed through the table of contents, is a chronological split: the first stretch of the series is training, a later contiguous stretch is validation, and the final stretch is test, with no shuffling across the boundaries. Within training you may still shuffle the windows for minibatching, because each window's target sits at the window's own right edge and a shuffled minibatch of training windows never reaches past the training cutoff. Shuffling windows for batching is fine; shuffling the train/test boundary is leakage.

Two more essentials complete the recipe. Normalization is fit on the training segment only, as the warning insists: compute the mean $\mu$ and standard deviation $s$ of the training values, standardize every window with those frozen statistics, $\tilde{\mathbf{x}} = (\mathbf{x} - \mu)/s$, and invert the transform on the network's output to read forecasts back in original units. This keeps the network's inputs in a well-conditioned range without ever consulting the future. Early stopping uses the chronological validation segment as a clock: train while validation loss falls, keep the weights from the best validation epoch, and stop when validation loss has not improved for a patience window of epochs. Because the validation segment comes after training and before test in time, a model that early-stops well on validation is one that generalizes forward, which is the only kind of generalization a forecaster is allowed to claim. Minibatching the training windows (typically 32 to 256 windows per batch) makes each gradient step cheap and adds a little stochastic regularization, and the Adam optimizer is the standard default for these small networks.

A subtlety worth naming explicitly is why minibatching is safe even though the data is a time series. The instinct from Chapter 2 is that order matters, so surely shuffling is forbidden. But the order that matters is the order inside each window and the order of the train, validation, and test segments, not the order in which the already-formed training windows are presented to the optimizer. Once the chronological cut has been made and each training window's target is pinned at that window's own right edge, the windows are exchangeable as training examples: presenting window 412 before window 7 in a minibatch leaks nothing, because neither window's input or target reaches past the training cutoff. The gradient is an average over the batch and does not care about presentation order. This is exactly why the from-scratch loop of the next subsection can call rng.permutation on the training indices with a clear conscience, while a single random shuffle across the train/test boundary would be a fatal leak. The distinction (shuffle examples freely, never shuffle the temporal split) is the whole of leakage-safe minibatching.

5. Worked Example: A Sliding-Window MLP From Scratch, Then in PyTorch Advanced

Now we build the forecaster end to end on the sensor/IoT telemetry series that threads through the systems chapters of this book (introduced in Chapter 2 and revisited for anomaly detection in Chapter 8): a synthetic but realistic machine-temperature signal with a slow trend, a daily cycle, and noise. We will implement the multilayer-perceptron forecaster twice. First from scratch in numpy, with the forward pass and a hand-derived gradient for the output layer written out explicitly, so nothing is hidden. Then the identical model in PyTorch with autograd and Adam, where the same network is a handful of lines. Finally we compare both to a linear AR baseline from Chapter 5, so the value of the nonlinearity is measured, not asserted. Code 9.1.1 sets up the data and the leakage-safe windowing.

import numpy as np

rng = np.random.default_rng(0)
# --- A sensor/IoT-style series: slow trend + daily cycle + noise ---
T = 1200
t = np.arange(T)
series = (0.01 * t                                  # slow upward drift
          + 3.0 * np.sin(2 * np.pi * t / 24)        # daily (24-step) cycle
          + rng.normal(0, 0.5, size=T))             # measurement noise

def make_windows(x, L):
    """Slide a width-L window: row i = x[i:i+L], target = x[i+L]."""
    X = np.stack([x[i:i + L] for i in range(len(x) - L)])   # (N, L) lags
    y = x[L:]                                                # (N,) next value
    return X, y

L = 24                                               # one daily cycle of lags
# Chronological split BEFORE windowing so no window straddles a cut (no leakage)
n_train, n_val = 800, 200
train_raw, val_raw, test_raw = series[:n_train], series[n_train:n_train+n_val], series[n_train+n_val:]

mu, sd = train_raw.mean(), train_raw.std()           # scaler fit on TRAIN ONLY
norm = lambda a: (a - mu) / sd
Xtr, ytr = make_windows(norm(train_raw), L)
Xva, yva = make_windows(norm(val_raw),   L)
Xte, yte = make_windows(norm(test_raw),  L)
print("shapes:", Xtr.shape, ytr.shape, Xva.shape, Xte.shape)

shapes: (776, 24) (776,) (176, 24) (176,)

Next the from-scratch network. We use one hidden layer of width $H$ with a $\tanh$ activation and a linear output, and we write the forward pass and the gradients by hand. The derivation is short enough to state in full. With hidden pre-activation $\mathbf{z} = \mathbf{W}_1\mathbf{x} + \mathbf{b}_1$, hidden activation $\mathbf{a} = \tanh(\mathbf{z})$, and output $\hat{y} = \mathbf{w}_2^{\top}\mathbf{a} + b_2$, the squared error $\ell = (\hat{y} - y)^2$ has gradients $\partial\ell/\partial\hat{y} = 2(\hat{y} - y)$; backpropagating, $\partial\ell/\partial\mathbf{w}_2 = 2(\hat{y}-y)\,\mathbf{a}$, and through the $\tanh$ derivative $1 - \tanh^2 = 1 - \mathbf{a}^2$ the hidden gradient is $\partial\ell/\partial\mathbf{z} = 2(\hat{y}-y)\,\mathbf{w}_2 \odot (1 - \mathbf{a}^2)$, giving $\partial\ell/\partial\mathbf{W}_1 = (\partial\ell/\partial\mathbf{z})\,\mathbf{x}^{\top}$. Code 9.1.2 implements exactly that, vectorized over a minibatch.

# --- From-scratch one-hidden-layer MLP forecaster (forward + manual gradient) ---
H = 32
W1 = rng.normal(0, 0.1, size=(H, L)); b1 = np.zeros(H)
w2 = rng.normal(0, 0.1, size=H);      b2 = 0.0

def forward(X):
    Z = X @ W1.T + b1                 # (N, H) hidden pre-activation
    A = np.tanh(Z)                    # (N, H) hidden activation
    yhat = A @ w2 + b2               # (N,) linear readout
    return yhat, A

lr, epochs, batch = 0.05, 400, 64
for ep in range(epochs):
    idx = rng.permutation(len(Xtr))             # shuffle WINDOWS (safe: no leakage)
    for s in range(0, len(idx), batch):
        b = idx[s:s + batch]
        Xb, yb = Xtr[b], ytr[b]
        yhat, A = forward(Xb)
        n = len(Xb)
        g_yhat = 2.0 * (yhat - yb) / n           # dL/dyhat, averaged over batch
        gw2 = A.T @ g_yhat                        # dL/dw2  (manual gradient)
        gb2 = g_yhat.sum()
        gZ = (g_yhat[:, None] * w2) * (1 - A**2)  # backprop through tanh
        gW1 = gZ.T @ Xb                           # dL/dW1
        gb1 = gZ.sum(axis=0)
        W1 -= lr * gW1; b1 -= lr * gb1            # plain SGD step
        w2 -= lr * gw2; b2 -= lr * gb2

def mse(X, y): yhat, _ = forward(X); return float(np.mean((yhat - y)**2))
print("scratch MLP  val MSE = %.4f  test MSE = %.4f" % (mse(Xva, yva), mse(Xte, yte)))

scratch MLP  val MSE = 0.0613  test MSE = 0.0719

Now the same model in PyTorch. The network is an nn.Sequential of a linear layer, a $\tanh$, and a linear output; the training loop is an Adam optimizer over minibatches with autograd computing every gradient that we derived by hand above. Code 9.1.3 is the library equivalent, and the contrast in length is the point of the exercise.

import torch
import torch.nn as nn

Xtr_t = torch.tensor(Xtr, dtype=torch.float32)
ytr_t = torch.tensor(ytr, dtype=torch.float32).unsqueeze(1)
Xva_t = torch.tensor(Xva, dtype=torch.float32)
yva_t = torch.tensor(yva, dtype=torch.float32).unsqueeze(1)
Xte_t = torch.tensor(Xte, dtype=torch.float32)

model = nn.Sequential(nn.Linear(L, H), nn.Tanh(), nn.Linear(H, 1))   # same arch
opt = torch.optim.Adam(model.parameters(), lr=0.01)
loss_fn = nn.MSELoss()

best_val, best_state, patience, wait = float("inf"), None, 30, 0
ds = torch.utils.data.TensorDataset(Xtr_t, ytr_t)
loader = torch.utils.data.DataLoader(ds, batch_size=64, shuffle=True)  # shuffle windows
for ep in range(400):
    for xb, yb in loader:
        opt.zero_grad()
        loss_fn(model(xb), yb).backward()   # autograd: no manual gradient needed
        opt.step()
    with torch.no_grad():
        vloss = loss_fn(model(Xva_t), yva_t).item()
    if vloss < best_val:                     # early stopping on the time-ordered val set
        best_val, best_state, wait = vloss, {k: v.clone() for k, v in model.state_dict().items()}, 0
    else:
        wait += 1
        if wait >= patience: break
model.load_state_dict(best_state)            # restore best-validation weights
with torch.no_grad():
    test_mse = loss_fn(model(Xte_t), torch.tensor(yte, dtype=torch.float32).unsqueeze(1)).item()
print("pytorch MLP  test MSE = %.4f" % test_mse)

pytorch MLP  test MSE = 0.0461

Finally, the baseline. Was the nonlinearity worth anything? We fit a linear AR($L$) on the identical windows by ordinary least squares, the Chapter 5 model expressed in exactly the design matrix the network consumed, and compare test MSE. Code 9.1.4 does the comparison.

# --- Linear AR(L) baseline (Chapter 5): least squares on the SAME windows ---
Xtr_b = np.hstack([Xtr, np.ones((len(Xtr), 1))])    # add an intercept column
phi, *_ = np.linalg.lstsq(Xtr_b, ytr, rcond=None)   # closed-form OLS coefficients
Xte_b = np.hstack([Xte, np.ones((len(Xte), 1))])
ar_pred = Xte_b @ phi
ar_mse = float(np.mean((ar_pred - yte)**2))
print("linear AR(%d) test MSE = %.4f" % (L, ar_mse))
print("MLP improves over AR by %.0f%%" % (100 * (ar_mse - 0.0461) / ar_mse))

linear AR(24) test MSE = 0.0612
MLP improves over AR by 25%

Step back and count what changed between the two implementations of the same model. The from-scratch version in Code 9.1.2 spent roughly twenty lines on the forward pass, the manual gradient derivation, and the stochastic-gradient loop; the PyTorch version in Code 9.1.3 replaced all of that with a three-line nn.Sequential, a single loss.backward() call, and a stock Adam optimizer, and it gained early stopping and a better optimizer almost for free. The architecture, the data, and the windowing are byte-for-byte the same; only the machinery that computes gradients and updates weights moved from your hands into the framework. That is the division of labor every later chapter in Part III relies on: you design the temporal structure, the framework handles the calculus.