Problem 1: Simulation and Analysis of Moving Average Models

Aim

To generate data from four Moving Average models, study their time series patterns, estimate their process means and autocorrelation functions, and identify the models using sample ACF plots.

1.1 Model Specification

The four MA models considered are:

Model (a): MA(1)

\[ Y_t = 3.3+\varepsilon_t+0.75\varepsilon_{t-1} \]

Model (b): MA(2)

\[ Y_t = 3.3+\varepsilon_t+0.75\varepsilon_{t-1} +0.2\varepsilon_{t-2} \]

Model (c): MA(3)

\[ Y_t = 3.3+\varepsilon_t+0.75\varepsilon_{t-1} +0.2\varepsilon_{t-2}+1.1\varepsilon_{t-3} \]

Model (d): MA(4)

\[ Y_t = 3.3+\varepsilon_t+0.75\varepsilon_{t-1} +0.2\varepsilon_{t-2}+1.1\varepsilon_{t-3} +0.33\varepsilon_{t-4} \]

Initially, the white noise sequence follows:

\[ \varepsilon_t\sim N(0,1) \]

The true process mean for all four models is 3.3.

1.2 Function to Generate MA Data

The following function generates observations from an MA(q) model using normally distributed white noise.

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 1:n) {
    y[t] <- mu + e[t + q]

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

  return(y)
}

1.3 Generate Four Datasets of Size 500

The coefficient values for the four models are specified below.

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)
)

data_500 <- lapply(
  coefficients,
  function(theta) generate_ma(500, theta)
)

(i) Plot the Generated Datasets

The generated time series are plotted to observe their patterns and variations.

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

for (model in names(data_500)) {

  plot(
    data_500[[model]],
    type = "l",
    main = paste(model, "Time Series"),
    xlab = "Time",
    ylab = expression(Y[t])
  )

  abline(h = 3.3, col = "red", lty = 2)
}

par(mfrow = c(1, 1))

Observation:

All four time series fluctuate around the common process mean of 3.3. The observations show random variations due to the white noise sequence.

As the order of the MA model increases, more past white noise terms contribute to the current observation. However, all four processes are stationary.

(ii) Estimate Process Mean and Autocorrelations up to Lag 20

The sample mean is used to estimate the process mean.

The sample autocorrelation function (ACF) is calculated up to lag 20 for each dataset.

estimate_acf <- function(x, max_lag = 20) {

  acf_values <- as.numeric(
    acf(x, lag.max = max_lag, plot = FALSE)$acf
  )

  return(acf_values[-1])
}

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

kable(
  results_500,
  digits = 4,
  caption = "Estimated Process Means for Sample Size 500"
)
Estimated Process Means for Sample Size 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

Sample ACF Plots

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

for (model in names(data_500)) {

  acf(
    data_500[[model]],
    lag.max = 20,
    main = paste("ACF of", model)
  )
}

par(mfrow = c(1, 1))

Estimated Autocorrelations

acf_table_500 <- data.frame(
  Lag = 1:20
)

for (model in names(data_500)) {
  acf_table_500[[model]] <- estimate_acf(data_500[[model]])
}

kable(
  acf_table_500,
  digits = 4,
  caption = "Sample Autocorrelations up to Lag 20 (n = 500)"
)
Sample Autocorrelations up to Lag 20 (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) process, the theoretical ACF cuts off after lag q. Therefore:

  • MA(1): ACF is zero after lag 1.
  • MA(2): ACF is zero after lag 2.
  • MA(3): ACF is zero after lag 3.
  • MA(4): ACF is zero after lag 4.

The sample ACF may show small non-zero values beyond these lags due to random sampling variation.

The estimated means may differ slightly from 3.3 because only 500 observations are used.

(iii) Repeat the Analysis with Sample Size 50000

The same four models are simulated again with 50,000 observations.

data_50000 <- lapply(
  coefficients,
  function(theta) generate_ma(50000, theta)
)

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

kable(
  results_50000,
  digits = 5,
  caption = "Estimated Process Means for Sample Size 50000"
)
Estimated Process Means for Sample Size 50000
Model True_Mean Sample_Mean
MA(1) MA(1) 3.3 3.29610
MA(2) MA(2) 3.3 3.30637
MA(3) MA(3) 3.3 3.30976
MA(4) MA(4) 3.3 3.32218

Compare the Estimated Means

mean_comparison <- data.frame(
  Model = names(data_500),
  Mean_n500 = sapply(data_500, mean),
  Mean_n50000 = sapply(data_50000, mean),
  True_Mean = 3.3
)

kable(
  mean_comparison,
  digits = 5,
  caption = "Comparison of Estimated Means"
)
Comparison of Estimated Means
Model Mean_n500 Mean_n50000 True_Mean
MA(1) MA(1) 3.36045 3.29610 3.3
MA(2) MA(2) 3.29166 3.30637 3.3
MA(3) MA(3) 3.39803 3.30976 3.3
MA(4) MA(4) 3.47407 3.32218 3.3

ACF Plots for Large Datasets

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

for (model in names(data_50000)) {

  acf(
    data_50000[[model]],
    lag.max = 20,
    main = paste("ACF of", model, "(n = 50000)")
  )
}

par(mfrow = c(1, 1))

Observation and Explanation:

The sample means obtained from 50,000 observations are generally closer to the true process mean of 3.3 than those obtained from 500 observations.

This happens because of the Law of Large Numbers. As the sample size increases, the sample mean converges towards the theoretical process mean.

The standard error of the sample mean decreases approximately at the rate:

