Problem 1: Moving Average (MA) Models

Aim

To simulate four MA models, study their time-series plots, estimate their means and autocorrelations, identify their orders using ACF plots, and examine the effect of increasing the white-noise variance.

1.1 Model definitions

The four models have the form

\[ Y_t = 3.3 + \varepsilon_t + \theta_1\varepsilon_{t-1} + \cdots + \theta_q\varepsilon_{t-q}, \]

where initially \(\varepsilon_t \sim N(0,1)\). Their orders are 1, 2, 3 and 4, respectively. The true mean of each process is 3.3.

1.2 Function to generate an MA series

generate_ma <- function(n, theta, mu = 3.3, noise_sd = 1) {
  q <- length(theta)
  e <- rnorm(n + q, mean = 0, sd = noise_sd)

  y <- numeric(n)

  for (t in seq_len(n)) {
    y[t] <- mu + e[t + q]

    for (j in seq_len(q)) {
      y[t] <- y[t] + theta[j] * e[t + q - j]
    }
  }

  y
}

1.3 Generate four datasets of size 500

ma_coefficients <- list(
  "MA(1)" = c(0.75),
  "MA(2)" = c(0.75, 0.2),
  "MA(3)" = c(0.75, 0.2, 1.1),
  "MA(4)" = c(0.75, 0.2, 1.1, 0.33)
)

ma_500 <- lapply(
  ma_coefficients,
  function(theta) generate_ma(500, theta)
)

(i) Plot the generated datasets

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

for (model in names(ma_500)) {
  plot(
    ma_500[[model]],
    type = "l",
    main = model,
    xlab = "Time",
    ylab = expression(Y[t])
  )
  abline(h = 3.3, col = "red", lty = 2)
}

par(old_par)

Observation: Each series fluctuates around 3.3. The observations are random, but the current value depends on a finite number of current and past white-noise terms. An MA process with finite order is stationary.

(ii) Estimate the mean and autocorrelations up to lag 20

ma_mean_500 <- data.frame(
  Model = names(ma_500),
  True_Mean = 3.3,
  Sample_Mean = sapply(ma_500, mean)
)

kable(ma_mean_500, digits = 4,
      caption = "Process mean estimates (n = 500)")
Process mean estimates (n = 500)
Model True_Mean Sample_Mean
MA(1) MA(1) 3.3 3.3605
MA(2) MA(2) 3.3 3.2917
MA(3) MA(3) 3.3 3.3980
MA(4) MA(4) 3.3 3.4741

The ACF function below returns sample autocorrelations at lags 1 to 20.

get_acf <- function(x, max_lag = 20) {
  as.numeric(acf(x, lag.max = max_lag, plot = FALSE)$acf)[-1]
}
old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

for (model in names(ma_500)) {
  acf(
    ma_500[[model]],
    lag.max = 20,
    main = paste("ACF:", model)
  )
}

par(old_par)
ma_acf_500 <- data.frame(Lag = 1:20)

for (model in names(ma_500)) {
  ma_acf_500[[model]] <- get_acf(ma_500[[model]])
}

kable(ma_acf_500, digits = 4,
      caption = "Sample ACF values (n = 500)")
Sample ACF values (n = 500)
Lag MA(1) MA(2) MA(3) MA(4)
1 0.4197 0.5537 0.3977 0.3662
2 -0.0842 0.0998 0.3802 0.2343
3 -0.0522 -0.0284 0.4187 0.4186
4 -0.0641 0.0196 -0.0061 0.0035
5 -0.0654 0.0294 0.0555 -0.1004
6 -0.0504 0.0435 -0.0233 -0.0556
7 0.0031 0.0150 -0.0449 -0.0908
8 0.0540 0.0052 -0.0117 -0.0493
9 0.0517 0.0486 -0.1001 -0.0077
10 0.0271 0.0489 -0.0401 -0.0492
11 -0.0448 0.0000 -0.0278 -0.0081
12 -0.1138 -0.0274 -0.0338 0.0168
13 -0.0428 0.0103 0.0272 0.0146
14 0.0691 -0.0052 0.0021 0.0619
15 0.0716 -0.0202 -0.0304 0.0605
16 -0.0150 -0.0187 -0.0213 0.0447
17 -0.0555 0.0086 -0.0664 -0.0295
18 -0.0538 0.0535 -0.0944 -0.0136
19 -0.0217 0.0623 -0.0716 -0.0980
20 -0.0235 0.0120 -0.1033 -0.2144

Observation: For an MA(q) model, the theoretical ACF is zero at all lags greater than q. In a finite sample, some sample autocorrelations beyond q may still be non-zero because of sampling variation. The dashed limits in the ACF plots are approximate 95% reference limits.

(iii) Repeat the mean and ACF estimation for datasets of size 50,000

