1 Problem 1: Simulation of Time Series Models

1.1 1.1 Generating the datasets

We generate 1000 observations from each of the seven given models. The first 500 observations are discarded to reduce the effect of initial values.

n <- 1000
burn <- 500

# Store the AR coefficients, MA coefficients and constant terms
models <- list(
  a = list(ar = 0.5, ma = 0.75, c = 9),
  b = list(ar = 0.5, ma = c(0.75, 0.2), c = 9),
  c = list(ar = c(-0.2, 0.7), ma = c(0.75, 0.2, 1.1), c = 0),
  d = list(ar = c(-0.4, 0.7), ma = c(0.75, 0.2, 1.1), c = 0),
  e = list(ar = numeric(0), ma = c(0.75, 0.2, 1.1), c = 3.3),
  f = list(ar = numeric(0), ma = numeric(0), c = 0),
  g = list(ar = numeric(0), ma = c(0.75, 0.2), c = 0)
)

# Check stationarity using the roots of the AR polynomial
is_stationary <- function(ar) {

  if (length(ar) == 0) return(TRUE)

  roots <- polyroot(c(1, -ar))

  all(Mod(roots) > 1)
}

# Generate observations using the given ARMA equation
simulate_model <- function(model) {

  ar <- model$ar
  ma <- model$ma
  constant <- model$c

  p <- length(ar)
  q <- length(ma)

  total <- n + burn

  # Generate independent random errors from N(0,1)
  errors <- rnorm(total + q)

  y <- numeric(total)

  # Calculate observations recursively
  for (t in seq_len(total)) {

    ar_part <- 0
    ma_part <- errors[t + q]

    # AR component
    if (p > 0) {
      for (j in seq_len(p)) {
        if (t > j) {
          ar_part <- ar_part + ar[j] * y[t-j]
        }
      }
    }

    # MA component
    if (q > 0) {
      for (j in seq_len(q)) {
        ma_part <- ma_part + ma[j] * errors[t + q-j]
      }
    }

    y[t] <- constant + ar_part + ma_part
  }

  # Remove the first 500 observations
  tail(y, n)
}

# Generate all seven datasets
simulated <- lapply(models, simulate_model)

names(simulated) <- names(models)

# Display the stationarity results
stationarity_table <- data.frame(
  Model = names(models),
  Stationary = sapply(
    models,
    function(m) is_stationary(m$ar)
  )
)

kable(stationarity_table, row.names = FALSE)
Model Stationary
a TRUE
b TRUE
c TRUE
d FALSE
e TRUE
f TRUE
g TRUE 
# Save the stationarity table
write.csv(
  stationarity_table,
  "PS2_Output/Tables/Stationarity.csv",
  row.names = FALSE
)

1.2 1.2 Time series plots

par(mfrow = c(4, 2), mar = c(4, 4, 3, 1))

# Plot all seven simulated datasets
for (name in names(simulated)) {

  x <- simulated[[name]]

  plot(
    seq_along(x), x,
    type = "l",
    main = paste("Model", toupper(name)),
    xlab = "Observation",
    ylab = "Value",
    col = "steelblue",
    xaxt = "n"
  )

  # Use integer observation numbers on the x-axis
  ticks <- pretty(c(1, length(x)), n = 6)
  ticks <- ticks[ticks >= 1 & ticks <= length(x)]

  axis(1, at = ticks, labels = as.integer(ticks))
}

par(mfrow = c(1, 1))

Observation: Models (a) and (b) fluctuate around a positive mean. Model (e) fluctuates around 3.3, while models (c), (f), and (g) fluctuate around zero. Model (d) shows nonstationary behaviour, with increasing variability over time.

1.3 1.3 Estimation of process mean

For a stationary ARMA process, the theoretical mean is

\[ \mu = \frac{c}{1-\phi_1-\phi_2-\cdots-\phi_p}. \]

# Select only stationary models
stationary_models <- names(models)[
  sapply(models, function(m) is_stationary(m$ar))
]

