Στατιστική Ανάλυση Time Series Analysis
Εργασία 11
Kostas Soantzoglou
2026-05-26
Εισαγωγή: Time Series Forecasting
Dataset: Johnson & Johnson Quarterly Earnings (1960–1980) Στόχος: Πρόβλεψη τριμηνιαίων κερδών για τα επόμενα 3 χρόνια
Αρχικά, φορτώνουμε και ελέγχουμε τα δεδομένα μας.
# Φόρτωση δεδομένων
data("JohnsonJohnson")
jj <- JohnsonJohnson## 📊 ΒΑΣΙΚΑ ΧΑΡΑΚΤΗΡΙΣΤΙΚΑ ΧΡΟΝΟΣΕΙΡΑΣ## • Τύπος Αντικειμένου: ts## • Έναρξη Χρονοσειράς: Έτος 1960 - Τρίμηνο 1## • Λήξη Χρονοσειράς: Έτος 1980 - Τρίμηνο 4## • Συχνότητα Δεδομένων: 4 (Τρίμηνα ανά έτος)## • Σύνολο Παρατηρήσεων: 84 τρίμηνα (21 χρόνια)## (Πρώτες 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: Οπτικοποίησε τη χρονοσειρά
Προβάλλουμε τα δεδομένα για να αποκτήσουμε μια αρχική εικόνα της συμπεριφοράς τους.
# Κύριο time series plot
autoplot(JohnsonJohnson) +
ggtitle("Quarterly Revenue Johnson & Johnson") +
xlab("Έτος") +
ylab("Έσοδα") +
theme_minimal()
#εχουμε τάση και εποχικότητα το Q3 είναι το μεγαλύτερο πάντα σε εσοδα και το q4 το χειρότερο
2. TODO 2: Εποχικός Έλεγχος
Ποιο τρίμηνο της χρονιάς είναι συστηματικά το ισχυρότερο; Διερευνούμε ποιο τρίμηνο της χρονιάς είναι συστηματικά το ισχυρότερο χρησιμοποιώντας εποχιακά γραφήματα.
# Seasonal plot: κάθε έτος ως χωριστή γραμμή
ggseasonplot(JohnsonJohnson, year.labels = TRUE) +
ggtitle("Seasonal plot — Johnson & Johnson")
# Subseries plot: ένα mini-plot ανά μήνα/τρίμηνο
ggsubseriesplot(JohnsonJohnson) +
ggtitle("Subseries plot — μεταβολή ανά Τρίμηνο (Quarter)")
Παρατηρούμε πως το Q2 ειναι ανεβασμένο ωστοσο το Q3 φαινεται τα τελευταια χρονια να ειναι το καλύτερο για κάθε χρόνο
TODO 3: Τύπος εποχικότητας
Teκμηρίωσε γιατί θεωρείς ότι είναι multiplicative ή additive
Παρατηρουμε οτι δεν ειναι σταθερη η αυξηση,Έχουμε multipicative μοντελο επειδη μεγαλωνει αποτομα το πλάτος των κυμάτων και η εποχικότητα επηρεάζεται απο την τάση, Οσο αυξανεται το trend αυξανεται και το πλατος του κυματος. Δεν ειναι σταθερή η αυξηση εποχιακά. Δηλαδη ενα Q3 που ειναι το καλύτερο quarter κάθε χρονο δεν αυξάνεται σταθερά, και αυτο φαίνεται αισθητά στα τελευται 2-3 χρόνια Λογο του trend, εκτοξεύεται.
TODO 4: Decomposition
Σχολίασε: τι δείχνει το trend component; ευθύγραμμο ή καμπυλωτό;
# Επειδή είδαμε ότι το πλάτος μεγαλώνει → multiplicative
dec_mult <- decompose(JohnsonJohnson, type = "multiplicative")
autoplot(dec_mult) +
ggtitle("Multiplicative Decomposition — J&J")
Το trend δεν αυξανεται γραμμικά.
TODO 5. Log Transformation & Έλεγχος Συσχέτισης
Εφάρμοσε log() στη σειρά και κάνε autoplot() δίπλα-δίπλα με την αρχική (grid.arrange + ncol=2). Τι αλλάζει μετά το λογάριθμο;
data_log <- log(JohnsonJohnson)
stl_fit <- stl(data_log, s.window = "periodic")
autoplot(stl_fit) +
ggtitle("STL Decomposition (σε log scale)")
Εδώ παρατηρούμε αφού μετατρέψαμε την χρονοσειρά σε additive με την χρήση του log(), η τάση αυξάνεται γραμμικα αντι για καμπυλωτά Εφαρμόζουμε λογάριθμο (log) στη σειρά για να σταθεροποιήσουμε τη διακύμανση και έπειτα ελέγχουμε τα ACF/PACF γραφήματα. Το STL δουλεύει μόνο σε additive σειρές, οπότε η log μετατροπή είναι απαραίτητη.
# Δύο plots δίπλα-δίπλα (log transformed)
p1 <- ggAcf(data_log) + ggtitle("ACF (Log Scale)")
p2 <- ggPacf(data_log) + ggtitle("PACF (Log Scale)")
grid.arrange(p1, p2, ncol = 2)
Απο το PACF plot παρατηρουμε πως το πρώτο lag ειναι στατιστικά σημαντικό αρα μπορουμε να προβλέψουμε με βάση την αμέσως προηγούμενη περίοδο. Το 2ο,3ο κτλπ lags δεν είναι στατιστικά τόσο σημαντικά. #Για αυτο το ACF plot φθίνει γεωμετρικά. σημαινει οτι η κάθε περίδοος εξαρταται απο τη προηγούμενη περίοδο σε μεγάλο βαθμό.
Συγκριτικά, ας δούμε τα ACF/PACF στην αρχική σειρά (χωρίς λογάριθμο):
# Multiplied (Αρχικά Δεδομένα)
p1 <- ggAcf(JohnsonJohnson) + ggtitle("ACF (Original)")
p2 <- ggPacf(JohnsonJohnson) + ggtitle("PACF (Original)")
grid.arrange(p1, p2, ncol = 2)
εδω πέρα το PACF diagram μας λεει οτι εχει αξία και το 2ο lag (2 περιόδους πίσω απο την τωρινή περίοδο).
Ας κανουμε την σειρά στάσιμη εδω πέρα..
air_diff <- diff(diff(log(JohnsonJohnson)), lag = 4)
autoplot(air_diff) +
ggtitle("Μετά από log + 1η διαφόριση + εποχιακή διαφόριση")
1. 📈 Τι παρατηρείτε για την τάση της J&J; Είναι γραμμική ή εκθετική;
Η τάση είναι εκθετική, καθως δεν αυξάνεται γραμμικα και το επιβεβαιωσαμε απο το decomposition που κάναμε.
2. 🔁 Σε ποιο τρίμηνο εμφανίζουν οι κερδοφορίες την ισχυρότερη εποχική ώθηση; Δώστε μια πιθανή επιχειρηματική εξήγηση.
Το Q3 ειναι το πιο profitable κάθε χρόνο. Η εταιρεία μάλλον κάνει τις περισσότερες πωλήσεις εκείνη τη περίοδο, με αποτέλεσμα η αξία της μετοχής να ανεβαίνει.
3. ➕ Multiplicative ή additive; Τεκμηριώστε την απάντηση με βάση τα plots.
ΜΕΡΟΣ Β — Modeling & Forecasting
TODO 6. Train/Test Split
Χωρίζουμε χρονολογικά τα δεδομέναμας για την εκπαίδευση (Train) και την αξιολόγηση (Test) των μοντέλων.
# 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] 8TODO 7: Έλεγχος στασιμότητας
Εφάρμοσε adf.test() στη σειρά Εάν δεν είναι στάσιμη, εφάρμοσε log + diff + diff(lag=4) και ξαναδοκίμασε Προσοχή: εδώ η εποχική περίοδος είναι 4 (τρίμηνα), όχι 12!
adf.test(air_diff)##
## Augmented Dickey-Fuller Test
##
## data: air_diff
## Dickey-Fuller = -6.8701, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryH σειρά ειναι στάσιμη. Το καταλαβαίνουμε απο το p-value < 0.05
TODO 8: ACF & PACF
Φτιάξε ggAcf() και ggPacf() της διαφορισμένης σειράς δίπλα-δίπλα Σχολίασε: ποια lags δείχνουν σημαντικές τιμές;
p1 <- ggAcf(air_diff) + ggtitle("ACF")
p2 <- ggPacf(air_diff) + ggtitle("PACF")
grid.arrange(p1, p2, ncol = 2)
Αν παρατηρήσουμε το PACF plot είναι στατιστικά σημαντικά το Lag 1 και Lag 4. Αυτό σημαίνει: Lag 1 -> Η ακριβώς προηγούμενη περίοδος (ένα τρίμηνο πριν είναι στατιστικά σημαντική) Lag 4 -> 4 περιόδους πριν δηλαδή ένα χρόνο πριν το ίδιο τρίμηνο είναι στατιστικά σημαντικό, άρα εποχικότητα;
TODO 9: Εκπαιδεύστε ΤΡΙΑ μοντέλα στο train set, ορίζοντας το horizon = 8
(α) Seasonal Naïve → snaive(train, h = 8) (β) Holt-Winters → hw(train, h = 8, seasonal = “???”) (γ) ARIMA → auto.arima(train, lambda = 0); forecast(…, h = 8) Εμφάνισε τα summary() των τριών μοντέλων
Πρόβλεψη: κάθε μήνας = ίδιος μήνας πέρυσι
fit_snaive <- snaive(train, h = 8)
autoplot(fit_snaive) +
ggtitle("Seasonal Naïve Forecast") +
autolayer(test, series = "Πραγματικά", color = "red")
Τρεις παράμετροι: α (επίπεδο), β (τάση), γ (εποχικότητα)
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") +
autolayer(test, series = "Πραγματικά", color = "red")
Η auto.arima δοκιμάζει αυτόματα συνδυασμούς (p,d,q)(P,D,Q) και επιλέγει αυτόν με το χαμηλότερο AIC
fit_arima <- auto.arima(train, lambda = 0) # lambda=0 → log transformation
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") +
autolayer(test, series = "Πραγματικά", color = "red")
TODO 10: Residual diagnostics
Εφάρμοσε checkresiduals() στο καλύτερό σου μοντέλο. Είναι τα residuals white noise; (Ljung-Box p-value > 0.05;)
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- Ναι το σφαλμα οφειλεται σε white noise (p > 0.05), βλεπουμε κανονικη κατανομη στο ιστογραμμα (bell-shaped),
- και τo ACF plot μας λεει οτι το σφάλμα σημερα δεν επηρεαζσεται απο καποιο παρελθοντικο σφάλμα (lag).
- Η χρονοσειρά με τα residuals περναει κοντα απο την τιμη 0.0 και φαινεται σαν τυχαιο white noise αντι για καποιο pattern.
TODO 11: Σύγκριση προβλέψεων — ένα plot με όλα
autoplot(JohnsonJohnson) +
autolayer(test, series = "Πραγματικά",PI = FALSE) +
autolayer(fit_snaive$mean, series = "Seasonal Naïve",PI = FALSE) +
autolayer(fit_hw$mean, series = "Holt-Winters",PI = FALSE) +
autolayer(fc_arima$mean, series = "ARIMA")
Παρατηρουμε πως το arima ειναι το πιo ακριβης μοντελο μας και μετα ερχεται το holt winters.
TODO 12: Accuracy metrics
Εφάρμοσε accuracy() σε κάθε μοντέλο με το test set Συγκέντρωσε τα RMSE, MAE, MAPE σε ένα πίνακα Ποιο μοντέλο νικά;
accuracy(fit_snaive, test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.5884722 0.8214427 0.5920833 14.27463 14.77554 1.000000
## Test set 2.5425000 2.7765401 2.5425000 17.89911 17.89911 4.294159
## ACF1 Theil's U
## Training set 0.6642405 NA
## Test set -0.2050122 0.7697759accuracy(fit_hw, test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.04204020 0.3595992 0.2602222 0.5285725 9.289055 0.4395027
## Test set 0.07397684 1.0865532 1.0427583 -0.4202121 7.763851 1.7611681
## ACF1 Theil's U
## Training set 0.1282904 NA
## Test set -0.8122651 0.290191accuracy(fc_arima, test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.001523642 0.3596867 0.2419193 -0.163280 6.770657 0.4085899
## Test set -0.142883118 0.8136629 0.7225012 -1.524645 5.419212 1.2202695
## ACF1 Theil's U
## Training set -0.1186641 NA
## Test set -0.6150056 0.1989989To arima ειναι το μοντελο με το μικροτερο σφάλμα
TODO 13: Τελική πρόβλεψη για τον portfolio manager
Επανεκπαίδευσε το νικητήριο μοντέλο σε όλα τα διαθέσιμα δεδομένα (jj) Κάνε forecast(h = 12) για 12 τρίμηνα = 3 χρόνια μπροστά Οπτικοποίησε με prediction intervals διακριτά
Τώρα που ξέρουμε το νικητή, εκπαιδεύουμε σε όλα τα δεδομένα
final_model <- auto.arima(JohnsonJohnson, lambda = 0)
final_forecast <- forecast(final_model, h = 12) # 3 χρόνια μπροστά
autoplot(final_forecast) +
ggtitle("Πρόβλεψη 1981-1983") +
ylab("Έσοδα") +
theme_minimal()
TODO 14 (BONUS): Επιχειρηματική ανάγνωση
Υπολόγισε το μέσο ετήσιο growth rate των προβλεπόμενων ετών Συγκρίνετέ το με το ιστορικό growth rate — είναι ρεαλιστικό;
#Φτιαξαμε ενα dataframe για να υπολογισουμε το μεσο ετησιο growth rate.
yearly_growth_df <- data.frame(
year = numeric(),
full_year_total = numeric(),
growth_rate = numeric()
)
#This function calculates the yearly
calculate_growth_rate <- function(start_year,df){
prev_year_total <- NULL
for (i in 1:21){
current_year <- start_year + i
year <- window(JohnsonJohnson,start=c(current_year,1),end=c(current_year,4))
current_year_total <- sum(year)
if (!is.null(prev_year_total)){
growth <- ((current_year_total - prev_year_total) / prev_year_total) * 100
growth <- round(growth, 2)
}
else {
growth <- NA # The first year (1960) has no 1959 data to compare to
}
new_row <- data.frame(
year=current_year,
full_year_total = current_year_total,
growth_rate = growth
)
# Append the new row to our temporary data frame
df <- rbind(df, new_row)
# Update the previous year total for the next loop iteration
prev_year_total <- current_year_total
}
return(df)
}
# 3. Run the function and assign the result back to your variable
yearly_growth_df <- calculate_growth_rate(1959, yearly_growth_df)
# 4. View dataset
#View(yearly_growth_df)
head(yearly_growth_df, 25)## year full_year_total growth_rate
## 1 1960 2.63 NA
## 2 1961 2.77 5.32
## 3 1962 3.01 8.66
## 4 1963 3.40 12.96
## 5 1964 4.16 22.35
## 6 1965 5.16 24.04
## 7 1966 6.06 17.44
## 8 1967 6.81 12.38
## 9 1968 8.19 20.26
## 10 1969 9.54 16.48
## 11 1970 13.50 41.51
## 12 1971 16.29 20.67
## 13 1972 19.35 18.78
## 14 1973 23.31 20.47
## 15 1974 25.20 8.11
## 16 1975 28.62 13.57
## 17 1976 31.77 11.01
## 18 1977 38.07 19.83
## 19 1978 45.00 18.20
## 20 1979 51.84 15.20
## 21 1980 58.50 12.85#Now let preview the predicted quarter values.
# 1. Convert the ts object into a matrix so kable handles the years as row names
forecast_matrix <- as.matrix(t(matrix(final_forecast$mean, nrow=4,
dimnames=list(c("Q1", "Q2", "Q3", "Q4"), 1981:1983))))
# 2. Render the beautiful table
kable(t(forecast_matrix),
digits = 2,
caption = "Προβλέψεις Εσόδων Johnson & Johnson (1981-1983)")| 1981 | 1982 | 1983 | |
|---|---|---|---|
| Q1 | 18.52 | 21.60 | 25.13 |
| Q2 | 17.06 | 19.83 | 23.10 |
| Q3 | 18.89 | 21.89 | 25.53 |
| Q4 | 13.47 | 15.69 | 18.27 |
τωρα θα συγκρινουμε το growth rate των ετων που προβλεψαμε για να δουμε αν ειναι ρεαλιστικό. * 1981 -> 11.23% * 1982 -> 16.64% * 1983 -> 16.64%
- 🏆 Ποιο μοντέλο νίκησε στο test set; Με πόση διαφορά από το baseline (Seasonal Naïve); ARIMA(p,d,q) Κέρδισε το -> ARIMA(2,0,0)(1,1,0)[4] with drift
Ο seasonal naive έχει:
Mean Error: 2.5425000 RMSE: 2.7765401 MAPE: 17.89911
το ARIMA model έχει: Mean Error: -0.142883118 RMSE: 0.8136629 MAPE: 5.419212
- 🔍 Τι μοντέλο διάλεξε η
auto.arima(); Ερμηνεύστε τα (p,d,q)(P,D,Q)₌.
Κέρδισε το -> ARIMA(2,0,0)(1,1,0)[4]
Το (2,0,0) ειναι η βραθυπρόθεσμη μνήμη χωρίς εποχικότητα.
΄\(p = 2\) (Autoregressive - AR): δηλαδή 2 παρελθοντικά lags. d παράμετρος (Differencing) -> ο αριθμός διαφόρησεων που κάναμε ΄\(q = 0\) (Moving Average - MA)΄ : Δεν υπάρχουν παρελθοντικά σφάλματα (past shocks). Το σημερινό τρίμηνο δεν επηρεάζεται από το τυχαίο “σφάλμα” ή τις εκπλήξεις του προηγούμενου τριμήνου.
Το (1,1,0) μας εξηγεί την μακρυπρόθεσμη επίδραση με βάση τα προηγούμενα χρόνια
p = 1 (δηλαδή η περσινή ίδια περίοδος είναι στατιστικά σημαντική) d = 1 (κάναμε εποχική διαφόριση) q = 0
το drift υποδηλώνει οτι υπάρχει μια ανοδική τάση μακροπρόθεσμα
- 📈 Πόσο μεγαλώνει το 95% prediction interval για το 1ο τρίμηνο vs το 12ο; Τι σας λέει αυτό για την αβεβαιότητα;
| Τρίμηνο | Lower | Upper | Lower | Upper |
|---|---|---|---|---|
| 1981 Q1 | 16.53479 | 20.73297 | 15.57364 | 22.01253 |
| 1981 Q2 | 15.17457 | 19.18062 | 14.26218 | 20.40766 |
| 1981 Q3 | 16.68561 | 21.37892 | 15.62610 | 22.82850 |
| 1981 Q4 | 11.87766 | 15.26528 | 11.11443 | 16.31355 |
| 1982 Q1 | 18.40188 | 25.34630 | 16.90664 | 27.58795 |
| 1982 Q2 | 16.84548 | 23.35265 | 15.45032 | 25.46140 |
| 1982 Q3 | 18.51659 | 25.88334 | 16.94573 | 28.28271 |
| 1982 Q4 | 13.25288 | 18.56606 | 12.12155 | 20.29887 |
| 1983 Q1 | 20.58554 | 30.67519 | 18.52304 | 34.09080 |
| 1983 Q2 | 18.87243 | 28.28574 | 16.95557 | 31.48350 |
| 1983 Q3 | 20.77142 | 31.37172 | 18.62382 | 34.98935 |
| 1983 Q4 | 14.84986 | 22.47328 | 13.30742 | 25.07811 |
- 💼 Σύσταση προς τον portfolio manager: βάση των προβλέψεών σας, θα συστήνατε long ή short θέση στη μετοχή; Τεκμηριώστε.
Θα προτείνουμε long θέση για την μετοχή καθώς υπάρχει ανοδική αύξηση και τα επόμενα χρόνια προβλέπεται αύξηση 11-16% ανα έτος.