Workbook 1 — What Is a Time Series?

STA 542 · Introduction to Time Series Analysis

\[\def\E{\mathbb{E}} \def\P{\mathbb{P}} \def\R{\mathbb{R}} \def\Z{\mathbb{Z}} \def\N{\mathbb{N}} \def\cH{\mathcal{H}} \def\cP{\mathcal{P}} \def\Var{\operatorname{Var}} \def\Cov{\operatorname{Cov}} \def\Corr{\operatorname{Corr}}\]

# Load the packages that supply the course data and dated series.
library(astsa)   # course data: gtemp_land, djia, soi, Hare, Lynx
library(xts)     # djia is stored as an xts (dated time series) object

Time Series Data

A finite time series record is an ordered collection of observed elements, \[x_1,x_2,\ldots,x_T.\] Here \(T\) is the number of observed time points and \(t\in\{1,\ldots,T\}\) is the time index. The symbol \(x_t\) means the element observed at time \(t\). The subscripts are part of the data: we will not be assuming the \(x_t\)’s to be exchangeable.

Yearly temperatures, daily stock prices, monthly sales, hourly traffic counts, and a voltage sampled a thousand times per second all have this form. What is observed at one time can be a single number or a richer object.

The simplest case records one real number at every time: \[x_t\in\R, \qquad t=1,\ldots,T.\] The complete record can then be displayed as the \(T\times1\) column \[x= \begin{pmatrix} x_1\\x_2\\\vdots\\x_T \end{pmatrix} \in\R^{T\times1}.\] A time plot graphs the points \((t,x_t)\), putting the time index on the horizontal axis and the observed real value on the vertical axis.

Sometimes \(d\) related numbers are observed at each time. We then write \[x_t= \begin{pmatrix} x_{t1}\\x_{t2}\\\vdots\\x_{td} \end{pmatrix} \in\R^d, \qquad t=1,\ldots,T.\] The first subscript \(t\) says when; the second subscript says which component.

The element indexed by time can be richer still. A satellite image observed at time \(t\) can be represented by a matrix \[x_t=(x_{t,ij})\in\R^{r\times c},\] with one entry for each of \(r\) rows and \(c\) columns of pixels. A complete ECG curve observed at time \(t\) can be written as a function \[x_t:\mathcal U\to\R,\] where \(u\in\mathcal U\) indexes location within the curve. In every case the subscript \(t\) still answers the same question: which image, curve, vector, or number was observed at that time? This course develops scalar and vector series directly; finite matrices can also be stacked into vectors, while function-valued series require additional tools built on the same indexing idea.

Example 1.1 (A tour of four series). The figure shows the four real series that accompany the course. Global land temperature: annual mean anomalies, °C relative to 1991–2020, 1850–2023 (astsa::gtemp_land, updated from Hansen et al. 2006). Dow Jones Industrial Average: daily closing values over 2518 trading days, April 2006 – April 2016 (astsa::djia). Southern Oscillation Index: a monthly standardized air-pressure difference between Tahiti and Darwin used to track El Niño–Southern Oscillation (ENSO) conditions, 1950–1987 (astsa::soi). Hare and lynx: pelts traded by the Hudson’s Bay Company, annually 1845–1935, in thousands (astsa::Hare, astsa::Lynx) — two series moving in linked cycles, and thus a first vector-valued series: one point of \(\R^2\) per year. Have a look at each of the plots and see if there are visible trends, wandering, cycles, and possibly interaction in each of the panels. ◊

# Plot four real series to compare their temporal patterns.
par(mfrow = c(2, 2), mar = c(2.5, 4, 2, 1))
plot(gtemp_land, col = 4, ylab = "Anomaly (°C)", main = "Land temperature")
plot(index(djia), as.numeric(djia$Close), type = "l", col = 4,
     xlab = "", ylab = "Close", main = "DJIA")
plot(soi, col = 4, ylab = "Index", main = "SOI")
plot(Hare, col = 4, ylab = "Pelts (thousands)", main = "Hare and lynx")
lines(Lynx, col = 2)

Figure 1: Four real series: annual land temperature anomaly (1850–2023), DJIA daily close (2006–2016), monthly SOI (1950–1987), and hare and lynx pelts (1845–1935). Real data: R package astsa.

What makes such data a subject of their own is the order: past values influence future values, and the index carries information. The contrast is with i.i.d. — more generally, exchangeable — data, where every reordering of the sample is statistically identical. Shuffle the temperature record and its histogram is untouched, while the trend and the clustering of warm years are destroyed.

# Compare the recorded temperature path with a shuffled version.
set.seed(542)
g <- as.numeric(gtemp_land)
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(1850:2023, g, type = "l", col = 4,
     xlab = "Year", ylab = "Anomaly (°C)", main = "As recorded")
plot(1850:2023, sample(g), type = "l", col = 4,
     xlab = "", ylab = "Anomaly (°C)", main = "Same numbers, shuffled")

Figure 2: The land temperature anomalies as recorded, and the same 174 numbers in random order. Real data: astsa::gtemp_land; right panel deliberately shuffled.

The temporal component need not be clock time: any ordered data invite the same techniques — the words of a paragraph arrive in order, and for a paragraph that makes sense they are anything but independent. Time is the leading case.

We will study methods to analyze and model time series data. The theme throughout: dependence changes what data can tell us — standard errors, confidence intervals, and train–test splits rest on independence assumptions, and under dependence they can be off by an order of magnitude.

Modelling timeseries: the difference with iid data

Statistical inference is inductive inference whose evidence is data generated by an unknown process involving randomness; because the data involve randomness, the conclusions are uncertain. Probability reasons forward — if we knew the process, what might we observe? — and statistics reasons backward: given what we observed, what can we learn about the process? The backward step needs candidates to reason among. A statistical model is that list: a collection of candidate probability distributions for the random object the data represent. Each member is a complete forward description; inference weighs the members against what was observed.

The first building block is the usual one. A random vector \(X\) with values in \(\R^d\) is specified by its cumulative distribution function, \[F(c) \;=\; \P\left(X \le c\right), \qquad c \in \R^d,\] the inequality read componentwise. The CDF is the object we actually touch: writing one down — or a density or probability mass function that yields it — is what ‘specifying a distribution’ means throughout the course.

