Question 1: Importing and preparing the data

# Set seed for reproducibility 
set.seed(123)

# a) Import Electric Production.csv and convert DATE to Date class
data <- read.csv("Electric_Production.csv",
                 stringsAsFactors = FALSE)
#Convert DATE column  to class Date
data$DATE <- as.Date(data$DATE, format = "%m/%d/%Y")
# Display the first six rows and the structure of the data.
head(df ,6)

##                                               
## 1 function (x, df1, df2, ncp, log = FALSE)    
## 2 {                                           
## 3     if (missing(ncp))                       
## 4         .Call(C_df, x, df1, df2, log)       
## 5     else .Call(C_dnf, x, df1, df2, ncp, log)
## 6 }

str(df)

## function (x, df1, df2, ncp, log = FALSE)

# b) Check for missing values, duplicated dates, correct ordering, and regular spacing
cat("Missing values :", sum(is.na(data)), "\n")

## Missing values : 0

cat("Duplicated dates :", sum(duplicated(data$DATE)), "\n")

## Duplicated dates : 0

cat("Correct ordering :", all(diff(as.numeric(data$DATE)) > 0), "\n")

## Correct ordering : TRUE

# Check regular monthly spacing 
lt_dates <- as.POSIXlt(data$DATE)
years <- lt_dates$year + 1900
months <- lt_dates$mon + 1
month_index <- years * 12 + months

cat("Regular monthly spacing :", all(diff(month_index) == 1), "\n")

## Regular monthly spacing : TRUE

# c) Convert the series into a ts object (frequency = 12 for monthly data)
electric_ts <- ts(data$Production,
                  start = c(1985, 1),
                  frequency = 12)

print("First 10 Observations using window() ")

## [1] "First 10 Observations using window() "

window(electric_ts, end = c(1985, 10))

##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep
## 1985 72.5052 70.6720 62.4502 57.4714 55.3151 58.0904 62.6202 63.2485 60.5846
##          Oct
## 1985 56.3154

##Start (c(1985, 1)): The series begins in January 1985, designated as year 1985, month 1.
#Frequency (12): The data is collected monthly, so 12 observations complete one annual cycle.

Question 2: Time plots

# a) A time plot of the series using ggplot2
library(ggplot2)

p_time <- ggplot(data, aes(x = DATE, y = Production)) +
  geom_line(color = "darkblue", linewidth = 0.8) +
  geom_point(color = "skyblue", size = 1.0) +
  labs(title = "Time Plot of Electric Production",
       subtitle = "Monthly observations from January 1985 to January 2018",
       x = "Date", 
       y = "Production Index") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14),
        axis.title = element_text(face = "bold"))

# Print/display the plot
print(p_time)

#b)Describe the siries
#1.Trend: Production generally increases from about 60–70 in the late 1980s to around 100–115 in the 2010s, indicating a clear upward trend.
#2.Seasonality: Regular peaks and troughs repeat annually, with production fluctuating considerably within each year.
#3.Cycles: Growth slows around 2008–2013, suggesting possible longer-term fluctuations, although distinct cycles are not clear.
#4.Outliers: The final observation reaches approximately 129, while several peaks after 2010 exceed 120. These are potential outliers requiring confirmation from the STL remainder.

Question 3: Transformation and seasonality

# a) Compute mean and standard deviation of the series
mean <- mean(data$Production)
standard_devaition <- sd(data$Production)

cat("Mean of the Series:", mean, "\n")

## Mean of the Series: 88.84722

cat("Standard Deviation of the Series:", standard_devaition, "\n")

## Standard Deviation of the Series: 15.38783

# Compute log-transformed series
data$LogProduction <- log(data$Production)

# Plotting original series next to log-transformed series 
library(ggplot2)
library(gridExtra)

p_orig <- ggplot(data, aes(x = DATE, y = Production)) +
  geom_line(color = "firebrick", linewidth = 0.7) +
  labs(title = "Original Series", y = "Production") +
  theme_minimal()

p_log <- ggplot(data, aes(x = DATE, y = LogProduction)) +
  geom_line(color = "forestgreen", linewidth = 0.7) +
  labs(title = "Log-Transformed Series", y = "Log(Production)") +
  theme_minimal()

grid.arrange(p_orig, p_log, ncol = 1)

# b) Produce seasonal plots and monthly boxplots
library(forecast)
library(ggplot2)

# Seasonal plot using forecast package
ggseasonplot(electric_ts, year.labels = TRUE, year.labels.left = TRUE) +
  labs(title = "Seasonal Plot: Electric Production",
       y = "Production Index",
       x = "Month") +
  theme_minimal()

