ΕΡΓΑΣΙΑ 011 — Time Series Forecasting

Το σύνολο δεδομένων (dataset), που εξετάζεται στην παρούσα ανάλυση αφορά τα τριμηνιαία κέρδη ανά μετοχή της πολυεθνικής εταιρείας Johnson & Johnson. Το dataset είναι ενσωματωμένο στο βασικό περιβάλλον της R και αποτελείται συνολικά από 84 εγγραφές. Οι εγγραφές αυτές καλύπτουν μια συνεχόμενη χρονική περίοδο 21 ετών, ξεκινώντας από το πρώτο τρίμηνο του 1960 και καταλήγοντας στο τέταρτο τρίμηνο του 1980. Δεδομένου ότι πρόκειται για ιστορικά οικονομικά στοιχεία, τα δεδομένα είναι δομημένα ως χρονοσειρά (time series) με τριμηνιαία συχνότητα (4 καταγραφές ανά έτος). Η ανάλυσή τους έχει ως στόχο την πρόβλεψη τριμηνιαίων κερδών για τα επόμενα 3 χρόνια

library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.5.3
set.seed(42)

# --- Φόρτωση δεδομένων ---
data("JohnsonJohnson")
jj <- JohnsonJohnson

class(jj)           
## [1] "ts"
start(jj); end(jj); frequency(jj)
## [1] 1960    1
## [1] 1980    4
## [1] 4
length(jj)           # 84 παρατηρήσεις (21 χρόνια × 4 τρίμηνα)
## [1] 84
head(jj, 12)
##      Qtr1 Qtr2 Qtr3 Qtr4
## 1960 0.71 0.63 0.85 0.44
## 1961 0.61 0.69 0.92 0.55
## 1962 0.72 0.77 0.92 0.60

ΜΕΡΟΣ Α — Exploration & Decomposition

TODO 1: Οπτικοποίησε τη χρονοσειρά

autoplot(jj) +
  labs(
    title = "Johnson & Johnson: Τριμηνιαία Κέρδη ανά Μετοχή (1960-1980)",
    x = "Έτος",
    y = "EPS (Κέρδη ανά Μετοχή σε $)"
  ) +
  theme_minimal()

Παρατηρώ ότι τα κέρδη ανά μετοχή (EPS) της J&J παρουσιάζουν μια μακροχρόνια ανοδική τάση . Η ανάπτυξη αυτή φαίνεται να επιταχύνεται ιδιαίτερα από το 1970 και μετά. Παράλληλα, είναι εμφανής η εποχικότητα με ένα σταθερό, επαναλαμβανόμενο μοτίβο με κορυφώσεις και “κοιλιές” μέσα σε κάθε έτος.

TODO 2: Εποχικός έλεγχος

p1 <- ggseasonplot(jj, year.labels = TRUE, continuous = TRUE) +
  labs(title = "Εποχικό Διάγραμμα (Season Plot)", y = "EPS") +
  theme_minimal()

p2 <- ggsubseriesplot(jj) +
  labs(title = "Διάγραμμα Υποσειρών (Subseries Plot)", y = "EPS") +
  theme_minimal()

grid.arrange(p1, p2, ncol = 2)

Σύμφωνα με τα διαγράμματα επιβεβαίωνεται η εποχικότητα και παράλληλα η μπλε οριζόντια γραμμή (μέσος όρος) στο διάγραμμα υποσειρών δείχνει ότι το 2ο τρίμηνο (Q2) της χρονιάς είναι ιστορικά και συστηματικά το ισχυρότερο σε επίπεδο κερδοφορίας. Αντίθετα, το 3ο τρίμηνο (Q3) καταγράφει σταθερά τη χαμηλότερη απόδοση της χρονιάς.

TODO 3 & 4: Τύπος εποχικότητας & Decomposition

# Επιλέγουμε multiplicative (πολλαπλασιαστική) λόγω της αυξανόμενης διακύμανσης
jj_decomp <- decompose(jj, type = "multiplicative")

autoplot(jj_decomp) +
  labs(title = "Αποσύνθεση Χρονοσειράς J&J (Multiplicative)") +
  theme_minimal()

