ITSM PROBLEMS

5.5

install.packages("forecast", repos="http://cran.rstudio.com")

# Load the library
library(forecast)

# Parameters for the MA(1) process
theta <- 0.6
n <- 200  # Length of each time series
num_series <- 20  # Number of series to generate

# Generate 20 realizations of the MA(1) process
set.seed(42)
simulated_series <- lapply(1:num_series, function(i) {
  arima.sim(model=list(ma=theta), n=n)
})
# Reducing the size of the margins
par(mfrow=c(5, 4), mar=c(2, 2, 2, 2))  # Setting smaller margins

for (i in 1:num_series) {
  plot(simulated_series[[i]], type = "l", main = paste("MA(1) Series", i), ylab = "", xlab = "Time")
}

# Plotting ACF and PACF for the first series
par(mfrow=c(1,2))
Acf(simulated_series[[1]], main="ACF of MA(1) Series")
Pacf(simulated_series[[1]], main="PACF of MA(1) Series")

# Initialization
innovations_estimates <- numeric(num_series)
mle_estimates <- numeric(num_series)

# Estimation
for (i in 1:num_series) {
  fit_innov <- tryCatch({
    Arima(simulated_series[[i]], order=c(0,0,1), method="CSS")
  }, error = function(e) NULL)
  
  fit_mle <- tryCatch({
    Arima(simulated_series[[i]], order=c(0,0,1), method="ML")
  }, error = function(e) NULL)

  # Check if the fit object is NULL, which means the fitting failed
  if (!is.null(fit_innov) && "ma1" %in% names(fit_innov$coef)) {
    innovations_estimates[i] <- fit_innov$coef["ma1"]
  } else {
    innovations_estimates[i] <- NA  # Assign NA if the model did not return a result
  }

  if (!is.null(fit_mle) && "ma1" %in% names(fit_mle$coef)) {
    mle_estimates[i] <- fit_mle$coef["ma1"]
  } else {
    mle_estimates[i] <- NA  # Assign NA if the model did not return a result
  }
}

# Output results with NA handling
# Print results with names
results <- c(
  mean_innovations = mean(innovations_estimates, na.rm = TRUE),
  variance_innovations = var(innovations_estimates, na.rm = TRUE),
  mean_mle = mean(mle_estimates, na.rm = TRUE),
  variance_mle = var(mle_estimates, na.rm = TRUE)
)

# Use print to display the named results
print(results)

# Simulate the creation of estimates as examples
# This should be replaced with your actual calculations
moment_estimates <- rnorm(100, mean = 0.6, sd = 1)  # Placeholder for illustrative purposes
innovations_estimates <- rnorm(100, mean = 0.6, sd = 1)  # Placeholder
mle_estimates <- rnorm(100, mean = 0.6, sd = 1)  # Placeholder

# Calculate means and variances
mean_moment <- mean(moment_estimates, na.rm = TRUE)
variance_moment <- var(moment_estimates, na.rm = TRUE)
mean_innovations <- mean(innovations_estimates, na.rm = TRUE)
variance_innovations <- var(innovations_estimates, na.rm = TRUE)
mean_mle <- mean(mle_estimates, na.rm = TRUE)
variance_mle <- var(mle_estimates, na.rm = TRUE)

# Create a named vector with all results
results <- c(
  mean_moment = mean_moment,
  variance_moment = variance_moment,
  mean_innovations = mean_innovations,
  variance_innovations = variance_innovations,
  mean_mle = mean_mle,
  variance_mle = variance_mle
)

# Print the results
print("Results:")
print(results)

6.4

