Time Series Analysis in R —

Comprehensive Learning Guide

For: MBA Students & Beginners in Time Series Analysis Tool: R Programming Language Coverage: 10 Lab Sessions (Labs 1–10,
with Lab 9 having two parts)

Table of Contents
1. Lab 1 — Introduction to Time Series & Decomposition
2. Lab 2 — Time Series Visualization Techniques
3. Lab 3 — Forecasting Benchmarks, Transformations & Residuals
4. Lab 4 — Forecast Accuracy & Time Series Regression
5. Lab 5 — Moving Averages & Decomposition Methods
6. Lab 6 — Exponential Smoothing Methods & ETS Models
7. Lab 7 — Stationarity, Differencing & ACF/PACF
8. Lab 8 — KPSS Test & ARIMA Models
9. Lab 9 — Seasonal ARIMA (SARIMA) Models
10. Lab 9 Part 2 — Stock Data & Financial Returns
11. Lab 10 — GARCH Models for Volatility Modelling
12. Glossary

Lab 1 — Introduction to Time Series &

Decomposition
Overview
This lab introduces the most fundamental building block of time series analysis: what is a time series, how to load one into R, how to
convert raw data into a time series object, and how to decompose it into its structural components.

Full Code
kings <- scan("[Link] skip=3)
kings
kingstimeseries <- ts(kings)
kingstimeseries

births <- scan("[Link]

births
birthstimeseries <- ts(births, frequency=12, start=c(1946,1))
birthstimeseries

[Link](birthstimeseries)

birthstimeseriescomponents <- decompose(birthstimeseries, type = 'multiplicative')

plot(birthstimeseriescomponents)
birthstimeseriescomponents$seasonal
birthstimeseriescomponents$trend

birthstimeseriescomponentsadditive <- decompose(birthstimeseries, type = 'additive')

plot(birthstimeseriescomponentsadditive)
birthstimeseriescomponentsadditive$seasonal
birthstimeseriescomponentsadditive$trend

1. Code Walkthrough — Line by Line

Loading the Kings Dataset

kings <- scan("[Link] skip=3)

scan() — A base R function that reads data from a file or URL, returning a numeric vector.
"[Link] — The URL of the dataset. The dataset contains the age at death of 42 successive kings of
England.
skip=3 — Skips the first 3 lines of the file (these are header/description lines, not data). Without this, R would try to read text as
numbers and throw an error.
Result: kings is now a plain numeric vector: [38, 43, 45, 36, ...] — just numbers, with no time information attached.

kingstimeseries <- ts(kings)

ts() — The workhorse function in R for creating a Time Series (ts) object. It takes a plain vector and wraps it with time metadata.
Here, no frequency or start is provided, so R assumes:
The data starts at time = 1
The frequency = 1 (one observation per "period" — here, one king per period)
Result: kingstimeseries is now a ts object with an implied time index (1, 2, 3, ..., 42).

Why convert to ts? R's time series functions ([Link], decompose, forecast, acf, etc.) all require a ts object. A plain vector
doesn't have any temporal structure.
 Loading the NY Births Dataset

births <- scan("[Link]

Loads monthly birth counts in New York City. No skip needed here — the file starts with data directly.
Result: births is a plain numeric vector of monthly values.

birthstimeseries <- ts(births, frequency=12, start=c(1946,1))

frequency=12 — Tells R there are 12 observations per year (monthly data). This is crucial for decomposition and seasonal
analysis.
start=c(1946,1) — Tells R the series starts in January (month 1) of 1946. The c(year, period) format is used.
Result: birthstimeseries is now a properly labelled monthly ts object, from Jan 1946 onwards.

Frequency Value Data Type

1 Annual
4 Quarterly
12 Monthly
52 Weekly
365 Daily

Plotting the Time Series

[Link](birthstimeseries)

[Link]() — A base R function specifically designed for plotting ts objects. It automatically uses the time index on the X-axis.
This gives you the first visual look at your data: trends, seasonality, and anomalies.

What the Graph Looks Like

X-axis: Time (years, from 1946 to approximately 1959)

Y-axis: Number of births per month
Shape: A wavy line that oscillates up and down each year (seasonal peaks in summer, dips in winter) while also showing a gentle
upward drift over time.
Pattern to note: The amplitude (height) of the waves might appear to grow over time — this is a signal that multiplicative
decomposition might be more appropriate.

Decomposition — Multiplicative

birthstimeseriescomponents <- decompose(birthstimeseries, type = 'multiplicative')

plot(birthstimeseriescomponents)

decompose() breaks a time series into three structural parts:

Component Symbol Meaning

Trend (T) Tₜ The long-run direction (up, down, flat)
Seasonal (S) Sₜ Regular, repeating pattern within each year
Random / Residual (R) Rₜ What's left over after removing T and S
 Multiplicative Model Formula:

Yₜ = Tₜ × Sₜ × Rₜ

Yₜ = Observed value at time t

Tₜ = Trend component at time t
Sₜ = Seasonal component at time t (expressed as a ratio, e.g., 1.15 means 15% above average)
Rₜ = Irregular/Random component (also a ratio, ideally close to 1.0)

When to use Multiplicative? Use multiplicative decomposition when the seasonal fluctuations grow proportionally with the level
of the series. For example, if sales double in December when the annual level is high, but only grow a little when annual sales are low
— that's multiplicative behaviour.

What the Decomposition Plot Looks Like

The plot() of the decomposed object produces 4 panels stacked vertically:

1. Top panel — Observed (Yₜ): The original data as-is.

2. Trend (Tₜ): A smooth, slow-moving curve showing the general direction. The jagged seasonality is removed. For births, this might
show a slow upward curve.
3. Seasonal (Sₜ): A repeating identical pattern for every year — like a template wave. Values above 1.0 = above-average months;
below 1.0 = below-average months. For births, summer months (July, August) likely show values > 1.0.
4. Random (Rₜ): Should look like "white noise" — no pattern, scattered around 1.0. If you see patterns here, your model is not
capturing all structure.

birthstimeseriescomponents$seasonal
birthstimeseriescomponents$trend

$seasonal — Extracts the seasonal component values (one per month over the full dataset)
$trend — Extracts the trend component values (Note: the first and last few values will be NA because the trend is computed using a
centred moving average which needs data on both sides)

Decomposition — Additive

birthstimeseriescomponentsadditive <- decompose(birthstimeseries, type = 'additive')

plot(birthstimeseriescomponentsadditive)

Additive Model Formula:

Yₜ = Tₜ + Sₜ + Rₜ

Here, Sₜ is measured in original units (e.g., +300 births in July, −200 births in February)
The seasonal component is the same fixed amount regardless of the trend level
 When to use Additive? Use additive decomposition when the seasonal fluctuations stay roughly constant in size regardless of
the overall level of the series.

Comparing Additive vs Multiplicative

Feature Additive Multiplicative

Formula Y=T+S+R Y=T×S×R
Seasonal size Constant over time Grows with trend
Seasonal values In original units (e.g., ±300) As ratios (e.g., 1.15 or 0.85)
Use when Stable seasonality Proportional seasonality
Example Temperature Retail sales, tourism

2. Statistical Concepts
What is a Time Series?

A time series is a sequence of data points collected at successive, equally spaced time intervals.

Examples: Monthly GDP, daily stock prices, annual rainfall, quarterly sales
Key property: Each observation is NOT independent — today's value depends on yesterday's value. This violates the assumption of
ordinary regression.

Decomposition

Decomposition is the process of separating a time series into its component parts to better understand and forecast it.

Real-world analogy: Think of your monthly electricity bill. It has:

A trend (slowly rising as you buy more appliances)

A seasonal pattern (spikes in summer from air conditioning, spikes in winter from heating)
A random component (the month you accidentally left a light on all day)

Decomposition helps you isolate each of these.

3. Inference & Conclusions

The kings dataset has no seasonality (annual data, one observation per king) — it only has a trend and random component.
The births dataset has clear seasonality (monthly data) — births peak in summer and dip in winter.
The multiplicative model is preferred when the seasonal amplitude increases over time.
After decomposition, examining the random (residual) component tells you if the model has captured all structure — ideally it
should look like white noise.
The trend component of the births data shows the long-run direction — useful for long-term planning.

Lab 2 — Time Series Visualization Techniques

Overview
This lab introduces a rich set of visualization tools from the fpp2 and GGally packages. Good visualization is the first step in any time
series analysis — it reveals patterns, relationships, and anomalies that guide all subsequent modelling decisions.

Full Code

[Link]("fpp2")
[Link]("GGally")

library(fpp2)
library(GGally)

help("melsyd")
melsyd

# Time plot
autoplot(melsyd[,"[Link]"]) + ggtitle("Economy class passengers: Melbourne-Sydney") +
xlab("Year") + ylab("Thousands")

help("a10")
a10
autoplot(a10) + ggtitle("Antidiabetic drug sales") + ylab("$million") + xlab("year")

# Seasonal plot
ggseasonplot(a10, [Link] = TRUE, [Link] = TRUE) + ylab("$million") +
ggtitle("Seasonal plot: antidiabetic drug sales")

ggseasonplot(a10, polar = TRUE) + ylab("$million") +

ggtitle("Seasonal plot: antidiabetic drug sales")

# Seasonal subseries plot

ggsubseriesplot(a10) + ylab("$million") +
ggtitle("Seasonal sub series plot: antidiabetic drug sales")

# Scatter plot
help("elecdemand")
elecdemand

autoplot(elecdemand[,c("Demand","Temperature")], facets=TRUE) +
xlab("Year:2014") + ylab("") +
ggtitle("Half-hourly electricity demand: Victoria, Australia")

qplot(Temperature, Demand, data = [Link](elecdemand)) +

ylab("Demand (GW)") + xlab("Temperature (Celsius)")

help("visnights")
visnights

autoplot(visnights[,1:5], facets=TRUE)
GGally::ggpairs([Link](visnights[,1:5]))
 1. Code Walkthrough
Package Setup

[Link]("fpp2")
library(fpp2)

fpp2 — "Forecasting: Principles and Practice, 2nd edition" — a package by Rob Hyndman that contains datasets and functions for
modern forecasting in R. It automatically loads ggplot2 and forecast packages.
GGally — An extension of ggplot2 that provides advanced plot types like pair plots (ggpairs).

The melsyd Dataset

melsyd
autoplot(melsyd[,"[Link]"]) + ggtitle(...) + xlab("Year") + ylab("Thousands")

melsyd — Melbourne–Sydney airline passengers dataset (weekly, from 1987–1992). Contains columns for Economy, Business, and
First Class passengers.
melsyd[,"[Link]"] — Extracts only the "Economy Class" column using matrix-style indexing.
autoplot() — The fpp2 version of plotting — it automatically detects the ts object and produces a ggplot2-style time plot.
ggtitle() — Adds a main title to the plot.
xlab() / ylab() — Labels for X and Y axes.

What the Graph Looks Like

X-axis: Years (1987–1992)

Y-axis: Passenger numbers in thousands
Pattern: The line shows weekly fluctuations, with a notable gap (missing data / interruption in service — possibly a strike or
grounding). A good analyst would immediately note and investigate this gap.

The a10 Dataset (Antidiabetic Drug Sales)

a10
autoplot(a10) + ggtitle("Antidiabetic drug sales") + ylab("$million") + xlab("year")

a10 — Monthly Australian government expenditure on antidiabetic drugs (in $AUD millions), from 1992–2008.
Pattern: Shows a strong upward trend with clear seasonal spikes every January (due to Australian government subsidy system —
prescriptions surge before the new year price policy kicks in).

Seasonal Plot

ggseasonplot(a10, [Link] = TRUE, [Link] = TRUE) + ylab("$million") +

ggtitle("Seasonal plot: antidiabetic drug sales")

ggseasonplot() — Creates a seasonal plot: each year is drawn as a separate coloured line, with months on the X-axis.
[Link] = TRUE — Prints the year label at the right end of each line.
 [Link] = TRUE — Also prints the year label at the left end of each line (visible on both sides).

What the Graph Looks Like

X-axis: Months (January through December)

Y-axis: Sales in $millions
Each coloured line: One year (e.g., a red line for 1992, blue for 1993, etc.)
Pattern: All lines show a dramatic spike in January, a dip in February, then gradual growth through the year. More recent years (lines)
sit higher than older years, confirming the upward trend.

How to read it: If all lines follow a similar seasonal shape (peak and trough in the same months), the seasonal pattern is consistent.
If lines cross randomly, seasonality is weak or changing.

Polar Seasonal Plot

ggseasonplot(a10, polar = TRUE) + ylab("$million") + ggtitle("Seasonal plot: antidiabetic drug sales")

polar = TRUE — Wraps the seasonal plot into a circular (radar) chart where the 12 months are arranged like a clock face.
Use: Visually emphasizes cyclic patterns. More intuitive for some users.

What the Graph Looks Like

The circle represents one full year. January is at the top (12 o'clock position).
Each year is a polygon drawn around the circle.
The January spike appears as a dramatic outward jut at the top.

Seasonal Subseries Plot

ggsubseriesplot(a10) + ylab("$million") + ggtitle("Seasonal sub series plot: antidiabetic drug sales")

ggsubseriesplot() — Creates a seasonal subseries plot: data is separated into 12 mini-panels (one per month), and each
panel shows the value for that month across all years.
A horizontal blue line in each panel shows the mean for that month.

What the Graph Looks Like

12 panels side by side (Jan, Feb, ..., Dec)

Each panel has a line trending upward across years (showing that every month has grown over time)
The January panel has the highest values; February is the lowest
This is useful for detecting whether the seasonal pattern is changing over time — if the January panel is steeper than others,
seasonality is increasing

Scatter Plots for Relationships

autoplot(elecdemand[,c("Demand","Temperature")], facets=TRUE) +
xlab("Year:2014") + ylab("") + ggtitle("Half-hourly electricity demand: Victoria, Australia")

elecdemand — Half-hourly electricity demand and temperature data for Victoria, Australia (year 2014).
facets=TRUE — Produces two separate panels stacked vertically: one for Demand, one for Temperature.
This lets you visually compare whether demand and temperature move together.

qplot(Temperature, Demand, data = [Link](elecdemand)) +

ylab("Demand (GW)") + xlab("Temperature (Celsius)")

qplot() — A "quick plot" function from ggplot2 for simple scatter plots.
[Link](elecdemand) — Converts the ts object to a data frame for qplot.
X-axis: Temperature (°C)
Y-axis: Electricity demand (GW)

What the Graph Looks Like

A U-shaped (non-linear) scatter plot: demand is high at very low temperatures (heating) AND at very high temperatures (air
conditioning), with minimum demand at mild temperatures (~18°C).
This is a classic example of a non-linear relationship — linear regression would not capture this well.

Pair Plots with visnights

autoplot(visnights[,1:5], facets=TRUE)
GGally::ggpairs([Link](visnights[,1:5]))

visnights — Quarterly visitor nights (in millions) in different Australian regions: NSW, VIC, QLD, SA, WA, etc.
visnights[,1:5] — Selects the first 5 region columns.
autoplot(..., facets=TRUE) — 5 separate mini time-plots, one per region, stacked for comparison.
GGally::ggpairs() — Creates a pair plot matrix (also called a scatterplot matrix):
Diagonal: Distribution (density plot) of each variable
Upper triangle: Correlation coefficient between each pair
Lower triangle: Scatter plot of each pair

What the Graph Looks Like

A 5×5 grid of plots

Strong positive correlations (e.g., NSW and VIC might both peak in the same quarters, being major tourism states)
Use: Understand which regions move together (useful if building regression models or studying how regional tourism is linked)

2. Statistical Concepts
Why Visualize First?

Before any modelling, visualization helps you answer:

Is there a trend? (upward, downward, flat)

Is there seasonality? (repeating pattern)
 Are there outliers or anomalies? (sudden spikes, missing data)
Is there a relationship between two variables?

Correlation

Correlation measures the strength of a linear relationship between two variables, ranging from −1 to +1.

+1: Perfect positive linear relationship

0: No linear relationship
−1: Perfect negative linear relationship

Formula (Pearson's r):

r = Σ[(xᵢ − x̄)(yᵢ − ȳ)] / √[Σ(xᵢ − x̄)² × Σ(yᵢ − ȳ)²]

3. Inference & Conclusions

The a10 dataset clearly shows increasing trend + consistent seasonal pattern → decomposition and exponential smoothing
would be appropriate next steps.
The melsyd data shows a disruption — this must be handled before forecasting.
The elecdemand scatter plot reveals a U-shaped non-linear relationship between temperature and demand — a linear model
would be misleading.
Seasonal plots are more informative than simple time plots because they reveal whether the shape of seasonality is changing.
Subseries plots reveal which specific months/quarters are driving the seasonal pattern.

Lab 3 — Forecasting Benchmarks,

Transformations & Residuals
Overview
This lab covers the simplest forecasting methods (used as benchmarks), data transformations (Box-Cox), and the critical concept of
residual diagnostics — how to check whether a forecasting model is good.

Full Code
library(fpp2)
beer2 <- window(ausbeer, start=1992, end=c(2007,4))
beer2
x <- meanf(beer2, h=11)
x
autoplot(x)

autoplot(beer2) +
autolayer(meanf(beer2, h=11), series="Mean", PI=FALSE) +
autolayer(naive(beer2, h=11), series="Naive", PI=FALSE) +
autolayer(snaive(beer2, h=11), series="Seasonal Naive", PI=FALSE) +
autolayer(rwf(beer2, h=11, drift=TRUE), series="Drift", PI=FALSE) +
ggtitle("Forecast for Quarterly Beer Production") +
xlab("Year") + ylab("Megalitres")

b <- rwf(beer2, h=11, drift=TRUE)

b
autoplot(b)

autoplot(elec)
(lambda <- [Link](elec))
autoplot(BoxCox(elec, lambda))

autoplot(goog200)
y <- naive(goog200)
res <- residuals(y)
autoplot(res)
gghistogram(res)
ggAcf(res)

[Link](res, lag=10, fitdf=0)

[Link](res, lag=10, fitdf=0, type="Lj")
checkresiduals(y)

1. Code Walkthrough
Subsetting the Data

beer2 <- window(ausbeer, start=1992, end=c(2007,4))

ausbeer — A built-in dataset in fpp2: quarterly Australian beer production (megalitres) from 1956 to 2010.
window() — Extracts a sub-window (subset) of a ts object between specified start and end dates.
start=1992 — Starts from Q1 of 1992.
end=c(2007,4) — Ends at Q4 of 2007 (the 4th quarter of 2007).
Result: beer2 contains quarterly beer production data from 1992–2007 (used for training models).

Benchmark Forecasting Methods

 x <- meanf(beer2, h=11)

meanf() — Mean Forecast: predicts all future values as the simple historical average of the training data.
h=11 — Forecasts 11 periods ahead (11 quarters = ~2.75 years).
Formula: ŷₜ₊ₕ = (1/T) × Σ yₜ for all t = constant horizontal line

autolayer(naive(beer2, h=11), series="Naive", PI=FALSE)

naive() — Naive Forecast (also called "Random Walk"): predicts that the next value equals the last observed value.
Formula: ŷₜ₊ₕ = yₜ (last known value)
PI=FALSE — Turns off Prediction Intervals (the shaded uncertainty bands) to keep the plot clean.

autolayer(snaive(beer2, h=11), series="Seasonal Naive", PI=FALSE)

snaive() — Seasonal Naive Forecast: predicts that the next value equals the value from the same season last year.
Formula: ŷₜ₊ₕ = yₜ₋ₘ₊ₕ where m = seasonal period (m=4 for quarterly data)
Example: Q3 2008 forecast = actual Q3 2007 value.

autolayer(rwf(beer2, h=11, drift=TRUE), series="Drift", PI=FALSE)

rwf() — Random Walk Forecast. When drift=TRUE, it becomes the Drift Method.
drift=TRUE — Allows the forecast to drift (increase or decrease over time) at the same average rate as in the historical data.
Formula:

ŷₜ₊ₕ = yₜ + h × [(yₜ − y₁) / (T − 1)]

This is equivalent to drawing a straight line from the first point to the last point and extending it.

Summary of Benchmark Methods

Method Formula Best for

Mean Average of all history Stationary, no trend/seasonality
Naive Last observed value Random walk processes (stocks)
Seasonal Naive Same season last year Strong seasonal, no trend
Drift Linear extrapolation of trend Trending data

What the Combined Plot Looks Like

The historical data (1992–2007) is shown as a black line.

Four coloured forecast lines extend beyond 2008, each representing a different method.
The Seasonal Naive line oscillates (correct for seasonal data).
The Mean line is flat.
The Drift line trends slightly downward (if the last period was lower than the first).
The Naive line is flat at the last observed value.

Box-Cox Transformation
 autoplot(elec)
(lambda <- [Link](elec))
autoplot(BoxCox(elec, lambda))

elec — Monthly electricity production in Australia (1956–1995).

[Link]() — Automatically selects the optimal λ (lambda) value for the Box-Cox transformation using maximum likelihood
estimation.

Box-Cox Transformation Formula:

w(t) = { (yₜ^λ − 1) / λ if λ ≠ 0
{ log(yₜ) if λ = 0

λ value Transformation applied

λ = 1 No transformation (original data)
λ = 0.5 Square root transformation
λ = 0 Natural log transformation
λ = −1 Inverse (1/y) transformation

BoxCox(elec, lambda) — Applies the Box-Cox transformation with the chosen λ.

Purpose: Stabilizes the variance of a time series. When the seasonal variation grows over time (heteroscedasticity), transforming
the data makes the variance more constant — which is an assumption of many forecasting models.

Analogy: If the waves in your data get bigger over time (like ocean waves growing near shore), the Box-Cox transformation shrinks
them back to a uniform size — making the data easier to model.

What the Graphs Look Like

Before (autoplot(elec)): The seasonal peaks grow taller as time progresses — classic multiplicative behaviour.
After (autoplot(BoxCox(elec, lambda))): The seasonal peaks are now more uniform in height — the variance has been
stabilized.

Residual Diagnostics

y <- naive(goog200)
res <- residuals(y)

goog200 — Daily closing prices of Google stock (first 200 trading days).
naive(goog200) — Applies the naive forecasting method (fitted values = previous day's price).
residuals() — Computes residuals: the difference between the actual value and the fitted (predicted) value.
Formula: eₜ = yₜ − ŷₜ

Residual Plot

autoplot(res)
 What the Graph Looks Like

X-axis: Time
Y-axis: Residual values (centred around 0)
Ideal pattern: Residuals should look like "white noise" — randomly scattered around zero with no visible pattern, no trend, no
seasonality, and no clustering.
Bad pattern: If you see a wave or trend in the residuals, the model has missed something important in the data.

Histogram of Residuals

gghistogram(res)

What the Graph Looks Like

A bar chart showing the frequency distribution of residuals.

Ideal: Bell-shaped (normal distribution), centred at 0.
If skewed: The model is systematically over- or under-forecasting.
Why it matters: Many forecasting models assume residuals are normally distributed — a non-normal distribution can affect
prediction interval accuracy.

ACF Plot of Residuals

ggAcf(res)

ggAcf() — Plots the ACF (AutoCorrelation Function) of the residuals.

ACF measures the correlation between a series and its own lagged values.

ACF Formula:

rₖ = Σ[(yₜ − ȳ)(yₜ₋ₖ − ȳ)] / Σ(yₜ − ȳ)²

Where k = lag (number of periods apart).

What the Graph Looks Like

X-axis: Lag number (1, 2, 3, ..., typically up to 20 or 30)

Y-axis: Correlation coefficient (−1 to +1)
Each lag has a vertical bar (spike).
Blue dashed lines: Confidence bounds (approximately ±1.96/√T). Spikes INSIDE these lines are NOT statistically significant.
Ideal (good) pattern: ALL bars inside the blue lines → no autocorrelation → residuals are white noise.
Bad pattern: Bars sticking outside the blue lines → significant autocorrelation → the model has missed a pattern.

Box-Pierce and Ljung-Box Tests

[Link](res, lag=10, fitdf=0)

[Link](res, lag=10, fitdf=0, type="Lj")

These are formal hypothesis tests for autocorrelation in residuals.

 Hypotheses:

H₀ (Null Hypothesis): No autocorrelation in the residuals (residuals are white noise)

H₁ (Alternative Hypothesis): Autocorrelation exists in the residuals

lag=10 — Tests the first 10 lags collectively.

fitdf=0 — Degrees of freedom used by the model. For a naive model, fitdf=0 (no parameters estimated).

type="Lj" — Uses the Ljung-Box version (more accurate for small samples than Box-Pierce).

Decision rule:

If p-value > 0.05 → Fail to reject H₀ → Residuals are white noise → Model is acceptable
If p-value < 0.05 → Reject H₀ → Autocorrelation remains → Model needs improvement

All-in-one Diagnostic

checkresiduals(y)

checkresiduals() — Produces three plots in one:

1. Residual time plot
2. ACF plot of residuals
3. Histogram of residuals
Also automatically runs the Ljung-Box test and prints the p-value.

3. Inference & Conclusions

For the Google stock data, the naive forecast is appropriate (stock prices behave like a random walk).
Residuals from the naive model should be white noise — if the ACF shows all bars within bounds, the model has captured all
predictable structure.
Box-Cox transformation is essential when seasonal variance is non-constant — it makes the data better suited for additive
decomposition or ARIMA modelling.
Benchmark methods (Mean, Naive, Seasonal Naive, Drift) are not necessarily good forecasters — but they serve as benchmarks: if
your complex model can't beat a simple naive forecast, something is wrong.

Lab 4 — Forecast Accuracy & Time Series

Regression
Overview
This lab covers two major topics:

1. How to measure forecast accuracy using error metrics

2. How to build regression models for time series — using other variables to explain and forecast a series

Full Code
library(fpp2)
autoplot(ausbeer)
beer2 <- window(ausbeer, start=1992, end=c(2007,4))
beerfit1 <- meanf(beer2, h=10)
beerfit2 <- rwf(beer2, h=10)
beerfit3 <- snaive(beer2, h=10)
beer3 <- window(ausbeer, start=2008)
beer3
accuracy(beerfit1, beer3)
accuracy(beerfit2, beer3)
accuracy(beerfit3, beer3)

uschange
autoplot(uschange[,c("Consumption","Income")]) + ylab("Percentage change") + xlab("year")
tslm(Consumption ~ Income, data = uschange)
uschange %>% [Link] %>% GGally::ggpairs()
[Link] <- tslm(Consumption ~ Income + Production + Savings, data=uschange)
summary([Link])
checkresiduals([Link])

beer2
fit_beer <- tslm(beer2 ~ trend + season)
summary(fit_beer)
autoplot(beer2)

1. Code Walkthrough
Training vs Test Set

beer2 <- window(ausbeer, start=1992, end=c(2007,4)) # Training set

beer3 <- window(ausbeer, start=2008) # Test set

Training set (beer2): Data used to build the model (1992–2007).

Test set (beer3): Data used to evaluate how well the model performs on unseen data (2008 onwards).
Critical rule: The test set must NEVER be used during model training — it must remain "unseen" until evaluation.

Analogy: It's like studying for an exam (training) and then taking a different exam (test) to see how well you truly learned.

Forecast Accuracy Metrics

accuracy(beerfit1, beer3)
accuracy(beerfit2, beer3)
accuracy(beerfit3, beer3)

accuracy(forecast_object, actual_test_data) — Computes multiple error metrics comparing the forecasts to the actual
test values.
 The output shows two rows:

Training set: How well the model fits the historical data it was trained on.
Test set: How well it forecasts the future (the important row!).

Key Metrics Explained:

Metric Full Name Formula Interpretation

Average bias; +ve = over-
ME Mean Error (1/n)Σeₜ
forecast, −ve = under-forecast
Root Mean Most common; penalizes large
RMSE √[(1/n)Σeₜ²]
Squared Error errors heavily
Average absolute
Mean Absolute
MAE (1/n)Σ eₜ error; easy to
Error
interpret
Mean Percentage
MPE (1/n)Σ(eₜ/yₜ)×100 Average % bias
Error
Mean Absolute Average % error;
MAPE (1/n)Σ eₜ/yₜ ×100
Percentage Error scale-independent
Mean Absolute MAE / Compares to naive forecast; <1
MASE
Scaled Error MAE_naive means better than naive
First-order ACF of Lag-1
ACF1 Should be close to 0
residuals autocorrelation

Where eₜ = yₜ − ŷₜ (actual minus forecast)

Best model = lowest RMSE and MAE on the test set (not the training set).

For the beer data:

beerfit3 (Seasonal Naive) should outperform the others because beer production has strong quarterly seasonality.

uschange Dataset — Multiple Variables

uschange
autoplot(uschange[,c("Consumption","Income")]) + ylab("Percentage change") + xlab("year")

uschange — Quarterly percentage changes in US consumption, income, production, savings, and unemployment (1970–2016).
Shows all variables are percentage changes (stationary — see Lab 7), making them suitable for regression.

Simple Linear Regression (TSLM)

tslm(Consumption ~ Income, data = uschange)

tslm() — Time Series Linear Model: fits an OLS regression to time series data. Similar to lm() but handles ts objects properly.
Formula: Consumption = β₀ + β₁ × Income + εₜ
~ — The tilde means "is modelled by" or "depends on".
This gives coefficient estimates (β₀ = intercept, β₁ = slope).

Pair Plot for Multiple Variables

uschange %>% [Link] %>% GGally::ggpairs()

%>% — The pipe operator (from magrittr/dplyr): passes the output of one function as the first input to the next. Reads as "then".
Converts to data frame and passes to ggpairs() for a full pair plot matrix.
Use: See correlations between ALL pairs of variables before building a multiple regression.

Multiple Regression

[Link] <- tslm(Consumption ~ Income + Production + Savings, data=uschange)

summary([Link])

Multiple Regression Formula:

Consumptionₜ = β₀ + β₁×Incomeₜ + β₂×Productionₜ + β₃×Savingsₜ + εₜ

summary() — Prints the full regression output including:

Coefficients (β): The estimated impact of each predictor
Standard Errors: Precision of each coefficient estimate
t-value: Coefficient / Standard Error — higher absolute value = more significant
p-value: If < 0.05, the predictor is statistically significant
R²: Proportion of variance in Consumption explained by the model (0 to 1)
Adjusted R²: R² penalized for number of predictors (use this for model comparison)
F-statistic: Tests if the overall model is significant

Interpreting Coefficients:

β₁ = 0.75 for Income → "For every 1% increase in income, consumption increases by 0.75%, holding other variables constant."

checkresiduals([Link])

Checks if the residuals of the regression model are white noise (no autocorrelation, normally distributed, centred at 0).
Why important: Autocorrelated residuals in regression violate OLS assumptions and make the standard errors unreliable.

Trend and Season in Regression

fit_beer <- tslm(beer2 ~ trend + season)

summary(fit_beer)

trend — A special variable in tslm() that creates a numeric time index (1, 2, 3, ..., T) to capture the linear trend.
season — A special variable in tslm() that creates dummy variables for each season (quarter in this case).
For quarterly data: it creates 3 dummy variables (Q2, Q3, Q4 — with Q1 as the base/reference category).
Full model:

Beerₜ = β₀ + β₁×t + β₂×D₂ + β₃×D₃ + β₄×D₄ + εₜ

Where D₂, D₃, D₄ are 1 if the observation is in Q2, Q3, Q4 respectively, 0 otherwise.

Interpretation of season coefficients:

 β₂ for Q2 = the average difference between Q2 and Q1 production (in megalitres), after removing the trend.
Negative β values for Q2 suggest Q2 production is lower than Q1.

3. Inference & Conclusions

Seasonal Naive beats Mean and Random Walk for beer production — confirming strong quarterly seasonality.
MASE < 1 for Seasonal Naive confirms it outperforms a naive benchmark.
The uschange regression shows that income is the most important predictor of consumption.
Savings is negatively correlated with consumption (people save less when they spend more).
The tslm(beer2 ~ trend + season) model confirms both a negative trend (declining production) and strong seasonal
differences between quarters.
Residual checks for regression models are as important as for pure forecasting models — autocorrelated residuals signal model
misspecification.

Lab 5 — Moving Averages & Decomposition

Methods
Overview
This lab covers moving averages (a data smoothing technique) and four different methods for decomposing a time series: Classical, X11,
SEATS, and STL. Each has different strengths, and the choice depends on the data characteristics.

Full Code
library(fpp2)
autoplot(elecsales)
ma(elecsales, 5)
autoplot(elecsales) + autolayer(ma(elecsales, 5))

beer2 <- window(ausbeer, start=1992)

autoplot(beer2)
ma4 <- ma(beer2, 4)
ma4 <- ma(beer2, 4, centre = FALSE)
ma2x4 <- ma(beer2, order = 4, centre = TRUE)
autoplot(beer2) + autolayer(ma4) + autolayer(ma2x4)

# Classical Decomposition
elecequip %>% decompose(type = "multiplicative") %>% autoplot()
y <- elecequip %>% decompose(type = "multiplicative")
y$trend
y$seasonal
y$figure

# X11
[Link]("seasonal")
library(seasonal)
elecequip %>% seas(x11="") -> fit
autoplot(fit)
trendcycle(fit)
seasonal(fit)

# SEATS
elecequip %>% seas() %>% autoplot()
y1 <- elecequip %>% seas()
trendcycle(y1)
seasonal(y1)

# STL
elecequip %>% stl([Link]=13, [Link]="periodic", robust=TRUE) %>% autoplot()
y2 <- elecequip %>% stl([Link]=13, [Link]="periodic", robust=TRUE)
y3 <- elecequip %>% stl([Link]=13, [Link]=7, robust=TRUE)
trendcycle(y2)
seasonal(y3)

forecasting <- stlf(elecequip, method='naive')

forecasting

1. Code Walkthrough
Moving Averages
ma(elecsales, 5)
autoplot(elecsales) + autolayer(ma(elecsales, 5))

elecsales — Annual electricity sales in South Australia (1989–1998).

ma(series, order) — Computes a Moving Average (MA) of the specified order.
Order = 5 means: each smoothed value is the average of 5 consecutive observations (the current one, 2 before, and 2 after —
centred).

MA(5) Formula:

MAₜ = (yₜ₋₂ + yₜ₋₁ + yₜ + yₜ₊₁ + yₜ₊₂) / 5

Effect: Smooths out short-term fluctuations to reveal the underlying trend.

Trade-off: Higher order → smoother but loses more data at the ends (first and last few values will be NA).

Analogy: Imagine you're tracking your daily weight. A 7-day moving average smooths out the day-to-day noise to show your actual
weight trend.

What the Graph Looks Like

Original series: Jagged line with ups and downs.

Moving average: A smoother curve that follows the general direction of the data.
The MA line doesn't extend to the very beginning or end (missing values at edges).

Moving Average for Seasonal Data

ma4 <- ma(beer2, 4, centre = FALSE)

ma2x4 <- ma(beer2, order = 4, centre = TRUE)

centre = FALSE (default) — Simple 4-period moving average. For quarterly data (m=4), each point is the average of 4 consecutive
quarters. BUT this produces a MA centred between two time points (not at an actual observation).
centre = TRUE — 2×4 Moving Average: applies a MA(4) twice, which re-centres it at actual time points. This is the standard
method for even-period seasonal data.

2×4 MA Formula:

MA(2×4) = (0.5yₜ₋₂ + yₜ₋₁ + yₜ + yₜ₊₁ + 0.5yₜ₊₂) / 4

(The endpoints get a weight of 0.5 instead of 1.)

Why 2×4? For data with an even seasonal period (4 quarters, 12 months), you need 2×m to avoid a phase shift in the smoothed
series.

Method 1: Classical Decomposition

elecequip %>% decompose(type = "multiplicative") %>% autoplot()
y$figure

elecequip — Monthly new orders for electrical equipment in the Euro area (1996–2012).
$figure — Extracts the seasonal indices — the average seasonal pattern (12 values for monthly data).

Classical Decomposition Algorithm:

1. Estimate trend using a centred moving average (MA)

2. Remove trend: For multiplicative: Sₜ × Rₜ = Yₜ / Tₜ
3. Average the de-trended values by season to get seasonal indices
4. Remove seasonal: Rₜ = Yₜ / (Tₜ × Sₜ)

Limitations of Classical Decomposition:

Trend estimates are unavailable for the first and last few periods (NA values)
Assumes the seasonal component is the same every year (cannot handle changing seasonality)
Not robust to outliers

Method 2: X11 Decomposition

library(seasonal)
elecequip %>% seas(x11="") -> fit
autoplot(fit)

seasonal — R package that interfaces with the US Census Bureau's X-13ARIMA-SEATS software.
seas(x11="") — Activates the X11 decomposition method (the empty string "" tells it to use X11 defaults).
-> — Alternative right-assignment operator (same as fit <- elecequip %>% seas(x11=""))

X11 Method:

Originally developed by the US Census Bureau in the 1960s, refined as X-11, X-12, X-13.
Uses iterative centred moving averages with various filters.
Advantages over Classical:
Trend estimates available at the start AND end of the series
Handles varying seasonal patterns over time
More robust to outliers and irregular observations
Can handle both additive and multiplicative models

trendcycle(fit)
seasonal(fit)

trendcycle() — Extracts the trend-cycle component (combined smooth trend + business cycle).
seasonal() — Extracts the seasonal component.

What the X11 Plot Looks Like

Similar 4-panel layout to classical decomposition (Observed, Trend, Seasonal, Irregular).

Key difference: The seasonal panel may not be perfectly regular — it can show slight variation from year to year (capturing changing
seasonality).
The trend extends all the way to the first and last observations (no NA at edges).
 Method 3: SEATS Decomposition

elecequip %>% seas() %>% autoplot()

seas() without x11="" uses SEATS (Seasonal Extraction in ARIMA Time Series) — the default method in the seasonal package.
SEATS is a model-based approach:
1. Fits an ARIMA model to the data (see Lab 8)
2. Derives the decomposition components from that ARIMA model analytically
Advantages: Statistically grounded, produces smooth components, gives measures of uncertainty.
Limitation: Only works if a suitable ARIMA model can be fitted.

Method 4: STL Decomposition

elecequip %>% stl([Link]=13, [Link]="periodic", robust=TRUE) %>% autoplot()

STL — Seasonal and Trend decomposition using Loess

LOESS (Locally Estimated Scatterplot Smoothing): A non-parametric smoothing technique that fits local polynomial regressions
in a moving window.

Parameters:

[Link]=13 — Controls the smoothness of the trend component. Larger = smoother trend. Must be odd. 13 means a 13-period
local regression window.
[Link]="periodic" — Controls the smoothness of the seasonal component.
"periodic" → The seasonal component is forced to be identical every year (perfectly periodic, no change over time). This
is the most restrictive setting.
A number (e.g., [Link]=7) → Allows the seasonal component to change slowly over time. Larger number = slower
change. Minimum value = 7.
robust=TRUE — Makes the decomposition robust to outliers by downweighting extreme observations in the LOESS fit.
Recommended for real-world data.

STL Advantages:

Can handle ANY type of seasonality (monthly, quarterly, weekly, daily)

The seasonal component CAN change over time (unlike classical)
Robust to outliers (with robust=TRUE)
Works only with additive decomposition (but data can be log-transformed first for multiplicative behaviour)

STL Disadvantages:

Does not handle trading day or calendar effects automatically (unlike X11/SEATS)

y3 <- elecequip %>% stl([Link]=13, [Link]=7, robust=TRUE)

[Link]=7 → The seasonal pattern is now allowed to change slowly over time (7-period local seasonal window).

Forecasting with STL

forecasting <- stlf(elecequip, method='naive')

 stlf() — STL Forecast: first decomposes using STL, seasonally adjusts the data, then forecasts the seasonally-adjusted data
using a specified method, and finally re-adds the seasonal component.
method='naive' — Uses the naive method on the seasonally-adjusted data.
Other options: method='ets' (Exponential Smoothing), method='arima'.

2. Decomposition Method Comparison

Feature Classical X11 SEATS STL
Additive or
Seasonal type Both Both (via ARIMA) Additive only*
Multiplicative
Changing Yes (with [Link]
No Yes Yes
seasonality number)
Robust to outliers No Somewhat No Yes (robust=TRUE)
Trend at edges Missing (NA) Available Available Available
Works with any Monthly/Quarterly Monthly/Quarterly
Yes Yes (any frequency)
frequency only only
Complexity Simple High High Moderate

*Use log transformation before STL for multiplicative data.

3. Inference & Conclusions

Moving averages reveal the trend by smoothing seasonal and random noise — higher order = smoother but less responsive.
For quarterly data, 2×4 MA is the correct choice to centre the trend estimate.
Classical decomposition is the simplest but assumes constant seasonality — often too restrictive for real data.
STL is the most flexible and recommended for most applications, especially when seasonality may be changing.
X11 and SEATS are industry standards (used by government statistical agencies) but limited to monthly/quarterly data.
The stlf() function makes it easy to forecast after decomposition — combine the best of decomposition and forecasting.

Lab 6 — Exponential Smoothing Methods & ETS

Models
Overview
Exponential Smoothing is one of the most widely used and practically successful forecasting approaches. Unlike the simple moving
average (which gives equal weight to all past observations), exponential smoothing gives more weight to recent observations and less
weight to older ones — the weights decay exponentially as you go further back in time.

This lab covers three flavours: Simple (SES), Holt's (Double), Holt-Winters (Triple), and the fully automatic ETS framework.

Full Code
library(fpp2)
oildata <- window(oil, start=1996)
autoplot(oildata)

# Simple Exponential Smoothing

fc <- ses(oildata, h=5)
summary(fc)
round(accuracy(fc), 2)
autoplot(oildata) + autolayer(fc)
autoplot(oildata) + autolayer(fc, PI=FALSE)
fc$fitted
autoplot(oildata) + autolayer(fitted(fc))

# Holt's Method (Double Exponential Smoothing - Additive)

air <- window(ausair, start=1990)
fc1 <- holt(air, h=5)
summary(fc1)
autoplot(ausair) + autolayer(fc1, PI=FALSE)
fc1$fitted
autoplot(ausair) + autolayer(fitted(fc1))

# Holt's Damped Method

fc2 <- holt(air, damped=TRUE, h=5)
summary(fc2)
autoplot(ausair) + autolayer(fc2, PI=FALSE)

fc3 <- holt(air, damped=TRUE, phi=0.9, h=5)

summary(fc3)
autoplot(ausair) + autolayer(fc3, PI=FALSE)

# Time Series Cross-Validation

autoplot(livestock)
e1 <- tsCV(livestock, ses, h=1)
e2 <- tsCV(livestock, holt, h=1)
e3 <- tsCV(livestock, holt, damped=TRUE, h=1)
mean(e1^2, [Link]=TRUE)
mean(e2^2, [Link]=TRUE)
mean(e3^2, [Link]=TRUE)

# Holt-Winters (Triple Exponential Smoothing)

aust <- window(austourists, start=2005)
autoplot(austourists)
fit1 <- hw(aust, seasonal="additive")
fit2 <- hw(aust, seasonal="multiplicative")
autoplot(austourists) + autolayer(fit1) + autolayer(fit2)
summary(fit1)
summary(fit2)
round(accuracy(fit1), 2)
round(accuracy(fit2), 2)
# ETS (Error, Trend, Seasonality)
fit <- ets(aust)
summary(fit)

gc()

1. Code Walkthrough
Simple Exponential Smoothing (SES)

oildata <- window(oil, start=1996)

fc <- ses(oildata, h=5)
summary(fc)

oil — Annual oil production (millions of tonnes), Saudi Arabia (1965–2013).

window(oil, start=1996) — Subsets from 1996 onward.
ses() — Simple Exponential Smoothing: the most basic form, suitable for data with no trend and no seasonality.
h=5 — Forecast 5 periods (5 years) ahead.

SES Formula — Level Equation:

Level: lₜ = α × yₜ + (1 − α) × lₜ₋₁

Forecast: ŷₜ₊ₕ = lₜ (flat forecast — same for all h)

α (alpha) — The smoothing parameter (0 < α < 1).

α close to 1: Heavy weight on the most recent observation → model reacts quickly to changes.
α close to 0: Heavy weight on older observations → model is slow to adapt, very smooth forecast.
α is automatically optimized by minimizing the Sum of Squared Errors (SSE).

Analogy: SES is like a person who mostly listens to what just happened (α close to 1) vs. someone who relies heavily on long-term
memory (α close to 0).

Expanding the formula:

lₜ = α × yₜ + α(1−α) × yₜ₋₁ + α(1−α)² × yₜ₋₂ + ...

The weights α, α(1−α), α(1−α)², ... sum to 1 and decay exponentially.

fc$fitted
autoplot(oildata) + autolayer(fitted(fc))

fc$fitted — The in-sample fitted values (the model's estimates of past observations).
fitted(fc) — Equivalent way to extract fitted values.
The fitted line will look like a smoothed version of the original, always lagging slightly behind.
 What the SES Forecast Plot Looks Like

Historical data is shown as a solid line.

The forecast (from the last observed year onward) is a flat horizontal line (because SES has no trend component).
Shaded regions show prediction intervals (80% and 95% confidence bands that widen as you forecast further).

Holt's Method (Double Exponential Smoothing)

air <- window(ausair, start=1990)

fc1 <- holt(air, h=5)
summary(fc1)

ausair — Annual passengers (millions) on Australian airlines.

holt() — Holt's Linear Method (also called Double Exponential Smoothing): extends SES to handle data with a linear trend.

Holt's Method Formula (Additive trend):

Level: lₜ = α × yₜ + (1 − α)(lₜ₋₁ + bₜ₋₁)

Trend: bₜ = β × (lₜ − lₜ₋₁) + (1 − β) × bₜ₋₁
Forecast: ŷₜ₊ₕ = lₜ + h × bₜ

lₜ = Level (smoothed estimate of the current value)

bₜ = Trend (smoothed estimate of the current slope)
α (alpha) = Level smoothing parameter (0 < α < 1)
β (beta) = Trend smoothing parameter (0 < β < 1)
h = Forecast horizon

Interpretation: The forecast is a straight line extending from the last known point, with slope bₜ.

What the Holt Forecast Plot Looks Like

Historical data shows an upward trend.

The forecast is a straight line continuing upward (linear extrapolation of the trend).
Prediction intervals widen significantly over time.

Holt's Damped Method

fc2 <- holt(air, damped=TRUE, h=5)

fc3 <- holt(air, damped=TRUE, phi=0.9, h=5)

damped=TRUE — Adds a damping parameter φ (phi) that gradually flattens the trend toward a horizontal line as the forecast
horizon increases.
φ (phi) — Damping parameter (0 < φ < 1).
φ = 1: No damping (equivalent to standard Holt's).
φ = 0: Trend immediately collapses to zero (equivalent to SES).
Typical range: 0.8 to 0.98

Damped Holt's Formula:

Level: lₜ = α × yₜ + (1 − α)(lₜ₋₁ + φ × bₜ₋₁)
Trend: bₜ = β × (lₜ − lₜ₋₁) + (1 − β) × φ × bₜ₋₁
Forecast: ŷₜ₊ₕ = lₜ + (φ + φ² + ... + φʰ) × bₜ

phi=0.9 — Manually sets φ = 0.9 (instead of letting R optimize it). This is a hyperparameter — it controls model behaviour but is
not estimated from the data in the usual sense.

What the Damped Forecast Looks Like

The forecast starts as a trend line but curves and flattens as it moves further into the future.
At large h, the forecast asymptotically approaches a constant level.
This is often more realistic than assuming a trend continues forever.

Real-world insight: Most things that have been growing don't grow forever at the same rate. The damped method embeds this
realism into the forecast.

Time Series Cross-Validation (k-fold CV)

e1 <- tsCV(livestock, ses, h=1)

e2 <- tsCV(livestock, holt, h=1)
e3 <- tsCV(livestock, holt, damped=TRUE, h=1)
mean(e1^2, [Link]=TRUE)
mean(e2^2, [Link]=TRUE)
mean(e3^2, [Link]=TRUE)

livestock — Annual sheep population in Asia (in millions).

tsCV() — Time Series Cross-Validation (also called rolling-origin evaluation or walk-forward validation):

Train the model on data up to time t, forecast h steps ahead, compute the error.
Move forward one step, retrain, re-forecast.
Repeat for all valid time points.
This is the time-series equivalent of k-fold cross-validation.

mean(e1^2, [Link]=TRUE) — Computes the Mean Squared Error (MSE) across all cross-validation errors. [Link]=TRUE
removes NA values (the first few periods where there isn't enough data to train).

Decision: The method with the lowest MSE from cross-validation is selected.

Why not just use training set accuracy? In-sample accuracy can be misleading — a complex model might perfectly fit the training
data but fail on new data (overfitting). Cross-validation gives a more honest estimate.

Holt-Winters' Method (Triple Exponential Smoothing)

aust <- window(austourists, start=2005)
fit1 <- hw(aust, seasonal="additive")
fit2 <- hw(aust, seasonal="multiplicative")

austourists — Quarterly international tourist visitor nights in Australia (2000–2015).

hw() — Holt-Winters' Method: extends Holt's method to also handle seasonality. Three components: Level (l), Trend (b), and
Seasonal (s).

Holt-Winters Additive Model:

Level: lₜ = α(yₜ − sₜ₋ₘ) + (1 − α)(lₜ₋₁ + bₜ₋₁)

Trend: bₜ = β(lₜ − lₜ₋₁) + (1 − β) bₜ₋₁
Seasonal: sₜ = γ(yₜ − lₜ₋₁ − bₜ₋₁) + (1 − γ) sₜ₋ₘ
Forecast: ŷₜ₊ₕ = lₜ + h×bₜ + sₜ₋ₘ₊ₕₘ₋

Holt-Winters Multiplicative Model:

Level: lₜ = α(yₜ / sₜ₋ₘ) + (1 − α)(lₜ₋₁ + bₜ₋₁)

Trend: bₜ = β(lₜ − lₜ₋₁) + (1 − β) bₜ₋₁
Seasonal: sₜ = γ(yₜ / (lₜ₋₁ + bₜ₋₁)) + (1 − γ) sₜ₋ₘ
Forecast: ŷₜ₊ₕ = (lₜ + h×bₜ) × sₜ₋ₘ₊ₕₘ₋

Where:

α (alpha): Level smoothing (0 < α < 1)

β (beta): Trend smoothing (0 < β < 1)
γ (gamma): Seasonal smoothing (0 < γ < 1)
m: Seasonal period (m=4 for quarterly, m=12 for monthly)
sₜ₋ₘ: The seasonal component from the same season last year

Additive vs Multiplicative Seasonality:

Additive: Seasonal effect is the same size regardless of the level (fixed number of extra visitors per season)
Multiplicative: Seasonal effect scales with the level (proportionally more visitors in peak season as total grows)

What the Holt-Winters Plot Looks Like

Two coloured forecast lines extend into the future.

Both show an upward trend with seasonal oscillations.
The multiplicative forecast has larger seasonal swings in the future (because the seasonal effect grows with the trend).
Both have prediction interval bands.

ETS Framework (Automatic Model Selection)

fit <- ets(aust)

summary(fit)

ets() — Error, Trend, Seasonality model selection. This is a systematic framework that:
1. Considers ALL possible combinations of Error type, Trend type, and Seasonality type.
2. Fits each model.
 3. Selects the best using AIC (Akaike Information Criterion).

ETS Notation: ETS(E, T, S)

Position Options Meaning

E — Error A (Additive), M (Multiplicative) How forecast errors enter the model
T — Trend N (None), A (Additive), Ad (Additive Damped) Trend component
S — Seasonal N (None), A (Additive), M (Multiplicative) Seasonal component

Common ETS models and their equivalents:

ETS Model Equivalent to

ETS(A,N,N) Simple Exponential Smoothing (SES)
ETS(A,A,N) Holt's Linear Method
ETS(A,Ad,N) Holt's Damped Method
ETS(A,A,A) Holt-Winters Additive
ETS(M,A,M) Holt-Winters Multiplicative

summary(fit) — Prints the selected model (e.g., "ETS(M,A,M)"), its parameters (α, β, γ, φ), initial states, and AIC/AICc/BIC values.

Comment from the code: "determined that error is Multiplicative, trend additive and seasonal multiplicative" — so the best model is
ETS(M,A,M), which is equivalent to a Holt-Winters multiplicative model.

gc()

gc() — Garbage Collection: frees up unused memory in R. Good practice after memory-intensive operations.

3. Inference & Conclusions

SES is appropriate for the oil data (no clear trend).
Holt's method captures the linear trend in airline passengers.
Damped trend is often more realistic for medium-to-long-term forecasts — trends rarely continue at the same rate indefinitely.
Cross-validation objectively identifies the best method: if mean(e3^2) is lowest, the damped Holt's model is best.
ETS(M,A,M) for the tourism data means: Multiplicative Error + Additive Trend + Multiplicative Seasonality — the seasonal pattern
grows with the trend (common in tourism data).
Comparing accuracy(fit1) vs accuracy(fit2): whichever has lower RMSE/MAE on the test set is preferred.

Lab 7 — Stationarity, Differencing & ACF/PACF

Overview
This lab is a critical gateway to ARIMA modelling. Before fitting an ARIMA model, you MUST ensure the data is stationary. This lab teaches
you what stationarity is, how to test for it visually, how to achieve it through differencing, and how to use ACF/PACF plots to guide model
selection.

Full Code
library(fpp2)

autoplot(lynx)
acf(goog200)
[Link](goog200)
y <- diff(goog200)
acf(y)
[Link](y)
ndiffs(goog200)
nsdiffs(goog200)

autoplot(a10)
y <- log(a10)
y1 <- diff(y)
cbind(a10, y, y1) %>% autoplot(facets=TRUE)
ndiffs(y1)

autoplot(usmelec)
d <- log(usmelec)
nsdiffs(usmelec)
nsdiffs(d)
d1 <- diff(d)
ndiffs(d1)
cbind(usmelec, d, d1) %>% autoplot(facets=TRUE)

par(mfrow=c(3,3))
plot(goog200, main="GOOG200")
plot(diff(goog200), main="Diff GOOG200")
plot(strikes, main="Strikes")
plot(hsales, main="House Sales")
plot(eggs, main="Eggs")
plot(pigs, main="Pigs")
plot(lynx, main="Lynx")
plot(beer, main="Beer")
plot(elec, main="Electricity")

1. Code Walkthrough
Visualizing the Lynx Data

autoplot(lynx)

lynx — Annual number of lynx trapped in Canada (1821–1934).

Pattern: Shows a remarkable cyclical pattern (peaks roughly every 10 years) — driven by the predator-prey cycle between lynx and
snowshoe hares.
Note: This is a cyclic pattern, NOT seasonal — cycles don't have a fixed period (they vary from 7–14 years).

ACF of Google Stock Price

acf(goog200)
[Link](goog200)

acf() — Base R version of the ACF plot (similar to ggAcf() from Lab 3).
For goog200 (Google stock PRICE), the ACF will show slowly decaying, significant correlations at all lags — a classic sign of
non-stationarity.
[Link](goog200) — Tests for autocorrelation in the raw series. With a very small p-value, we reject H₀ and conclude there IS
significant autocorrelation → the series is not white noise → likely non-stationary.

Differencing to Achieve Stationarity

y <- diff(goog200)
acf(y)
[Link](y)

diff() — Computes first-order differences: yₜ' = yₜ − yₜ₋₁

This transforms a series of prices into a series of changes in price (returns).
After differencing goog200:
acf(y) — Should show NO significant spikes beyond lag 0 → white noise.
[Link](y) — p-value should be > 0.05 → fail to reject H₀ → no autocorrelation.
This confirms that Google STOCK PRICE is non-stationary, but STOCK RETURNS (differenced prices) are stationary.

How Many Differences Are Needed?

ndiffs(goog200)
nsdiffs(goog200)

ndiffs() — Uses statistical tests (KPSS or ADF) to determine the number of regular (non-seasonal) differences needed to
make the series stationary.
Returns 1 for goog200 → one first difference is sufficient.
nsdiffs() — Determines the number of seasonal differences needed.
Returns 0 for goog200 (annual data — no seasonality).

Types of Differencing:

Regular (first-order) difference: yₜ' = yₜ − yₜ₋₁ (removes trend)

Second-order difference: yₜ'' = yₜ' − yₜ'₋₁ = yₜ − 2yₜ₋₁ + yₜ₋₂ (removes quadratic trend)
Seasonal difference (lag m): yₜ' = yₜ − yₜ₋ₘ (removes seasonal pattern)

Transforming and Differencing the a10 Drug Sales Data

autoplot(a10)
y <- log(a10)
y1 <- diff(y)
cbind(a10, y, y1) %>% autoplot(facets=TRUE)
ndiffs(y1)

a10 has 3 issues: trend (upward), seasonality (monthly), and growing variance.
 Step 1: log(a10) — Log transformation stabilizes the growing variance (as noted in Lab 3 with Box-Cox, log is Box-Cox with λ=0).
Step 2: diff(y) (default lag=1) — Takes the first difference to remove the trend AND, since this is monthly data with lag=12
implicitly captured in step 1, here one difference step makes it more stationary.

Note: For seasonal data, you typically need:

1. A seasonal difference (diff(y, lag=12) for monthly) to remove seasonality.

2. A regular difference (diff(...)) to remove remaining trend.

cbind(a10, y, y1) — Combines multiple ts objects side by side. Used here because all three have different units/scales.
autoplot(facets=TRUE) — Plots each column in a separate panel (since they have different Y-axis scales). This is the right
choice when comparing original, log, and differenced data.

Key comment in code: "when y axes are different use cbind" — this is the correct approach for comparing series at different scales.

US Monthly Electricity (usmelec)

autoplot(usmelec)
d <- log(usmelec)
nsdiffs(usmelec) # returns 1 → needs 1 seasonal difference
nsdiffs(d)
d1 <- diff(d)
ndiffs(d1)
cbind(usmelec, d, d1) %>% autoplot(facets=TRUE)

usmelec — US monthly electricity generation (billion kWh, 1973–2013).

Has trend, monthly seasonality, and growing variance.
nsdiffs(usmelec) → 1 (needs a seasonal difference)
d <- log(usmelec) — Stabilizes variance.
d1 <- diff(d) — Regular differencing after logging.
ndiffs(d1) — Checks if further differencing is needed after one seasonal+regular difference sequence.

Multi-Panel Plot of 9 Time Series

par(mfrow=c(3,3))
plot(goog200, main="GOOG200")
...
plot(elec, main="Electricity")

par(mfrow=c(3,3)) — Sets up the plotting area as a 3-row × 3-column grid for 9 plots. par() controls base R graphical
parameters; mfrow = "multiple frames by row".
This is the base R equivalent of facets=TRUE in ggplot2.
When to use par(mfrow) vs cbind(...) + facets=TRUE:
par(mfrow) → when X-axes are different (different time periods)
cbind() + facets=TRUE → when Y-axes are different (different scales but same time range)
 2. Statistical Concepts
What is Stationarity?

A time series is stationary if its statistical properties do NOT change over time:

1. Constant mean: The series fluctuates around a fixed level

2. Constant variance: The spread of values doesn't change over time
3. Constant autocovariance: The correlation between yₜ and yₜ₋ₖ depends only on k, not on t

Why does stationarity matter?

ARIMA models (and most time series models) assume stationarity.

Non-stationary series can produce spurious regressions — apparent relationships that are actually meaningless.

Types of non-stationarity:

Trend non-stationarity: Mean changes over time → fix with regular differencing
Seasonal non-stationarity: Seasonal pattern present → fix with seasonal differencing
Variance non-stationarity (heteroscedasticity): Variance changes over time → fix with log or Box-Cox transformation

Differencing — The Mathematical View

First difference (removes linear trend):

yₜ' = yₜ − yₜ₋₁ (d=1 in ARIMA)

Seasonal difference (removes seasonal pattern):

yₜ' = yₜ − yₜ₋ₘ (D=1 in ARIMA, where m = seasonal period)

Second difference (removes quadratic trend — rarely needed):

yₜ'' = yₜ' − yₜ'₋₁ (d=2 in ARIMA)

Rule of thumb: Use the minimum number of differences needed. Over-differencing can introduce spurious structure.

3. Inference & Conclusions

Google stock price (goog200) is non-stationary — one regular difference is sufficient to make it stationary.
a10 drug sales require: log transformation (for variance) + differencing (for trend).
usmelec requires: log transformation + seasonal differencing.
The ndiffs() and nsdiffs() functions provide algorithmic guidance for the number of differences — this feeds directly into
ARIMA model specification (the d and D parameters).
The 9-panel comparison plots in par(mfrow=c(3,3)) give a rich visual overview of how different economic and natural
phenomena behave over time.
 Lab 8 — KPSS Test & ARIMA Models
Overview
This lab introduces a formal statistical test for stationarity (the KPSS test) and the powerful ARIMA model family — the most widely used
class of models for univariate time series forecasting. ARIMA models are built on three components: AutoRegression (AR), Integration (I =
differencing), and Moving Average (MA).

Full Code

library(fpp2)
library(urca)

autoplot(goog)
goog %>% [Link]() %>% summary()
goog %>% diff() %>% [Link]() %>% summary()

autoplot(uschange[,"Consumption"])
fit <- [Link](uschange[,"Consumption"])
fit
fit2 <- [Link](uschange[,"Consumption"], seasonal=FALSE)
fit2
fit2 %>% forecast(h=10) %>% autoplot()
fit %>% forecast(h=10) -> y
y

ggAcf(uschange[,"Consumption"])
ggPacf(uschange[,"Consumption"])
fit3 <- Arima(uschange[,"Consumption"], order=c(3,0,0))
fit3
checkresiduals(fit3)
autoplot(fit3)

1. Code Walkthrough
KPSS Test for Stationarity

library(urca)
goog %>% [Link]() %>% summary()

urca — Unit Root and Cointegration Analysis package.

[Link]() — Performs the KPSS Test (Kwiatkowski-Phillips-Schmidt-Shin Test).

KPSS Test Hypotheses (OPPOSITE to ADF):

H₀ (Null): The series IS stationary (trend-stationary)

H₁ (Alternative): The series is NOT stationary (has a unit root)
 Decision rule:

Test statistic > Critical value → Reject H₀ → Series is NON-STATIONARY

Test statistic < Critical value → Fail to reject H₀ → Series IS STATIONARY

Comment from code: "test statistic value is much much higher than the critical values so we will reject the null hypothesis and the
dataset is non-stationary"

goog %>% diff() %>% [Link]() %>% summary()

Comment from code: "test statistic has become lower than critical values hence not to reject the null hypothesis and so it is now
trend stationary" → After one difference, goog becomes stationary.

Comparison of Stationarity Tests:

Test H₀ H₁ Reject H₀ means

KPSS Stationary Non-stationary Series is non-stationary
ADF (Augmented Dickey-Fuller) Non-stationary (unit root) Stationary Series is stationary
PP (Phillips-Perron) Non-stationary Stationary Series is stationary

Best practice: Use multiple tests together for confirmation. KPSS and ADF have opposite null hypotheses, so you want: KPSS fails
to reject (stationary) AND ADF rejects (stationary) simultaneously.

ARIMA Background

ARIMA stands for AutoRegressive Integrated Moving Average.

An ARIMA(p, d, q) model has three parameters:

p = order of the AR (AutoRegressive) part

d = degree of differencing (Integration)
q = order of the MA (Moving Average) part

AR(p) — AutoRegressive Model:

yₜ = c + φ₁yₜ₋₁ + φ₂yₜ₋₂ + ... + φₚyₜ₋ₚ + εₜ

The current value depends on its own p previous values.

φ₁, φ₂, ..., φₚ = AR coefficients
εₜ = white noise error
Real-world analogy: Today's temperature depends on yesterday's, two days ago's, etc.

MA(q) — Moving Average Model:

yₜ = c + εₜ + θ₁εₜ₋₁ + θ₂εₜ₋₂ + ... + θqεₜ₋q

 The current value depends on q previous forecast errors (residuals).
θ₁, θ₂, ..., θq = MA coefficients
Analogy: Today's sales are influenced not just by past sales, but by past "surprises" (unexpected demand spikes from q periods
ago).

Integrated (I(d)): The data has been differenced d times to achieve stationarity.

Combined ARIMA(p,d,q):

Φ(B) × (1-B)^d × yₜ = c + Θ(B) × εₜ

(Using backshift operator notation — not required for MBA level)

Automatic ARIMA

fit <- [Link](uschange[,"Consumption"])

fit

[Link]() — Automatically selects the best ARIMA model by:

1. Determining d using KPSS tests (how many differences needed)
2. Searching over combinations of p and q
3. Selecting the model with the lowest AICc (corrected AIC)
uschange[,"Consumption"] — Quarterly percentage change in US consumption.

AIC (Akaike Information Criterion):

AIC = −2 × log(Likelihood) + 2k

Where k = number of parameters. Lower AIC = better model. AIC penalizes complexity.

AICc = AIC corrected for small sample sizes (preferred over AIC).

BIC (Bayesian Information Criterion):

BIC = −2 × log(Likelihood) + k × log(n)

BIC penalizes complexity more heavily than AIC (especially for large n).

fit2 <- [Link](uschange[,"Consumption"], seasonal=FALSE)

fit2

seasonal=FALSE — Restricts the search to non-seasonal ARIMA models only (ARIMA(p,d,q) without seasonal components).

fit2 %>% forecast(h=10) %>% autoplot()

forecast(h=10) — Generates forecasts 10 periods ahead.

autoplot() — Plots the forecast with prediction intervals.

Manual ARIMA Using ACF and PACF

 ggAcf(uschange[,"Consumption"])
ggPacf(uschange[,"Consumption"])

Reading ACF and PACF for Model Identification:

The ACF and PACF plots are the key diagnostic tools for identifying p and q.

Pattern Model suggested

ACF cuts off after lag q, PACF tails off MA(q)
PACF cuts off after lag p, ACF tails off AR(p)
Both tail off gradually ARMA(p,q) — mixed model
ACF shows slow decay (large significant values at many lags) Non-stationary → difference first
Spike at seasonal lag in ACF Seasonal component → use SARIMA

"Cuts off" = drops abruptly to near zero after lag k (inside the confidence bounds)

"Tails off" = gradually decreases toward zero over many lags

ggPacf() — Partial AutoCorrelation Function (PACF): measures the correlation between yₜ and yₜ₋ₖ after removing the effect of
all intermediate lags (1, 2, ..., k−1). The PACF directly reveals the order of the AR component.

PACF Formula (conceptually):

PACF at lag k = correlation between yₜ and yₜ₋ₖ | yₜ₋₁, ..., yₜ₋(k-1)

fit3 <- Arima(uschange[,"Consumption"], order=c(3,0,0))

fit3

Arima() (capital A) — Fits a manually specified ARIMA model.

order=c(3,0,0) — Specifies ARIMA(3,0,0) = AR(3):
p=3 → Uses 3 lagged values as predictors
d=0 → No differencing needed (Consumption is already stationary as % changes)
q=0 → No MA component

Comment from code: "AIC and BIC are coming lower than [Link] method, hence manual method with ACF and PACF plot is
better performing"

This shows that manual inspection of ACF/PACF can sometimes outperform the automated selection.

Residual Diagnostics for ARIMA

checkresiduals(fit3)

Checks: (1) residual time plot, (2) ACF of residuals, (3) histogram.
p-value > 0.05 from the Ljung-Box test → fail to reject H₀ of no autocorrelation → residuals are white noise → ARIMA model is well-
specified.
 Comment from code: "p-value > alpha which is 0.05 hence failed to reject the null hypothesis, hence exist no correlation"

Inverse Characteristic Root Plot

autoplot(fit3)

Plots the inverse characteristic roots (also called inverse AR roots and MA roots) on the complex unit circle.
AR roots (φ): Must lie INSIDE the unit circle (|root| < 1) for the model to be stationary.
MA roots (θ): Must lie INSIDE the unit circle for the model to be invertible.

Comment from code: "phi1, phi2, phi3 are lying within −1 to +1" — all three AR coefficients are within bounds → the model is
stationary and valid.

What the Root Plot Looks Like

A circle of radius 1 is drawn.

Dots represent the roots (inverse roots actually).
All dots inside the circle → model is valid (stationary + invertible).
A dot on or outside the circle → model is problematic (unit root or explosive behaviour).

3. Inference & Conclusions

KPSS test confirms goog is non-stationary; one difference makes it stationary.
[Link]() is a reliable starting point, but manual inspection of ACF/PACF can sometimes yield a better model (lower AIC/BIC).
For uschange[,"Consumption"]: ARIMA(3,0,0) was selected manually and performed better — this is a pure AR(3) model
(today's consumption depends on the last 3 quarters).
Residual diagnostics are non-negotiable — always check that residuals are white noise.
Inverse root plot confirms model validity — all roots inside the unit circle.

Lab 9 — Seasonal ARIMA (SARIMA) Models

Overview
This lab extends ARIMA modelling to seasonal data using the SARIMA model (Seasonal ARIMA). SARIMA adds seasonal AR and MA
components to the regular ARIMA framework, making it one of the most comprehensive and flexible models for seasonal time series.

Full Code
 library(fpp2)
autoplot(euretail) + ylab("Retail index") + xlab("Year")
nsdiffs(euretail)
euretail %>% ggtsdisplay()
euretail %>% diff(lag=4) %>% ggtsdisplay()
euretail %>% diff(lag=4) %>% ndiffs()
euretail %>% diff(lag=4) %>% diff() %>% ggtsdisplay()

fit1 <- euretail %>% Arima(order=c(0,1,1), seasonal=c(0,1,1))

fit1
fit2 <- euretail %>% Arima(order=c(1,1,0), seasonal=c(1,1,0))
fit2
fit3 <- euretail %>% Arima(order=c(0,1,1), seasonal=c(1,1,0))
fit3
fit4 <- euretail %>% Arima(order=c(1,1,0), seasonal=c(0,1,1))
fit4
fit5 <- euretail %>% Arima(order=c(1,1,1), seasonal=c(1,1,1))
fit5

euretail %>% Arima(order=c(0,1,1), seasonal=c(0,1,1)) %>% residuals() %>% ggtsdisplay()

euretail %>% Arima(order=c(0,1,2), seasonal=c(0,1,1)) %>% residuals() %>% ggtsdisplay()
euretail %>% Arima(order=c(0,1,3), seasonal=c(0,1,1)) %>% residuals() %>% ggtsdisplay()
fit6 <- euretail %>% Arima(order=c(0,1,3), seasonal=c(0,1,1))
fit6
[Link](euretail)
fit6 %>% forecast(h=12) %>% autoplot()

1. Code Walkthrough
Dataset and Initial Visualization

autoplot(euretail) + ylab("Retail index") + xlab("Year")

nsdiffs(euretail)

euretail — Quarterly retail trade index for the Euro area (1996–2011).
nsdiffs(euretail) → Returns 1 → needs 1 seasonal difference (D=1, seasonal period m=4 for quarterly).

SARIMA Notation

Full SARIMA model notation: ARIMA(p,d,q)(P,D,Q)[m]

Parameter Full Name Meaning

p Non-seasonal AR order Number of lagged regular values used
d Non-seasonal differencing How many times to regularly difference
q Non-seasonal MA order Number of lagged regular errors used
P Seasonal AR order Number of lagged seasonal values used
D Seasonal differencing How many seasonal differences needed
Q Seasonal MA order Number of lagged seasonal errors used
m Seasonal period 4 for quarterly, 12 for monthly
 SARIMA(p,d,q)(P,D,Q)[m] General Formula:

Φ_P(Bᵐ) × φ_p(B) × (1−Bᵐ)^D × (1−B)^d × yₜ = c + Θ_Q(Bᵐ) × θ_q(B) × εₜ

Where B = backshift operator (Byₜ = yₜ₋₁).

Building the SARIMA Model Step-by-Step

Step 1: Assess the original series

euretail %>% ggtsdisplay()

ggtsdisplay() — Produces 3 panels simultaneously: (1) time plot, (2) ACF, (3) PACF. Essential for ARIMA identification.
Original euretail shows trend and seasonal pattern → needs differencing.

Step 2: Apply seasonal differencing

euretail %>% diff(lag=4) %>% ggtsdisplay()

diff(lag=4) — Applies a seasonal difference with lag=4 (for quarterly data).

Formula: yₜ' = yₜ − yₜ₋₄ (current quarter minus same quarter last year)
Removes the seasonal pattern. After this, check ACF/PACF for remaining autocorrelation.

Step 3: Check if regular differencing is also needed

euretail %>% diff(lag=4) %>% ndiffs()

Returns 1 → yes, one regular difference is also needed.

Step 4: Apply both differences and analyze

euretail %>% diff(lag=4) %>% diff() %>% ggtsdisplay()

diff(lag=4) %>% diff() — First seasonal difference, THEN regular difference.

This gives us: d=1, D=1 → the "I" part of SARIMA is determined.
Now read ACF and PACF of the doubly-differenced series to choose p, q, P, Q.

Comment from code: "p=1, d=1, q=1; P=1, D=1, Q=1" — the ACF/PACF suggest starting candidates.

Fitting Multiple SARIMA Models

fit1 <- euretail %>% Arima(order=c(0,1,1), seasonal=c(0,1,1))

order=c(p,d,q) — Non-seasonal ARIMA parameters: p=0, d=1, q=1

seasonal=c(P,D,Q) — Seasonal parameters: P=0, D=1, Q=1
This is SARIMA(0,1,1)(0,1,1)[4] — a common starting model for quarterly data.
fit2 <- euretail %>% Arima(order=c(1,1,0), seasonal=c(1,1,0))

SARIMA(1,1,0)(1,1,0)[4] — AR model instead of MA for both non-seasonal and seasonal parts.

Multiple models (fit1 through fit5) are compared by their AIC/BIC values printed by the output.

Residual Analysis to Refine the Model

euretail %>% Arima(order=c(0,1,1), seasonal=c(0,1,1)) %>% residuals() %>% ggtsdisplay()

euretail %>% Arima(order=c(0,1,2), seasonal=c(0,1,1)) %>% residuals() %>% ggtsdisplay()
euretail %>% Arima(order=c(0,1,3), seasonal=c(0,1,1)) %>% residuals() %>% ggtsdisplay()

After fitting SARIMA(0,1,1)(0,1,1), plot the residuals' ACF/PACF to check for remaining autocorrelation.
If ACF/PACF of residuals still shows significant spikes → the model hasn't captured all the structure → increase p or q.
Progressively increasing q from 1 → 2 → 3 until all residual spikes disappear.

fit6 <- euretail %>% Arima(order=c(0,1,3), seasonal=c(0,1,1))

fit6

SARIMA(0,1,3)(0,1,1)[4] — The final selected model: no AR, 3 MA terms, 1 regular difference + 1 seasonal difference with 1
seasonal MA term.

Automatic SARIMA

[Link](euretail)

The automatic algorithm may select the same or a similar model.

Comparison: If the manually identified model (fit6) has lower AIC than [Link], the manual approach wins.

Forecast

fit6 %>% forecast(h=12) %>% autoplot()

h=12 — Forecast 12 quarters ahead = 3 years (since m=4 quarters per year).
Comment: "next 3 years data h=12 means 4 quarters each year, h=horizon"

What the Forecast Plot Looks Like

Historical data with clear seasonal pattern and trend.

Forecasts continue the seasonal oscillation with an upward trend.
Prediction intervals widen over the 3-year horizon.
The seasonal amplitude in the forecast mirrors the historical seasonal pattern.

3. Inference & Conclusions

euretail requires: D=1 (seasonal difference) + d=1 (regular difference) → total integration = I(1,1).
 Manual identification process: Start from ACF/PACF of doubly-differenced data → try multiple candidate models → check residuals
→ refine.
SARIMA(0,1,3)(0,1,1)[4] (fit6) is selected as the best model.
The quarterly forecast for the next 3 years shows continuing seasonal patterns and trend.
[Link] serves as a validation check — agreement between manual and automatic selection builds confidence.

Lab 9 Part 2 — Stock Data & Financial Returns

Overview
This supplementary lab introduces tools for fetching and visualizing real-time financial data using the quantmod and
PerformanceAnalytics packages. It sets up the foundation for GARCH modelling (Lab 10) by showing how to compute financial returns.

Full Code

[Link]("quantmod")
[Link]("rugarch")
[Link]("xts")
[Link]("PerformanceAnalytics")
library(quantmod)
library(rugarch)
library(xts)
library(PerformanceAnalytics)

df <- getSymbols("TSLA", from="2010-01-01", to="2024-12-31")

head(TSLA)
chartSeries(TSLA)
chartSeries(TSLA["2020-12"])
return <- CalculateReturns(TSLA$[Link])
return
return <- return[-c(1),]
chart_Series(return)

1. Code Walkthrough
Package Overview
Package Purpose
quantmod Download financial data (stocks, FX, crypto) from Yahoo Finance
rugarch Fit GARCH-family models for volatility
xts eXtensible Time Series — a powerful time series format for financial data
PerformanceAnalytics Portfolio performance and risk metrics

Downloading Tesla Stock Data

df <- getSymbols("TSLA", from="2010-01-01", to="2024-12-31")

 getSymbols() — Downloads daily OHLCV (Open, High, Low, Close, Volume) + Adjusted price data from Yahoo Finance.
"TSLA" — The stock ticker symbol for Tesla Inc.
Important behaviour: getSymbols() does NOT return data to df directly. Instead, it creates an object named TSLA in the global
environment (the variable df just stores the string "TSLA").

head(TSLA)

Shows the first 6 rows of the TSLA xts object. Columns:

[Link] — Opening price of the day
[Link] — Highest price of the day
[Link] — Lowest price of the day
[Link] — Closing price of the day
[Link] — Number of shares traded
[Link] — Closing price adjusted for dividends and stock splits (most important for analysis)

Visualizing Stock Price

chartSeries(TSLA)
chartSeries(TSLA["2020-12"])

chartSeries() — Creates a financial chart including:

Candlestick chart (OHLC prices)
Volume bars at the bottom
Optional technical indicators
TSLA["2020-12"] — Uses xts time-based subsetting: extracts only December 2020 data. This is a feature of the xts format —
you can use date strings directly.

What the Chart Looks Like

Top panel: Candlestick chart. Each bar represents one day. Green = price went up; Red = price went down. The top and bottom of
the rectangle show open and close; the lines (wicks) show high and low.
Bottom panel: Volume bars for each day.
December 2020 was a period of extreme volatility for Tesla (it was being added to the S&P 500).

Computing Financial Returns

return <- CalculateReturns(TSLA$[Link])

return <- return[-c(1),]

TSLA$[Link] — Extracts only the Adjusted Close price column.

CalculateReturns() — Computes simple (or log) returns from price data.

Simple Return Formula:

Rₜ = (Pₜ − Pₜ₋₁) / Pₜ₋₁ = Pₜ/Pₜ₋₁ − 1

Log Return Formula:

 rₜ = log(Pₜ/Pₜ₋₁) = log(Pₜ) − log(Pₜ₋₁)

return[-c(1),] — Removes the first row (which is NA because there is no previous price to compute the first return).
-c(1) means "exclude row 1". The result is returns from the second trading day onwards.

Why returns instead of prices?

Stock prices are non-stationary (random walk).

Stock returns are approximately stationary (though they have fat tails and volatility clustering — addressed by GARCH).
Returns are comparable across different stocks regardless of price level.

chart_Series(return)

Plots the return series.

What the Returns Plot Looks Like

X-axis: Time (2010–2024)

Y-axis: Daily return (as a fraction or percentage)
Pattern: Returns fluctuate randomly around 0 (mean ≈ 0), BUT the spread (volatility) is not constant — some periods have very
large swings (high volatility), others are calm (low volatility).
This volatility clustering phenomenon is the motivation for GARCH models (Lab 10).

3. Inference & Conclusions

Tesla's adjusted price shows a massive upward trend with explosive growth especially 2019–2021.
The return series is approximately stationary (mean ~0) but shows clear volatility clustering — periods of high turbulence followed
by calm.
This heteroscedasticity (non-constant variance) in returns violates the ARIMA assumptions → GARCH models are needed.

Lab 10 — GARCH Models for Volatility Modelling

Overview
GARCH (Generalized AutoRegressive Conditional Heteroscedasticity) models are the gold standard for modelling time-varying
volatility in financial time series. While ARIMA models the mean (level) of a series, GARCH models the variance (volatility). This is critical
for risk management, option pricing, and portfolio optimization.

Full Code
library(dplyr)
library(tidyverse)
library(tseries)
library(rugarch)
library(xts)
library(PerformanceAnalytics)
library(quantmod)

df <- getSymbols("TSLA", from="2010-01-01", to="2024-12-31")

return <- CalculateReturns(TSLA$[Link])
return <- return[-c(1),]

# Model 1: sGARCH(1,1) with Normal distribution

mod_specify <- ugarchspec(
[Link] = list(armaorder=c(0,0)),
[Link] = list(model="sGARCH", garchorder=c(1,1)),
[Link]='norm'
)
mod_fitting <- ugarchfit(data=return, spec=mod_specify, [Link]=20)
mod_fitting
plot(mod_fitting, which='all')

# Model 2: sGARCH(1,1) with Skewed t-distribution

mod_specify2 <- ugarchspec(
[Link] = list(armaorder=c(0,0)),
[Link] = list(model="sGARCH", garchorder=c(1,1)),
[Link]='sstd'
)
mod_fitting2 <- ugarchfit(data=return, spec=mod_specify2, [Link]=20)
mod_fitting2
plot(mod_fitting2, which='all')

# Model 3: GJR-GARCH(1,1) with Skewed t-distribution

mod_specify3 <- ugarchspec(
[Link] = list(armaorder=c(0,0)),
[Link] = list(model="gjrGARCH", garchorder=c(1,1)),
[Link]='sstd'
)
mod_fitting3 <- ugarchfit(data=return, spec=mod_specify3, [Link]=20)
mod_fitting3
plot(mod_fitting3, which='all')

# Forecasting
fore <- ugarchforecast(fitORspec=mod_fitting2, [Link]=20)
fore
plot(fitted(fore))
plot(sigma(fore))

1. Code Walkthrough
 Understanding Conditional Heteroscedasticity

Before diving into GARCH, understand the problem it solves:

Heteroscedasticity = non-constant variance over time. Conditional Heteroscedasticity = the variance at time t depends (is conditional on)
the variance and shocks from previous periods.

Why does this matter for financial data?

After a large market shock (crash), volatility stays high for several days/weeks.
During quiet periods, volatility is low and stays low.
This is volatility clustering — large changes are followed by large changes (of either sign).

ARCH effect: If the squared residuals (εₜ²) are autocorrelated, there is an ARCH effect — and GARCH is the right tool.

Specifying a GARCH Model

mod_specify <- ugarchspec(

[Link] = list(armaorder=c(0,0)),
[Link] = list(model="sGARCH", garchorder=c(1,1)),
[Link]='norm'
)

ugarchspec() — UGARCH Specification: defines the structure of the GARCH model (from the rugarch package). Think of it as
writing a blueprint before building.

Three components to specify:

1. [Link] = list(armaorder=c(0,0))

Specifies the ARMA model for the mean (level) of the return series.
armaorder=c(p,q) → ARMA(0,0) = no AR, no MA → the mean is just a constant (μ).
This means: we assume returns have a constant mean (approximately zero for financial returns).
Could also be ARMA(1,0) or ARMA(1,1) if mean has autocorrelation structure.

2. [Link] = list(model="sGARCH", garchorder=c(1,1))

model="sGARCH" — Standard GARCH model.

garchorder=c(p,q) — GARCH(p,q):
p=1 (q in GARCH notation) = order of the ARCH term (lagged squared residuals)
q=1 (p in GARCH notation) = order of the GARCH term (lagged conditional variance)

GARCH(1,1) Equations:

Mean equation:

rₜ = μ + εₜ, εₜ = σₜ × zₜ, zₜ ~ iid(0,1)

Variance equation:

σₜ² = ω + α × εₜ₋₁² + β × σₜ₋₁²

Where:
 σₜ² = Conditional variance (today's volatility, estimated by the model)
ω (omega) = Long-run average variance (baseline variance; must be > 0)
α (alpha) = ARCH coefficient: how strongly yesterday's shock (εₜ₋₁²) affects today's volatility
β (beta) = GARCH coefficient: how strongly yesterday's volatility (σₜ₋₁²) persists into today
εₜ₋₁² = Squared residual from yesterday (the "shock")
σₜ₋₁² = Estimated variance from yesterday

Key constraint: α + β < 1 → ensures the variance is stationary (returns to a long-run mean).

Analogy: Imagine σₜ ² is your daily stress level:

ω = your baseline resting stress

α = how much yesterday's bad event (shock) raises today's stress
β = how much yesterday's stress level carries forward to today
If α + β is close to 1, stress takes a long time to dissipate (high persistence).

3. [Link]='norm'

'norm' — Assumes the standardized residuals (zₜ) follow a Normal (Gaussian) distribution.
Financial returns often have fat tails (extreme events happen more frequently than Normal predicts), so 'norm' may underestimate
risk.

Distribution options:

Code Distribution Notes

'norm' Normal (Gaussian) Standard; thin tails
'std' Student's t Fat tails; symmetric
'sstd' Skewed Student's t Fat tails + asymmetric (left-skewed returns)
'ged' Generalized Error Distribution Flexible tail behaviour
'snorm' Skewed Normal Asymmetric; thin tails

Fitting the GARCH Model

mod_fitting <- ugarchfit(data=return, spec=mod_specify, [Link]=20)

mod_fitting

ugarchfit() — Fits the specified GARCH model to the data using Maximum Likelihood Estimation (MLE).
data=return — The Tesla daily return series.
[Link]=20 — Keeps the last 20 observations as an out-of-sample test set (not used for fitting, used for evaluation).

Output interpretation:

The mod_fitting printout shows:

Optimal Parameters:
mu (μ) = estimated constant mean return
omega (ω) = baseline variance
alpha1 (α) = ARCH coefficient
beta1 (β) = GARCH coefficient
α + β = persistence: if close to 1.0, volatility is highly persistent
Information Criteria: AIC, BIC (lower = better model)
 Ljung-Box Tests: On residuals and squared residuals — both should have p > 0.05 for a good model.
Sign Bias Test: Tests whether positive and negative shocks have different effects on volatility (asymmetry).

plot(mod_fitting, which='all')

which='all' — Produces ALL available diagnostic plots (typically 12 plots for rugarch).

Key Diagnostic Plots for GARCH

1. Series with 2 Conditional SD Superimposed: Original return series with the estimated volatility bands (±2σₜ) overlaid. Good model
should have most observations within bands.
2. Conditional Standard Deviation: The time-varying σₜ estimated by the model. Should spike during market crises (2020 COVID
crash, 2022 bear market).
3. Standardized Residuals: rₜ/σₜ — should look like white noise.
4. ACF of Standardized Residuals: Should show no significant spikes.
5. ACF of Squared Standardized Residuals: Should show no significant spikes (confirms ARCH effects are captured).
6. QQ Plot of Standardized Residuals: Points should fall along the diagonal. Deviations at the tails indicate fat-tail behaviour
(suggests using 'sstd' instead of 'norm').

Model 2: Skewed t-Distribution

mod_specify2 <- ugarchspec(

[Link] = list(armaorder=c(0,0)),
[Link] = list(model="sGARCH", garchorder=c(1,1)),
[Link]='sstd'
)
mod_fitting2 <- ugarchfit(data=return, spec=mod_specify2, [Link]=20)

Same GARCH(1,1) structure but now with 'sstd' (Skewed Student's t-distribution).
This adds two more parameters: degrees of freedom (ν, nu — controls tail thickness) and skewness (ξ, xi).
For Tesla returns, negative skewness is expected (crashes are sharper than rallies).
Better model if: AIC/BIC lower AND QQ plot improves AND Ljung-Box tests still pass.

Model 3: GJR-GARCH (Asymmetric GARCH)

mod_specify3 <- ugarchspec(

[Link] = list(armaorder=c(0,0)),
[Link] = list(model="gjrGARCH", garchorder=c(1,1)),
[Link]='sstd'
)

model="gjrGARCH" — GJR-GARCH (Glosten-Jagannathan-Runkle GARCH), also known as Threshold GARCH (TGARCH).

GJR-GARCH Variance Equation:

σₜ² = ω + α × εₜ₋₁² + γ × εₜ₋₁² × I(εₜ₋₁ < 0) + β × σₜ₋₁²

Where:
 I(εₜ₋₁ < 0) = Indicator function: equals 1 if yesterday's return was negative (bad news), 0 otherwise.
γ (gamma) = Asymmetry/leverage coefficient: If γ > 0, negative shocks (losses) increase volatility MORE than positive shocks of
the same magnitude.

Why this matters — The Leverage Effect:

In stock markets, negative returns (price drops) tend to increase volatility MORE than positive returns (price rises) of the same size.
This is called the leverage effect — companies have more debt (leverage) when prices fall, making them riskier.
Standard GARCH assumes symmetric responses; GJR-GARCH captures the asymmetry.

Model Comparison:

Model Distribution Captures

sGARCH + norm Normal Basic volatility clustering
sGARCH + sstd Skewed t Clustering + fat tails + skewness
gjrGARCH + sstd Skewed t Clustering + fat tails + skewness + leverage effect

Forecasting with GARCH

fore <- ugarchforecast(fitORspec=mod_fitting2, [Link]=20)

fore
plot(fitted(fore))
plot(sigma(fore))

ugarchforecast() — Generates forecasts from a fitted GARCH model.

fitORspec=mod_fitting2 — Uses the fitted Model 2 as the base.
[Link]=20 — Forecasts 20 periods (trading days) ahead.

Output:

fitted(fore) — Forecasted mean returns (μ) for the next 20 days. For ARMA(0,0), this is a constant flat line.
sigma(fore) — Forecasted conditional standard deviation (σₜ) for the next 20 days.

What the Forecasted Sigma Plot Looks Like

X-axis: Forecast horizon (1 to 20 days ahead)

Y-axis: Forecasted volatility (σ)
Pattern: The volatility forecast starts at the current level and gradually converges toward the long-run unconditional volatility
level (ω / (1−α−β)).
If current volatility is high, the forecast will start high and decay.
If current volatility is low, it will start low and slightly converge upward.
This mean-reversion of volatility is a key feature of GARCH.

3. Inference & Conclusions

Tesla daily returns show clear ARCH effects (volatility clustering) — GARCH is appropriate.
GARCH(1,1) with 'norm' is the baseline; check if α+β is close to 1 (high persistence = volatility lasts long).
GARCH(1,1) with 'sstd' should improve the QQ plot fit (better tail modelling) with lower AIC/BIC.
GJR-GARCH with 'sstd' captures the leverage effect — if γ > 0 and statistically significant, negative returns cause higher volatility
than positive ones (this is almost always true for equities).
The volatility forecast reverts to the long-run average — useful for options pricing and Value-at-Risk (VaR) calculation over a 20-day
horizon.
 Model selection: Compare AIC, BIC, residual diagnostics, and QQ plots across all three models. The best model (likely GJR-
GARCH with sstd) should be used for risk management.

Comprehensive Glossary
All technical terms, abbreviations, notations, and symbols used across all 10 labs, in alphabetical order.

Term / Symbol Full Form / Definition

A dataset in fpp2: monthly Australian antidiabetic drug sales ($
a10
million, 1992–2008)
AutoCorrelation Function — measures the correlation between
ACF a time series and its own lagged values at various lags k. Plot
shows spikes at each lag.
Augmented Dickey-Fuller Test — tests for a unit root (non-
ADF Test
stationarity). H₀: Non-stationary. Reject H₀ → stationary.
Akaike Information Criterion — model selection criterion: AIC =
AIC −2×log(L) + 2k. Lower AIC = better model balancing fit and
parsimony.
Corrected AIC — AIC adjusted for small sample size. Preferred
AICc
over AIC in small samples.
AutoRegressive Conditional Heteroscedasticity — a model for
ARCH time-varying variance where current variance depends on past
squared errors.
AutoRegressive — a model where the current value depends
AR
on its own lagged values: yₜ = φ₁yₜ₋₁ + ... + φₚyₜ₋ₚ + εₜ
AutoRegressive Integrated Moving Average — a class of
ARIMA
models for non-seasonal time series. ARIMA(p,d,q).
AutoRegressive Moving Average — combines AR(p) and MA(q)
ARMA
without differencing.
Dataset in fpp2: quarterly Australian beer production
ausbeer
(megalitres, 1956–2010)
ausair Dataset in fpp2: annual Australian air passengers (millions)
Dataset in fpp2: quarterly international tourist visitor nights in
austourists
Australia
Function from fpp2/ggplot2 that automatically creates
autoplot()
appropriate time series plots
Adds additional layers (forecast lines, fitted values) to an
autolayer()
existing autoplot
B Backshift Operator — Byₜ = yₜ₋₁. Powers: Bᵏyₜ = yₜ₋ₖ. Used in
ARIMA algebra.
1. Trend smoothing parameter in Holt/Holt-Winters (0 < β < 1).
β (beta) 2. GARCH coefficient (lagged variance). 3. Regression slope
coefficient.
Bayesian Information Criterion — similar to AIC but penalizes
BIC
complexity more: BIC = −2log(L) + k×log(n). Lower = better.
A family of power transformations to stabilize variance: wₜ =
Box-Cox
(yₜ^λ−1)/λ if λ≠0; log(yₜ) if λ=0
Box-Pierce Test A test for autocorrelation in residuals. H₀: No autocorrelation.
Less accurate than Ljung-Box for small samples.
Column bind — combines multiple time series side by side
cbind()
when they have the same time range but different Y-axis scales
Function from fpp2 that produces residual time plot, ACF,
checkresiduals()
histogram, and runs the Ljung-Box test
 Term / Symbol Full Form / Definition
Classical Simplest decomposition: uses centred moving averages for
Decomposition trend, then averages by season. Assumes constant seasonality.
Non-seasonal differencing order in ARIMA(p,d,q). Number of
d
regular differences applied.
Seasonal differencing order in SARIMA(p,d,q)(P,D,Q)[m].
D
Number of seasonal differences applied.
A modification to Holt's method where the trend gradually
damped flattens toward a constant (controlled by φ). Prevents over-
forecasting.
Base R function for classical additive or multiplicative time
decompose()
series decomposition into Trend + Seasonal + Random
Takes differences of a time series. diff(y) = yₜ−yₜ₋₁. diff(y,
diff()
lag=m) = yₜ−yₜ₋ₘ (seasonal difference).
A forecasting method (Drift method / rwf with drift=TRUE) that
drift
extrapolates the average historical trend into the future
Error, Trend, Seasonality — a unified framework for exponential
ETS smoothing models. Each component can be None (N), Additive
(A), Multiplicative (M), or Additive Damped (Ad).
Dataset in fpp2: monthly new orders for electrical equipment in
elecequip
the Euro area (1996–2012)
Dataset in fpp2: half-hourly electricity demand and temperature
elecdemand
in Victoria, Australia (2014)
Dataset in fpp2: annual electricity sales in South Australia
elecsales
(1989–1998)
White noise error term at time t — random, mean-zero,
εₜ (epsilon)
constant variance, uncorrelated across time
Dataset in fpp2: quarterly retail trade index for the Euro area
euretail
(1996–2011)
Function in fpp2 that automatically selects the best ETS model
ets()
using AIC
1. Seasonal smoothing parameter in Holt-Winters (0 < γ < 1). 2.
γ (gamma)
GJR-GARCH asymmetry coefficient.
Generalized AutoRegressive Conditional Heteroscedasticity —
GARCH models time-varying conditional variance: σₜ² = ω + αεₜ₋₁² +
βσₜ₋₁²
In ugarchspec(): c(p,q) for GARCH order — p = ARCH terms
garchorder
(lagged ε²), q = GARCH terms (lagged σ²)
Function from quantmod that downloads historical financial data
getSymbols()
from Yahoo Finance
R package extending ggplot2 with functions like ggpairs() for
GGally
pair plots
Creates a scatterplot matrix showing all pairwise scatter plots,
ggpairs()
distributions, and correlations
ggAcf() ggplot2-style ACF plot from fpp2
ggPacf() ggplot2-style PACF plot from fpp2
Creates a seasonal plot where each year is a separate coloured
ggseasonplot()
line, months on X-axis
Creates a seasonal subseries plot: separate mini-panels for
ggsubseriesplot()
each season across years
Displays time plot + ACF + PACF in one panel — essential for
ggtsdisplay()
ARIMA identification
Glosten-Jagannathan-Runkle GARCH — asymmetric GARCH
gjrGARCH
model: σₜ² = ω + αεₜ₋₁² + γεₜ₋₁²×I(εₜ₋₁<0) + βσₜ₋₁²
goog200 Dataset in fpp2: first 200 daily closing prices of Google stock
 Term / Symbol Full Form / Definition
Null Hypothesis — the default assumption being tested.
H₀
Rejected if p-value < significance level (α = 0.05).
H₁ Alternative Hypothesis — the conclusion if H₀ is rejected.
Forecast horizon — the number of periods ahead to forecast.
h
h=12 = 12 periods ahead.
Function in fpp2 implementing Holt's Linear (Double
holt()
Exponential Smoothing) method for trending data
Double Exponential Smoothing — extends SES with a trend
Holt's Method
component. Two smoothing parameters: α and β.
Triple Exponential Smoothing — extends Holt's with a seasonal
Holt-Winters
component. Three parameters: α, β, γ.
Function in fpp2 for Holt-Winters' method (additive or
hw()
multiplicative seasonal)
Non-constant variance in a time series or regression residuals.
heteroscedasticity
The opposite of homoscedasticity (constant variance).
Integrated of order d — means the series needs to be
I(d)
differenced d times to become stationary.
Independent and Identically Distributed — a sequence of
iid random variables with no autocorrelation, all from the same
distribution.
Box-Cox transformation parameter. λ=0 → log transformation;
λ (lambda)
λ=0.5 → square root; λ=1 → no transformation.
The number of time periods between two observations being
Lag
compared. Lag k means comparing yₜ with yₜ₋ₖ.
In finance: the tendency for stock price declines to increase
Leverage Effect volatility more than equivalent price increases. Captured by
GJR-GARCH (γ > 0).
livestock Dataset in fpp2: annual sheep population in Asia (millions)
An improved version of Box-Pierce test for autocorrelation in
Ljung-Box Test residuals. H₀: No autocorrelation. p > 0.05 → model is
acceptable.
lₜ Level component at time t in exponential smoothing models
Log Return rₜ = log(Pₜ/Pₜ₋₁) — continuously compounded return. Preferred
in econometrics; additively compoundable.
Locally Estimated Scatterplot Smoothing — a non-parametric
LOESS
regression technique used in STL decomposition
Seasonal period — number of observations per year. m=4
m
(quarterly), m=12 (monthly), m=52 (weekly).
Moving Average — in ARIMA context: current value depends
MA
on past errors. yₜ = εₜ + θ₁εₜ₋₁ + ... + θqεₜ₋q
Function in fpp2 computing a simple or centred moving average
ma()
of specified order
. Easy to
interpret;
MAE Mean Absolute Error — (1/n)Σ eₜ
same units as
data.
×100%.
Scale-free;
MAPE Mean Absolute Percentage Error — (1/n)Σ eₜ/yₜ cannot be
used when
yₜ=0.
Mean Absolute Scaled Error — MAE / MAE_naïve. MASE<1
MASE
means better than naive forecast.
ME Mean Error — (1/n)Σeₜ. Measures bias; positive = over-
forecast; negative = under-forecast.
 Term / Symbol Full Form / Definition
Function in fpp2 for the Mean forecast method (constant at
meanf()
historical mean)
Dataset in fpp2: weekly economy class passengers on
melsyd
Melbourne-Sydney airline (1987–1992)
Maximum Likelihood Estimation — method of parameter
MLE estimation that finds values most likely to have generated the
observed data
MPE Mean Percentage Error — (1/n)Σ(eₜ/yₜ)×100%. Measures
average percentage bias.
MSE Mean Squared Error — (1/n)Σeₜ². Penalizes large errors
heavily. Used in cross-validation.
Mean — the average value of a distribution or series. In
μ (mu)
GARCH, μ = constant mean return.
Multiplicative Decomposition where Y = T × S × R. Used when seasonal
Model variation grows proportionally with the trend level.
Function in fpp2 for the Naive forecast method (next value =
naive()
last observed value)
Function in fpp2/forecast that determines the number of regular
ndiffs()
differences needed for stationarity
Function in fpp2/forecast that determines the number of
nsdiffs()
seasonal differences needed
Degrees of freedom parameter in Student's t-distribution.
ν (nu)
Controls tail thickness; lower ν = fatter tails.
Long-run average variance in GARCH — the baseline variance
ω (omega)
intercept (must be > 0).
Ordinary Least Squares — regression method minimizing the
OLS
sum of squared residuals
In ugarchfit(): the number of observations held out from fitting
[Link]
for out-of-sample evaluation
1. Non-seasonal AR order in ARIMA. 2. ARCH term order in
p
GARCH.
The probability of observing a test statistic at least as extreme
p-value
as the one computed, assuming H₀ is true. If p < 0.05, reject H₀.
P Seasonal AR order in SARIMA(p,d,q)(P,D,Q)[m]
Partial AutoCorrelation Function — correlation between yₜ and
PACF yₜ₋ₖ after controlling for intermediate lags. Used to identify AR
order (p).
Base R parameter setting for multi-panel plotting.
par(mfrow)
par(mfrow=c(r,c)) creates an r×c grid of plots.
φ (phi) 1. AR coefficient (φ₁, φ₂, ..., φₚ). 2. Damping parameter in Holt's
damped method (0 < φ < 1).
Prediction Interval — a range within which a future observation
PI is expected to fall with a specified probability (e.g., 80% or
95%).
[Link]() Base R function for plotting time series objects
1. Non-seasonal MA order in ARIMA. 2. GARCH term order in
q
GARCH.
Q Seasonal MA order in SARIMA(p,d,q)(P,D,Q)[m]
qplot() Quick plot function from ggplot2 for simple scatter plots
R package for downloading and analysing financial data
quantmod
(Quantitative Financial Modelling)
R-squared — proportion of variance in the dependent variable
R²
explained by the regression model (0 to 1). Higher = better fit.
 Term / Symbol Full Form / Definition
1. Random/Irregular component in decomposition. 2. Financial
Rₜ
return at time t.
A process where each value is the previous value plus a
Random Walk random shock: yₜ = yₜ₋₁ + εₜ. Non-stationary; requires
differencing.
residuals() Computes forecast errors (actual minus fitted): eₜ = yₜ − ŷₜ
Root Mean Squared Error — √[(1/n)Σe₮²]. The most commonly
RMSE
used accuracy metric; penalizes large errors.
rugarch R package for Univariate GARCH modelling
Random Walk Forecast — with drift=TRUE implements the Drift
rwf()
method
R The residual component after removing Trend and Seasonal
(Random/Irregular) from a time series. Should look like white noise.
The repeating periodic pattern in a time series (yearly, quarterly,
S (Seasonal)
monthly cycle)
Seasonal ARIMA — ARIMA(p,d,q)(P,D,Q)[m]. Extends ARIMA
SARIMA
with seasonal AR and MA components.
Base R function to read data from a file or URL into a numeric
scan()
vector
Function from the seasonal package for X11 or SEATS
seas()
decomposition
R package interfacing with the US Census Bureau's X-
seasonal
13ARIMA-SEATS software
Seasonal Extraction in ARIMA Time Series — a model-based
SEATS
decomposition method using ARIMA models
Simple Exponential Smoothing — for data with no trend or
SES
seasonality. One parameter: α. ŷₜ₊ₕ = lₜ.
ses() Function in fpp2 for Simple Exponential Smoothing
sGARCH Standard GARCH — the basic symmetric GARCH model
Conditional standard deviation at time t. σₜ² = conditional
σₜ (sigma)
variance modelled by GARCH.
Simple Return Rₜ = (Pₜ − Pₜ₋₁)/Pₜ₋₁ — percentage price change per period
Seasonal Naive forecast — next value equals the value from
snaive()
the same season last year
Sum of Squared Errors — Σ(yₜ − ŷₜ)². Minimized in exponential
SSE
smoothing parameter estimation.
Skewed Student's t-distribution — used in GARCH to model fat
sstd
tails and asymmetry in financial returns
Property of a time series whose statistical properties (mean,
Stationarity variance, autocovariance) do not change over time. Required
for ARIMA.
Seasonal and Trend decomposition using Loess — a flexible,
STL robust decomposition method. Parameters: [Link], [Link],
robust.
STL Forecast — decomposes with STL, forecasts seasonally-
stlf()
adjusted data, re-adds seasonal component
STL parameter controlling seasonal component smoothness.
[Link]
"periodic" = constant seasonality; integer = allowed to change.
The long-run direction (upward, downward, or flat) of a time
T (Trend)
series
> 2 roughly
t-statistic Test statistic = coefficient / standard error. t corresponds
to p < 0.05.
 Term / Symbol Full Form / Definition
θ (theta) MA coefficient in ARIMA (θ₁, θ₂, ..., θq). Applied to lagged error
terms.
Creates a time series object in R from a numeric vector, with
ts()
attributes: start, frequency
Time Series Cross-Validation — rolling-origin evaluation of
tsCV()
forecast methods. Gives honest out-of-sample error estimates.
Time Series Linear Model — fits an OLS regression to time
tslm()
series, supporting trend and season terms
TSLA Tesla Inc. stock ticker symbol on NASDAQ
STL parameter controlling trend component smoothness.
[Link]
Larger = smoother trend. Must be odd.
Fits a GARCH model specification to data using MLE (from
ugarchfit()
rugarch package)
Generates forecasts of mean and variance from a fitted
ugarchforecast()
GARCH model
Specifies the structure of a GARCH model (mean model,
ugarchspec()
variance model, distribution)
A non-stationarity condition where the AR polynomial has a root
Unit Root
equal to 1. A random walk has a unit root.
R package for Unit Root and Cointegration Analysis. Contains
urca
[Link](), [Link](), etc.
[Link]() KPSS stationarity test from the urca package
Dataset in fpp2: quarterly % changes in US consumption,
uschange
income, production, savings, unemployment
Dataset in fpp2: US monthly electricity generation (billion kWh,
usmelec
1973–2013)
Variance Equation In GARCH: σₜ² = ω + αεₜ₋₁² + βσₜ₋₁². Describes how conditional
variance evolves over time.
Dataset in fpp2: quarterly visitor nights in Australian regions
visnights
(millions)
The empirical observation that large changes in financial
Volatility
returns tend to be followed by large changes (of either sign),
Clustering
and small changes by small changes.
A sequence of independent, identically distributed random
White Noise variables with mean 0 and constant variance. Ideal residuals
from a good model.
Extracts a sub-window of a ts object between specified start
window()
and end dates
A decomposition method developed by the US Census Bureau
X11 using iterative centred moving averages. More robust than
Classical.
eXtensible Time Series — a flexible time series class in R used
xts
for financial data, supporting date-based indexing
Forecast/fitted value at time t (y-hat) — the model's prediction
ŷₜ
of yₜ
yₜ Observed value at time t
Skewness parameter in skewed distributions (sstd, snorm) used
ξ (xi)
in GARCH
 Final Note: Time series analysis is a journey from raw data → understanding structure (decomposition, stationarity) → building
models (ARIMA, SARIMA, GARCH) → evaluating and forecasting. Each lab in this guide represents one step in that journey. Revisit
any section when working on real data — the concepts build on each other.

Guide compiled for MBA Time Series Analysis students. All R code is based on lab sessions using R packages: fpp2, GGally, forecast,
urca, seasonal, rugarch, quantmod, xts, PerformanceAnalytics, dplyr, tidyverse, tseries.