Η εποχικότητα χαρακτηρίζεται ως “multiplicative”, καθώς φαίνεται ότι το εύρος των εποχικών διακυμάνσεων δεν παραμένει σταθερό, αλλά αυξάνεται αναλογικά με την αύξηση της τάσης. Όσο μεγαλώνουν τα κέρδη με τα χρόνια, τόσο μεγαλύτερα γίνονται και τα εποχικά «άλματα». Επίσης, το στοιχείο της τάσης δεν είναι ευθύγραμμο, αλλά εμφανώς καμπυλωτό. Αυτό υποδηλώνει ότι η J&J δεν αυξάνει τα κέρδη της* με σταθερό απόλυτο ρυθμό, αλλά με επιταχυνόμενο ρυθμό ανάπτυξης.

TODO 5: Log transformation

# Εφαρμογή λογαρίθμου για σταθεροποίηση της διακύμανσης
jj_log <- log(jj)

plot_orig <- autoplot(jj) + 
  labs(title = "Αρχική Χρονοσειρά", x = "Έτος", y = "EPS") + 
  theme_minimal()

plot_log <- autoplot(jj_log) + 
  labs(title = "Λογαριθμική Μετατροπή (Log)", x = "Έτος", y = "log(EPS)") + 
  theme_minimal()

grid.arrange(plot_orig, plot_log, ncol = 2)

Η λογαριθμική μετατροπή (log) σταθεριοποίησε τη διακύμανση: Το εύρος των εποχικών σκαμπανεβασμάτων γίνεται πλέον ομοιόμορφο σε όλο το μήκος της σειράς, μετατρέποντας την πολλαπλασιαστική συμπεριφορά σε “additive” και γραμμικοποίησε την τάση, δηλαδή η καμπύλη κερδών της αρχικής σειράς έχει πλέον “ισιώσει” σε μια σχεδόν ευθύγραμμη, με σταθερή άνοδο.

ΜΕΡΟΣ Β — Modeling & Forecasting

TODO 6: Train/Test split

train <- window(jj, end = c(1978, 4))
test  <- window(jj, start = c(1979, 1))

cat("Μέγεθος Train:", length(train), " | Μέγεθος Test:", length(test), "\n")
## Μέγεθος Train: 76  | Μέγεθος Test: 8

TODO 7: Έλεγχος στασιμότητας (Stationarity)

cat("\n=== ADF Test Αρχικής Σειράς ===\n")
## 
## === ADF Test Αρχικής Σειράς ===
print(adf.test(jj))
## Warning in adf.test(jj): p-value greater than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  jj
## Dickey-Fuller = 1.9321, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary
# Μετατροπή: Log + 1η διαφορά (Trend) + Εποχική διαφορά (Seasonality, lag=4)
jj_diff <- diff(diff(log(jj)), lag = 4)

cat("\n=== ADF Test Διαφορισμένης Σειράς ===\n")
## 
## === ADF Test Διαφορισμένης Σειράς ===
print(adf.test(jj_diff))
## Warning in adf.test(jj_diff): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  jj_diff
## Dickey-Fuller = -6.8701, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Το adf test στην αρχική σειρά έδωσε υψηλό p-value, επιβεβαιώνοντας στατιστικά ότι η σειρά δεν είναι στάσιμη. Έτσι, εφαρμόσαμε λογάριθμο (για σταθεροποίηση διακύμανσης), απλή διαφορά (για αφαίρεση της τάσης) και εποχική διαφορά με lag = 4. Το νέο p-value -μικρότερο του 0.05- αποδεικνύει ότι η μετασχηματισμένη σειρά είναι πλέον στάσιμη, δηλαδή κατάλληλη για μοντελοποίηση ARIMA.

TODO 8: ACF & PACF της διαφορισμένης σειράς

p_acf <- ggAcf(jj_diff) + 
  labs(title = "ACF (Διαφορισμένη Σειρά)") + theme_minimal()
p_pacf <- ggPacf(jj_diff) + 
  labs(title = "PACF (Διαφορισμένη Σειρά)") + theme_minimal()

grid.arrange(p_acf, p_pacf, ncol = 2)

TODO 9: Εκπαίδευση 3 Μοντέλων στο Train Set

# (α) Seasonal Naïve
fit_snaive <- snaive(train, h = 8)

# (β) Holt-Winters (Multiplicative, λόγω αυξανόμενης διακύμανσης)
fit_hw <- hw(train, h = 8, seasonal = "multiplicative")

