Time Series Analysis Lab PS1

1 Problem 1: Moving Average (MA) Models

1.1 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.2 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.3 1.1 Function to generate an MA series

# Function to generate an MA(q) series
generate_ma <- function(n, theta, mu = 3.3, noise_sd = 1) {

  q <- length(theta)

  # Generate white-noise errors
  e <- rnorm(n + q, mean = 0, sd = noise_sd)

  y <- numeric(n)

  for (t in seq_len(n)) {

    y[t] <- mu + e[t + q]

    # Add the lagged error terms
    for (j in seq_len(q)) {
      y[t] <- y[t] + theta[j] * e[t + q - j]
    }
  }

  y
}

1.4 1.2 Generate four datasets of size 500

# Define the coefficients of the four MA models
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)
)

# Generate 500 observations from each model
ma_500 <- lapply(
  ma_coefficients,
  function(theta) generate_ma(500, theta)
)

# Save the generated datasets
write.csv(
  data.frame(Time = seq_len(500), ma_500),
  "PS1_Output/Tables/MA_Datasets_n500.csv",
  row.names = FALSE
)

1.5 (i) Plot the generated datasets

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

# Plot all four MA series
for (model in names(ma_500)) {

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

  # Red line indicates the theoretical mean
  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.

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

# Compare theoretical and sample means
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

# Save the mean estimates
write.csv(
  ma_mean_500,
  "PS1_Output/Tables/MA_Mean_n500.csv",
  row.names = FALSE
)

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

# Extract sample ACF values
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))

# Plot ACF for each MA model
for (model in names(ma_500)) {

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

par(old_par)

# Create a table of ACF values for lags 1 to 20
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

# Save the ACF table
write.csv(
  ma_acf_500,
  "PS1_Output/Tables/MA_ACF_n500.csv",
  row.names = FALSE
)

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.

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

# Generate 50,000 observations from each MA model
ma_50000 <- lapply(
  ma_coefficients,
  function(theta) generate_ma(50000, theta)
)

# Save the generated datasets
write.csv(
  data.frame(Time = seq_len(50000), ma_50000),
  "PS1_Output/Tables/MA_Datasets_n50000.csv",
  row.names = FALSE
)

# Compare means for both sample sizes
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

# Save the comparison table
write.csv(
  ma_mean_comparison,
  "PS1_Output/Tables/MA_Mean_Comparison.csv",
  row.names = FALSE
)

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

# Plot ACF for the larger datasets
for (model in names(ma_50000)) {

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

par(old_par)

# Calculate ACF values for the larger datasets
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

# Save the ACF table
write.csv(
  ma_acf_50000,
  "PS1_Output/Tables/MA_ACF_n50000.csv",
  row.names = FALSE
)

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.

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

# Randomly arrange the four series
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))

# Plot ACFs with anonymous labels
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.

# Reveal the actual models after visual identification
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

# Save the identification table
write.csv(
  identification_check,
  "PS1_Output/Tables/MA_Identification_Check.csv",
  row.names = FALSE
)

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.

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

# Generate MA series with white-noise variance 10
ma_high_variance <- lapply(
  ma_coefficients,
  function(theta) {

    generate_ma(
      n = 500,
      theta = theta,
      mu = 3.3,
      noise_sd = sqrt(10)
    )
  }
)

# Save the high-variance datasets
write.csv(
  data.frame(Time = seq_len(500), ma_high_variance),
  "PS1_Output/Tables/MA_Datasets_Variance10.csv",
  row.names = FALSE
)

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

# Plot the high-variance MA series
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)

# Estimate means for the high-variance datasets
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

# Save the mean estimates
write.csv(
  high_variance_means,
  "PS1_Output/Tables/MA_High_Variance_Means.csv",
  row.names = FALSE
)

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

# Plot ACF for high-variance datasets
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.10 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.

2 Problem 2: Autoregressive (AR) Models

2.1 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.2 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.3 2.2 Functions to generate AR data

# Generate an AR(1) process
generate_ar1 <- function(n, intercept, phi, burn = 200) {

  total <- n + burn
  e <- rnorm(total)
  y <- numeric(total)

  # Random initial value
  y[1] <- rnorm(1)

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

  # Return only observations after burn-in
  y[(burn + 1):total]
}

# Generate an AR(2) process
generate_ar2 <- function(n, phi1, phi2, burn = 200) {

  total <- n + burn
  e <- rnorm(total)
  y <- numeric(total)

  # Random initial values
  y[1:2] <- rnorm(2)

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

  # Return only observations after burn-in
  y[(burn + 1):total]
}

2.4 2.3 Generate four datasets of size 500

set.seed(2026)

# Generate the four AR datasets
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)
)

# Save all simulated AR datasets
write.csv(
  data.frame(Time = seq_len(500), ar_data),
  "PS1_Output/Tables/AR_Datasets_n500.csv",
  row.names = FALSE
)

2.5 (i) Plot the generated datasets

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

# Plot the four AR series
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.

2.6 (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.

# Compare sample means with theoretical means
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	10.0776
AR(1), phi = 1.01	AR(1), phi = 1.01	Non-stationary	Not defined as a stationary mean	63490.0850
AR(2), phi1 = 0.74, phi2 = -0.9	AR(2), phi1 = 0.74, phi2 = -0.9	Stationary	0	-0.0323
AR(2), phi1 = 0.74, phi2 = 0.9	AR(2), phi1 = 0.74, phi2 = 0.9	Non-stationary	Not defined as a stationary mean	-13472818230779573829828884044622260600464866820802266006468280644480400864008242804824468808620404.0000

# Save the mean estimates
write.csv(
  ar_mean_results,
  "PS1_Output/Tables/AR_Mean_Estimates.csv",
  row.names = FALSE
)

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.

2.7 (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 of model (a)
acf(
  ar_data[["AR(1), phi = 0.7"]],
  lag.max = 20,
  main = "Model (a): AR(1), phi = 0.7"
)

# ACF of model (b)
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.

2.8 (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 of model (c)
acf(
  ar_data[["AR(2), phi1 = 0.74, phi2 = -0.9"]],
  lag.max = 20,
  main = "Model (c): AR(2), phi2 = -0.9"
)

# ACF of model (d)
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.9 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.

3 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.