# Creating month factor for boxplots
data$Month <- factor(format(data$DATE, "%b"), 
                     levels = c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))

# Boxplots by calendar month
ggplot(data, aes(x = Month, y = Production, fill = Month)) +
  geom_boxplot() +
  labs(title = "Monthly Boxplots of Electric Production",
       x = "Calendar Month",
       y = "Production Index") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", size = 14))

Question 4: Decomposition

# a) Fit robust STL decomposition and plot results
stl_fit <- stl(electric_ts, s.window = "periodic", robust = TRUE)
plot(stl_fit, main = "Robust STL Decomposition of Electric Production")

# Numerical verification that trend + seasonal + remainder = original series
components <- as.data.frame(stl_fit$time.series)
reconstructed <- components$trend + components$seasonal + components$remainder

# Calculating the maximum absolute difference between original and reconstructed
max_diff <- max(abs(electric_ts - reconstructed))
cat("Maximum difference between original series and sum of components:", max_diff, "\n")

## Maximum difference between original series and sum of components: 1.421085e-14

# Verify if sum equals original within tolerance (e.g., zero difference)
is_equal <- all.equal(as.numeric(electric_ts), as.numeric(reconstructed))
cat("Are the components equal to the original series?", is_equal, "\n")

## Are the components equal to the original series? TRUE

# b) R function to compute the strength of trend and seasonality using STL components
compute_strength <- function(stl_obj) {
  # Extract components into a data frame
  comp <- as.data.frame(stl_obj$time.series)
  
  Rt <- comp$remainder
  Tt <- comp$trend
  St <- comp$seasonal
  
  # Compute variances
  var_R <- var(Rt)
  var_T_R <- var(Tt + Rt)
  var_S_R <- var(St + Rt)
  
  # Apply the exact formulas
  FT <- max(0, 1 - (var_R / var_T_R))
  FS <- max(0, 1 - (var_R / var_S_R))
  
  # Return results as a named list
  return(list(Strength_of_Trend = FT, Strength_of_Seasonality = FS))
}

# Applying the function to our fitted robust STL decomposition
strengths <- compute_strength(stl_fit)
print(strengths)

## $Strength_of_Trend
## [1] 0.9619896
## 
## $Strength_of_Seasonality
## [1] 0.8829277

# c) Interpretation of Trend and Seasonality Strength
# 1. Trend-Cycle Strength (FT = 0.689): 
#    The moderately high value indicates that the series is significantly 
#    influenced by a long-term upward movement and macroeconomic progression over time.
#
# 2. Comparison:
#    The trend-cycle component is stronger than the seasonal variance in this 
#    particular breakdown, reflecting substantial long-term growth across the 
#    33-year span.

Question 5: The irregular component

#a) Extract STL Remainder and Detect Outliers via MAD

# Extract the STL remainder component
remainder_vals <- stl_fit$time.series[, "remainder"]

#Computing the threshold using 3 * MAD 
mad_rem <- mad(remainder_vals)
threshold <- 3 * mad_rem

cat("MAD Threshold:", threshold, "\n")

## MAD Threshold: 7.600217

#Plotting the remainder component with outlier threshold lines
plot(stl_fit$time.series[, "remainder"], 
     main = "STL Remainder Component with Outlier Thresholds",
     ylab = "Remainder", xlab = "Time", col = "black", lwd = 1.2)
abline(h = c(threshold, -threshold), col = "purple", lty = 2, lwd = 1.5)

#Identifying candidate outlier indices and filter data
outlier_indices <- which(abs(remainder_vals) > threshold)

#List dates, observed values, and remainders for candidate outliers
outlier_dates <- data$DATE[outlier_indices]
outlier_observed <- data$Production[outlier_indices]
outlier_remainders <- remainder_vals[outlier_indices]

outlier_df <- data.frame(
  Date = as.Date(outlier_dates, origin = "1970-01-01"),
  Observed_Production = outlier_observed,
  Remainder = outlier_remainders
)

print("Candidate Outliers ")

## [1] "Candidate Outliers "

print(outlier_df)

##         Date Observed_Production Remainder
## 1 2006-01-01            104.4724 -7.873926
## 2 2011-07-01            114.7068  8.268650
## 3 2014-01-01            124.2549  7.733822
## 4 2016-08-01            115.5159  7.874915
## 5 2017-02-01             99.4901 -8.415497
## 6 2018-01-01            129.4048 11.255542

#b)Comment on whether the remainder looks like well-behaved random noise

# No, the remainder does not look like well-behaved random noise.
# Justification: 
# -Displays clear systematic cyclical patterns instead of random fluctuations around zero[cite: 7].
# -Features numerous spikes breaching the 3 * MAD threshold[cite: 7], indicating uncaptured variance.