Προετοιμασία Δεδομένων

class(jj)
## [1] "ts"
start(jj); end(jj); frequency(jj)
## [1] 1960    1
## [1] 1980    4
## [1] 4
length(jj)  
## [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

Μέρος Α - Exploration & Description

Α.1 Οπτικοποιήση Χρονοσειράς

autoplot(jj) +
  ggtitle("Κέρδη τριμήνου Ανά Μετοχή της Johnson & Johnson 1960-1980") +
  xlab("Έτος") +
  ylab("Τιμή σε δολάρια") +
  theme_minimal()

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

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

ggseasonplot(jj, year.labels = TRUE) +
  ggtitle("Εποχιακό Διάγραμμα - Johnson & Johnson")

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

ggsubseriesplot(jj) +
  ggtitle("Subseries plot")

Από το subseries plot φαίνεται πως το τρίμηνο με τις καλύτερες επιδώσεις είναι το 3ο, με το 2ο και 1ο να ακολουθούν. Το 4ο εξάμηνο παρουσιάζει σημαντικά χαμηλότερα κέρδη σε σχέση με τα άλλα τρία.

A.3 Τύπος εποχικότητας

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

A.4 Decomposition

dec <- decompose(jj, type = "multiplicative")
  autoplot(dec) +
    ggtitle("Multiplicative Decomposition - Johnson & Johnson")

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

A.5 Λογαριθμικός Μετασχηματισμός

jjlog <- log(jj)

series1 <- autoplot(jj) + 
  ggtitle("Κανονική Χρονοσειρά") +
  theme_minimal()

series2 <- autoplot(jjlog) + 
  ggtitle("Log Transformation") +
  theme_minimal()

grid.arrange(series1, series2, ncol=2, top = "Διαφορά Log Transformation & Κανονικής Χρονοσειράς")

Εδώ παρατηρείται πως η εφαρμογή log transformation στην χρονοσειρά κάνει την εποχικότητα να παραμένει περισσότερο σταθερή. Αυτό γιατί ο λογάριθμος μετατρέπει τον πολλαπλασιασμό σε άθροισμα. Συνεπώς, με την λογαρίθμηση της χρονοσειράς αυτής, μετατρέπεται από multiplicative σε additive.

Μέρος Β - Modeling & Forecasting

B.1 Train/Test Split

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

cat("Αριθμός εγγραφών στο test είναι",  length(test),"ενώ στο train είναι",length(train))
## Αριθμός εγγραφών στο test είναι 8 ενώ στο train είναι 76

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

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

Εδώ το p-value είναι τρομερά υψηλό με τιμή που προσεγγίζει την μονάδα. Παρόλα αυτά, από το διάγραμμα γνωρίζουμε πως η χρονοσειρά δεν είναι στάσιμη, συνεπώς το αποτελέσμα αυτό είναι αφερέγγιο. Ως εκ τούτου, αφαιρούμε την τάση και την εποχικότητα, μέσω διαφοροποιήσης και λογαρίθμησης.

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

autoplot(jj_diff)+
  ggtitle("Μετά από log + 1η διαφόριση + εποχιακή διαφόριση")

adf.test(jj_diff)
## Warning in adf.test(jj_diff): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  jj_diff
## Dickey-Fuller = -6.8701, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

B.3 ACF & PACF

p1 <- ggAcf(jj_diff)  + ggtitle("ACF")
p2 <- ggPacf(jj_diff) + ggtitle("PACF")
grid.arrange(p1, p2, ncol = 2)

Εδώ παρατηρούνται spikes στο 1ο και στο 7ο spike στο ACF και στο 1ο και το 4ο spike στο PACF. Μετά από το τελευταίο spike στο κάθε ένα φαίνεται σταδιακή μείωση και σύγκλιση προς την τιμή 0, με το ACF να έχει σταθερά πιο υψηλές τιμές. Δεδομένου πως το ACF κόβει ξαφνικά μετά το lag 7 και το PACF φθίνει σταδιακά τότε πρέπει να γίνει χρήση κανονικά MA(7).

B.4 Εκπαίδευση Μοντέλων

B.4.1 Seasonal Naive

fit_snaive <- snaive(train, h = 8)

autoplot(fit_snaive) +
  ggtitle("Seasonal Naïve Forecast") +
  autolayer(test, series = "Πραγματικά", color = 'red')

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

Στο Season Naive παρατηρείται πως το test βρίσκεται σχεδόν εξολοκλήρου εκτός του διαστήματος εμπιστοσύνης της πρόβλεψης. Αυτό, υποδικνυεί χαμηλή προγνωστική ικανότητα, και κακή ποιότητα προβλέψεων. To σφάλμα είναι υψηλό και η πρόβληψη απέχει πολύ από την θεμιτή απόδοση.

