library(astsa)

sim <- 1000 # number of realizations to be simulated
T_values <- c( 50, 200, 500)
true_phi <- 0.9
true_theta <- 0.5

# Initiailizing the required matrices with empty values
phi_matrix <- matrix(NA, nrow = sim, ncol = 3)
theta_matrix <- matrix(NA, nrow = sim, ncol = 3)
sigma2_matrix <- matrix(NA, nrow = sim, ncol = 3)

# Looping through the 3 cases one by one
for( j in 1:3) {
  T = T_values[j]
  
  for( i in 1:sim) {
    # simulating ARMA(1,1) process
    x = arima.sim( n = T, list(ar = true_phi, ma = true_theta), sd = sqrt(1) )
    
    # Fitting ARIMA model
    fit = arima(x, order = c(1, 0, 1))
    
    # Storing MLE estimates
    phi_matrix[i, j] = fit$coef[1]
    theta_matrix[i, j] = fit$coef[2]
    sigma2_matrix[i ,j] = fit$sigma2
  }
}

# Computing mean square error, mean absolute deviation, and coverage of the 95% confidence intervals 
mse_phi = apply((phi_matrix - true_phi)^2, 2, mean)
mse_theta = apply((theta_matrix - true_theta)^2, 2, mean)
mse_sigma2 = apply((sigma2_matrix - 1)^2, 2, mean)

mad_phi = apply(abs(phi_matrix - true_phi), 2, mean)
mad_theta = apply(abs(theta_matrix - true_theta), 2, mean)
mad_sigma2 = apply(abs(sigma2_matrix - 1), 2, mean)

# Computing 95% confidence interval coverage
phi_coverage = rep(NA, 3)
theta_coverage = rep(NA, 3)
sigma2_coverage = rep(NA, 3)

for(j in 1:3) {
  se_phi = sd(phi_matrix[ , j])
  se_theta = sd(theta_matrix[ , j])
  
  # Coverage for phi
  phi_coverage[j] = mean(abs(phi_matrix[ , j] - true_phi) <= 1.96 * se_phi)
  
  # Coverage for theta 
  theta_coverage[j] = mean(abs(theta_matrix[ , j] - true_theta) <= 1.96 * se_theta)
  
  # For sigma2, we assign NA in case standard error is not available
  if(!is.null(sd(sigma2_matrix[ , j]))) {
    se_sigma2 = sd(sigma2_matrix[ , j])
    sigma2_coverage[j] = mean(abs(sigma2_matrix[ , j] - 1) <= 1.96 * se_sigma2)
  } else {
    sigma2_coverage[j] = NA
  }
}

# combining the above results into required defined matrices
mse_matrix = rbind(mse_phi, mse_theta, mse_sigma2)
mad_matrix = rbind(mad_phi, mad_theta, mad_sigma2)
coverage_matrix = rbind(phi_coverage, theta_coverage, sigma2_coverage)

# Displaying the required output
print("MSE Matrix: ")

## [1] "MSE Matrix: "

print(mse_matrix)

##                  [,1]        [,2]         [,3]
## mse_phi    0.01697729 0.001679379 0.0004918123
## mse_theta  0.02096685 0.004595710 0.0016932093
## mse_sigma2 0.04409070 0.010687945 0.0041272778

print("MAD Matrix: ")

## [1] "MAD Matrix: "

print(mad_matrix)

##                  [,1]       [,2]       [,3]
## mad_phi    0.09448723 0.03040329 0.01719714
## mad_theta  0.10998433 0.05451957 0.03256178
## mad_sigma2 0.17002784 0.08236150 0.05095508

print("Coverage Matrix: ")

## [1] "Coverage Matrix: "

print(coverage_matrix)

##                  [,1]  [,2]  [,3]
## phi_coverage    0.869 0.913 0.930
## theta_coverage  0.934 0.947 0.952
## sigma2_coverage 0.937 0.949 0.947