Εβδομαδιαία Εργασία 11:
Agapi Serasidou
2026-05-31
Johnson & Johnson Dataset:
Το dataset JohnsonJohnson περιλαμβάνει
84 τριμηνιαίες παρατηρήσεις των κερδών ανά μετοχή
(Earnings Per Share - EPS) της ομώνυμης εταιρείας, καλύπτοντας την
περίοδο από το 1ο τρίμηνο του 1960 έως και το 4ο τρίμηνο του
1980. Πρόκειται για μια κλασική χρονοσειρά (συχνότητας 4 ανά
έτος) που καταγράφει την οικονομική πορεία της εταιρείας σε δολάρια, με
τις τιμές να κυμαίνονται από $0.71 έως $16.20. Το σύνολο δεδομένων
παρουσιάζει ιδιαίτερο στατιστικό ενδιαφέρον καθώς συνδυάζει μια έντονη,
εκθετικά αυξητική μακροχρόνια τάση με ένα σταθερά
επαναλαμβανόμενο μοτίβο εποχικότητας μεταξύ των
τεσσάρων τριμήνων κάθε έτους.
#Εγκατάσταση & Φόρτωση πακέτων:
library(forecast)
library(tseries)
library(ggplot2)
library(gridExtra)
set.seed(42)
# --- Φόρτωση δεδομένων ---
data("JohnsonJohnson")
jj <- JohnsonJohnsonΜε την παρακάτω εντολή επιβεβαίωνουμε ότι τα δεδομένα είναι αποθηκευμένα ως ts (Time Series):
class(jj) # έλεγχος: είναι ήδη ts object;## [1] "ts"Για τον υπολογισμό των χρονικών ορίων και της συχνότητας, χρησιμοποιώ τις εξής εντολές:
start(jj); end(jj); frequency(jj)## [1] 1960 1## [1] 1980 4## [1] 4Πιο συγκεκριμένα: * start(jj) -> 1960 1: Τα δεδομένα ξεκινούν από το 1ο τρίμηνο του 1960.
end(jj) -> 1980 4: Τα δεδομένα σταματούν στο 4ο τρίμηνο του 1980.
frequency(jj) -> 4: Η συχνότητα είναι 4, που σημαίνει ότι έχουμε 4 παρατηρήσεις ανά έτος (τριμηνιαία δεδομένα).
Απο την παρακάτω εντολή, υπολογίζεται το πλήθος των δεδομένων του dataset που είναι 84.
length(jj) # 84 παρατηρήσεις (21 χρόνια × 4 τρίμηνα)## [1] 84Τέλος, δημιουργία ενός πίνακα που θα περίεχει τις πρώτες 12 τιμές του αρχικού πίνακα δεδομένων.
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Μεταβλητές του Dataset:
- Έτος (Year): Η χρονική μεταβλητή που ξεκινά από το 1960 και φτάνει μέχρι το 1980 (21 έτη συνολικά).
- Τρίμηνο (Quarter): Δείχνει σε ποια περίοδο του έτους αντιστοιχεί η κάθε μέτρηση (κάθε έτος χωρίζεται σε 4 τρίμηνα).
- Κέρδη ανά Μετοχή (EPS): Είναι η βασική αριθμητική τιμή (π.χ. 0.71, 0.63, 16.20). Δείχνει πόσα δολάρια κέρδους αντιστοιχούν σε κάθε μία μετοχή της εταιρείας για το συγκεκριμένο έτος και τρίμηνο.
# ============================================================
# ΜΕΡΟΣ Α — Exploration & Decomposition
# ============================================================
#TODO 1: Οπτικοποίηση της χρονοσειράς:
# Κύριο time series plot
autoplot(JohnsonJohnson) +
ggtitle("Τριμηνιαία κέρδη μετοχών της εταιρείας Johnson & Johnson (1960–1980)") +
xlab("Έτος") +
ylab("κέρδη") +
theme_minimal()
Το γράφημα δείχνει μια αυξητική τάση στσα κέρδη ανά μετοχή της Johnson
& Johnson, τα οποία από περίπου $0.71 το 1960, ξεπερνούν τα $15 το
1980.Παρατηρερίται επίσης, ένα σταθερό μοτίβο
εποχικότητας. Κάθε χρόνο, τα κέρδη της εταιρείας κάνουν
μια “κορυφή” στο 1ο και 2ο τρίμηνο, ενώ παρουσιάζουν μια σταθερή
“βουτιά” στο 3ο τρίμηνο . Αυτό το σκαμπανέβασμα επαναλαμβάνεται με τον
ίδιο ακριβώς τρόπο σε όλη τη διάρκεια των 21 ετών, δείχνοντας ότι η
κερδοφορία της εταιρείας επηρεάζεται συστηματικά από την εποχή του
χρόνου.. Επιπλέον, η διακύμανση των δεδομένων αυξάνεται με την πάροδο
του χρόνου, γεγονός που υποδηλώνει πολλαπλασιαστική δομή (multiplicative
time series).
# TODO 2: Εποχικός έλεγχος:
# Seasonal plot: κάθε έτος ώς χωριστή γραμμή
ggseasonplot(JohnsonJohnson, year.labels = TRUE) +
ggtitle("Seasonal plot — Johnson & Johnson")
Το διάγραμμα αυτό, δείχνει καθαρά τη συνεχή μακροχρόνια άνοδο των κερδών και επιβεβαιώνει ότι η αυξητική ταση κυριαρχεί σε σχέση με την εποχικότητα.
# Subseries plot: ένα mini-plot ανά μήνα
ggsubseriesplot(JohnsonJohnson) +
ggtitle("Subseries plot — μεταβολή ανά μήνα")
Το συγκεκριμένο διάγραμμα, απομονώνει κάθε τρίμηνο χχωριστά για να δείξει την εξέλιξή του μέσα στα 21 χρόνια. Οι μαύρες γραμμές δείχνουν ότι και τα τέσσερα τρίμηνα ακολουθούν την ίδια εκθετική ανάπτυξη, ενώ οι μπλε οριζόντιες γραμμές αποκαλύπτουν ότι το Q3 έχει τον υψηλότερο συνολικό μέσο όρο, με το Q4 να είναι το πιο αδύναμο.
TODO 3: Τύπος εποχικότητας:
Η χρονοσειρά έχει πολλαπλασιαστική (multiplicative) εποχικότητα. Αυτό σημαίνει ότι όσο μεγαλώνουν τα κέρδη της εταιρείας με τα χρόνια, τόσο μεγαλώνουν και τα σκαμπανέβασματά της. Στο ξεκίνημα (1960), που τα κέρδη είναι χαμηλά, οι διαφορές ανάμεσα στα τρίμηνα είναι πολύ μικρές. Προς το τέλος όμως (1980), που τα κέρδη έχουν εκτοξευθεί, οι ίδιες εποχικές διαφορές γίνονται τεράστιες, με πολύ ψηλές κορυφές και βαθιές βουτιές. Η εποχικότητα δηλαδή δεν προσθέτει ένα σταθερό ποσό κάθε χρόνο, αλλά λειτουργεί ως ποσοστό πάνω στα κέρδη της εκάστοτε περιόδου.
# TODO 4: Decomposition
dec_mult <- decompose(JohnsonJohnson, type = "multiplicative")
autoplot(dec_mult) +
ggtitle("Multiplicative Decomposition — Johnson & Johnson")
Trend Component: Το Trend Component δείχνει τη μακροχρόνια κατεύθυνση των κερδών της εταιρείας, έχοντας αφαιρέσει την εποχικότητα και τον τυχαίο θόρυβο.
Μορφή: Είναι καθαρά καμπυλωτό και όχι ευθύγραμμο.
Ερμηνεία: Τα κέρδη ανά μετοχή της Johnson & Johnson αυξάνονται με ολοένα και ταχύτερο ρυθμό όσο περνούν τα χρόνια. Συγκεκριμένα, από το 1960 έως το 1970 η άνοδος είναι σχετικά ήπια, ενώ από το 1970 έως το 1980 η καμπύλη παίρνει πολύ πιο απότομη κλίση προς τα πάνω.
# TODO 5: Log transformation
# Εφαρμογή λογαρίθμου
log_jj <- log(jj)
# Δημιουργία των δύο γραφημάτων
p1 <- autoplot(jj) + ggtitle("Αρχική Χρονοσειρά") + ylab("EPS")
p2 <- autoplot(log_jj) + ggtitle("Μετά το Λογάριθμο (Log)") + ylab("Log(EPS)")
# Εμφάνιση δίπλα-δίπλα
grid.arrange(p1, p2, ncol = 2)
Στα δύο γραφήματα φαίνεται ξεκάθαρα η επίδραση του λογαριθμικού μετασχηματισμού:
Αριστερό Γράφημα (Αρχική Χρονοσειρά): Η τάση είναι καμπυλωτή (εκθετική) και το μέγεθος των διακυμάνσεων μεγαλώνει όσο περνάει ο χρόνος (πολλαπλασιαστική δομή). Το 1960 οι αυξομειώσεις είναι ανεπαίσθητες, ενώ το 1980 γίνονται τεράστιες.
Δεξί Γράφημα (Μετά το Λογάριθμο): Η καμπύλη έχει ισιώσει και έχει γίνει μια σταθερή ευθεία γραμμή (γραμμική τάση). Το πιο σημαντικό είναι ότι το ύψος των σκαμπανεβασμάτων έχει πλέον το ίδιο σταθερό μέγεθος από την αρχή μέχρι το τέλος του γραφήματος.
# ============================================================
# ΜΕΡΟΣ Β — Modeling & Forecasting
# ============================================================
# TODO 6: Train/Test split — χρονολογικά!
# Train: 1960 Q1 — 1978 Q4 (76 παρατηρήσεις)
# Test: 1979 Q1 — 1980 Q4 ( 8 παρατηρήσεις = 2 χρόνια)
train <- window(JohnsonJohnson, end = c(1978, 4) )
test <- window(JohnsonJohnson, start = c(1979, 1))
length(train) ## [1] 76length(test) ## [1] 8# TODO 7: Έλεγχος στασιμότητας
adf.test(JohnsonJohnson)## Warning in adf.test(JohnsonJohnson): p-value greater than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: JohnsonJohnson
## Dickey-Fuller = 1.9321, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary# 2. Επειδή ΔΕΝ είναι στάσιμη, εφαρμόζουμε: Λογάριθμο -> Πρώτη Διαφορά -> Εποχική Διαφορά (lag=4)
# Σημείωση: Χρησιμοποιούμε το train set που φτιάξαμε στο TODO 6
log_train <- log(train)
stat_train <- diff(diff(log_train, differences = 1), lag = 4)
# 3. Επανέλεγχος στασιμότητας στη μετασχηματισμένη σειρά
adf.test(stat_train)## Warning in adf.test(stat_train): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: stat_train
## Dickey-Fuller = -6.3584, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryp-value = 0.01 :Η τιμή αυτή είναι πολύ μικρότερη από το κλασικό όριο του 0.05 (5%). Αυτό σημαίνει ότι απορρίπτουμε την αρχική υπόθεση ότι η σειρά έχει μοναδιαία ρίζα (δεν είναι στάσιμη).
alternative hypothesis: stationary: Μας επιβεβαιώνει την εναλλακτική υπόθεση, δηλαδή ότι η χρονοσειρά stat_train είναι πλέον στατιστικά στάσιμη.
Dickey-Fuller = -6.3584: Είναι η τιμή της στατιστικής ελέγχου. Όσο πιο αρνητική είναι (εν προκειμένω το -6.35 είναι αρκετά χαμηλό), τόσο πιο ισχυρή είναι η απόρριψη της μη-στασιμότητας.
# TODO 8: ACF & PACF
# Δύο plots δίπλα-δίπλα
p1 <- ggAcf(stat_train) + ggtitle("ACF")
p2 <- ggPacf(stat_train) + ggtitle("PACF")
grid.arrange(p1, p2, ncol = 2)
Διάγραμμα ACF: Στο διάγραμμα ACF της διαφορισμένης χρονοσειράς, οι τιμές που ξεπερνούν τις οριζόντιες μπλε διακεκομμένες γραμμές θεωρούνται στατιστικά σημαντικές. Παρατηρείται μια έντονα σημαντική αρνητική τιμή στο Lag 1, η οποία υποδηλώνει την ύπαρξη μιας ισχυρής μη-εποχικής εξάρτησης και μας κατευθύνει στην προσθήκη ενός όρου κινητού μέσου όρου πρώτης τάξης, MA(1). Επιπλέον, εμφανίζεται μια οριακά σημαντική αρνητική τιμή στο Lag 4 και μια θετική τιμή στο Lag 7, γεγονός που επιβεβαιώνει ότι υπάρχουν ακόμη υπολειμματικές εποχικές συσχετίσεις στα δεδομένα.
Διάγραμμα PACF: Στο διάγραμμα PACF, η στατιστική σημαντικότητα ελέγχεται με τον ίδιο τρόπο βάσει των μπλε ορίων. Εδώ ξεχωρίζει μια έντονα αρνητική τιμή στο Lag 1, η οποία αποτελεί κλασική ένδειξη για την ανάγκη εισαγωγής ενός αυτοπαλίνδρομου όρου πρώτης τάξης, AR(1). Το πιο αξιοσημείωτο εύρημα όμως είναι η πολύ έντονη αρνητική τιμή στο Lag 4. Επειδή το dataset αφορά τριμηνιαία δεδομένα (συχνότητα 4), η σημαντικότητα στο Lag 4 αποδεικνύει μια ισχυρή εποχική αυτοπαλίνδρομη εξάρτηση, υποδεικνύοντας την ανάγκη για έναν εποχικό όρο SAR(1).
#TODO 9: Εκπαιδεύστε ΤΡΙΑ μοντέλα στο train set, ορίζοντας το horizon = 8
# Πρόβλεψη: κάθε μήνας = ίδιος μήνας πέρυσι
fit_snaive <- snaive(train, h = 8)
autoplot(fit_snaive) +
ggtitle("Seasonal Naïve Forecast") +
autolayer(test, series = "Πραγματικά", colour = TRUE)
## Μοντέλο 2 — Holt-Winters
fit_hw <- hw(train, h = 8, seasonal = "multiplicative")
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.59970autoplot(fit_hw) +
ggtitle("Holt-Winters Forecast (h = 8)") +
autolayer(test, series = "Πραγματικά") +
scale_colour_manual(values = c("Πραγματικά" = "red", "Forecast" = "blue", "Σειρά" = "black"))
## Μοντέλο 3 — ARIMA (με auto.arima)
fit_arima <- auto.arima(train, lambda = 0)
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.1186641fc_arima <- forecast(fit_arima, h = 8)
autoplot(fc_arima) +
ggtitle("ARIMA Forecast (h = 8)") +
autolayer(test, series = "Πραγματικά") +
scale_colour_manual(values = c("Πραγματικά" = "red", "Forecast" = "blue", "Σειρά" = "black"))
# 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Με βάση τα αποτελέσματα της συνάρτησης για το βέτιστο μοντέλο ARIMA, προκύπτουν τα εξής συμπεράσματα:
Στατιστικός Έλεγχος Ljung-Box:Το τεστ Ljung-Box επιστρέφει p-value = 0.5788. Επειδή η τιμή αυτή είναι σημαντικά μεγαλύτερη από το κρίσιμο όριο του 0.05, αποτυγχάνουμε να απορρίψουμε την αρχική υπόθεση. Αυτό σημαίνει ότι δεν υπάρχει στατιστικά σημαντική αυτοσυσχέτιση στα κατάλοιπα του μοντέλου.
Γραφική Απεικόνιση (ACF & Ιστόγραμμα):
Στο διάγραμμα ACF, καμία από τις ράβδους των lags δεν ξεπερνάει τα μπλε διακεκομμένα όρια στατιστικής σημαντικότητας, επιβεβαιώνοντας την απουσία συσχετίσεων.Το ιστόγραμμα των καταλοίπων εμφανίζει μια ικανοποιητική κανονική κατανομή με μέσο όρο πολύ κοντά στο μηδέν.
Τα κατάλοιπα του μοντέλου συμπεριφέρονται ως λευκός θόρυβος. Αυτό αποδεικνύει ότι το ARIMA μοντέλο μας είναι απόλυτα επαρκές, καθώς κατάφερε να αποσπάσει όλη τη συστηματική πληροφορία (τάση και εποχικότητα) από τη χρονοσειρά, καθιστώντας τις προβλέψεις του αξιόπιστες.
# TODO 11: Σύγκριση προβλέψεων — ένα plot με όλα
# Όλες οι προβλέψεις σε ένα plot
autoplot(train, series = "Training") +
autolayer(test, series = "Test (Πραγματικά)") +
autolayer(fit_snaive$mean, series = "Seasonal Naïve") +
autolayer(fit_hw$mean, series = "Holt-Winters") +
autolayer(fc_arima$mean, series = "ARIMA") +
ggtitle("Σύγκριση μοντέλων στο test set") +
ylab("Επιβάτες") +
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)
# Απομόνωση των γραμμών που αφορούν το Test set (η 2η γραμμή κάθε αποτελέσματος)
metrics_table <- rbind(
"Seasonal Naive" = acc_snaive[2, c("RMSE", "MAE", "MAPE")],
"Holt-Winters" = acc_hw[2, c("RMSE", "MAE", "MAPE")],
"ARIMA" = acc_arima[2, c("RMSE", "MAE", "MAPE")]
)
# Εμφάνιση του πίνακα στην κονσόλα
print(metrics_table)## RMSE MAE MAPE
## Seasonal Naive 2.7765401 2.5425000 17.899110
## Holt-Winters 1.0865532 1.0427583 7.763851
## ARIMA 0.8136629 0.7225012 5.419212# Εμφάνιση σε όμορφο πίνακα για το R Markdown αναφοράς
knitr::kable(metrics_table, caption = "Σύγκριση Δεικτών Ακρίβειας στο Test Set", digits = 4)| RMSE | MAE | MAPE | |
|---|---|---|---|
| Seasonal Naive | 2.7765 | 2.5425 | 17.8991 |
| Holt-Winters | 1.0866 | 1.0428 | 7.7639 |
| ARIMA | 0.8137 | 0.7225 | 5.4192 |
Το μοντέλο ARIMA νικάει, καθώς επιτυγχάνει τις χαμηλότερες τιμές σε όλους τους δείκτες σφάλματος.
# TODO 13: Τελική πρόβλεψη για τον portfolio manager
# Τώρα που ξέρουμε το νικητή, εκπαιδεύουμε σε όλα τα δεδομένα
final_model <- auto.arima(JohnsonJohnson, lambda = 0)
final_forecast <- forecast(final_model, h = 36) # 3 χρόνια μπροστά
autoplot(final_forecast) +
ggtitle("Πρόβλεψη 1980-1983") +
ylab("Μετοχές (χιλιάδες)") +
theme_minimal()
Συνέχιση της Τάσης (Trend): Το μοντέλο αναγνωρίζει και επεκτείνει με επιτυχία την έντονη εκθετικά αυξητική τάση των προηγούμενων ετών. Η μπλε γραμμή των προβλέψεων συνεχίζει να κινείται με επιταχυνόμενο ρυθμό προς τα πάνω, ξεπερνώντας τα 60 κοντά στο 1983.
Διατήρηση της Εποχικότητας (Seasonality): Οι προβλέψεις διατηρούν το σταθερά επαναλαμβανόμενο κυκλικό μοτίβο με τις έντονες κορυφές και τις βαθιές βουτιές ανά έτος, αποδεικνύοντας ότι το μοντέλο έμαθε σωστά την τριμηνιαία συμπεριφορά της σειράς.
Πολλαπλασιαστική Συμπεριφορά: Το εύρος των διακυμάνσεων στις προβλέψεις μεγαλώνει όσο αυξάνεται το επίπεδο της σειράς (οι κορυφές γίνονται ψηλότερες και οι βουτιές βαθύτερες), ακολουθώντας την πολλαπλασιαστική φύση των ιστορικών δεδομένων.
Διαστήματα Εμπιστοσύνης (Μπλε Ζώνες): Οι σκούρες και ανοιχτές μπλε ζώνες γύρω από την κεντρική πρόβλεψη αντιπροσωπεύουν την αβεβαιότητα (διαστήματα εμπιστοσύνης 80% και 95%). Όσο ο ορίζοντας πρόβλεψης μακραίνει προς το 1983, η ζώνη ανοίγει και μεγαλώνει, γεγονός απόλυτα φυσιολογικό καθώς η αβεβαιότητα για το μέλλον αυξάνεται.