ma_50000 <- lapply(
  ma_coefficients,
  function(theta) generate_ma(50000, theta)
)

ma_mean_comparison <- data.frame(
  Model = names(ma_500),
  True_Mean = 3.3,
  Mean_n500 = sapply(ma_500, mean),
  Mean_n50000 = sapply(ma_50000, mean)
)

kable(ma_mean_comparison, digits = 5,
      caption = "Mean comparison for two sample sizes")
Mean comparison for two sample sizes
Model True_Mean Mean_n500 Mean_n50000
MA(1) MA(1) 3.3 3.36045 3.29610
MA(2) MA(2) 3.3 3.29166 3.30637
MA(3) MA(3) 3.3 3.39803 3.30976
MA(4) MA(4) 3.3 3.47407 3.32218 
old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

for (model in names(ma_50000)) {
  acf(
    ma_50000[[model]],
    lag.max = 20,
    main = paste("ACF:", model, "(n = 50000)")
  )
}

par(old_par)
ma_acf_50000 <- data.frame(Lag = 1:20)

for (model in names(ma_50000)) {
  ma_acf_50000[[model]] <- get_acf(ma_50000[[model]])
}

kable(ma_acf_50000, digits = 4,
      caption = "Sample ACF values (n = 50000)")
Sample ACF values (n = 50000)
Lag MA(1) MA(2) MA(3) MA(4)
1 0.4792 0.5586 0.3957 0.5031
2 -0.0018 0.1216 0.3654 0.3712
3 -0.0007 0.0013 0.3902 0.4633
4 -0.0004 -0.0034 0.0048 0.1159
5 -0.0030 -0.0111 0.0054 0.0053
6 -0.0065 -0.0134 0.0033 0.0077
7 -0.0138 -0.0142 0.0107 0.0058
8 -0.0131 -0.0078 0.0043 -0.0023
9 -0.0066 0.0020 -0.0010 -0.0005
10 -0.0032 0.0042 -0.0021 -0.0051
11 -0.0084 0.0025 -0.0082 -0.0114
12 -0.0134 -0.0046 -0.0105 -0.0025
13 -0.0066 -0.0056 -0.0090 -0.0016
14 -0.0024 -0.0031 -0.0069 -0.0007
15 -0.0018 0.0018 -0.0037 0.0037
16 0.0047 0.0063 -0.0013 0.0040
17 0.0078 0.0066 -0.0004 0.0038
18 0.0061 0.0010 0.0012 0.0040
19 0.0011 -0.0022 0.0058 0.0007
20 -0.0017 -0.0032 0.0053 0.0005

Observation and explanation: The sample means for 50,000 observations should generally be closer to 3.3 than the means for 500 observations. This follows from the law of large numbers. The standard error of the sample mean decreases as the sample size increases. For a stationary MA process,

\[ \operatorname{Var}(\bar Y) \approx \frac{\gamma(0)}{n} \left(1+2\sum_{k=1}^{q}\rho(k)\right). \]

Thus, increasing the sample size makes the estimate more stable.

(iv) Identify the MA models using their sample ACF plots

The series are given anonymous labels in a random display order. Their time order is kept unchanged, because shuffling the observations within a series would destroy the autocorrelation pattern.

set.seed(456)

display_order <- sample(names(ma_500))
unknown_series <- ma_500[display_order]

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

for (i in seq_along(unknown_series)) {
  acf(
    unknown_series[[i]],
    lag.max = 20,
    main = paste("Unknown series", i)
  )
}

par(old_par)

Use the approximate cutoff in the ACF to propose the order:

  • ACF cuts off after lag 1: candidate MA(1).
  • ACF cuts off after lag 2: candidate MA(2).
  • ACF cuts off after lag 3: candidate MA(3).
  • ACF cuts off after lag 4: candidate MA(4).

The sample ACF may not show a perfectly sharp cutoff, so use the overall pattern and the approximate confidence limits.

identification_check <- data.frame(
  Anonymous_Series = paste("Unknown", seq_along(display_order)),
  Actual_Model = unname(display_order),
  Actual_Order = sapply(ma_coefficients[display_order], length)
)

kable(identification_check,
      caption = "Cross-check of the anonymous series")
Cross-check of the anonymous series
Anonymous_Series Actual_Model Actual_Order
MA(1) Unknown 1 MA(1) 1
MA(4) Unknown 2 MA(4) 4
MA(3) Unknown 3 MA(3) 3
MA(2) Unknown 4 MA(2) 2

Observation: The ACF is useful for identifying an MA order because its theoretical autocorrelations become zero after the model order. The table above reveals the true model labels after the visual identification exercise.

(v) Change the white-noise distribution to N(0,10)

Here \(N(0,10)\) is interpreted as mean 0 and variance 10. Therefore, the standard deviation supplied to rnorm() is \(\sqrt{10}\).