An i.i.d. model constructs this collection economically. Start with a collection \(\cP\) of candidate CDFs for one \(\R^d\)-valued observation. For each \(F\in\cP\), independence and identical distribution determine the joint CDF of the entire random sample \(X_1,\ldots,X_T\): \[F_{1:T}(c_1,\ldots,c_T) :=\P\left(X_1\le c_1,\ldots,X_T\le c_T\right) =\prod_{t=1}^T F(c_t).\] The statistical model is still a collection of joint distributions. The i.i.d. assumption means that each member of that collection is determined by one candidate CDF for a single observation. Estimation selects a member of \(\cP\) (by maximum likelihood, for example, or through a Bayesian posterior over \(\cP\)); prediction simulates fresh draws from the selected member(s) of the model. In this precise sense, to write down an i.i.d. model is to write down a collection of one-observation CDFs.

For time series, specifying marginals at each time point won’t result in a model that captures dependence across time — as settled by the following example.

Example 1.2 (Same marginals, different processes). Two processes built from standard normal ingredients. Process A is i.i.d.: \(X_t \overset{\text{iid}}{\sim} N(0,1)\). Process B draws a single \(X_0 \sim N(0,1)\) once and repeats it: \(X_t := X_0\) for all \(t\). Under either process every \(X_t\) is exactly \(N(0,1)\), so any model that speaks only about single observations cannot tell them apart — yet \(\Cov(X_s, X_t) = 0\) for \(s \neq t\) under A and \(\Cov(X_s, X_t) = 1\) for every pair under B (verified in Exercise 1.1). One realization of each:

# Simulate two processes with identical marginals and different dependence.
set.seed(1)
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), bg = "gray92")
plot(rnorm(200), type = "l", col = 4, ylim = c(-3, 3),
     xlab = "t", ylab = expression(x[t]), main = "A: i.i.d. N(0,1)")
plot(rep(rnorm(1), 200), type = "l", col = 4, ylim = c(-3, 3),
     xlab = "t", ylab = expression(x[t]), main = "B: one draw, repeated")

Figure 3: One realization each of Processes A and B of Example 1.2. Every observation in both panels is a draw from N(0,1). Simulated data.

The marginal distribution does not determine the process. A time series model therefore cannot be a collection of distributions for one data point: it must specify the relationship between all past and future values. ◊

So how do we derive a model for time series data? For \(T\) observations with values in \(\R^d\), the random object is the complete random data matrix \[\mathbf X_{1:T} := \begin{pmatrix} X_1^\top\\X_2^\top\\\vdots\\X_T^\top \end{pmatrix} \in\R^{T\times d},\] of which the observed matrix, with lowercase rows \(x_1^\top,\ldots,x_T^\top\), is one realization. It too is specified by a CDF, \[F_{1:T}(c_1,\ldots,c_T) :=\P\left(X_1\le c_1,\ldots,X_T\le c_T\right), \qquad c_1,\ldots,c_T\in\R^d,\] a single function of \(Td\) real arguments, describing how all \(T\) observations vary together. Just specifying a CDF for the \(T\) observations is typically insufficient. A time series model is asked for more: From the observed path we want to describe the dependence along it, forecast \(X_{T+1},X_{T+2},\ldots\), reconstruct values missing from the record or preceding it, and attach uncertainty to each answer (the task list of the section “What We Want From a Time Series”). Each of these concerns how the coordinates relate to one another — observed to observed, observed to unobserved. This brings us to the following definitions.

The Formal Objects

Definition 1.3 (Time series). A time series (or stochastic process) is a family of random vectors \((X_t)_{t \in \Z}\) with values in \(\R^d\), one for each time point, specified through its finite-dimensional distributions: for every dimension \(k\) and every choice of time points \(t_1 < t_2 < \cdots < t_k\), the joint distribution of the random vector \[(X_{t_1}, X_{t_2}, \dots, X_{t_k}) \in (\R^d)^k,\] that is, the function \[F_{t_1, \dots, t_k}(c_1, \dots, c_k) := \P\left(X_{t_1} \le c_1, \dots, X_{t_k} \le c_k\right), \qquad c_1, \dots, c_k \in \R^d ,\] the inequalities read componentwise when \(d > 1\).

In words: to specify a process is to specify the joint distribution of every finite window of it. Four remarks:

• The index set. We take \(t \in \Z\), a doubly infinite sequence, as the default; a dataset occupies the stretch \(t = 1, \dots, T\). (Oddly, a time series is thus a sequence, not a series.) Why \(\Z\) is taken up below.
• Windows are random vectors. For each fixed choice of times, \((X_{t_1}, \dots, X_{t_k})\) is an ordinary random vector — a point of \(\R^{dk}\), specified by one CDF: exactly the \(F_{t_1, \dots, t_k}\) the definition lists.
• Consistency. The distributions must agree with one another: integrating the last coordinate out of \(F_{t_1,\dots,t_k}\) returns \(F_{t_1,\dots,t_{k-1}}\). A family with this property is called consistent, and consistent families are all we will ever write down.
• The value space. The definition costs nothing extra at general \(d\), and several series studied jointly form one vector-valued series — hare and lynx form a single series in \(\R^2\); images and curves per moment stretch the same idea. The theory is developed at \(d = 1\), exceptions flagged; the letter \(d\) is reserved for this dimension.

Remark 1.4 (♠ Does a consistent family define a process?). Does a consistent family of finite-dimensional distributions determine a single random object \((X_t)_{t\in\Z}\)? It does — Kolmogorov’s extension theorem, which justifies the phrase “specified through” in Definition 1.3. The proof needs measure theory this course does not assume; we take the theorem on faith (Brockwell and Davis 1991, Ch. 1). Nothing later depends on the construction, only on the finite-dimensional distributions.

The object demanded at the end of Example 1.2 can now be named.

Definition 1.5 (Time series model). A time series model is a collection \(\cP\) of time series in the sense of Definition 1.3 — equivalently, a collection of consistent families of finite-dimensional distributions — each member regarded as a candidate mechanism for the observed series.

The i.i.d. models of the earlier courses are the special case in which every member’s family is built by multiplication from a single marginal. In general each member carries the joint law of every window — the relationships Example 1.2 showed the marginals cannot see. How a member of \(\cP\) can be written down, given that a family of distributions in every dimension is unwieldy, is answered below: through generating equations.

The Index Set

Why index by \(\Z\)? The displays in Definition 1.3 never use the ordering of the time points: the family makes sense for random variables indexed by an arbitrary set. (Indexing by \(\Z^2\) gives a random field, the object of spatial statistics.) The interpretation of the index as time enters through structure that \(\Z\) carries, and each piece of structure buys a concept.

The total order makes ‘past’ and ‘future’ meaningful: at each \(t\) the remaining observations split into before and after, and prediction — saying something about \(X_{t+1}\) from the variables up to \(t\) — becomes a well-posed question. On \(\Z^2\) no such split exists.

