---
title: "Workbook 8 — Choosing a Transformation"
subtitle: "STA 542 · Introduction to Time Series Analysis"
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 2
    number-sections: false
    code-copy: true
    fig-width: 9
    fig-height: 3.8
    fig-dpi: 150
execute:
  warning: false
  message: false
editor:
  markdown:
    wrap: 72
---

::: hidden
$$\def\E{\mathbb{E}}
\def\P{\mathbb{P}}
\def\R{\mathbb{R}}
\def\Z{\mathbb{Z}}
\def\Var{\operatorname{Var}}
\def\Cov{\operatorname{Cov}}
\def\Corr{\operatorname{Corr}}$$
:::

```{=html}
<style>
.exercise { border-left: 3px solid #00539B; padding: 0.7em 1.1em; background: #f6f8fa; margin: 1.4em 0; border-radius: 0 4px 4px 0; }
.exercise p:first-child { margin-top: 0; }
figure figcaption { color: #555; font-size: 0.9em; }
</style>
```

```{r}
#| label: setup
#| include: false
knitr::opts_chunk$set(comment = "#>")
```

## What Changes in the Observed Record?

In earlier workbooks, we saw that weak stationarity allows us a common mean and ACVF let us combine observations taken at different times. An observed record with a changing level or seasonal pattern raises a prior
question: how do we transform a timeseries such that it could reasonably be modeled as weakly stationary?

**Example 8.1 (Monthly airline passengers).** The series
`datasets::AirPassengers` records monthly international airline passenger
totals, in thousands, from January 1949 through December 1960. We seek an
adjustment after which a common mean and ACVF might describe the record.
The time plot and three annual profiles show the observations before any
adjustment.

```{r}
#| label: fig-wb8-air-raw
#| code-fold: true
#| code-summary: "R code for the plots"
#| fig-cap: "Monthly international airline passenger totals, 1949–1960. The three annual profiles compare the same calendar months at different levels. Real data: datasets::AirPassengers."
air <- datasets::AirPassengers
air_years <- 1949:1960
air_matrix <- matrix(as.numeric(air), nrow = 12)
profile_years <- c(1949, 1954, 1960)
profile_cols <- c("#00539B", "#AD5B00", "#40764B")
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), bg = "white")
plot(air, col = "#00539B", xlab = "Year",
     ylab = "Passengers (thousands)", main = "Observed record")
matplot(1:12, air_matrix[, match(profile_years, air_years)],
        type = "l", lty = 1, lwd = 2, col = profile_cols, xaxt = "n",
        xlab = "Calendar month", ylab = "Passengers (thousands)",
        main = "Three annual profiles")
axis(1, at = 1:12, labels = month.abb, cex.axis = 0.75)
legend("topleft", legend = profile_years, col = profile_cols,
       lty = 1, lwd = 2, bty = "n", cex = 0.85)
```

The later years have a higher level and a larger seasonal range. The annual
profiles also share a summer peak. A change of scale, removal of a trend,
and adjustment for calendar month address different parts of this pattern.
$\diamond$

For a weakly stationary process $(Y_t)$, the three requirements are
$$
\E[Y_t]=\mu_Y,\qquad
\Var(Y_t)=\gamma_Y(0)<\infty,\qquad
\Cov(Y_t,Y_{t+h})=\gamma_Y(h).
$$
Each right-hand side is independent of $t$. The covariance may be large at
some lags. Weak stationarity requires the same covariance at a given lag
across time.

Our first objective is a series for which these stationarity assumptions
are plausible. Further, we will see that sometimes multiple transformations
result in weakly stationary series. Depending on their autocorrelation structure
and implied noise innovations, different transformations may be preferred.

## Would a Different Scale Help?

The larger seasonal range at higher passenger levels suggests comparing
relative changes. A log transformation expresses a fixed ratio as a fixed
difference. A log transformation is an example of a point transformation.

**Definition 8.2 (Point transformation).** Let $(X_t)_{t\in\Z}$ be a
process taking values in $D\subseteq\R$. A *point transformation* applies
the same function $h:D\to\R$ at every time, giving
$$
L_t=h(X_t),\qquad t\in\Z.
$$
Each transformed value depends only on the observation at that time.

For example, if the process takes positive values, choosing $h(x)=\log x$
gives the log-transformed process $L_t=\log X_t$.

Suppose
$$
X_t=\exp\{m(t)+Y_t\},
$$
where $m$ is a deterministic function of time and $(Y_t)$ is weakly
stationary with mean zero. Then
$$
L_t=\log X_t=m(t)+Y_t.
$$
The log converts this multiplicative model into an additive mean plus
stationary fluctuations.

If $m$ is constant, $(L_t)$ is already weakly stationary. A changing
$m(t)$ still has to be addressed. Subtracting a known $m$ yields
$L_t-m(t)=Y_t$. If $m$ is unknown, we estimate it first and subtract
the fitted values. Differencing is a second option: it replaces each
observation by its change from the previous time.

The *backshift operator* $B$ is defined by
$$
BX_t:=X_{t-1}.
$$
Powers of $B$ are defined by $B^kX_t:=X_{t-k}$ for $k=0,1,2,\ldots$,
with $B^0$ the identity. The *first-difference operator* is
$\nabla:=1-B$, where $1$ denotes that identity, so
$$
\nabla X_t:=(1-B)X_t=X_t-X_{t-1}.
$$

Applied to the logged series,
$$
\begin{aligned}
\nabla L_t
&=L_t-L_{t-1}\\
&=\bigl(m(t)-m(t-1)\bigr)+(Y_t-Y_{t-1})\\
&=\nabla m(t)+\nabla Y_t.
\end{aligned}
$$
If $m(t)=a+bt$, then $\nabla m(t)=b$, and
$$
\nabla L_t=b+\nabla Y_t.
$$
The mean of $\nabla L_t$ is $b$. The section "How Does Filtering Change
Stationary Dependence?" calculates $\Cov(\nabla Y_t,\nabla Y_{t+h})$,
which does not depend on $t$.
Subtracting the specified linear mean recovers $Y_t$.

The long-run variance workbook formed DJIA log returns by this
construction. If $c_t$ denotes a closing level, the log return is
$\nabla\log c_t$. Modeling $\log c_t=m(t)+Y_t$ with linear $m$ makes
that return equal to $b+\nabla Y_t$.

The calculations above use the displayed multiplicative model. A
nonlinear point transformation of an arbitrary weakly stationary
process need not preserve weak stationarity.

