library(forecast)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
library(gridExtra)

Εισαγωγή και παρουσίαση Dataset

Στην παρούσα οικονομική ανάλυση, αναλαμβάνουμε τον ρόλο ενός χρηματοοικονομικού αναλυτή με σκοπό να μελετήσυομε την δυναμική τριμηναίψν κερδών ανά μετοχή της εταιρίας Jonhson & Jonhson για την περίοδο 1960-1980.

Στόχος μα είναι η δημιουργία μιας τεκμηριωμένης και αξιόπιστης πρόβληψης για τα επόμενα 3 έτη , προκειμένου να ενημερώσουμε τον Manger του πορτοφολίου της τράπεζας για την απόφαση της θέσης που θα λάβει το fund .

Exploration & Decomposition

data("JohnsonJohnson")
jj_data <- JohnsonJohnson

class(jj_data)
## [1] "ts"
start(jj_data)
## [1] 1960    1
end(jj_data)
## [1] 1980    4
frequency(jj_data)
## [1] 4
length(jj_data)
## [1] 84
# 1. Οπτικοποίηση Χρονοσειράς 
autoplot(jj_data) +
  labs(title = "Τριμηνιαία κέρδη ανά μετοχή  - Johnson & Johnson",
       x = "Έτος",
       y = "Κέρδη ανά μετοχή ") +
  theme_minimal()

# 2. Εποχικός Έλεγχος 
ggseasonplot(jj_data, year.labels = TRUE, year.labels.left = TRUE) +
  labs(title = "Εποχικό διάγραμμα κερδών ανά Έτος",
       x = "Τρίμηνο",
       y = "Κέρδη ανά μετοχή ") +
  theme_minimal()

ggsubseriesplot(JohnsonJohnson) +
  ggtitle("Subseries Plot - Μεταβολή ανά Τρίμηνο (Johnson & Johnson)") +
  xlab("Τρίμηνο") +
  ylab("Κέρδη ανά μετοχή ") +
  theme_minimal()

# 4. Αποσύνθεση Χρονοσειράς

jj_decomp <- decompose(jj_data, type = "multiplicative")

autoplot(jj_decomp) +
  labs(title = "Πολλαπλασιαστική Αποσύνθεση (Multiplicative Decomposition) της Χρονοσειράς") +
  theme_minimal()

# 5.  Αποσύνθεση με λογαριθμικό μετασχηματισμό 
jj_log <- log(JohnsonJohnson)
stl_fit <- stl(jj_log, s.window = "periodic")
autoplot(stl_fit) + 
  ggtitle("STL Αποσύνθεση") +
  theme_minimal()

Απαντήσεις ερωτημάτων Α

1. Γραμμική ή Εκθετική τάση ;

Η τάση της χρονοσειράς είναι καθαρά εκθετική. Όπως παρατηρείται στο αρχικό διάγραμμα, τα κέρδη ανά μετοχή αυξάνονται με επιταχυνόμενο ρυθμό όσο περνούν τα χρόνια (ειδικά κατά τη δεκαετία του 1970), αντί να παρουσιάζουν μια σταθερή, γραμμική ετήσια αύξηση.

2. Ποιο τρίμηνο έχει την ισχυρότερη εποχιακή ώθηση ;

Βάση των γραφιμάτων, οι κερδοφορίες εμφανίζουν συστηματικά την ισχυρότερη εποχιακή ώθηση το 3ο τρίμηνο κάθε έτους. Μια πιθανή εξήγηση αποτελεί η αύξηση τους τους καλοκαιρινούς μήνες , όπου παρατηρείται παραδοσιακά υψηλή ζήτηση για καταναλωτικά προϊόντα υγείας.

3. Multiplicative ή additive ;

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

Modeling & Forecasting

# 6. Train / Test split
train_set <- window(jj_data, end = c(1978, 4))
test_set <- window(jj_data, start = c(1979, 1))

# 7. Έλεγχος στασιμότητας 

print("--- Αρχικός Έλεγχος ADF ---")
## [1] "--- Αρχικός Έλεγχος ADF ---"
initial_adf <- adf.test(train_set)
## Warning in adf.test(train_set): p-value greater than printed p-value
print(initial_adf)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_set
## Dickey-Fuller = 0.85296, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary
# Απαιτήται διαμόρφωση για να γίνει στάσιμη 

train_seasonal_diff <- diff(train_set, lag = 4)
train_stationary <- diff(train_seasonal_diff, lag = 1)

print("--- Τελικός Έλεγχος ADF μετά τις Διαφορίσεις ---")
## [1] "--- Τελικός Έλεγχος ADF μετά τις Διαφορίσεις ---"
final_adf <- adf.test(train_stationary)
## Warning in adf.test(train_stationary): p-value smaller than printed p-value
print(final_adf) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_stationary
## Dickey-Fuller = -5.6292, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
# 8. ACF & PACF

 ggAcf(train_stationary) +
  ggtitle("ACF") + 
  theme_minimal()

 ggPacf(train_stationary)+
  ggtitle("PACF")+
  theme_minimal()

 # 9. Εκπαίδευση μοντέλων 
 
 fit_snaive <- snaive(train_set, h= 8)
 
 fit_hw <- hw(train_set, seasonal = "multiplicative", h = 8)
 
 fit_arima <- auto.arima(train_set, lambda = 0, stepwise = FALSE, approximation = FALSE)
