Time Series and Forecasting: Quick recap

01 Time series and decomposition

A time series is a sequence $y_1, \dots, y_n$ in time order; its values are dependent, so never shuffle them or split them at random.
Look at a time plot, a seasonal plot (one line per period) and a lag plot ( $y_t$ against $y_{t-k}$ ) before any modelling.
Additive decomposition: $y_t = T_t + S_t + R_t$ ; multiplicative: $y_t = T_t \cdot S_t \cdot R_t$ , which becomes additive after a log transform.
A seasonal pattern has a fixed calendar period $m$ ; a cycle has no fixed period and belongs to the trend-cycle.
A trailing $k$ -MA lags a linear trend by $(k-1)/2$ steps; a symmetric MA has no lag; any MA over exactly one period removes a pattern with that period.
For even $m$ the centred trend is the $2 \times m$ MA, with weights $\tfrac{1}{2m}, \tfrac1m, \dots, \tfrac1m, \tfrac{1}{2m}$ .
Classical decomposition: trend by MA, seasonal indices as centred averages of the detrended values, remainder as the rest; it has no trend at the ends and a fixed seasonal pattern. STL is the robust default.
Strength of seasonality $F_S = \max(0, 1 - \operatorname{Var}(R_t)/\operatorname{Var}(S_t + R_t))$ , and $F_T$ in the same way with $T_t$ .

02 Stationarity and autocorrelation

White noise has mean 0, variance $\sigma^2$ and no autocorrelation; a random walk $y_t = y_{t-1} + \varepsilon_t$ has $\operatorname{Var}(y_t) = t\sigma^2$ and is not stationary.
Weak stationarity: constant mean, finite variance, and $\operatorname{Cov}(y_t, y_{t+k}) = \gamma(k)$ depends only on the lag $k$ . Trends and seasonal patterns break it.
A stationary AR(1) process $y_t = \phi y_{t-1} + \varepsilon_t$ ( $|\phi| < 1$ ) has variance $\sigma^2/(1-\phi^2)$ and ACF $\rho(k) = \phi^k$ ; $\phi = 1$ is a unit root (a random walk).
The sample ACF divides by $n$ , which keeps $|\hat\rho(k)| \le 1$ and the autocovariance matrix positive semidefinite; for white noise about 95 % of the values lie in $\pm 1.96/\sqrt n$ , one lag at a time.
The Ljung–Box statistic $Q = n(n+2)\sum_{k=1}^{\ell} \hat\rho(k)^2/(n-k)$ is approximately $\chi^2_\ell$ for white noise ( $\ell - p - q$ degrees of freedom for ARMA residuals).
The PACF $\phi_{kk}$ is the last coefficient of the best AR( $k$ ) predictor; $\phi_{22} = (\rho(2) - \rho(1)^2)/(1 - \rho(1)^2)$ ; for AR( $p$ ) it is 0 after lag $p$ .
Differencing $\nabla = 1 - B$ and seasonal differencing $\nabla_m = 1 - B^m$ (after a log if needed) make many series stationary; $\hat\rho(1) \approx -0.5$ signals over-differencing.
ADF has the null "unit root", KPSS the null "stationary"; a large p-value proves nothing. Regressing one random walk on an independent one gives spurious significance.

03 Baselines and backtesting

Benchmarks: mean $\bar y$ ; naive $y_n$ ; seasonal naive $y_{n+h-m(k+1)}$ with $k = \lfloor (h-1)/m\rfloor$ ; drift $y_n + h\,(y_n - y_1)/(n-1)$ . Each is the best forecast for one simple model (white noise, random walk, seasonal random walk, random walk with drift).
Under squared error the best forecast is the conditional mean; for a random walk this is the last value, and the $h$ -step error variance is $h\sigma^2$ (the fan).
Good residuals are uncorrelated and have mean zero (essential); constant variance and normality are needed only for the usual intervals. Residuals are in-sample and do not measure accuracy.
Evaluate with a rolling origin: many forecast origins, an expanding or sliding training window, a gap when data arrive late or features overlap the future, and errors reported per horizon $h$ .
Random K-fold splits of a time series interpolate and see the future, so they are usually too optimistic; they are valid only for autoregressive models with uncorrelated errors on stationary data.
MSE is minimised by the mean, MAE by the median, and MAPE by a value below the median ( $e^{\mu - \sigma^2}$ for a log-normal); choose the metric from the decision.
MASE divides the test MAE by the in-sample MAE of the one-step seasonal naive forecast, $\frac{1}{n-m}\sum_{t=m+1}^{n}\lvert y_t - y_{t-m}\rvert$ ; it is scale-free, and MASE < 1 means better than that in-sample benchmark.
Avoid leakage: fit every preprocessing step inside each training window, never fill from the future, use only data and covariates available at each origin.