```{r}
#| label: fig-wb8-air-log
#| code-fold: true
#| code-summary: "R code for the plots"
#| fig-cap: "Log passenger totals and the same three annual profiles. Seasonal amplitudes are more comparable on this scale, while the levels of the annual profiles still differ. Real data: datasets::AirPassengers."
log_air <- log(air)
log_air_matrix <- log(air_matrix)
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), bg = "white")
plot(log_air, col = "#00539B", xlab = "Year",
     ylab = "Log passengers", main = "Log record")
matplot(1:12, log_air_matrix[, match(profile_years, air_years)],
        type = "l", lty = 1, lwd = 2, col = profile_cols, xaxt = "n",
        xlab = "Calendar month", ylab = "Log passengers",
        main = "Annual profiles on the log scale")
axis(1, at = 1:12, labels = month.abb, cex.axis = 0.75)
legend("topleft", legend = profile_years, col = profile_cols,
       lty = 1, lwd = 2, bty = "n", cex = 0.85)
```

The profiles support investigating an additive seasonal mean on the log
scale. A plot of within-year standard deviation against annual mean offers
another description of the changing range. That standard deviation includes
the variation in seasonal means, so it does not estimate an innovation
variance.

The DJIA log returns from the long-run variance workbook illustrate
the same distinction. Volatility clustering in those returns is a
separate feature and can coexist with weak stationarity.

## Subtract a Trend or Take Differences?

The log passenger series still rises over time. Subtracting a mean trend
and taking successive differences are two possible adjustments. We first
compare them in simulations where the mean trend is known.

**Example 8.3 (Same Mean Trend, Different Fluctuations).** Records A and B
contain 300 observations from the respective processes
$$
X_t^{(D)}=bt+Z_t,\qquad
X_t^{(S)}=bt+\sum_{j=1}^{t}U_j,\qquad t\geq1,
$$
with $b=0.02$ and independent i.i.d. standard normal sequences $(Z_t)$
and $(U_t)$. Record A has a linear trend plus independent noise. Record B
is a random walk with drift: $X_t^{(S)}=X_{t-1}^{(S)}+b+U_t$, with
$X_0^{(S)}=0$.

Both processes have mean $bt$. In the first, a shock affects only its own
observation. In the second, each shock remains in all subsequent levels.
We first compare the adjustments assuming that the mean trend is known.
For differences within the observed records, $t=2,\ldots,300$.

| Adjustment | Record A: trend plus independent noise | Record B: random walk with drift |
|---|---|---|
| Subtract the known mean $bt$ | $Z_t$: stationary, variance $1$ | $\sum_{j=1}^{t}U_j$: nonstationary, variance $t$ |
| Take first differences | $b+Z_t-Z_{t-1}$: stationary, with dependent successive changes | $b+U_t$: stationary, with independent successive changes |

For Record B, the independent shocks each have variance one, so their
sum has variance $t$. Even subtracting the exact mean leaves a variance
that grows with time. Taking differences instead leaves the independent
increments $b+U_t$, with constant mean $b$ and variance one.

For Record A, both adjustments give stationary outputs. Mean subtraction
retains deviations from the expected level; differencing retains
successive changes. Those changes share shocks: $Z_t$ enters one change
with a plus sign and the next with a minus sign, introducing negative
lag-one correlation. Stationarity alone does not select between these
quantities; their interpretation and dependence structure matter.

The orange line below is the known mean $bt$ in both rows. The middle
column subtracts that line; the last column takes first differences.

```{r}
#| label: fig-wb8-trend-candidates
#| fig-height: 6.1
#| code-fold: true
#| code-summary: "Simulation and adjustment code"
#| fig-cap: "The same known mean trend is subtracted from both records. Record A leaves independent fluctuations; Record B leaves accumulated shocks. First differences are stationary under both models, with dependent changes for A and independent increments for B. Simulated data, seed 542."
set.seed(542)
n_trend <- 300  # T in the mathematical notation.
time_trend <- seq_len(n_trend)
x_deterministic <- 0.02 * time_trend + rnorm(n_trend)
x_accumulated <- cumsum(0.02 + rnorm(n_trend))
known_trend <- 0.02 * time_trend
trend_records <- list("Record A" = x_deterministic,
                      "Record B" = x_accumulated)
par(mfrow = c(2, 3), mar = c(3.5, 4, 3, 1),
    mgp = c(2.2, 0.7, 0), bg = "gray92")
for (name in names(trend_records)) {
  x <- trend_records[[name]]
  plot(time_trend, x, type = "l", col = "#00539B",
       ylim = range(x, known_trend),
       xlab = "Observation", ylab = "Level", main = name)
  lines(time_trend, known_trend, col = "#AD5B00", lwd = 2)
  plot(time_trend, x - known_trend, type = "l", col = "#00539B",
       xlab = "Observation", ylab = "Deviation",
       main = "Subtract known mean bt", cex.main = 0.9)
  abline(h = 0, col = "gray60")
  plot(time_trend[-1], diff(x), type = "l", col = "#00539B",
       xlab = "Observation", ylab = "Change", main = "First difference")
}
```

The growing variance calculated above establishes nonstationarity after
mean subtraction for Record B. The single simulated path illustrates
the accumulated fluctuations. $\diamond$

In observations, the mean trend is usually unknown. We next consider
estimating it and how fitted residuals differ from the underlying
fluctuations.

## Fitting a Mean with Time Regressors

In the stationary causal AR(1) family, every fixed parameter choice with
$|\phi|<1$ gives a weakly stationary process. Estimation learns its mean,
variance, and dependence within that family. Here the unknown coefficients
describe a changing mean. We estimate them to construct residuals that
approximate the unobserved stationary component. In either setting, a
fitted model's properties still need to be assessed against the data.

To allow correlated fluctuations around a linear mean, extend Record A
to $X_t=a+bt+Y_t$, where $(Y_t)$ is weakly stationary with mean zero.
Subtracting a fixed candidate line $u+vt$ gives
$$
X_t-(u+vt)=(a-u)+(b-v)t+Y_t.
$$
This process is weakly stationary exactly when $v=b$: a wrong slope
leaves a changing mean. A wrong intercept only changes the constant
mean. Estimating the slope therefore matters to removing the source of
nonstationarity.

The straight trend $a+bt$ uses two known regressors: the constant $1$
and time $t$. A curved trend $a+bt+ct^2$ adds the known regressor $t^2$.
Both are linear regression models because the unknown coefficients enter
linearly. The curve need not be a straight line in time.