Mathematically we obtain an appealing property: closure under the shifts \(t \mapsto t + h\) is what lets the course’s central regularity assumption, stationarity — invariance of the finite-dimensional distributions under all shifts — be written at all. A finite index set admits no such shift, one reason the model extends beyond any finite record.

Against \(\N\), the advantage of \(\Z\) is an infinite past at every time point: constructions expressing today as an accumulation of all past shocks have no starting boundary, and the generality is free — a stationary process indexed by \(\N\) always extends to \(\Z\) (a fact we shall take on faith). Against \(\R\): this would correspond to continuous time modeling. In practice, most data arrives sampled. Furthermore, continuous time models require heavier mathematical machinery for a satisfactory treatment.

One Process, Many Realizations

A process is a random object, and each draw of it is an entire path — a realization, one value per time point. The data are something else: numbers recorded in the world, carrying no probability of their own. Modeling is the act of treating the observed series as if it were part of one realization of a candidate process from the model, and every probability statement this course makes about data is made under that stance. The notation carries this stance: capital \(X_t\) for the random variable, lowercase \(x_t\) for the observed value, so that \(x_1, \dots, x_T\) is read as a realization of the window \((X_1, \dots, X_T)\). The sample size is written \(T\) throughout, observations \(t = 1, \dots, T\); the letter \(i\) stay reserved for cross-sectional. (One clash: in R code, T abbreviates TRUE, so the code in these workbooks writes n where the text writes \(T\).) The figure shows three realizations of one process.

# Simulate three realizations of the same moving-average process.
set.seed(542)
par(mfrow = c(3, 1), mar = c(2, 4, 1.5, 1), bg = "gray92")
for (i in 1:3) {
  z <- rnorm(202)
  v <- stats::filter(z, rep(1/3, 3))[2:201]
  plot(v, type = "l", col = 4, ylim = c(-1.6, 1.6),
       xlab = "", ylab = expression(x[t]), main = paste("Realization", i))
}

Figure 4: Three realizations of the same process (a three-point moving average of Gaussian noise, defined in Example 1.8). A data analyst receives one panel. Simulated data — shown here and throughout on a gray background, where figures with real data have a white one.

The panels differ point by point, yet they look alike: comparable typical size, comparable speed of variation. Which features of one path reflect the law that produced it — and why a single realization displays them at all — begins with stationarity and mean ergodicity in Workbook 2, then becomes the central question of Workbook 5, Learning Covariances from One Path. The program meanwhile is to learn properties of the process while observing part of one realization. Hence the insistence on notation: Shumway and Stoffer (2025) write \(x_t\) for both process and realization, but many statements in this course are true for one reading and false for the other.

Finitely Observed

The process has infinitely many coordinates; the dataset has \(T\). The observed series occupies the window \(\{1, \dots, T\} \subset \Z\), and the random vector behind it, \((X_1, \dots, X_T)\), has law \(F_{1, \dots, T}\) — one member of the specifying family. A model in the sense of Definition 1.5 induces a model for the data by evaluating each candidate law at the window; for irregularly spaced or missing observations, the relevant law is \(F_{t_1, \dots, t_k}\) at the observed times.

Why carry an infinite model for a finite window? First, no generality is lost: every law for the window arises from some process. Second, the questions of interest live outside the window in both directions — a forecast concerns \(X_{T+1}\), and the past behind an archive’s first page concerns \(X_0, X_{-1}, \dots\); only a process-level model can state either. Third, the window law says nothing about the coordinates beyond it: two process laws can agree on the window and disagree completely about the future, so whatever links the observed stretch to the unobserved one must be assumed — and the assumptions that do this work, stationarity first among them, are process-level statements, as is the course’s asymptotic vocabulary (\(T \to \infty\), ergodicity, stability).

Remark 1.6 (One realization). In the i.i.d. case, one observation of \(T\) variables is really \(T\) observations of one distribution, and classical inference rests on that reduction. Here the dataset is one observation of the single random vector \((X_1, \dots, X_T)\): without further assumptions, one draw from an arbitrary unknown distribution on \(\R^T\), from which nothing can be learned. The subject of time series exists because regularity across time — the same mechanism operating throughout — can be assumed, letting one long realization do the work of many observations. Workbook 2 makes stationarity and mean ergodicity precise; Workbook 5 develops process ergodicity and the broader question of what one path can reveal.

Models from Generating Equations

Definition 1.5 asks for joint distributions in every finite dimension. We rarely write them down one by one. Instead, we introduce an auxiliary process \((Z_t)\), called a driving sequence, and describe \((X_t)\) through an equation such as \[X_t=g_t(Z_t,X_{t-1},X_{t-2},\ldots), \qquad t\in\Z.\]

The displayed equation is only one part of the specification. We must also state the process law of \((Z_t)\), including any independence assumptions, and the parameter restrictions. If the equation contains feedback through earlier \(X\)’s, we must further say which solution on \(\Z\) is intended; for a finite forward simulation, an initial state plays that role. A complete driving law can therefore accompany an equation that has no solution or several solutions.

Sometimes we impose only moment conditions and relations to the past. The resulting model then contains every process law consistent with those restrictions. Such a partial specification can determine means and covariances, and a conditional mean restriction can support least squares, but it does not supply a likelihood. One can instead specify conditional CDFs directly, but those CDFs still need a compatible marginal law on \(\Z\), or an initial-block law for a finite forward model. The workbook on solutions of generating equations formalizes both routes.

White Noise as a Baseline

To see what temporal structure an equation creates, we first need a driving sequence with no linear dependence across time.

Definition 1.7 (White noise). A process \((Z_t)_{t \in \Z}\) is white noise with variance \(\sigma^2\), written \((Z_t) \sim \mathrm{WN}(0, \sigma^2)\), if \[\E[Z_t] = 0, \qquad \Var(Z_t) = \sigma^2, \qquad \Cov(Z_s, Z_t) = 0 \ \text{ for } s \neq t .\]

The notation \(\mathrm{WN}(0,\sigma^2)\) is deliberately incomplete: it fixes means, variances, and covariances, but neither the marginal distribution nor independence. This is the building block used in the first model equations below. The stronger statement \[Z_t\overset{\mathrm{iid}}{\sim}N(0,\sigma^2), \qquad t\in\Z.\] selects one member of the white-noise class and fixes its full process law; we call it Gaussian white noise.

Dependence in the Level

White noise has zero covariance between distinct times. Reusing the same driving variables in neighboring observations creates nonzero covariance over a finite range.

