Η ανεργία αποτελεί ένα χρόνιο πρόβλημα της χώρας μας, μιας και ιστορικά τα τελευταία 25 χρόνια είναι σε υψηλότερα επίπεδα του ευρωπαϊκού μέσου και τις χώρες του ΟΟΣΑ. Το φαινόμενο επιδεινώθηκε τα χρόνια της οικονομικής κρίσης όπου στις χειρότερες στιγμές της το ένα τέταρτο του ενεργού πληθυσμού δεν μπορούσε να βρει εργασία. Ακόμα χειρότερη ήταν η κατάσταση για τους νέους της χώρας μιας και η νεανική ανεργία ανήλθε στο 46% , η οποία είναι και χειρότερη επίδοση στην ΕΕ.
Πριν ξεκινήσουμε να αναλύουμε κρίνω σημαντικό να δούμε τι σημαίνει ο όρος ανεργία και πώς προκύπτει αυτό το ποσοστό. Όταν αναφερόμαστε σε άνεργους αναφερόμαστε σε άτομα τα οποία δεν εργάζονται, αλλά δεν περιλαμβάνονται σε αυτά σπουδαστές, συνταξιούχοι, δηλαδή άτομα τα οποία έχουν την ηλικία (18-67), αλλά και τη δυνατότητα να εργαστούν. Για να υπολογιστεί το ποσοστό ανεργίας:
\text{Ποσοστό ανεργίας} = \frac{\text{Αριθμός ανέργων}}{\text{Εργατικό Δυναμικό}}
Σε αυτό το άρθρο έχω ως σκοπό τη πρόβλεψη της πορείας της ανεργίας τους επόμενους μήνες. Έτσι λοιπόν πήρα κάποια ιστορικά δεδομένα για την ανεργία στην ΕΕ, τον ΟΟΣΑ και την χώρα μας. Θα χρησιμοποιήσω ένα απλό μοντέλο (S)ARIMA, προκειμένου να κάνω μία εκτίμηση του μεγέθους της τους επόμενους μήνες. Τα δεδομένα που χρησιμοποιώ κυμαίνονται από τη περίοδο του 1998 μέχρι και το 2022. Αν θέλετε μία γρήγορη απάντηση, στη συγκεκριμένη ανάλυση προβλέπω ότι η πτωτική τάση της ανεργίας αναμένεται να συνεχιστεί τους επόμενους μήνες. Τον Φεβρουάριο του 2023, αυτή θα κυμαίνεται μεταξύ του 10% - 13%.
Για την ανάλυση αυτή θα χρειαστούμε τυπικές βιβλιοθήκες για την εισαγωγή και την επεξεργασία των δεδομένων μου, όπως οι readr (Wickham, Hester, & Bryan, 2024) και dplyr (Wickham, François, Henry, Müller, & Vaughan, 2023). Το πακέτο kableExtra (Zhu, 2024) χρησιμοποιήθηκε για την αποτύπωση των αποτελεσμάτων σε μορφή πίνακα, ενώ το πακέτο flextable (Gohel & Skintzos, 2024) χρησιμοποιήθηκε για τους πίνακες των αποτελεσμάτων των τεστ Dickey-Fuller και KPSS.
Στη συνέχεια, λόγω της φύσης των δεδομένων (χρονοσειρές) κρίθηκε απαραίτητη η χρήση σχετικών βιβλιοθηκών όπως τα πακέτα lubridate(Spinu, Grolemund, & Wickham, 2024), tseries(Trapletti & Hornik, 2024) & forecast(R. Hyndman κ.ά., 2025).
Τέλος, χρησιμοποιήθηκε το highcharter, ώστε να δημιουργήσω βασικά τα διαγράμματα pacf, acf μέσω ραβδογραμμάτων και κάποια διαγράμματα πρόβλεψης τάσης της ανεργίας με απλά γραφήματα γραμμών.
# General purpose R libraries
library(dplyr)
library(readr)
library(kableExtra)
library(flextable)
library(reactable)
library(purrr)
library(tibble)
# Graphs
library(highcharter)
library(glue)
# Time Series
library(urca)
library(lubridate)
library(tseries)
library(forecast)
library(tsibble)
library(feasts)
library(fable)
library(strucchange)Αφού συμπεριλάβω τις απολύτως απαραίτητες βιβλιοθήκες, μπορώ να εισάγω τα δεδομένα μου, μέσω της βιλιοθήκης readr. Το αρχείο δεδομένων μου είναι ένα τυπικό αρχείο csv και για την χρήση αυτού θα χρησιμοποιήσω την εντολή read_csv(). Το αρχείο περιέχει δεδομένα ανεργίας
Προκειμένου να αποφύγω την επανάληψη του κώδικά μου θα ορίσω και μερικές συναρτήσεις. Αρχικά, θέλω να έχω πίνακες που βασίζονται στη βιβλιοθήκη reactable με ορισμένες αλλαγές στην εμφάνιση για να φαίνονται ωραία. Για αυτό σε πολλούς πίνακες θα δείτε χρήση της reactable_custom αντί απλώς της reactable εντολής.
reactable_custom = function(data, head_max = 5){
reactable(
head(data, head_max),
bordered = FALSE, # No borders (booktabs style)
striped = FALSE, # No stripes for a clean look
highlight = FALSE, # No hover highlighting
defaultColDef = colDef(align = "center"),
style = list(fontSize = "14px", border = "none"),
theme = reactableTheme(
borderColor = "transparent", # Remove border
cellStyle = list(
borderBottom = "transparent" # Subtle line between rows
),
headerStyle = list(
borderBottom = "2px solid rgb(117, 117, 117)", # Thick top rule
borderTop = "2px solid rgb(117, 117, 117)",
fontWeight = "bold"
)
)
)
}Επιπλέον, έχω γράψει και άλλες δύο συναρτήσεις για τη κατασκευή διαγραμμάτων αυτοσυσχέτισης (ACF) και μερικής αυτοσυσχέτισης (pACF). Τις έχω ονομάσει acf_plot_function και pacf_plot_function αντίστοιχα. Και οι δύο συναρτήσεις βασίζονται στη βιβλιοθήκη highcharter για τη δημιουργία των αντίστοιχων διαγραμμάτων.
acf_plot_function = function(data, min = -1, max = 1){
ts_data <- data
N <- length(ts_data)
# Compute confidence limit
conf_limit <- qnorm(0.975) / sqrt(N)
# Use forecast::Acf (biased estimator, recommended for ARIMA)
acf_obj <- Acf(ts_data, plot = FALSE, lag.max = 30)
# Extract lags and values
# Note: Acf() returns acf_obj$acf as a matrix; acf_obj$lag as a matrix
lags_acf <- as.numeric(acf_obj$lag[-1]) # skip lag 0
acf_values <- as.numeric(acf_obj$acf[-1]) # skip lag 0
# Build the highcharter ACF plot
acf_chart <- highchart() %>%
hc_chart(type = "column") %>%
hc_title(text = "Autocorrelation Plot") %>%
hc_xAxis(categories = lags_acf, title = list(text = "Lag")) %>%
hc_yAxis(min = min, max = max,
title = list(text = "ACF"),
plotLines = list(
list(value = 0, color = "#000000", width = 3),
list(value = conf_limit, color = "#FF0000", dashStyle = "ShortDash", width = 1),
list(value = -conf_limit, color = "#FF0000", dashStyle = "ShortDash", width = 1)
)
) %>%
hc_legend(enabled = FALSE) %>%
hc_add_series(name = "ACF", data = round(acf_values, digits = 4))
# Show chart
acf_chart
}
pacf_plot_function <- function(data, min = -1, max = 1) {
ts_data <- data
N <- length(ts_data)
# Compute confidence limit (same formula as for ACF)
conf_limit <- qnorm(0.975) / sqrt(N)
# Compute PACF using forecast::Pacf (biased estimator, recommended for ARIMA)
pacf_obj <- forecast::Pacf(ts_data, plot = FALSE, lag.max = 20)
# Extract lags and PACF values
lags_pacf <- as.numeric(pacf_obj$lag) # already excludes lag 0
pacf_values <- as.numeric(pacf_obj$acf)
# Build the highcharter PACF plot
pacf_chart <- highcharter::highchart() %>%
highcharter::hc_chart(type = "column") %>%
highcharter::hc_title(text = "Partial Autocorrelation Plot") %>%
highcharter::hc_xAxis(categories = lags_pacf, title = list(text = "Lag")) %>%
highcharter::hc_yAxis(
min = min, max = max,
title = list(text = "PACF"),
plotLines = list(
list(value = 0, color = "#000000", width = 3),
list(value = conf_limit, color = "#FF0000", dashStyle = "ShortDash", width = 1),
list(value = -conf_limit, color = "#FF0000", dashStyle = "ShortDash", width = 1)
)
) %>%
highcharter::hc_legend(enabled = FALSE) %>%
highcharter::hc_add_series(name = "PACF", data = round(pacf_values, digits = 4))
# Show chart
pacf_chart
}Το αρχείο περιέχει δεδομένα ανεργίας για διάφορες χώρες ή οντότητες χωρών όπως η Ευρωπαϊκή Ένωση και οι χώρες του ΟΟΣΑ (Οργανισμός Οικονομικής Συνεργασίας και Ανάπτυξης), γ
reactable_custom(unemployment)Τα δεδομένα μου αποτελούνται από 3 μεταβλητές (στήλες). Πιο συγκεκριμένα οι στήλες μου είναι οι εξής:
| Μεταβλητή | Τύπος Μεταβλητής | Περιγραφή |
|---|---|---|
LOCATION |
Ποιοτική (κατηγορική) |
Χώρα ή Οντότητα χωρών |
TIME |
Ποιοτική (διατάξιμη) |
Μήνας-Έτος που αναφέρεται η μέτρηση |
Value |
Ποσοτική (συνεχής) |
Ύψος ανεργίας (βάσει περιοχής και μήνα) |
Συνεπώς, το δείγμα μου αποτελείται από 3 μεταβλητές, εκ των οποίων οι δύο είναι ποιοτικές και μία ποσοτική που είναι και η τιμή που θέλω να προβλέψω (ανεργία). Σε αυτό το σημείο ίσως να πρέπει να τονιστεί ότι στη μεταβλητή LOCATION έχει τρεις τιμές, την Ελλάδα, τις χώρες του ΟΟΣΑ και τις χώρες των 27 Ευρωπαϊκών χωρών. Τέλος, όσον αφορά τη κατηγοριοποίηση της μεταβλητής του χρόνου, με δεδομένο ότι οι τιμές μου έχουν τη μορφή μήνα-έτος (μμ/ετος), δεν είναι ξεκάθαρο το είδος της. Θα μπορούσαμε να τη χωρίσουμε σε δύο επιπλέον μεταβλητές όπου η μία να είναι το έτος και να τη χαρακτηρίσουμε ως μία ποσοτική μεταβλητή και ο μήνας μία ποιοτική διατάξιμη μεταβλητή.
Η μεταβλητή που δηλώνει το μήνα για τον οποίο αναφέρεται η αντίστοιχη ανεργία (TIME) αναγνωρίζεται αυτόματα ως ένα διάνυσμα χαρακτήρων. Ένα από τα πρώτα πράγματα που πρέπει να κάνουμε όταν χειριζόμαστε δεδομένα που δηλώνουν διάστημα χρόνου είναι να τα μετατρέψουμε στο αντίστοιχο είδος μεταβλητής, που στην R, αυτό το είδος καλείται Date. Ο παρακάτω πίνακας δηλώνει τις υπάρχουσες μεταβλητές, καθώς και το είδος το οποίο τους έχει αποδοθεί αυτόματα με βάση τις τιμές που περιέχουν. Η R έκανε καλή δουλειά και εντόπισε ότι η μεταβλητή Value πρόκειται για ένα διάνυσμα αριθμών μιας και αντιπορσωπεύει το ύψος της ανεργίας. Η αντιστοίχιση των τριών ονοτήτων είναι και αυτή σωστή μιας και αναφερόμαστε σε ονόματα αυτών και άρα θα είναι ένα διάνυσμα χαρακτήρων.
Παραπάνω επισημάνθηκε ότι οι ημερομηνίες έχουν την μορφή “ΕΕΕΕ-ΜΜ” (Έτος-Μήνας) και από το λογισμικό αναγνωρίστηκαν αυτόματα ως χαρακτήρες. Με τη βοήθεια του πακέτου lubridate θα μετατρέψω τη μεταβλητή του χρόνου σε τύπου Date.
Και αφού κάναμε τη μετατροπή, αν ελέγξουμε άλλη μία φορά θα δούμε ότι η αλλαγή ήταν επιτυχημένη.
Και τώρα που μετέτρεψα τη κύρια μεταβλητή μου σε τύπου Date είμαι έτοιμος να συνεχίσω την ανάλυσή μου.
Ούφ! Έχουμε κάποια καλά νέα. Σε αυτό το σύνολο δεδομένων υπάρχουν συνολικά 0 ελλείπουσες τιμές. Σε περίπτωση που το σύνολό μου έλειπαν παρατηρήσεις θα έπρεπε να ερευνήσω σε πρώτη φάση ποια από τις μεταβλητές αυτές παρατηρήθηκαν. Σε δεύτερη φάση και αναλόγως του τύπου της μεταβλητής θα έπρεπε είτε να διώξω εντελώς εκείνες τις σειρές - παρατηρήσεις ή θα μπορούσα να προσπαθήσω με διάφορες μεθόδους να προβλέψω τις τιμές τους.
Τα δεδομένα της ανεργίας για την Ελλάδα αναφέρονται στη περίοδο του Απριλίου του 1998 μέχρι και τον Αύγουστο του 2022. Όσον αφορά τα δεδομένα της ΕΕ, ξεκινάνε από τον Ιανουάριο του 2000 μέχρι και τον Αύγουστο το 2022. Τα τελευταία 20 χρόνια έχουν γίνει σοκαριστικά μεγάλες αλλαγές όσον αφορά το ποσοστό ανεργίας στη χώρα μας με την πιο απότομη μεταβολή να σημειώνεται τις περιόδους τις οικονικής κρίσης (μετά το 2009). Ενδεικτικά μπορούμε να παρατηρήσουμε τη διαφορά στο ύψος της ανεργίας μεταξύ του Σεπτέμβρίου του 2010, που ήταν στο 10% και τρία χρόνια αργότερα, τον Σεπτέμβριο του 2013, έφτασε στο 28.1%. Για λόγους πληρότητας παρακάτω επισυνάπτονται πίνακες που δείχνουν τους 5 μήνες με τη μεγαλύτερη και τη μικρότερη ανεργία στην Ελλάδα και στην Ευρώπη τα τελευταία 20 χρόνια.
unemploymentType$TIMEchr <- format(unemploymentType$TIME, "%Y %b")unemploymentType %>%
as.tibble() %>%
dplyr::filter(LOCATION == "GRC") %>%
select(TIMEchr, Value) %>%
arrange(Value) %>%
tail(5) %>%
setNames(c("Μήνας", "Ποσοστό ανεργίας (%)")) %>%
reactable_custom()
unemploymentType %>%
as.tibble() %>%
dplyr::filter(LOCATION == "EU27_2020") %>%
select(TIMEchr, Value) %>%
arrange(Value) %>%
tail(5) %>%
setNames(c("Μήνας", "Ποσοστό ανεργίας (%)")) %>%
reactable_custom()Από την άλλη μεριά έχει ενδιαφέρον να παρατηρήσουμε και τις περιόδους με τη χαμηλότερη παρατηρούμενη ανεργία. Στην Ελλάδα αυτή η περίοδος ήταν λίγο πριν την οικονομική κρίση, το 2008
unemploymentType %>%
as.tibble() %>%
dplyr::filter(LOCATION == "GRC") %>%
select(TIMEchr, Value) %>%
arrange(Value) %>%
head(5) %>%
setNames(c("Μήνας", "Ποσοστό ανεργίας (%)")) %>%
reactable_custom()
unemploymentType %>%
as.tibble() %>%
dplyr::filter(LOCATION == "EU27_2020") %>%
select(TIMEchr, Value) %>%
arrange(Value) %>%
head(5) %>%
setNames(c("Μήνας", "Ποσοστό ανεργίας (%)")) %>%
reactable_custom()On the other hand, at the moment in Europe we are at the lowest levels of the last 20 years, with unemployment hovering around 6%.
Finally, it is evident that in the case of Greece the changes are more intense than in Europe.
All of the above can be summarized by the following graph:
`summarise()` has grouped output by 'LOCATION'. You can override using the
`.groups` argument.country_colors <- c("#000","#1f77b4","#2ca02c")
highchart() %>%
hc_chart(type = "line") %>% # clean white background
hc_title(text = "Ποσοστό ανεργίας", style = list(fontSize = "20px", fontWeight = "bold")) %>%
hc_subtitle(text = "Ανά χώρα / οντότητα (μέση ανεργία ανά έτος)", style = list(color = "#666")) %>%
hc_xAxis(
title = list(text = "Έτος"),
gridLineWidth = 0,
labels = list(style = list(fontSize = "12px"))
) %>%
hc_yAxis(
title = list(text = "Ανεργία (%)"),
gridLineDashStyle = "ShortDot", # subtle dashed grid
labels = list(format = "{value} %", style = list(fontSize = "12px"))
) %>%
hc_tooltip(
shared = TRUE,
valueDecimals = 2,
valueSuffix = " %",
borderColor = "#333333",
backgroundColor = "#f9f9f9",
style = list(fontSize = "13px")
) %>%
hc_add_series(
data = g,
type = "line",
hcaes(x = Year, y = mean_unemployment, group = LOCATION),
) %>%
hc_plotOptions(
series = list(
lineWidth = 3,
marker = list(enabled = FALSE, symbol = "circle"), # hide cluttered markers
dataLabels = list(
enabled = TRUE,
formatter = JS("
function() {
if (this.point.index === this.series.data.length - 1) {
return Highcharts.numberFormat(this.y, 1) + ' %';
}
return null;
}"
),
style = list(fontWeight = "bold", color = "#333")
)
)
) %>%
hc_legend(
enabled = TRUE,
layout = "horizontal",
align = "center",
verticalAlign = "bottom"
)Σύγκριση χρονοσειρών ανεργίας μεταξύ της Ελλάδας, των χωρών του ΟΟΣΑ και των 27 χωρών της Ευρώπης (λαμβάνοντας υπόψιν και την έξοδο του Ηνωμένου Βασιλείου από την ΕΕ)
Στην ανάλυση χρονοσειρών είναι σημαντικό να ξεχωρίσουμε τις πηγές της διασποράς μιας χρονοσειράς και να διαπιστώσουμε από που προέρχεται αυτή. Οι χρονοσειρές έχουν τρία βασικά στοιχεία, τη τάση (T), την εποχικότητα (S) και την τυχαιότητα (E). Η τάση αναφέρεται στην The trend refers to whether the series shows a specific direction (upward/downward). Seasonality refers to recurring patterns at regular intervals.
y_t = S_t + T_t + E_t Όπου:
Παρόμοια, το πολλαπλασιαστικό μοντέλο:
y_t = S_t \cdot T_t \cdot E_t
όπου τα στοιχεία που συνθέτουν τη χρονοσειρά πολλαπλασιάζονται, αντί να προσθέτονται.
library(forecast)
library(highcharter)
decomp <- unemploymentType %>%
model(stl = STL(Value)) %>%
components() %>%
as_tibble() %>%
dplyr::filter(LOCATION == "GRC") %>%
mutate(TIME = as.Date(TIME),
trend = round(trend, 2),
season_year = round(season_year, 4),
remainder = round(remainder, 4)
)
# Build highchart with multiple y-axes (stacked)
highchart(type = "stock") %>%
hc_title(text = "Κατάτμηση χρονοσειράς σε τάση, εποχικότητα και θόρυβο") %>%
hc_subtitle(text = "Τα ιστορικά δεδομένα της ανεργίας στην Ελλάδα δείχνουν ότι δεν υπάρχει κάποια ξεκάθαρη τάση. Επιπλέον, παρατηρούμε μία ανεμική εποχικότητα μέχρι το 2012, ενώ τα επόμενα χρόνια παρατηρείται εντονότερα το φαινόμενο αυτής.") %>%
hc_caption(text="stesiam, 2023") %>%
# Observed
hc_add_series(decomp, "line", hcaes(x = datetime_to_timestamp(TIME), y = Value), name = "Παρατήρηση", yAxis = 0) %>%
# Trend
hc_add_series(decomp, "line", hcaes(x = datetime_to_timestamp(TIME), y = trend), name = "Τάση", yAxis = 1) %>%
# Seasonal
hc_add_series(decomp, "line", hcaes(x = datetime_to_timestamp(TIME), y = season_year), name = "Εποχικότητα", yAxis = 2) %>%
# Random
hc_add_series(decomp, "line", hcaes(x = datetime_to_timestamp(TIME), y = remainder), name = "Θόρυβος", yAxis = 3) %>%
#
# # Define 4 stacked y-axes
hc_yAxis_multiples(
list(
title = list(text = "Δεδομένα"),
height = "25%", top = "0%",
offset=0, labels = list(enabled = FALSE), gridLineWidth = 0),
list(
title = list(text = "Τάση"),
height = "25%", top = "25%",
offset=0, labels = list(enabled = FALSE), gridLineWidth = 0),
list(
title = list(text = "Εποχικότητα"),
height = "25%", top = "50%",
offset=0, labels = list(enabled = FALSE), gridLineWidth = 0),
list(
title = list(text = "Υπόλοιπο"),
height = "25%", top = "75%",
offset=0, labels = list(enabled = FALSE), gridLineWidth = 0)
) %>%
hc_tooltip(shared = TRUE, valueDecimals = 2) %>%
hc_rangeSelector(enabled = FALSE) %>%
hc_scrollbar(enabled = FALSE) %>%
hc_navigator(enabled = FALSE)Μία σημαντική έννοια στις χρονοσειρές είναι η στασιμότητα. Μία χρονοσειρά καλείται στάσιμη (‘Applied Time Series Analysis’, χ.χ.) αν:
Από αυτό το διάγραμμα είναι εμφανέστατο ότι η σειρά μας δεν κινείται γύρω από κάποια συγκεκριμένη τιμή, παραβιάζοντας την πρώτη προϋπόθεση για να θεωρηθεί μία χρονοσειρά στατική. Αυτό μας υποδεικνύει την ανάγκη χρήσης διαφορών πρώτης τάξης για την ανεργία της Ελλάδας.
highchart() %>%
hc_chart(type = "line") %>%
hc_title(text = "Ποσοστό ανεργίας") %>%
hc_subtitle(text = "Η κορύφωση της ανεργίας στην Ελλάδα συνέβη το 2013 όπου εκείνη τη περίοδο άνω του ένα τέταρτου του εργατικού δυναμικού δεν μπορούσε να βρει εργασία.") %>%
hc_xAxis(title = list(text = "Έτος")) %>%
hc_yAxis(title = list(text = "Ποσοστό ανεργίας (%)")) %>%
hc_tooltip(shared = TRUE, valueDecimals = 2, valueSuffix = " %") %>%
hc_add_series(
data = g,
type = "line",
hcaes(x = Year, y = mean_unemployment, group = LOCATION)
) %>%
hc_plotOptions(
series = list(
dataLabels = list(
enabled = TRUE,
formatter = JS("function() {
if (this.point.index === this.series.data.length - 1) {
return this.y.toFixed(2) + ' %';
}
return null;
}")
)
)
)Από τη πρώτη διαφορά παρατηρώ μεγάλη βελτίωση μιας και δεν έχουμε τις τεράστιες αποκλίσεις του προηγούμενου διαγράμματος. Οι τιμές ως επί το πλείστον δεν παρουσιάζουν κάποια τάση και κινούνται σε τιμές σχετικά κοντά στο μηδέν. Αυτό είναι ένα καλό στοιχείο, όμως έχω έναν ελαφρύ προβληματισμό καθώς υπάρχουν δύο σημεία στη χρονοσειρά με σχετικά μεγάλες αποκλίσεις από το μηδέν. Η πρώτη είναι στο σημείο ~150ο όπου η σειρά έχει και μία ελαφριά ανοδική τάση και η δεύτερη είναι στο 250ο σημείο.
grc_unemployment = unemployment %>% dplyr::filter(LOCATION == "GRC")
grc_unemployment_diff1 <- diff(grc_unemployment$Value, differences = 1)
highchart() %>%
hc_chart(type = "line") %>%
hc_title(text = "Πρώτες διαφορές ανεργίας στην Ελλάδα") %>%
hc_xAxis(title = list(text = "Έτος")) %>%
hc_yAxis(title = list(text = "Παρατήρηση")) %>%
hc_tooltip(shared = TRUE, valueDecimals = 2, valueSuffix = " %") %>%
hc_add_series(
data = grc_unemployment_diff1
)Λαμβάνοντας υπόψιν τους παραπάνω προβληματισμούς μου, έλαβα και τις δεύτερες διαφορές και τις οπτικοποίησα. Το διάγραμμα είναι σχεδόν το ίδιο, με τις τιμές να κυμαίνονται στο μηδέν ακόμα και στο προβληματικό σημείο κοντά στο 150ο σημείο.
grc_unemployment_diff2<- diff(grc_unemployment$Value, differences = 2)
highchart() %>%
hc_chart(type = "line") %>%
hc_title(text = "Δεύτερες διαφορές ανεργίας") %>%
hc_xAxis(title = list(text = "Έτος")) %>%
hc_yAxis(title = list(text = "Παρατήρηση")) %>%
hc_tooltip(shared = TRUE, valueDecimals = 2, valueSuffix = " %") %>%
hc_add_series(
data = grc_unemployment_diff2
)Ο γραφικός έλεγχος της στασιμότητας είναι ένας αρκετά εύκολος τρόπος για να διαπιστώσουμε την ύπαρξη τάσεων ή αν η σειρά μας έχει γενικότερα σταθερή συμπεριφορά. Αν εξαιρέσουμε κάποιες αρκετά ξεκάθαρες περιπτώσεις, θα υπάρχουν φορές που πολλοί μπορεί να διαφωνήσουν ως προς τη στασιμότητα της σειράς απλώς από τη πορεία της. Αυτό είναι λογικό μιας και ως μέτρο είναι κατά κάποιο τρόπο υποκειμενικό μιας και βασίζεται στην άποψη / ερμηνεία του καθενός που θα δώσει στη κίνηση της χρονοσειράς. Έτσι λοιπόν η χρήση αυτού του τρόπου μπορεί να οδηγήσει σε μη συνεπείς ή σταθερές αποφάσεις μιας και ένα άτομο μπορεί να θεωρήσει μία ελαφριά αυξητική τάση που επανέρχεται ως στάσιμη σε και κάποιο άλλο άτομο να μεταφράσει το μοτίβο ως ισχνή αλλά υπαρκτή τάση. Σε αναλογία των ελέγχων κανονικότητας που έχουμε γραφικούς (quantile-quantile γράφημα) αλλά και στατιστικούς ελέγχους (έλεγχος Kolmogorov-Smirnov, έλεγχος Shapiro-Wilk), έτσι έχουμε ανάλογες εναλλακτικές και για τον έλεγχο στασιμότητας. Με αυτό το τρόπο μπορούμε να έχουμε ένα πιο αντικειμενικό κριτήριο για το αν οι σειρές μας είναι στάσιμες ή όχι. Κάποιοι από τους πιο γνωστούς ελέγχους στασιμότητας, οι οποίοι είναι γνωστοί και ως έλεγχοι μοναδιαίας ρίζας, είναι οι εξής:
Και κάπου εδώ ξεκινάει το χάος που άρχισα να συνειδητοποιώ όταν ξεκίνησα να μαθαίνω την R. Χρησιμοποιώντας μία γλώσσα προγραμματισμού δεν έχουμε ωραία μενού και κουτάκια με επιλογές για τον κάθε έλεγχο. Στην R αλλά και σε άλλες γλώσσες, οπψς στη Python, υπάρχουν πακέτα που προσθέτουν λειτουργίες και δυνατότητες κάθε γλώσσας. Για τους ελέγχους στασιμότητας έχουν κατασκευαστεί αρκετά πακέτα της R, που επιτελούν παρόμοιο σκοπό. Κάποια αρκετά διάσημα πακέτα που προσφέρουν ελέγχους σταστιμότητας είναι τα:
Οκ. Ποιο από όλα όμως να διαλέξω; Έχουν διαφορά; Θα πρέπει να μελετήσουμε προσεκτικά τις επιλογές μας. Στα έτοιμα προγράμματα το λογισμικό λαμβάνει κάποιες λογικές αποφάσεις για εμάς ή μας δίνει τη δυνατότητα να διαλέξουμε κάποιες τιμές ή ελέγχους. Όταν κάνουμε αυτή τη μετάβαση σε μία γλώσσα προγραμματισμού έχουμε την ελευθερία να δούμε πώς προέκυψε αυτό το αποτέλεσμα (μελετώντας τον αντίστοιχο κώδικα του πακέτου). Αυτή η αυξημένη ελευθερία κινήσεων συνοδεύεται με αυξημένες ευθύνες από τη μεριά του αναλυτή αφού θα πρέπει να ελέγξει τι υπολογίζει το κάθε πακέτο, ποιες είναι οι δυνατότητές του και οι περιορισμοί του. Σε αυτή τη περίπτωση και τα τρία πακέτα και ειδικά τα δύο πρώτα, το tseries και το urca προωθούνται ως πακέτα που παρέχουν ελέγχους στασιμότητας, δηλαδή παρόμοια λειτουργικότητα . Εν συντομία, το πακέτο tseries με το οποίο έχω τη μεγαλύτερη εξοικείωση είναι αυτό που είναι αρκετά περιοριστικό. Παρά το γεγονός ότι πολλοί οδηγοί χρησιμοποιούν αυτό το πακέτο, συνειδητοποίησα ότι δεν μπορούσα να θέσω τα lags για τους ελέγχους ή να θέσω χαρακτηριστικά της χρονοσειράς, όπως αν έχει (αυξητική ή μειούμενη) τάση. Από την άλλη μεριά έχουμε το πακέτο urca, το οποίο απαντάει σε αυτούς τους περιορισμούς επιτρέποντας στους χρήστες να θέτουν τον αριθμό των lags καθώς και στοιχεία της χρονοσειράς. Το μοναδικό μειονέκτημα του urca πακέτου είναι η μη παροχή ενός p-value στα αποτελέσματα των ελέγχων του, το οποίο είναι και το πιο γνωστό μέτρο ερμηνείας των ελέγχων. Για τον έλεγχο της υπόθεσης υπολογίζεται η στατιστική τιμή και ελέγχεται με την αντίστοιχη κριτική τιμή και αν η στατιστική τιμή είναι μεγαλύτερη της τελευταίας τότε ο έλεγχος απορρίπτεται στο αντίστοιχο επίπεδο σημαντικότητας.
summary_stationarity_results <- data.frame(
"Τύπος" = c("Αρχ. Παρατηρήσεις", "Δ(GRC)", "Δ(GRC)"),
"ADF test" = c("Μη στάσιμη", "Στάσιμη", "Στάσιμη"),
"PP test" = c("Μη στάσιμη", "Στάσιμη", "Στάσιμη"),
"KPSS test" =c("Μη στάσιμη", "Μη στάσιμη", "Στάσιμη"),
check.names = FALSE
)
summary_stationarity_results %>% reactable_custom()Ο έλεγχος Dickey - Fuller είναι ένας από τους πιο απλούς ελέγχους μοναδιαίας ρίζας για να διαπιστώσουμε τη στασιμότητα ή μη μίας χρονοσειράς. Αυτός ο έλεγχος βασίζεται στο αυτοπαλίνδρομο μοντέλο πρώτης τάξης, \text{AR}(1):
y_t = \phi y_{t-1} + e_t
Όπου y_t είναι η τιμή της χρονοσειράς και το e_t ο όρος του σφάλματος. Δηλαδή είναι μία χρονοσειρά που η τιμές της επηρεάζονται - εξαρτώνται από τις προηγούμενες τιμές της.
\begin{equation} \begin{split} y_t & = \rho y_{t-1} + e_t \\ \\ y_t - y_{t-1}& = \rho y_{t-1} - y_{t-1} + e_t \\ \\ \Delta y_t & = (\rho - 1) y_{t-1} + e_t \\ \\ \Delta y_t & = \gamma y_{t-1} + e_t \end{split} \end{equation}
Ο έλεγχος έχει συγκεκριμένες παραλλαγές οι οποίες βασίζονται στη συμπεριφορά της εκάστοτε χρονοσειράς. Πιο συγκεκριμένα υπάρχουν τρεις παραλλαγές:
Ένας απο τους πιο χρησιμοποιούμενους ελέγχους μοναδιαίας ρίζας - στασιμότητας είναι ο επαυξημένος έλεγχος ADF (Augmented Dickey - Fuller) και αποτελεί γενίκευση του απλού ελέγχου Dickey - Fuller, μιας και εξαρτάτει από ένα αυτοπαλίνδρομο μοντέλο μεγαλύτερης τάξης (\rho >1). Δηλαδή εννοούμε ότι η τιμή της χρονοσειράς δεν εξαρτάται μόνο από τη προηγούμενη τιμή της αλλά και από άλλες \rho προηγούμενες. Ένα αυτοπαλίνδρομο μοντέλο \rho τάξης έχει την εξής μορφή:
y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + \epsilon_t Αντίστοιχα με τον απλό έλεγχο Dickey-Fuller υπάρχουν ορίζονται παραλλαγές οι οποίες βασίζονται στη συμπεριφορά της εκάστοτε χρονοσειράς:
Ο έλεγχος ADF έχει την εξής μορφή υποθέσεων: H_0 : \gamma = 0 \textsf{, η σειρά δεν είναι στασιμη}\\ H_1 : \gamma \neq 0 \text{, η σειρά είναι στασιμη} Δηλαδή, αν τα αποτελέσματά μου απορρίψουν τη μηδενική υπόθεση τότε αυτό μας υποδεικνύει ότι η χρονοσειρά είναι στάσιμη. Ας εφαρμόσουμε τον παραπάνω έλεγχο για την ανεργία της Ελλάδας και της Ευρώπης. Όμως πριν το κάνω αυτό θα πρέπει να διαπιστώσω πόσες χρονικές υστερήσεις (lags) πρέπει να λάβω υπόψιν μου. Αρχικά, θα πρέπει να θέσω το μέγιστο αριθμό τον οποίο μπορούν να λάβουν οι υστερήσεις μου, p_{max}. Αυτό μπορεί να προκύψει εμπειρικά από τη μορφή των δεδομένων μου, αν είναι ανά τετράμηνο θα θέσω p_{max} = 4 ή αν είναι μηνιαία ενδείκνυται να θέσω p_{max} = 12. Ένας πιο διαδεδομένος τρόπος είναι ο κανόνας του Schwert, ο οποίος βασίζεται στο πλήθος των παρατηρήσεων:
\begin{equation} \begin{split} p_{max} & = \left[12 \cdot \left( \frac{T}{100}\right)^{0.25} \right] \\ \\ & = \left[ 12 \cdot \left( \frac{293}{100}\right)^{0.25} \right]\\ \\ & = \left[ 15.7 \right] \\ \\ & = 15 \text{ υστερήσεις} \end{split} \end{equation}
Στη βιβλιογραφία υπάρχει και μία πιο συντηρητική μέθοδος καθορισμού μέγιστων υστερήσεων με παράγοντα 4 αντί για 12,
\begin{equation} \begin{split} p_{max} & = \left[4 \cdot \left( \frac{T}{100}\right)^{0.25} \right] \\ \\ & = \left[ 4 \cdot \left( \frac{293}{100}\right)^{0.25} \right]\\ \\ & = \left[ 5.23 \right] \\ \\ & = 5 \text{ υστερήσεις} \end{split} \end{equation}
Προς αποφυγή οποιασδήποτε σύμπτωσης, αυτή τη στιγμή δεν βρήκα τις υστερήσεις που ενδείκνυται να λάβω για τον επαυξημένο έλεγχο Dickey-Fuller. Έθεσα το άνω όριο αυτών που με βάση τον κανόνα του Schwerz είναι p_{max}. Δηλαδή ο αριθμός των υστερήσεων \rho είναι ένας ακέραιας αριθμός ανάμεσα μεταξύ του 1 και του 15, το οποίο μπορεί να εκφραστεί ως εξής:
\rho \in \mathbb{Z}, 1 \leq \rho \leq \rho_{max} = 15
Αρκετά όμως με αυτά. Ώρα να βρούμε τις ιδανικές υστερήσεις. Υπάρχουν δύο τρόποι για να διευκρινίσω τις υστερήσεις μου:
Πλέον είμαι σε θέση να εκτελέσω τον επαυξημένο έλεγχο Dickey - Fuller μέσω του πακέτου urca με την εντολή ur.df, θέτοντας τη παράμετρο lags ίσης με τη μέγιστη τιμή των υπερβάσεων p_{max} = 12. Μη ξεγελιέσται από το όνομα της παραμέτρου, δεν θέτετε τις υπερβάσεις, αλλά μέγιστο αριθμό αυτών.
types <- c("none", "drift", "trend")
# run ur.df for each case
results <- lapply(types, function(t) {
ur.df(grc_unemployment$Value, type = t, lags = 15, selectlags = "AIC")
})
result_table <- data.frame(
"Έλεγχος" = c("ADF (none)", "ADF (drift)", "ADF (trend)"),
"Στατιστική Τιμή" = c(round(results[[1]]@teststat[1], 3),
round(results[[2]]@teststat[1], 3),
round(results[[3]]@teststat[1], 3)),
"Υστέρηση" = c(results[[1]]@testreg$coefficients[,1] %>% length() - 1,
results[[2]]@testreg$coefficients[,1] %>% length() - 2,
results[[3]]@testreg$coefficients[,1] %>% length() - 3),
"Κρίσιμη Τιμή (1%)" = c(results[[1]]@cval[1,1],
results[[2]]@cval[1,1],
results[[3]]@cval[1,1]),
"Κρίσιμη Τιμή (5%)" = c(results[[1]]@cval[1,2],
results[[2]]@cval[1,2],
results[[3]]@cval[1,2]),
check.names = FALSE
)
result_table %>% reactable_custom()Με δεδομένο ότι τα αποτελέσματα του πίνακα Πίνακας 7 δεν είναι στατιστικά σημαντικά μπορώ να συμπεράνω ότι η χρονοσειρά μου δεν είναι στάσιμη. Συνεπώς η χρήση των διαφορών των παρατηρήσεων έχει νόημα και ο επανέλεγχος αυτών ως προς τη στασιμότητα. Αξίζει να υπενθυμίοσουμε ότι αυτές τις διαφορές τις έχουμε στο Σχήμα [].. και σε αυτό το σχήμα παρατηρούμε ότι εξαλείφθηκε τόσο η τάση όσο και η τάση να κυμαίνεται η παρατήρησή μας σε τιμές διάφορες τιυ μηδενός. Έτσι το πιο ακριβές σε αυτή τη περίπτωση είναι να κρίνουμε τον έλεγχο ως προς το πιο απλό μοντέλο του όπου δεν έχει ούτε σταθερό όρο (drift), ούτε και κάποιου είδους τάση. Στα αποτελέσματα για λόγους πληρότητας θα συμπεριλάβουν όλες τις περιπτώσεις.
types <- c("none", "drift", "trend")
# run ur.df for each case
results <- lapply(types, function(t) {
ur.df(diff(grc_unemployment$Value, 1), type = t, lags = 10, selectlags = "Fixed")
})
result_table <- data.frame(
"Έλεγχος" = c("ADF (none)", "ADF (drift)", "ADF (trend)"),
"Στατιστική Τιμή" = c(glue(round(results[[1]]@teststat[1], 3), "*"),
round(results[[2]]@teststat[1], 3),
round(results[[3]]@teststat[1], 3)),
"Υστέρηση" = c(results[[1]]@testreg$coefficients[,1] %>% length() - 1,
results[[2]]@testreg$coefficients[,1] %>% length() - 2,
results[[3]]@testreg$coefficients[,1] %>% length() - 3),
"Κρίσιμη Τιμή (1%)" = c(results[[1]]@cval[1,1],
results[[2]]@cval[1,1],
results[[3]]@cval[1,1]),
"Κρίσιμη Τιμή (5%)" = c(results[[1]]@cval[1,2],
results[[2]]@cval[1,2],
results[[3]]@cval[1,2]),
check.names = FALSE
)
result_table %>% reactable_custom()Αφού η στατιστική τιμή (-2.119) είναι μικρότερη της αντίστοιχης κρίσιμης τιμής (-1.95), τότε μπορούμε να απορρίψουμε την μηδενική υπόθεση μη στασιμότητας για τις διαφορές των αρχικών παρατηρήσεων της ανεργίας στην Ελλάδα, σε επίπεδο σημαντικότητας 5%.
Άλλος ένας έλεγχος στασιμότητας είναι ο έλεγχος Phillips - Perron, ο οποίος βασίζεται στον απλό Dickey Fuller. Η λογική είναι ίδια με τον ADF, μιας και η υπόθεση έχει την ίδια μορφή, αλλά διαφέρουν σημαντικά ως προς τον τρόπο που αντιμετωπίζει την αυτοσυσχέτιση των σφαλμάτων. O έλεγχος δεν προσθέτει lags για να περιορίσει την αυτοσυσχέτιση των σφαλμάτων (όπως ο ADF), αντιθέτως προσπαθεί να τροποποιήσει τη στατιστική τιμή εκτιμώντας τη μακροχρόνια διακύμανση \hat{\lambda}. Η στατιστική t του ελέγχου Dickey Fuller:
t = \frac{\hat{\rho}}{se(\hat{\rho})}
αλλάζει στο και υπολογίζει μία τροποποιημένη στατιστική τιμή του t η οποία βασίζεται στη μακροχρόνια () και τη βραχυχρόνια διακύμανση
t_{adj} = t \frac{s}{\hat{\lambda}} Όπου, t είναι η στατιστική τιμή του απλού ελέγχου Dickey-Fuller, s είναι η βραχυχρόνια διακύμανση:
s =
και τέλος, η \hat{\lambda} είναι η μακροχρόνια διακύμανση που ορίζεται ως εξής:
\hat{\lambda} =
Όπως είναι εμφανές ο μοναδικός όρος´που δεν μπορεί να διευκρινιστεί εκ των προτέρων είναι η μακροχρόνια διακύμανση καθώς το άθροισμά της εξαρτάται από έναν παράγοντα q. Δηλαδή δεν έχουμε αποφασίσει πόσα αθροίσματα θα λάβουμε υπόψιν μας στον υπολογισμό - εκτίμηση της μακροχρόνιας διακύμανσης \hat{\lambda}. Για τον παράγοντα εύρους της μακροχρόνιας διακύμανσης q υπάρχουν κάποιοι απλοί κανόνες υπολογισμού που συμπίπτουν με αυτούς των μέγιστων αριθμών υστερήσεων του ελέγχου ADF. Για μικρά δείγματα,
q_{\text{μικρό}} = \left[4\cdot \left( \frac{T}{100}\right)^{0.25} \right]
ενώ για μεγαλύτερα δείγματα:
q_{μεγάλο} = \left[12 \cdot \left( \frac{T}{100}\right)^{0.25} \right]
Ο έλεγχος υποθέσεων έχει την ίδια μορφή με τους προαναφερόμενους ελέγχους, με τη μηδενική υπόθεση να υποθέτει τη μη στασιμότητα της χρονοσειράς μας.
H_0 : \text{Η σειρά δεν είναι στάσιμη} \\ H_1 : \text{Η σειρά είναι στάσιμη}
Ο έλεγχος Phillips-Perron μπορεί να εκτελεστεί μέσω του urca πακέτου με την εντολή ur.pp. Ως παραμέτρους δέχεται την χρονοσειρά μας, τη στατιστική τιμή που μας ενδιαφέρει να υπολογίσουμε, το είδος της χρονοσειράς (τάση, χωρίς τάση). Τέλος, ένας παράγοντας που μπορεί να καθοριστεί από την εντολή μας είναι ο μέγιστος αριθμός αθροισμάτων που θα λάβω υπόψιν για την εκτίμηση της μακροχρόνιας διασποράς όπου στο πακέτο αυτή η εντολή ονομάζεται lags. Αυτό δεν είναι άλλο από το μέγεθος q στο οποίο αναφερθήκαμε προηγουμένως. Το πακέτο μας δίνει δύο επιλογές, είτε να θέσουμε εμείς τη δική μας τιμή στη παράμετρο use.lag, είτε να χρησιμοποιήσουμε μία εκ των δύο επιλογών (short ή long) στη παράμετρο lags. Αν διαβάσουμε τον κώδικα της εντολής ur.pp θα διαπιστώσουμε ότι αν επιλέξουμε τον αυτόματο τρόπο προσδιορισμό του q, μέσω των δύο επιλογών που δίνονται, οι τύποι είναι ίδιοι με αυτούς που προσδιορίστηκαν παραπάνω για τα q_{\text{μικρό}} και q_{\text{μεγάλο}}.
types <- c("constant", "trend")
# run ur.df for each case
results <- lapply(types, function(t) {
ur.pp(grc_unemployment$Value, type = "Z-tau", model = t, lags = "long")
})
result_table <- data.frame(
"Έλεγχος" = c("PP (constant)", "PP (trend)"),
"Στατιστική Τιμή" = c(round(results[[1]]@teststat[1], 3),
round(results[[2]]@teststat[1], 3)),
"Όρος αθροίσματος μακροχρ. διακύμανσης" = c(results[[1]]@lag,
results[[2]]@lag),
"Κρίσιμη Τιμή (1%)" = c(round(results[[1]]@cval[1,1], 3),
round(results[[2]]@cval[1,1],3)),
"Κρίσιμη Τιμή (5%)" = c(round(results[[1]]@cval[1,2],3),
round(results[[2]]@cval[1,2],3)),
check.names = FALSE
)
result_table %>% reactable_custom()types <- c("constant", "trend")
# run ur.df for each case
results <- lapply(types, function(t) {
ur.pp(diff(grc_unemployment$Value,1), type = "Z-tau", model = t,
lags = "long")
})
result_table <- data.frame(
"Έλεγχος" = c("PP (constant)", "PP (trend)"),
"Στατιστική Τιμή" = c(glue(round(results[[1]]@teststat[1], 3), "**"),
glue(round(results[[2]]@teststat[1], 3), "**")),
"Υστέρηση" = c(results[[1]]@lag,
results[[2]]@lag),
"Κρίσιμη Τιμή (1%)" = c(round(results[[1]]@cval[1,1], 3),
round(results[[2]]@cval[1,1],3)),
"Κρίσιμη Τιμή (5%)" = c(round(results[[1]]@cval[1,2],3),
round(results[[2]]@cval[1,2],3)),
check.names = FALSE
)
result_table %>% reactable_custom()Συνεπώς, έχουμε απόρριψη της μηδενικής (H_0) υπόθεσης, στις διαφορές των αρχικών δεδομένων μου, αλλά μη απόρριψη στις πρώτες διαφορές των παραρτηρήσεων. Αυτό σημαίνει ότι οι παρατηρήσεις τις ανεργίας στην Ελλάδα δεν ήταν στάσιμες, όμως οι διαφορές αυτών ήταν. Συνεπώς, ο έλεγχος PP επιβεβαίωσε τα ευρήματα του ελέγχου ADF.
Ένας ακόμα έλεγχος μαναδιαίας ρίζας, είναι ο έλεγχος KPSS. Αν και ο στόχος του ελέγχου είναι ο ίδιος, έχει αντίθετη λογική σε σχέση με τους προηγούμενους ελέγχους αφού η μηδενική του υπόθεση είναι η στασιμότητα, έναντι της μη στασιμότητας.
H_0 : \text{Η σειρά είναι στάσιμη} \\ H_1 : \text{Η σειρά δεν είναι στάσιμη}
results = ur.kpss(grc_unemployment$Value, type = "mu", lags = "long")
results_diff = ur.kpss(diff(grc_unemployment$Value,1), type = "mu", lags = "long")
results_diff2 = ur.kpss(diff(grc_unemployment$Value,2), type = "mu", lags = "long")
result_table <- data.frame(
"Μεταβλητή" = c("GR_Unemp", "Δ Gr_Unemp", "Δ2 Gr_Unemp"),
"Στατιστική Τιμή" = c(glue(round(results@teststat[1], 3), "**"),
round(results_diff@teststat[1],3),
round(results_diff2@teststat[1],3)),
"Όρος αθροίσματος μακροχρ. διακύμανσης" = c(results@lag,
results_diff@lag,
results_diff2@lag),
"Κρίσιμη Τιμή (1%)" = c(round(results@cval[,4], 3),
round(results_diff@cval[,4], 3),
round(results_diff2@cval[,4], 3)),
"Κρίσιμη Τιμή (5%)" = c(round(results@cval[,2],3),
round(results_diff@cval[,2], 3),
round(results_diff2@cval[,2], 3)),
check.names = FALSE
)
result_table %>% reactable_custom()Η μηδενική μου υπόθεση υπέθεσε στασιμότητα η οποία απορρίπτεται με βάση τον έλεγχο KPSS σε επίπεδο σημαντικότητας 1%. Eπιπλέον ο έλεγχος δεν είναι σε θέση να απορρίψει τη στασιμότητα στις διαφορές των παρατηρήσεων. Αυτά τα αποτελέσματα επιβεβαιώνουν και τους ελέγχους PP καθώς και τον έλεγχο ADF.
Τέλος, ένας άλλος έλεγχος για τη στασιμότητα είναι ο έλεγχος Zivott - Andrews (έλεγχος ZA). Η συγκεκριμένη στατιστική συνάρτηση διαφέρει από τις προηγούμενες μιας και λαμβάνει υπόψιν της ορισμένα σημεία (structural breaks) από τα οποία η χρονοσειρά αλλάζει συμπεριφορά. Αυτά τα χρονικά σημεία θα μπορούσαν να είναι διάφορες ημερομηνίες όπου συνέβησαν σημαντικά γεγονότα που μπορεί να επηρέασαν τη συμπεριφορά της χρονοσειράς. Στην δική μας περίπτωση αναλύουμε την ανεργία της Ελλάδας η οποία εκτοξεύτηκε μετά την οικονομική κρίση του 2009 και στα μέσα αυτής είχαμε ιστορικό υψηλό. Είναι εμφανές ότι η κρίση χρέους αποτελεί ένα χρονικό σημείο από το οποίο επηρεάστηκε η συμπεριφορά της χρονοσειράς. Η χρήση ενός εξειδικευμένων ελέγχων με διαρθρωτικές τομές (structural breaks) κρίνεται απαραίτητη για τη περίπτωσή μας.
Οι πιο γνωστοί έλεγχοι για τη στασιμότητα σε χρονοσειρές με διαρθρωτικές τομές είναι οι:
Πόσες όμως είναι οι τομές (breaks) σε μία χρονοσειρά;
breakPoints = summary(breakpoints(grc_unemployment$Value ~ 1, h = 0.15))
summaryBreakPoints = t(breakPoints[["RSS"]]) %>% as.data.frame() %>%
mutate(across(everything()), round(., 2))
highchart() %>%
hc_chart(type = "spline") %>% # spline = smooth line
hc_title(text = "Απόδοση ανά αριθμό διαρθρωτικών τομών") %>%
hc_subtitle(text = "Το με βάση το Μπεύζιανό πληροφοριακό κριτήριο μας δηλώνει την ύπαρξη δύο διαρθρωτικών τομών. Όσο μικρότερος ο αριθμός, τόσο το καλύτερο.") %>%
# custom colors (muted, professional palette)
hc_colors(c("#0072B2", "#E69F00")) %>%
# series
hc_add_series(data = summaryBreakPoints$RSS, name = "RSS", lineWidth = 3) %>%
hc_add_series(data = summaryBreakPoints$BIC, name = "BIC", lineWidth = 3) %>%
# axes
hc_xAxis(
type = "integer",
lineColor = "#e6e6e6",
tickColor = "#e6e6e6",
title = list(text = "(#) Αριθμός διαρθρωτικών τομών"),
labels = list(style = list(color = "#333333", fontSize = "12px"))
) %>%
hc_yAxis(
title = list(text = "Τιμή"),
gridLineColor = "#f0f0f0",
labels = list(style = list(color = "#333333", fontSize = "12px"))
) %>%
# tooltip
hc_tooltip(
shared = TRUE,
crosshairs = TRUE,
backgroundColor = "#ffffff",
borderColor = "#cccccc",
style = list(color = "#333333", fontSize = "12px"),
xDateFormat = "%b %Y"
) %>%
# plot options for styling
hc_plotOptions(
series = list(
marker = list(enabled = FALSE, radius = 4, symbol = "circle"),
states = list(hover = list(lineWidthPlus = 1))
)
) %>%
# legend
hc_legend(
align = "center",
verticalAlign = "bottom",
layout = "horizontal",
itemStyle = list(fontSize = "13px", color = "#333333")
) %>%
# export button
hc_exporting(enabled = TRUE)Υπάρχει αρκετά μεγάλη μείωση του Μπεϋζιανού κριτηρίου πληροφορίας (BIC) όταν μεταβαίνω από ένα μοντέλο με καμία διαρθρωτική τομή, σε ένα άλλο με δύο διαρθρωτικές τομές. Αυτό είναι μία σημαντική ένδειξη ότι η ανεργία μου όντως επηρεάστηκε από ξαφνικούς παράγοντες. Φυσικά το υποψιαζόμασταν αυτό καθώς έχουμε τη περίοδο της κρίσης που συνετέλεσε σε υψηλά ποσοστά ανεργίας. Αφού καθορίσαμε τον αριθμό των τομών, είναι η στιγμή να καθορίσουμε ποια είναι τα κομμάτια που πρέπει να αναλυθούν, κοινώς να διαπιστώσουμε ποια είναι αυτά τα εύρη ημερομηνιών που εντοπίστηκε σημαντική διαφοροποίηση στη συμπεριφορά της χρονοσειράς. Σύμφωνα με τα αποτελέσματα έχω:
data.frame("Αριθμός τομών" = c(1,2,3,4,5),
"Τομή A" = breakPoints$breakpoints[,1] %>% set_names(NULL),
"Τομή B" = breakPoints$breakpoints[,2] %>% set_names(NULL),
"Τομή Γ" = breakPoints$breakpoints[,3] %>% set_names(NULL),
"Τομή Δ" = breakPoints$breakpoints[,4] %>% set_names(NULL),
"Τομή Ε" = breakPoints$breakpoints[,5] %>% set_names(NULL),
check.names = FALSE) %>%
reactable_custom()και οι ημερομηνίες των αντίστοιχων τομών:
data.frame("Αριθμός τομών" = c(1,2,3,4,5),
"Τομή Α" = min(unemployment$TIME) %m+% months(breakPoints$breakpoints[,1]) %>% set_names(NULL),
"Τομή Β" = min(unemployment$TIME) %m+% months(breakPoints$breakpoints[,2]) %>% set_names(NULL),
"Τομή Γ" = min(unemployment$TIME) %m+% months(breakPoints$breakpoints[,3]) %>% set_names(NULL),
"Τομή Δ" = min(unemployment$TIME) %m+% months(breakPoints$breakpoints[,4]) %>% set_names(NULL),
"Τομή Ε" = min(unemployment$TIME) %m+% months(breakPoints$breakpoints[,5]) %>% set_names(NULL),
check.names = FALSE) %>%
reactable_custom()Παραπάνω κατέληξα ότι οι παρατηρήσεις τις ανεργίας δεν είναι στάσιμες, ωστόσο οι πρώτες διαφορές τους αποτελούν στάσιμη χρονοσειρά. Αυτή τη στιγμή θα ήθελα να μελετήσω ποιο μοντέλο (S)ARIMA(p, d, q) ενδείκνυται στη περίπτωσή μου. Για αυτό το λόγο θα κατασκευάσω τα διαγράμματα αυτοσυσχέτισης και μερικής αυτοσυσχέτισης για να αναγνωρίσω τις ιδανικές τιμές p, q. Τέλος, υπενθυμίζω ότι σε μοντέλο ARIMA, το d συμβολίζει την τάξη της διαφοράς, όπου στη περίπτωσή μας είναι 1.
Για τυπικούς λόγους θα ξεκινήσουμε από τις δοσμένες παρατηρήσεις τις ανεργίας. Τυπικά ήδη γνωρίζουμε ότι δεν είναι στάσιμες, αλλά υπάρχουν συγκεκριμένα μοτίβα που πρέπει να παρατηρήσουμε στα διαγράμματα αυτοσυσχέτισης προκειμένου να το επιβεβαιώσουμε. Το διάγραμμα αυτοσυσχέτισης πέφτει με έναν εξαιρετικά αργό ρυθμό που είναι ισχυρή ένδειξη μη στασιμότητας της σειράς.
acf_plot_function(grc_unemployment_diff1, min=-1, max = 1)pacf_plot_function(grc_unemployment_diff1, min=-1, max=1)Η μη στασιμότητα των δεδομένων μου που έχει χιλιοεπιβεβαιωθεί μας υποχρεώνει να λάβουμε τις πρώτες διαφορές και να λάβουμε τα αντίστοιχα διαγράμματα αυτοσυσχέτισης. Παρατηρώ μία αρκετά δυσνόητη θα έλεγα κατάσταση από τα διαγράμματα καθώς υπάρχουν σημαντικές αυξομειώσεις και στα δύο διαγράμματα που κάνουν δύσκολη την απόφαση για το κατάλληλο μοντέλο ARIMA. Συνήθως για να λάβω τα κατάλληλα μέτρα βλέπω που έχω απότομη μείωση εντός της κόκκινης περιοχής, αλλά έχω αρκετές στατιστικά σημαντικές υπερβάσεις. Μετά τη 10η υστέρηση αυτή η τάση μειώνεται ή εξαλείφεται σημαντικά, επομένως το μοντέλο μου έχει μάλλον παραμέτρους p και q κάπου ανάμεσα στο 1 και στο 10.
acf_plot_function(grc_unemployment_diff1, min=-0.1, max = 1)pacf_plot_function(grc_unemployment_diff1, min=-0.2, max=0.5)Για το απλό ARIMA μοντέλο θα πρέπει να καθορίσω τις τρεις παραμέτρους μου. Γνωρίζω ότι το d = 1, μιας και στις πρώτες διαφορές πέτυχα τη στασιμότητα. Για τον καθορισμό του p υποψιάζομαι το 1, το 6 και το 10 καθώς μετά από αυτές τις τιμές είχα μεγάλη απότομη μείωση στο διάγραμμα μερικής αυτοσυσχέτισης των πρώτων διαφορών. Τέλος, για το q υποψιάζομαι το 1, το 3, το 6 και το 10. Γενικότερα από το σχήμα δεν μπορώ να βγάλω κάποιο ξεκάθαρο συμπέρασμα για τις τιμές τους μιας και υπάρχουν αρκετά σημαντικές αυτοσυσχετίσεις τουλάχιστον μέχρι τη 10η υστέρηση. Για αυτό το λόγο πέρα από τους χειροκίνητους ελέγχους που καθορίζω τις τιμές των παραμέτρων του μοντέλου, θα αξιοποιήσω την εντολή auto.arima, ώστε αυτόμαστα να εξεταστεί όλο το εύρος τιμών που υποψιάζομαι βρίσκεται το καλύτερο μοντέλο.
auto_model <- auto.arima(grc_unemployment$Value, seasonal = TRUE, d = 1, trace = TRUE, stepwise = FALSE, max.p = 11, max.q = 11,max.order = 30)
arimaModel_1=arima(grc_unemployment$Value, order=c(0,1,1))
arimaModel_2=arima(grc_unemployment$Value, order=c(1,1,0))
arimaModel_3=arima(grc_unemployment$Value, order=c(2,1,1))
arimaModel_4=arima(grc_unemployment$Value, order=c(1,1,2))
arimaModel_5=arima(grc_unemployment$Value, order=c(5,1,5))
arimaModel_6=arima(grc_unemployment$Value, order=c(9,1,1))accuracy_table <- data.frame(
"Όνομα μοντέλου" = c("Αυτόματο μοντέλο", "Υποψήφιο Μοντέλο #1", "Υποψήφιο Μοντέλο #2", "Υποψήφιο Μοντέλο #3", "Υποψήφιο Μοντέλο #4", "Υποψήφιο Μοντέλο #5", "Υποψήφιο Μοντέλο #6"),
Model = c("ARIMA(5,1,6)", "ARIMA(0,1,1)", "ARIMA(1,1,0)", "ARIMA(2,1,1)",
"ARIMA(1,1,2)", "ARIMA(5,1,5)", "ARIMA(9,1,1)"),
AIC =c(auto_model$aic, arimaModel_1$aic, arimaModel_2$aic, arimaModel_3$aic,
arimaModel_4$aic, arimaModel_5$aic, arimaModel_6$aic),
check.names = FALSE
) %>%
mutate(AIC = round(AIC, 2))
accuracy_table %>% reactable_custom(head_max = 8)Αν και βρήκα ένα κάποια καλά μοντέλα με αρκετά καλή απόδοση μόνος μου, με βάση τον αυτόματο έλεγχο βρήκα ότι το ARIMA(5,1,6) περιγράφει καλύτερα τη χρονοσειρά μου. Αυτό τεκμαίρεται συγκρίνοντας το κριτήριο πληροφορίας των διάφορων μοντέλων ARIMA, όπου στο βέλτιστο μοντέλο έχω το μικρότερο αριθμό AIC, με 271 μονάδες. Αυτό είναι 6 μονάδες χαμηλότερα της δεύτερης επιλογής μου και προφανώς αποτελεί την κύριά μου επιλογή για τον καθορισμό προβλέψεων.
Έφτασε η στιγμή όλα αυτά να βγάζουν κάποιο νόημα. Βρήκα ποιο είναι το καλύτερο μοντέλο, χρησιμοποιωντας το θα προβλέψω το ύψος της ανεργίας στην Ελλάδα για τους επόμενους 12 μήνες. Τα δεδομένα μου σταματάνε στον Αύγουστο του 2022 (ανεργία 12.2%), άρα η όποια εκτίμηση θα γίνει για τους μήνες Σεπτέμβριο του 2022 μέχρι και τον Αύγουστο του 2023.
s = forecast(auto_model,20)
month = seq(max(grc_unemployment$TIME) %m+% months(1), # start at next month
by = "1 month",
length.out = 20)
highchart() %>%
hc_title(text = "Πρόβλεψη ανεργίας για το επόμενο έτος") %>%
hc_subtitle(text = "Η ανεργία θα βαίνει μειούμενη για το επόμενο διάστημα με βάση το μοντέλο μου.") %>%
hc_xAxis(type = "datetime") %>%
hc_yAxis(title = list(text="Ανεργία (%)")) %>%
hc_add_series(name="Ιστορικά δεδομένα",
data = list_parse2(data.frame(x = datetime_to_timestamp(grc_unemployment$TIME), y = grc_unemployment$Value)),
type = "line"
) %>%
hc_add_series(name = "Πρόβλεψη",
data = list_parse2(data.frame(
x = datetime_to_timestamp(month),
y = round(s$mean, digits=2)
)),
type = "line",
color = "#d62728") %>%
hc_add_series(name = "Confidence Interval",
data = list_parse2(data.frame(
x = datetime_to_timestamp(month),
low = round(s$lower[, "80%"], digits=2),
high = round(s$upper[, "80%"], digits=2)
)),
type = "arearange",
color = hex_to_rgba("#d62728", 0.2),
linkedTo = ":previous")Και ο αντίστοιχος πίνακας με τις προβλέψεις:
s %>%
as.data.frame() %>%
dplyr::mutate(Date = seq(max(grc_unemployment$TIME) %m+% months(1), # start at next month
by = "1 month",
length.out = 20)) %>%
setNames(c("Πρόβλεψη", "Κάτω διάστημα 80%", "Άνω διάστημα 80%", "Κάτω διάστημα 95%", "Άνω διάστημα 95%", "Μήνας")) %>%
dplyr::relocate("Μήνας", .before = "Πρόβλεψη") %>%
dplyr::mutate(across(where(is.numeric), ~ round(.x, 1))) %>%
reactable_custom(head_max = 7)Με δεδομένο το διάγραμμα αλλά και το πίνακα των προβλέψεων του μοντέλου με τη καλύτερη επίδοση, συμπεραίνω ότι η μείωση της ανεργίας στην Ελλάδα θα συνεχιστεί και την επόμενη περίοδο (τους επόμενους έξι μήνες). Πιο συγκεκριμένα, η ανεργία στην Ελλάδα τον Φεβρουάριο του 2023 θα είναι 11.4%, κυμαίνεται μεταξύ του 10.1% (9.4%) και 12.8% (13.5%) σμε πιθανότητα 80% (95%).
Φωτογραφία από Rosy / Bad Homburg / Germany από το Pixabay
@online{2022,
author = {, stesiam},
title = {Προβλέποντας την Ανεργία στην Ελλάδα},
date = {2022-10-22},
url = {https://stesiam.com/el/posts/forecasting-greek-unemployment/},
langid = {el}
}