More generally, choose known deterministic functions
$r_0(t),\ldots,r_p(t)$, with $r_0(t)=1$, and write a candidate mean as
$$
m_\beta(t)=\sum_{j=0}^{p}\beta_jr_j(t).
$$
We propose $X_t=m(t)+Y_t$, where the actual mean $m$ belongs to this
class and $(Y_t)$ is weakly stationary with mean zero. Subtracting the
actual $m(t)$ recovers $Y_t$. This extends the known-mean subtraction
used for Record A in Example 8.3.

For observations $x_1,\ldots,x_T$, ordinary least squares (OLS) chooses
the coefficients to minimize
$$
Q_T(\beta):=\sum_{t=1}^{T}\{x_t-m_\beta(t)\}^2.
$$
The regressors form the columns of an ordinary regression design matrix.
If these columns are linearly independent on the observed times, the
minimizing coefficient vector is unique.

Why use this criterion when the fluctuations are correlated? For any
fixed candidate $\beta$, expanding the squares gives
$$
\E\left[\sum_{t=1}^{T}\{X_t-m_\beta(t)\}^2\right]
=\sum_{t=1}^{T}\{m(t)-m_\beta(t)\}^2+T\gamma_Y(0).
$$
The cross terms vanish because $\E[Y_t]=0$. Thus the actual mean
minimizes the expected criterion. This calculation uses no independence
between observations. Accuracy of the observed OLS fit still depends on
the dependence and regressor assumptions; this identity alone proves
neither consistency nor the validity of the usual independent-error
standard errors.

Let $\widehat m(t)=m_{\widehat\beta}(t)$ be the fitted mean and
$y_t=x_t-m(t)$ the unobserved realized fluctuation. Then
$$
e_t:=x_t-\widehat m(t)=y_t+\{m(t)-\widehat m(t)\}.
$$
For a straight trend, the extra term is
$(a-\hat a)+(b-\hat b)t$; for a quadratic it also contains
$(c-\hat c)t^2$.

The fitted coefficients are random functions of the same observations,
so their errors are dependent on the fluctuations being estimated. The
fixed-line calculation above cannot be applied by treating $\hat a$ and
$\hat b$ as fixed. Under the correctly specified OLS mean model, the
residuals have expectation zero, but their covariance can still depend
on position within the finite record. They approximate the stationary
component; they need not themselves have its exact population properties.

We now fit a line to each record from Example 8.3. The left panels compare
the fitted line with the known mean. The right panels compare fitted
residuals with the deviations obtained by subtracting that known mean.

```{r}
#| label: fig-wb8-fitted-trend-comparison
#| code-fold: true
#| code-summary: "Fit lines to the same two records"
#| fig-height: 5.8
#| fig-cap: "Fitted lines and residuals for the same records as Example 8.3. For A, fitted residuals approximate the independent fluctuations. For B, the fit absorbs part of the realized accumulation, although the population mean is also linear. A fitted line alone does not justify a stationary remainder. Simulated data, seed 542."
par(mfrow = c(2, 2), mar = c(3.5, 4, 3, 1),
    mgp = c(2.2, 0.7, 0), bg = "gray92")
for (name in names(trend_records)) {
  x <- trend_records[[name]]
  trend_fit <- lm(x ~ time_trend)
  known_remainder <- x - known_trend
  plot(time_trend, x, type = "l", col = "gray65",
       ylim = range(x, known_trend, fitted(trend_fit)),
       xlab = "Observation", ylab = "Level", main = name)
  lines(time_trend, known_trend, col = "#AD5B00", lwd = 2)
  lines(time_trend, fitted(trend_fit), col = "#00539B", lwd = 1, lty = 2)
  legend("topleft", c("Known mean", "Fitted line"),
         col = c("#AD5B00", "#00539B"), lwd = c(2, 1), lty = c(1, 2),
         bty = "n", cex = 0.7)
  plot(time_trend, known_remainder, type = "l", col = "#AD5B00", lwd = 2,
       ylim = range(known_remainder, residuals(trend_fit)),
       xlab = "Observation", ylab = "Deviation / residual",
       main = "Subtract known or fitted mean")
  abline(h = 0, col = "gray60")
  lines(time_trend, residuals(trend_fit), col = "#00539B", lwd = 0.8)
  legend("topleft", c("Known-mean deviation", "Fitted residual"),
         col = c("#AD5B00", "#00539B"), lwd = c(2, 0.8),
         bty = "n", cex = 0.7)
}
```

For Record A, fitted residuals estimate the stationary fluctuations
$Z_t$. Record B also has a linear population mean, but its deviations
from that mean accumulate shocks. Its fitted line can absorb some of
this realized wandering; fitting the line supplies no justification for
a stationary remainder.

To see why the regressors matter, add $0.0001t^2$ to Record A.
Its mean is now $0.02t+0.0001t^2$, with the same i.i.d. fluctuations.
The two fits below estimate all their coefficients from this record.

```{r}
#| label: fig-wb8-quadratic-trend
#| code-fold: true
#| code-summary: "Fit straight and quadratic trends"
#| fig-height: 3.8
#| fig-cap: "A quadratic trend with the same Gaussian fluctuations as Record A. The orange curve in the first panel is the actual mean. A straight-line fit leaves a curved residual pattern; the quadratic fit estimates deviations around the specified mean. The last panel overlays the unobserved fluctuations used in this simulation. Simulated data, seed 542."
x_quadratic <- x_deterministic + 0.0001 * time_trend^2
straight_fit <- lm(x_quadratic ~ time_trend)
quadratic_fit <- lm(x_quadratic ~ time_trend + I(time_trend^2))
true_quadratic_mean <- 0.02 * time_trend + 0.0001 * time_trend^2
true_fluctuation <- x_deterministic - 0.02 * time_trend
residual_limits <- range(residuals(straight_fit),
                         residuals(quadratic_fit), true_fluctuation)
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1),
    mgp = c(2.2, 0.7, 0), bg = "gray92")
plot(time_trend, x_quadratic, type = "l", col = "#00539B",
     xlab = "Observation", ylab = "Level", main = "Quadratic mean")
lines(time_trend, true_quadratic_mean, col = "#AD5B00", lwd = 2)
plot(time_trend, residuals(straight_fit), type = "l", col = "#00539B",
     xlab = "Observation", ylab = "Residual", ylim = residual_limits,
     main = "Subtract fitted line")
abline(h = 0, col = "gray60")
plot(time_trend, true_fluctuation, type = "l", col = "#AD5B00", lwd = 2,
     xlab = "Observation", ylab = "Fluctuation / residual",
     ylim = residual_limits, main = "Subtract fitted quadratic")
lines(time_trend, residuals(quadratic_fit), col = "#00539B", lwd = 0.8)
legend("topleft", c("True fluctuation", "Fitted residual"),
       col = c("#AD5B00", "#00539B"), lwd = c(2, 0.8),
       bty = "n", cex = 0.65)
```