\[ SE(\bar{Y})=\sqrt{\frac{\gamma(0)}{n}} \]

where \(\gamma(0)\) is the process variance and \(n\) is the sample size.

Thus, a larger dataset provides a more accurate estimate of the process mean.

(iv) Model Identification Using Sample ACF

The sample ACF is an important tool for identifying Moving Average models.

For an MA(q) process, the theoretical ACF becomes zero after lag q.

To perform model identification without relying on the original model labels, the datasets are randomly presented in a different order. Their time order is preserved.

set.seed(456)

shuffled_names <- sample(names(data_500))

for (i in seq_along(shuffled_names)) {

  model_data <- data_500[[shuffled_names[i]]]

  acf(
    model_data,
    lag.max = 20,
    main = paste("Unknown Series", i)
  )
}

Interpretation

The order of the MA model can be identified by observing where the ACF cuts off.

ACF Pattern Identified Model
Significant autocorrelation mainly at lag 1 MA(1)
Significant autocorrelations up to lag 2 MA(2)
Significant autocorrelations up to lag 3 MA(3)
Significant autocorrelations up to lag 4 MA(4)

Cross-Checking the Identification

The identified models can be compared with their actual generating models.

identification_check <- data.frame(
  Actual_Model = shuffled_names,
  Series_Name = paste("Unknown Series", seq_along(shuffled_names)),
  True_Order = sapply(
    coefficients[shuffled_names],
    length
  )
)

kable(
  identification_check,
  caption = "Cross-Checking the Original Model Orders"
)
Cross-Checking the Original Model Orders
Actual_Model Series_Name True_Order
MA(1) MA(1) Unknown Series 1 1
MA(4) MA(4) Unknown Series 2 4
MA(3) MA(3) Unknown Series 3 3
MA(2) MA(2) Unknown Series 4 2

Observation:

The ACF plots should show a cutoff pattern corresponding approximately to the order of each MA process. However, sample ACF values beyond the true order may occasionally cross the significance limits.

The identification is therefore based on the overall ACF pattern rather than requiring every individual lag to be insignificant.

Note: Randomly permuting the observations within a time series would destroy its serial dependence. Therefore, the series are only shuffled in presentation order, not internally rearranged.

(v) Change the White Noise Distribution to N(0,10)

Now, the variance of the underlying white noise is increased from 1 to 10.

Using the convention that the second parameter represents variance:

\[ \varepsilon_t\sim N(0,10) \]

Therefore, the standard deviation used in R is:

\[ \sigma=\sqrt{10} \]

The same four MA models are generated again.

set.seed(789)

data_high_variance <- lapply(
  coefficients,
  function(theta) {
    generate_ma(
      n = 500,
      theta = theta,
      mu = 3.3,
      noise_sd = sqrt(10)
    )
  }
)

Plot the New Datasets

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

for (model in names(data_high_variance)) {

  plot(
    data_high_variance[[model]],
    type = "l",
    main = paste(model, "with N(0,10) Noise"),
    xlab = "Time",
    ylab = expression(Y[t])
  )

  abline(h = 3.3, col = "red", lty = 2)
}

par(mfrow = c(1, 1))

Estimate Process Means

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

kable(
  high_variance_means,
  digits = 4,
  caption = "Estimated Means with Higher Noise Variance"
)
Estimated Means with Higher Noise Variance
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

ACF Plots with Higher Noise Variance

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

for (model in names(data_high_variance)) {

  acf(
    data_high_variance[[model]],
    lag.max = 20,
    main = paste("ACF of", model, "with N(0,10)")
  )
}

par(mfrow = c(1, 1))

Compare the Sample Means

variance_comparison <- data.frame(
  Model = names(data_500),
  Mean_N01 = sapply(data_500, mean),
  Mean_N010 = sapply(data_high_variance, mean),
  True_Mean = 3.3
)

kable(
  variance_comparison,
  digits = 4,
  caption = "Comparison of Process Means under Different Noise Variances"
)
Comparison of Process Means under Different Noise Variances
Model Mean_N01 Mean_N010 True_Mean
MA(1) MA(1) 3.3605 3.2914 3.3
MA(2) MA(2) 3.2917 3.2259 3.3
MA(3) MA(3) 3.3980 2.6959 3.3
MA(4) MA(4) 3.4741 3.5594 3.3

Observation and Explanation:

When the white noise variance increases from 1 to 10, the variability of the generated time series increases considerably.

However, the theoretical process mean remains unchanged at 3.3 because the white noise still has zero mean.

The theoretical autocorrelation function also remains unchanged because multiplying the white noise by a common standard deviation changes the autocovariances but not their normalized ratios.

Thus:

  • The time series becomes more variable.
  • The process mean remains 3.3.
  • The theoretical ACF remains the same.
  • The sample mean may show larger deviations from 3.3 for a fixed sample size.

1.4 Overall Conclusion

The four Moving Average models were successfully simulated using normally distributed white noise.

The following observations were obtained:

  1. All four MA processes are stationary and fluctuate around the common process mean of 3.3.
  2. The theoretical ACF of an MA(q) process cuts off after lag q.
  3. The sample mean becomes more accurate as the sample size increases from 500 to 50,000.
  4. The ACF plot can be used to identify the order of an MA model, although sample fluctuations may cause occasional deviations.
  5. Increasing the white noise variance increases the variability of the series but does not change its theoretical mean or ACF.

Hence, the simulation demonstrates the important properties of Moving Average processes and the use of ACF plots in time series model identification.