Time Series Analysis - Johnson & Johnson

Μέρος Α - Exploration & Decomposition

Πρώτη Οπτικοποίηση Χρονοσειράς

Σκοπός αυτού του μέρους είναι να αποκτήσουμε μία εικόνα των δεδομένων μας πριν εφαρμόσουμε οποιοδήποτε περίπλοκο μαθηματικό μοντέλο.

# Σχεδίαση της αρχικής σειράς
autoplot(jj, color = "#8B5E3C", linewidth = 1) +
  labs(title = "Johnson & Johnson Quarterly Earnings (1960-1980)",
       subtitle = "Ιστορική εξέλιξη των τriμηνιαίων κερδών ανά μετοχή (EPS)",
       x = "Έτος", y = "Κέρδη ανά Μετοχή (EPS)") +
  theme_fintech

Ευρήματα:

• Εξετάζοντας την ιστορική διαδρομή της J&J, διαπιστώνουμε μια ισχυρή, μακροχρόνια ανοδική τάση, η οποία δεν είναι γραμμική αλλά παρουσιάζει εκθετική μορφή (ανάπτυξη με ανατοκισμό).

• Εμφανής είναι μια έντονη εποχικότητα, με τα κέρδη να εμφανίζουν τακτικά σκαμπανεβάσματα εντός του έτους.

• Ακόμη πιο κρίσιμο είναι ότι το εύρος αυτών των διακυμάνσεων διευρύνεται με την πάροδο του χρόνου, ακολουθώντας την άνοδο της τάσης.
  

Εποχικός Έλεγχος

Ας εμβαθύνουμε τώρα στη διάσταση της εποχικότητας. Θέλουμε να απομονώσουμε τους μήνες ή τα τρίμηνα για να καταλάβουμε ποια περίοδος του έτους φέρνει συστηματικά το περισσότερο χρήμα στα ταμεία της εταιρείας.

p_season <- ggseasonplot(jj, year.labels = TRUE, year.labels.left = TRUE) +
  labs(title = "Εποχικό Διάγραμμα", x = "Τρίμηνο", y = "EPS") +
  theme_fintech

p_sub <- ggsubseriesplot(jj) +
  labs(title = "Διάγραμμα Υποσειρών", x = "Τρίμηνο", y = "Μέσο EPS") +
  theme_fintech

grid.arrange(p_season, p_sub, ncol = 2)

Ευρήματα:

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

• Η κορύφωση στο Q3 (Ιούλιος - Σεπτέμβριος) σχετίζεται άμεσα με τον προϊοντικό κύκλο των καταναλωτικών αγαθών υγείας και φροντίδας της J&J (consumer health). Οι λιανεμπορικές πωλήσεις προϊόντων που συνδέονται με το καλοκαίρι και την προστασία από τον ήλιο, καθώς και η αναπλήρωση αποθεμάτων από τις φαρμακαποθήκες ενόψει του φθινοπώρου, καταγράφονται και οριστικοποιούνται στα λογιστικά βιβλία κατά τη διάρκεια αυτού του τριμήνου. Η απότομη πτώση του Q4 αντανακλά τη φυσιολογική κάμψη της ζήτησης αυτών των κατηγοριών μετά το πέρας της θερινής περιόδου.

Ανάλυση Τύπου Εποχικότητας

Με βάση την οπτική ανάλυση των διαγραμμάτων του, στην περίπτωσή μας ισχύει ο πολλαπλασιαστικός τύπος (multiplicative seasonality) για τη συγκεκριμένη χρονοσειρά των δεδομένων μας.Συγκεκριμένα, το εύρος των εποχικών διακυμάνσεων (η απόσταση μεταξύ των κορυφών και των πυθμένων) δεν παραμένει σταθερό στο χρόνο, αλλά διευρύνεται αναλογικά με την άνοδο του επιπέδου της τάσης. Στις αρχές της περιόδου (κάτω γραμμές), όπου τα κέρδη είναι χαμηλά, οι εποχικές αυξομειώσεις είναι μικρές. Στα τελευταία έτη (πάνω γραμμές του 1979-1980), όπου η εταιρεία έχει αναπτυχθεί, οι εποχικές διακυμάνσεις γίνονται πολύ μεγαλύτερες. Αυτό σημαίνει ότι η εποχικότητα δεν προσθέτει ένα σταθερό ποσό κερδών κάθε χρόνο (additive), αλλά λειτουργεί ως σταθερό ποσοστό (multiplicative) επί της συνεχώς αυξανόμενης βάσης των κερδών.