set.seed(789)

ma_high_variance <- lapply(
  ma_coefficients,
  function(theta) {
    generate_ma(
      n = 500,
      theta = theta,
      mu = 3.3,
      noise_sd = sqrt(10)
    )
  }
)
old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

for (model in names(ma_high_variance)) {
  plot(
    ma_high_variance[[model]],
    type = "l",
    main = paste(model, "with noise variance 10"),
    xlab = "Time",
    ylab = expression(Y[t])
  )
  abline(h = 3.3, col = "red", lty = 2)
}

par(old_par)
high_variance_means <- data.frame(
  Model = names(ma_high_variance),
  True_Mean = 3.3,
  Sample_Mean = sapply(ma_high_variance, mean)
)

kable(high_variance_means, digits = 4,
      caption = "Mean estimates with white-noise variance 10")
Mean estimates with white-noise variance 10
Model True_Mean Sample_Mean
MA(1) MA(1) 3.3 3.2914
MA(2) MA(2) 3.3 3.2259
MA(3) MA(3) 3.3 2.6959
MA(4) MA(4) 3.3 3.5594 
old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

for (model in names(ma_high_variance)) {
  acf(
    ma_high_variance[[model]],
    lag.max = 20,
    main = paste("ACF:", model, "with variance 10")
  )
}

par(old_par)

Observation and explanation: Increasing the white-noise variance increases the variability of the observations. The process mean remains 3.3. The theoretical ACF does not change, because multiplying all autocovariances by the same variance factor leaves their ratios unchanged. Sample ACF values will differ slightly because of random sampling variation.

1.4 Conclusion for Problem 1

The four MA models were simulated and examined. Their series fluctuate around the common mean of 3.3. The theoretical ACF of an MA(q) process cuts off after lag q, which makes the ACF useful for identifying the MA order. A larger sample size gives a more stable estimate of the process mean. Increasing the white-noise variance increases the spread of the series but does not change its theoretical mean or ACF.

Problem 2: Autoregressive (AR) Models

Aim

To simulate four AR models using random initialization and a burn-in period of 200 observations, examine their time-series plots, estimate the means of stationary processes, and compare their sample ACF plots.

2.1 Model definitions

The four models are:

\[ \begin{aligned} (a)\quad Y_t &= 3+0.7Y_{t-1}+\varepsilon_t,\\ (b)\quad Y_t &= 3+1.01Y_{t-1}+\varepsilon_t,\\ (c)\quad Y_t &= 0.74Y_{t-1}-0.9Y_{t-2}+\varepsilon_t,\\ (d)\quad Y_t &= 0.74Y_{t-1}+0.9Y_{t-2}+\varepsilon_t. \end{aligned} \]

For each model, \(\varepsilon_t\sim N(0,1)\). We generate 700 observations, use the first 200 as burn-in, and retain the remaining 500.

2.2 Functions to generate AR data

generate_ar1 <- function(n, intercept, phi, burn = 200) {
  total <- n + burn
  e <- rnorm(total)
  y <- numeric(total)

  y[1] <- rnorm(1)

  for (t in 2:total) {
    y[t] <- intercept + phi * y[t - 1] + e[t]
  }

  y[(burn + 1):total]
}

generate_ar2 <- function(n, phi1, phi2, burn = 200) {
  total <- n + burn
  e <- rnorm(total)
  y <- numeric(total)

  y[1:2] <- rnorm(2)

  for (t in 3:total) {
    y[t] <- phi1 * y[t - 1] +
      phi2 * y[t - 2] + e[t]
  }

  y[(burn + 1):total]
}

2.3 Generate four datasets of size 500

set.seed(2026)

ar_data <- list(
  "AR(1), phi = 0.7" = generate_ar1(500, 3, 0.7),
  "AR(1), phi = 1.01" = generate_ar1(500, 3, 1.01),
  "AR(2), phi1 = 0.74, phi2 = -0.9" =
    generate_ar2(500, 0.74, -0.9),
  "AR(2), phi1 = 0.74, phi2 = 0.9" =
    generate_ar2(500, 0.74, 0.9)
)

(i) Plot the generated datasets

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

for (model in names(ar_data)) {
  plot(
    ar_data[[model]],
    type = "l",
    main = model,
    xlab = "Time",
    ylab = expression(Y[t])
  )
}

par(old_par)

Observation:

  • Model (a) has a stable pattern and fluctuates around its theoretical mean.
  • Model (b) has an explosive pattern because its AR coefficient is greater than 1.
  • Model (c) is stationary and may show oscillatory behaviour.
  • Model (d) is non-stationary and may show an explosive pattern.

The stationarity of the AR(2) models depends on the roots of their AR characteristic polynomial, not simply on whether the individual coefficients are smaller than 1.

