For every week: the concepts, the right R code, how to read the plots, and the conclusions on the patterns. Built on the course lectures (Krispin / TSstudio) and the official worksheet & homework solutions. Press Print → Save as PDF to keep it offline.
ts the axis "1,2,3" means lags 12, 24, 36 (×frequency).ccf(x, y), a peak at a negative lag means x leads y; never judge two series by lag-0 cor() alone.log).p > 0.05 and lowest validation MASE.A ts = values + start + frequency. Components: trend, seasonal, cycle, irregular. Additive vs multiplicative. Inspect, create, subset (window), change frequency (aggregate), and lag a series.
library(TSstudio); library(dplyr); library(lubridate) load("NGC_monthly.RData") # instructor's NGC series ts_info(NGC); frequency(NGC); deltat(NGC) # 12 ; 1/12 head(cycle(NGC)); head(time(NGC)); start(NGC); end(NGC) y <- ts(x, start = c(2020,1), frequency = 12) # build one yourself window(NGC, start=c(2010,1), end=c(2012,12)) # subset aggregate(NGC, nfrequency = 4, FUN = "mean") # monthly -> quarterly stats::lag(NGC, k = -12) # seasonal lag ts_plot(NGC, title="US Natural Gas", slide=TRUE)
ts_info() confirms the frequency is right.Moving-average smoothing (order = n_left + n_right + 1; one-sided vs two-sided; symmetric only for odd order). White noise & its ACF / Ljung–Box test. Classical decomposition into trend + seasonal + random.
library(TSstudio); library(forecast) ts_ma(USVSales, n_left=11, plot=FALSE) # one-sided order-12 (trailing) ts_ma(USVSales, n=6) # symmetric order-13 set.seed(1234); wn <- ts(rnorm(120)) acf(wn); Box.test(wn, lag=24, type="Ljung-Box") plot(decompose(USVSales)) # additive (default) plot(decompose(AirPassengers, type="multiplicative"))
decompose(), the random panel is the judge: flat noise = the model fits; a leftover wave / megaphone = wrong model.Single vs multiple seasonality; descriptive stats by frequency unit; density plots; ggseasonplot (normal + polar); ts_seasonal (normal/cycle/box); heatmap; quantile plots. Seasonal (fixed period) vs cycle (variable).
library(forecast); library(TSstudio); library(dplyr) df <- data.frame(month=factor(month.abb[cycle(USgas)], levels=month.abb), y=as.numeric(USgas)) df %>% group_by(month) %>% summarize(mean=mean(y), sd=sd(y)) # bar-chart data ggseasonplot(USgas, year.labels=TRUE); ggseasonplot(USgas, polar=TRUE) ts_seasonal(USgas, type="all") # normal + cycle + box ts_heatmap(USgas); ts_quantile(UKgrid)
Autocorrelation (ACF) and partial autocorrelation (PACF); lag plots; causation vs correlation; cross-correlation (CCF) for two series.
acf(nottem, lag.max=60); pacf(nottem, lag.max=60) ts_lags(nottem, lags=c(1,6,12,24,36)) # scatter vs each lag ccf(x=advertising, y=sales, lag.max=12) # lead-lag of two series cor(sales, advertising) # lag-0 only (can mislead)
ccf(x,y) peak at negative lag = x leads y (advertising leads sales by 3 weeks at lag −3).Workflow: prepare → train → test → evaluate → refit. Train/test split (train ≥70%, test = horizon); backtesting (expanding vs sliding); residual analysis; scores MSE/RMSE/MAE/MAPE/MASE; benchmarks; finalize on full data; prediction intervals.
library(TSstudio); library(forecast) sp <- ts_split(USgas, sample.out=12); train <- sp$train; test <- sp$test md <- auto.arima(train); fc <- forecast(md, h=12) accuracy(fc, test) # score vs actuals + benchmark checkresiduals(md) fc_final <- forecast(auto.arima(USgas), h=12) # refit on ALL data
checkresiduals() → residuals ~0, flat ACF, bell histogram, Ljung–Box p > 0.05 = white noise. Accuracy table → compare test RMSE/MAPE/MASE to the naïve; big test-vs-train gap = overfit; MASE < 1 beats naïve. Backtest error stable across partitions = robust.OLS review and the time-series problem (autocorrelated errors); features trend = 1:n and season = factor(month); trend-only, season-only, trend+season, polynomial trend.
train$trend <- 1:nrow(train) train$seasonal <- factor(month(train$ds, label=TRUE), ordered=FALSE) md1 <- lm(y ~ trend + seasonal, data=train) # workhorse md2 <- lm(y ~ trend + I(trend^2) + seasonal, data=train) # curved trend md <- tslm(train.ts ~ trend + season) # on the ts directly
I(trend^2) only if the trend curves (don't extrapolate it far); autocorrelated residuals → regression with ARIMA errors (Week 13).tslm; modeling outliers, structural breaks (H_t = max(0, t−k) = slope change), recurring non-seasonal events; engineering wday, month, lag365 for several seasonal cycles; anova selection; forecasting with future regressors; Box–Cox.
k <- which(time(USgas) >= 2010)[1]; trend <- 1:length(USgas) s_break <- pmax(trend - k, 0) # slope-change feature md <- tslm(USgas ~ season + trend + s_break) md3 <- tslm(train_ts ~ season + trend + wday + month + lag365, data=train_df) anova(md3) # which terms are significant fc <- forecast(final_md, h=h, newdata=UK_fc_df) # needs FUTURE predictors
wday/month profile is a separate seasonal cycle. anova p-value small = the term carries signal. Even high R²/low MAPE can hide autocorrelated residuals — check the ACF.Simple MA forecast (average the last m; m=1 = naïve); window-length trade-off; limitations (no trend/season, flattens); weighted MA (weights ≥ 0, sum to 1); seasonal weighting.
# all weight on lag 12 == seasonal naive fc_a <- sma_forecast(USgas_df, h=24, m=12, w=c(1, rep(0,11))) # blend seasonal lag (0.8) + recent (0.2) fc_b <- sma_forecast(USgas_df, h=24, m=12, w=c(0.8, rep(0,10), 0.2))
m reacts fast but is noisy; long m is smooth but flat; the forecast converges to a constant → it can't hold a trend and averages seasonality away. Weighted MA with all weight on lag 12 reproduces the seasonal oscillation.SES (level only, flat forecast; weight decays by (1−α); α=1 = naïve); Holt (+trend); Holt–Winters (+season, additive vs multiplicative); ets() auto-selects.
fc_ses <- ses(train, h=12) # flat fc_holt <- holt(train, h=8) # + trend md_hw <- HoltWinters(train); fc_hw <- forecast(md_hw, h=12) fit <- ets(train); fc <- forecast(fit, h=12) # auto E,T,S
Weak stationarity (constant mean/variance/autocovariance); differencing (1st, 2nd, seasonal); log for variance; random walk (non-stationary; its first difference is white noise); AR(1)/AR(p); AR(1) property ρₘ = φᵐ.
diff(y); diff(log(AirPassengers)) # difference / log-then-difference library(tseries); adf.test(y); kpss.test(y) # stationarity tests ndiffs(y); nsdiffs(y) # how many differences ar2 <- arima.sim(list(order=c(2,0,0), ar=c(0.9,-0.3)), n=500)
MA(1)/MA(q) (a weighted sum of recent shocks); identify MA; MA(1) has only lag-1 ACF nonzero; ARMA(p,q); order selection by AIC/BIC; residual check; forecasting.
ma2 <- arima.sim(list(order=c(0,0,2), ma=c(0.5,-0.3)), n=500) md_ma <- arima(ma2, order=c(0,0,2), method="ML") arma_md<- arima(arma, order=c(1,0,2)); AIC(md_a, md_b); BIC(md_a) checkresiduals(arma_md); forecast(ar(ar2), h=100)
checkresiduals Ljung–Box p > 0.05 = adequate. Forecast reverts to the mean; interval widens then stabilizes.ARIMA(p,d,q) (I = differencing); special cases (000 = white noise, 010 = random walk); choosing d; SARIMA(p,d,q)(P,D,Q)ₛ; seasonal AR/MA/differencing; tuning; grid-search AIC; USgas case study.
arima(robusta_price, order=c(1,1,0)); lmtest::coeftest(md, df=Inf) d12 <- diff(train, 12); d12_1 <- diff(d12, 1) # seasonal then ordinary best <- arima(train, order=c(1,1,1), seasonal=list(order=c(0,1,2))) fc <- forecast(best, h=12); accuracy(fc, test); checkresiduals(best)
coeftest → drop a clearly non-significant term. Ljung–Box p > 0.05 = adequate.OLS on time series leaves autocorrelated errors → understated standard errors; fix = regression with ARIMA errors (yₜ = β′xₜ + nₜ, nₜ ~ ARIMA); forecasting needs future predictor values.
X <- cbind(trend=seq_along(train), fourier(train, K=3)) Xf <- cbind(trend=length(train)+1:h, fourier(train, K=3, h=h)) fit <- auto.arima(train, xreg=X, lambda=0) # regression + ARIMA errors fc <- forecast(fit, xreg=Xf, h=h) # MUST pass future xreg checkresiduals(fit)
tslm residual ACF with spikes = autocorrelated errors (t-values untrustworthy). After the ARIMA-errors fit, the residual ACF should be flat → the ARIMA part absorbed the dynamics and the coefficients now have correct standard errors.No new plots — chain the reads: time plot → ACF/PACF → candidate fit → checkresiduals → forecast, each feeding the next.
y <- ... # your ts sp <- ts_split(y, sample.out=h); train <- sp$train; test <- sp$test models <- list(snaive=snaive(train,h), ets=forecast(ets(train),h), arima=forecast(auto.arima(train, lambda=0),h)) sapply(models, function(m) accuracy(m, test)["Test set","MASE"]) final <- auto.arima(y, lambda=0); forecast(final, h=24) # refit on all data
| You see | Conclusion → move |
|---|---|
| ACF decays slowly, stays outside band | Non-stationary → diff() |
| Tall ACF spikes at 12, 24… | Seasonal → diff(y,12) (D=1) |
| PACF cuts off at p · ACF tails off | AR(p) → c(p,d,0) |
| ACF cuts off at q · PACF tails off | MA(q) → c(0,d,q) |
| Both tail off | ARMA → compare AIC |
| Residual ACF flat + Ljung–Box p > 0.05 | Adequate → forecast |
| Decompose random panel flat / wavy | Model fits / switch add↔mult |
STAT 411 · Study Guide · built from the course lectures (Krispin / TSstudio) and the official worksheet & homework solution keys.