In the R formula, `I(time_trend^2)` supplies the squared time values as
another predictor. The fit estimates $a$, $b$, and $c$ by the same OLS
criterion as the straight-line fit.

For $X_t=a+bt+ct^2+Y_t$ with stationary, mean-zero $(Y_t)$,
$$
\nabla X_t=b+c(2t-1)+\nabla Y_t,\qquad
\nabla^2X_t=2c+\nabla^2Y_t.
$$
One difference leaves a changing mean when $c\ne0$. Two differences leave
a constant mean and a stationary filtered fluctuation, as the finite-filter
result below verifies. Subtracting the actual quadratic mean instead
retains $Y_t$ itself.

Likewise, cubic and higher-degree trends use the additional regressors
$t^3,t^4,\ldots$. Known seasonal functions can enter the same OLS fit.
The same model $X_t=m(t)+Y_t$ with a stationary remainder and the same
fitted-residual identity apply.

## How Does Filtering Change Stationary Dependence?

For Record A in Example 8.3, differencing leaves $b+\nabla Z_t$.
We now calculate how such weighted combinations change a stationary
process's mean and ACVF.

**Definition 8.4 (Finite linear filter).** For fixed coefficients
$a_0,\ldots,a_q$, define
$$
a(B):=\sum_{j=0}^{q}a_jB^j,\qquad
Y_t:=a(B)X_t=\sum_{j=0}^{q}a_jX_{t-j}.
$$
The first-difference filter has coefficients $1,-1$. The same coefficients
are applied at every time.

**Proposition 8.5 (Moments after finite filtering).** If $(X_t)$ is weakly
stationary with mean $\mu$ and ACVF $\gamma_X$, then the process in
Definition 8.4 is weakly stationary, with
$$
\E[Y_t]=\mu\sum_{j=0}^{q}a_j,\qquad
\gamma_Y(h)=
\sum_{j=0}^{q}\sum_{k=0}^{q}a_ja_k\gamma_X(h+j-k).
$$

*Proof.* Finite linear combinations have finite second moments. Linearity
gives the mean formula. For the covariance,
$$
\begin{aligned}
\Cov(Y_t,Y_{t+h})
&=\sum_{j=0}^{q}\sum_{k=0}^{q}
a_ja_k\Cov(X_{t-j},X_{t+h-k})\\
&=\sum_{j=0}^{q}\sum_{k=0}^{q}
a_ja_k\gamma_X(h+j-k).
\end{aligned}
$$
Neither expression depends on $t$. $\square$

This result starts with a stationary input. For a nonstationary process,
we first use its model equation to determine what a filter removes.

For the general trend model
$$
X_t=a+bt+Y_t,
$$
where $(Y_t)$ is weakly stationary with mean zero and ACVF $\gamma_Y$,
subtracting the mean gives $X_t-(a+bt)=Y_t$. Differencing gives
$$
\nabla X_t=b+Y_t-Y_{t-1}.
$$
Proposition 8.5 applied to $(Y_t)$ shows that the differenced process is
also weakly stationary, with mean $b$ and covariance
$$
\Cov(\nabla X_t,\nabla X_{t+h})
=2\gamma_Y(h)-\gamma_Y(h-1)-\gamma_Y(h+1).
$$
This expression depends only on the lag $h$. Thus both operations can
give stationary outputs even when the original fluctuations are
correlated. Mean subtraction retains $Y_t$; differencing changes its
covariance. Record A is the special case $a=0$ and $Y_t=Z_t$.

**Example 8.6 (Differencing white noise).** If
$(Y_t)\sim\mathrm{WN}(0,\sigma^2)$ with $\sigma^2>0$, then
$\nabla Y_t=Y_t-Y_{t-1}$ has
$$
\gamma_{\nabla Y}(0)=2\sigma^2,\qquad
\gamma_{\nabla Y}(1)=-\sigma^2,\qquad
\gamma_{\nabla Y}(h)=0\quad(|h|>1).
$$
Its lag-one correlation is $-1/2$ and its long-run variance is
$2\sigma^2+2(-\sigma^2)=0$.
Differencing has changed the dependence and the scale relevant to mean
inference. A large negative sample lag-one correlation can motivate
checking for an unnecessary difference; it does not identify its cause.
$\diamond$

**Remark 8.7 (Fixed coefficients and estimated parameters).** Fixed
coefficients in Definition 8.4 may be unknown and require estimation.
The proposition describes the population filter with those coefficients;
substituting estimates computed from the same record introduces fitting
error. Subtracting a fitted trend is a different operation from the
fixed-coefficient lag sum in Definition 8.4: its adjustment depends
explicitly on time and on a fit to the whole record. Both operations can
require parameter estimation.

A moving average can estimate a slowly varying trend, after which the
candidate remainder is $x_t-\hat m_t$. Inspecting the smoothed curve itself
addresses a different object. For example, the two-point average of a
random walk remains nonstationary: Exercise 8.1 derives its growing
variance. A smoother-looking record is therefore insufficient evidence
for a stationary model.

## What Does the Seasonal Pattern Require?

The annual profiles in Example 8.1 show why removing a changing level may
leave a recurring calendar pattern. We can estimate a fixed seasonal mean
or compare observations one full season apart.

For monthly observations, a smooth annual mean can be modeled by
$$
S_t=c\cos(2\pi t/12)+d\sin(2\pi t/12),\qquad
X_t=a+S_t+Y_t,
$$
where $(Y_t)$ is weakly stationary with mean zero. The period of twelve
months is specified; the coefficients $a,c,d$ are unknown. This is
*harmonic regression*: the sine and cosine values are known regressors,
so the same OLS criterion estimates their coefficients.

Including both terms lets the data determine the amplitude and the time
of the peak. Indeed,
$$
c\cos(2\pi t/12)+d\sin(2\pi t/12)
=A\cos(2\pi t/12-\delta),
$$
with $c=A\cos\delta$, $d=A\sin\delta$, and $A=\sqrt{c^2+d^2}$.
We can therefore estimate the two linear coefficients rather than
optimize over amplitude and phase directly. The period must be fixed
for this OLS construction; estimating it is an additional problem.

For a monthly record stored in `x`, the fit is:

```{r}
#| label: wb8-harmonic-fit-template
#| eval: false
seasonal_time <- seq_along(x)
seasonal_fit <- lm(as.numeric(x) ~ cos(2*pi*seasonal_time/12) +
                    sin(2*pi*seasonal_time/12))
seasonal_residual <- residuals(seasonal_fit)
```

The residuals estimate deviations from the mean $a+S_t$. If the proposed
mean also contains a trend, we include its regressors in the same fit.
The mean-estimation error discussed above still applies.

One sine--cosine pair describes a single smooth wave per year. Adding
terms such as $\cos(4\pi t/12)$ and $\sin(4\pi t/12)$ permits a more
complicated shape. Separate calendar-month indicators allow an arbitrary
fixed monthly pattern. These are choices of regressors for the same
mean-fitting method; each assumes that the seasonal mean repeats across
years.

**Definition 8.8 (Seasonal difference).** If one seasonal cycle contains
$s$ observations, the *seasonal-difference filter* is
$$
\nabla_s:=1-B^s,\qquad
\nabla_sX_t=X_t-X_{t-s}.
$$
For monthly observations with an annual cycle, $s=12$.

**Example 8.9 (Two seasonal structures).** Suppose
$$
X_t=c_t+Y_t,\qquad c_{t+s}=c_t,
$$
with deterministic $(c_t)$ and weakly stationary $(Y_t)$. Subtracting
$c_t$ recovers $Y_t$, whereas seasonal differencing produces
$Y_t-Y_{t-s}$. The latter is stationary but has a different ACVF.

For example, if $(Y_t)\sim\mathrm{WN}(0,\sigma^2)$ with $\sigma^2>0$,
subtracting $c_t$ leaves variance $\sigma^2$ and zero correlations at
nonzero lags. Seasonal differencing instead gives
$$
\Var(Y_t-Y_{t-s})=2\sigma^2,\qquad
\Corr(Y_t-Y_{t-s},Y_{t+s}-Y_t)=-\tfrac12.
$$
Mean subtraction retains deviations from the expected seasonal level;
seasonal differencing retains changes from one season to the next.
Either quantity may be useful, although their covariance structures differ.

If instead
$$
X_t=X_{t-s}+Y_t
$$
with weakly stationary $(Y_t)$, then $\nabla_sX_t=Y_t$. This equation
provides a different justification for the same filter. With independent
nondegenerate shocks and fixed starting values, the variance accumulates
within each seasonal subsequence; subtracting a fixed seasonal curve
cannot remove it. Exercise 8.3 compares the two structures, including
the effect of estimating the seasonal mean. $\diamond$

A changing seasonal pattern alone does not establish seasonal
accumulation or justify differencing. For either adjustment, we assess
the remaining mean, spread, and lag relationships.

**Example 8.10 (Competing airline adjustments).** On the log scale, one
candidate model is
$$
L_t=a+bt+c_{j(t)}+Y_t,\qquad L_t:=\log X_t,
$$
where $j(t)\in\{1,\ldots,12\}$ is calendar month, $c_1=0$ fixes the
intercept convention, and $(Y_t)$ is weakly stationary with mean zero.
Removing the specified mean would recover $Y_t$. Since $a$, $b$, and
the monthly effects are unknown, we fit them by OLS with regressors $1$,
$t$, and eleven calendar-month indicators. Omitting the January indicator
implements $c_1=0$ and avoids a redundant column.

```{r}
#| label: wb8-air-calendar-fit
air_index <- seq_along(log_air)
air_month <- factor(cycle(log_air))
air_mean_fit <- lm(as.numeric(log_air) ~ air_index + air_month)
air_mean_residual <- ts(residuals(air_mean_fit),
                        start = start(air), frequency = 12)
```

The fit removes a linear trend and a separate mean effect for each
calendar month. Its residuals provide one candidate record. We use
regression to estimate the mean here; its usual independent-error
standard errors would require additional assumptions.

Two further candidates are
$$
W_t:=\nabla_{12}L_t,\qquad
V_t:=\nabla\nabla_{12}L_t.
$$
A model with stationary year-over-year log changes motivates $W_t$.
A model in which those year-over-year changes have stationary increments
motivates $V_t$. A linear deterministic log trend can also be removed by a
seasonal difference, leaving a constant; observed success of a difference
does not establish seasonal accumulation.

```{r}
#| label: wb8-air-difference-candidates
air_seasonal_difference <- diff(log_air, lag = 12)
air_both_differences <- diff(air_seasonal_difference)
```

The product filter expands as
$$
\nabla\nabla_{12}L_t
=L_t-L_{t-1}-L_{t-12}+L_{t-13}.
$$
It compares this month's log change with the log change in the same month
one year earlier. Ordinary and seasonal differences commute because their
finite polynomials in $B$ commute.

The three outputs have different interpretations:

| Candidate | Quantity being modeled | Parameters estimated for the adjustment |
|---|---|---|
| Fitted mean residual | Deviation of log passengers from a linear trend and fixed monthly effects | Intercept, slope, and monthly effects |
| Seasonal difference $W_t$ | Year-over-year log growth, $\log(X_t/X_{t-12})$ | None |
| Both differences $V_t$ | Month-to-month change in year-over-year log growth, $W_t-W_{t-1}$ | None |

A stationary model for $V_t$ describes changes in annual log growth.
It does not directly describe the deviations represented by the fitted
mean residuals. We assess stationarity for each proposed quantity before
choosing which to model. $\diamond$

## Can We Recover the Original Observations?

Taking logarithms can be undone: if $y_t=\log x_t$ with $x_t>0$, then
$x_t=\exp(y_t)$ recovers each original observation.

Subtracting a known trend can also be undone. If $y_t=x_t-m(t)$, then
adding that same trend back gives $x_t=y_t+m(t)$.

First differencing cannot be undone uniquely. The records $(2,5,4)$
and $(12,15,14)$ both give differences $(3,-1)$. More generally, adding
any constant to every observation leaves the differences unchanged.
Differencing removes the overall level.

**Definition 8.11 (Inverse transformations).** An *inverse transformation*
is a rule that recovers every original series from its transformed
series. A transformation is *invertible* if an inverse exists. If
distinct original series produce the same transformed series, the
transformation is not invertible.

Knowing the first observation resolves the ambiguity. Given $x_1$ and
all successive differences $d_t=x_t-x_{t-1}$, $t=2,\ldots,T$, repeated
addition gives
$$
x_t=x_1+\sum_{j=2}^{t}d_j,\qquad t=2,\ldots,T.
$$
The differences determine the changes; the initial observation supplies
the missing level.