(ii) Estimate the process means for stationary models

For a stationary AR(1) process,

\[ \mu=\frac{c}{1-\phi}. \]

For model (a),

\[ \mu=\frac{3}{1-0.7}=10. \]

For an AR(2) process without an intercept,

\[ \mu=\phi_1\mu+\phi_2\mu, \]

so its stationary mean is zero, provided the process is stationary.

The stationarity checks give:

  • Model (a): stationary, with mean 10.
  • Model (b): non-stationary, because \(|1.01|>1\).
  • Model (c): stationary; its characteristic roots are outside the unit circle.
  • Model (d): non-stationary; one characteristic root is inside the unit circle.
ar_mean_results <- data.frame(
  Model = c(
    "AR(1), phi = 0.7",
    "AR(1), phi = 1.01",
    "AR(2), phi1 = 0.74, phi2 = -0.9",
    "AR(2), phi1 = 0.74, phi2 = 0.9"
  ),
  Stationarity = c(
    "Stationary",
    "Non-stationary",
    "Stationary",
    "Non-stationary"
  ),
  True_Mean = c(
    "10",
    "Not defined as a stationary mean",
    "0",
    "Not defined as a stationary mean"
  ),
  Sample_Mean = sapply(ar_data, mean)
)

kable(ar_mean_results, digits = 4,
      caption = "AR process mean estimates and stationarity")
AR process mean estimates and stationarity
Model Stationarity True_Mean Sample_Mean
AR(1), phi = 0.7 AR(1), phi = 0.7 Stationary 10 1.007760e+01
AR(1), phi = 1.01 AR(1), phi = 1.01 Non-stationary Not defined as a stationary mean 6.349008e+04
AR(2), phi1 = 0.74, phi2 = -0.9 AR(2), phi1 = 0.74, phi2 = -0.9 Stationary 0 -3.230000e-02
AR(2), phi1 = 0.74, phi2 = 0.9 AR(2), phi1 = 0.74, phi2 = 0.9 Non-stationary Not defined as a stationary mean -1.347282e+97

Observation and explanation: Only models (a) and (c) have stationary process means that can be compared with their sample means. Models (b) and (d) are non-stationary, so their sample averages should not be interpreted as estimates of a constant stationary process mean.

(iii) Compare the sample ACF plots for models (a) and (b)

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

acf(
  ar_data[["AR(1), phi = 0.7"]],
  lag.max = 20,
  main = "Model (a): AR(1), phi = 0.7"
)

acf(
  ar_data[["AR(1), phi = 1.01"]],
  lag.max = 20,
  main = "Model (b): AR(1), phi = 1.01"
)

par(old_par)

Observation:

For model (a), the ACF generally decreases gradually as the lag increases, which is the usual pattern for a stationary AR(1) process.

For model (b), the coefficient is greater than 1, so the process is explosive. Its sample ACF may remain very high across many lags. The usual stationary AR(1) ACF interpretation does not apply to this model.

(iv) Compare the sample ACF plots for models (c) and (d)

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

acf(
  ar_data[["AR(2), phi1 = 0.74, phi2 = -0.9"]],
  lag.max = 20,
  main = "Model (c): AR(2), phi2 = -0.9"
)

acf(
  ar_data[["AR(2), phi1 = 0.74, phi2 = 0.9"]],
  lag.max = 20,
  main = "Model (d): AR(2), phi2 = 0.9"
)

par(old_par)

Observation:

Model (c) is stationary and can show a damped oscillating ACF pattern because of its AR coefficients.

Model (d) is non-stationary and may show a persistent or increasing dependence pattern. The stationary ACF interpretation is not valid for this model.

The ACF of an AR process generally tails off rather than cutting off after a fixed lag. Therefore, ACF plots alone are not usually sufficient to identify the order of an AR process.

The partial autocorrelation function (PACF) is commonly used as an additional identification tool for stationary AR processes. For an AR(p) process, the theoretical PACF cuts off after lag p.

2.3 Conclusion for Problem 2

The four AR models were simulated using random initial values and a burn-in period of 200 observations. Models (a) and (c) are stationary, whereas models (b) and (d) are non-stationary. The sample mean was compared with the theoretical mean only for stationary models. The ACF plots show that AR processes generally have a gradual tailing-off pattern, unlike the cutoff pattern of MA processes. Hence, the PACF is a useful additional tool for identifying stationary AR model orders.

Overall Conclusion

The simulations illustrate the main differences between Moving Average and Autoregressive time series models. MA processes have an ACF that cuts off after their order, while stationary AR processes generally have an ACF that tails off. Increasing the sample size improves mean estimation, and increasing white-noise variance changes the spread of a series without changing its theoretical mean or normalized autocorrelation structure.