This is a problem I am trying to work on.

• Suppose I have some random data
• I fit a time series model to this data
• After fitting a model to this data, I want to simulate 100 points from this model

Here is what I tried.

First, simulate the data:

library(forecast)
library(ggplot2)
    set.seed(123)


n <- 1000  
random_data <- rnorm(n, mean = 0, sd = 1)
random_ts <- ts(random_data, frequency = 12)  

Next, I fit the model

    fit <- auto.arima(random_ts)

Series: random_ts 
ARIMA(0,0,0)(0,0,1)[12] with zero mean 

Coefficients:
         sma1
      -0.0792
s.e.   0.0305

sigma^2 = 0.9771:  log likelihood = -1406.89
AIC=2817.79   AICc=2817.8   BIC=2827.6

Here is where I am stuck: How do I simulate 100 trajectories from this model and record the average? Do I need arima.sim() for this?

I thought this could be done the same way as forecasting?

Thanks!

My own attempt:

library(forecast)
library(ggplot2)

set.seed(123)


fit <- auto.arima(random_ts)

simulate_forecast <- function(model, h) {
    sim <- simulate(model, nsim = h)
    return(as.numeric(sim))
}

n_sims <- 100
h <- 100
simulations <- replicate(n_sims, simulate_forecast(fit, h))

avg_prediction <- rowMeans(simulations)

time <- seq_len(n + h)
original_data <- c(as.numeric(random_ts), rep(NA, h))
sim_data <- rbind(matrix(NA, nrow = n, ncol = n_sims), simulations)
avg_data <- c(rep(NA, n), avg_prediction)

plot_data <- data.frame(
    time = rep(time, n_sims + 2),
    value = c(original_data, as.vector(sim_data), avg_data),
    type = factor(rep(c("Original", rep("Simulation", n_sims), "Average"), each = n + h))
)

 ggplot(plot_data, aes(x = time, y = value, group = interaction(type, rep(1:(n_sims+2), each = n + h)), color = type)) +
    geom_line(aes(alpha = type)) +
    scale_color_manual(values = c("Original" = "blue", "Simulation" = "red", "Average" = "black")) +
    scale_alpha_manual(values = c("Original" = 1, "Simulation" = 0.02, "Average" = 1)) +  #
    theme_minimal() +
    labs(title = "Original Data, 100 Simulations, and Average Prediction",
         x = "Time",
         y = "Value",
         color = "Type") +
    theme(legend.position = "bottom")