Αποσύνθεση (Decomposition)

Σε αυτή τη φάση θα αποσυνθέσουμε τα επιμέρους στοιχεία που απαρτίζουν την χρονοσειρά μας, στηριζόμενοι στον πολλαπλασιαστικό τύπο που τεκμηριώσαμε παραπάνω:

# Εκτέλεση της πολλαπλασιαστικής αποσύνθεσης
fit_decomp <- decompose(jj, type = "multiplicative")

autoplot(fit_decomp) +
  labs(title = "Πολλαπλασιαστική Αποσύνθεση (Decomposition) της Χρονοσειράς J&J",
       x = "Έτος") +
  theme_fintech

Ευρήματα

• Η μακροχρόνια τάση της Johnson & Johnson είναι ξεκάθαρα καμπυλωτή (εκθετική) και όχι γραμμική, επιβεβαιώνοντας τη σταθερή ικανότητα της εταιρείας να επιτυγχάνει σύνθετη ετήσια ανάπτυξη.

• Ιδιαίτερο ενδιαφέρον παρουσιάζει η χρονική περίοδος μετά το 1970, όπου η καμπύλη παρουσιάζει μια έντονη αλλαγή κλίσης και επιταχύνεται εντυπωσιακά. Για τον portfolio manager, αυτό αποτελεί ισχυρή ένδειξη ότι η J&J κατά τη δεκαετία του ’70 πέρασε σε μια νέα φάση επιχειρηματικής κλίμακας, αυξάνοντας γεωμετρικά τα κέρδη της.

• Επιπλέον, το Seasonal Component επιβεβαιώνει τη σταθερότητα του τριμηνιαίου μοτίβου, δείχνοντας ότι η εποχική κορύφωση ενισχύει τα κέρδη κατά \(\sim 12\%\) πάνω από τη βασική τάση, ενώ η κάμψη του Q4 τα περιορίζει κατά \(\sim 15\%\).

Μετατροπή σε log()

Έχοντας αναγνωρίσει την εκθετική τάση και την αυξανόμενη διακύμανση της σειράς στα προηγούμενα βήματα, μεταβαίνουμε στον λογαριθμητικό μετασχηματισμό για να σταθεροποιήσουμε τα δεδομένα μας πριν από τη μοντελοποίηση.

# Γράφημα 1: Αρχικά Δεδομένα
p1 <- autoplot(jj, color = "#8B5E3C", linewidth = 0.8) + 
  labs(title = "Αρχικά Δεδομένα", 
       subtitle = "Εκθετική τάση & πολλαπλασιαστική διακύμανση", 
       x = "Έτος", y = "Κέρδη ανά Μετοχή (EPS)") + 
  theme_fintech

# Γράφημα 2: Λογαριθμημένα Δεδομένα
p2 <- autoplot(log(jj), color = "#5C4033", linewidth = 0.8) + 
  labs(title = "Λογαριθμημένα Δεδομένα [Log(J&J)]", 
       subtitle = "Γραμμικοποίηση τάσης & σταθερή διακύμανση", 
       x = "Έτος", y = "Λογάριθμος Κερδών [Log(EPS)]") + 
  theme_fintech

# Εμφάνιση των δύο γραφημάτων δίπλα-δίπλα σε 2 στήλες
grid.arrange(p1, p2, ncol = 2)

Ευρήματα:

• Στο αριστερό διάγραμμα, τα αρχικά δεδομένα χαρακτηρίζονται από εκθετική αύξηση και μια διακύμανση που διευρύνεται επικίνδυνα όσο αυξάνεται η αξία της μετοχής. Αυτή η μη-σταθερότητα στη διακύμανση καθιστά την απευθείας μοντελοποίηση αναξιόπιστη.

