mod_ar <-Arima(WWWusage, order = c(2,0,0)) # AR(2) model
mod_ma <-Arima(WWWusage, order = c(0,0,2)) # MA(2) model
mod_arma <-Arima(WWWusage, order = c(1,0,1)) # ARMA(1,1) model
mod_arima <-auto.arima(WWWusage) # ARIMA model
mod_arima <-auto.arima(WWWusage)
# model summary
summary(mod_arima)
## Series: WWWusage
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.6504 0.5256
## s.e. 0.0842 0.0896
##
## sigma^2 = 9.995: log likelihood = -254.15
## AIC=514.3 AICc=514.55 BIC=522.08
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3035616 3.113754 2.405275 0.2805566 1.917463 0.5315228
## ACF1
## Training set -0.01715517
# model parameters
mod_arima$coef
## ar1 ma1
## 0.6503760 0.5255959
# 12-step-ahead forecast
forecast(mod_arima, h = 12)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 101 218.8805 214.8288 222.9322 212.6840 225.0770
## 102 218.1524 208.4496 227.8552 203.3133 232.9915
## 103 217.6789 202.3128 233.0449 194.1786 241.1792
## 104 217.3709 196.6302 238.1115 185.6508 249.0910
## 105 217.1706 191.4321 242.9091 177.8069 256.5343
## 106 217.0403 186.6844 247.3962 170.6149 263.4657
## 107 216.9556 182.3341 251.5771 164.0066 269.9046
## 108 216.9005 178.3266 255.4744 157.9068 275.8942
## 109 216.8646 174.6121 259.1172 152.2449 281.4844
## 110 216.8413 171.1478 262.5348 146.9592 286.7235
## 111 216.8262 167.8979 265.7544 141.9969 291.6555
## 112 216.8163 164.8325 268.8001 137.3140 296.3187
forecast(mod_arima, h = 12) %>% autoplot()
library(forecast)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.4.3
library(datasets)
# Load data
data("WWWusage")
www_data <- WWWusage
# Fit models (corrected syntax)
mod_ar <- Arima(WWWusage, order = c(2,0,0)) # AR(2) model
mod_ma <- Arima(WWWusage, order = c(0,0,2)) # MA(2) model
mod_arma <- Arima(WWWusage, order = c(1,0,1)) # ARMA(1,1) model
mod_arima <- auto.arima(WWWusage) # Auto ARIMA model
mod_arima_manual <- Arima(WWWusage, order = c(1,1,1)) # Manual ARIMA(1,1,1)
# Print model summaries
cat("=== MODEL SUMMARIES ===\n")
## === MODEL SUMMARIES ===
cat("\nAR(2) Model:\n")
##
## AR(2) Model:
print(summary(mod_ar))
## Series: WWWusage
## ARIMA(2,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 mean
## 1.8107 -0.8297 140.1117
## s.e. 0.0554 0.0571 16.3143
##
## sigma^2 = 11.47: log likelihood = -265.47
## AIC=538.94 AICc=539.36 BIC=549.36
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1177152 3.335831 2.722642 0.01575721 2.132211 0.6016552
## ACF1
## Training set 0.2008619
cat("\nMA(2) Model:\n")
##
## MA(2) Model:
print(summary(mod_ma))
## Series: WWWusage
## ARIMA(0,0,2) with non-zero mean
##
## Coefficients:
## ma1 ma2 mean
## 1.7426 0.9547 137.4306
## s.e. 0.0407 0.0427 4.2075
##
## sigma^2 = 136.1: log likelihood = -389.23
## AIC=786.47 AICc=786.89 BIC=796.89
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.1616367 11.48822 9.249899 -2.35186 7.465587 2.044063 0.8213117
cat("\nARMA(1,1) Model:\n")
##
## ARMA(1,1) Model:
print(summary(mod_arma))
## Series: WWWusage
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.9927 0.7984 149.3662
## s.e. 0.0089 0.0459 48.4678
##
## sigma^2 = 14.78: log likelihood = -278.24
## AIC=564.49 AICc=564.91 BIC=574.91
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.6382458 3.786451 3.002197 0.3772227 2.336051 0.6634319 0.4253221
cat("\nAuto ARIMA Model:\n")
##
## Auto ARIMA Model:
print(summary(mod_arima))
## Series: WWWusage
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.6504 0.5256
## s.e. 0.0842 0.0896
##
## sigma^2 = 9.995: log likelihood = -254.15
## AIC=514.3 AICc=514.55 BIC=522.08
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3035616 3.113754 2.405275 0.2805566 1.917463 0.5315228
## ACF1
## Training set -0.01715517
cat("\nManual ARIMA(1,1,1) Model:\n")
##
## Manual ARIMA(1,1,1) Model:
print(summary(mod_arima_manual))
## Series: WWWusage
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.6504 0.5256
## s.e. 0.0842 0.0896
##
## sigma^2 = 9.995: log likelihood = -254.15
## AIC=514.3 AICc=514.55 BIC=522.08
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3035616 3.113754 2.405275 0.2805566 1.917463 0.5315228
## ACF1
## Training set -0.01715517
# Generate forecasts
h <- 20
forecast_ar <- forecast(mod_ar, h = h)
forecast_ma <- forecast(mod_ma, h = h)
forecast_arma <- forecast(mod_arma, h = h)
forecast_arima <- forecast(mod_arima, h = h)
forecast_arima_manual <- forecast(mod_arima_manual, h = h)
# PLOT 1: Individual forecast plots
cat("\n=== CREATING FORECAST PLOTS ===\n")
##
## === CREATING FORECAST PLOTS ===
p1 <- autoplot(forecast_ar) +
ggtitle(paste("AR(2) Forecast\nAIC:", round(AIC(mod_ar), 2))) +
theme_minimal()
p2 <- autoplot(forecast_ma) +
ggtitle(paste("MA(2) Forecast\nAIC:", round(AIC(mod_ma), 2))) +
theme_minimal()
p3 <- autoplot(forecast_arma) +
ggtitle(paste("ARMA(1,1) Forecast\nAIC:", round(AIC(mod_arma), 2))) +
theme_minimal()
p4 <- autoplot(forecast_arima) +
ggtitle(paste("Auto ARIMA(",
paste(mod_arima$arma[1:3], collapse = ","),
") Forecast\nAIC:", round(AIC(mod_arima), 2))) +
theme_minimal()
p5 <- autoplot(forecast_arima_manual) +
ggtitle(paste("Manual ARIMA(1,1,1) Forecast\nAIC:", round(AIC(mod_arima_manual), 2))) +
theme_minimal()
# Arrange all forecast plots
grid.arrange(p1, p2, p3, p4, p5, ncol = 2)
# PLOT 2: Comparison of all forecasts
cat("\n=== FORECAST COMPARISON ===\n")
##
## === FORECAST COMPARISON ===
plot(forecast_ar, main = "WWWusage - All Model Forecasts Comparison",
ylab = "Number of Users", xlab = "Time")
lines(forecast_ma$mean, col = "blue", lwd = 2, lty = 2)
lines(forecast_arma$mean, col = "green", lwd = 2, lty = 2)
lines(forecast_arima$mean, col = "purple", lwd = 2, lty = 2)
lines(forecast_arima_manual$mean, col = "orange", lwd = 2, lty = 2)
legend("topleft",
legend = c("AR(2)", "MA(2)", "ARMA(1,1)", "Auto ARIMA", "ARIMA(1,1,1)"),
col = c("red", "blue", "green", "purple", "orange"),
lwd = 2, lty = c(1, 2, 2, 2, 2), bty = "n")
# PLOT 3: Fitted values comparison
cat("\n=== FITTED VALUES COMPARISON ===\n")
##
## === FITTED VALUES COMPARISON ===
par(mfrow = c(2, 3))
plot(WWWusage, main = "Original Data", ylab = "Users", type = "l", col = "black")
plot(fitted(mod_ar), main = "AR(2) Fitted", ylab = "Users", type = "l", col = "red")
plot(fitted(mod_ma), main = "MA(2) Fitted", ylab = "Users", type = "l", col = "blue")
plot(fitted(mod_arma), main = "ARMA(1,1) Fitted", ylab = "Users", type = "l", col = "green")
plot(fitted(mod_arima), main = "Auto ARIMA Fitted", ylab = "Users", type = "l", col = "purple")
plot(fitted(mod_arima_manual), main = "ARIMA(1,1,1) Fitted", ylab = "Users", type = "l", col = "orange")
par(mfrow = c(1, 1))
# PLOT 4: Residuals analysis
cat("\n=== RESIDUALS ANALYSIS ===\n")
##
## === RESIDUALS ANALYSIS ===
par(mfrow = c(3, 2))
# AR(2) residuals
res_ar <- residuals(mod_ar)
plot(res_ar, main = "AR(2) Residuals", type = "l", ylab = "Residuals")
abline(h = 0, col = "red")
acf(res_ar, main = "AR(2) Residuals ACF")
# MA(2) residuals
res_ma <- residuals(mod_ma)
plot(res_ma, main = "MA(2) Residuals", type = "l", ylab = "Residuals")
abline(h = 0, col = "red")
acf(res_ma, main = "MA(2) Residuals ACF")
# ARMA(1,1) residuals
res_arma <- residuals(mod_arma)
plot(res_arma, main = "ARMA(1,1) Residuals", type = "l", ylab = "Residuals")
abline(h = 0, col = "red")
acf(res_arma, main = "ARMA(1,1) Residuals ACF")
par(mfrow = c(1, 1))
# PLOT 5: Model coefficients comparison
cat("\n=== MODEL COEFFICIENTS ===\n")
##
## === MODEL COEFFICIENTS ===
models <- list(
"AR(2)" = mod_ar,
"MA(2)" = mod_ma,
"ARMA(1,1)" = mod_arma,
"Auto ARIMA" = mod_arima,
"ARIMA(1,1,1)" = mod_arima_manual
)
# Create a summary table
model_summary <- data.frame(
Model = names(models),
Order = c("(2,0,0)", "(0,0,2)", "(1,0,1)",
paste("(", mod_arima$arma[1], ",", mod_arima$arma[6], ",", mod_arima$arma[2], ")", sep = ""),
"(1,1,1)"),
AIC = sapply(models, AIC),
BIC = sapply(models, BIC),
RMSE = sapply(models, function(m) sqrt(mean(residuals(m)^2))),
Coefficients = sapply(models, function(m) {
coefs <- round(m$coef, 4)
paste(names(coefs), "=", coefs, collapse = ", ")
})
)
print(model_summary)
## Model Order AIC BIC RMSE
## AR(2) AR(2) (2,0,0) 538.9395 549.3602 3.335831
## MA(2) MA(2) (0,0,2) 786.4656 796.8863 11.488221
## ARMA(1,1) ARMA(1,1) (1,0,1) 564.4870 574.9077 3.786451
## Auto ARIMA Auto ARIMA (1,1,1) 514.2995 522.0848 3.113754
## ARIMA(1,1,1) ARIMA(1,1,1) (1,1,1) 514.2995 522.0848 3.113754
## Coefficients
## AR(2) ar1 = 1.8107, ar2 = -0.8297, intercept = 140.1117
## MA(2) ma1 = 1.7426, ma2 = 0.9547, intercept = 137.4306
## ARMA(1,1) ar1 = 0.9927, ma1 = 0.7984, intercept = 149.3662
## Auto ARIMA ar1 = 0.6504, ma1 = 0.5256
## ARIMA(1,1,1) ar1 = 0.6504, ma1 = 0.5256
# PLOT 6: Actual vs Fitted for all models
cat("\n=== ACTUAL VS FITTED VALUES ===\n")
##
## === ACTUAL VS FITTED VALUES ===
plot(WWWusage, type = "l", lwd = 2, col = "black",
main = "Actual vs Fitted Values - All Models",
ylab = "Number of Users", xlab = "Time")
lines(fitted(mod_ar), col = "red", lty = 2)
lines(fitted(mod_ma), col = "blue", lty = 2)
lines(fitted(mod_arma), col = "green", lty = 2)
lines(fitted(mod_arima), col = "purple", lty = 2)
lines(fitted(mod_arima_manual), col = "orange", lty = 2)
legend("topleft",
legend = c("Actual", "AR(2)", "MA(2)", "ARMA(1,1)", "Auto ARIMA", "ARIMA(1,1,1)"),
col = c("black", "red", "blue", "green", "purple", "orange"),
lwd = c(2, 1, 1, 1, 1, 1), lty = c(1, 2, 2, 2, 2, 2), bty = "n")
# Final model comparison
cat("\n=== FINAL MODEL RANKING BY AIC ===\n")
##
## === FINAL MODEL RANKING BY AIC ===
ranked_models <- model_summary[order(model_summary$AIC), ]
print(ranked_models)
## Model Order AIC BIC RMSE
## Auto ARIMA Auto ARIMA (1,1,1) 514.2995 522.0848 3.113754
## ARIMA(1,1,1) ARIMA(1,1,1) (1,1,1) 514.2995 522.0848 3.113754
## AR(2) AR(2) (2,0,0) 538.9395 549.3602 3.335831
## ARMA(1,1) ARMA(1,1) (1,0,1) 564.4870 574.9077 3.786451
## MA(2) MA(2) (0,0,2) 786.4656 796.8863 11.488221
## Coefficients
## Auto ARIMA ar1 = 0.6504, ma1 = 0.5256
## ARIMA(1,1,1) ar1 = 0.6504, ma1 = 0.5256
## AR(2) ar1 = 1.8107, ar2 = -0.8297, intercept = 140.1117
## ARMA(1,1) ar1 = 0.9927, ma1 = 0.7984, intercept = 149.3662
## MA(2) ma1 = 1.7426, ma2 = 0.9547, intercept = 137.4306
cat("\nBest model:", ranked_models$Model[1], "with AIC =", round(ranked_models$AIC[1], 2))
##
## Best model: Auto ARIMA with AIC = 514.3