Untitled

Dataset

JohnsonJohnson (R base) — τριμηνιαία κέρδη ανά μετοχή (EPS) της Johnson & Johnson — 84 παρατηρήσεις, 1960–1980 (συχνότητα = 4)

Φόρτωση δεδομένων

set.seed(42)

data("JohnsonJohnson")
jj <- JohnsonJohnson

class(jj)            # έλεγχος: είναι ήδη ts object;

## [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

Παρατήρηση: Τις μέγιστες τιμές έχει το 3ο τετράμηνο

ΜΕΡΟΣ Α — Exploration & Decomposition

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

autoplot(jj) +
  ggtitle("Κέρδη ανά μετοχή της Johnson & Johnson (1949–1960)") +
  xlab("Τρίμηνο") +
  ylab("Κέρδη ανά μετοχή") +
  theme_minimal()

Με βάση το διαγραμμα φαίνεται πως:

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

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

# Seasonal plot: κάθε έτος ώς χωριστή γραμμή
ggseasonplot(jj, year.labels = TRUE) +
  ggtitle("Seasonal plot — Johnson & Johnson")

# Subseries plot: ένα mini-plot ανά τρίμηνο
ggsubseriesplot(jj) +
  ggtitle("Subseries plot — μεταβολή ανά τρίμηνο")

Είναι φανερό ότι υπάρχει εποχικότητα με την σείρα τους αυξανόμενου κέρδους ανα μετοχή: Q4 < Q1 < Q2 < Q3. Άρα το Q3 είναι ισχυρότερο από όλα τα τρίμηνα.

Η επιχειρηματική εξήγηση μπορεί να είναι οτι το καλοκαίρι ο κόσμος χρησιμοποιέι περισσότερα είδοι περιποίησης. Για παράδειγμα, στις δεκαετίες 60’ - 80’ το μαύρισμα ήταν στην μόδα, όποτε συχνά αγόραζαν Johnson’s Baby Oil. Ακόμα ένα παράδειγμα είναι Johnson’s Baby Powder που το χρησιμοποιούσαν για απορρόφηση ιδρώτα.

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

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

TODO 4: Decomposition

# Επειδή είδαμε ότι το πλάτος μεγαλώνει → multiplicative
dec_mult <- decompose(jj, type = "multiplicative")
autoplot(dec_mult) +
  ggtitle("Multiplicative Decomposition — Johnson & Johnson")

Είναι εμφανές ότι το trend component είναι καμπυλωτό με κοίλη προς τα πάνω κάτι που επιβεβαιώνει την αυξηση τιμής μετοχών με αύξοντα ρυθμό.

TODO 5: Log transformation

# STL (Seasonal-Trend-Loess) - αλλά δουλεύει μόνο additive,
# οπότε εφαρμόζουμε log() πρώτα για να γίνει additive
air_log <- log(jj)
stl_fit <- stl(air_log, s.window = "periodic")
grid.arrange(autoplot(dec_mult) +
  ggtitle("Multiplicative Decomposition"), autoplot(stl_fit) +
  ggtitle("STL Decomposition (σε log scale)"), ncol=2)

Μετά το λογάριθμο, αλλάζει:

Το μέγεθος των κύματων γίνονται ομοιόμορφα στο data component. Είναι προφανές αποτέλεσμα καθώς η εφαρμογή του log() μετατρέπει τον εποχιακό τύπο από multiplicative σε additive.
Το trend component τείνει να είναι ευθύγραμμη σε σχέση με την προηγούμενη καμπυλωτή εκδοχή.

ΜΕΡΟΣ Β — Modeling & Forecasting

TODO 6: Train/Test split — χρονολογικά!

# Train: 1960 Q1 — 1978 Q4  (76 παρατηρήσεις)
# Test:  1979 Q1 — 1980 Q4  ( 8 παρατηρήσεις = 2 χρόνια)

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

length(train)  # 76

## [1] 76

length(test)   # 8

## [1] 8

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

# Έλεγχος στασιμότητας
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

Παρατήρηση: Και από το plot στο TODO 1 και από το p-value φαίνεται πως η σειρά jj δεν είναι στάσιμη. Άρα πρέπει να εφαρμοστεί το log + diff + diff(lag=4)

air_diff <- diff(diff(log(jj)), lag = 4)

# Ξανά ελέγχουμε
adf.test(air_diff)

## Warning in adf.test(air_diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  air_diff
## Dickey-Fuller = -6.8701, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Τωρα το p-value είναι μικρό και Dickey-Fuller αρνητικό. Οπότε, η σειρά είναι στάσιμη.

TODO 8: ACF & PACF

# Δύο plots δίπλα-δίπλα
p1 <- ggAcf(air_diff)  + ggtitle("ACF")
p2 <- ggPacf(air_diff) + ggtitle("PACF")
grid.arrange(p1, p2, ncol = 2)

#### Σημαντικότερα lag

Στο ACF είναι:

lag 1
lag 7

Στο PACF είναι:

lag 1
lag 4

TODO 9: Εκπαίδευση ΤΡΙΩΝ μοντέλων στο train set, ορίζοντας το horizon = 8

Μοντέλο 1 — Seasonal Naïve (baseline)