• Στο δεξί διάγραμμα, μετά την εφαρμογή του λογαρίθμου, η εκθετική τάση μετατρέπεται επιτυχώς σε γραμμική. Το σημαντικότερο εύρημα είναι η σταθεροποίηση της διακύμανσης, όπου το εύρος των εποχικών διακυμάνσεων παραμένει πλέον απόλυτα ομοιόμορφο σε όλο το ιστορικό εύρος (1960–1980), μετατρέποντας την εποχικότητα σε προσθετική (additive).

Μέρος Β - Modeling & Forecasting

Χρονολογική Διαίρεση του Dataset

• Κρατάμε το μεγαλύτερο μέρος της ιστορίας (1960 Q1 έως 1978 Q4 = 76 παρατηρήσεις) για να εκπαιδεύσουμε τα μοντέλα μας (Train set).

• Κρύβουμε τα τελευταία 2 χρόνια (1979 Q1 έως 1980 Q4 = 8 παρατηρήσεις) και τα ορίζουμε ως το πεδίο εξέτασης (Test set).

Το μοντέλο θα κάνει πρόβλεψη για αυτά τα 2 χρόνια και εμείς θα συγκρίνουμε τις προβλέψεις του με τα πραγματικά δεδομένα που κρύψαμε, για να δούμε ποιο μοντέλο είναι το πιο μάχιμο.

# 1. Χρονολογικός διαχωρισμός
train <- window(jj, end = c(1978, 4))
test  <- window(jj, start = c(1979, 1))

# 2. Δημιουργία δομημένου πίνακα για την αναφορά
split_summary <- data.frame(
  `Subset` = c("Train Set", "Test Set", "Total Dataset"),
  `Timeline` = c("1960 Q1 – 1978 Q4", "1979 Q1 – 1980 Q4", "1960 Q1 – 1980 Q4"),
  `Observations` = c(length(train), length(test), length(jj)),
  `Objective` = c("Εκπαίδευση μοντέλων & αναγνώριση μοτίβων", 
                  "Προσομοίωση & αξιολόγηση ακρίβειας", 
                  "Πλήρες ιστορικό εύρος")
)

# Εμφάνιση του πίνακα
knitr::kable(split_summary, align = "cccl")

Subset	Timeline	Observations	Objective
Train Set	1960 Q1 – 1978 Q4	76	Εκπαίδευση μοντέλων & αναγνώριση μοτίβων
Test Set	1979 Q1 – 1980 Q4	8	Προσομοίωση & αξιολόγηση ακρίβειας
Total Dataset	1960 Q1 – 1980 Q4	84	Πλήρες ιστορικό εύρος

Απομονώνουμε τις τελευταίες 8 παρατηρήσεις (έτη 1979–1980) ως Test set, ενώ οι προηγούμενες 76 παρατηρήσεις (1960–1978) αποτελούν το Train set. Η διατήρηση της χρονολογικής αλληλουχίας είναι κρίσιμη, καθώς διασφαλίζει ότι δεν θα υπάρξει διαρροή μελλοντικής πληροφορίας στα μοντέλα μας, προσομοιώνοντας ακριβώς τις συνθήκες κάτω από τις οποίες θα κληθεί να πάρει αποφάσεις ο portfolio manager στην πραγματική αγορά.

Έλεγχος Στασιμότητας

# Έλεγχος στασιμότητας στην αρχική σειρά (Train Set)
adf_original <- adf.test(train)

train_stationary <- diff(diff(log(train)), lag = 4)

adf_transformed <- adf.test(train_stationary)

# Μορφοποιημένος πίνακας αποτελεσμάτων
adf_summary <- data.frame(
  `Σειρά Δεδομένων` = c("Αρχικό Train Set", "Μετασχηματισμένο Train Set (Log + Diff + Seasonal Diff)"),
  `ADF Statistic` = c(round(adf_original$statistic, 3), round(adf_transformed$statistic, 3)),
  `p-value` = c(round(adf_original$p.value, 3), round(adf_transformed$p.value, 3)),
  `Κατάσταση` = c("Μη Στάσιμη (p > 0.05)", "Στάσιμη (p < 0.05)")
)

