Ergasia011

============================================================

ΕΡΓΑΣΙΑ 011 — Time Series Forecasting

Dataset: Johnson & Johnson Quarterly Earnings (1960–1980)

Στόχος: Πρόβλεψη τριμηνιαίων κερδών για τα επόμενα 3 χρόνια

============================================================

# --- Φόρτωση δεδομένων ---
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

Οπτικοποιηση Χρονοσειράς

autoplot(jj) +
  ggtitle("Johnson & Johnson Quarterly Earnings") +
  xlab("Έτος") + ylab("Κέρδη") +
  theme_minimal()

Ως αναλυτής, βλέπω μια σαφή ανοδική τάση στα κέρδη της εταιρείας. Παράλληλα, παρατηρώ έντονη εποχικότητα που επαναλαμβάνεται κάθε έτος. Είναι σημαντικό να σημειωθεί ότι το πλάτος των εποχικών διακυμάνσεων αυξάνεται όσο αυξάνεται η τάση, γεγονός που υποδηλώνει ότι η εποχικότητα εξαρτάται από το μέγεθος των κερδών (multiplicative decomposition).

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

# Εποχικό plot (σύγκριση ετών)
ggseasonplot(jj, polar = FALSE) + ggtitle("Εποχικό Plot")

# Subseries plot (μέσος όρος ανά τρίμηνο)
ggsubseriesplot(jj) + ggtitle("Subseries Plot")

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

Είναι multiplicative γιατί οι διακυμάνσεις έχουν μεγαλύτερο πλάτος καθώς αυξάνεται η τάση. Σε ένα addictive μοντέλο οι διακυμάνσεις θα έπρεπε να έχουν σταθερό πλάτος.

Decomposition

decomp <- decompose(jj, type = "multiplicative")
autoplot(decomp) + theme_minimal()

Log Transformation

p1 <- autoplot(jj) + ggtitle("Αρχική Σειρά")
p2 <- autoplot(log(jj)) + ggtitle("Με Log transformation")

grid.arrange(p1, p2, ncol=2)

Μετά τον λογάριθμο, το πλάτος των εποχικών διακυμάνσεων σταθεροποιείται. Το log() μετατρέπει μια πολλαπλασιαστική σχέση σε μια προσθετική (additive) σχέση. Αυτό κάνει τα δεδομένα πιο εύκολα στην ανάλυση και στην πρόβλεψη.

TODO 6: Train/Test Split

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

Έλεγχος Στασιμοτητας

adf.test(train) 

## Warning in adf.test(train): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train
## Dickey-Fuller = 0.85296, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary

train_diff <- diff(diff(log(train), lag = 4))
adf.test(train_diff)

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

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

ACF & PACF

ggAcf(train_diff)

ggPacf(train_diff)

Εκπαίδευση μοντέλων

fit_snaive <- snaive(train, h = 8)

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

# ARIMA
fit_arima <- auto.arima(train, lambda = 0) 
fc_arima <- forecast(fit_arima, h = 8)

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

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

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

Residual diagnostics

# Έλεγχος για το καλύτερο μοντέλο (συνήθως το ARIMA ή το HW)
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

Σύγκριση προβλέψεων

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") +
  ggtitle("Σύγκριση Προβλέψεων με Πραγματικά Δεδομένα") +
  xlab("Έτος") + ylab("Κέρδη") +
  guides(colour=guide_legend(title="Μοντέλο")) +
  theme_minimal()

#Accurasy Metrics

acc_snaive <- accuracy(fit_snaive, test)[2,]
acc_hw <- accuracy(fit_hw, test)[2,]
acc_arima <- accuracy(fc_arima, test)[2,]


comparison_table <- rbind(acc_snaive, acc_hw, acc_arima)
rownames(comparison_table) <- c("SNaive", "Holt-Winters", "ARIMA")

# Εμφάνιση του πίνακα
print(comparison_table[, c("RMSE", "MAE", "MAPE")])

##                   RMSE       MAE      MAPE
## SNaive       2.7765401 2.5425000 17.899110
## Holt-Winters 1.0865532 1.0427583  7.763851
## ARIMA        0.8136629 0.7225012  5.419212

Το μοντέλο που νικά είναι αυτό με τη μικρότερη τιμή RMSE δηλαδή το ARIMA

Τελική πρόβλεψη

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

final_fc <- forecast(final_model, h = 12)

autoplot(final_fc) +
  ggtitle("Τελική Πρόβλεψη κερδών J&J για τα επόμενα 3 χρόνια") +
  xlab("Έτος") + ylab("Κέρδη") +
  theme_minimal()