# Calculate theoretical and sample means
mean_table <- data.frame(
  Model = stationary_models,

  True_Mean = sapply(
    models[stationary_models],
    function(m) m$c / (1 - sum(m$ar))
  ),

  Sample_Mean = sapply(
    simulated[stationary_models],
    mean
  )
)

# Calculate the difference between the two means
mean_table$Difference <-
  mean_table$Sample_Mean - mean_table$True_Mean

kable(mean_table, digits = 4, row.names = FALSE)
Model True_Mean Sample_Mean Difference
a 18.0 18.0415 0.0415
b 18.0 17.9199 -0.0801
c 0.0 -0.1162 -0.1162
e 3.3 3.2873 -0.0127
f 0.0 0.0235 0.0235
g 0.0 -0.0140 -0.0140 
# Save the mean comparison table
write.csv(
  mean_table,
  "PS2_Output/Tables/Mean_Comparison.csv",
  row.names = FALSE
)

Observation: The sample means are generally close to their theoretical values. Small differences occur due to random variation in the simulated observations.

1.4 1.4 Sample ACF plots

par(mfrow = c(4, 2), mar = c(4, 4, 3, 1))

# ACF shows correlation between observations at different lags
for (name in names(simulated)) {

  a <- acf(
    simulated[[name]],
    lag.max = 30,
    plot = FALSE
  )

  plot(
    a,
    main = paste("ACF - Model", toupper(name)),
    xaxt = "n"
  )

  # Display integer lag values
  axis(1, at = seq(0, 30, by = 5),
       labels = seq(0, 30, by = 5))
}

par(mfrow = c(1, 1))

1.5 1.5 Sample PACF plots

par(mfrow = c(4, 2), mar = c(4, 4, 3, 1))

# PACF measures correlation after removing intermediate lag effects
for (name in names(simulated)) {

  p <- pacf(
    simulated[[name]],
    lag.max = 30,
    plot = FALSE
  )

  plot(
    p,
    main = paste("PACF - Model", toupper(name)),
    xaxt = "n"
  )

  # Display integer lag values
  axis(1, at = seq(1, 30, by = 5),
       labels = seq(1, 30, by = 5))
}

par(mfrow = c(1, 1))

1.6 1.6 Comparison of ACF and PACF

Model Process Expected pattern
(a) ARMA(1,1) Both generally tail off
(b) ARMA(1,2) Both generally tail off
(c) ARMA(2,3) Both generally tail off
(d) Nonstationary ARMA(2,3) Persistent or growing dependence
(e) MA(3) ACF cuts off after lag 3
(f) White noise ACF and PACF near zero
(g) MA(2) ACF cuts off after lag 2

Observation: MA models generally show a cutoff in the ACF after their MA order. ARMA models show gradual decay in both ACF and PACF. White noise has no meaningful autocorrelation at non-zero lags. The nonstationary model may show persistent dependence.

2 Problem 2: Identification of Models from Given Datasets

2.1 2.1 Importing the CSV datasets

All five CSV files should be kept in the same folder as this R Markdown file. Each file contains one column of observations without a header.

# Import the five datasets
y1 <- read.csv("y1 - y1.csv", header = FALSE)[[1]]
y2 <- read.csv("y2 - y2.csv", header = FALSE)[[1]]
y3 <- read.csv("y3 - y3.csv", header = FALSE)[[1]]
y4 <- read.csv("y4 - y4.csv", header = FALSE)[[1]]
y5 <- read.csv("y5 - y5.csv", header = FALSE)[[1]]

# Convert the datasets into time series
datasets <- list(
  y1 = ts(as.numeric(y1)),
  y2 = ts(as.numeric(y2)),
  y3 = ts(as.numeric(y3)),
  y4 = ts(as.numeric(y4)),
  y5 = ts(as.numeric(y5))
)

# Display basic information about the datasets
data_summary <- data.frame(
  Dataset = names(datasets),
  Observations = sapply(datasets, length),
  Mean = sapply(datasets, mean),
  SD = sapply(datasets, sd)
)