B.4.2 Holt-Winters

fit_hw <- hw(train, h = 8, seasonal = "multiplicative")

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

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

Εδώ το Holt-Winters φαίνεται να βρίσκεται όχι μόνο εντός του 95% διαστήματος εμπιστοσύνης, αλλά σε πολλά σημεία εφάπτεται με την γραμμή των πραγματικών τιμών. Συνολικά το σφάλμα είναι αρκετά χαμηλό και ειδικά σε σύγκριση με το Seasonal Naive.

B.4.3 ΑRIMA

fit_arima <- auto.arima(train, lambda = 0)

fc_arima <- forecast(fit_arima, h = 8)

autoplot(fc_arima) +
  ggtitle("ARIMA Forecast") +
  autolayer(test, series = "Πραγματικά", color = "red")

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 παρουσιάζει πολύ καλές επιδόσεις, με υψηλη ταύτηση της πραγματικης και προβλεπόμενης γραμμής. Οι τιμές του σφάλματος είναι οι χαμηλότερες όλων των μεθόδων κάτι το οποίο υποδεικνυεί ότι η μέθοδος αυτή λογικά είναι η καλύτερη, κάτι το οποίο χρήζει περισσότερης μελέτης.

B.5 Residual Diagnostics

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

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

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

Το Holt-Winter’s ακολουθεί την κανονική κατανομή και η χρονοσειρά των καταλοίπων παρουσιάζει μια ήπια εποχικότητα, κάτι το οποίο δεν θα έπρεπε, όπως και μετά από το 1971 και έπειτα η απόκλιση από το μηδέν μικραίνει. Όσον αφορά το ACF εχουμε σημαντικό αριθμό spikes εκτός των μπλε γραμμών. Συνολικά η HW ενώ δεν είναι βέλτιστη, βρίσκεται στην μέση μεταξύ των άλλων δύο μεθόδων.

B.5 Σύγκριση Προβλέψεων σε ένα Διάγραμμα

autoplot(train) +
  autolayer(test, series = "Πραγματικά") +
  autolayer(fit_snaive$mean, series = "Seasonal Naïve") +
  autolayer(fit_hw$mean, series = "Holt-Winters") +
  autolayer(fc_arima$mean, series = "ARIMA")

Εδώ φαίνεται πως η ARIMA είναι αυτή που εφάπτεται με την αρχική πραγματική χρονοσειρά, ενώ η Holt-Winters είναι για λίγο χειρότερη. Η Season Naive απέχει πολύ από την πραγματική τιμή.

Β.6 Accuracy Metrics

accuracy(fit_arima)
##                       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
accuracy(fit_hw)
##                     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
accuracy(fit_snaive)
##                     ME      RMSE       MAE      MPE     MAPE MASE      ACF1
## Training set 0.5884722 0.8214427 0.5920833 14.27463 14.77554    1 0.6642405

Η μέθοδος με τις καλύτερες τιμές είναι η ARIMA λόγω των χαμηλότερων τιμών σφάλματος. Ακολουθεί η HW και στο τέλος με μεγάλη διαφορά βρίσκεται η Seasonal Naive.

B.7 Εκπαίδευση σε όλα τα δεδομένα & Πρόβλεψη

fit_arima <- auto.arima(jj, lambda = 0)

fc_arima <- forecast(fit_arima, h = 12)

autoplot(fc_arima) +
  ggtitle("ARIMA Forecast")

summary(fit_arima)
## Series: jj 
## ARIMA(2,0,0)(1,1,0)[4] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1     ar2     sar1   drift
##       0.2686  0.2855  -0.2695  0.0382
## s.e.  0.1137  0.1214   0.1212  0.0042
## 
## sigma^2 = 0.007793:  log likelihood = 82.47
## AIC=-154.95   AICc=-154.14   BIC=-143.04
## 
## Training set error measures:
##                       ME      RMSE      MAE        MPE     MAPE     MASE
## Training set -0.00516599 0.4009399 0.273459 -0.1504525 6.505327 0.389751
##                    ACF1
## Training set -0.2525095

Ερωτήσεις