# Πρόβλεψη: κάθε τρίμηνο = ίδιο τρίμηνο πέρυσι
fit_snaive <- snaive(train, h = 8)
autoplot(fit_snaive) +
  ggtitle("Seasonal Naive Forecast") +
  autolayer(test, series = "Πραγματικά", color = "red")

Μοντέλο 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.59970

autoplot(fit_hw) +
  ggtitle("Holt-Winters Forecast") +
  autolayer(test, series = "Πραγματικά")

Μοντέλο 3 — ARIMA (με auto.arima)

# Η 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.1186641

fc_arima <- forecast(fit_arima, h = 8)
autoplot(fc_arima) +
  ggtitle("ARIMA Forecast") +
  autolayer(test, series = "Πραγματικά")

# Seasonal Naïve (baseline)
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

# Holt-Winters
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

# ARIMA
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) να είναι το καλύτερο.

Η auto.arima() επέλεξε το ARIMA(2,0,0)(1,1,0)[4] with drift. Αυτό σημαίνει εξάρτηση από τα 2 προηγούμενα τρίμηνα (p=2) και το ίδιο τρίμηνο πέρυσι (P=1), με μία εποχική διαφόριση (D=1). Το drift καλύπτει την ανοδική τάση.

TODO 10: Residual diagnostics

# Ελέγχουμε αν τα residuals είναι white noise με μοντέλο 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

Συμπέρασμα: Τα residuals είναι white noise αφού p-value = 0.5788 > 0.05. Αυτό σημαίνει πως ένα μεγάλο ποσοστό residuals ήταν άχρηστη πληροφορία (white noise). Άρα,το μοντέλο που επιλέξαμε (ARIMA) φιλτράρει καλα την πληροφορία.

TODO 11: Σύγκριση προβλέψεων — ένα plot με όλα

# Όλες οι προβλέψεις σε ένα plot
autoplot(train, series = "Training") +
  autolayer(test,           series = "Πραγματικά") +
  autolayer(fit_snaive$mean, series = "Seasonal Naive") +
  autolayer(fit_hw$mean,     series = "Holt-Winters") +
  autolayer(fc_arima$mean,   series = "ARIMA") +
  ggtitle("Σύγκριση μοντέλων στο test set") +
  ylab("Κέρδη ανά μετοχή") +
  theme_minimal() +
  guides(colour = guide_legend(title = "Μοντέλο"))

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

TODO 12: Accuracy metrics

# Seasonal Naïve (baseline)
temp <- accuracy(fit_snaive, test)
temp

##                     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.7697759

acc_snaive <- temp ["Test set", ]

# Holt-Winters
temp <- accuracy(fit_hw, test)
temp

##                      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.290191

acc_hw <- temp ["Test set", ]

# ARIMA
temp <- accuracy(fc_arima, test)
temp

##                        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.1989989

acc_arima <- temp ["Test set", ]

Aν τα βάλουμε σε ένα πίνακα θα είναι:

table <- data.frame(
  Model = c("Seasonal Naive", "Holt-Winters", "ARIMA"),
  RMSE  = c(acc_snaive["RMSE"], acc_hw["RMSE"], acc_arima["RMSE"]),
  MAE   = c(acc_snaive["MAE"], acc_hw["MAE"], acc_arima["MAE"]),
  MAPE  = c(acc_snaive["MAPE"], acc_hw["MAPE"], acc_arima["MAPE"])
)

print(table)

##            Model      RMSE       MAE      MAPE
## 1 Seasonal Naive 2.7765401 2.5425000 17.899110
## 2   Holt-Winters 1.0865532 1.0427583  7.763851
## 3          ARIMA 0.8136629 0.7225012  5.419212

To μοντέλο που νικά είναι το ARIMA καθώς έχει τα μικρότερα RMSE, MAE και MAPE που σημαίνει ότι έχει τα ελάχιστα λάθη. Συγκριτικά με το Seasonal Naïve (baseline), το ARIMA υπερτερεί με τεράστια διαφορά, ρίχνοντας το σφάλμα MAPE στο 6.77%, καθώς λαμβάνει υπόψη την ανάπτυξη της εταιρείας αντί να αντιγράφει απλώς τις περσινές τιμές.

TODO 13: Τελική πρόβλεψη για τον portfolio manager

# Τώρα που ξέρουμε το νικητή, εκπαιδεύουμε σε όλα τα δεδομένα
final_model <- auto.arima(jj, lambda = 0)
final_forecast <- forecast(final_model, h = 12)  # 3 χρόνια μπροστά

autoplot(final_forecast) +
  ggtitle("Πρόβλεψη 1981-1983") +
  ylab("Κέρδη ανά μετοχή") +
  theme_minimal()

Είναι φανερό πως η πρόγνωση φαίνεται ομαλή και δείχνει αύξηση στα κερδη. Παρατηρούμε επίσης ότι η αβεβαιότητα (γκρι σκιά / 95% interval) είναι μικρή στο 1ο τρίμηνο αλλά αυξάνεται αισθητά στο 12ο. Η τελική σύσταση στον portfolio manager είναι ξεκάθαρα Long θέση (αγορά), καθώς η εκθετική ανοδική τάση κάνει οποιαδήποτε short θέση υπερβολικά ρισκοδόρα.