Example 1.8 (Moving average). Let \((Z_t)\sim\mathrm{WN}(0,\sigma^2)\) and define \[V_t := \tfrac{1}{3}\left(Z_{t-1} + Z_t + Z_{t+1}\right).\] Every white-noise law produces a process \((V_t)\) through this explicit equation. Its unspecified higher-order features can change with the law of \((Z_t)\), but its covariances cannot: consecutive values share two terms, lag-two values share one, and values three or more time points apart are uncorrelated (Exercise 1.1). For the figure only, we set \(\sigma^2=1\) and select Gaussian white-noise for \((Z_t)\) (you can try it with another WN process as an exercise). Averaging makes that simulated series vary more slowly than the noise from which it is constructed. ◊

# Compare Gaussian noise with its three-point moving average.
set.seed(3)
par(mfrow = c(2, 1), mar = c(2.5, 4, 2, 1), bg = "gray92")
z <- rnorm(252)
v <- stats::filter(z, rep(1/3, 3))[2:251]
plot(z[2:251], type = "l", col = 4, ylim = c(-3, 3),
     ylab = expression(z[t]), main = "I.i.d. Gaussian noise")
plot(v, type = "l", col = 4, ylim = c(-3, 3),
     ylab = expression(v[t]), main = "Three-point moving average of the same noise")

Figure 5: I.i.d. Gaussian noise, and the three-point moving average of the same draws. Simulated data.

The moving average covariance vanishes once two observations no longer share driving variables. Feedback can carry a second-order effect further forward.

Example 1.9 (Autoregression, and the random walk). Fix \(\phi\in\R\) and \(\sigma^2>0\), let \((Z_t)\sim\mathrm{WN}(0,\sigma^2)\), and consider \[X_t = \phi X_{t-1}+Z_t, \qquad t\in\Z.\] The white-noise condition partially specifies the driving sequence \((Z_t)\), and the feedback equation leaves a separate question: for a given \(\phi \in \R\) and an admissible white-noise law, which processes on \(\Z\) satisfy it? A later workbook on solutions of generating equations answers that question. The figure created below simulates according to the generating equation for finitely many steps: it takes \(X_0=0\), sets \(\sigma^2=1\), selects Gaussian white noise, and runs forward.

At \(\phi=1\), the equation becomes \[X_t=X_{t-1}+Z_t, \qquad X_t=Z_1+\cdots+Z_t \quad (t\geq 1)\] from the stated zero start. This so called random walk has mean \(0\) and variance \(t\sigma^2\); each driving variable remains in every later value (Exercise 1.2). Adding a constant \(\delta\) to each step gives the random walk with drift, \(X_t=\delta t+\sum_{j=1}^t Z_j\). The temperature record of Example 1.1 poses the later modeling question of whether an apparent trend is a deterministic change or accumulated noise. ◊

# Compare a stable recursion and a random walk driven by the same noise.
set.seed(542)
z <- rnorm(200)
par(mfrow = c(2, 1), mar = c(2.5, 4, 2, 1), bg = "gray92")
plot(as.numeric(stats::filter(z, 0.9, method = "recursive")), type = "l", col = 4,
     ylab = expression(x[t]), main = expression(phi == 0.9))
plot(cumsum(z), type = "l", col = 4,
     ylab = expression(x[t]), main = expression(phi == 1 ~ "(random walk)"))

Figure 6: Two finite recursions from X₀ = 0 driven by the same i.i.d. N(0,1) draws: φ = 0.9 (top) and φ = 1, the random walk (bottom). Simulated data.

Both preceding constructions use a WN process as driving variable, indexed by every time point. A process can instead be random because a finite set of quantities is drawn once.

Example 1.10 (A deterministic trigonometric series). Fix \(\lambda\in(0,\pi)\), let \(A,B\) be uncorrelated random variables with mean zero and variance \(\sigma^2\), and set \[X_t=A\cos(\lambda t)+B\sin(\lambda t), \qquad t\in\Z.\] These moment conditions determine the mean and covariance of \((X_t)\) but leave the joint distribution of \((A,B)\) open. Every realization is nevertheless a pure sinusoid, and any two consecutive observations determine \(A\), \(B\), and hence the entire path (Exercise 1.5). For the figure, we select independent \(N(0,1)\) coefficients. The process is random, but its future contains no new randomness once two observations are known. The workbook on dependence by frequency returns to these sine-and-cosine coordinates. ◊

# Simulate three sinusoidal paths from random coefficient pairs.
set.seed(542)
par(bg = "gray92")
tt <- 1:72
plot(NULL, xlim = c(1, 72), ylim = c(-3.2, 3.2), xlab = "t", ylab = expression(x[t]),
     main = "Three draws of (A, B), three sinusoids")
for (col in c(4, 2, 1)) {
  ab <- rnorm(2)
  lines(tt, ab[1] * cos(pi * tt / 6) + ab[2] * sin(pi * tt / 6), col = col)
}

Figure 7: Three realizations of the trigonometric process (λ = π/6, period 12): each pair of independent Gaussian coefficients (A, B) fixes one sinusoid. Simulated data.

Dependence in the Variability

The preceding examples place dependence in the level of the series. Financial returns often show a different pattern: little visible dependence in the signed values, but persistence in their magnitudes.

Example 1.11 (Stock returns with drift). From the DJIA closings of Example 1.1, form the daily log returns \(r_t = \log(x_t) - \log(x_{t-1})\), which for small moves approximate percentage changes. At the process level, a simple second-order model is \[R_t=\delta+Z_t, \qquad (Z_t)\sim\mathrm{WN}(0,\sigma^2).\] Equivalently, \(\log X_t=\log X_{t-1}+\delta+Z_t\), so \(\delta\) is the drift in the log index. A drift term makes sense because a stock index can grow on average over time, giving its log level a nonzero average slope and its log returns a nonzero mean. After that mean is removed, the signed returns have little visible serial pattern, but their magnitude changes in clusters: quiet days tend to follow quiet days, and turbulent days tend to follow turbulent days.

# Compute and plot DJIA log returns after estimating their drift.
r <- diff(log(as.numeric(djia$Close)))
delta_hat <- mean(r)
z_hat <- r - delta_hat
round(c("estimated drift" = delta_hat,
        "innovation variance" = var(z_hat),
        "lag-1 innovation covariance" =
          cov(z_hat[-1], z_hat[-length(z_hat)])), 6)

            estimated drift         innovation variance 
                   0.000186                    0.000146 
lag-1 innovation covariance 
                  -0.000015 

plot(index(djia)[-1], r, type = "l", col = 4,
     xlab = "", ylab = "Log return", main = "DJIA daily log returns")
abline(h = delta_hat, lty = 2)