Ερωτήσεις Μέρους Α

  1. Τι παρατηρείτε για την τάση της J&J; Είναι γραμμική ή εκθετική; Η τάση της J&J είναι εκθετική καθώς ενώ ξεκινά μικρή πολύ γρήγορα μεγάλωνε και παρουσιάζεται καμπύλη στο κάτω μέρος της γραμμής.

  2. Σε ποιο τρίμηνο εμφανίζουν οι κερδοφορίες την ισχυρότερη εποχική ώθηση; Δώστε μια πιθανή επιχειρηματική εξήγηση. Στο τρίμηνο το οποίο εμφανίζεται η μεγαλύτερη εποχιακή ώθηση είναι το 3ο τρίμηνο (δηλαδή Ιούλιο, Αύγουστο και Σεπτέμβρη), και αυτό πιθανώς γιατί κατά το τρίτο τρίμηνο αυξάνονται οι παραγγελίες για παιδικά φάρμακα και είδη πρώτων βοθειών για σχολέια. Επίσης πολλά φαρμακεία παραγγέλνουν φάρμακα και εμβολια για να προλάβουν τις εποχικές αρρώστιες του φθινωπόρου.

  3. Multiplicative ή additive; Τεκμηριώστε την απάντηση με βάση τα plots. Είναι Multiplicative γιατί είναι εκθετική, και το έυρος της εποχικότητας δεν είναι σταθερό (πχ είναι πολλαπλάσιο το 19080 σε σχέση με το 1960).

Ερωτήσεις Μέρους Β

  1. Ποιο μοντέλο νίκησε στο test set; Με πόση διαφορά από το baseline (Seasonal Naïve) Το μοντέλο το οποίο νίκησε στο test set ηταν το ARIMA και μάλιστα με μεγάλη διαφορα σε σχέση με το Seasonal Naive.

  2. Τι μοντέλο διάλεξε η auto.arima(); Ερμηνεύστε τα (p,d,q)(P,D,Q). Το μοντέλο που διάλεξε το auto.arima() είναι το (2,0,0)(1,1,0)[4]. Αυτό σημαίνει για (p,d,q) = (2,0,0), οτι χρειάστηκαν 2 lags η δύο προηγούμενες παρατηρήσεις για μια νέα πρόβλεψη (p), 0 φορές έγινε αφαίρεση των ακατέργαστων δεδομένων από τις προηγούμενες τιμές του για την αφαίρεση τάσης (d) και 0 μέγεθος moving average παραθύρου (q), το οποίο υπολογίζει το ποσό white noise που θα χρειαστεί στο παρελθόν για μελλοντικές προβλέψεις. Για (P,D,Q)[m] = (1,1,0)[4], οπου το m=4 δείχνει τον αριθμό διαχωρισμού της εποχής (δηλαδή εδώ τα τρίμηνα). Αυτή η παρένθεση χρησιμοποιείται όταν παρουσιάζεται εποχικότητα. Συγκεκιρμένα το P = 1 σημαίνει πως το μοντέλο κοιτάει μόνο μια τιμή του ίδιου τριμήνου παλιότερου έτους για την νέα πρόβλεψη. Το D = 1 σημαίνει πως χρειάστηκε μια αφαίρεση εποχιακής τιμής, ενώ το Q = 0 σημαίνει πως για κάθε εποχή χρειάστηκε μηδενικό μέγεθος moving average παραθύρου.

  3. Συμπεριφέρονται ως white noise τα residuals; Τι σημαίνει αυτό; Τα Residuals του ARIMA συμπεριφέρονται ως white noise γιατί έχουν τυχαιότητα, δεν υπάρχει τάση ούτε και εποχικότητα.

  4. Πόσο μεγαλώνει το 95% prediction interval για το 1ο τρίμηνο vs το 12ο; Τι σας λέει αυτό για την αβεβαιότητα; Το 95% του 12ου τριμήνου είναι αρκετά μεγαλύτερο του 1ου στις προβλέψεις, καθώς όσο προχωράμε χρονικά στο μέλλον η αβεβαιότητα αυξάνεται. Είναι αναμενόμενο να είμαστε πιο βέβαιαοι για τις προβλέψεις ένα τρίμηνο μπροστά, επειδή λογικά οι συγκυρίες θα είναι ίδιες με το τωρινό, ενώ όσο πιο μελλοντικά πηγαίνουμε τόσο περισσότερο ανοίγει το prediction interval, καθώς δεν μπορούμε ρεαλιστικά να γνωρίζουμε τις ακριβείς συγκυρίες της αγοράς.

  5. Σύσταση προς τον portfolio manager: βάση των προβλέψεών σας, θα συστήνατε long ή short θέση στη μετοχή; Τεκμηριώστε. Με δεδομένη την πρόβλεψη αύξησης των κέρδων της μετοχής θα πρότεινα την long θέση, καθώς η long απαιτεί την άνοδο της μετοχής για την κερδοφορία, ενώ η sort προϋποθέτει την μειώση της, και άρα της πτώσης των κερδών. Η απόφαση αυτή κάνει την παραδοχή πως η πρόβλεψη θα παραμένει στο άνω κομμάτι του prediction interval. Γνωρίζοντας από το διάγραμμα πως υπάρχει μεγάλη μακροχρόνια αβεβαιότητα, υπάρχει περίπτωση να αξίζει και η short, αν η τιμή των κερδών γίνει όσο του πατώματος του prediction interval.