kable(data_summary, digits = 4, row.names = FALSE)
Dataset Observations Mean SD
y1 800 -0.0984 1.7366
y2 300 0.1102 1.4593
y3 500 -0.2209 10.1554
y4 1000 3.3454 2.8311
y5 800 108361716055497466280884004204002664648208.0000 660255097369126116160262482682442880602682.0000 
# Save the dataset summary
write.csv(
  data_summary,
  "PS2_Output/Tables/Dataset_Summary.csv",
  row.names = FALSE
)

2.2 2.2 Time series plots

par(mfrow = c(5, 1), mar = c(4, 4, 3, 1))

# Plot all five datasets
for (name in names(datasets)) {

  x <- as.numeric(datasets[[name]])

  plot(
    seq_along(x), x,
    type = "l",
    main = paste("Time Series -", name),
    xlab = "Observation",
    ylab = "Value",
    col = "steelblue",
    xaxt = "n"
  )

  # Use actual integer observation numbers
  ticks <- pretty(c(1, length(x)), n = 6)
  ticks <- ticks[ticks >= 1 & ticks <= length(x)]

  axis(1, at = ticks, labels = as.integer(ticks))
}

par(mfrow = c(1, 1))

Observation: The time series plots provide an initial understanding of each dataset. They help identify trends, fluctuations, and possible changes in variability.

2.3 2.3 ACF and PACF plots

par(mfrow = c(5, 2), mar = c(4, 4, 3, 1))

# Plot ACF and PACF for each dataset
for (name in names(datasets)) {

  x <- datasets[[name]]

  # ACF plot
  a <- acf(x, lag.max = 30, plot = FALSE)

  plot(
    a,
    main = paste(name, "- ACF"),
    xaxt = "n"
  )

  axis(1, at = seq(0, 30, by = 5),
       labels = seq(0, 30, by = 5))

  # PACF plot
  p <- pacf(x, lag.max = 30, plot = FALSE)

  plot(
    p,
    main = paste(name, "- PACF"),
    xaxt = "n"
  )

  axis(1, at = seq(1, 30, by = 5),
       labels = seq(1, 30, by = 5))
}

par(mfrow = c(1, 1))

Observation: ACF and PACF plots help identify possible model orders. An ACF cutoff suggests an MA model, while a PACF cutoff suggests an AR model. If both gradually decay, an ARMA model may be suitable.

2.4 2.4 Fitting candidate ARIMA models

We fit candidate ARIMA models and compare them using AICc. A smaller AICc indicates a better model among the candidates considered.

# Calculate corrected AIC
calculate_aicc <- function(fit, n) {

  k <- length(fit$coef) + 1

  if (n <= k + 1) return(Inf)

  fit$aic + (2 * k * (k + 1)) / (n - k - 1)
}

# Fit ARIMA models with different orders
find_models <- function(x) {

  results <- data.frame()
  fits <- list()

  # Try differencing orders 0 and 1
  for (d in 0:1) {

    for (p in 0:3) {

      for (q in 0:3) {

        fit <- tryCatch(
          suppressWarnings(
            arima(
              x,
              order = c(p, d, q),
              include.mean = (d == 0),
              method = "ML",
              optim.control = list(maxit = 2000)
            )
          ),
          error = function(e) NULL
        )

        # Keep models that were fitted successfully
        if (!is.null(fit) &&
            is.finite(fit$aic) &&
            all(is.finite(fit$coef))) {

          model_name <- paste0(
            "ARIMA(", p, ",", d, ",", q, ")"
          )

          value <- calculate_aicc(fit, length(x))

          if (is.finite(value)) {

            results <- rbind(
              results,
              data.frame(
                Model = model_name,
                p = p,
                d = d,
                q = q,
                AICc = value
              )
            )

            fits[[model_name]] <- fit
          }
        }
      }
    }
  }

  # Stop if no model could be fitted
  if (nrow(results) == 0) {
    stop("No candidate model could be fitted. Check the dataset.")
  }

  # Arrange models from lowest to highest AICc
  results <- results[order(results$AICc), ]
  rownames(results) <- NULL

  list(
    table = results,
    fits = fits,
    best = fits[[results$Model[1]]]
  )
}

