Time Series Analysis

The following notes come from the course Time Series Analysis, taught at ETH Zurich by Professor Fadoua Balabdaoui.

Feel free to reach out to me with any comments, clarifications, or corrections.

1. Introduction to Time Series

Time Series Definition

A time series is a stochastic process indexed by time $t \in T \subseteq R$ .

Examples (they are all continuous!):

$T = [0, 1]$ : Brownian bridge from $(0, 0)$ to $(1, 0)$ , where the pair notation denotes $(time, value) = (t, X_{t})$
$T = [0, + \infty]$ : Brownian motion starting at $(0, 0)$
$T = R$ : two-sided Brownian motion
$T = N_{0}$ , or $T = Z$ : discrete random walk.

We use the notation ${X_{t}}_{t \in T}$ to denote a random variable and the notation ${x_{t}}_{t \in T}$ to denote realizations of that variable.

A time series model is either the specification of:

all joint distributions $f ({X_{t}}_{t \in T})$ (i.e., the model provides $Pr (X_{t_{1}} \leq x_{t_{1}}, \dots, X_{t_{k}} \leq x_{t_{k}})$ for all $k \in N$ , $t_{1}, \dots, t_{k}$ , and $x_{1}, \dots, x_{k} \in R$ )
all their means and covariances $(E (X_{t}), Cov (X_{t}, X_{p}))$ for all $t, p \in T$

White noise

The i.i.d. noise ${X_{t}}_{t \in T}$ is a sequence of independent and identically distributed random variables such that $E (X_{t}) = 0$ for all $t$ .
In general, we also assume that $E [X_{t}^{2}] = Var (X_{t}) < \infty$ for all $t$ . Note that in this case $E [X_{t}^{2}]$ does not depend on $t$ .

Rademacher Process

Consider ${X_{t}}_{t \in T}$ such that $Pr (X_{t} = - 1) = Pr (X_{t} = 1) = 1/2$ for all $t$ .

If ${B_{t}}_{t} \sim_{iid} Bernoulli (1/2)$ , then:

${X_{t}}_{t} = 2 {B_{t}}_{t} - 1$
$Var (X_{t}) = 4 Var (B_{t}) = 1$ (note that $X_{t}^{2} = 1$ and $E (X_{t}) = 0$ )

White noise

${X_{t}}_{t}$ is a sequence of centered and uncorrelated random variables. More specifically, $E (X_{t}) = 0$ , $E (X_{t}^{2}) < \infty$ , and $Cov (X_{t}, X_{s}) = 0$ for all $t \neq = s$ .

Note that i.i.d. noise such that $E (X_{t}^{2}) < \infty$ is also white noise.
The term "white" is used by analogy with white light. Indeed, white noise with constant variance $v^{2}$ has a flat spectral density given by $f (λ) = v^{2} /2 π$ , with $λ \in [- π, π]$ (this means that all time-series frequencies are present).

Random walk

Let ${X_{t}}_{t}$ be i.i.d. noise such that $E (X_{t}^{2}) < \infty$ for all $t \in N$ .

${S_{t}}_{t \in N_{0}}$ is a random walk if $S_{0} = 0$ and $S_{t} = \sum_{i = 1}^{t} X_{i}$ for all $t \geq 1$ .

If ${X_{t}}_{t \in N_{0}}$ is a Rademacher process, then ${S_{t}}_{t \in N_{0}}$ is called a symmetric random walk, with $E (S_{t}) = \sum_{i = 1}^{t} E (X_{i}) = 0$ .

Time Series with trend

Let ${Y_{t}}_{t}$ be a zero-mean time series and ${m_{t}}_{t}$ be a deterministic (non-random) sequence.
Consider $X_{t} = Y_{t} + m_{t}$ for all $t \in T$ .

Then $E [X_{t}] = m_{t}$ is called the trend of ${X_{t}}_{t}$ .

Time Series with seasonality

Given a zero-mean time series ${Y_{t}}_{t}$ and a function $t \mapsto p_{t}$ such that $p_{t}$ is periodic with known period $d$ (i.e., $p_{t + d} = p_{t}$ ).

Consider the time series $X_{t} = Y_{t} + p_{t}$ for all $t \in T$ . Then ${p_{t}}_{t}$ is called the seasonal component of ${Y_{t}}_{t}$ .

For example, we could have the seasonal component in the following forms:

p_{t} = a_{0} + j = 1 \sum k (a_{j} cos (\frac{2 π α _{j} t}{d}) + b_{j} sin (\frac{2 π β _{j} t}{d}))

with known integers $k \geq 1$ , $α_{j}, β_{j} > 0$ , and unknown $a_{0}, a_{j}, b_{j}$ for all $1 \leq j \leq k$ .

Time series with both trend and seasonal component

$X_{t} = m_{t} + p_{t} + Y_{t}, \forall t \in T$

with ${m_{t}}_{t}$ the trend, ${p_{t}}_{t}$ the seasonal component, and ${Y_{t}}_{t}$ a zero-mean time series.

Stationarity

A time series ${X_{t}}_{t}$ is said to be strictly (strongly) stationary if, for all $n \geq 1$ and $h \in Z$ , it holds that:

(X_{1}, \dots, X_{n}) = d (X_{1 + h}, \dots, X_{n + h})

i.e., the stochastic behavior of ${X_{t}}_{t}$ remains unchanged for any time shift.

Let ${X_{t}}_{t}$ be a time series.