Figure 8: DJIA daily log returns, April 2006–April 2016. The dashed line is the estimated daily drift; variability occurs in clusters. Real data: astsa::djia.

The sample mean \(\hat\delta=\bar r\) estimates the drift. After subtracting it, the code estimates the variance and lag-one covariance of the zero-mean innovations. One lag-one covariance cannot establish the white-noise conditions at every lag. Moreover, white noise is only a second-order specification, so uncorrelated drift-adjusted returns can still have dependent magnitudes. Whether a sample covariance is distinguishable from zero requires the inference-under-dependence tools developed later. ◊

One way to represent the volatility clustering of the DJIA daily log returns is to let the scale of the next observation depend on recent squared drift-adjusted returns.

Example 1.12 (GARCH). Fix \(\delta\in\R\), \(\omega>0\), and \(\alpha,\beta\geq0\). Let \((Z_t)\) be i.i.d. with \[\E[Z_t]=0, \qquad \Var(Z_t)=1,\] its common distribution otherwise unspecified. Require \(Z_t\) to be independent of \((X_s,\sigma_s^2)_{s<t}\), and consider \[X_t = \delta+\sigma_t Z_t, \qquad \sigma_t^2 = \omega + \alpha (X_{t-1}-\delta)^2 + \beta \sigma_{t-1}^2.\] The unit variance and the independence condition imply (check this!) \[\E[X_t\mid (X_s,\sigma_s^2)_{s<t}]=\delta, \qquad \Var(X_t\mid (X_s,\sigma_s^2)_{s<t})=\sigma_t^2.\] Thus, the second equation gives “feedback” in the conditional variance: there is dependence on the volatility in the past. For the figure below we set \((\delta,\omega,\alpha,\beta)=(0.02,0.05,0.15,0.8)\) and \(\sigma_1^2=1\), then select the standard normal distribution for the \(Z_t\)’s. These extra choices completely specify the finite forward simulation distribution. ◊

# Simulate a GARCH path with recursively changing conditional variance.
set.seed(542)
n <- 400
delta <- 0.02; omega <- 0.05; alpha <- 0.15; beta <- 0.8
z <- rnorm(n); s2 <- numeric(n); x <- numeric(n)
s2[1] <- 1; x[1] <- delta + sqrt(s2[1]) * z[1]
for (t in 2:n) {
  s2[t] <- omega + alpha * (x[t-1] - delta)^2 + beta * s2[t-1]
  x[t] <- delta + sqrt(s2[t]) * z[t]
}
par(bg = "gray92")
plot(x, type = "l", col = 4, ylab = expression(x[t]), main = "GARCH(1,1)")
abline(h = delta, lty = 2)

Figure 9: A finite GARCH(1,1) recursion with drift δ = 0.02 and i.i.d. Gaussian innovations (ω = 0.05, α = 0.15, β = 0.8, σ₁² = 1): calm and turbulent stretches alternate. Simulated data. Are you satisfied with how this compares to the DJIA daily log returns? Try changing the parameter values until you get a more comparable picture.

In a GARCH model, the variance state is determined by its previous value and the preceding drift-adjusted observation. Stochastic volatility gives that state its own driving noise.

Example 1.13 (Stochastic volatility). Fix \(\delta,\mu\in\R\), \(|\phi|<1\), and \(\sigma_\eta^2>0\). Let \((Z_t)\) and \((\eta_t)\) be independent i.i.d. sequences satisfying \[\E[Z_t]=\E[\eta_t]=0, \qquad \Var(Z_t)=1, \qquad \Var(\eta_t)=\sigma_\eta^2,\] with both common marginal distributions otherwise unspecified. Let \((h_t)\) satisfy \[h_t = \mu + \phi(h_{t-1}-\mu)+\eta_t\]

(as before; in later workbooks we will worry about verifying whether a time series \((h_t)\) exists that satisfies the above equation — for now we take this on faith).

Define the observed process by \[X_t=\delta+e^{h_t/2}Z_t.\] The assumptions define a family of process laws. In every member, \(h_t\) is an unobserved, or latent, log-variance process and \(\E[X_t\mid h_t]=\delta\) and \(\Var(X_t\mid h_t)=e^{h_t}\). For the figure we use \(\delta=0.02\), \(\mu=-1\), \(\phi=0.95\), and \(\sigma_\eta=0.3\), select Gaussian distributions for both driving sequences \((Z_t)\) and \((\eta_t)\), and draw \(h_1\) from the resulting Gaussian stationary marginal distribution. The state-space workbooks develop models with this observed/latent structure. ◊

# Simulate returns driven by a latent autoregressive log variance.
set.seed(542)
n <- 400
delta <- 0.02; mu <- -1; phi <- 0.95; sigma_eta <- 0.3
h <- numeric(n)
h[1] <- rnorm(1, mean = mu, sd = sigma_eta / sqrt(1 - phi^2))
eta <- rnorm(n, sd = sigma_eta)
z <- rnorm(n)
for (t in 2:n) h[t] <- mu + phi * (h[t-1] - mu) + eta[t]
par(bg = "gray92")
plot(delta + exp(h/2) * z, type = "l", col = 4,
     ylab = expression(x[t]), main = "Stochastic volatility")
abline(h = delta, lty = 2)

Figure 10: A stationary stochastic-volatility simulation with drift δ = 0.02 (μ = −1, φ = 0.95, ση = 0.3): the visible series is driven by an unobserved log-variance process. Simulated data.

The model statements above retain only the distributional assumptions needed for their claims: white-noise moments for the moving average and autoregression, coefficient moments for the trigonometric process, innovation independence for conditional variances, and a selected stationary solution for the latent recursion. The figures choose Gaussian members of those model classes so that concrete paths can be drawn, but we could alternatively have chosen them e.g. Laplace or even uniform on \([-1,1]\) — try either of those choices for the WN process and see how it changes the observed series in the examples.

What We Want From a Time Series

Facing series like those of Example 1.1, we typically want some of the following.

1. Description. Summarize the dependence: how strongly is today related to yesterday, to last month? Which regularities — trend, seasonality, cycles — are present?
2. Modeling. Propose a mechanism, in the sense of Definition 1.5 and its generating equations, that could plausibly have generated the observed realization; estimate its parameters; check its fit.
3. Forecasting. Say something about \(X_{T+1}, X_{T+2}, \dots\) given the observed \(x_1, \dots, x_T\), with uncertainty attached — a question only a process-level model can pose. In case we wish to infer the unobserved past, say \(X_{-1}\) or \(X_0\), we call it “backcasting”, but the mechanism is the same.
4. Inference under dependence. Answer questions like “is the warming trend real?” with tests and standard errors that remain valid under correlation; the independence assumption enters such procedures silently, through the standard errors.
5. Validation. Judge models and forecasts on data they have not seen; when data share a timeline, the future must not leak into the training past.