04 Exponential smoothing

SES forecast: $\hat y_{t+h\mid t} = \ell_t$ (flat), with $\ell_t = \alpha y_t + (1-\alpha)\ell_{t-1} = \ell_{t-1} + \alpha e_t$ ; equivalently a weighted average with weights $\alpha(1-\alpha)^j$ plus $(1-\alpha)^t\ell_0$ , summing to 1.
Small $\alpha$ means long memory (average age $(1-\alpha)/\alpha$ ); $\alpha = 1$ is the naive forecast. Fit $\alpha$ (and $\ell_0$ ) by minimising the SSE of one-step errors.
Holt: $\hat y_{t+h\mid t} = \ell_t + h b_t$ , a straight line; the damped trend uses $\ell_t + (\phi + \dots + \phi^h) b_t$ , which tends to $\ell_t + \frac{\phi}{1-\phi} b_t$ .
Holt-Winters adds $m$ seasonal states; the forecast uses $s_{t+h-m(k+1)}$ with $k = \lfloor (h-1)/m \rfloor$ . Additive for swings of fixed size, multiplicative (or logs) for swings that grow with the level.
Error-correction form: each state = its prediction + gain × one-step error, with gains $\alpha$ , $\alpha\beta^*$ , $\gamma$ .
ETS (Error, Trend, Seasonal) adds an error model to the recursions: ETS(A,N,N) is $y_t = \ell_{t-1} + \varepsilon_t$ , $\ell_t = \ell_{t-1} + \alpha\varepsilon_t$ ; its point forecasts equal SES.
The ETS model gives a likelihood (model choice by AICc) and prediction intervals; for ETS(A,N,N) the $h$ -step variance is $\sigma^2[1 + \alpha^2(h-1)]$ .
On the CO₂ backtest, additive and multiplicative Holt-Winters have MASE about 0.27, against 1.00 for the seasonal naive forecast; non-seasonal methods lose to the benchmark.

05 AR, MA and ARMA

AR( $p$ ): $\phi(B)y_t = c + \varepsilon_t$ with $\phi(z) = 1 - \phi_1 z - \dots - \phi_p z^p$ ; a stationary AR( $p$ ) has mean $\mu = c/(1 - \phi_1 - \dots - \phi_p)$ .
An AR( $p$ ) is stationary (and causal) when all roots of $\phi(z)$ lie outside the unit circle, equivalently all eigenvalues of the companion matrix lie inside it; eigenvalues are the reciprocals of the roots.
AR(2) is stationary inside the triangle $\phi_1 + \phi_2 < 1$ , $\phi_2 - \phi_1 < 1$ , $|\phi_2| < 1$ ; below the parabola $\phi_1^2 + 4\phi_2 < 0$ it has a quasi-cycle with period $2\pi/\omega$ , $\cos\omega = \phi_1/(2\sqrt{-\phi_2})$ .
Yule–Walker: $\rho(k) = \sum_{i=1}^p\phi_i\rho(k-i)$ for $k \ge 1$ , a Toeplitz system that gives the ACF from $\phi$ or estimates $\phi$ from the sample ACF.
MA( $q$ ) is always stationary and its ACF cuts off after lag $q$ ; for MA(1), $\rho(1) = \theta/(1 + \theta^2)$ , so $|\rho(1)| \le 0.5$ .
$\theta$ and $1/\theta$ give the same ACF; choose the invertible MA (roots of $\theta(z)$ outside the unit circle), which is an AR( $\infty$ ), so its PACF decays.
Identification: AR( $p$ ) has a PACF cut-off at $p$ , MA( $q$ ) an ACF cut-off at $q$ , ARMA neither; then compare candidates by AIC/AICc/BIC on the same data and check residuals with Ljung–Box.
AR(2) fitted to the yearly sunspots gives $\hat\phi = (1.392, -0.690)$ and a quasi-period of 10.9 years.

06 ARIMA and SARIMA

