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 = 100 observed_data = ts(data[1:k]) # Fit model with fixed parameter values true_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 model true_model$sigma2 = 1 #---------------------------------------------------- # Forecast future states using perfect knowledge of model #---------------------------------------------------- n_forecast = 100 # number of points to be forecast forecast_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 interval lower_true <- forecast_true$lower[50,2] # [50,2] gives 95% # Test if true X_{k+50} lies in interval data[k+50] >= lower_true & data[k+50] <= upper_true # Illustrate forecast with plot # Create plot of forecast output using autoplot p_true<-autoplot(forecast_true,legend=TRUE) + ggtitle("Forecast - known coefficients") + xlab("Time") + ylab("Data points") #Ploting the forecast p_true #---------------------------------------------------- # Q4 is done at same time as Q5 #---------------------------------------------------- #---------------------------------------------------- # Forecast future states using perfect knowledge of model #---------------------------------------------------- # Fit the model to the data estimated_model <- arima(observed_data, order = c(2,1,2))# Fits a ARIMA(2,1,2) model to the data set coef<-as.array(estimated_model$coef) print(coef) # Print the estimated values of the coefficients # Forecast forecast_est= forecast(estimated_model, h=n_forecast) p_est<-autoplot(forecast_est, legend=TRUE) + ggtitle("Forecast - estimated coefficients") + xlab("Time") + ylab("Data points") p_est grid.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 interval lower_est <- forecast_est$lower[50,2] # [50,2] gives 95% # Test if true X_{k+50} lies in predicted interval data[k+50] >= lower_est & data[k+50] <= upper_est # Compare to interval for true model print(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 them sims <- list() # Store the simulations in a list true_accuracy <- 0 # Variables to count the proportion of simulations where est_accuracy <- 0 # simulated value lies in 95% prediction interval # Generate nsim simulated future paths for (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 interval if (sims[[i]][50] >= forecast_true$lower[50,2] & sims[[i]][50] <= forecast_true$upper[50,2]){ # add 1/nsim to proportion if so true_accuracy <- true_accuracy + 1/nsim } # Similarly for estimated if (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 intervals true_accuracy est_accuracy # Adding 10 of the simulated paths to forecast plots for (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 #----------------------------------------------------