knitr::kable(adf_summary, align = "lccc")

Σειρά.Δεδομένων	ADF.Statistic	p.value	Κατάσταση
Αρχικό Train Set	0.853	0.99	Μη Στάσιμη (p > 0.05)
Μετασχηματισμένο Train Set (Log + Diff + Seasonal Diff)	-6.358	0.01	Στάσιμη (p < 0.05)

Παρατηρήσεις

• Το 0.99 στο αρχικό σετ επιβεβαιώνει με απόλυτη μαθηματική βεβαιότητα τη μη στασιμότητα.

• Το 0.01 (που στην R σημαίνει πρακτικά \(p \le 0.01\)) στο μετασχηματισμένο σετ δείχνει ότι η διαδικασία πέτυχε 100%.

Αυτοσυσχέτιση

# Παραγωγή των plots με το στυλ σου
p_acf  <- ggAcf(train_stationary) + 
  labs(title = "Συνάρτηση Αυτοσυσχέτισης (ACF)", x = "Lag", y = "ACF") +
  theme_fintech

p_pacf <- ggPacf(train_stationary) + 
  labs(title = "Μερική Συνάρτηση Αυτοσυσχέτισης (PACF)", x = "Lag", y = "PACF") +
  theme_fintech

# Εμφάνιση δίπλα-δίπλα
grid.arrange(p_acf, p_pacf, ncol = 2)

Ευρήματα:

• Βραχυπρόθεσμη εξάρτηση (Lag 1): Τόσο στο διάγραμμα ACF όσο και στο PACF, παρατηρείται ένα πολύ ισχυρό και στατιστικά σημαντικό αρνητικό spike στο Lag 1, το οποίο ξεπερνά καθαρά τα όρια σημαντικότητας (μπλε διακεκομμένες γραμμές). Αυτό υποδηλώνει έντονη συσχέτιση ανάμεσα σε διαδοχικά τρίμηνα και επιβεβαιώνει την ανάγκη εισαγωγής μη-εποχικών όρων AR(1) ή MA(1) στο μοντέλο.

• Εποχική εξάρτηση (Lag 4): Στο Lag 4 (που αντιπροσωπεύει την επίδραση του αντίστοιχου τριμήνου του προηγούμενου έτους) εντοπίζουμε επίσης σημαντικά αρνητικά spikes και στις δύο συναρτήσεις. Η παρουσία αυτών των spikes δείχνει ότι η διαφοροποίηση που εφαρμόσαμε δεν εξάλειψε πλήρως την εποχική εξάρτηση, γεγονός που καθιστά απαραίτητη την προσθήκη εποχικών παραμέτρων (Seasonal AR ή Seasonal MA) για την πλήρη εξομάλυνση των καταλοίπων.

• Αυτός ο συνδυασμός (σημαντικά lags τόσο χαμηλά όσο και στα πολλαπλάσια του 4) συνιστά ένδειξη ότι η σειρά χρειάζεται ένα Seasonal ARIMA μοντέλο.

Μοντέλα Εκπαίδευσης

Σε αυτό το στάδιο εκπαιδεύουμε τρεις διαφορετικές μεθοδολογίες πρόβλεψης στο Train set, ορίζοντας ορίζοντα πρόβλεψης 8 τριμήνων (h = 8), ώστε να καλύψουμε πλήρως την περίοδο του Test set.

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

# (β) Holt-Winters (Με πολλαπλασιαστική εποχικότητα βάσει του Μέρους Α)
fit_hw <- hw(train, h = 8, seasonal = "multiplicative")

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

# Εμφάνιση των summary() των τριών μοντέλων
summary(fit_snaive)