Two of the five — inference and forecasting — can be made precise with this workbook’s framework alone; doing so now fixes what success will mean later.

Estimating a Parameter

A parametric model is a model in the sense of Definition 1.5 whose members are indexed by a parameter set: \(\{P_\theta : \theta \in \Theta\}\). The estimand is \(\theta\), or a function of it; an estimator is a function of the data window, \[\hat\theta_T = f(X_1, \dots, X_T) ;\] and a level \(1-\alpha\) confidence set is a set-valued function of the data \(C_T \equiv C_T(X_1,\dots,X_T)\) (i.e. a random subset of \(\Theta\), depending on the observed data) with \[P_\theta\left(\theta \in C_T\right) \ge 1 - \alpha \qquad \text{for all } \theta \in \Theta .\] These are the typical i.i.d. setting definitions unchanged; the novelty is that \(P_\theta\) is now the law of a process, so the coverage probability is computed under dependence.

Example 1.14 (Hare and lynx). The pelt records of Example 1.1 move in linked cycles of roughly ten years. Predator-prey feedback suggests two directions: a larger lynx population may reduce the next hare population, while a larger hare population may support the next lynx population. The pelt records are only proxies for those populations, but they motivate a paired model.

We model the original hare and lynx pelt records by positive-valued processes \((X_{H,t})\) and \((X_{L,t})\). The label \(H\) or \(L\) identifies the species, while \(t\) identifies the year. Define their log transforms by \[Y_{H,t}:=\log X_{H,t}, \qquad Y_{L,t}:=\log X_{L,t}.\] The \(X\) variables are measured in thousands of pelts; the \(Y\) variables are their logs. We consider candidate process laws satisfying the one-lag restrictions \[\begin{aligned} Y_{H,t+1} &=a_H+b_{HH}Y_{H,t}+b_{HL}Y_{L,t}+Z_{H,t+1},\\ Y_{L,t+1} &=a_L+b_{LH}Y_{H,t}+b_{LL}Y_{L,t}+Z_{L,t+1}. \end{aligned}\] Write \(Z_{t+1}:=(Z_{H,t+1},Z_{L,t+1})^\top\) and bundle the history of both processes through year \(t\) as \[\cH_t:=\bigl(\ldots, Y_{H,t-1},Y_{L,t-1},Y_{H,t},Y_{L,t}\bigr).\] We assume \[\E[Z_{t+1}\mid\cH_t]=0, \qquad \Cov(Z_{t+1}\mid\cH_t)=\Sigma,\] where \(\Sigma\) is a fixed \(2\times 2\) positive semidefinite matrix. The first restriction makes the two equations conditional-mean equations: after the current pair is included, earlier records do not change the predicted means. The second restriction keeps the conditional innovation covariance fixed over time.

Take \(b_{LH}\) as the estimand. The first subscript names the response, and the second names the lagged predictor: \(b_{LH}\) multiplies the current log hare record in the equation for the next log lynx record. Informally, it describes how the one-step lynx prediction changes with the current hare record, holding the current lynx record fixed.

For the exact interpretation, fix current records \(x_H,x_L>0\). Replacing \(x_H\) by \(c x_H\), for \(c>0\), changes the conditional mean of the next log lynx record by \[\begin{aligned} &\{a_L+b_{LH}\log(c x_H)+b_{LL}\log x_L\} -\{a_L+b_{LH}\log x_H+b_{LL}\log x_L\}\\ &\qquad=b_{LH}\log c. \end{aligned}\] If \(b_{LH}=0\), the current hare record adds no one-step conditional-mean information once the current lynx record is included. However, if the expected predator–prey relationship is present in these records, we would expect \(b_{LH}>0\): holding the current lynx record fixed, a larger current hare record predicts a larger lynx record next year.

The restriction \(\E[Z_{L,t+1}\mid\cH_t]=0\) supplies the population regression condition used by least squares. However, the standard errors provided assume independence of the errors, meaning the standard confidence interval will only be a nominal one. Calling the interval below nominally 95% means that the ordinary OLS confidence-interval formula would attain 95% coverage under its usual assumptions (which includes independent errors). Here independence is doubtful; the residual calculation below gives evidence against it. We should therefore not expect the interval’s actual coverage to equal 95%.

# Fit the paired models and summarize the lynx equation and residual dependence.
X_H <- as.numeric(Hare)
X_L <- as.numeric(Lynx)
Y_H <- log(X_H)
Y_L <- log(X_L)
n <- length(X_H)  # R uses n here; the text uses T for the window length

lagged <- data.frame(
  Y_H_next = Y_H[-1], Y_L_next = Y_L[-1],
  Y_H_now = Y_H[-n], Y_L_now = Y_L[-n]
)
fit_H <- lm(Y_H_next ~ Y_H_now + Y_L_now, data = lagged)
fit_L <- lm(Y_L_next ~ Y_H_now + Y_L_now, data = lagged)

round(coef(fit_H), 3)

(Intercept)     Y_H_now     Y_L_now 
      2.002       0.740      -0.372 

round(coef(fit_L), 3)

(Intercept)     Y_H_now     Y_L_now 
      0.230       0.206       0.700 

b_LH_hat <- unname(coef(fit_L)["Y_H_now"])
b_LH_ci <- unname(confint(fit_L, "Y_H_now"))
round(b_LH_ci, 3)  # ordinary 95% regression interval

      [,1]  [,2]
[1,] 0.113 0.298

resid_acf1 <- c(
  hare = acf(resid(fit_H), plot = FALSE)$acf[2],
  lynx = acf(resid(fit_L), plot = FALSE)$acf[2]
)
round(resid_acf1, 2)

hare lynx 
0.03 0.52 

Least squares gives \(\hat b_{LH}=0.206\). For a doubling of the current hare record, exponentiating the two fitted log-scale predictions gives the ratio \[2^{\hat b_{LH}}=1.153.\] This fit says roughly that, holding the current lynx record fixed, a doubling of the current hare record therefore multiplies the fitted one-step lynx prediction by \(1.153\).

The residual correlations between \(t\) and \(t+1\) are \(0.03\) in the hare equation and \(0.52\) in the lynx equation. Under the conditional-mean restriction, each innovation is uncorrelated with earlier innovations. The lynx residual correlation supports the concern about independent regression errors. ◊

Predicting the Next Value