The mean function of ${X_{t}}_{t}$ is defined as:

μ_{x} (t) = E (X_{t}), \forall t \in T

The covariance function of ${X_{t}}_{t}$ is defined as:

γ_{x} (r, p) = Cov (X_{r}, X_{p}) = = E [(X_{r} - μ_{x} (r)) (X_{p} - μ_{x} (p))] = E [X_{r} X_{p}] - μ_{x} (r) μ_{x} (p)

provided that $E [∣ X_{r} X_{p} ∣] < \infty$ and $E [X_{r}] < \infty$ for all $r, p$ .

On the other hand, ${X_{t}}_{t}$ is said to be (weakly) stationary if $E [X_{t}^{2}] < \infty$ for all $t$ , and:

$μ_{x} (t)$ is constant in $t$
$γ_{x} (t, t + h)$ depends only on $h$

1. Week 2

Weak Stationarity

A time series ${X_{t}}_{t \in Z}$ is said to be weakly stationary if:

the mean function is constant over time: $μ_{x} (t) = E (X_{t}) = μ \in R$
the covariance $γ_{X} (r, s) = Cov (X_{r}, X_{s})$ depends only on $h = s - r$

Intuition: Weak stationarity means a series’ average level and dependence pattern stay stable over time.

Note also that the covariance between two time points depends only on their time distance, not on their absolute values (i.e., not on their calendar dates). On the other hand, strict stationarity is a strong concept since the full distribution must be constant over time.

For example, financial prices often show a trend and are not stationary, so we usually work with log returns, such as $r_{t} = lo g P_{t} - lo g P_{t - 1}$ , which can be closer to weak stationarity and supports tools like ARMA models, autocorrelation analysis, and covariance-based risk estimates.

However, financial returns are not reliable in this sense since volatility clustering, market regimes, and structural breaks can strongly change their distribution.

Given the definition, we can understand why the strict stationarity of ${X_{t}}$ implies that all the $X_{t}$ have the same marginal distribution.

To show this, set $n = 1$ . By the definition of strict stationarity, for any lag $h \in Z$ we have:

$X_{1} = d X_{1 + h}$

Furthermore, since ${1 + h : h \in Z} = Z$ , it follows that $X_{1} = d X_{t}$ for all $t \in Z$ .

On the other hand, we can now spend some time understanding what a moment means in this context.

In the table below, you can find a summary comparing the first and second moments for single and paired random variables.

W : (W_{1}, W_{2}) : 1st moment: E (W), 2nd moment: E (W^{2}) / Var (W) . 1st moments: E (W_{1}), E (W_{2}), 2nd moments: E (W_{1}^{2}), E (W_{2}^{2}), E (W_{1} W_{2}) / Var (W_{1}), Var (W_{2}), Cov (W_{1}, W_{2}) .

We also recall that moments can be generated by the moment-generating function, defined as:

$M_{W} (t) = E [e^{t W}]$

M_{W_{1}, W_{2}} (s, t) = E [e^{s W_{1} + t W_{2}}]

Autocovariance and autocorrelation functions

Let ${X_{t}}_{t}$ be a stationary time series. The Autocovariance Function (ACVF) at lag $h \in Z$ is defined as:

$γ_{X} (h) = Cov (X_{t}, X_{t + h})$

Similarly, the Autocorrelation Function (ACF) at lag $h$ is defined as:

ρ_{X} (h) = \frac{Cov ( X _{t} , X _{t + h} )}{Var ( X _{t} ) Var ( X _{t + h} )} = \frac{γ _{X} ( h )}{γ _{X} ( 0 )}

The second equality holds because, if ${X_{t}}_{t}$ is stationary, the following hold:

$Var (X_{t}) = Cov (X_{t}, X_{t}) = γ_{X} (0), \forall t \in Z .$

Var (X_{t}) Var (X_{t + h}) = γ_{X} (0) γ_{X} (0) = γ_{X} (0) .

Key Properties of ACVF and ACF

Boundedness: by the Cauchy–Schwarz inequality, we have that:

∣ Cov (X_{t}, X_{t + h}) ∣ \leq Var (X_{t}) Var (X_{t + h}) ⟺ ∣ γ_{X} (h) ∣ \leq ∣ γ_{X} (0) ∣ ⟺ ∣ ρ_{X} (h) ∣ \leq 1.

Recall that the Cauchy–Schwarz inequality states that, for random variables $X$ and $Y$ with finite second moments:

∣ E [X Y] ∣ \leq E [X^{2}] E [Y^{2}] .

Symmetry: since $γ_{X}$ and $ρ_{X}$ are even functions (i.e., $f (x) = f (- x)$ for all $x$ ), for any $h \in Z$ we have:

γ_{X} (h) = Cov (X_{t}, X_{t + h}) = Cov (X_{t + h}, X_{t}) .