## 
## Forecast method: Seasonal naive method
## 
## Model Information:
## Call: snaive(y = train, h = 8) 
## 
## Residual sd: 0.8214 
## 
## Error measures:
##                     ME      RMSE       MAE      MPE     MAPE MASE      ACF1
## Training set 0.5884722 0.8214427 0.5920833 14.27463 14.77554    1 0.6642405
## 
## Forecasts:
##         Point Forecast     Lo 80     Hi 80     Lo 95    Hi 95
## 1979 Q1          11.88 10.827279 12.932721 10.270002 13.49000
## 1979 Q2          12.06 11.007279 13.112721 10.450002 13.67000
## 1979 Q3          12.15 11.097279 13.202721 10.540002 13.76000
## 1979 Q4           8.91  7.857279  9.962721  7.300002 10.52000
## 1980 Q1          11.88 10.391227 13.368773  9.603119 14.15688
## 1980 Q2          12.06 10.571227 13.548773  9.783119 14.33688
## 1980 Q3          12.15 10.661227 13.638773  9.873119 14.42688
## 1980 Q4           8.91  7.421227 10.398773  6.633119 11.18688

summary(fit_hw)

## 
## Forecast method: Holt-Winters' multiplicative method
## 
## Model Information:
## Holt-Winters' multiplicative method 
## 
## Call:
## hw(y = train, h = 8, seasonal = "multiplicative")
## 
##   Smoothing parameters:
##     alpha = 0.197 
##     beta  = 0.1082 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 0.5902 
##     b = 0.0072 
##     s = 0.8267 1.0468 1.0827 1.0439
## 
##   sigma:  0.1337
## 
##      AIC     AICc      BIC 
## 174.6172 177.3444 195.5938 
## 
## Error measures:
##                     ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.0420402 0.3595992 0.2602222 0.5285725 9.289055 0.4395027
##                   ACF1
## Training set 0.1282904
## 
## Forecasts:
##         Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 1979 Q1       12.83495 10.636472 15.03343  9.472669 16.19723
## 1979 Q2       13.76326 11.303396 16.22312 10.001224 17.52529
## 1979 Q3       13.74580 11.122117 16.36949  9.733221 17.75839
## 1979 Q4       11.19950  8.872544 13.52645  7.640730 14.75826
## 1980 Q1       14.57307 11.235132 17.91102  9.468132 19.67802
## 1980 Q2       15.56606 11.610776 19.52135  9.516974 21.61515
## 1980 Q3       15.48923 11.116368 19.86209  8.801514 22.17695
## 1980 Q4       12.57631  8.637829 16.51479  6.552923 18.59970

summary(fit_arima)

## Series: train 
## ARIMA(2,0,0)(1,1,0)[4] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1     ar2     sar1   drift
##       0.3131  0.2432  -0.2970  0.0386
## s.e.  0.1244  0.1300   0.1285  0.0046
## 
## sigma^2 = 0.008377:  log likelihood = 71.78
## AIC=-133.56   AICc=-132.65   BIC=-122.18
## 
## Training set error measures:
##                       ME      RMSE       MAE      MPE     MAPE      MASE
## Training set 0.001523642 0.3596867 0.2419193 -0.16328 6.770657 0.4085899
##                    ACF1
## Training set -0.1186641

Καλύτερη Επίδοση: ARIMA

Έχει το χαμηλότερο MAPE (6.77%), που σημαίνει ότι οι προβλέψεις του πέφτουν έξω κατά μέσο όρο μόνο 6.7%, έναντι 9.2% του Holt-Winters και 14.7% του Seasonal Naive. Επίσης, έχει το χαμηλότερο MAE (0.241). Το RMSE του είναι ουσιαστικά ολόιδιο με του Holt-Winters (0.359).

Residual Diagnostics

Έλεγχος Καταλοίπων για το ARIMA

# Έλεγχος καταλοίπων για το μοντέλο 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

Ευρήματα:

• Time Plot (Πάνω διάγραμμα): Τα κατάλοιπα διακυμαίνονται τυχαία γύρω από το μηδέν καθ’ όλη τη διάρκεια της περιόδου 1960–1978. Δεν παρατηρείται κάποια συστηματική μεταβολή στη διακύμανση ούτε κάποιο εμφανές κυκλικό μοτίβο, γεγονός που δείχνει σταθερότητα.

• ACF Plot (Κάτω αριστερά): Όλα τα κάθετα spikes για όλες τις χρονικές υστερήσεις (Lags) βρίσκονται αυστηρά μέσα στα μπλε διακεκομμένα όρια σημαντικότητας. Αυτό αποτελεί την ισχυρότερη οπτική απόδειξη ότι δεν έχει απομείνει καμία στατιστικά σημαντική αυτοσυσχέτιση στα σφάλματα.

