[SOLVED] CS # R Lab 3

# R Lab 3
# Forecasting usin an ARIMA(2,1,2) model

# Load in required packages
library(forecast)
library(ggplot2)
library(gridExtra) # Used to display plots in grid

#-
# Simulate and plot the data
#-

n = 1000
phi1=0.7
phi2=0.2
psi1=0.5
psi2=0.8
model <- list(ar = c(phi1, phi2), ma = c(psi1, psi2),order = c(2, 1, 2))data <- arima.sim(model = model, n = n)autoplot(data)#—————————————————-# Fit true model to first 100 datapoints#—————————————————-k = 100observed_data = ts(data[1:k])# Fit model with fixed parameter valuestrue_model <-arima(observed_data, order = c(2,1,2), fixed = c(phi1, phi2, psi1, psi2)) #Fitting an ARIMA(2,1,2) model for simulation# Set variance of modeltrue_model$sigma2 = 1#—————————————————-# Forecast future states using perfect knowledge of model#—————————————————-n_forecast = 100 # number of points to be forecastforecast_true = forecast(true_model, h = n_forecast) # Forecast using true_model# 95% prediction interval for X_{k+50}upper_true <- forecast_true$upper[50,2] # [50,1] gives 80% confidence intervallower_true <- forecast_true$lower[50,2] # [50,2] gives 95%# Test if true X_{k+50} lies in intervaldata[k+50] >= lower_true & data[k+50] <= upper_true# Illustrate forecast with plot# Create plot of forecast output using autoplotp_true<-autoplot(forecast_true,legend=TRUE) +ggtitle(“Forecast – known coefficients”) + xlab(“Time”) +ylab(“Data points”) #Ploting the forecastp_true#—————————————————-# Q4 is done at same time as Q5#—————————————————-#—————————————————-# Forecast future states using perfect knowledge of model#—————————————————-# Fit the model to the dataestimated_model <- arima(observed_data, order = c(2,1,2))# Fits a ARIMA(2,1,2) model to the data setcoef<-as.array(estimated_model$coef)print(coef) # Print the estimated values of the coefficients# Forecastforecast_est= forecast(estimated_model, h=n_forecast)p_est<-autoplot(forecast_est, legend=TRUE) +ggtitle(“Forecast – estimated coefficients”) +xlab(“Time”) +ylab(“Data points”)p_estgrid.arrange(p_true, p_est,nrow=2) # Predicted interval for X_{k+50}upper_est <- forecast_est$upper[50,2] # [50,1] gives 80% confidence intervallower_est <- forecast_est$lower[50,2] # [50,2] gives 95%# Test if true X_{k+50} lies in predicted intervaldata[k+50] >= lower_est & data[k+50] <= upper_est# Compare to interval for true modelprint(c(lower_est, upper_est))print(c(lower_true, upper_true))#—————————————————-# Investigate accuracy of two models using simulation#—————————————————-nsim = 10000 #Number of simulations- try changing themsims <- list() # Store the simulations in a listtrue_accuracy <- 0 # Variables to count the proportion of simulations whereest_accuracy <- 0# simulated value lies in 95% prediction interval# Generate nsim simulated future pathsfor (i in 1:nsim ){# Generate data for n_forecast=100 more timesteps sims[[i]]<-simulate(true_model,nsim=n_forecast)# Test if simulated value lies in true prediction intervalif (sims[[i]][50] >= forecast_true$lower[50,2] & sims[[i]][50] <= forecast_true$upper[50,2]){# add 1/nsim to proportion if sotrue_accuracy <- true_accuracy + 1/nsim}# Similarly for estimatedif (sims[[i]][50] >= forecast_est$lower[50,2] & sims[[i]][50] <= forecast_est$upper[50,2]){est_accuracy <- est_accuracy + 1/nsim}}# Alternative method# for (i in 1:nsim ){# sims[[i]]<-simulate(true_model,nsim=n_forecast)# ifelse(sims[[i]][50] >= forecast_true$lower[50,2] & sims[[i]][50] <= forecast_true$upper[50,2], true_accuracy <- true_accuracy + 1/nsim, true_accuracy <- true_accuracy)# ifelse(sims[[i]][50] >= forecast_est$lower[50,2] & sims[[i]][50] <= forecast_est$upper[50,2], est_accuracy <- est_accuracy + 1/nsim, est_accuracy <- est_accuracy)# }# }# Compare accuracy of prediction intervalstrue_accuracyest_accuracy# Adding 10 of the simulated paths to forecast plotsfor (i in 1:10){p_true <- p_true + autolayer(sims[[i]])p_est <- p_est + autolayer(sims[[i]])}grid.arrange(p_true, p_est,nrow=2) #—————————————————-# To alter value of k, return to top and change value#—————————————————-