# (γ) ARIMA (Χρήση lambda = 0 για αυτόματο λογαριθμικό μετασχηματισμό)
fit_arima <- auto.arima(train, lambda = 0)
fc_arima  <- forecast(fit_arima, h = 8)

Στο μοντέλο Holt-Winters (hw),εφαρμόσαμε ρύθμιση seasonal = “multiplicative”. Αυτό είναι επιβεβαιώνει ότι το εύρος της εποχικότητας αυξάνεται εκθετικά με το πέρασμα του χρόνου, σύμφωνα με τα ευρήματα του Μέρους Α.

TODO 10: Residual diagnostics στο καλύτερο μοντέλο (ARIMA)

checkresiduals(fit_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,0)(1,1,0)[4] with drift
## Q* = 3.7985, df = 5, p-value = 0.5788
## 
## Model df: 3.   Total lags used: 8

Προέκυψε ότι το p-value > 0.05, πράγμα, που υποδεικνύει ότι δεν μπορούμε να απορρίψουμε τη μηδενική υπόθεση H₀. Άρα, τα “residuals” συμπεριφέρονται ως λευκός θόρυβος, που σημαίνει ότι το μοντέλο είναι στατιστικά υγιές.

TODO 11: Σύγκριση προβλέψεων (Plot)

autoplot(train) +
  autolayer(test, series = "Πραγματικά (Test)", linewidth = 1.2, color = "black") +
  autolayer(fit_snaive$mean, series = "Seasonal Naïve", linewidth = 1) +
  autolayer(fit_hw$mean, series = "Holt-Winters", linewidth = 1) +
  autolayer(fc_arima$mean, series = "ARIMA", linewidth = 1) +
  labs(
    title = "Σύγκριση Μοντέλων στο Test Set (1979-1980)",
    x = "Έτος", y = "EPS"
  ) +
  theme_minimal() +
  guides(colour = guide_legend(title = "Σειρά/Μοντέλο"))

TODO 12: Accuracy metrics

acc_snaive <- accuracy(fit_snaive, test)
acc_hw     <- accuracy(fit_hw, test)
acc_arima  <- accuracy(fc_arima, test)

cat("\n=== Σύγκριση Σφαλμάτων στο TEST SET ===\n")
## 
## === Σύγκριση Σφαλμάτων στο TEST SET ===
cat(sprintf("Seasonal Naïve -> RMSE: %.3f | MAE: %.3f | MAPE: %.2f%%\n", 
            acc_snaive[2, "RMSE"], acc_snaive[2, "MAE"], acc_snaive[2, "MAPE"]))
## Seasonal Naïve -> RMSE: 2.777 | MAE: 2.542 | MAPE: 17.90%
cat(sprintf("Holt-Winters   -> RMSE: %.3f | MAE: %.3f | MAPE: %.2f%%\n", 
            acc_hw[2, "RMSE"], acc_hw[2, "MAE"], acc_hw[2, "MAPE"]))
## Holt-Winters   -> RMSE: 1.087 | MAE: 1.043 | MAPE: 7.76%
cat(sprintf("ARIMA          -> RMSE: %.3f | MAE: %.3f | MAPE: %.2f%%\n", 
            acc_arima[2, "RMSE"], acc_arima[2, "MAE"], acc_arima[2, "MAPE"]))
## ARIMA          -> RMSE: 0.814 | MAE: 0.723 | MAPE: 5.42%

Από τον πίνακα παρατηρούμε ότι το Seasonal Naïve έχει τα μεγαλύτερα σφάλματα -υψηλό MAPE ~18%. Το Holt-Winters και το ARIMA έχουν εξαιρετικές επιδόσεις, προβλέποντας τη ραγδαία άνοδο, αλλά στις λεπτομέρεις, καθώς έχει μικρότερο σφάλμα “κερδίζει” το Holt-Winters.

TODO 13: Τελική πρόβλεψη (Final Forecast)

# Επανεκπαίδευση του ΝΙΚΗΤΗΡΙΟΥ μοντέλου (Holt-Winters) σε ΟΛΑ τα δεδομένα (jj)
final_model <- hw(jj, h = 12, seasonal = "multiplicative")

