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.
Figure 1: Monthly international airline passenger totals, 1949–1960. The three annual profiles compare the same calendar months at different levels. Real data: datasets::AirPassengers.
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 code for the plots
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)
Figure 2: 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.
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\).
\(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.
Simulation and adjustment code
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_trendtrend_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 innames(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")}
Figure 3: 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.
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.
Fit lines to the same two records
par(mfrow =c(2, 2), mar =c(3.5, 4, 3, 1),mgp =c(2.2, 0.7, 0), bg ="gray92")for (name innames(trend_records)) { x <- trend_records[[name]] trend_fit <-lm(x ~ time_trend) known_remainder <- x - known_trendplot(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)}
Figure 4: 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.
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.
Fit straight and quadratic trends
x_quadratic <- x_deterministic +0.0001* time_trend^2straight_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^2true_fluctuation <- x_deterministic -0.02* time_trendresidual_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)
Figure 5: 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.
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.
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.
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.
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.
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 code for the candidate comparison
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 inseq_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")}
Figure 6: 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.
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 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.")
Descriptive summaries within two chronological portions.
Candidate
Period
n
Mean
SD
ACF(1)
ACF(12)
Fitted mean residual
Feb 1950–Jun 1955
65
0.006
0.053
0.725
0.017
Fitted mean residual
Jul 1955–Dec 1960
66
0.003
0.058
0.811
0.458
Seasonal difference
Feb 1950–Jun 1955
65
0.131
0.069
0.667
-0.391
Seasonal difference
Jul 1955–Dec 1960
66
0.110
0.050
0.756
-0.119
Both differences
Feb 1950–Jun 1955
65
0.002
0.056
-0.323
-0.454
Both differences
Jul 1955–Dec 1960
66
-0.002
0.033
-0.404
-0.196
The selected lags give a compact comparison; the following plots show the intermediate lags as well.
Figure 7: 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.
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 code for the calendar profiles
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)
Figure 8: 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.
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 0.056 in the earlier portion to 0.033 in the later portion; the lag-12 sample correlation changes from -0.454 to -0.196. 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, discusses seasonal indicators and harmonic regression; Section 9.4 treats infinite-past invertibility for moving-average models.
Exercises
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.
Derive \(\E[Y_t]\) and \(\gamma_Y(h)\) directly.
Specialize the result to \(Y_t=X_t-X_{t-1}\). Express \(\gamma_Y(0)\) and \(\gamma_Y(1)\) in terms of \(\gamma_X\).
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.
\]
Compare \(M_t\) with \(\nabla X_t\). Explain why a smoother time plot need not describe a stationary process.
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.
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.
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.
Derive the mean and ACVF of \(\nabla X_t^{(D)}\). Explain its lag-one correlation without fitting a time-series model.
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 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}\).
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.
For \(t\geq13\), derive \(X_t-S_t\) and \(\nabla_{12}X_t\). Calculate the variance and lag-12 correlation of the latter.
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?
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\).
For each model, which adjustment recovers or estimates \(Z_t\)? Explain using parts (a)–(b). Does a stationary output have to be uncorrelated?
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.
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.
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.
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?
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.
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 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.
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.
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.
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.
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.