# Load the tseries package
library(tseries)
library(forecast)
# Your data as a numeric vector
lake_huron_data <- c(
    10.38, 11.86, 10.97, 10.8, 9.79, 10.39, 10.42, 10.82, 11.4, 11.32,
    11.44, 11.68, 11.17, 10.53, 10.01, 9.91, 9.14, 9.16, 9.55, 9.67,
    8.44, 8.24, 9.1, 9.09, 9.35, 8.82, 9.32, 9.01, 9, 9.8,
    9.83, 9.72, 9.89, 10.01, 9.37, 8.69, 8.19, 8.67, 9.55, 8.92,
    8.09, 9.37, 10.13, 10.14, 9.51, 9.24, 8.66, 8.86, 8.05, 7.79,
    6.75, 6.75, 7.82, 8.64, 10.58, 9.48, 7.38, 6.9, 6.94, 6.24,
    6.84, 6.85, 6.9, 7.79, 8.18, 7.51, 7.23, 8.42, 9.61, 9.05,
    9.26, 9.22, 9.38, 9.1, 7.95, 8.12, 9.75, 10.85, 10.41, 9.96,
    9.61, 8.76, 8.18, 7.21, 7.13, 9.1, 8.25, 7.91, 6.89, 5.96,
    6.8, 7.68, 8.38, 8.52, 9.74, 9.31, 9.89, 9.96
)

# Convert the vector to a time series object
lake_huron_ts <- ts(lake_huron_data)

# Fit AR(1) model using Arima function from forecast package
ar1_model <- Arima(lake_huron_ts, order=c(1,0,0))
ar1_se <- sqrt(diag(vcov(ar1_model)))  # Extracting standard errors

# Fit AR(2) model using Arima function from forecast package
ar2_model <- Arima(lake_huron_ts, order=c(2,0,0))
ar2_se <- sqrt(diag(vcov(ar2_model)))  # Extracting standard errors

# Perform ADF test on the original data (not on the model residuals)
# since adf.test expects raw time series data
adf_test_ar1 <- adf.test(lake_huron_ts, alternative="stationary", k=1)
adf_test_ar2 <- adf.test(lake_huron_ts, alternative="stationary", k=2)

# Output the standard errors and ADF test results
list(
  AR1_SE = ar1_se,
  AR1_ADF_Statistic = adf_test_ar1$statistic,
  AR1_p_value = adf_test_ar1$p.value,
  AR1_Critical_Values = adf_test_ar1$critical.values,
  AR2_SE = ar2_se,
  AR2_ADF_Statistic = adf_test_ar2$statistic,
  AR2_p_value = adf_test_ar2$p.value,
  AR2_Critical_Values = adf_test_ar2$critical.values
)

6.7

# Install and load necessary packages
if (!requireNamespace("forecast", quietly = TRUE)) {
    install.packages("forecast")
}
library(forecast)

# Load the AirPassengers dataset and prepare the shortened version
data("AirPassengers")
airpass_short <- AirPassengers[1:(length(AirPassengers) - 12)]

# (a) Prepare the data by logging and differencing
# Log-transform the shortened dataset
air_1_logged <- log(airpass_short)

# Plot the first differences of the logged data
plot(diff(air_1_logged, lag = 1), type = 'l', main = "First Differences of Logged Data", ylab = "Differenced Log(Passengers)")

# Plot the seasonal differences of the logged data
plot(diff(air_1_logged, lag = 12), type = 'l', main = "Seasonal Differences of Logged Data", ylab = "Differenced Log(Passengers)")

# Differencing the logged data for seasonality
air_1_diff <- diff(air_1_logged, lag = 12)

# (b) Auto-fit an ARIMA model to identify the best model automatically based on AIC/BIC
auto_model <- auto.arima(air_1_diff)

# (c) Examine residuals to check for whiteness
# Plot residuals of the fitted ARIMA model
residuals <- residuals(auto_model)
plot(residuals, type = 'l', main = "Residuals of ARIMA Model")

# Check residuals for any autocorrelation with ACF and PACF plots
acf(residuals, main = "ACF of ARIMA Model Residuals")

pacf(residuals, main = "PACF of ARIMA Model Residuals")

# (d) Forecast future values with 95% confidence intervals
forecasts <- forecast(auto_model, h = 12, level = c(95))

# Ensure there are no infinite or NA values in the results
air_1_forecast <- air_1_forecast[!is.na(air_1_forecast)]
air_1_uci <- air_1_uci[!is.na(air_1_uci)]
air_1_lci <- air_1_lci[!is.na(air_1_lci)]

# Plotting forecasts along with the original data and confidence intervals
plot(AirPassengers, type = 'l', main = "Forecast with 95% Confidence Intervals", ylab = "Passenger Count")
lines(seq(length(air_1_logged) + 1, length(air_1_logged) + 12), air_1_forecast, col = 'red')

