6 Recurrent Neural Networks

6.1 Overview

Feed-forward neural networks (FNNs) are flexible nonlinear maps from an input vector to an output, but they do not have a built-in state variable. If we want to use them for time-series forecasting, we usually have to create lagged inputs manually, such as the L most recent observations (y_t, y_{t-1}, \ldots, y_{t-L+1}). That approach is often useful, but it fixes the memory length L in advance and treats temporal dependence as a feature-engineering problem.

Recurrent neural networks (RNNs) instead maintain a hidden state that is updated as new observations arrive. This makes them closer in spirit to econometric models with recursively updated states. Examples include autoregressive moving-average (ARMA) models, generalized autoregressive conditional heteroskedasticity (GARCH) models, and state-space models, which represent dynamics through a state that is updated or filtered over time. Both classical recursions and RNNs estimate coefficients within a chosen function family. ARMA and GARCH models use parsimonious recursions for quantities such as conditional means and variances; an RNN uses a neural-network update to learn predictive features of the history. Its architecture, including the dimension of the state, the connections, and the activation function, still restricts the admissible update functions. Estimation selects the weights and biases within that family.

This chapter introduces the basic RNN architecture in a forecasting setting, interprets it through an econometric state-variable lens, explains how training works when the data are one long time series, and then turns to the central difficulty of RNNs: vanishing and exploding gradients. The basic recurrent architecture is closely associated with Elman (1990), while training through an unrolled sequence uses the backpropagation ideas of Rumelhart, Hinton, and Williams (1986) and backpropagation through time as described by Werbos (1990). The vanishing-gradient difficulty is analyzed in Bengio, Simard, and Frasconi (1994). That difficulty is the main reason why the next chapter moves to Long Short-Term Memory (LSTM) networks. Recurrent forecasting models are in routine industrial use: Uber’s forecasting team reports using an LSTM to predict rider demand around holidays and other special events, the dates at which forecast errors are most costly (Smyl, Shanmugam, and Laptev 2017).

6.2 Roadmap

1. We first define the forecasting setup and ask how a model can use the available history.
2. We contrast a lagged feed-forward network with a recurrent state update.
3. We introduce the RNN equations, parameter sharing—the same weights reused at every date—and the unrolled computational graph, in which the recurrence is drawn out over time as a chain.
4. We connect RNNs to econometric state-variable models such as GARCH.
5. We discuss deterministic versus probabilistic outputs for forecasting applications.
6. We derive backpropagation through time and analyze the vanishing and exploding gradient problems.
7. We finish with training on short sequences, contrasting resetting the memory with carrying it forward, and motivate LSTMs.

6.3 Forecasting Setup and the Information Set

Suppose we want to forecast a scalar outcome Y_{t+1} using information available at time t. Recall the forecast-origin information set \mathcal{F}_t from the time-series cross-validation discussion: it contains the variables in the form in which they were actually available at the forecast origin. In a macroeconomic or financial application, \mathcal{F}_t may contain past values of the target variable, lagged predictors, realized returns, volatility measures, survey variables, text-based indicators, and other signals known at time t.

A lagged feed-forward network uses a fixed window of this information. Write f_\theta for its prediction map, where \theta collects that network’s trainable weights and biases. For example,

\hat{y}_{t+1\mid t} = f_\theta(y_t, y_{t-1}, \ldots, y_{t-L+1}, \mathbf{x}_t, \mathbf{x}_{t-1}, \ldots, \mathbf{x}_{t-L+1}).

Here \hat{y}_{t+1\mid t} denotes a forecast of Y_{t+1} formed at the origin t; the two-index subscript records both dates, and the forecast itself is built from the realized history, following the book’s case convention. This is a reasonable baseline. It is close to how econometricians build distributed-lag or autoregressive forecasting models. The limitation is that the window length L must be chosen before estimation, and the model does not update an internal state as new observations arrive.

From this point onward, \mathbf{x}_t denotes the complete per-period input to the RNN. It may include y_t together with any other predictors known at time t. An RNN computes a hidden state \mathbf{h}_t by processing the inputs (\mathbf{x}_1, \ldots, \mathbf{x}_t) in order from a fixed initial state. The hidden state replaces the fixed lag vector with a recursively updated summary.

The hidden state \mathbf{h}_t is neither an observed economic variable nor a latent stochastic state that must be filtered. Once the network parameters and initial state are fixed, \mathbf{h}_t is computed exactly from the input history. It is a learned vector summary of the past that is useful for prediction.

Econometric Interpretation

An RNN carries information forward through a state variable, as a state-space model or a GARCH recursion does. Its state coordinates are learned predictive features rather than quantities assigned a specific econometric meaning. The neural architecture specifies the family of state-update functions, and estimation fits its weights and biases.

6.4 From a Lagged FNN to a Recurrent State

The previous section described the hidden state abstractly, as a recursively updated summary of the past. We now make the recursion explicit. A feed-forward network with lagged inputs can learn nonlinear relationships among all the lags supplied to it, but cannot access additional history unless that history is summarized in an included feature. An RNN instead processes inputs sequentially, compressing them into a hidden state with the same update rule at every step. If the state is carried forward without resets, all earlier inputs can in principle influence the current prediction through that state. The fixed-window training scheme introduced later truncates this history in estimation.

Write g_\theta for the RNN’s state-update map and q_\theta for its output map. Here \theta collects the RNN’s trainable weights and biases, a separate parameter collection from that of the lagged FNN. At time t, the state-update map combines the current input \mathbf{x}_t and the previous hidden state \mathbf{h}_{t-1} to produce the new hidden state:

\mathbf{h}_t = g_\theta(\mathbf{h}_{t-1}, \mathbf{x}_t).

The output map q_\theta then turns the state into a forecast, classification probability, or distributional parameter. For a scalar output,

\hat{y}_{t+1\mid t} = q_\theta(\mathbf{h}_t).

The essential point is that the model learns both the update rule and the predictive mapping.

6.5 RNN Architecture and Notation

The recurrence can be visualized in two ways: a compact form with a self-loop, representing the recurrence, and an “unrolled” form, which shows the flow of information through time.

(a) Compact Form

(b) Unrolled Form