## Which Candidate Supports a Stationary Approximation?

For each candidate we compare chronological portions of the record:

- Is a common mean plausible, allowing for persistent runs?
- Is the spread comparable across time?
- Are the relationships at the same lag comparable across time?

Time plots and local means and standard deviations address the first two
questions. Local ACFs address the third after normalizing by variance,
so they must be read alongside the spread comparison. Sample summaries
can differ under stationarity; these checks are descriptive, with no
calibrated rejection threshold.

Differencing drops initial observations, whereas fitting a mean uses the
whole record. We compare the three outputs on their common dates,
February 1950 through December 1960, retaining 131 observations in each.

```{r}
#| label: wb8-air-align
air_candidates <- ts.intersect(
  "Fitted mean residual" = air_mean_residual,
  "Seasonal difference" = air_seasonal_difference,
  "Both differences" = air_both_differences
)
stopifnot(nrow(air_candidates) == 131)
```

The time plots show when deviations occur; the ACFs summarize lag
relationships across each candidate's whole record. The horizontal axes
of the ACFs below count observations, so lag 12 means twelve months.
Reference bands are omitted because these plots describe dependence
rather than implement a calibrated test.

```{r}
#| label: fig-wb8-air-candidate-comparison
#| code-fold: true
#| code-summary: "R code for the candidate comparison"
#| fig-height: 8
#| fig-cap: "Three candidate adjustments on the same 131 months. The fitted-mean residuals and seasonal differences retain persistent stretches; both differences leave visible short-lag and annual-lag dependence. All ACF lags are in months. Real data: datasets::AirPassengers."
air_candidate_labels <- c("Log residual", "12-month log change",
                          "Change in 12-month log change")
par(mfrow = c(3, 2), mar = c(3.5, 4, 3, 1),
    mgp = c(2.2, 0.7, 0), bg = "white")
for (j in seq_len(ncol(air_candidates))) {
  plot(as.numeric(time(air_candidates)), as.numeric(air_candidates[, j]), type = "l",
       col = "#00539B", xlab = "Year", ylab = air_candidate_labels[j],
       main = colnames(air_candidates)[j])
  abline(h = mean(air_candidates[, j]), col = "gray60", lty = 2)
  acf(as.numeric(air_candidates[, j]), lag.max = 24, ci = 0,
      xlab = "Lag (months)", main = "Whole-record sample ACF")
}
```

The fitted-mean residuals sum to zero over the original regression window.
They can nevertheless have sustained stretches above or below zero.
Likewise, a sample ACF pooled over all dates can conceal different lag
relationships in different portions of the record.

We now use the first 65 common observations, February 1950–June 1955, and
the remaining 66, July 1955–December 1960. In each portion we compute its
own sample mean, standard deviation, and sample autocorrelations at lags 1
and 12, using the divisor equal to that portion's length for the sample
ACVF. Pairs crossing the split do not enter either portion's ACF.

```{r}
#| label: wb8-air-portion-table
#| code-fold: true
#| code-summary: "R code for the descriptive summaries"
air_cut <- floor(nrow(air_candidates) / 2)
air_portions <- list(
  "Feb 1950–Jun 1955" = seq_len(air_cut),
  "Jul 1955–Dec 1960" = (air_cut + 1):nrow(air_candidates)
)
air_portion_summary <- do.call(rbind, lapply(
  seq_len(ncol(air_candidates)), function(j) {
    do.call(rbind, lapply(names(air_portions), function(period) {
      x <- as.numeric(air_candidates[air_portions[[period]], j])
      r <- as.numeric(acf(x, lag.max = 12, plot = FALSE)$acf)
      data.frame(Candidate = colnames(air_candidates)[j], Period = period,
                 n = length(x), Mean = mean(x), SD = sd(x),
                 ACF1 = r[2], ACF12 = r[13])
    }))
  }
))
knitr::kable(air_portion_summary, digits = 3, row.names = FALSE,
             col.names = c("Candidate", "Period", "n", "Mean", "SD",
                           "ACF(1)", "ACF(12)"),
             caption = "Descriptive summaries within two chronological portions.")
```

The selected lags give a compact comparison; the following plots show
the intermediate lags as well.

```{r}
#| label: fig-wb8-air-portion-acfs
#| code-fold: true
#| code-summary: "R code for the ACF comparison"
#| fig-height: 7.6
#| fig-cap: "Sample ACFs computed separately in the two chronological portions, with all lags measured in months. Each estimate uses only its own portion, so differences between panels include sampling variation. Real data: datasets::AirPassengers."
par(mfrow = c(3, 2), mar = c(3.5, 4, 3.5, 1),
    mgp = c(2.2, 0.7, 0), bg = "white")
for (j in seq_len(ncol(air_candidates))) {
  for (period in names(air_portions)) {
    x <- as.numeric(air_candidates[air_portions[[period]], j])
    acf(x, lag.max = 18, ci = 0, ylim = c(-0.6, 1),
        xlab = "Lag (months)",
        main = paste(colnames(air_candidates)[j],
                     if (period == names(air_portions)[1]) "earlier" else "later",
                     sep = ": "),
        cex.main = 0.85)
  }
}
```

For the fitted-mean candidate, we can also compare calendar-month
residual means within each portion. The regression sets each month's
residual mean to zero over its full 144-month fitting window. It does not
impose that constraint separately in these two portions.

```{r}
#| label: fig-wb8-air-calendar-residuals
#| code-fold: true
#| code-summary: "R code for the calendar profiles"
#| fig-height: 3.6
#| fig-cap: "Calendar-month means of fitted-mean residuals in the two chronological portions. The different profiles motivate checking the fixed seasonal shape; each plotted mean uses only five or six observations. Real data: datasets::AirPassengers."
air_calendar_profiles <- sapply(air_portions, function(ii) {
  tapply(as.numeric(air_candidates[ii, 1]),
         cycle(air_candidates)[ii], mean)
})
par(mfrow = c(1, 1), mar = c(4, 4, 2, 1), bg = "white")
matplot(1:12, air_calendar_profiles, type = "b", lty = c(1, 2),
        pch = c(16, 17), col = c("#00539B", "#AD5B00"), xaxt = "n",
        xlab = "Calendar month", ylab = "Mean log residual",
        main = "Calendar profiles after fitting one seasonal shape")
axis(1, at = 1:12, labels = month.abb)
abline(h = 0, col = "gray60")
legend("topleft", legend = names(air_portions), bty = "n",
       col = c("#00539B", "#AD5B00"), lty = c(1, 2), pch = c(16, 17),
       cex = 0.8)
```