# Check the datasets for missing or infinite values
for (name in names(datasets)) {

  if (any(!is.finite(datasets[[name]]))) {
    stop(paste(name, "contains missing or infinite values."))
  }

  if (sd(datasets[[name]]) == 0) {
    stop(paste(name, "has zero variance."))
  }
}

# Fit candidate models for all five datasets
selection <- lapply(datasets, find_models)

# Display and save the top five candidate models
for (name in names(selection)) {

  cat("\n\nDataset:", name, "\n\n")

  print(
    kable(
      head(selection[[name]]$table, 5),
      digits = 3,
      row.names = FALSE
    )
  )

  write.csv(
    selection[[name]]$table,
    paste0(
      "PS2_Output/Tables/",
      name,
      "_Model_Comparison.csv"
    ),
    row.names = FALSE
  )
}
## 
## 
## Dataset: y1 
## 
## 
## 
## |Model        |  p|  d|  q|     AICc|
## |:------------|--:|--:|--:|--------:|
## |ARIMA(3,0,1) |  3|  0|  1| 2238.185|
## |ARIMA(2,0,0) |  2|  0|  0| 2239.683|
## |ARIMA(3,0,0) |  3|  0|  0| 2240.066|
## |ARIMA(2,0,1) |  2|  0|  1| 2240.135|
## |ARIMA(3,0,2) |  3|  0|  2| 2241.440|
## 
## 
## Dataset: y2 
## 
## 
## 
## |Model        |  p|  d|  q|    AICc|
## |:------------|--:|--:|--:|-------:|
## |ARIMA(1,0,2) |  1|  0|  2| 872.634|
## |ARIMA(0,0,1) |  0|  0|  1| 873.062|
## |ARIMA(2,0,1) |  2|  0|  1| 873.188|
## |ARIMA(0,0,3) |  0|  0|  3| 873.759|
## |ARIMA(1,0,1) |  1|  0|  1| 874.090|
## 
## 
## Dataset: y3 
## 
## 
## 
## |Model        |  p|  d|  q|     AICc|
## |:------------|--:|--:|--:|--------:|
## |ARIMA(2,0,2) |  2|  0|  2| 3737.412|
## |ARIMA(0,1,1) |  0|  1|  1| 3739.640|
## |ARIMA(3,1,3) |  3|  1|  3| 3739.891|
## |ARIMA(0,0,0) |  0|  0|  0| 3739.967|
## |ARIMA(3,0,2) |  3|  0|  2| 3740.200|
## 
## 
## Dataset: y4 
## 
## 
## 
## |Model        |  p|  d|  q|     AICc|
## |:------------|--:|--:|--:|--------:|
## |ARIMA(3,0,3) |  3|  0|  3| 4397.542|
## |ARIMA(3,0,2) |  3|  0|  2| 4426.897|
## |ARIMA(2,0,3) |  2|  0|  3| 4427.733|
## |ARIMA(3,1,3) |  3|  1|  3| 4428.046|
## |ARIMA(3,0,0) |  3|  0|  0| 4449.607|
## 
## 
## Dataset: y5 
## 
## 
## 
## |Model        |  p|  d|  q|     AICc|
## |:------------|--:|--:|--:|--------:|
## |ARIMA(2,1,3) |  2|  1|  3| 146159.5|
## |ARIMA(2,1,2) |  2|  1|  2| 146260.2|
## |ARIMA(3,1,2) |  3|  1|  2| 146306.1|
## |ARIMA(2,1,1) |  2|  1|  1| 146811.3|
## |ARIMA(2,1,0) |  2|  1|  0| 147196.0|

2.5 2.5 Selecting the best model

