Load packages
library(tidyverse)
library(astsa)
library(fpp3)
library(epidatasets)
Periodic processes
n = 100
t = 1:n
x1 = sqrt(2^2 + 3^2)*cos(2*pi*t*3/n + atan(-3/2))
y1 = 2*cos(2*pi*t*3/n) + 3*sin(2*pi*t*3/n)
x2 = sqrt(4^2 + 5^2)*cos(2*pi*t*6/n + atan(-5/4))
y2 = 4*cos(2*pi*t*6/n) + 5*sin(2*pi*t*6/n)
par(mar = c(4.5, 4.5, 0.5, 0.5))
matplot(t, cbind(x1, y1), type = "l", lty = 1:2, col = 1:2, xlab = "Time",
ylab = "", ylim = c(-10, 10))
matplot(t, cbind(x2, y2), type = "l", lty = 1:2, col = c(8,4), add = TRUE)
rect(-atan(-3/2)*n/3*1/(2*pi), -20, n/3 - atan(-3/2)*n/3*1/(2*pi), 20,
border = NA, col = rgb(1, 0, 0, 0.2))
rect(3*n/6 - atan(-5/4)*n/6*1/(2*pi), -20, 4*n/6 - atan(-6/4)*n/6*1/(2*pi), 20,
border = NA, col = rgb(0, 0, 1, 0.2))
Mixtures
par(mfrow = c(2, 1), mar = c(4.5, 4.5, 0.5, 0.5))
plot(t, y1 + y2, type = "l", xlab = "Time", ylab = "")
y3 = 6*cos(2*pi*t*18/100) + 7*sin(2*pi*t*18/100)
plot(t, y1 + y2 + y3, type = "l", xlab = "Time", ylab = "")
Periodogram, mixture data
P = Mod(fft(y1 + y2))^2 / n
Q = Mod(fft(y1 + y2 + y3))^2 / n
par(mfrow = c(2, 1), mar = c(4.5, 4.5, 0.5, 0.5))
plot(0:(n-1)/n, P, type = "h", xlab = "Frequency", ylab = "Periodogram",
xlim = c(0, 0.5), lwd = 3)
abline(v = 0.5, lty = 2, col = 2)
plot(0:(n-1)/n, Q, type = "h", xlab = "Frequency", ylab = "Periodogram",
xlim = c(0, 0.5), lwd = 3)
abline(v = 0.5, lty = 2, col = 2)
Periodogram, star data
n = length(star)
Per = Mod(fft(star - mean(star)))^2 / n
Freq = 0:(n-1)/n
par(mfrow = c(2, 1), mar = c(4.5, 4.5, 0.5, 0.5))
plot(star, ylab="Star magnitude", xlab = "Day")
plot(Freq, Per, type = "h", xlab = "Frequency", ylab = "Periodogram",
xlim = c(0, 0.08), lwd = 3)
text(0.05, 7000, "24 day cycle")
text(0.027, 9000, "29 day cycle")
ST decomposition
us_retail_employment <- us_employment |>
filter(year(Month) >= 1990, Title == "Retail Trade") |>
select(-Series_ID, -Title)
x = as.Date(us_retail_employment$Month)
y = us_retail_employment$Employed
# Compute and plot HP filter solution at decently large lambda value
n = nrow(us_retail_employment)
D = diag(rep(-2,n)) # -2s on the diagonal
D[row(D) == col(D)-1] = 1 # 1s above the diagonal
D[row(D) == col(D)+1] = 1 # 1s below the diagonal
D = D[-c(1,n), ] # Drop first and last row
I = diag(n) # n x n identity matrix
lam = 1000
hp = solve(I + lam * t(D) %*% D, y)
par(mfrow = c(2, 1), mar = c(4.5, 4.5, 0.5, 0.5))
plot(x, y, col = 8, type = "l", xlab = "Date", ylab = "Employed")
lines(x, hp, type = "l", lty = 1, lwd = 2, col = 2)
# Compute and plot periodogram and pick out 4 largest components (but ignoring
# ones that bunched up together)
d = fft(y - hp)
Per = Mod(d)^2 / n
Freq = 0:(n-1)/n
plot(Freq, Per, type = "h", xlab = "Frequency", ylab = "Periodogram",
xlim = c(0, 0.5), lwd = 2)
ord = order(Per[Freq <= 1/2], decreasing = TRUE)
cbind(Per[ord][1:5], Freq[ord][1:5])
## [,1] [,2]
## [1,] 3092160.9 0.08403361
## [2,] 1845802.0 0.16806723
## [3,] 1672142.4 0.16526611
## [4,] 921830.6 0.24929972
## [5,] 369734.6 0.33333333
ind = ord[c(1, 2, 4, 5)]
text(Freq[ind] + 0.005, Per[ind], adj = c(0, 0.5),
paste(round(1/Freq[ind], 2), "month cycle"))
Spectral density, moving average
theta = seq(0, 0.9, length = 8)
omega = seq(0, 1/2, length = 500)
fmat = matrix(NA, length(omega), length(theta))
for (j in 1:ncol(fmat)) {
fmat[,j] = (1 + theta[j])^2 + 2* theta[j] * cos(2*pi*omega)
}
par(mar = c(4.5, 4.5, 0.5, 0.5))
matplot(omega, fmat, type = "l", lty = 1, col = 1:8,
xlab = "Frequency", ylab = "Spectral density")
legend("topright", lty = 1, col = 1:8,
legend = paste("theta =", round(theta, 2)))
Spectral density, autoregressive
phi1 = c(1, 0.6, 0.5, -0.5)
phi2 = c(-0.9, -0.9, 0.4, 0.4)
omega = seq(0, 1/2, length = 500)
fmat = matrix(NA, length(omega), length(theta))
for (j in 1:ncol(fmat)) {
fmat[,j] = 1 / ((1 + phi1[j]^2 + phi2[j]^2) -
2 * phi1[j] * (1 - phi2[j]) * cos(2*pi*omega) -
2 * phi2[j] * cos(4*pi*omega))
}
par(mar = c(4.5, 4.5, 0.5, 0.5))
matplot(omega, fmat, type = "l", lty = 1, col = 1:8,
xlab = "Frequency", ylab = "Spectral density")
legend("topright", lty = 1, col = 1:8,
legend = paste0("(phi1, phi2) = (", round(phi1, 2), ", ",
round(phi2, 2), ")"))