fc_arima  <- forecast(fit_arima, h = 8)

summary(fit_snaive)
## 
## Forecast method: Seasonal naive method
## 
## Model Information:
## Call: snaive(y = train_set, 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_set, 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_set 
## 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
# 10. Residual Diagnostics 

checkresiduals(fit_snaive)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 123.73, df = 8, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 8
checkresiduals(fit_hw)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 71.136, df = 8, p-value = 2.919e-12
## 
## Model df: 0.   Total lags used: 8
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
# 11. Σύγκριση προβλέψεων 

autoplot(train_set) +
  autolayer(test_set,        series = "Πραγματικά") +
  autolayer(fit_snaive,      series = "Seasonal Naïve", PI = FALSE) +
  autolayer(fit_hw,          series = "Holt-Winters",   PI = FALSE) +
  autolayer(fc_arima,        series = "ARIMA",          PI = FALSE) +
  labs(title = "Σύγκριση Μοντέλων Πρόβλεψης — Johnson & Johnson",
       x = "Έτος",
       y = "Κέρδη ανά μετοχή",
       colour = "Μοντέλο") +
  theme_minimal()

# 12. Accuracy metrics 

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

rbind(
  "Seasonal Naive" = acc_snaive[2,],
  "Holt_Winters"  = acc_hw[2,],
  "Arima" = acc_arima[2,]
)[, c("RMSE", "MAE", "MAPE")]
##                     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
# 13. Τελική πρόβλεψη 

final_arima <- auto.arima(jj_data, lambda = 0, stepwise = FALSE, approximation = FALSE)

fc_final <- forecast(final_arima, h = 12)

print(fc_final)
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1981 Q1       18.50781 16.51885 20.73625 15.55400 22.02257
## 1981 Q2       17.21349 15.31420 19.34832 14.39514 20.58362
## 1981 Q3       18.68227 16.50472 21.14713 15.45669 22.58098
## 1981 Q4       13.53940 11.96128 15.32573 11.20176 16.36487
## 1982 Q1       21.58352 18.24246 25.53649 16.68856 27.91423
## 1982 Q2       20.07410 16.92950 23.80280 15.46945 26.04936
## 1982 Q3       21.78697 18.28444 25.96045 16.66436 28.48428
## 1982 Q4       15.78943 13.25108 18.81403 12.07698 20.64310
## 1983 Q1       25.17036 20.42500 31.01821 18.28664 34.64535
## 1983 Q2       23.41010 18.96299 28.90012 16.96179 32.30984
## 1983 Q3       25.40762 20.49945 31.49096 18.29758 35.28048
## 1983 Q4       18.41339 14.85634 22.82210 13.26060 25.56844
autoplot(fc_final) +
  labs(title = "Πρόβλεψη Τριμηνιαίων Κερδών J&J — Επόμενα 3 Χρόνια ",
       x = "Έτος",
       y = "Κέρδη ανά μετοχή",
       caption = "Μοντέλο: ARIMA με λογαριθμικό μετασχηματισμό ς") +
  theme_minimal()

Απαντήσεις ερωτήσεων Β

1. Ποιο μοντέλο νίκηση το test set ;

Το ARIMA (2,0,0)(1,1,0)[4] with drift ανέδειξε σαφώς καλύτερη προβλεπτική ικανότητα στο test set, επιτυγχάνοντας MAPE 5.4% έναντι του 17.9% του Seasonal Naive .

2. Τι μοντέλο διάλεξε το auto.arima ;

Η auto.arima επέλεξε το ARIMA(2,0,0)(1,1,0)[4] with drift:

p=2, d=0, q=0: το μοντέλο χρησιμοποιεί τις 2 προηγούμενες τιμές (AR) χωρίς ανάγκη κανονικής διαφόρισης P=1, D=1, Q=0: 1 εποχικός AR όρος και 1 εποχική διαφόριση (lag=4) που αφαιρεί την εποχική μη-στασιμότητα with drift: σταθερά ανοδικής τάσης, συνεπής με την εκθετική ανάπτυξη της J&J

3. Συμπεριφέρονται ως white noise τα residuals ;

Μόνο τα residuals του ARIMA πέρασαν τον έλεγχο Ljung-Box(p-value = 0.578>> 0.05), επιβεβαιώνοντας white noise. Αυτό σημαίνει ότι το μοντέλο εξήντλησε όλη την πληροφορία της χρονοσειράς.

4. Πόσο μεγαλώνει το 95% prediction interval από το 1ο τρίμηνο vs το 12ο ; Τι λέει αυτό για την αβεβαιότητα;

Το 95% διάστημα εμπιστοσύνης διευρύνεται σημαντικά με τον χρόνο : 1981 Q1: εύρος ~6.5 (15.6 – 22.0) 1983 Q3: εύρος ~17.0 (18.3 – 35.3)

Αυτό αντικατροπτρίζει τη φυσική αύξηση της αβεβαιότητας όσο προβλέπουμε μακρύτερα. Στην πράξη, οι προβλέψεις του 1ου έτους είναι αξιόπιστες, ενώ του 3ου πρέπει να αντιμετωπίζονται με μεγαλύτερη επιφύλαξη.

5. Σύσταση προσ τον portfolio manager: long ή short;

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