Figure 6.1: RNN architecture in compact form (a) and as a three-step excerpt of the unrolled sequence (b). Circles are hidden states, the nodes below them are inputs, the nodes above them are outputs, and arrows show the flow of information. The labels ‘state update’ refer to the same rule combining the previous hidden state with the current input at every step. The loop represents recurrence; the ellipses and continuation arrows indicate omitted earlier and later steps. Each output’s two-index subscript records the target date before the bar and the forecast origin after it.

The left panel in Figure 6.1 shows the compact representation: a recurrent link feeds the previous hidden state into the next state update. The right panel draws three consecutive updates from a longer sequence, with the ellipses showing that the recurrence continues on both sides. Unrolling is mainly a computational device for training; the same update rule and parameters are reused at every date.

The output carries the two-index subscript introduced above: \hat{\mathbf{y}}_{t+1\mid t} is formed at origin t and concerns date t+1. In a scalar one-step-ahead application it is the point forecast \hat{y}_{t+1\mid t}; more generally it collects the parameters of the predictive distribution of Y_{t+1} given \mathcal{F}_t. As everywhere in the book, an unhatted lowercase y_s denotes the outcome realized at calendar time s, and uppercase Y_s the corresponding random variable.

At each time step t, define:

• \mathbf{x}_t \in \mathbb{R}^{p}: The input vector at the current time step (e.g., today’s stock returns and trading volume); p is the per-period input dimension, as in the feed-forward chapter.
• \mathbf{h}_{t-1} \in \mathbb{R}^{H}: The hidden state from the previous time step, which acts as the network’s memory of the past.
• \mathbf{h}_t \in \mathbb{R}^{H}: The new hidden state, calculated by combining the current input \mathbf{x}_t with the previous state \mathbf{h}_{t-1}.
• \hat{\mathbf{y}}_{t+1\mid t} \in \mathbb{R}^{O}: The output produced at forecast origin t for date t+1 (e.g., the prediction for the next period’s volatility), obtained from the current hidden state. The hat marks it as a model output rather than an observed quantity.

For the basic RNN, the maps g_\theta and q_\theta take the following concrete forms. Write g for the fixed hidden-state activation, applied componentwise, and g_{\mathrm{out}} for the output activation, chosen to match the forecasting object:

\begin{aligned} \mathbf{h}_t &= g_\theta(\mathbf{h}_{t-1}, \mathbf{x}_t) \\ &= g(\mathbf{W}_{hh}\mathbf{h}_{t-1} + \mathbf{W}_{xh}\mathbf{x}_t + \mathbf{b}_h), \\[0.5em] \hat{\mathbf{y}}_{t+1\mid t} &= q_\theta(\mathbf{h}_t) \\ &= g_{\mathrm{out}}(\mathbf{W}_{hy}\mathbf{h}_t + \mathbf{b}_y). \end{aligned}

The activation g is one component of the full state-update map g_\theta; it does not include the weights and biases. The output activation may be the identity for an unrestricted point forecast, a sigmoid for a binary-event probability, softmax for probabilities over several classes, or softplus for a positive variance forecast. These choices are the output activations introduced in the feed-forward networks chapter.

We fix \mathbf{h}_0 by convention, typically \mathbf{h}_0 = \mathbf{0}, since the state recursion is otherwise undefined at t=1. The trainable parameters are:

• \mathbf{W}_{hh} \in \mathbb{R}^{H \times H} (hidden-to-hidden weights)
• \mathbf{W}_{xh} \in \mathbb{R}^{H \times p} (input-to-hidden weights)
• \mathbf{W}_{hy} \in \mathbb{R}^{O \times H} (hidden-to-output weights)
• \mathbf{b}_h \in \mathbb{R}^{H}, \mathbf{b}_y \in \mathbb{R}^{O} (bias terms)

The letter pairs in the matrix names give the source then the destination: xh means input-to-hidden and hy means hidden-to-output. Matrix entries retain the convention from the feed-forward chapter: rows index receiving units and columns index sending units.

The RNN’s parameter collection is \theta = \{\mathbf{W}_{hh}, \mathbf{W}_{xh}, \mathbf{W}_{hy}, \mathbf{b}_h, \mathbf{b}_y\}. The state-update map g_\theta uses the first two matrices and \mathbf{b}_h; the output map q_\theta uses \mathbf{W}_{hy} and \mathbf{b}_y.

Optional direct input-to-output connection. A further weight matrix \mathbf W_{xy}\in\mathbb R^{O\times p} lets the current input affect the output directly. The output map then takes both the hidden state and the current input as arguments:

\begin{aligned} \hat{\mathbf y}_{t+1\mid t} &=q_\theta(\mathbf h_t,\mathbf x_t) \\ &=g_{\mathrm{out}}(\mathbf W_{hy}\mathbf h_t+\mathbf W_{xy}\mathbf x_t+\mathbf b_y). \end{aligned}

The additional matrix joins the trainable parameter collection \theta. The hidden state still provides the recurrent summary of the past. Setting \mathbf W_{xy}=\mathbf0 recovers the basic architecture used in this chapter’s figures and exercises.

Why use tanh? We use \tanh as the hidden-state activation g, recalling its shape from the activation-function discussion. Its range (-1,1) keeps every hidden-state component bounded through repeated updates. Its symmetry around zero lets a unit represent both positive and negative values. Near zero, \tanh(z)\approx z and \tanh'(0)=1, so the activation approximately preserves small changes in its argument while allowing nonlinear responses away from zero (Goodfellow, Bengio, and Courville 2016, sec. 6.3.2).

The tradeoff is saturation: for large absolute arguments, tanh flattens near \pm1, and its derivative \tanh'(z)=1-\tanh^2(z) approaches zero. Repeated multiplication by small derivatives can weaken the learning signal from distant inputs. Bounded hidden-state values do not guarantee bounded gradients, either: gradient propagation also involves \mathbf W_{hh}. We derive both effects in Section 6.9.

Before listing the properties of this recursion, we follow it through two dates with numbers.

Example: Two recurrent updates

Consider a single hidden unit with a scalar input and an identity output activation, g_{\mathrm{out}}(z)=z, so that every quantity is a scalar and we drop the boldface. The supplied parameters are W_{xh}=1, W_{hh}=0.5, b_h=0, W_{hy}=1.2, and b_y=0; the initial state is h_0=0; the inputs are x_1=0.4 and x_2=0.2. Use \tanh(0.4)\approx 0.380 and \tanh(0.39)\approx 0.371. The two updates are