**Remark 8.12 (A finite-record comparison).** The two portions need not
produce equal estimates under stationarity. They contain only about five
annual cycles each, and dependence further limits information about a
mean or covariance. These comparisons do not attach a rejection
probability to a discrepancy. A stationarity assumption also supplies
neither the ergodicity conditions for learning moments nor a CLT for
their estimators.

For the airline series, we choose the combined log differences
$$
V_t=(1-B)(1-B^{12})\log X_t
$$
as our starting point for a stationary dependence model. This series
measures month-to-month changes in annual log growth.

The time plot shows fewer sustained movements than the fitted-mean
residuals or annual log growth. Its ACF has a relatively simple pattern,
with the most visible dependence at short lags and around twelve months.
We retain this dependence for modeling. Nonzero autocorrelations alone
give no reason to take another difference.

The choice leaves a concern about stability. The sample standard
deviation falls from `r sprintf("%.3f", air_portion_summary$SD[5])`
in the earlier portion to `r sprintf("%.3f", air_portion_summary$SD[6])`
in the later portion; the lag-12 sample correlation changes from
`r sprintf("%.3f", air_portion_summary$ACF12[5])` to
`r sprintf("%.3f", air_portion_summary$ACF12[6])`.
The quieter later stretch is visible in the time plot. We therefore
use stationarity as a working approximation whose adequacy remains to
be assessed.

The comparison does not show that the extra ordinary difference was
necessary. The slower movements in annual log growth or the fitted-mean
residuals could themselves be stationary fluctuations. Seasonal
differences remain a candidate for modeling annual growth directly;
the fitted residuals remain a candidate for modeling deviations from
the specified trend and monthly means. Exercise 8.5 asks for a defended
choice using the same evidence.

## What Remains to Be Modeled?

For the chosen airline series, the short-lag and annual-lag correlations
are now features for a dependence model to explain. Once a stationary
approximation is adopted, fitting that dependence can support a further
filter aimed at recovering white noise.

To illustrate this step on a simpler model, suppose a stationary series
$(Y_t)$ with mean $\mu$ follows the causal AR(1) model
$$
Y_t-\mu=\phi(Y_{t-1}-\mu)+Z_t,\qquad |\phi|<1,
$$
where $(Z_t)\sim\mathrm{WN}(0,\sigma^2)$. For any fixed coefficient
$\psi$, consider the residual filter
$$
R_t(\psi):=(Y_t-\mu)-\psi(Y_{t-1}-\mu)
=Z_t+(\phi-\psi)(Y_{t-1}-\mu).
$$
Proposition 8.5 makes $R_t(\psi)$ weakly stationary for every fixed
$\psi$. Choosing $\psi=\phi$ recovers $Z_t$; other choices can leave
autocorrelation. Here coefficient accuracy concerns recovery of the
white-noise sequence. With trend subtraction, a wrong fixed slope can
leave a changing mean.

Both trend removal and dependence filtering can require fitting unknown
parameters. Replacing $\mu$ and $\phi$ by estimates gives residuals intended
to estimate the white-noise sequence; they need not have its population
properties exactly. The ARMA material develops these dependence models
and their residuals more generally.

To return from changes to levels, we retain the initial values needed
for reconstruction. For a first-difference series
$D_t:=X_t-X_{t-1}$, one level anchor suffices:
$$
X_{T+h}=X_T+\sum_{j=1}^{h}D_{T+j}.
$$
Seasonal differences require one anchor for each seasonal subsequence.
For log differences, we reconstruct the log levels and then exponentiate.

All adjustments in the airline comparison were estimated or computed for
describing the observed record. For forecast evaluation, estimated trends
and seasonal effects must be fitted using only the training period.
The forecasting material develops reconstruction of predictions and their
uncertainty on the original scale.

## Fixed Transformations of a Stationary Process

The following result concerns a fixed rule applied to a strictly
stationary input. 

**Proposition 8.13 (Preservation of strict stationarity).** Suppose
$(X_t)$ is strictly stationary. For a fixed integer $q\ge0$ and a fixed
function $g:\R^{q+1}\to\R$ for which the random variables are well
defined, let
$$
Y_t:=g(X_t,X_{t-1},\ldots,X_{t-q}).
$$
Then $(Y_t)$ is strictly stationary.

*Proof.* For any $t_1,\ldots,t_k$ and shift $r$, the vector
$(Y_{t_1+r},\ldots,Y_{t_k+r})$ is obtained by applying the same rule to
the shifted collection of input coordinates. Strict stationarity makes
the distribution of that input collection independent of $r$, hence
also the distribution of the output vector. $\square$

To conclude weak stationarity of this output, finite second moments are
also needed. Neither a general nonlinear transformation of a merely
weakly stationary input nor a fitted rule depending on the entire
observed record is covered by Proposition 8.13.

## Sources

Shumway and Stoffer (2025), *Time Series Analysis and Its Applications:
With R Examples*, 5th ed., Sections 2.2–2.3, develops detrending,
differencing, and smoothing through data examples. The moment calculations
here use the course definitions of weak stationarity and the ACVF.
The airline data are `datasets::AirPassengers`; the temperature exercise
uses `astsa::gtemp.month`.