• Ιστόγραμμα (Κάτω δεξιά): Η κατανομή των σφαλμάτων ακολουθεί με πολύ μεγάλη ακρίβεια την καμπύλη της κανονικής κατανομής (orange line), με μέση τιμή πολύ κοντά στο μηδέν, ικανοποιώντας τις βασικές υποθέσεις της γραμμικής μοντελοποίησης.

Συμπέρασμα Λευκού Θορύβου: Με βάση τα παραπάνω οπτικά ευρήματα και το Ljung-Box test (το οποίο παράγει \(p\text{-value} > 0.05\)), επιβεβαιώνεται ότι τα κατάλοιπα του μοντέλου \(ARIMA(2,0,0)(1,1,0)_4\) με drift αποτελούν Λευκό Θόρυβο (White Noise). Το μοντέλο έχει εξαγάγει το 100% της διαθέσιμης πληροφορίας από τα ιστορικά δεδομένα της J&J.

Σύγκριση Προβλέψεων

# Δημιουργία του κοινού γραφήματος σύγκρισης
autoplot(train, color = "#8B5E3C", linewidth = 0.8) +
  autolayer(test,            series = "Πραγματικά Δεδομένα (Test Set)", linewidth = 1) +
  autolayer(fit_snaive$mean, series = "Seasonal Naïve (Benchmark)",     linewidth = 0.8) +
  autolayer(fit_hw$mean,     series = "Holt-Winters (Multiplicative)",  linewidth = 0.8) +
  autolayer(fc_arima$mean,   series = "Auto ARIMA (2,0,0)(1,1,0)[4]",   linewidth = 0.8) +
  labs(title = "Σύγκριση Προβλέψεων Μοντέλων έναντι Πραγματικών Δεδομένων",
       subtitle = "Αξιολόγηση εκτός δείγματος (Out-of-sample Forecast) για τα έτη 1979-1980",
       x = "Έτος", y = "Κέρδη ανά Μετοχή (EPS)") +
  scale_color_manual(values = c(
    "Πραγματικά Δεδομένα (Test Set)" = "#222222", 
    "Seasonal Naïve (Benchmark)"     = "#D9534F", 
    "Holt-Winters (Multiplicative)"  = "#5CB85C", 
    "Auto ARIMA (2,0,0)(1,1,0)[4]"   = "#0275D8"
  )) +
  theme_fintech

Μέσα από το γράφημα της σύγκρισης, φαίνεται πεντακάθαρα πώς το ARIMA (μπλε γραμμή) και το Holt-Winters (πράσινη γραμμή) ακολουθούν την πραγματική πορεία (μαύρη γραμμή), ενώ το Seasonal Naïve (κόκκινη γραμμή) μένει πίσω επειδή δεν ακολουθεί την τάση.

Accuracy Metrics

# Υπολογισμός δεικτών ακρίβειας με βάση το πραγματικό test set
accuracy_snaive <- accuracy(fit_snaive, test)
accuracy_hw     <- accuracy(fit_hw, test)
accuracy_arima  <- accuracy(fc_arima, test)

# Απομόνωση των δεικτών για το Test Set (2η γραμμή των αποτελεσμάτων)
metrics_summary <- data.frame(
  `Μοντέλο` = c("Seasonal Naïve", "Holt-Winters", "Auto ARIMA"),
  `RMSE`    = c(accuracy_snaive[2, "RMSE"], accuracy_hw[2, "RMSE"], accuracy_arima[2, "RMSE"]),
  `MAE`     = c(accuracy_snaive[2, "MAE"],  accuracy_hw[2, "MAE"],  accuracy_arima[2, "MAE"]),
  `MAPE`    = c(accuracy_snaive[2, "MAPE"], accuracy_hw[2, "MAPE"], accuracy_arima[2, "MAPE"])
)

# Εμφάνιση του τελικού πίνακα
knitr::kable(metrics_summary, digits = 3, align = "lccc")