\begin{aligned} h_1 &= \tanh(W_{hh}h_0 + W_{xh}x_1 + b_h) \\ &= \tanh(0.5\times 0 + 0.4) \approx 0.380, \\[0.5em] h_2 &= \tanh(W_{hh}h_1 + W_{xh}x_2 + b_h) \\ &= \tanh(0.5\times 0.380 + 0.2) \\ &= \tanh(0.39) \approx 0.371, \end{aligned}

and the output at origin 2 is \hat{y}_{3\mid 2} = W_{hy}h_2 + b_y \approx 0.445, the forecast of Y_3. The same three weights and two biases are used at both dates; only the input and the previous state change. The first input x_1 affects h_2, and hence \hat{y}_{3\mid 2}, only through h_1; the training section differentiates exactly this route.

Key Properties

• Parameter sharing: The same weight matrices (\mathbf{W}_{hh}, \mathbf{W}_{xh}, \mathbf{W}_{hy}) and biases are used at every time step. The total parameter count is therefore independent of sequence length, so the same network can be applied to sequences of different lengths.
• Memory: The hidden state \mathbf{h}_t carries information from past inputs forward, so observations at t-d for d\geq 1 can still influence \hat{y}_{t+1\mid t} through the state.

6.6 Econometric State-Variable Interpretation

An RNN models time-series dynamics through a recursively updated state. A useful econometric benchmark is the GARCH model of Bollerslev (1986).

In the GARCH(1,1) model, the return at t is decomposed as R_t = \mu_t + \varepsilon_t, where \mu_t=\mathbb{E}[R_t\mid\mathcal{F}_{t-1}] and \varepsilon_t=\sigma_t Z_t. The innovations Z_t are independent and identically distributed, independent of \mathcal{F}_{t-1}, and normalized so that \mathbb{E}[Z_t]=0 and \mathbb{E}[Z_t^2]=1. Under these conditions, \sigma_t^2=\operatorname{Var}(\varepsilon_t\mid\mathcal{F}_{t-1}). GARCH indexes the variance by the date of the innovation whose variance it describes. At forecast origin t, its next-period update is

\sigma_{t+1}^2 = \omega + \alpha \varepsilon_t^2 + \beta \sigma_t^2,

where \omega>0 and \alpha,\beta\geq 0. Under the stated normalization, covariance stationarity with finite unconditional innovation variance requires \alpha+\beta<1; that variance is \omega/(1-\alpha-\beta). The update combines the previous variance state \sigma_t^2 with the newly observed squared innovation \varepsilon_t^2 to produce \sigma_{t+1}^2.

An RNN-style volatility model could instead use

\begin{aligned} \mathbf{h}_t &= \tanh(\mathbf{W}_{hh}\mathbf{h}_{t-1} + \mathbf{W}_{xh}\mathbf{x}_t + \mathbf{b}_h) \\ \sigma_{t+1}^2 &= \text{softplus}(\mathbf{W}_{hy}\mathbf{h}_t + \mathbf{b}_y), \end{aligned}

where \mathbf{x}_t might contain lagged returns, squared returns, realized volatility measures, or macro-financial predictors known at time t. Both models use information available at t to forecast the variance at t+1. In the RNN, the hidden state \mathbf{h}_t learns a vector summary of the history, and the separate output map turns that summary into a positive variance forecast.

The analogy has clear limits. A GARCH conditional variance is an observation-driven state (Cox 1981): given the parameters, mean specification, and initial variance, it is computed recursively from past observations. Its probabilistic role is the conditional variance of the innovation; specifying the distribution of Z_t additionally gives a full conditional distribution. An RNN hidden state is likewise observation-driven, a deterministic function of past inputs given its parameters and initialization, but its components carry no assigned probabilistic role: the hidden state is a predictive representation. The RNN becomes a probabilistic econometric model once we specify how its output maps into a conditional distribution. A proper scoring rule then determines how the parameters of that model are estimated.

RNNs and Classical Time-Series Models

The hidden state \mathbf{h}_t plays a role similar to a state vector in a state-space model or a volatility state in GARCH. Like the GARCH variance, \mathbf{h}_t is exactly computable from past inputs once the model parameters and initial state are fixed. A latent state in a stochastic state-space model instead requires filtering: inference about its conditional distribution given the observations. The RNN’s state coordinates are learned predictive features, without a direct one-to-one econometric interpretation.

6.7 Outputs and Loss Functions

The output layer determines what kind of forecasting object the RNN produces. For a conditional mean forecast, we can use a linear output

\hat{y}_{t+1\mid t} = \mathbf{W}_{hy}\mathbf{h}_t + \mathbf{b}_y

and train by minimizing squared loss. For a binary event, such as recession prediction, we can pass the output through a sigmoid and train with cross-entropy. For a probabilistic forecast, the output should parameterize a distribution and be trained with a proper scoring rule. For example, suppose we impose a zero conditional mean for returns, \mu_{t+1}=0. A Gaussian RNN model then specifies

\begin{aligned} R_{t+1} \mid \mathcal{F}_t &\sim \mathcal{N}(0, \sigma_{t+1}^2), \\ \sigma_{t+1}^2 &= \text{softplus}(\mathbf{W}_{hy}\mathbf{h}_t + \mathbf{b}_y). \end{aligned}

The training score is evaluated on the realized returns r_{t+1}. The zero-mean restriction is a modeling choice; a nonzero mean can instead be estimated jointly with the variance using an additional unrestricted output, as developed in the distributional networks chapter.

A deterministic RNN trained by squared loss gives a point forecast. A distributional RNN gives a predictive distribution; negative log-likelihood is one familiar training score, but other differentiable proper scores, such as the continuous ranked probability score (CRPS), are also valid. The resulting forecasts can be evaluated with the scoring rules from the predictive-distribution chapter.

From next-token prediction to foundation models

Text is another ordered sequence. Consider continuing a central-bank statement beginning “Inflation is expected to”. We first divide the text into tokens: units such as words, word fragments, or punctuation. The prediction steps now correspond to positions in the text. An RNN takes a numerical encoding of each token as input and updates its hidden state using the same recurrence as above. A softmax output assigns probabilities to the possible next tokens in its vocabulary, the set of tokens it can predict. Training uses the observed preceding tokens and minimizes cross-entropy: the negative log probability assigned to the actual next token, summed or averaged across positions. Recurrent language models use precisely this forecasting setup (Mikolov et al. 2010).