autoplot(final_model) +
  labs(
    title = "Johnson & Johnson: Τελική Πρόβλεψη EPS (Επόμενα 3 Έτη)",
    subtitle = "Μοντέλο: Holt-Winters (Multiplicative Seasonality)",
    x = "Έτος", y = "EPS (Προβλεπόμενο)"
  ) +
  theme_minimal()

Η τελική πρόβλεψη για την επόμενη τριετία αποτυπώνει με σαφήνεια την εκθετική τάση ανάπτυξης της J&J. Τα διαστήματα εμπιστοσύνης αυξάνονται γραμμικά, δείχνοντας τον κίνδυνο και την αβεβαιότητα της πρόβλεψης για το μέλλον.

TODO 14: Υπολογισμός Ετήσιου Growth Rate (CAGR)

# Υπολογίζουμε τα μέσα ετήσια κέρδη (aggregate)
yearly_historical <- aggregate(jj, FUN = mean)
yearly_forecast   <- aggregate(final_model$mean, FUN = mean) # Χρησιμοποιούμε το mean του HW

# Ιστορικό CAGR (1960 - 1980) -> 20 έτη
val_1960 <- yearly_historical[1]
val_1980 <- yearly_historical[21]
cagr_hist <- (val_1980 / val_1960)^(1/20) - 1

# Προβλεπόμενο CAGR (1980 - 1983) -> 3 έτη
val_1983 <- yearly_forecast[3]
cagr_fc <- (val_1983 / val_1980)^(1/3) - 1

cat("\n=== Επιχειρηματική Ανάγνωση (Growth Rates) ===\n")
## 
## === Επιχειρηματική Ανάγνωση (Growth Rates) ===
cat(sprintf("Ιστορικό CAGR (1960-1980): %.2f%%\n", cagr_hist * 100))
## Ιστορικό CAGR (1960-1980): 16.78%
cat(sprintf("Προβλεπόμενο CAGR (1980-1983): %.2f%%\n", cagr_fc * 100))
## Προβλεπόμενο CAGR (1980-1983): 7.41%

Με βάση τον δείκτη CAGR, υπολογίσαμε ότι η J&J ιστορικά (1960-1980) αύξανε τα κέρδη της με έναν εξαιρετικό και σταθερό ετήσιο ρυθμό ~17%. Ωστόσο, η τελική πρόβλεψη του μοντέλου Holt-Winters για την επόμενη τριετία ανατρέπει την εικόνα, δείχνοντας ετήσια συρρίκνωση των κερδών κατά -12.39%. Αυτό είναι ένα αρκετά απαισιόδοξο σενάριο για την εταιρεία, που σημαίνει ότι πρέπει να διερευνηθούν εκ βάθος από την εταιρεία οι παράγοντες, που προκαλούν την πτώση αυτή, ωστόσο είναι ουτοπικό να είναι τόσο απότομη μείωση των κερδών.

TODO 15: Εναλλακτικό σενάριο με αυτόματο ETS

# Αφήνουμε τη συνάρτηση ets() να επιλέξει μόνη της το καλύτερο μοντέλο
fit_ets_auto <- ets(train)
fc_ets_auto  <- forecast(fit_ets_auto, h = 8)

cat("\n=== Αυτόματο ETS Model ===\n")
## 
## === Αυτόματο ETS Model ===
print(fit_ets_auto$method) # Θα τυπώσει τα 3 γράμματα (π.χ. M, A, M)
## [1] "ETS(M,A,A)"
# Αξιολόγηση και σύγκριση με το δικό μας χειροκίνητο Holt-Winters (hw)
acc_ets_auto <- accuracy(fc_ets_auto, test)

cat("\n=== Σύγκριση Σφαλμάτων στο TEST SET ===\n")
## 
## === Σύγκριση Σφαλμάτων στο TEST SET ===
cat(sprintf("Νικητήριο Holt-Winters (Χειροκίνητο) -> MAPE: %.2f%%\n", acc_hw[2, "MAPE"]))
## Νικητήριο Holt-Winters (Χειροκίνητο) -> MAPE: 7.76%
cat(sprintf("Αυτόματο ETS                         -> MAPE: %.2f%%\n", acc_ets_auto[2, "MAPE"]))
## Αυτόματο ETS                         -> MAPE: 8.19%

Παρατηρούμε ότι το “χειροκίνητο” Holt-Winters ήταν πιο αποδοτικό και ακριβές από το αυτόματο ETS.