16. Time series¶
Almost every method in this book so far assumed independent observations: shuffle the data and nothing changes. Data recorded in time order break that rule. Today's temperature is close to yesterday's, this month's sales echo last year's, the air pressure on Mars rises and falls with its seasons. This chapter asks how to describe and predict such data honestly. You will split a series into trend, season and noise, measure how strongly a series remembers its past with the autocorrelation function, smooth with moving averages, study the simplest model with memory, the AR(1), and prove its variance and correlation formulas, make forecasts with error bars that grow the right way, and end with a real example from Mars: the seasonal pressure cycle first measured by the Viking landers.
In this chapter
You will learn:
- how to describe a series as trend plus seasonal pattern plus noise, additively or multiplicatively
- the autocovariance and autocorrelation functions, how to estimate them, and what "significant" lags mean
- moving averages and exponential smoothing, and how smoothing can create cycles that are not there
- stationarity, the AR(1) model, and proofs of its variance \(\sigma^2/(1-\phi^2)\) and correlations \(\phi^h\)
- how to forecast with an AR(1), why forecast intervals widen and then level off, and how to test forecasts honestly
- how the seasonal pressure cycle on Mars was measured and how a time-series analyst would approach it
You need: Chapter 4 (correlation), Chapter 11 (least squares), and from the Probability book expectation and variance and joint distributions (covariance of sums). Geometric series are used in 16.3 to 16.5 (series). Sines and cosines appear in 16.6.
Levels: 16.1 to 16.3 are intermediate; 16.4 and 16.5 are advanced (full proofs); 16.6 is intermediate, with an advanced part on harmonic regression.
16.1 Trend, season, noise¶
A picture first¶
Here is four years of made-up monthly data, say the number of visitors to a museum, in thousands:
Three things are visible at once. The values drift upwards over the years: a trend. They rise and fall in the same pattern every twelve months: a seasonal component. And around that regular shape they jitter unpredictably: noise, or the irregular component. (The made-up series was built exactly that way: a trend rising 0.8 per month, a sine wave of height 12 with period 12, and normal noise with standard deviation 4.)
The decomposition¶
A time series is a sequence of observations \(x_1, x_2, \dots, x_n\) indexed by time, usually at equal spacing. The classical additive decomposition writes
with a slowly varying trend \(T_t\), a seasonal component \(S_t\) that repeats with period \(s\) (\(S_{t+s} = S_t\), and by convention \(S_1 + \dots + S_s = 0\) so that the season does not secretly carry part of the level), and an irregular component \(E_t\). When the seasonal swings grow in proportion to the level, as they often do for sales or populations, the multiplicative decomposition \(x_t = T_t \times S_t \times E_t\) fits better; taking logarithms turns it into the additive one.
The decomposition is a description, not a law of nature. "Trend" and "season" are names for the parts that change slowly and the parts that repeat, and different methods split the same series differently. Still, it is the first thing to do with any series, because each component answers a different question: is it growing, when in the year is it high, and how unpredictable is it?
A real example in general terms¶
The best-known real example is the record of carbon dioxide in the air measured at Mauna Loa Observatory in Hawaii, begun by Charles David Keeling in 1958 and often called the Keeling curve. It shows a steady upward trend, from burning fossil fuels, and a seasonal cycle of a few parts per million each year: concentrations fall during the northern growing season, when plants take up carbon dioxide, and rise again in autumn and winter. Trend and season have completely different causes, and the decomposition separates them.
Why independence fails¶
In every earlier chapter the observations were independent. In a time series, neighbouring values are related: a month above trend tends to be followed by another above trend. This changes everything that depended on independence. The standard error \(s/\sqrt n\) of a mean, for example, is wrong for correlated data, usually too small, because \(n\) correlated observations carry less information than \(n\) independent ones. Measuring that dependence is the job of the next section.
16.2 Autocorrelation¶
Covariance with yourself¶
Definition. For a series \(X_t\) whose mean \(\mu\) and variance do not change with time, the autocovariance at lag \(h\) is
and the autocorrelation at lag \(h\) is \(\rho(h) = \gamma(h)/\gamma(0)\), the ordinary correlation between the series and itself shifted by \(h\) steps. Always \(\rho(0) = 1\), \(\rho(-h) = \rho(h)\) and \(|\rho(h)| \le 1\) (it is a correlation, Chapter 4). The function \(h \mapsto \rho(h)\) is the autocorrelation function (ACF).
Estimates. From data \(x_1, \dots, x_n\) with mean \(\bar x\),
The divisor is \(n\), not \(n - h\), on purpose: with \(n\) the estimated autocovariances always form a valid covariance structure (the matrix they make is positive semi-definite), and the estimates are less noisy at long lags.
Worked example (made-up data). Ten daily counts: \(3, 5, 4, 6, 5, 3, 2, 4, 5, 3\). The mean is 4, and the deviations are \(-1, 1, 0, 2, 1, -1, -2, 0, 1, -1\).
- \(\hat\gamma(0) = (1 + 1 + 0 + 4 + 1 + 1 + 4 + 0 + 1 + 1)/10 = 14/10 = 1.4\).
- Lag 1 products of neighbours: \((-1)(1) + (1)(0) + (0)(2) + (2)(1) + (1)(-1) + (-1)(-2) + (-2)(0) + (0)(1) + (1)(-1) = 1\), so \(\hat\gamma(1) = 0.1\) and \(r_1 = 0.1/1.4 = 0.071\).
- Lag 2: \((-1)(0) + (1)(2) + (0)(1) + (2)(-1) + (1)(-2) + (-1)(0) + (-2)(1) + (0)(-1) = -4\), so \(\hat\gamma(2) = -0.4\) and \(r_2 = -0.286\).
Is a sample autocorrelation real?¶
If the series is white noise (independent, identically distributed values with finite variance), then for each fixed \(h \ge 1\), \(r_h\) is approximately normal with mean 0 and variance \(1/n\) for large \(n\) (a result due to Bartlett, 1946). So values outside \(\pm 1.96/\sqrt n\) are unusual at a single lag. For the ten counts the bands are \(\pm 0.62\), and neither \(r_1\) nor \(r_2\) is outside them: no evidence of memory, which with ten points is not surprising.
Two cautions. With 20 lags plotted, about one will cross the 5% bands by chance; a joint test such as the Ljung-Box statistic \(n(n+2)\sum_{h=1}^{H} r_h^2/(n-h)\), compared with a chi-square with \(H\) degrees of freedom (Ljung and Box, 1978), is better than staring at single bars. And the bands assume white noise; they are not the right yardstick for lags beyond the first in a series that does have memory.
What the ACF shows¶
- Trend. A trending series has \(r_h\) close to 1 for many lags, decaying very slowly, because values far apart are both above or both below the overall mean. The ACF of a trending series says "trend", and little else.
- Season. A seasonal monthly series has peaks in its ACF at lags 12, 24, ..., and troughs at 6, 18, ....
- Short memory. A stationary series with short memory (16.4) has an ACF that dies away quickly.
Nonsense correlations¶
In 1926 George Udny Yule published "Why do we sometimes get nonsense-correlations between time-series?" (Journal of the Royal Statistical Society). His point, still ignored daily: two series that both trend, or that both wander like random walks, will often show large correlations with each other even when they are completely unrelated, and the usual significance test, which assumes independent observations, will declare the correlation significant far more often than 5% of the time. A correlation between two time series means little until the trends and the autocorrelation have been dealt with, for example by correlating the changes \(x_t - x_{t-1}\) instead of the levels.
16.3 Moving averages and smoothing¶
Centred moving averages¶
To see the trend, average away the season and the noise. A centred moving average of order \(2k+1\) replaces each value by the mean of itself and its \(k\) neighbours on each side:
For monthly data the natural window is 12 months, which is even, so there is no middle month. The standard fix is the 2×12 moving average: average two consecutive 12-month averages, which is the same as giving weight \(1/24\) to the two end months and \(1/12\) to the eleven in between:
The weights add up to \(\frac{1}{12}(\frac12 + 11 + \frac12) = 1\), and the window covers each calendar month with total weight exactly \(1/12\) (the month at both ends shares its weight). So a fixed seasonal pattern \(S_t\) with period 12 and sum zero is averaged to exactly zero, and a linear trend passes through unchanged (the weights are symmetric). That is why the orange line in the picture in 16.1 is smooth and rises at about 0.8 per month. Its price: it is undefined for the first and last six months, exactly where people most want to know the trend.
Seasonal adjustment. Subtracting the trend estimate leaves \(S_t + E_t\); averaging these over all Januaries, all Februaries, and so on (then shifting so that the twelve add to zero) estimates the seasonal component. The series minus its seasonal component is the seasonally adjusted series that official statistics publish. Agencies use refined versions of this idea (the X-13ARIMA-SEATS program of the US Census Bureau, for example), but the logic is the one above.
Exponential smoothing¶
A moving average needs future values. For forecasting we need a smoother that uses only the past. Simple exponential smoothing keeps a running level
and forecasts every future value by the current level. Unrolling the recursion,
a weighted average of the past with weights that shrink geometrically (hence "exponential"). The weights add to 1: \(\alpha\sum_{j=0}^{t-1}(1-\alpha)^j + (1-\alpha)^t = \alpha \cdot \frac{1 - (1-\alpha)^t}{\alpha} + (1-\alpha)^t = 1\), by the geometric series. A large \(\alpha\) follows the data closely; a small one smooths heavily.
Worked example (made-up data). Weekly orders \(12, 15, 11, 14, 16, 13\), with \(\alpha = 0.3\) and \(\ell_1 = 12\):
The forecast for week 7 (and, with this simple method, every later week) is about 13.5 orders.
Robert G. Brown developed exponential smoothing for inventory control in the 1950s; Charles Holt (1957) added a smoothed trend, and Peter Winters (1960) a smoothed seasonal component. The Holt-Winters method is still a strong baseline. These methods were later shown to be the optimal forecasts for certain statistical models (state space models), which gave them error bars; see Hyndman and Athanasopoulos in the further reading.
Smoothing can invent cycles¶
In the 1920s Eugen Slutsky, working in Moscow, showed (his paper appeared in English in Econometrica in 1937 as "The summation of random causes as the source of cyclic processes") that taking moving sums or averages of pure noise produces series with smooth, wave-like swings that look like business cycles. The reason is the ACF: averaging \(k\) neighbours makes values \(k\) apart share terms, so they become correlated, and correlated noise drifts in waves. The lesson is a warning. A smoothed series with regular-looking cycles is not evidence of a cycle in the world; the smoother may have made them.
16.4 AR(1) and stationarity¶
Level: advanced.
Stationarity¶
To use one series to learn about its own future, we need something that does not change over time. The standard assumption is:
Definition. A series \(X_t\) (with finite variances) is weakly stationary if \(E[X_t] = \mu\) is the same for all \(t\), and \(\operatorname{Cov}(X_t, X_{t+h}) = \gamma(h)\) depends only on the lag \(h\), not on \(t\).
A series with a trend is not stationary (its mean changes); a series with a season is not stationary (its mean depends on the month). The decomposition of 16.1 and differencing (\(x_t - x_{t-1}\) removes a linear trend; \(x_t - x_{t-12}\) removes a fixed monthly season) are ways to get from a raw series to something stationary.
White noise is the simplest stationary series: \(\varepsilon_t\) with mean 0, variance \(\sigma^2\), and \(\operatorname{Cov}(\varepsilon_s, \varepsilon_t) = 0\) for \(s \ne t\). It has no memory at all: \(\rho(h) = 0\) for \(h \ne 0\).
The AR(1) model¶
The simplest model with memory says that each value is a fraction of the previous one plus a new shock:
where \(\varepsilon_t\) is white noise and \(\varepsilon_t\) is uncorrelated with all earlier \(X_s\) (\(s < t\)). This is the autoregressive model of order 1, AR(1): the series is regressed on itself. (A non-zero mean is added by writing \(X_t - \mu = \phi(X_{t-1} - \mu) + \varepsilon_t\).) The constant \(\phi\) measures memory: \(\phi = 0\) is white noise, \(\phi\) near 1 is a series with long swings, \(\phi\) negative makes it zigzag.
Here are three series built from the same 150 shocks \(\varepsilon_t\): white noise (\(\phi = 0\)), an AR(1) with \(\phi = 0.9\), and a random walk (\(\phi = 1\)):
The AR(1) series wanders in long swings but keeps coming back to 0; the random walk drifts away and does not come back. The next theorem explains the difference.
The variance and correlation formulas¶
Theorem. Let \(\varepsilon_t\) be white noise with variance \(\sigma^2\) and \(|\phi| < 1\).
- The series \(X_t = \sum_{j=0}^{\infty}\phi^j \varepsilon_{t-j}\) is well defined, satisfies \(X_t = \phi X_{t-1} + \varepsilon_t\), and is weakly stationary with mean 0.
- Its variance is \(\gamma(0) = \dfrac{\sigma^2}{1 - \phi^2}\).
- Its autocorrelations are \(\rho(h) = \phi^{|h|}\).
- Conversely, any weakly stationary solution of the AR(1) equation has variance \(\sigma^2/(1-\phi^2)\); and when \(|\phi| \ge 1\) no weakly stationary solution with \(\varepsilon_t\) uncorrelated with the past exists.
Proof. Part 1 and 2. Consider the partial sums \(Y_N = \sum_{j=0}^{N}\phi^j \varepsilon_{t-j}\). The terms are uncorrelated, so for \(M < N\),
by the geometric series with ratio \(\phi^2 < 1\). So the partial sums form a Cauchy sequence in mean square and converge to a limit \(X_t\) (this uses the completeness of the space of random variables with finite variance; for this book it is enough to know that the limit exists, with the mean and variance that the formulas below give). Its mean is 0, and
Splitting off the \(j = 0\) term, \(X_t = \varepsilon_t + \phi\sum_{j=0}^{\infty}\phi^j\varepsilon_{t-1-j} = \varepsilon_t + \phi X_{t-1}\), so it solves the equation; and \(X_{t-1}\) is built from shocks up to time \(t-1\), which are uncorrelated with \(\varepsilon_t\).
Part 3. For \(h \ge 1\), the equation gives \(X_{t+h} = \phi X_{t+h-1} + \varepsilon_{t+h}\), and \(\varepsilon_{t+h}\) is uncorrelated with \(X_t\) (which uses only shocks up to time \(t\)). Taking covariance with \(X_t\):
Repeating, \(\gamma(h) = \phi^h\gamma(0)\), so \(\rho(h) = \phi^h\); by symmetry \(\rho(-h) = \rho(h)\). These depend only on \(h\), which completes the proof of stationarity.
Part 4. Suppose \(X_t\) is weakly stationary and solves the equation, with \(\varepsilon_t\) uncorrelated with \(X_{t-1}\). Take variances of both sides of \(X_t = \phi X_{t-1} + \varepsilon_t\); the cross term vanishes, so
If \(|\phi| < 1\) this gives \(\gamma(0) = \sigma^2/(1-\phi^2)\). If \(|\phi| \ge 1\) the left side is \(\le 0\) while \(\sigma^2 > 0\), a contradiction, so no such stationary solution exists. \(\blacksquare\)
Two things are worth noticing in the formula \(\sigma^2/(1-\phi^2)\). Memory inflates variance: with \(\phi = 0.9\) the series has \(1/(1 - 0.81) = 5.26\) times the variance of its shocks, because shocks pile up before they fade. And as \(\phi \to 1\) the variance blows up: the stationary picture breaks down.
The random walk: \(\phi = 1\)¶
With \(\phi = 1\) and starting value \(X_0 = 0\), \(X_t = \varepsilon_1 + \dots + \varepsilon_t\), so \(\operatorname{Var}(X_t) = t\sigma^2\), growing without bound: the random walk of Probability, Chapter 15. It is not stationary, but its differences \(X_t - X_{t-1} = \varepsilon_t\) are. Such a series is called integrated or said to have a unit root. Many economic and financial series behave like this, which is why analysts difference them first.
Starting from a fixed value. If an AR(1) with \(|\phi| < 1\) starts at \(X_0 = 0\) instead of in its stationary state, then \(X_t = \sum_{j=0}^{t-1}\phi^j\varepsilon_{t-j}\) and \(\operatorname{Var}(X_t) = \sigma^2(1 - \phi^{2t})/(1 - \phi^2)\), which rises towards \(\sigma^2/(1-\phi^2)\) geometrically fast. That is why simulations throw away a start-up stretch.
Estimating \(\phi\)¶
Because \(\rho(1) = \phi\), a natural estimate is the lag-1 sample autocorrelation \(r_1\) (the Yule-Walker estimate). Regressing \(x_t\) on \(x_{t-1}\) by least squares gives almost the same number. For a stationary AR(1), \(\hat\phi\) is approximately normal with mean \(\phi\) and variance \((1-\phi^2)/n\) for large \(n\) (stated here without proof; see Brockwell and Davis, chapter 8).
Try it in Python¶
import random, statistics
random.seed(1)
phi, n = 0.7, 500
x = [0.0]
for t in range(1, n + 200):
x.append(phi * x[-1] + random.gauss(0, 1))
x = x[200:] # drop the start-up stretch
m = statistics.mean(x)
c0 = sum((v - m) ** 2 for v in x) / n
c1 = sum((x[t] - m) * (x[t - 1] - m) for t in range(1, n)) / n
print(round(c0, 2), round(1 / (1 - phi ** 2), 2)) # sample variance vs theory
print(round(c1 / c0, 3)) # estimate of phi
It prints a sample variance of 2.11 against the theoretical 1.96, and \(\hat\phi = 0.732\) against the true 0.7, within about one standard error (\(\sqrt{(1 - 0.49)/500} = 0.032\)).
Testing for a unit root is harder than it looks: when \(\phi = 1\), \(\hat\phi\) is not approximately normal, and the usual \(t\) statistic for \(\phi = 1\) has a different distribution, tabulated by David Dickey and Wayne Fuller (1979). Using ordinary \(t\) tables there is a classic error.
Beyond AR(1)¶
The AR(\(p\)) model \(X_t = \phi_1 X_{t-1} + \dots + \phi_p X_{t-p} + \varepsilon_t\) allows longer memory and cycles; Yule used an AR(2) in 1927 to model sunspot numbers, finding a quasi-period of roughly a decade, and Gilbert Walker extended the method in 1931. Adding moving-average terms in past shocks gives ARMA models, and differencing first gives ARIMA; George Box and Gwilym Jenkins (1970) turned these into a complete method of identification, estimation and checking.
16.5 Forecasting honestly¶
Level: advanced.
The best forecast for an AR(1)¶
Suppose \(X_t - \mu = \phi(X_{t-1} - \mu) + \varepsilon_t\) with the \(\varepsilon_t\) independent \(N(0, \sigma^2)\) and \(\phi\), \(\mu\), \(\sigma\) known. Having seen the series up to time \(n\), what should we predict for time \(n + h\)?
Apply the equation \(h\) times:
The first term is known at time \(n\); the sum involves only future shocks, which have mean 0 and are independent of everything up to time \(n\). The forecast that minimises the expected squared error is the conditional mean given the past (the fact from Chapter 3 that the mean minimises squared error, applied to the conditional distribution), so
The forecast decays geometrically from the last value back to the long-run mean. And the forecast error is the sum of future shocks, with variance
At \(h = 1\) this is \(\sigma^2\) (only one new shock); as \(h \to \infty\) it rises to \(\sigma^2/(1-\phi^2)\), the variance of the series itself. Far enough ahead, knowing today's value is useless, and the best forecast is the mean with the series' full spread around it.
Worked example. An AR(1) with \(\mu = 20\), \(\phi = 0.6\), \(\sigma = 2\), and last value \(x_n = 26\):
| steps ahead \(h\) | forecast \(20 + 0.6^h \times 6\) | error s.d. \(2\sqrt{(1 - 0.36^h)/0.64}\) | 95% interval |
|---|---|---|---|
| 1 | 23.60 | 2.00 | 19.7 to 27.5 |
| 2 | 22.16 | 2.33 | 17.6 to 26.7 |
| 3 | 21.30 | 2.44 | 16.5 to 26.1 |
| 5 | 20.47 | 2.49 | 15.6 to 25.4 |
The limit of the error standard deviation is \(2/\sqrt{0.64} = 2.5\). Check: at \(h = 2\), \(0.6^2 \times 6 = 2.16\) and \(2\sqrt{1 + 0.36} = 2.332\).
Why real forecast intervals are too narrow¶
The intervals above assume we know \(\phi\), \(\mu\) and \(\sigma\), and that the model is right. In practice:
- Parameter uncertainty. Estimated parameters add error, roughly of order \(1/n\) relative to the shock variance for stationary models; more for long horizons and near unit roots.
- Model uncertainty. The chosen model is one of many that fit the past about equally well, and they disagree about the future. Intervals computed from one model ignore this, and are the main reason published intervals cover less often than they claim.
- Change. The future may not resemble the past at all: a new policy, a pandemic, a broken sensor. No model fitted to the past can see a change that has not happened yet.
Testing forecasts honestly¶
The only honest test of a forecasting method is how well it forecast data it had not seen.
- Hold out the end. Fit to the first part of the series, forecast the rest, and compare. Never choose the model by how well it fits the data it was fitted to; that rewards overfitting (Chapter 12's lesson, in time order).
- Rolling origin. Repeat: fit to data up to time \(t\), forecast \(t + 1, \dots, t + h\), move \(t\) forward, and average the errors. This is cross-validation that respects time; ordinary random-split cross-validation leaks the future into the past.
- Beat a benchmark. Compare with the naive forecast (tomorrow equals today) and the seasonal naive forecast (this January equals last January). A method that cannot beat them is not worth its complexity. Scaled error measures such as the mean absolute scaled error of Hyndman and Koehler (2006) build this comparison in.
The forecasting competitions organised by Spyros Makridakis and colleagues since 1982 put these ideas to large-scale tests on thousands of real series. Their recurring lessons have been that sophisticated methods do not automatically beat simple ones, and that combinations of several methods are hard to beat; in the M4 competition (results published in 2020) several pure machine-learning entries did worse than simple statistical benchmarks, while in M5 (retail sales, 2022) tree-based machine-learning methods did best. What works depends on the data, which is exactly why out-of-sample testing matters.
16.6 A Mars example: the seasonal pressure cycle¶
Level: intermediate, with an advanced part on harmonic regression.
What was measured¶
When Viking Lander 1 (in Chryse Planitia) and Viking Lander 2 (in Utopia Planitia) set down on Mars in 1976, each carried a meteorology instrument that recorded air pressure, temperature and wind. Viking Lander 1 kept sending data until 1982, covering several Mars years. The pressure record turned out to be one of the most striking time series in planetary science (Hess and colleagues, 1980).
Mars' air is mostly carbon dioxide, and at the surface its pressure is less than one percent of Earth's. In the polar winter it is cold enough for carbon dioxide to freeze out of the atmosphere onto the polar cap as frost and ice. So as each pole goes into winter, it pulls a large part of the atmosphere down onto the ground, and the pressure everywhere falls; in spring the ice turns back to gas and the pressure rises. Over a Mars year (687 Earth days, about 669 Mars days or sols) a substantial fraction of the whole atmosphere, often quoted as roughly a quarter, cycles between the air and the caps.
The seasonal curve is not a simple sine wave. It has two minima and two maxima per Mars year. Mars' orbit is noticeably eccentric, and southern winter happens near the far point of the orbit, so it is longer and colder than northern winter. The southern cap therefore takes up more carbon dioxide, and the deepest pressure minimum comes in southern winter, with a shallower one in northern winter. On top of the seasonal cycle there are daily oscillations driven by the Sun's heating (atmospheric thermal tides), weather systems passing through, and changes during dust storms, which Viking observed in 1977.
Later missions repeated the measurement: Mars Pathfinder (1997), Phoenix (2008), the Curiosity rover's REMS instrument (since 2012) and the InSight lander (2018 to 2022), among others. Each site sees the same basic shape, at a different overall level because pressure falls with height (the atmosphere's scale height, the height over which pressure falls by a factor \(e\), is about 11 km).
How a time-series analyst approaches it¶
Everything in this chapter applies. We use no numbers from the missions here, only the methods.
- Choose the time variable. For seasons on Mars the natural clock is not the Earth date but the solar longitude \(L_s\), the angle of Mars around its orbit, from \(0°\) (northern spring) to \(360°\). Using \(L_s\) lines up the seasons of different Mars years exactly.
- Separate time scales. Average each sol's readings to remove the daily tide (a moving average over one sol, as in 16.3), or model the daily cycle separately if it is the object of interest.
- Decompose. Daily-mean pressure = seasonal component + weather and noise. Because different sites have different levels, a multiplicative description, pressure divided by its annual mean, puts them on a common scale (16.1).
- Fit the season by harmonic regression (below), with at least two harmonics to capture the two minima.
- Look at the residuals. They are not white noise: weather systems make them autocorrelated over days, so the ACF (16.2) and models like AR(1) (16.4) are needed before attaching standard errors to anything.
- Compare years and missions. The same season in different Mars years, at the same site, can be compared directly. Comparing missions decades apart is a way to ask whether the amount of carbon dioxide stored in the caps is changing. It requires correcting for elevation, sensor calibration and the weather of particular years, which is exactly the kind of careful, honest time-series work this chapter describes.
Harmonic regression¶
Level: advanced.
A periodic seasonal component with period one year can be written as a Fourier series in the season angle \(s\) (for Mars, \(s = L_s\) in radians):
The \(k = 1\) terms give one rise and fall per year; \(k = 2\) adds a second, which is what makes two minima possible. This is a linear regression (Chapter 12) on the predictors \(\cos s, \sin s, \cos 2s, \sin 2s, \dots\), so least squares applies directly.
When the observations are equally spaced over whole cycles, the least squares solution has a simple form, because the predictors are orthogonal. Take \(N\) observations at angles \(s_m = 2\pi m/N\), \(m = 0, \dots, N-1\). For any whole number \(k\) that is not a multiple of \(N\), let \(z = e^{2\pi \mathrm{i} k/N}\), where \(\mathrm{i} = \sqrt{-1}\); then \(z \ne 1\) and
since \(z^N = e^{2\pi \mathrm{i} k} = 1\). Taking real and imaginary parts, \(\sum_m \cos ks_m = \sum_m \sin ks_m = 0\). Writing products as sums (\(\cos A\cos B = \frac12[\cos(A-B) + \cos(A+B)]\) and similar), it follows that for \(1 \le j, k < N/2\) the columns \(\cos js\), \(\sin js\), \(\cos ks\), \(\sin ks\) are orthogonal to each other and to the constant, with \(\sum_m \cos^2 ks_m = \sum_m \sin^2 ks_m = N/2\). The normal equations then separate, one unknown at a time, and give
These are the discrete Fourier coefficients of the data.
Try it in Python (made-up data, not Mars data). The series below is invented: a seasonal shape with two harmonics, sampled 24 times per cycle for three cycles, plus noise. The fit recovers the shape.
import math, random
random.seed(3)
N, K = 72, 2 # 72 equally spaced points over 3 cycles
true = lambda s: 8 - 0.9 * math.cos(s) + 0.5 * math.sin(s) + 0.6 * math.cos(2 * s)
season = [2 * math.pi * (i % 24) / 24 for i in range(N)]
y = [true(s) + random.gauss(0, 0.2) for s in season]
a0 = sum(y) / N
for k in range(1, K + 1):
a = 2 / N * sum(v * math.cos(k * s) for v, s in zip(y, season))
b = 2 / N * sum(v * math.sin(k * s) for v, s in zip(y, season))
print(k, round(a, 2), round(b, 2))
print(round(a0, 2))
It prints \(\hat a_1 = -0.92\), \(\hat b_1 = 0.50\), \(\hat a_2 = 0.63\), \(\hat b_2 = -0.02\) and \(\hat a_0 = 7.95\), close to the true \(-0.9\), \(0.5\), \(0.6\), \(0\) and \(8\). With real data the samples are not perfectly equally spaced and some sols are missing, so one solves the normal equations in general, but the idea is the same. And the residuals must be checked for autocorrelation before trusting any standard error that the regression prints.
Check yourself¶
1. For the ten counts \(3, 5, 4, 6, 5, 3, 2, 4, 5, 3\), verify \(r_1 = 0.071\) and \(r_2 = -0.286\), and say whether either is significant at 5%.
Deviations from the mean 4: \(-1, 1, 0, 2, 1, -1, -2, 0, 1, -1\); sum of squares 14, so \(\hat\gamma(0) = 1.4\). Lag-1 products sum to 1 (\(\hat\gamma(1) = 0.1\), \(r_1 = 0.071\)); lag-2 products sum to \(-4\) (\(\hat\gamma(2) = -0.4\), \(r_2 = -0.286\)). The white-noise bands are \(\pm 1.96/\sqrt{10} = \pm 0.62\); neither is outside them.
2. An AR(1) has \(\phi = 0.8\) and shock variance 1. Find its variance and its autocorrelation at lag 3.
\(\gamma(0) = 1/(1 - 0.64) = 2.78\). \(\rho(3) = 0.8^3 = 0.512\).
3. Describe the ACF of an AR(1) with \(\phi = -0.5\), and the look of the series.
\(\rho(h) = (-0.5)^h\): \(-0.5, 0.25, -0.125, \dots\), alternating in sign and halving. The series tends to jump from above its mean to below it and back: a zigzag.
4. For a random walk \(X_t = \varepsilon_1 + \dots + \varepsilon_t\) with shock variance \(\sigma^2\), find \(\operatorname{Cov}(X_s, X_t)\) for \(s \le t\) and the correlation between \(X_{10}\) and \(X_{20}\).
\(X_t = X_s + (\varepsilon_{s+1} + \dots + \varepsilon_t)\), and the bracket is uncorrelated with \(X_s\), so \(\operatorname{Cov}(X_s, X_t) = \operatorname{Var}(X_s) = s\sigma^2\). Correlation: \(10\sigma^2/\sqrt{10\sigma^2 \cdot 20\sigma^2} = \sqrt{10/20} = 0.707\). It depends on \(s\) and \(t\), not only on \(t - s\): not stationary.
5. With simple exponential smoothing and \(\alpha = 0.3\), what weight does the forecast put on the most recent observation, the one before it, and the one five steps back?
\(\alpha = 0.3\), \(\alpha(1-\alpha) = 0.21\), and \(\alpha(1-\alpha)^5 = 0.3 \times 0.16807 = 0.050\).
6. Show that the 2×12 moving average turns any seasonal pattern with period 12 and sum zero into exactly zero.
Over the 13 terms from \(t - 6\) to \(t + 6\), months \(t - 6\) and \(t + 6\) are the same calendar month and get weight \(\frac1{24}\) each, a total of \(\frac1{12}\); the eleven months in between each get \(\frac1{12}\). So every calendar month gets total weight \(\frac1{12}\), and the average of \(S\) is \(\frac1{12}(S_1 + \dots + S_{12}) = 0\).
7. Using the AR(1) with \(\mu = 20\), \(\phi = 0.6\), \(\sigma = 2\) and last value 26, give the 3-step-ahead forecast and a 95% forecast interval.
\(\hat X = 20 + 0.6^3 \times 6 = 20 + 1.296 = 21.30\). Error variance \(4(1 + 0.36 + 0.1296) = 5.958\), s.d. 2.441. Interval \(21.30 \pm 1.96 \times 2.441 = (16.5, 26.1)\).
8. A series has a linear trend, \(x_t = a + bt + \varepsilon_t\) with white noise \(\varepsilon_t\). What does differencing give? Is the result white noise?
\(x_t - x_{t-1} = b + \varepsilon_t - \varepsilon_{t-1}\). The trend becomes a constant \(b\), but the noise becomes \(\varepsilon_t - \varepsilon_{t-1}\), which has lag-1 correlation \(-\sigma^2/(2\sigma^2) = -0.5\): not white noise. Differencing a series that did not need it introduces negative correlation (over-differencing).
9. For the moving-average model \(X_t = \varepsilon_t + \theta\varepsilon_{t-1}\), find \(\rho(1)\) and \(\rho(h)\) for \(h \ge 2\). What is the largest possible \(|\rho(1)|\)?
\(\gamma(0) = (1 + \theta^2)\sigma^2\), \(\gamma(1) = \operatorname{Cov}(\varepsilon_t + \theta\varepsilon_{t-1}, \varepsilon_{t+1} + \theta\varepsilon_t) = \theta\sigma^2\), and \(\gamma(h) = 0\) for \(h \ge 2\) (no shared shocks). So \(\rho(1) = \theta/(1 + \theta^2)\). Since \(1 + \theta^2 \ge 2|\theta|\), \(|\rho(1)| \le \frac12\), with equality at \(\theta = \pm1\).
10. Why does the ordinary test for a correlation coefficient mislead when applied to two independent random walks?
The test assumes the pairs are independent across time. Random walks are strongly autocorrelated: each wanders in long excursions, so the effective number of independent observations is tiny, and two unrelated walks often share long upward or downward stretches by chance. Large sample correlations appear far more often than the test's 5%, Yule's nonsense correlation. Correlating the differences, which are independent, gives a valid test.
11. Show that \(\sum_{m=0}^{N-1}\cos(2\pi k m/N) = 0\) when \(k\) is a whole number that is not a multiple of \(N\), and check it for \(N = 4\), \(k = 1\).
With \(z = e^{2\pi \mathrm{i} k/N} \ne 1\), \(\sum_{m=0}^{N-1} z^m = (z^N - 1)/(z - 1) = 0\) because \(z^N = 1\); the cosine sum is the real part. For \(N = 4\), \(k = 1\): \(\cos 0 + \cos\frac\pi2 + \cos\pi + \cos\frac{3\pi}2 = 1 + 0 - 1 + 0 = 0\).
12. (Harder.) An AR(1) started at \(X_0 = 0\) with \(|\phi| < 1\). Show that \(\operatorname{Var}(X_t) = \sigma^2(1 - \phi^{2t})/(1 - \phi^2)\), and find how many steps it takes to reach 99% of the stationary variance when \(\phi = 0.9\).
\(X_t = \sum_{j=0}^{t-1}\phi^j\varepsilon_{t-j}\), a sum of \(t\) uncorrelated terms, so the variance is \(\sigma^2\sum_{j=0}^{t-1}\phi^{2j} = \sigma^2(1 - \phi^{2t})/(1 - \phi^2)\). It is 99% of the limit when \(\phi^{2t} \le 0.01\), that is \(t \ge \log 0.01/(2\log 0.9) = 21.9\): 22 steps.
What you learned¶
- A time series is often described as trend + season + noise (or a product of them); the decomposition separates parts with different causes, like the trend and the growing-season cycle in the Keeling curve.
- The autocorrelation function measures memory; for white noise, sample autocorrelations are roughly \(N(0, 1/n)\), and trending series produce nonsense correlations.
- Centred moving averages estimate trends and remove seasons; exponential smoothing forecasts with geometrically decaying weights; smoothing noise can create cycles that are not there.
- A weakly stationary AR(1) needs \(|\phi| < 1\), and then has variance \(\sigma^2/(1-\phi^2)\) and autocorrelations \(\phi^h\); at \(\phi = 1\) it becomes a random walk whose variance grows like \(t\).
- AR(1) forecasts decay to the mean as \(\phi^h\), and their error variance rises to the series' variance; real intervals must also allow for parameter and model uncertainty, and methods must be judged out of sample against naive benchmarks.
- Mars' surface pressure has a strong seasonal cycle, with two minima per Mars year, because carbon dioxide freezes onto and sublimes from the polar caps; harmonic regression in the solar longitude is the natural way to model it.
Next: Chapter 17, from patterns to causes: designed experiments, potential outcomes and what can be learned from observational data.
Further reading¶
- G. U. Yule, "Why do we sometimes get nonsense-correlations between time-series?", Journal of the Royal Statistical Society 89, 1926, and "On a method of investigating periodicities in disturbed series, with special reference to Wolfer's sunspot numbers", Philosophical Transactions of the Royal Society A 226, 1927.
- E. Slutsky, "The summation of random causes as the source of cyclic processes", Econometrica 5, 1937.
- G. E. P. Box and G. M. Jenkins, Time Series Analysis: Forecasting and Control, Holden-Day, 1970 (later editions with G. C. Reinsel and G. M. Ljung).
- P. J. Brockwell and R. A. Davis, Time Series: Theory and Methods, 2nd edition, Springer, 1991, and their Introduction to Time Series and Forecasting, 3rd edition, Springer, 2016.
- R. J. Hyndman and G. Athanasopoulos, Forecasting: Principles and Practice, 3rd edition, OTexts, 2021. Free online, practical and honest about forecast evaluation.
- R. H. Shumway and D. S. Stoffer, Time Series Analysis and Its Applications, 4th edition, Springer, 2017.
- D. A. Dickey and W. A. Fuller, "Distribution of the estimators for autoregressive time series with a unit root", Journal of the American Statistical Association 74, 1979.
- S. L. Hess, R. M. Henry, C. B. Leovy, J. A. Ryan and J. E. Tillman, "Meteorological results from the surface of Mars: Viking 1 and 2", Journal of Geophysical Research 82, 1977, and S. L. Hess and colleagues, "The annual cycle of pressure on Mars measured by Viking Landers 1 and 2", Geophysical Research Letters 7, 1980.
- C. B. Leovy, "Weather and climate on Mars", Nature 412, 2001.
- On this site: Series (geometric series), Random walks and martingales and Markov chains in the Probability book, Correlation and Correlation and causation (gentler), and the Statistics book's contents.