ARIMA( $p,d,q$ ): $\phi(B)(1 - B)^d y_t = c + \theta(B)\varepsilon_t$ ; the $d$ -times differenced series is a stationary, invertible ARMA.
Random walk = ARIMA(0,1,0); with a constant it has drift; SES with parameter $\alpha$ = ARIMA(0,1,1) with $\theta = \alpha - 1$ .
Forecast recursion: write the model at time $n + h$ , replace future values by forecasts, future shocks by 0 and past shocks by residuals.
Forecast error variance: $\operatorname{Var}(e_{n+h}) = \sigma^2\sum_{j=0}^{h-1}\psi_j^2$ ; bounded for a stationary model, grows like $h$ (width like $\sqrt h$ ) for a random walk.
The constant and $d$ fix the long-run forecast: mean ( $d = 0$ ), constant or line ( $d = 1$ ), line or quadratic ( $d = 2$ ).
SARIMA $(p,d,q)(P,D,Q)_m$ multiplies seasonal polynomials in $B^m$ ; the airline model $(0,1,1)(0,1,1)_{12}$ has MA terms at lags 1, 12, 13 and only two MA coefficients.
Box–Jenkins: transform, difference, identify, compare by AICc on the same differencing, check residuals (Ljung–Box), forecast, back-transform ( $\exp$ of a log forecast is the median).
Regression with ARIMA errors needs the future covariates; Fourier terms model long seasonal periods.

07 State space and Kalman filter

A linear Gaussian state-space model is $x_t = Fx_{t-1} + w_t$ , $w_t \sim \mathcal N(0, Q)$ (hidden state) and $y_t = Hx_t + v_t$ , $v_t \sim \mathcal N(0, R)$ (observation). Local level, trend, seasonal and ARMA models are all special cases.
The Kalman filter predicts $\hat x_{t\mid t-1} = F\hat x_{t-1\mid t-1}$ , $P_{t\mid t-1} = FP_{t-1\mid t-1}F^\top + Q$ , then updates with the innovation $\nu_t = y_t - H\hat x_{t\mid t-1}$ , its covariance $S_t = HP_{t\mid t-1}H^\top + R$ and the gain $K_t = P_{t\mid t-1}H^\top S_t^{-1}$ .
One scalar update is Bayes' rule for two Gaussians: precisions add, and the new mean is the prediction plus $K = P/(P + R)$ times the innovation. $K \to 0$ when $R \to \infty$ and $K \to 1$ when $R \to 0$ .
Under the model assumptions the filter gives the exact posterior $x_t \mid y_{1:t} \sim \mathcal N(\hat x_{t\mid t}, P_{t\mid t})$ (Gaussian conditioning, by induction).
For the local level model the gain converges to $K_\infty = \tfrac12(\sqrt{q^2 + 4q} - q)$ with $q = Q/R$ , and the filter becomes simple exponential smoothing with $\alpha = K_\infty$ .
Forecasts repeat the predict step; missing observations skip the update, so the variance grows (by $Q$ per step for the local level).
The log-likelihood is $-\tfrac12\sum_t(\log\det(2\pi S_t) + \nu_t^\top S_t^{-1}\nu_t)$ (prediction-error decomposition); maximise it to estimate $Q$ , $R$ and other parameters. Nile: $\hat R \approx 15\,099$ , $\hat Q \approx 1\,469$ .
The RTS smoother uses all data and has smaller variances; use it to describe the past, never at a forecast origin.

08 ML for forecasting

Reduction to supervised learning: each row has features known at the origin (lags $y_{s-h}, \dots, y_{s-h-L+1}$ , rolling statistics, calendar, covariates) and the target $y_s$ ; the table has $n - L - h + 1$ rows.
Leakage is any feature that uses information after the forecast origin; compute window features after shift(h), and fit encodings and scales on training data only.
Covariates are static, known in advance, or observed only; in the direct strategy a lag $k$ of the target or of an observed-only covariate is usable at horizon $h$ only if $k \ge h$ .
Tree ensembles predict means of training targets and are constant beyond the largest split thresholds, so they cannot follow a trend; model differences, ratios or residuals from a trend instead.
Recursive forecasting uses one model and feeds forecasts back (errors accumulate); direct forecasting uses one model per horizon; multi-output predicts all horizons at once.
With one lag, $\mathrm{MSE}_{\text{rec}}(h) - \mathrm{MSE}_{\text{dir}}(h) = \gamma(0)(\rho(1)^h - \rho(h))^2 \ge 0$ : direct is never worse under misspecification, and equal for a true AR(1).
A global model trained on many series shares patterns and has more rows than a local model; it needs per-series scaling computed from training data.
Validate gradient boosting with a rolling origin and a time-based early-stopping block, never with a random split.

09 Deep learning for forecasting