To generate text, first process a prompt, the supplied starting text. Then sample a next token from the output distribution, feed that token into the RNN, and repeat. With the fitted weights held fixed, the hidden state changes as the generated sequence grows. This is autoregressive generation: each new token is predicted conditional on the preceding tokens.

A foundation model is trained on broad data at scale and can be adapted to many subsequent tasks (Bommasani et al. 2021). Next-token prediction can provide this initial training, called pretraining, as in Brown et al. (2020). The architecture need not be recurrent. A Transformer uses attention to combine numerical vectors representing tokens, with weights that depend on the text (Vaswani et al. 2017). For next-token prediction, it can access only positions before the token being predicted, within a limited span called the context window. An RNN carries the preceding text through its recurrent hidden state; a Transformer can access vectors at multiple positions within that window. The forecasting objective can therefore be the same even though the mechanism for using the past differs.

Forecast-Origin Discipline

For time-series applications, \mathbf{x}_t and \mathbf{h}_t must only use information available at the forecast origin. The recurrent architecture does not remove the usual econometric concern about look-ahead bias, revised macroeconomic data, publication lags, or temporal leakage in validation.

Question for Reflection

If a macroeconomic predictor is published with a one-month delay and later revised, what should be included in \mathcal{F}_t when training an RNN for real-time forecasting?

Suggested Answer

\mathcal{F}_t should contain only the data vintage available at forecast origin t. With a one-month publication delay, the current-month value is not yet available, and revisions released after t are also unavailable. Each training row should therefore use the vintage that had been released by its own forecast origin, rather than values from the revised historical series available later.

6.8 Training Through Time

Training an RNN raises a problem absent from feed-forward networks: because the same weight matrices (\mathbf{W}_{hh}, \mathbf{W}_{xh}) are applied at every time step, the gradient with respect to each matrix is the sum of contributions from every time step. The algorithm that handles this is backpropagation through time (BPTT).

The core idea of BPTT is to unroll the network through time, as shown in Figure 6.1. This conceptually transforms the RNN into a deep feed-forward network where each time step becomes a layer, but with the critical constraint that the weights are shared across all layers. Once unrolled, standard backpropagation can be applied.

For clarity, the chapter specializes to a single terminal loss evaluated at the terminal step T, written in the forecasting convention: the step-T output is scored against the outcome it forecasts,

L_T = \ell(y_{T+1}, \hat{\mathbf{y}}_{T+1\mid T}),

with the observed target first and the network output second, as in the feed-forward networks chapter. For an objective with a loss at every step, the gradient is obtained by summing analogous contributions over the per-step losses \sum_{t=1}^{T}\ell(y_{t+1},\hat{\mathbf{y}}_{t+1\mid t}).

Notation: Gradients and Jacobians

All vectors are columns. We use the conventions from feed-forward backpropagation, with i and j indexing hidden units (matrix rows and columns):

• The loss gradient \nabla_{\mathbf{h}_t}L_T is an H\times1 column vector; its i-th entry is \partial L_T/\partial(\mathbf{h}_t)_i.
• The loss gradient \nabla_{\mathbf{W}_{hh}}L_T is an H\times H matrix; its (i,j)-th entry is \partial L_T/\partial(\mathbf{W}_{hh})_{ij}. A matrix gradient has the same shape as the parameter matrix it updates.
• The vector derivative \partial\mathbf{h}_t/\partial\mathbf{h}_s is an H\times H Jacobian; its (i,j)-th entry is \partial(\mathbf{h}_t)_i/\partial(\mathbf{h}_s)_j. Rows index the resulting state, and columns index the state being perturbed.

We use \nabla for scalar-loss gradients with respect to vectors or matrices, and \partial for ordinary scalar derivatives and vector-to-vector Jacobians. The chain rule propagates column gradients backward by multiplying by the transpose of the forward Jacobian.

From one weight to a matrix gradient. Define the pre-activation \mathbf{z}_t=\mathbf{W}_{hh}\mathbf{h}_{t-1}+\mathbf{W}_{xh}\mathbf{x}_t+\mathbf{b}_h. As in feed-forward backpropagation, let \boldsymbol{\delta}_t be the gradient of the terminal loss with respect to this pre-activation:

\begin{aligned} \boldsymbol{\delta}_t &= \nabla_{\mathbf{z}_t}L_T \\ &= \operatorname{diag}\big(g'(\mathbf{z}_t)\big) \nabla_{\mathbf{h}_t}L_T. \end{aligned}

The vector g'(\mathbf{z}_t) contains the componentwise activation derivatives, and \operatorname{diag} places them on the diagonal of an H\times H matrix. Thus \boldsymbol{\delta}_t is an H\times1 column vector that includes all downstream effects on L_T.

For a single recurrent weight (\mathbf{W}_{hh})_{ij}, its direct effect on (\mathbf{z}_t)_i at step t is (\mathbf{h}_{t-1})_j, holding the previous state fixed. Summing the contributions from its repeated use gives

\frac{\partial L_T}{\partial(\mathbf{W}_{hh})_{ij}} = \sum_{t=1}^{T}(\boldsymbol{\delta}_t)_i(\mathbf{h}_{t-1})_j.

Stacking these scalar derivatives yields the matrix gradient

\nabla_{\mathbf{W}_{hh}}L_T = \sum_{t=1}^{T}\boldsymbol{\delta}_t\mathbf{h}_{t-1}^{\top}.

Each summand is an outer product: an H\times1 column times a 1\times H row, giving an H\times H matrix. The previous state is held fixed only when differentiating the local step-t update; its own dependence on the shared weights enters through the other time-step contributions. Including that dependence again in the local factor would double-count paths. This componentwise construction avoids treating \partial\mathbf{h}_t/\partial\mathbf{W}_{hh}, which would have three coordinate indices, as an ordinary matrix.

Propagating the state gradient backward. The column gradient \nabla_{\mathbf{h}_t}L_T measures how an adjustment to the hidden state at time t affects the terminal loss. The matrix chain rule assembles the ordinary scalar chain rule into a matrix product. With the inputs and parameters held fixed, the terminal loss depends on \mathbf h_t through \mathbf h_T. For coordinate j, the scalar chain rule gives

\frac{\partial L_T}{\partial(\mathbf h_t)_j} = \sum_{i=1}^{H} \frac{\partial L_T}{\partial(\mathbf h_T)_i} \frac{\partial(\mathbf h_T)_i}{\partial(\mathbf h_t)_j}.

The last factor is the (i,j)-th entry of the forward Jacobian \partial\mathbf h_T/\partial\mathbf h_t. Thus the sum is the j-th entry of the following matrix product:

\nabla_{\mathbf{h}_t}L_T = \left(\frac{\partial\mathbf{h}_T}{\partial\mathbf{h}_t}\right)^{\top} \nabla_{\mathbf{h}_T}L_T.

The dimensions are (H\times H)(H\times1)=H\times1. The transpose makes the sum run over the later-state coordinates i, producing one gradient entry for each earlier-state coordinate j. In general, for a scalar loss depending on a column vector through a differentiable vector map, the gradient with respect to the input vector equals the map’s Jacobian transpose times the gradient with respect to its output vector.

At the terminal step, \nabla_{\mathbf{h}_T}L_T differentiates through the output map q_\theta, including its output activation g_{\mathrm{out}}; an identity output activation is not required.

For a composition of vector maps, the matrix chain rule multiplies their forward Jacobians in order. For an earlier date s<t, the forward state Jacobian is therefore the ordered product

\begin{aligned} \frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_s} &= \frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_{t-1}}\, \frac{\partial \mathbf{h}_{t-1}}{\partial \mathbf{h}_{t-2}}\cdots \frac{\partial \mathbf{h}_{s+1}}{\partial \mathbf{h}_{s}}, \\[0.5em] \frac{\partial \mathbf{h}_k}{\partial \mathbf{h}_{k-1}} &= \operatorname{diag}\big(g'(\mathbf{z}_k)\big)\,\mathbf{W}_{hh}, \end{aligned}

where k indexes time and each factor is H\times H. Each one-step Jacobian follows from \mathbf{h}_k = g(\mathbf{z}_k) with \mathbf{z}_k affine in \mathbf{h}_{k-1}. The factors need not commute, so their order matters. For backward propagation, the one-step transpose is \mathbf{W}_{hh}^{\top}\operatorname{diag}(g'(\mathbf{z}_k)): transposing reverses the multiplication order.

Substituting the ordered Jacobian product into the gradient formula and reversing the factor order under transposition gives:

\begin{aligned} \nabla_{\mathbf h_t}L_T &= \left(\frac{\partial \mathbf{h}_{t+1}}{\partial \mathbf{h}_{t}}\right)^{\top} \cdots \left(\frac{\partial \mathbf{h}_{T}}{\partial \mathbf{h}_{T-1}}\right)^{\top} \nabla_{\mathbf h_T}L_T \\[0.5em] &= \mathbf W_{hh}^{\top}\operatorname{diag}\big(g'(\mathbf z_{t+1})\big) \cdots \mathbf W_{hh}^{\top}\operatorname{diag}\big(g'(\mathbf z_{T})\big)\, \nabla_{\mathbf h_T}L_T . \end{aligned}

Read from the right, the two factors nearest the terminal gradient belong to date T: multiplying by \operatorname{diag}(g'(\mathbf z_T)) gives \boldsymbol\delta_T, and multiplying by \mathbf W_{hh}^{\top} gives \nabla_{\mathbf h_{T-1}}L_T. Each further pair of factors moves the gradient one date earlier. The error signal at date t therefore contains the whole product from t to T,

\boldsymbol\delta_t = \operatorname{diag}\big(g'(\mathbf z_t)\big) \left(\frac{\partial\mathbf h_T}{\partial\mathbf h_t}\right)^{\top} \nabla_{\mathbf h_T}L_T,

where for t=T the Jacobian product is empty and equals the identity, and the matrix gradient \sum_{t=1}^{T}\boldsymbol\delta_t\mathbf h_{t-1}^{\top} adds these carried gradients over the dates at which the shared weights were used.

The backward calculation can therefore be implemented one step at a time, and the product itself is never formed. Starting from \nabla_{\mathbf h_T}L_T, compute \boldsymbol\delta_t as above, accumulate the weight-gradient contribution, and propagate to the previous state using

\nabla_{\mathbf h_{t-1}}L_T = \mathbf W_{hh}^{\top}\boldsymbol\delta_t, \qquad t=T,T-1,\ldots,1.

When H=1, gradients and Jacobians reduce to scalar derivatives and the transposes have no effect; Exercise 6.1 develops that scalar case in detail.

6.9 Vanishing and Exploding Gradients

The long product of Jacobians is the source of the vanishing gradient problem (Bengio, Simard, and Frasconi 1994). Here \lVert\cdot\rVert denotes the Euclidean norm for vectors and its induced matrix norm, the spectral norm (the largest singular value). If each one-step Jacobian along the trajectory satisfies \lVert \operatorname{diag}(g'(\mathbf{z}_k))\,\mathbf{W}_{hh} \rVert \leq \kappa for some 0\leq\kappa < 1, the norm of the product is at most \kappa^{T-t}. This follows from submultiplicativity: the norm of a matrix product is at most the product of the matrix norms. For tanh, \lVert\mathbf W_{hh}\rVert<1 is a sufficient condition, since the activation derivatives lie between zero and one. For a fixed terminal loss, the gradient propagated to an early state is therefore multiplied by a geometrically small factor. When the remaining local derivative factors are bounded, early-step contributions to the parameter gradient can become numerically negligible.

For econometric forecasting, this is the central limitation of a plain RNN. Credit assignment means using the loss to determine how earlier state updates should change. Training can fail when this learning signal must cross a long time gap.

The same recurrent-state Jacobian product can also grow, in which case gradients explode rather than vanish. Since g' \leq 1 for the \tanh activation, each one-step Jacobian is bounded in norm by \lVert\mathbf{W}_{hh}\rVert, so \lVert\mathbf{W}_{hh}\rVert > 1 is necessary for this upper bound to permit exponential pathwise growth. It is not sufficient: whether the product actually grows depends on the entire trajectory of factors \operatorname{diag}(g'(\mathbf{z}_k))\,\mathbf{W}_{hh}. A per-factor bound below one forces decay, but individual factors with norms above one do not by themselves force growth. Nor does \lVert\mathbf{W}_{hh}\rVert\leq 1 uniformly bound every total parameter gradient, because contributions from many loss terms can still accumulate with sequence length.

In the scalar case, the pathwise picture is transparent: the product is \prod_{k=t+1}^{T} |W_{hh}|\,g'(z_k) and grows geometrically whenever the factors stay above some \gamma > 1 along the trajectory. Exploding gradients show up as numerical overflow or parameter updates that jump erratically between iterations, rather than as a silent failure to learn long-run dependence.

The standard remedy is gradient clipping (Pascanu, Mikolov, and Bengio 2013): if the norm of the stacked raw gradient vector \mathbf{v} exceeds a threshold c>0, replace \mathbf{v} by c\,\mathbf{v}/\lVert\mathbf{v}\rVert before the optimizer sees it. This transformation caps the raw gradient norm and preserves the raw gradient direction. Under plain gradient descent, it also caps the norm of the immediate parameter step, up to the learning rate; with momentum or adaptive optimizers, the eventual update also depends on stored optimizer state. Clipping stabilizes the exploding case but does nothing for vanishing gradients, whose problem is the loss of signal across a long credit-assignment path.

6.10 Training on Short Sequences

Backpropagating through an entire long series is costly, so we train on shorter sequences. A window of \tau inputs ending at date t is paired with the target y_{t+1}. For \tau=3, the first two training pairs are

\begin{aligned} (\mathbf x_1,\mathbf x_2,\mathbf x_3)&\longrightarrow y_4,\\ (\mathbf x_2,\mathbf x_3,\mathbf x_4)&\longrightarrow y_5. \end{aligned}

The inputs remain in chronological order within each window.

In fixed-window training, the hidden state starts at zero for each window. Alternatively, we process consecutive, non-overlapping segments of up to \tau steps and carry the final state into the next segment. At each boundary, we detach that state: retain its numerical value, but treat it as a constant in the next backward pass.

The state horizon is how far back the forward recursion extends; the gradient horizon is how far back we differentiate:

Scheme	State horizon	Gradient horizon
Reset each window	\tau steps	At most \tau steps
Carry and detach the state	Can extend beyond \tau	Stops at the segment boundary

Cutting the backward path at these boundaries is truncated backpropagation through time (Williams and Peng 1990). The forecast can still depend on earlier inputs through a carried state. In the current segment’s backward pass, however, shared weights receive contributions only from their uses within that segment. With fitted weights fixed, a reset window can use older information only if its inputs already contain summaries of that information.

Overlapping windows are dependent observations; forecast evaluation still requires chronological splits. Neither scheme by itself resolves a structural break absent from the inputs or state, since both reuse the same maps across dates. The empirical chapter uses reset windows and provides the window diagram and full Python workflow.

Bridge to LSTM Networks

Shorter backward passes reduce computation, but do not repair vanishing gradients. The next chapter introduces LSTMs, which modify the state update to preserve information and gradient signals over longer sequences.

6.11 Summary

Key Takeaways

1. RNNs use shared parameters to update a hidden state that summarizes past inputs.
2. Backpropagation through time applies the chain rule to the unrolled sequence of state updates.
3. Repeated state transitions can attenuate or amplify gradients, making long-range learning difficult.
4. Reset windows restrict both forward and backward paths; detachment cuts the backward path while retaining the carried state.

Common Pitfalls

• Hidden-state coordinates do not automatically correspond to identified economic states.
• A deterministic RNN output is not a full predictive distribution without a probabilistic output model.
• Keep observations ordered within windows and keep future targets out of the training split.
• Gradient clipping limits large gradients but does not repair vanishing gradients.

6.12 Exercises

Exercise 6.1: Vanishing Gradient Problem in RNNs

Consider a scalar RNN with fixed initial state h_0=0:

\begin{aligned} z_t &= W_{hh}h_{t-1}+W_{xh}x_t+b_h, \\ h_t &= \tanh(z_t), \\ \hat y_{t+1\mid t} &= W_{hy}h_t+b_y. \end{aligned}

Let L_T=\ell(y_{T+1},\hat y_{T+1\mid T}) be a differentiable terminal loss. In Parts 1–2, hold the observed sequence, terminal date T, and fitted parameter values fixed, and assume \partial L_T/\partial h_T is finite.

1. Derive the scalar BPTT gradient. Express \partial L_T/\partial W_{hh} as a sum of contributions from the repeated uses of W_{hh}. Evaluate the local derivative of h_t with respect to W_{hh}, holding h_{t-1} fixed, and express \partial h_T/\partial h_t as a product of one-step derivatives.

2. Bound contributions from earlier steps. Prove that |\partial h_T/\partial h_t|\leq|W_{hh}|^{T-t} and use h_0=0 and the range of tanh to bound each summand from Part 1. Explain what |W_{hh}|<1 implies for contributions across long time gaps, and why it does not require the total parameter gradient to vanish.

3. Local versus delayed signals. Suppose the RNN is trained with squared loss at each forecast origin. Consider two data-generating processes with finite second moments. In each process, the errors U_{t+1} are independent and identically distributed with mean zero, and independent of all inputs available at time t.

   • First-order autoregressive (AR(1)) process: Y_{t+1}=\phi Y_t+U_{t+1} with |\phi|<1. The RNN’s input at date t is x_t=y_t.
   • Delayed signal: Y_{t+1}=\psi X_{t-d}+U_{t+1} with \psi\ne0 and d\geq1. The inputs X_s are independent and identically distributed with positive variance; the RNN receives only their realizations x_s.

   Work at t\geq d+1, with no state reset from t-d through t. For each process, write the conditional mean of Y_{t+1} given the inputs up to t and identify the input appearing in it. Then count the state updates after that input has entered the hidden state, up to forecast origin t. Each such update contributes a factor W_{hh}\tanh'(z_k) to the derivative of the forecast with respect to that input. Which process requires the training signal to survive a long product when d is large, and why does \phi close to one not by itself create such a product?

Exam level. The three parts cover the gradient calculation, its contraction bound, and the forecasting interpretation.

Hint for Part 1

Temporarily regard the recurrent weight at each time step as a separate copy. Differentiate with respect to each copy, holding that step’s previous state fixed, then add the contributions because the copies share one parameter.

Hint for Part 3

The input that determines the conditional mean enters the state at one specific date. Start counting with the next state update, and stop at date t. The final step from h_t to \hat y_{t+1\mid t} is the output map and adds no state update.

Solution

Part 1: Scalar backpropagation through time

The local weight derivative and the state sensitivity are

\begin{aligned} \left.\frac{\partial h_t}{\partial W_{hh}}\right|_{h_{t-1}\,\mathrm{fixed}} &= \tanh'(z_t)h_{t-1}, \\[0.5em] \frac{\partial h_T}{\partial h_t} &= \prod_{k=t+1}^{T} W_{hh}\tanh'(z_k), \end{aligned}

because each one-step state derivative is W_{hh}\tanh'(z_k). The empty product equals one when t=T. Applying the chain rule to every use of the shared weight gives

\frac{\partial L_T}{\partial W_{hh}} = \sum_{t=1}^{T} \frac{\partial L_T}{\partial h_T} \frac{\partial h_T}{\partial h_t} \tanh'(z_t)h_{t-1}.

The local factor holds h_{t-1} fixed; the dependence of that state on earlier uses of W_{hh} is already included in the other summands. Including it again would double-count paths.

Part 2: Geometric attenuation

Since 0\leq\tanh'(z)=1-\tanh^2(z)\leq1,

\begin{aligned} \left|\frac{\partial h_T}{\partial h_t}\right| &= \prod_{k=t+1}^{T}|W_{hh}|\tanh'(z_k) \\ &\leq |W_{hh}|^{T-t}. \end{aligned}

Also, h_0=0 and subsequent states are tanh outputs, so |h_{t-1}|\leq1 and the local weight derivative has magnitude at most one. The absolute value of the t-th summand is therefore bounded by

\left|\frac{\partial L_T}{\partial h_T}\right|\,|W_{hh}|^{T-t}.

The terminal loss derivative is a common finite factor under the stated fixed-sequence convention. If |W_{hh}|<1, the bound is geometrically smaller for earlier steps, weakening the learning signal across long gaps. Recent-step contributions need not be small, so the sum of all contributions need not vanish.

Part 3: Local versus delayed signals

Under squared loss, the optimal forecast given the available inputs is the conditional mean. Since U_{t+1} has mean zero and is independent of those inputs, the AR(1) and delayed-signal processes give, respectively,

\begin{aligned} \mathbb E[Y_{t+1}\mid Y_1,\ldots,Y_t] &=\phi Y_t,\\ \mathbb E[Y_{t+1}\mid X_1,\ldots,X_t] &=\psi X_{t-d}. \end{aligned}

The realized optimal forecasts are therefore \phi y_t and \psi x_{t-d}. The inputs appearing in them are x_t=y_t and x_{t-d}. In the delayed task, the independence of the X_s and the error means that the intervening inputs x_{t-d+1},\ldots,x_t provide no information about Y_{t+1}.

In the AR(1) task, the input enters at date t, following the path x_t\to h_t\to\hat y_{t+1\mid t}. There are zero state updates after input entry. To see the corresponding derivative, hold the network parameters and all other inputs fixed:

\frac{\partial\hat y_{t+1\mid t}}{\partial x_t} =W_{hy}\tanh'(z_t)W_{xh}.

The input-to-state and output maps contribute factors, but there is no W_{hh} factor.

In the delayed task, x_{t-d} first enters h_{t-d}. The subsequent state updates occur at dates t-d+1,\ldots,t, giving d state updates after input entry. Under the same derivative convention,

\begin{aligned} \frac{\partial\hat y_{t+1\mid t}}{\partial x_{t-d}} &=W_{hy}\left[\prod_{k=t-d+1}^{t}W_{hh}\tanh'(z_k)\right]\\ &\qquad\times\tanh'(z_{t-d})W_{xh}. \end{aligned}

The bracket contains exactly d recurrent factors. By Part 2, its magnitude satisfies

\left|\prod_{k=t-d+1}^{t}W_{hh}\tanh'(z_k)\right| \leq |W_{hh}|^d.

When |W_{hh}|<1 and d is large, this factor attenuates the learning signal from the delayed input, even though the forecast itself is only one step ahead. A coefficient \phi close to one makes the AR(1) series persistent, but its conditional mean depends only on the current observation, already supplied as x_t=y_t.

Exercise 6.2: Manual RNN Forward Pass

Consider the simple RNN shown in Figure 6.2: a single hidden unit processing a sequence of length 3, with hidden update h_t = \tanh(W_{hh} h_{t-1} + W_{xh} x_t + b_h) and linear output \hat{y}_{t+1\mid t} = W_{hy} h_t + b_y (identity output activation).

Figure 6.2: The single-unit RNN of this exercise, unrolled over t = 0, 1, 2, 3. Blue circles are the scalar hidden states h_t, initialized at h_0 = 0; green circles are the inputs x_t with their numerical values printed below; red circles are the outputs \hat{y}_{t+1\mid t}, the forecasts formed at origin t. Arrow labels give the scalar weight applied along each edge: W_{xh} (input-to-hidden), W_{hh} (hidden-to-hidden), and W_{hy} (hidden-to-output).

The parameters are:

• Weight for input-to-hidden: W_{xh} = 0.5
• Weight for hidden-to-hidden: W_{hh} = 0.8
• Bias for hidden unit: b_h = 0.1
• Weight for hidden-to-output: W_{hy} = 1.2
• Bias for output: b_y = 0.0

The input sequence is (x_1, x_2, x_3) = (0.2, -0.5, 0.3) and the initial hidden state is h_0 = 0. Use the approximations \tanh(0.2) \approx 0.197, \tanh(0.008) \approx 0.008, and \tanh(0.256) \approx 0.251.

1. Compute the forward pass. For t = 1, 2, 3, compute the hidden state h_t and the output \hat{y}_{t+1\mid t}.
2. Separate value from sensitivity. Using your results from Part 1, show numerically how the memory term W_{hh} h_1 and the bias b_h offset the input contribution at t=2. Then compute the sensitivity \partial h_2/\partial x_1 and state what the value of h_2 does and does not imply about the memory of x_1, connecting your answer to the geometric decay factor from Exercise 6.1.

Exam level, approximately 10–15 minutes. Part 1 checks the recursion, and Part 2 separates state value from sensitivity.

Hint for Part 2

For the sensitivity, apply the chain rule along x_1 \to h_1 \to h_2 and recall that \tanh'(z)=1-\tanh^2(z).

Solution

Part 1: Forward pass

At t = 1:

\begin{aligned} h_1 &= \tanh(0.8 \cdot 0 + 0.5 \cdot 0.2 + 0.1) \\ &= \tanh(0.2) \approx 0.197, \\[0.5em] \hat{y}_{2\mid 1} &= 1.2 \cdot 0.197 \approx 0.236 \end{aligned}

At t = 2:

\begin{aligned} h_2 &= \tanh(0.8 \cdot 0.197 + 0.5 \cdot (-0.5) + 0.1) \\ &= \tanh(0.158 - 0.25 + 0.1) \\ &= \tanh(0.008) \approx 0.008, \\[0.5em] \hat{y}_{3\mid 2} &= 1.2 \cdot 0.008 \approx 0.010 \end{aligned}

At t = 3:

\begin{aligned} h_3 &= \tanh(0.8 \cdot 0.008 + 0.5 \cdot 0.3 + 0.1) \\ &= \tanh(0.006 + 0.15 + 0.1) \\ &= \tanh(0.256) \approx 0.251, \\[0.5em] \hat{y}_{4\mid 3} &= 1.2 \cdot 0.251 \approx 0.301 \end{aligned}

Part 2: A near-zero state value is not the same as lost memory

At t = 2, the input contribution W_{xh} x_2 = 0.5 \times (-0.5) = -0.25 is almost exactly offset by the memory term W_{hh} h_1 = 0.8 \times 0.197 = 0.158 plus the bias b_h = 0.1, giving a pre-activation near zero and thus h_2 \approx 0. This is cancellation: two sizable contributions of opposite sign happen to sum to nearly zero.

Cancellation of the value does not mean the memory of x_1 is gone. Memory in the gradient sense is the sensitivity

\begin{aligned} \frac{\partial h_2}{\partial x_1} &= \underbrace{\tanh'(z_2)\, W_{hh}}_{\partial h_2/\partial h_1} \cdot \underbrace{\tanh'(z_1)\, W_{xh}}_{\partial h_1/\partial x_1} \\[0.5em] &\approx (1.000 \times 0.8) \times (0.961 \times 0.5) \\ &\approx 0.3844 \approx 0.38, \end{aligned}

using \tanh'(z) = 1 - \tanh^2(z) with z_1 = 0.2 and z_2 = 0.008. A sensitivity of 0.38 after one recurrence step is substantial—in fact, because z_2 \approx 0 places \tanh in its steepest region, h_2 is close to maximally sensitive to perturbations of h_1. Had x_1 been slightly different, h_2 would have been noticeably different: the state still remembers x_1 even though its value happens to sit near zero.

The genuine connection to Exercise 6.1 runs through the rate of decay, not the value of the state. Each additional recurrence step multiplies the sensitivity by at most W_{hh}\,\max_z \tanh'(z) = 0.8, so the influence of x_1 fades geometrically over many steps (0.8^{10} \approx 0.11, 0.8^{30} \approx 0.001). The two-part conclusion: after one step, the memory of x_1 is still strong despite the near-zero state value; over long horizons, it is the geometric factor from Exercise 6.1—not value cancellation—that makes plain RNNs forget.

6.13 References

Bengio, Yoshua, Patrice Simard, and Paolo Frasconi. 1994. “Learning Long-Term Dependencies with Gradient Descent Is Difficult.” IEEE Transactions on Neural Networks 5 (2): 157–66. https://doi.org/10.1109/72.279181.

Bollerslev, Tim. 1986. “Generalized Autoregressive Conditional Heteroskedasticity.” Journal of Econometrics 31 (3): 307–27. https://doi.org/10.1016/0304-4076(86)90063-1.

Bommasani, Rishi, Drew A. Hudson, Ehsan Adeli, Russ Altman, Simran Arora, et al. 2021. “On the Opportunities and Risks of Foundation Models.” Stanford Center for Research on Foundation Models (CRFM). https://arxiv.org/abs/2108.07258.

Brown, Tom B., Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, et al. 2020. “Language Models Are Few-Shot Learners.” In Advances in Neural Information Processing Systems 33, 1877–1901.

Cox, D. R. 1981. “Statistical Analysis of Time Series: Some Recent Developments.” Scandinavian Journal of Statistics 8 (2): 93–115.

Elman, Jeffrey L. 1990. “Finding Structure in Time.” Cognitive Science 14 (2): 179–211. https://doi.org/10.1207/s15516709cog1402_1.

Goodfellow, Ian, Yoshua Bengio, and Aaron Courville. 2016. Deep Learning. MIT Press. https://www.deeplearningbook.org.

Mikolov, Tomáš, Martin Karafiát, Lukáš Burget, Jan Černocký, and Sanjeev Khudanpur. 2010. “Recurrent Neural Network Based Language Model.” In Proceedings of Interspeech 2010, 1045–48. https://doi.org/10.21437/Interspeech.2010-343.

Pascanu, Razvan, Tomas Mikolov, and Yoshua Bengio. 2013. “On the Difficulty of Training Recurrent Neural Networks.” In Proceedings of the 30th International Conference on Machine Learning, 1310–18. Atlanta.

Rumelhart, David E., Geoffrey E. Hinton, and Ronald J. Williams. 1986. “Learning Representations by Back-Propagating Errors.” Nature 323: 533–36. https://doi.org/10.1038/323533a0.

Smyl, Slawek, Santhosh Shanmugam, and Nikolay Laptev. 2017. “Engineering Extreme Event Forecasting at Uber with Recurrent Neural Networks.” Uber Engineering Blog. https://www.uber.com/blog/neural-networks/.

Vaswani, Ashish, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017. “Attention Is All You Need.” In Advances in Neural Information Processing Systems 30, 5998–6008.

Werbos, Paul J. 1990. “Backpropagation Through Time: What It Does and How to Do It.” Proceedings of the IEEE 78 (10): 1550–60. https://doi.org/10.1109/5.58337.

Williams, Ronald J., and Jing Peng. 1990. “An Efficient Gradient-Based Algorithm for on-Line Training of Recurrent Network Trajectories.” Neural Computation 2 (4): 490–501. https://doi.org/10.1162/neco.1990.2.4.490.