Hyndman and Athanasopoulos, *Forecasting: Principles and Practice*, 3rd ed.,
[Section 7.4](https://otexts.com/fpp3/useful-predictors.html), discusses
seasonal indicators and harmonic regression;
[Section 9.4](https://otexts.com/fpp3/MA.html) treats infinite-past
invertibility for moving-average models.

## Exercises

::: {.exercise}
**Exercise 8.1 (A filter is not automatically a stationarizer).** Let
$(X_t)$ be weakly stationary with mean $\mu$ and ACVF $\gamma_X$, and define
$Y_t=\sum_{j=0}^q a_jX_{t-j}$ for fixed coefficients $a_0,\ldots,a_q$.

*Method and report.* Work by hand. Derive the two filter moments and the
random-walk variance, then explain what each calculation establishes.

(a) Derive $\E[Y_t]$ and $\gamma_Y(h)$ directly.

(b) Specialize the result to $Y_t=X_t-X_{t-1}$. Express $\gamma_Y(0)$ and
$\gamma_Y(1)$ in terms of $\gamma_X$.

(c) Now let $X_0=0$ and $X_t=\sum_{j=1}^tZ_j$ for $t\geq1$, where
$(Z_t)$ is i.i.d. with mean zero and variance $\sigma^2>0$. For
$M_t=(X_t+X_{t-1})/2$, show that
$$
\Var(M_t)=\left(t-\frac34\right)\sigma^2.
$$

(d) Compare $M_t$ with $\nabla X_t$. Explain why a smoother time plot need
not describe a stationary process.
:::

::: {.exercise}
**Exercise 8.2 (Two trends, two treatments).** Generate $T=300$ observations
with seed 542 from
$$
X_t^{(D)}=0.02t+Z_t,
\qquad
X_t^{(S)}=X_{t-1}^{(S)}+0.02+U_t,
\qquad X_0^{(S)}=0,
$$
where $(Z_t)$ and $(U_t)$ are independent i.i.d. $N(0,1)$ sequences.

*Method and report.* Use R for the plots and descriptive summaries. Work
by hand for part (c). Keep the difference between a population component
and fitted residual data explicit.

(a) Plot both observed series. For each, fit and subtract a linear trend,
and also construct its first differences. Compare the four outputs on the
common indices $t=2,\ldots,300$ using time plots and ACFs through lag 20
observations.

(b) Split each output into its first 149 and last 150 values. Report the
mean, variance with divisor equal to the block length, and lag-one sample
autocorrelation within each block. Treat the comparisons as descriptive.

(c) Derive the mean and ACVF of $\nabla X_t^{(D)}$. Explain its lag-one
correlation without fitting a time-series model.

(d) For each model, does subtracting its specified mean $0.02t$ give a
stationary process? Does first differencing? When both work, explain
what each output measures and which retains the fluctuations around
the trend. Why do fitted residuals from $X^{(D)}$ differ from $Z_t$?
Why does their overall mean being zero not establish stability across time?
:::

::: {.exercise}
**Exercise 8.3 (Two seasonal structures).** Compare two models for monthly
observations. In the first, the observations fluctuate around a fixed
mean that repeats each year. In the second, each observation equals the
value in the same month of the preceding year plus a new shock.
Both models use the same i.i.d. $N(0,1)$ sequence $(Z_t)_{t\geq1}$.

**Fixed seasonal mean:**
$$
X_t=S_t+Z_t,\qquad S_t=2\cos(2\pi t/12),\qquad t\geq1.
$$

**Seasonal accumulation:**
$$
W_t=W_{t-12}+Z_t,\qquad t\geq1,
$$
with the twelve values before the observed record set to
$W_{-11}=\cdots=W_0=0$.

*Method and report.* Derive the population adjustments by hand, then
compare fitted seasonal residuals with seasonal differences in simulated
records.

(a) For $t\geq13$, derive $X_t-S_t$ and $\nabla_{12}X_t$. Calculate the variance and
lag-12 correlation of the latter.

(b) Write $t=12k+r$ with $k\in\{0,1,\ldots\}$ and $r\in\{1,\ldots,12\}$.
Find $\Var(W_t)$ and $\nabla_{12}W_t$. Why can subtracting a fixed
seasonal curve not make $(W_t)$ weakly stationary?

(c) Simulate both series for 30 years using seed 542. For each record, fit
$$
m(t)=a+c\cos(2\pi t/12)+d\sin(2\pi t/12)
$$
by OLS, treating its coefficients as unknown. Compare the fitted
residuals and seasonal differences using time plots and ACFs through
lag 36 months, on the common indices $13,\ldots,360$.

(d) For each model, which adjustment recovers or estimates $Z_t$?
Explain using parts (a)--(b).
Does a stationary output have to be uncorrelated?
:::

::: {.exercise}
**Exercise 8.4 (An estimated annual mean in global temperature).** The
data frame `astsa::gtemp.month` contains monthly global average surface
temperatures in degrees Celsius for 1975--2023. Its rows are months and
its columns are years in reverse chronological order.

```{r}
#| eval: false
data("gtemp.month", package = "astsa")
g <- gtemp.month[, ncol(gtemp.month):1]
x <- ts(as.vector(as.matrix(g)), start = c(1975, 1), frequency = 12)
t <- seq_along(x)
```

*Method and report.* Use R to fit the stated adjustment and inspect its
output. Report descriptive summaries; do not use the usual independent-error
standard errors printed by `lm`.

(a) Plot the observations and fit
$$
X_t=a+bt+c\cos(2\pi t/12)+d\sin(2\pi t/12)+U_t.
$$
As a working model to assess, suppose $(U_t)$ has mean zero and is weakly
stationary. Explain why fixing the annual frequency makes this a linear regression.
Report the fitted trend in degrees Celsius per decade.

(b) Define
$\widehat S_t=\hat c\cos(2\pi t/12)+\hat d\sin(2\pi t/12)$.
Plot both $x_t-\widehat S_t$ and
$e_t=x_t-\hat a-\hat bt-\widehat S_t$. What does each adjustment remove,
and in what units is each output measured?

(c) Compare the residuals in 1975--1998 and 1999--2023 using their
within-period means, variances, and sample ACFs through lag 24 months.
Report the sample correlations at lags 1 and 12 months.

(d) Identify one feature still requiring attention before adopting a
stationary approximation. Explain why remaining serial dependence is not,
by itself, a reason to difference again. Distinguish the proposed process
$U_t$ from the fitted residual data.
:::

::: {.exercise}
**Exercise 8.5 (Choose a process to model).** Use
`datasets::AirPassengers`, monthly international airline passenger totals
in thousands from January 1949 through December 1960. Let
$\ell_t=\log x_t$.

*Method and report.* Reconstruct the three candidates in R. Compare them
on the common February 1950--December 1960 window. End with a short
recommendation that records a remaining limitation.

(a) Fit $\ell_t$ to an intercept, a linear time trend, and calendar-month
indicators, omitting one month indicator to avoid redundancy. Construct
its residuals, $\ell_t-\ell_{t-12}$, and
$(\ell_t-\ell_{t-12})-(\ell_{t-1}-\ell_{t-13})$.
State the proposed structure behind each adjustment and what its output
measures.

(b) Plot each candidate and its ACF through lag 24 months. Compare the
first 65 months (February 1950--June 1955) with the remaining 66 months
(July 1955--December 1960). Report the mean, variance with divisor equal
to the block length, and sample correlations at lags 1 and 12 months.

(c) Recommend one candidate for a stationary model, or explain why none
is yet satisfactory. State the quantity that would be modeled and refer
to level, spread, and dependence across time. Do not select a candidate
solely because its sample ACF is smallest.

(d) State one further check that could change the recommendation.
Explain why these plots cannot determine whether the original series has
a deterministic trend or accumulated shocks.
:::