A neural forecaster maps a context of length $L$ to $H$ future values or to a distribution; each window is scaled with its own context statistics (mean scaling, RevIN), and the scaling is undone at the output. Scaling with future values is leakage.
DeepAR is a global autoregressive RNN with a likelihood head: training with true inputs (teacher forcing) is maximum likelihood, and ancestral sampling gives exact samples of the joint forecast distribution.
A stack of causal convolutions with kernel $k$ and dilations $d_i$ has receptive field $r = 1 + (k-1)\sum_i d_i$ and never uses future inputs; doubling dilations make $r$ grow exponentially with depth. A TCN block has two convolutions.
DLinear, $\hat y = W_s(x - Ax) + W_tAx + b$ , has exactly the function class of a linear map on the window, and such linear models beat several long-horizon transformers on standard benchmarks: always include a linear baseline.
PatchTST turns patches of $P$ values with stride $S$ into $N = \lfloor (L-P)/S \rfloor + 2$ tokens, cutting attention cost by about $S^2$ , and shares one model across channels (channel independence).
Time-series foundation models are pretrained on many series and used zero-shot (no training on the target data) or fine-tuned; they differ in tokens (quantised values, patches), heads (bins, points, distributions, quantiles, flow matching) and backbone.
Chronos divides by the mean absolute value $s$ and quantises into about 4,000 bins on $[-15, 15]$ ; the rounding error is at most $15s/(B-1)$ and values beyond $\pm 15s$ are clipped.
Evaluate deep and pretrained models like any forecaster: rolling-origin backtests on your own series, against seasonal naive, ETS and boosting, and watch for benchmark contamination.

10 Probabilistic forecasting

A probabilistic forecast is a predictive distribution, a set of quantiles, an interval or a set of sample paths; only sample paths answer questions about sums or maxima over several horizons.
The pinball loss $\rho_\tau(u) = \max\{\tau u, (\tau - 1)u\}$ with $u = y - q$ has expected value minimised at the $\tau$ -quantile, because its derivative is $F(q) - \tau$ ; training any model with it gives quantile regression.
The best stock level under linear costs is the quantile at the critical ratio $\tau = c_u/(c_u + c_o)$ .
Calibration means that the PIT values $u_t = F_t(y_t)$ are uniform (central intervals have their nominal coverage); a U-shaped PIT histogram means intervals that are too narrow, a hump means too wide.
Calibration is necessary but not sufficient: maximise sharpness subject to calibration.
The CRPS $\int (F(x) - \mathbf 1\{x \ge y\})^2\,dx = \mathbb E|X - y| - \frac12\mathbb E|X - X'| = 2\int_0^1 \rho_\tau(y - F^{-1}(\tau))\,d\tau$ is a strictly proper score in the units of $y$ ; for a point forecast it is the absolute error.
Split conformal uses the $\lceil (n + 1)(1 - \alpha)\rceil$ -th smallest calibration score as the half-width; with exchangeable data its coverage is at least $1 - \alpha$ (marginal, not conditional).
Time series break exchangeability; adaptive conformal inference $\alpha_{t+1} = \alpha_t + \gamma(\alpha - \mathrm{err}_t)$ keeps the long-run miss rate within $(\max(\alpha_1, 1 - \alpha_1) + \gamma)/(\gamma T)$ of $\alpha$ for any sequence.

11 Anomaly detection

Time-series anomalies come in four kinds: point anomalies, contextual anomalies (normal value, wrong time), collective anomalies (unusual shape) and change points (a lasting change of the distribution).
The standard recipe is expected value (decomposition, one-step forecast or filter), residual, robust score and threshold; the expected value must include the season, or contextual anomalies stay hidden.
The robust z-score $(r_t - \operatorname{med} r)/(1.4826\cdot\mathrm{MAD})$ is on the scale of an ordinary z-score for Gaussian data and is not pulled by the anomalies; the mean and standard deviation suffer from masking.
Choose the threshold from a false-alarm budget: $N$ checks with per-check rate $p$ give $Np$ false alarms, for example 14.4 per day for one check per minute at $p = 0.01$ .
The CUSUM $C_t = \max(0, C_{t-1} + z_t - k)$ with alarm at $C_t > h_c$ adds up evidence for a shift; $k = \delta/2$ makes it the maximised log-likelihood ratio for a shift of $\delta$ , and it detects small shifts much faster than a single-value rule at the same false-alarm rate.
Detectors of changes are compared by the average run length: large $\mathrm{ARL}_0$ (no change) and small $\mathrm{ARL}_1$ (delay after a shift).
The matrix profile gives each window the z-normalised distance to its nearest non-overlapping neighbour; its maximum, the discord, is a collective anomaly.
Point-adjusted F1 lets random detectors score high; evaluate with event recall, alarm precision, delay and false alarms per day, against a random and a simple robust baseline.