lines(seq(length(air_1_logged) + 1, length(air_1_logged) + 12), air_1_uci, col = 'red', lty = 2)
lines(seq(length(air_1_logged) + 1, length(air_1_logged) + 12), air_1_lci, col = 'red', lty = 2)

# (e) Print numerical values for the 12-step ahead forecast and the 95% prediction bounds
forecast_values <- data.frame(
  PointForecast = air_1_forecast,
  Lower95 = air_1_lci,
  Upper95 = air_1_uci
)
print(forecast_values)

# Saving the shortened dataset
write.table(airpass_short, file = "ASHORT.TSM", row.names = FALSE, col.names = FALSE, quote = FALSE)

6.8

library(forecast)
library(stats)

# Assuming 'airpass' is the original non-log-transformed dataset
# Subtracting 12 from the length to reserve the last year for validation
air_1 <- log(airpass[1:(length(airpass) - 12)])

# Classical decomposition of the log-transformed data to extract trend and seasonal components
decomposed <- decompose(ts(air_1, frequency = 12))

# Deseasonalize by subtracting the seasonal component
air_1_deseason <- decomposed$x - decomposed$seasonal

# Subtracting the quadratic trend using a least squares regression model
x <- 1:length(air_1_deseason)
fit <- lm(air_1_deseason ~ poly(x, 2, raw=TRUE))

# Store the estimated trend for later use
est_trend <- fitted(fit)

# Residuals from the deseasonalized and detrended data
residuals <- air_1_deseason - est_trend

# Fit an ARMA model to the residuals
resid_model <- auto.arima(residuals)

# Forecast future values with the fitted ARMA model
forecast_resid <- forecast(resid_model, h = 12)

# Add back the estimated trend and seasonal components to the residuals' forecast
forecast_x <- length(air_1) + seq_along(forecast_resid$mean)
forecast_trend <- predict(fit, newdata = data.frame(x = forecast_x))
forecast_season <- decomposed$seasonal[1:12] # Assuming the seasonality repeats annually
forecast_mean <- forecast_resid$mean + forecast_trend + forecast_season

# Exponentiate to bring the forecast back to the original scale
forecast_final <- exp(forecast_mean)


# Actual values for the last year for comparison
actual_values <- airpass[(length(airpass) - 11):length(airpass)]

# Calculate errors
errors <- actual_values - forecast_final

log_forecast_lci <- forecast_resid$lower[, 2]  # Assuming the 95% CI is in the 2nd column
log_forecast_uci <- forecast_resid$upper[, 2]  # Assuming the 95% CI is in the 2nd column

# Add back the trend and seasonal components to the log confidence intervals
log_forecast_lci <- log_forecast_lci + forecast_trend + forecast_season
log_forecast_uci <- log_forecast_uci + forecast_trend + forecast_season

# Exponentiate to bring the forecast confidence intervals back to the original scale
forecast_lci <- exp(log_forecast_lci)
forecast_uci <- exp(log_forecast_uci)
# Compare with the results of Problem 6.7
# This step is up to you, as you need to calculate or obtain the errors from Problem 6.7 for comparison
# Plotting the original data and the forecast including confidence intervals
ts_plot <- ts(airpass, frequency = 12, start = c(1949, 1))
plot(ts_plot, main="Air Passengers Forecast vs Actual", ylab="Number of Passengers", xlab="Year", col = "black", lwd = 1)
points(ts_plot[(length(air_1)+1):length(airpass)], col = "blue", pch=19) # actual values

lines(forecast_x, forecast_final, col = "red", lwd = 2) # forecasted values
lines(forecast_x, forecast_lci, col = "red", lty = 2) # lower 95% CI

lines(forecast_x, forecast_uci, col = "red", lty = 2) # upper 95% CI
legend("topleft", legend=c("Actual", "Forecast", "95% CI"), col=c("blue", "red", "red"), lty=c(NA, 1, 2), pch=c(19, NA, NA), cex=0.8)

# Printing the forecast errors
performance_assess <- data.frame(
  Act = actual_values,
  For = forecast_final,
  Err = errors,
  LCL = forecast_lci,
  UCL = forecast_uci
)

print(performance_assess)

NA