Μοντέλο	RMSE	MAE	MAPE
Seasonal Naïve	2.777	2.542	17.899
Holt-Winters	1.087	1.043	7.764
Auto ARIMA	0.814	0.723	5.419

Με MAPE 5.419% (έναντι 7.764% του Holt-Winters και 17.899% του Seasonal Naïve) και το χαμηλότερο RMSE (0.814), το ARIMA αναδεικνύεται ως το ικανότερο μοντέλο.

Τελική Πρόβλεψη

Σε αυτό το σημείο, επανεκπαιδεύουμε το βέλτιστο μοντέλο ARIMA στο πλήρες εύρος των δεδομένων και παράγουμε προβλέψεις για τα επόμενα 12 τρίμηνα (3 χρόνια), παρέχοντας στον portfolio manager τα διαστήματα εμπιστοσύνης.

# 1. Επανεκπαίδευση στο πλήρες dataset
final_arima_fit <- auto.arima(jj, lambda = 0)

# 2. Πρόβλεψη για 12 τρίμηνα (3 χρόνια) με επίπεδα εμπιστοσύνης 80% και 95%
final_forecast <- forecast(final_arima_fit, h = 12, level = c(80, 95))

# 3. Χειροκίνητη σχεδίαση για απόλυτο έλεγχο των χρωμάτων
# Σχεδιάζουμε την ιστορική σειρά και την κεντρική πρόβλεψη χωρίς το default shading
plot(final_forecast, shaded = FALSE, fcol = "#0275D8", flwd = 2,
     main = "Τελική Πρόβλεψη Κερδών Johnson & Johnson (1981-1983)",
     xlab = "Έτος", ylab = "Κέρδη ανά Μετοχή (EPS)")

# Προσθήκη ορίων 80% με Πράσινο χρώμα (στήλη 1 στο final_forecast$lower/upper)
lines(final_forecast$lower[,1], col = "green", lty = 2, lwd = 1.5)
lines(final_forecast$upper[,1], col = "green", lty = 2, lwd = 1.5)

# Προσθήκη ορίων 95% με Κόκκινο χρώμα (στήλη 2 στο final_forecast$lower/upper)
lines(final_forecast$lower[,2], col = "red", lty = 3, lwd = 1.5)
lines(final_forecast$upper[,2], col = "red", lty = 3, lwd = 1.5)

# Προσθήκη υπομνήματος για τον Portfolio Manager
legend("topleft", legend = c("Αναμενόμενη Πρόβλεψη", "Όρια Εμπιστοσύνης 80%", "Όρια Εμπιστοσύνης 95%"),
       col = c("#0275D8", "green", "red"), lty = c(1, 2, 3), lwd = 2, bty = "n")

Ευρήματα:

• Συνέχιση Τάσης: Η κεντρική πρόβλεψη (μπλε γραμμή) δείχνει σταθερή, ρεαλιστική άνοδο των κερδών, χωρίς ακραίες εξάρσεις, χάρη στο drift και τον λογάριθμο.

• Εποχικότητα: Το μοντέλο διατηρεί με ακρίβεια το ιστορικό μοτίβο, με τις κορυφές να εμφανίζονται σταθερά στο 3ο τρίμηνο (Q3) και τις πτώσεις στο 4ο τρίμηνο (Q4).

• Διαχείριση Ρίσκου: * Η πράσινη ζώνη (80%) δείχνει την πιο πιθανή πορεία των κερδών.

• Η κόκκινη ζώνη (95%) ορίζει τα ακραία σενάρια. Η κάτω κόκκινη γραμμή λειτουργεί ως «δίχτυ ασφαλείας» (worst-case scenario), δείχνοντας στον portfolio manager το ελάχιστο αναμενόμενο κέρδος ακόμη και σε κακές συνθήκες αγοράς.

Επιχειρηματική Ουσία

# 1. Απομόνωση των σημειακών προβλέψεων (Point Forecasts)
forecast_values <- final_forecast$mean

# 2. Υπολογισμός συνολικού ετήσιου EPS
eps_1980 <- sum(tail(jj, 4))          # Ιστορικό 1980
eps_1981 <- sum(forecast_values[1:4])   # Πρόβλεψη 1981
eps_1982 <- sum(forecast_values[5:8])   # Πρόβλεψη 1982
eps_1983 <- sum(forecast_values[9:12])  # Πρόβλεψη 1983