If we then set $t^{'} = t + h$ , we have

Cov (X_{t^{'}}, X_{t^{'} - h}) = γ_{X} (- h), so γ_{X} (- h) = γ_{X} (h) .

ρ_{X} (- h) = \frac{γ _{X} ( - h )}{γ _{X} ( 0 )} = \frac{γ _{X} ( h )}{γ _{X} ( 0 )} = ρ_{X} (h) .

I.I.D. Noise Process

Let ${X_{t}}_{t}$ be an i.i.d. noise process such that $E [X_{t}] = 0$ and $E [X_{t}^{2}] = Var (X_{t}) = σ^{2}$ . In this case, we use the notation ${X_{t}}_{t} \sim IID (0, σ^{2})$ .

For any $h \in Z$ , we also have

γ_{X} (t, t + h) = Cov (X_{t}, X_{t + h}) = {Var (X_{t}) = σ^{2}, 0, if h = 0, if h \neq = 0.

The zero covariance for $h \neq = 0$ is implied by independence.

Furthermore, $μ_{X} (t) = 0$ for all $t$ , so the process ${X_{t}}_{t}$ is stationary, with ACVF and ACF given by:

$γ_{X} (h) = σ^{2} 1_{h = 0}, ρ_{X} (h) = 1_{h = 0} .$

Also, note that an i.i.d. noise process is also strictly stationary.

White Noise Process

On the other hand, let ${X_{t}}_{t}$ be a white noise process with variance $σ^{2}$ . In this case, we say ${X_{t}}_{t} \sim WN (0, σ^{2})$ if the process satisfies the following:

$E [X_{t}] = 0, E [X_{t}^{2}] = σ^{2} < \infty, Cov (X_{t}, X_{s}) = 0 \forall t \neq = s .$

Note also that the process is stationary with ACVF $γ_{X} (h) = σ^{2} 1_{h = 0}$ .

Note that in the case of i.i.d. noise, the $X_{t}$ are independent and identically distributed, while for white noise we just require uncorrelatedness, zero mean, and finite constant variance (independence is not required!)

Example: Random Walk Process

We recall that ${S_{t}}_{t}$ is a random walk if $S_{0} = 0$ and $S_{t} = \sum_{i = 1}^{t} X_{i}$ , where ${X_{t}}_{t} \sim IID (0, σ^{2})$ .

Also, recall that a random walk has mean zero, $Var (S_{0}) = 0$ , and, for all $t \geq 1$ , we have:

Var (S_{t}) = Var (i = 1 \sum t X_{i}) = i = 1 \sum t Var (X_{i}) = t σ^{2} .

We immediately notice that, since $Var (S_{t})$ depends on time $t$ , the process is not stationary.

We now spend some time computing the autocovariance for a random walk. Consider $t \geq 1$ and $h \geq 1$ , so we have:

γ_{S} (t, t + h) = Cov (S_{t}, S_{t + h}) = Cov (i = 1 \sum t X_{i}, i = 1 \sum t X_{i} + i = t + 1 \sum t + h X_{i}) = Cov (i = 1 \sum t X_{i}, i = 1 \sum t X_{i}) + Cov (i = 1 \sum t X_{i}, i = t + 1 \sum t + h X_{i}) = Var (S_{t}) + 0 = t σ^{2} .

The last equality holds because

Cov (i = 1 \sum t X_{i}, j = t + 1 \sum t + h X_{j}) = i = 1 \sum t j = t + 1 \sum t + h Cov (X_{i}, X_{j}) = 0

by independence.

If we also consider the special case $t = 0$ , for which $γ_{S} (0, h) = 0$ , we can conclude that, for all $t, h \in N_{0}$ , we have

$γ_{S} (t, t + h) = t σ^{2} .$

Secondly, we can do the same for the autocorrelation of a random walk. For all $h \in N_{0}$ and $t \geq 1$ , we can compute:

ρ_{S} (t, t + h) = \frac{Cov ( S _{t} , S _{t + h} )}{Var ( S _{t} ) Var ( S _{t + h} )} = \frac{t σ ^{2}}{t σ ^{2} ( t + h ) σ ^{2}} = \frac{t}{t + h} .

We can now characterize the asymptotic behavior of random walk autocorrelation.

If we fix lag $h$ and let $t \to \infty$ , we have stronger correlation as the length of shared history grows:

$lim_{t \to \infty} ρ_{S} (t, t + h) = 1.$

If we fix $t$ and let $h \to \infty$ , we have weaker correlation for larger time gaps:

$lim_{h \to \infty} ρ_{S} (t, t + h) = 0.$

First-Order Moving Average Process | MA(1)

Let ${Z_{t}}_{t} \sim WN (0, σ^{2})$ and $θ \in R$ . Consider the time series ${X_{t}}_{t}$ defined as:

$X_{t} = Z_{t} + θ Z_{t - 1}$

Intuition: An MA(1) process combines today's random shock $Z_{t}$ with the previous period's shock $Z_{t - 1}$ , scaled by $θ$ . The sign of $θ$ determines whether consecutive shocks reinforce or offset each other.

Given this, we define the mean function as:

μ_{X} (t) = E [X_{t}] = E [Z_{t} + θ Z_{t - 1}] = E [Z_{t}] + θ E [Z_{t - 1}] = 0, \forall t \in Z .

Secondly, we can define the autocovariance function (ACVF) of an MA(1) process for all $h \in Z$ :

γ_{X} (t, t + h) = Cov (X_{t}, X_{t + h}) = Cov (Z_{t} + θ Z_{t - 1}, Z_{t + h} + θ Z_{t + h - 1}) = Cov (Z_{t}, Z_{t + h}) + θ Cov (Z_{t}, Z_{t + h - 1}) + θ Cov (Z_{t - 1}, Z_{t + h}) + θ^{2} Cov (Z_{t - 1}, Z_{t + h - 1}) = σ^{2} 1_{h = 0} + θ σ^{2} 1_{h = 1} + θ σ^{2} 1_{h = - 1} + θ^{2} σ^{2} 1_{h = 0} = σ^{2} (1 + θ^{2}) 1_{h = 0} + θ σ^{2} (1_{h = 1} + 1_{h = - 1}) = σ^{2} (1 + θ^{2}) 1_{h = 0} + θ σ^{2} 1_{∣ h ∣ = 1} .

The equalities follow from the fact that white noise has zero covariance across different times ( $Cov (Z_{i}, Z_{j}) = 0$ ) and variance $σ^{2}$ at each time. Each indicator marks when two shocks in the expansion are the same shock.

So, since this depends only on the lag $h$ , the MA(1) process is stationary, with the following ACVF:

γ_{X} (h) = ⎩ ⎨ ⎧ σ^{2} (1 + θ^{2}), θ σ^{2}, 0, if h = 0, if h \in {- 1, 1}, otherwise .

Similarly, we can compute the autocorrelation function of MA(1), given by:

ρ_{X} (h) = \frac{γ _{X} ( h )}{γ _{X} ( 0 )} = ⎩ ⎨ ⎧ 1, \frac{θ}{1 + θ ^{2}}, 0, if h = 0, if h \in {- 1, 1}, otherwise .

First-Order Autoregressive Process | AR(1)

Let ${Z_{t}}_{t} \sim WN (0, σ^{2})$ and consider a stationary time series ${X_{t}}_{t}$ satisfying:

$X_{t} = Φ X_{t - 1} + Z_{t}, t \in Z$

where $Φ \in (- 1, 1)$ and $Cov (X_{s}, Z_{t}) = 0$ for all $s < t$ .

Intuition: The current value is the previous value, scaled by $Φ$ , plus a shock. $Φ$ controls how much of the past carries forward. Note that AR(1) is useful when past values help predict future values, while MA(1) is useful when the dependence comes from recent shocks; these processes will soon be combined in ARMA models.

The mean function of ${X_{t}}_{t}$ is then:

Since the time series is stationary, we have that $μ_{X} (t) = E [X_{t}] = μ \in R$ . Taking expectations on both sides yields:

E [X_{t}] = Φ E [X_{t - 1}] + E [Z_{t}] ⟺ μ = Φ μ + 0 ⟺ μ (1 - Φ) = 0.

But since $Φ \in (- 1, 1)$ , we cannot have $Φ = 1$ , so it follows that $μ = 0$ .

As before, the autocovariance function for AR(1) can be computed as follows for $h \geq 1$ :

γ_{X} (h) = Cov (X_{t - h}, X_{t}) = Cov (X_{t - h}, Φ X_{t - 1} + Z_{t}) = Φ Cov (X_{t - h}, X_{t - 1}) + Cov (X_{t - h}, Z_{t}) .

Since $t - h < t$ , according to the condition in the AR(1) definition, $Cov (X_{t - h}, Z_{t}) = 0$ , so we write:

γ_{X} (h) = Φ γ_{X} (h - 1), ⟹ γ_{X} (h) = Φ^{h} γ_{X} (0), \forall h \geq 1.

However, we still have to compute $γ_{X} (0)$ , which in this case is:

γ_{X} (0) = Var (X_{t}) = Var (Φ X_{t - 1} + Z_{t}) = Φ^{2} Var (X_{t - 1}) + 2Φ Cov (X_{t - 1}, Z_{t}) + Var (Z_{t}) .

So we can conclude that $γ_{X} (0) = Φ^{2} γ_{X} (0) + σ^{2}$ , which gives $γ_{X} (0) = \frac{σ ^{2}}{1 - Φ ^{2}}$ .

And, since the function $h \to γ_{X} (h)$ is even, we finally state:

$γ_{X} (h) = \frac{σ ^{2}}{1 - Φ ^{2}} Φ^{∣ h ∣}, \forall h \in Z .$

Similarly, the autocorrelation function of AR(1) is simply given by:

$ρ_{X} (h) = \frac{γ _{X} ( h )}{γ _{X} ( 0 )} = Φ^{∣ h ∣}, \forall h \in Z .$

3. Sample Estimation and Testing

We now consider $x_{1}, \dots, x_{n}$ to be observations from a stationary time series with:

sample mean $\overset{x}{ˉ} = \frac{1}{n} \sum_{t = 1}^{n} x_{t}$
sample ACVF at lag $h$ given by: $\overset{γ}{^} (h) = \frac{1}{n} \sum_{t = 1}^{n - ∣ h ∣} (x_{t + ∣ h ∣} - \overset{x}{ˉ}) (x_{t} - \overset{x}{ˉ}), \forall∣ h ∣ \leq n - 1$
sample ACF at lag $h$ given by: $\overset{ρ}{^} (h) = \frac{γ ^ ( h )}{γ ^ ( 0 )} = \frac{\sum _{t = 1}^{n - ∣ h ∣} ( x _{t + ∣ h ∣} - x ˉ ) ( x _{t} - x ˉ )}{\sum _{t = 1}^{n} ( x _{t} - x ˉ ) ^{2}}, \forall∣ h ∣ \leq n - 1$

Note that we divide by $n$ and not by $n - ∣ h ∣$ to obtain a positive semidefinite matrix.

In particular, let $y_{t} = x_{t} - \overset{x}{ˉ}$ be the centered observation and define, for $h \geq 0$ :

$\overset{γ}{^} (h) = \frac{1}{n} \sum_{t = 1}^{n - h} y_{t} y_{t + h} .$

Then build the Toeplitz matrix $\hat{Γ}$ by setting $\hat{Γ}_{ij} = \overset{γ}{^} (∣ i - j ∣)$ .

Recall also that "positive semidefinite" means that, for all vectors $a$ , we have:

$a^{T} \hat{Γ} a \geq 0.$

To see the reason, imagine the data $y_{t}$ outside the observed range (i.e., we only observed $(y_{1}, \dots, y_{n})$ , so define $y_{t} = 0$ for all other integer times). The quadratic form can be written as

$a^{T} \hat{Γ} a = \frac{1}{n} \sum_{k \in Z} (\sum_{i = 1}^{n} a_{i} y_{k - i})^{2},$

which is nonnegative because it is a sum of squares.

Intuition: Using $n$ as the denominator preserves the structure and allows the matrix to represent squared combinations of the data.

Note also that, for all time series such that the seasonal component $s_{t} = 0$ , we can consider the sample estimation of the trend as

$x_{t} = m_{t} + y_{t}, t \in Z, with E [Y_{t}] = 0 for all t .$

Estimation with Two-Sided Moving Average

We now spend some time presenting two different estimation approaches. Firstly, we begin with the two-sided moving average approach.

Suppose there exist $a, b \in R$ and an integer $q \geq 1$ such that, near $t$ , the following approximation for the trend term $m_{t}$ holds (i.e., a two-sided moving average), assuming the underlying model is $X_{t} = m_{t} + Y_{t}$ :

m_{t^{'}} \approx a + b t^{'}

for all times $t^{'}$ in the interval $[t - q, t + q]$ . Note that the trend does not have to be a line globally, just to look linear on a small window of width $2 q + 1$ around $t$ .

The estimate then averages the $2 q + 1$ observations centered on $t$ and is given by:

\overset{m}{^}_{t} = \frac{1}{2 q + 1} i = t - q \sum t + q X_{i} = \frac{1}{2 q + 1} j = - q \sum q X_{t - j}

If we substitute $X_{t - j} = m_{t - j} + Y_{t - j}$ and split the sum, we get:

\overset{m}{^}_{t} = signal part \frac{1}{2 q + 1} j = - q \sum q m_{t - j} + noise part \frac{1}{2 q + 1} j = - q \sum q Y_{t - j}

If we use the local linearity approximation, i.e., $m_{t - j} \approx a + b (t - j)$ , we get:

j = - q \sum q (a + b (t - j)) = j = - q \sum q (a + b t) - b j = - q \sum q j = (2 q + 1) (a + b t)

We notice now that the offsets cancel, i.e., $\sum_{j = - q}^{q} j = 0$ , so averaging preserves the trend's value at the center while averaging the noise.

The condition $q + 1 \leq t \leq n - q$ just says the whole window must fit within the observed data, from time $1$ to $n$ .

So we can conclude that:

\overset{m}{^}_{t} \approx a + b t + \frac{1}{2 q + 1} j = - q \sum q Y_{t - j} \approx m_{t} + \frac{1}{2 q + 1} j = - q \sum q Y_{t - j}

That translates into saying that $\overset{m}{^}_{t}$ equals the true trend plus a small, averaged-out noise term.

Consider also that:

Large $q$ averages out more noise, so the estimate is smoother, but the window should be small enough for the trend to be linear.
Small $q$ , on the other hand, keeps the linearity assumption more plausible but leaves more noise (and this is the famous bias-variance trade-off!)

Furthermore, for a large $q$ , we can apply the Weak Law of Large Numbers (WLLN) to conclude that:

\frac{1}{N} i = 1 \sum N Y_{i} P 0 as N \to \infty

Intuition: for $N = 2 q + 1$ , as $q$ grows, the noise term is concentrated near $0$ .

Furthermore, for uncorrelated noise with variance $σ^{2}$ , Chebyshev's inequality allows us to be more precise and state that:

P (\frac{1}{2 q + 1} j = - q \sum q Y_{t - j} > ε) \leq \frac{σ ^{2}}{( 2 q + 1 ) ε ^{2}} \to 0

We now spend some more time on the bias-variance trade-off. Let's start by considering:

{Y_{t}}_{t} \sim WN (0, σ^{2}), Δ_{t^{'}} = m_{t^{'}} - (a + b t^{'}), t^{'} \in {t - q, \dots, t + q}

Note that we assumed ${Y_{t}}_{t}$ to be white noise, so it has mean $0$ , variance $σ^{2}$ , and no correlation between samples. Also, note that $Δ_{t^{'}}$ is the linearization error (measuring how far the true trend is from the approximation line $a + b t^{'}$ )

We now calculate the expectation of the estimator. The white-noise term vanishes because it has zero mean:

E (\overset{m}{^}_{t}) = \frac{1}{2 q + 1} j = - q \sum q E (X_{t - j}) = \frac{1}{2 q + 1} j = - q \sum q m_{t - j}

If we subtract and add the line $a + b t^{'}$ , we have:

E [\overset{m}{^}_{t}] = \frac{1}{2 q + 1} j = - q \sum q (m_{t - j} - (a + b (t - j))) + \frac{1}{2 q + 1} j = - q \sum q (a + b (t - j)) = \frac{1}{2 q + 1} j = - q \sum q Δ_{t - j} + a + b t

where the second equality holds since the average of the line equals $a + b t$ , as shown before.

We now define the bias of the estimator as:

Bias (q) = E (\overset{m}{^}_{t}) - m_{t} = \frac{1}{2 q + 1} j = - q \sum q Δ_{t - j} + a + b t - m_{t} = \frac{1}{2 q + 1} j = - q \sum q Δ_{t - j} - Δ_{t}

Intuition: the bias is the average linearization error over the window, minus the linearization error at the center. Therefore, if the trend is truly linear, both terms are $0$ and the estimator is unbiased.
On the other hand, if the trend is curved, the bias is non-zero, and it usually grows with $q$ (window size) because a wider window includes points where the line fits worse.

Furthermore, we can compute the mean squared error (MSE) using its standard decomposition:

MSE (q) = E [(\overset{m}{^}_{t} - m_{t})^{2}] = Bias (q)^{2} + Var (\overset{m}{^}_{t})

Lastly, if we compute its variance, we immediately notice that, since $m_{t}$ is deterministic, it only comes from the noise, where, since the errors are uncorrelated, the variance "simply" adds up:

Var (\frac{1}{2 q + 1} j = - q \sum q Y_{t - j}) = \frac{( 2 q + 1 ) σ ^{2}}{( 2 q + 1 ) ^{2}} = \frac{σ ^{2}}{2 q + 1}

It is still interesting to point out that the plain average defined before is a special case of a broader family of estimators. Indeed, so far we used one specific way to average the window, equal weights, but this is not the only choice:

\overset{m}{^}_{t} = j = - q \sum q \frac{1}{2 q + 1} X_{t - j}

For example, we could put more weight near the center or more at the edges. In general, we can then say that:

\overset{m}{^}_{t} = j = - q \sum q w_{j} X_{t - j}, j = - q \sum q w_{j} = 1, j = - q \sum q j w_{j} = 0

where each condition can be written as follows:

$\sum_{j = - q}^{q} w_{j} = 1$ implies each constant trend passes through unchanged.
$\sum_{j = - q}^{q} j w_{j} = 0$ implies a linear trend passes through unchanged.

Note also that any symmetric weights satisfy the second condition, and that $w_{j} = 1/ (2 q + 1)$ is the simplest choice.

Parametric Estimation via Least Squares

We now spend some time with parametric estimation. Consider a model $X_{t} = m_{t} + Y_{t}$ , where $Y_{t}$ is noise, with the parametric assumption that the trend belongs to a finite-dimensional family of functions. That is, we can consider a polynomial of degree $k$ such that $m_{t} = α_{0} + α_{1} t + \dots + α_{k} t^{k}$ , where the unknowns are the $k + 1$ coefficients.

For the least-squares estimator, we want to find the coefficients that make the curve as close as possible to the data under the squared-error sense, i.e., we want to find:

(\overset{α}{^}_{0}, \overset{α}{^}_{1}, \dots, \overset{α}{^}_{k}) = α_{0}, \dots, α_{k} arg min C (α_{0}, \dots, α_{k}),

C (α_{0}, \dots, α_{k}) = t = 1 \sum n (X_{t} - α_{0} - α_{1} t - \dots - α_{k} t^{k})^{2} .

Intuition: $C$ is the total squared vertical distance between the data points and the polynomial, and note that we square it since we want a smooth, convex function. We therefore want to minimize it by finding the unique point where the slope in every direction is $0$ ; this means the least-squares estimator is a stationary point of $C$ .

We can then set $\partial C / \partial α_{i} = 0$ , which gives the following linear system of $k + 1$ equations in $k + 1$ unknowns:

t = 1 \sum n t^{i} (X_{t} - j = 0 \sum k α_{j} t^{j}) = 0, i = 0, \dots, k .

Intuition: the residuals are orthogonal to every regressor $1, t, t^{2}, \dots, t^{k}$ .

We can proceed with illustrative examples.

If $k = 0$ (i.e., we have a constant trend), the equation is simply $\sum_{t = 1}^{n} (X_{t} - α_{0}) = 0$ , which gives $α_{0} = \overset{ˉ}{X}$ ; the best constant fit is the sample mean.

If $k = 1$ (i.e., we have a linear trend), the two equations are:

t = 1 \sum n (X_{t} - \overset{α}{^}_{0} - \overset{α}{^}_{1} t) = 0, t = 1 \sum n t (X_{t} - \overset{α}{^}_{0} - \overset{α}{^}_{1} t) = 0.

If we divide by $n$ and do some algebra, we get:

{\overset{α}{^}_{0} + \frac{n + 1}{2} \overset{α}{^}_{1} = \overset{ˉ}{X}, \frac{n + 1}{2} \overset{α}{^}_{0} + \frac{( n + 1 ) ( 2 n + 1 )}{6} \overset{α}{^}_{1} = \frac{1}{n} \sum_{t = 1}^{n} t X_{t} .

{\overset{α}{^}_{1} = \frac{12}{n ( n ^{2} - 1 )} \sum_{t = 1}^{n} t (X_{t} - \overset{ˉ}{X}), \overset{α}{^}_{0} = \overset{ˉ}{X} - \frac{n + 1}{2} \overset{α}{^}_{1} .

Note that this is the simple-regression formula, with $\overset{α}{^}_{1} = Cov (t, X) / Var (t)$ .

We now consider the matrix formula by stacking the data into $X = (X_{1}, \dots, X_{n})^{T}$ and the coefficients into $α = (α_{0}, \dots, α_{k})^{T}$ .

Let $A$ be the $n \times (k + 1)$ matrix whose row $t$ is $(1, t, t^{2}, \dots, t^{k})$ , called the Vandermonde matrix.

Then the model becomes $X = A α + Y$ , which gives:

\overset{α}{^} = (A^{T} A)^{- 1} A^{T} X .

(The normal equation in matrix form) is $\partial C / \partial α = - 2 A^{T} (X - A α) = 0$ , so $A^{T} A \overset{α}{^} = A^{T} X$ , given that $n \geq k + 1$ (i.e., we need at least as many data points as parameters, so that $A$ has full column rank and $A^{T} A$ is invertible).

In this case, if the noise has zero mean, $E [\overset{α}{^}] = α$ , i.e., the estimator is unbiased.

Furthermore, in the case of white noise, $Var (\overset{α}{^}) = σ^{2} (A^{T} A)^{- 1}$ .

Estimating trend and seasonality together (nonparametric approach)

We now consider the following model:

X_{t} = m_{t} + S_{t} + Y_{t}

where $m_{t}$ is the slowly moving trend, $S_{t}$ is a seasonal component with known period $d$ such that $S_{t + d} = S_{t}$ and sums to zero over one period, $\sum_{t = 1}^{d} S_{t} = 0$ (without it, you could move a constant between $m_{t}$ and $S_{t}$ and get the same $X_{t}$ ).

Again, we assume that $m_{t}$ is locally linear.

The starting idea is that, to estimate the trend, we have to remove the seasonal component; averaging over exactly one full period does exactly this.

So, given that the seasonal period has length $d$ , we can set $q$ to be the half-width of the interval, so the simple equal-weight averaging works as follows:

\overset{m}{^}_{t} = \frac{1}{d} j = - q \sum q X_{t - j} (d odd, q = (d - 1) /2)

\overset{m}{^}_{t} = \frac{1}{d} \frac{1}{2} X_{t - q} + \frac{1}{2} X_{t + q} + j = - (q - 1) \sum q - 1 X_{t - j} (d even, q = d /2)

We now focus on the case where $d$ is even (even if the opposite case can be solved similarly), so we can say:

\overset{m}{^}_{t} = =: A_{t} \frac{1}{d} \frac{1}{2} m_{t - q} + \frac{1}{2} m_{t + q} + j = - (q - 1) \sum q - 1 m_{t - j} + =: B_{t} \frac{1}{d} \frac{1}{2} S_{t - q} + \frac{1}{2} S_{t + q} + j = - (q - 1) \sum q - 1 S_{t - j} + =: C_{t} \frac{1}{d} \frac{1}{2} Y_{t - q} + \frac{1}{2} Y_{t + q} + j = - (q - 1) \sum q - 1 Y_{t - j}

(Note that this looks very similar to a moving average...)

Indeed, the weights $w_{- q} = w_{q} = \frac{1}{2 d}$ and $w_{j} = \frac{1}{d}$ for $j \in {- (q - 1), \dots, (q - 1)}$ satisfy

j = - q \sum q w_{j} = 1 and j = - q \sum q j w_{j} = 0,

so we can use the same argument as before to show that $A_{t} \approx a + b t \approx m_{t}$ .

Furthermore, the window covers exactly one period, and $S_{t - q} = S_{t + q}$ since they are $d$ steps apart, so the two half-weights add up to a single term and $B_{t}$ is the sum of $S$ over one full period, i.e., $B_{t} = 0$ .

Lastly, using the WLLN, we can show that $C_{t} \approx 0$ with high probability for $q$ large enough, so $\overset{m}{^}_{t} \approx m_{t}$ .

Remark that $q$ is large enough if and only if $d$ is large enough. Note also that we always take $d$ to be large enough because if $d_{0}$ is the smallest period of the seasonal component, then $d = k d_{0}$ is also a period for $k \in N$ .

Nonparametric estimation of the seasonal component

Suppose now that we know the trend $m_{t}$ exactly, so we can consider it as a “new observation” and consider:

$Z_{t} = X_{t} - m_{t} = s_{t} + Y_{t}$

Let $r \geq 1$ be the largest integer such that $1 + r d \leq n$ . Then we have that:

$s_{1} = s_{1 + d} = \dots = s_{1 + r d} = a$

where $a$ is unknown. So we can write:

$Z_{1 + (k - 1) d} = a + Y_{1 + (k - 1) d}$

To estimate $a$ , we then naturally use the empirical sample average, that is:

\overset{a}{^} = \frac{1}{r + 1} k = 1 \sum r + 1 Z_{1 + (k - 1) d} = \frac{1}{r + 1} k = 1 \sum r + 1 (X_{1 + (k - 1) d} - m_{1 + (k - 1) d}) .

Instead, suppose now the trend $m_{t}$ is unknown, but we can replace it with $\overset{m}{^}_{t}$ for $q + 1 \leq t \leq n - q$ (since the two-sided window must fit inside the data). Then $W_{j}$ is the average deviation from the trend at seasonal position $j$ , defined as:

W_{j} = \frac{1}{N _{j}} k : q + 1 \leq j + k d \leq n - q \sum (X_{j + k d} - \overset{m}{^}_{j + k d})

N_{j} = # {k : q + 1 \leq j + k d \leq n - q} = # {k : \frac{q + 1 - j}{d} \leq k \leq \frac{n - q - j}{d}}

Then, in order to enforce the constraint $\sum_{i = 1}^{d} s_{i} = 0$ , the final estimator is given by:

\overset{s}{^}_{j} = {W_{j} - \frac{1}{d} \sum_{i = 1}^{d} W_{i}, \overset{s}{^}_{j - d}, j \in {1, \dots, d}, j > d .

Intuition: the estimator averages across periods (to kill noise), subtracts the trend estimate (to isolate the season), and then re-centres the result (to fix the identifiability ambiguity between $m$ and $S$ ).

Parametric Estimation of the seasonal component

We now introduce the parametric estimation of $s$ , with the assumption that the seasonal component follows a harmonic function, that is:

s_{t} = a_{0} + j = 1 \sum k (a_{j} cos (\frac{2 π α _{j} t}{d}) + b_{j} sin (\frac{2 π β _{j} t}{d}))

where $k \in N$ , $α_{j}, β_{j} \in N$ are known, and $a_{0}, a_{j}$ , and $b_{j}$ are unknown. To avoid trivialities, we assume that $α_{j}, β_{j}$ are not multiples of $d$ for all $j$ .

Furthermore, note that

t = 1 \sum d cos (\frac{2 π α _{j} t}{d}) = t = 1 \sum d sin (\frac{2 π β _{j} t}{d}) = 0.

To prove that, simply set $c = e^{i 2 πa / d}$ , where $a \in N$ is not a multiple of $d$ , so we have $c \neq = 1$ . Then $c^{d} = 1$ , and we have that:

t = 1 \sum d c^{t} = c i = 0 \sum d - 1 c^{i} = c \frac{c ^{d} - 1}{c - 1} = 0.

Note that in order to enforce the constraint $\sum_{i = 1}^{d} s_{i} = 0$ , we must impose $a_{0} = 0$ .

To estimate $s_{t}$ , we can use the least-squares approach based on the de-trended data $\tilde{X}_{t} = X_{t} - \overset{m}{^}_{t}$ , with $q + 1 \leq t \leq n - q$ . Then, the LSE estimator of $(a_{1}, \dots, a_{k}, b_{1}, \dots, b_{k})$ is defined as:

(\overset{a}{^}_{1}, \dots, \overset{a}{^}_{k}, \hat{b}_{1}, \dots, \hat{b}_{k}) = \tilde{a}_{1}, \dots, \tilde{a}_{k}, \tilde{b}_{1}, \dots, \tilde{b}_{k} arg min t = q + 1 \sum n - q (\tilde{X}_{t} - j = 1 \sum k (\tilde{a}_{j} cos (\frac{2 π α _{j} t}{d}) + \tilde{b}_{j} sin (\frac{2 π β _{j} t}{d})))^{2} .

Estimating by differencing

Suppose now the model is $X_{t} = m_{t} + Y_{t}$ , with $E [Y_{t}] = 0$ . We want to eliminate $m_{t}$ (the trend component) instead of estimating it, so we apply an operator that sends a polynomial trend into a constant (or to zero). In particular, we then define:

$B$ , the backward shift operator, which moves the series back by one step, defined as:

$B X_{t} = X_{t - 1}$
$\nabla$ , the lag-1 differencing operator, which compares each value with the previous one (the discrete analogue of a derivative measuring the change per step), defined as:

$\nabla X_{t} = X_{t} - X_{t - 1} = (1 - B) X_{t}$

where $1 = B^{0}$ is the identity operator.

For $j \geq 2$ , we can define recursively:

$B^{j} = B \circ B^{j - 1} ⟹ B^{j} X_{t} = X_{t - j}$
$\nabla^{j} = \nabla \circ \nabla^{j - 1} ⟹ \nabla^{j} X_{t} = (1 - B)^{j} X_{t}$

We can also show that:

\nabla^{j} = (1 - B)^{j} = k = 0 \sum j (- 1)^{k} (k j) B^{k} ⟹ \nabla^{j} X_{t} = k = 0 \sum j (- 1)^{k} (k j) X_{t - k} .

Intuition: each application of $\nabla$ lowers the degree by one, so after $j$ applications, a degree- $j$ polynomial is left as a coefficient, and only the leading coefficient survives, multiplied by $j!$ (the same factor that appears when we differentiate $t^{j}$ $j$ times).

Note also that:

Stationarity is preserved, i.e., if $Y_{t}$ is stationary, then $\nabla^{j} Y_{t}$ is a linear combination of shifted copies of $Y_{t}$ and is therefore stationary too. $⟹$ The differenced series can therefore be managed with the usual stationary-friendly tools.
Zero mean, i.e., since $E [Y_{t}] = 0$ , we get $E [\nabla^{k} X_{t}] = k! a_{k}$ , so $a_{k}$ can be estimated from the sample mean of the differenced series.

We also have the following properties:

$\nabla X_{t} = X_{t} - X_{t - 1} = \frac{X _{t} - X _{t - 1}}{t - ( t - 1 )}$ .
If $m (t) = a_{0} + a_{1} t + \dots + a_{k} t^{k}$ with $t \in Z$ , then $\nabla^{k} m (t) = k! a_{k}$ .
If $m_{t} = a_{0} + a_{1} t + \dots + a_{k} t^{k}$ and $X_{t} = m_{t} + Y_{t}$ with $t \in Z$ , then $\nabla^{k} X_{t} = k! a_{k} + \nabla^{k} Y_{t}$ .
If ${Y_{t}}_{t \in Z}$ is stationary, then ${\nabla^{k} Y_{t}}_{t \in Z}$ is also stationary.