# Select the model with the lowest AICc for each dataset
best_models <- do.call(
  rbind,
  lapply(names(selection), function(name) {

    result <- selection[[name]]
    fit <- result$best

    p <- result$table$p[1]
    q <- result$table$q[1]

    # Extract model residuals
    res <- na.omit(as.numeric(residuals(fit)))

    # Perform the Ljung-Box test
    lag_value <- min(20, floor(length(res) / 5))

    test <- Box.test(
      res,
      lag = lag_value,
      type = "Ljung-Box",
      fitdf = p + q
    )

    data.frame(
      Dataset = name,
      Best_Model = result$table$Model[1],
      AICc = result$table$AICc[1],
      Ljung_Box_p = test$p.value
    )
  })
)

rownames(best_models) <- NULL

kable(
  best_models,
  digits = 4,
  row.names = FALSE,
  caption = "Selected models and residual checks"
)
Selected models and residual checks
Dataset Best_Model AICc Ljung_Box_p
y1 ARIMA(3,0,1) 2238.1848 0.0350
y2 ARIMA(1,0,2) 872.6338 0.4253
y3 ARIMA(2,0,2) 3737.4121 0.4415
y4 ARIMA(3,0,3) 4397.5417 0.0000
y5 ARIMA(2,1,3) 146159.5040 0.0000 
# Save the final model comparison
write.csv(
  best_models,
  "PS2_Output/Tables/Best_Models.csv",
  row.names = FALSE
)

2.6 2.6 Residual diagnostics

Residuals should behave approximately like white noise if the selected model adequately captures the dependence in the data.

par(mfrow = c(5, 2), mar = c(4, 4, 3, 1))

# Plot residuals and their ACF for each selected model
for (name in names(selection)) {

  res <- as.numeric(residuals(selection[[name]]$best))

  # Residual time series
  plot(
    seq_along(res), res,
    type = "l",
    main = paste(name, "- Residuals"),
    xlab = "Observation",
    ylab = "Residual",
    col = "steelblue",
    xaxt = "n"
  )

  ticks <- pretty(c(1, length(res)), n = 5)
  ticks <- ticks[ticks >= 1 & ticks <= length(res)]

  axis(1, at = ticks, labels = as.integer(ticks))

  # Residual ACF
  acf(
    res,
    lag.max = 30,
    main = paste(name, "- Residual ACF")
  )
}

par(mfrow = c(1, 1))

2.7 2.7 Final observations

# Print a short observation for each selected model
for (i in seq_len(nrow(best_models))) {

  name <- best_models$Dataset[i]
  model <- best_models$Best_Model[i]
  p_value <- best_models$Ljung_Box_p[i]

  cat(
    "**", name, ":** The selected model is ",
    model, " based on the lowest AICc. ",
    sep = ""
  )

  if (p_value > 0.05) {

    cat(
      "The Ljung-Box test does not indicate significant remaining autocorrelation.\n\n"
    )

  } else {

    cat(
      "The Ljung-Box test indicates remaining autocorrelation, so the model may need improvement.\n\n"
    )
  }
}

y1: The selected model is ARIMA(3,0,1) based on the lowest AICc. The Ljung-Box test indicates remaining autocorrelation, so the model may need improvement.

y2: The selected model is ARIMA(1,0,2) based on the lowest AICc. The Ljung-Box test does not indicate significant remaining autocorrelation.

y3: The selected model is ARIMA(2,0,2) based on the lowest AICc. The Ljung-Box test does not indicate significant remaining autocorrelation.

y4: The selected model is ARIMA(3,0,3) based on the lowest AICc. The Ljung-Box test indicates remaining autocorrelation, so the model may need improvement.

y5: The selected model is ARIMA(2,1,3) based on the lowest AICc. The Ljung-Box test indicates remaining autocorrelation, so the model may need improvement.

3 Conclusion

Seven time series models were simulated and studied using time series plots, sample means, ACF, and PACF. The estimated means were compared with their theoretical values for stationary processes.

For the five given datasets, ACF and PACF plots were used to understand their dependence structures. Candidate ARIMA models were compared using AICc, and residual diagnostics were used to assess the selected models.