Once \(x_1,\ldots,x_T\) have been observed, a one-step forecast is a number \(\hat x_{T+1\mid T}\) computed using only this observed history. Its target is the future value \(X_{T+1}\), which remains random under the model. Under squared-error loss, its quality is measured by the conditional mean squared prediction error \[\E\!\left[ (X_{T+1}-\hat x_{T+1\mid T})^2 \mathrel{\Big|} X_1=x_1,\ldots,X_T=x_T \right].\]

A one-step prediction interval is a set, computed from the observed window, intended to contain a future observation. For example, under a known model we might require a level \(1-\alpha\) interval to satisfy \[\P\!\left(X_{T+1} \in C_T \mid X_1, \dots, X_T\right) \ge 1 - \alpha.\] Workbook 3 develops prediction under a known process law; Workbooks 15 and 16 return to interval construction after fitting and to empirical coverage across forecast origins.

A prediction interval concerns a future random variable and must reflect whatever uncertainty remains about that value under the model. It need not collapse even when the model is known, although it can collapse when the observed past determines the future. This is unlike the typical confidence interval, which concerns a fixed unknown number and should shrink as information accumulates.

Example 1.15 (Forecasting a random walk). For the driftless random walk of Example 1.9 with Gaussian noise, \[X_{T+h}=X_T+\sum_{j=1}^h Z_{T+j},\] so the conditional distribution of \(X_{T+h}\) given the window is \(N(x_T, h\sigma^2)\). The point forecast is therefore \(\hat x_{T+h \mid T} = x_T\) at every horizon, and \[x_T \pm 1.96\, \sqrt{h}\, \sigma\] is an exact 95% prediction interval for each fixed horizon \(h\). Across horizons, these intervals form a pointwise 95% prediction fan, not a simultaneous 95% region for the entire future path. The point forecast is flat — the walk’s best guess for every future is wherever it stands today — and the uncertainty grows like \(\sqrt{h}\) without bound: both consequences of the permanent memory computed in Example 1.9, and both different for the mean-reverting processes studied in Workbook 2. ◊

# Simulate a random walk and plot its pointwise prediction intervals.
set.seed(11)
z <- rnorm(120)
x <- cumsum(z)
h <- 1:20
par(bg = "gray92")
plot(1:100, x[1:100], type = "l", col = 4, xlim = c(1, 120),
     ylim = range(c(x, x[100] + 1.96 * sqrt(20), x[100] - 1.96 * sqrt(20))),
     xlab = "t", ylab = expression(x[t]), main = "Random walk: pointwise 95% prediction fan")
lines(100:120, x[100:120], col = 2)
lines(100 + h, rep(x[100], 20), lty = 3)
lines(100 + h, x[100] + 1.96 * sqrt(h), lty = 2)
lines(100 + h, x[100] - 1.96 * sqrt(h), lty = 2)
legend("topleft", legend = c("Observed", "Continuation", "Point forecast", "Pointwise 95%"),
       col = c(4, 2, 1, 1), lty = c(1, 1, 3, 2), bty = "n")

Figure 11: A Gaussian random walk observed to T = 100 (blue), its actual continuation (red), the flat point forecast (dotted), and pointwise 95% prediction intervals (dashed). The fan is not a simultaneous 95% region for the entire future path. Simulated data.

Looking Ahead

Workbook 2, Solutions of Generating Equations, asks what it means to solve an equation, when a solution exists, and what happens when many do; it then introduces stationarity, the autocovariance function, and mean ergodicity. Workbook 5, Learning Covariances from One Path, studies which averages can recover the covariance features of a process. Workbook 6, Estimating Uncertainty from One Path, estimates the uncertainty of the sample mean and compares its precision with independent sampling. Workbook 7, From Dependence to a Central Limit Theorem, supplies conditions for Gaussian inference.

Exercises

Exercise 1.1 (The marginal is not the process). Three processes, all built from i.i.d. standard normal ingredients: Process A is i.i.d. \(N(0,1)\); Process B repeats a single draw, \(X_t := X_0 \sim N(0,1)\) for all \(t\); Process C is the scaled moving average \(W_t := \sqrt{3}\, V_t\) with \(V_t\) as in Example 1.8 (noise variance \(\sigma^2 = 1\)).

1. Verify that under all three processes, every observation has exactly the \(N(0,1)\) distribution. Hint: for C, a linear combination of independent normals is normal; compute its variance.

2. Compute \(\Cov(X_s, X_t)\) for all \(s \neq t\) under A and under B, verifying the claims of Example 1.2.

3. For Process C, compute \(\Cov(W_t, W_{t+h})\) for \(h = 1, 2, 3\), and conclude that the correlation between neighboring values is \(2/3\) while values three or more time points apart are independent. Hint: expand and count shared noise terms; only matching \(Z\)’s contribute.

4. Simulate a realization of length \(T = 300\) from each process (set the seed 542) and plot the three side by side on a common vertical scale.

# Set up simulations of three processes with common marginals.
set.seed(542)
# Process A:
# Process B:
# Process C:   (stats::filter(z, rep(1/3, 3)) computes the moving average)

5. Interpret in complete sentences. All three processes in (d) have \(N(0,1)\) marginals. State, for each process, what it asserts about the relationship between today and tomorrow, and what this implies about how much a histogram from one realization can reveal about a time series.

Exercise 1.2 (The random walk never settles). Let \((Z_t)\) be i.i.d. noise with mean zero and variance \(\sigma^2\), and let \(X_t = \delta t + \sum_{j=1}^{t} Z_j\) be the random walk with drift of Example 1.9.

1. Show that \(\E[X_t] = \delta t\) and \(\Var(X_t) = t\,\sigma^2\).

2. Show that \(\Cov(X_s, X_t) = \min(s, t)\,\sigma^2\); the drift does not enter. Hint: for \(s \le t\), split \(X_t\) into \(X_s\) plus the increment \(X_t - X_s\), and note the increment is built from fresh noise, uncorrelated with everything in \(X_s\).

3. Conclude that \(\Corr(X_t, X_{2t}) = 1/\sqrt{2}\), whatever the value of \(t\). Contrast with Process C of Exercise 1.1, where the correlation is zero at lag 3 and beyond: in what sense does the walk not forget?

4. Simulate 10,000 replications of a driftless walk with \(\sigma = 1\) (seed 542) and estimate \(\Var(X_{100})\) and \(\Var(X_{400})\) empirically. Report the two estimates next to the values predicted by (a) in a small table.

# Set up repeated simulations of the random-walk variance.
set.seed(542)
R <- 10000
# hint: a matrix of noise with cumsum applied columnwise