# 3. Υπολογισμός ετήσιων μεταβολών %
growth_81 <- (eps_1981 - eps_1980) / eps_1980 * 100
growth_82 <- (eps_1982 - eps_1981) / eps_1981 * 100
growth_83 <- (eps_1983 - eps_1982) / eps_1982 * 100

# Μέσος όρος προβλεπόμενης ανάπτυξης
mean_forecast_growth <- mean(c(growth_81, growth_82, growth_83))

cat("Ο μέσος ετήσιος ρυθμός ανάπτυξης των προβλέψεων είναι:", round(mean_forecast_growth, 2), "%\n")

## Ο μέσος ετήσιος ρυθμός ανάπτυξης των προβλέψεων είναι: 16.3 %

• Αξιολόγηση Ρεαλισμού: Ο προβλεπόμενος ρυθμός ανάπτυξης 16.3% κρίνεται απόλυτα ρεαλιστικός και συμβαδίζει με την εντυπωσιακή ιστορική πορεία της Johnson & Johnson κατά τις δεκαετίες ’60 και ’70.

• Financial Logic: Το ARIMA κατάφερε να αποτυπώσει την εκθετική τάση χωρίς να δημιουργήσει «φούσκα» (hyper-inflation). Λόγω του λογαρίθμου (lambda = 0) και του drift, η R μετέτρεψε την ανοδική τροχιά σε μια σταθερή, ισχυρή ετήσια μεγέθυνση της τάξης του 16.3%.

• Investment Decision: Για τον portfolio manager, αυτό το ποσοστό επιβεβαιώνει ότι η J&J παραμένει μια εξαιρετικά ελκυστική low-risk / high-yield επένδυση, προσφέροντας σταθερά διψήφιο ρυθμό ανάπτυξης για την επόμενη τριετία (1981-1983).

Σύγκριση με το ETS

# Εκπαίδευση αυτόματου μοντέλου ETS
fit_ets <- ets(train)
summary(fit_ets)

## ETS(M,A,A) 
## 
## Call:
## ets(y = train)
## 
##   Smoothing parameters:
##     alpha = 0.16 
##     beta  = 0.1018 
##     gamma = 0.3989 
## 
##   Initial states:
##     l = 0.6253 
##     b = 4e-04 
##     s = -0.1882 0.1782 -0.0081 0.0181
## 
##   sigma:  0.092
## 
##      AIC     AICc      BIC 
## 117.2365 119.9638 138.2131 
## 
## Training set error measures:
##                     ME     RMSE       MAE      MPE     MAPE      MASE
## Training set 0.0492013 0.418255 0.2614132 1.321354 7.079831 0.4415143
##                    ACF1
## Training set 0.02747311

Με βάση το επίσημο output, το μοντέλο που επιλέχθηκε αυτόματα είναι το ETS(M,A,A).

• Τι σημαίνουν τα 3 γράμματα ETS(M,A,A):

  • M (Error = Multiplicative): Το σφάλμα του μοντέλου είναι πολλαπλασιαστικό.

  • A (Trend = Additive): Η τάση της J&J αναγνωρίζεται ως προσθετική (γραμμική άνοδος).

  • A (Seasonality = Additive): Η εποχικότητα επιλέχθηκε ως προσθετική στο train set.

• Νικάει το Holt-Winters;

  • Όχι, το Holt-Winters κερδίζει οριακά. Αν συγκρίνουμε το RMSE στο Training set, το Holt-Winters έχει χαμηλότερο σφάλμα (0.359) έναντι του ETS (0.418).

  • Ωστόσο, το ETS(M,A,A) πέτυχε πολύ καλύτερο (χαμηλότερο) AICc (119.96) σε σχέση με το Holt-Winters (177.34), πράγμα που σημαίνει ότι στατιστικά είναι ένα πιο ισορροπημένο μοντέλο με λιγότερες παραμέτρους, παρόλο που το Holt-Winters έχει ελαφρώς καλύτερη ακρίβεια στα ιστορικά δεδομένα.