5.5

# Install necessary library
install.packages("forecast")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/forecast_8.22.0.tgz'
Content type 'application/x-gzip' length 2484322 bytes (2.4 MB)
==================================================
downloaded 2.4 MB

The downloaded binary packages are in
    /var/folders/d_/mt52wgnd2rnbdnw9k1hp2t5r0000gn/T//RtmplbAQ3n/downloaded_packages

# Load the library
library(forecast)

Warning: package ‘forecast’ was built under R version 4.2.3Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
This is forecast 8.22.0 
  Want to stay up-to-date? Read the Hyndsight blog:
  https://robjhyndman.com/hyndsight/

# Parameters for the MA(1) process
theta <- 0.6
n <- 200  # Length of each time series
num_series <- 20  # Number of series to generate

# Generate 20 realizations of the MA(1) process
set.seed(42)
simulated_series <- lapply(1:num_series, function(i) {
  arima.sim(model=list(ma=theta), n=n)
})
# Reducing the size of the margins
par(mfrow=c(5, 4), mar=c(2, 2, 2, 2))  # Setting smaller margins

for (i in 1:num_series) {
  plot(simulated_series[[i]], type = "l", main = paste("MA(1) Series", i), ylab = "", xlab = "Time")
}

NA

# Plotting ACF and PACF for the first series
par(mfrow=c(1,2))
Acf(simulated_series[[1]], main="ACF of MA(1) Series")
Pacf(simulated_series[[1]], main="PACF of MA(1) Series")

# Initialization
innovations_estimates <- numeric(num_series)
mle_estimates <- numeric(num_series)

# Estimation
for (i in 1:num_series) {
  fit_innov <- tryCatch({
    Arima(simulated_series[[i]], order=c(0,0,1), method="CSS")
  }, error = function(e) NULL)
  
  fit_mle <- tryCatch({
    Arima(simulated_series[[i]], order=c(0,0,1), method="ML")
  }, error = function(e) NULL)

  # Check if the fit object is NULL, which means the fitting failed
  if (!is.null(fit_innov) && "ma1" %in% names(fit_innov$coef)) {
    innovations_estimates[i] <- fit_innov$coef["ma1"]
  } else {
    innovations_estimates[i] <- NA  # Assign NA if the model did not return a result
  }

  if (!is.null(fit_mle) && "ma1" %in% names(fit_mle$coef)) {
    mle_estimates[i] <- fit_mle$coef["ma1"]
  } else {
    mle_estimates[i] <- NA  # Assign NA if the model did not return a result
  }
}

# Output results with NA handling
# Print results with names
results <- c(
  mean_innovations = mean(innovations_estimates, na.rm = TRUE),
  variance_innovations = var(innovations_estimates, na.rm = TRUE),
  mean_mle = mean(mle_estimates, na.rm = TRUE),
  variance_mle = var(mle_estimates, na.rm = TRUE)
)

# Use print to display the named results
print(results)

    mean_innovations variance_innovations             mean_mle         variance_mle 
         0.590236635          0.003113181          0.592986705          0.003742935 

# Calculating and printing results with names
results <- c(
  mean_moment = mean(moment_estimates, na.rm = TRUE),
  variance_moment = var(moment_estimates, na.rm = TRUE),
  mean_innovations = mean(innovations_estimates, na.rm = TRUE),
  variance_innovations = var(innovations_estimates, na.rm = TRUE),
  mean_mle = mean(mle_estimates, na.rm = TRUE),
  variance_mle = var(mle_estimates, na.rm = TRUE)
)

# Print the results with names
print(results)

         mean_moment      variance_moment     mean_innovations variance_innovations 
         0.363182282          0.001352328          0.590236635          0.003113181 
            mean_mle         variance_mle 
         0.592986705          0.003742935