5. Interpret. By (a), the uncertainty about the level of the walk grows without bound, although each increment is perfectly well behaved. Return to the temperature record of Example 1.1 and describe, in two or three sentences, what is at stake in deciding whether such a series is trend plus stable fluctuation or a drifting random walk.

Exercise 1.3 (One path or many?). Averages can be taken in two directions: along time, over one realization (\(T^{-1}\sum_{t=1}^T x_t\)), or across realizations, at a fixed time (the average of \(x_{50}\) over many independent runs). In the i.i.d. case the law of large numbers makes the first behave like the second. Under dependence this equivalence becomes a real question.

1. For Process C of Exercise 1.1 (the scaled moving average): simulate one long realization of length \(T = 10{,}000\) and compute its time average. Then simulate 10,000 independent short realizations, record the value at \(t = 50\) of each, and compute their average. Use seed 542 for each experiment.

2. Repeat both computations for Process B (one draw, repeated).

3. Report the four numbers in a 2 × 2 table: process (C, B) against direction of averaging (along time, across realizations).

# Set up time-average and ensemble-average simulations.
set.seed(542)
# (a) one path of length 10000; then 10000 paths, keep t = 50
# (b) same for the constant process

4. Interpret. Both processes have \(\E[X_t] = 0\), and averaging across realizations recovers it both times; averaging along time recovers it only once. Explain, in terms of what a single path of each process explores, why the time average of Process B fails, and state informally what property of a process makes time averages work for Process C. Workbook 2 names this target-specific property mean ergodicity; Workbook 5 develops full process ergodicity.

Exercise 1.4 (Reading real series). Consider the same data as in Example 1.1: astsa::gtemp_land, astsa::gtemp_ocean (annual global mean temperature anomalies, °C relative to 1991–2020, 1850–2023) and astsa::djia (DJIA daily closings, 2518 trading days, April 2006 – April 2016; an xts object, so load xts).

1. Plot the two temperature series on one set of axes, with a legend and labeled axes including units.

# Set up a comparison of land and ocean temperature anomalies.
plot(gtemp_land, col = 4, ylab = "Anomaly (°C vs 1991–2020)")
# add the ocean series and a legend

2. Compute the DJIA daily log returns, estimate their drift by \(\hat\delta=\bar r\), and add the estimated drift to the return plot as a horizontal line. Then form the drift-adjusted returns \(\hat z_t=r_t-\hat\delta\) and plot \(\hat z_t^2\) against date.

# Compute drift-adjusted returns for plotting their variability.
r <- diff(log(as.numeric(djia$Close)))
delta_hat <- mean(r)
z_hat <- r - delta_hat
# plot r with a line at delta_hat, then plot z_hat^2

3. For each of the two series (temperature; returns), describe in complete sentences three features visible in your plots — trend, seasonality, cycles, changing variability, persistence, outliers — and for each feature name the task from “What We Want From a Time Series” that it calls for.

4. Shuffle the drift-adjusted returns (sample(z_hat)) and plot the shuffled series next to the real one; do the same for their squares. Two of these four panels look alike and two do not. Which, and why? Hint: shuffling preserves exactly the marginal distribution and destroys exactly the time structure — Example 1.2 in reverse.

5. Interpret (d) in two or three sentences: what does the survival of one pattern and the loss of the other indicate about where the dependence in drift-adjusted stock returns lives? Reconcile this with the distinction between white noise and i.i.d. noise in Definition 1.7 and Example 1.11.

Exercise 1.5 (Destiny in two observations). Use the trigonometric process of Example 1.10, with \(\lambda \in (0, \pi)\) known and \(A,B\) uncorrelated, mean-zero, variance-\(\sigma^2\) random variables.

1. Show that \(\E[X_t] = 0\) and, using the cosine addition formula, that \[\Cov(X_t, X_{t+h}) = \sigma^2 \cos(\lambda h) :\] the covariance depends only on the lag, and it never decays.

2. Show that \((A, B)\) — and hence the entire path, for all \(t \in \Z\) — is determined by any two observations \((x_1, x_2)\). Hint: two linear equations in two unknowns; check the determinant is nonzero for \(\lambda \in (0, \pi)\).

3. Simulate five realizations (\(T = 72\), \(\lambda = \pi/6\), seed 542) and overlay them; then, for one realization, reconstruct \((A, B)\) from its first two values and verify the reconstruction reproduces the path.

# Set up recovery of sinusoid coefficients from two observations.
set.seed(542)
# solve a 2x2 linear system for (A, B) from x[1], x[2]

4. Interpret. The series is random, and by (a) its dependence reaches every lag undiminished; yet by (b) its future is decided by two numbers — prediction error zero. Contrast this with i.i.d. noise, where each future value contains a fresh draw, and place the moving average of Example 1.8 between the two poles. In one closing sentence, explain why series of this kind are called ‘signals’ rather than ‘noise’ and how their sine-and-cosine form anticipates Workbook 10’s frequency-domain coordinates.

Exercise 1.6 (One-step memory). Restrict the generating equation to one lag: \(X_t = g(Z_t, X_{t-1})\), with \((Z_t)\) i.i.d. and \(Z_t\) independent of \((X_s)_{s < t}\).

1. Explain why the conditional distribution of \(X_t\) given the entire past \((X_{t-1}, X_{t-2}, \dots)\) depends on \(X_{t-1}\) only. Hint: given \(X_{t-1} = x\), the value \(X_t = g(Z_t, x)\) is built from \(x\) and fresh noise, and the fresh noise knows nothing of the deeper past.

2. Which equations of the gallery have this one-lag form as written, and which do not? For the GARCH equation of Example 1.12, exhibit a state whose evolution is one-lag. Hint: track the pair \((X_t, \sigma_t^2)\) together.

3. Simulate a nonlinear one-lag chain, \(X_t = \sin(X_{t-1}) + 0.5\, Z_t\) with standard normal noise (\(T = 500\), start \(x_0 = 0\), seed 542), and plot the path.

# Set up simulation of the nonlinear one-lag process.
set.seed(542)
x <- numeric(500)
# x[1] <- sin(0) + 0.5 * rnorm(1); then iterate

4. The property established in (a) — conditionally on the present, the deeper past is irrelevant — has a name, Markov. Workbook 2 formalizes Markov processes before turning to stationarity. State in one sentence why this restriction might make estimation easier.

---

Sources for this workbook: van der Vaart (2010), Time Series, lecture notes, VU Amsterdam, §1.1; Shumway and Stoffer (2025), Time Series Analysis and Its Applications, 5th ed., Springer, Ch. 1; Brockwell and Davis (1991), Time Series: Theory and Methods, 2nd ed., Springer, Ch. 1. Data: R package astsa (D. Stoffer).