סדרות עיתיות הן סוג אחר של בעיות ממה שטיפלנו בו עד כה. במובן מסוים הן שייכות לבעיות “supervised”, אך בסדרות עיתיות המטרה אינה לחזות ערך \(y\) בודד, אלא לחזות סדרה של ערכים בהתבסס על סדרה שקדמה להם.

בסדרות עיתיות אנחנו מנסים להבין תופעות כגון עונתיות, ומגמות, כדי לנסות לחזות דברים שונים כגון:

ביחידה זו נשתמש בנתונים מנורמלים של חיפושי גוגל, ננתח, וננסה לפרש סדרות עיתיות עבור מונחי חיפוש כגון:

כמו כן נשתמש בנתוני פתיחה וסגירה של מניית Tesla

פרק זה לא ידון רבות בתיאוריה של סדרות עיתיות, אך הקורא המעוניין מופנה לספר הבא, ספר המציג קצת תיאוריה והרבה פרקטיקה בנושא של סדרות עיתיות:

Forecasting: Principles and Practice, OTexts (2018?). J. Hyndman and G. Athanasopoulos, Available online: https://otexts.org/fpp2/, (feched 12/11/2018).

מה היא סדרה עיתית (Time Series)

סדרה עיתית היא אוסף של תצפיות בעלי יחס סדר (התצפיות באות אחת אחרי השנה), במרווחי זמן קבועים (לדוגמה פעם בחודש).

חיזוי סדרות עיתיות מעניין אותנו בכמה רמות:

תהליך העבודה עליו דיברנו תופס גם בעבודה עם סדרות עיתיות (יבוא הנתונים, סידור הנתונים, ויז’ואליזציה, מידול, חיזוי, והצגת הממצאים בפני מקבלי החלטות).

הפעולה הראשונה בעבודה עם סדרה עיתית היא קריאת הסדרה והכנסתה לפורמט של סדרה עיתית ts.

library(tidyverse)
zombie_raw <- read_csv("data-files/zombieTimeline.csv", skip = 3,
                       col_names = c("Month", "search_volume"))
glimpse(zombie_raw)
## Observations: 179
## Variables: 2
## $ Month         <chr> "2004-01", "2004-02", "2004-03", "2004-04", "200...
## $ search_volume <int> 6, 6, 8, 6, 6, 6, 6, 6, 7, 10, 7, 7, 6, 7, 7, 7,...
# most of the time the data format will not be so nice, and you will have to consolidate data.
# in this case we already get a year-month data ready. Just need to turn it into date.
zombie_raw <- zombie_raw %>%
  mutate(Month = lubridate::ymd(paste0(Month, "-01")))

zombie_ts <- ts(zombie_raw$search_volume, start = 2004, frequency = 12)
# frequency is super important, it will be used later on in our forecast algorithms.
# Frequency can have values:
# 1=Annual
# 4=Quarterly
# 12=Monthly
# 52=Weekly
# 60=Seconds
# 24=Hourly
# 7=Daily
# You get the drift...

# You can see that R prints the data in a nice table:
zombie_ts
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2004   6   6   8   6   6   6   6   6   7  10   7   7
## 2005   6   7   7   7   8  10  14   9   9  13  12  11
## 2006  10  11  12  11   9  11  10  11  11  15   9   9
## 2007   9   9   9  12  14  13  13  16  17  20  14  15
## 2008  14  13  15  17  17  16  14  14  15  25  21  21
## 2009  21  21  22  23  25  23  22  27  27  44  28  25
## 2010  21  23  22  23  25  22  22  23  24  41  39  35
## 2011  37  36  35  34  45  34  34  34  37  64  37  34
## 2012  30  32  36  35  69 100  43  36  41  67  47  44
## 2013  34  37  35  35  32  34  35  34  37  65  33  28
## 2014  25  28  29  26  27  28  27  27  30  62  28  25
## 2015  25  26  26  25  25  24  24  26  28  56  30  26
## 2016  23  24  23  26  22  22  21  21  25  40  19  19
## 2017  18  18  17  17  20  18  18  16  18  38  18  16
## 2018  18  17  20  20  19  18  18  23  18  34  20
# Here is the timeplot describing search volume of zombie
ggplot(zombie_raw, aes(x = Month, y = search_volume)) + 
  geom_line() + 
  ggtitle("Google search volume for the word Zombie") + 
  ylab("Google coefficient for search volume") + 
  scale_x_date(date_breaks = "6 months") + 
  theme(axis.text.x = element_text(angle = 90))

התרשים הבא מציג את ה“עונות השונות” (כל שנה מהווה עונה). הוא מוצג כגרף ggplot2, אך בעצם מופעל מחבילת forecast שבה נשתמש בהמשך הפרק לויז’ואליזציה וכמו כן לבניית מודלים.

library(forecast)
## Warning: package 'forecast' was built under R version 3.4.4
ggseasonplot(zombie_ts, year.labels = TRUE, year.labels.left = TRUE) + 
  ylab("Search volume (google trend coefficient)") + 
  ggtitle("Seasonal plot: the search for Zombies")

# polar view
ggseasonplot(zombie_ts, polar = TRUE) + 
  ylab("Search volume (google trend coefficient)") + 
  ggtitle("Seasonal plot: the search for Zombies")

מעיון בתרשימים עולה תמונה של עליה של החיפושים של המונח עד שנת 2011, והחל משנת 2012, עונתיות שנתית (שיא בחיפושים סביב חודש אוקטובר). למה?

כמו כן, יש שיא מקומי בחודש יוני של שנת 2012.

כמו כן, ישנה היחלשות מסוימות בעוצמת החיפושים (הפיקים הולכים וקטנים בשנים האחרונות).

מקובל גם לבחון את הקורלציה בין תצפיות לאורך זמן.

ggAcf(zombie_ts, lag = 48)

התרשים מציג קורלציה גבוהה יחסית בין חודש לחודשים שקדמו לו, וגם בין חודש לבין אותו החודש שנה לפני (12 חודשים אחורה). מעבר לקשר של 12 חודשים אחורה, הקשר הולך ונחלש, וכמעט נעלם לאחר 25 חודשים.

כעת נרצה להשתמש בפקודה שתתאים מודל שיתפוס את העונתיות מחד, ואת המגמה מאידך.

המודלים הנאייביים ביותר הם:

autoplot(zombie_ts) +
  autolayer(meanf(zombie_ts, h=24),
    series="Mean", PI=FALSE) +
  autolayer(naive(zombie_ts, h=24),
    series="Naive", PI=FALSE) +
  autolayer(snaive(zombie_ts, h=24),
    series="Seasonal naive", PI=FALSE) +
  ggtitle("Forecasts for zombie searches") +
  xlab("Year") + ylab("Zombie searches") +
  guides(colour=guide_legend(title="Forecast")) 

# With confidence intervals using the seasonal naive model:
autoplot(zombie_ts) +
  autolayer(snaive(zombie_ts, h=24),
    series="Seasonal naive", PI=TRUE) +
  ggtitle("Forecasts for zombie searches - with confidence intervals") +
  xlab("Year") + ylab("Zombie searches") +
  guides(colour=guide_legend(title="Forecast"))

מעבר למודלים הבסיסיים שהוצגו לעיל, ישנן עוד שתי משפחות של מודלים מתקדמים יותר:

משפחת המודלים מסוג החלקה אקספוננציאלית מבוססים על הרעיון שכל תצפית עתידית מבוססת על ממוצע משוקלל (החלקה) של תצפיות קודמות. ייתכן שממוצע זה מורכב מרכיבים מיוחדים (כגון מגמה, או עונתיות).

\[ \hat{y}_{T+1|T}=\alpha y_T + \alpha(1-\alpha)y_{T-1}+\alpha(1-\alpha)^2y_{T-1}+\ldots \]

הפרמטר \(\alpha\) הוא פרמטר ההחלקה.

פקודת ets מבצעת התאמה של הסדרה העיתית לאחד מ-18 משפחות הכוללות סוגים שונים של מודלים דומים (טרנדים, עונתיות, ושגיאות חיבוריות וכפליות)

zombie_ets <- ets(zombie_ts)
summary(zombie_ets)
## ETS(M,N,M) 
## 
## Call:
##  ets(y = zombie_ts) 
## 
##   Smoothing parameters:
##     alpha = 0.9382 
##     gamma = 3e-04 
## 
##   Initial states:
##     l = 7.9029 
##     s = 0.9539 1.0513 1.6211 0.9665 0.9472 0.9522
##            1.0004 0.9908 0.8864 0.8987 0.8669 0.8647
## 
##   sigma:  0.1609
## 
##      AIC     AICc      BIC 
## 1364.966 1367.911 1412.777 
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.09152998 5.943424 2.726626 -0.756665 10.63041 0.4627506
##                     ACF1
## Training set -0.09240921
autoplot(zombie_ets)

המודל המיטבי הינו מודל בעל שגיאה כפלית, עונתיות כפלית, וללא טרנד. שימו לב שהמודל עשוי להשתנות, כתלות באורך הסדרות (התצפיות) שמזינים לתוכו.

zombie_ets_3yrs <- ets(tail(zombie_ts, 36))
summary(zombie_ets_3yrs)
## ETS(M,N,A) 
## 
## Call:
##  ets(y = tail(zombie_ts, 36)) 
## 
##   Smoothing parameters:
##     alpha = 0.3661 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 24.4847 
##     s = -1.3401 16.7692 -1.3084 -1.0286 -2.6781 -2.1174
##            -0.9071 -0.6359 -1.2091 -1.9143 -1.9626 -1.6675
## 
##   sigma:  0.1159
## 
##      AIC     AICc      BIC 
## 205.7536 229.7536 229.5064 
## 
## Training set error measures:
##                      ME     RMSE     MAE       MPE    MAPE      MASE
## Training set -0.3272774 1.953499 1.64154 -2.094269 7.86075 0.4863822
##                     ACF1
## Training set -0.05379723
autoplot(zombie_ets_3yrs)

כדי לספק את המודל, פונקציית ets מבצעת Best Fit בהתבסס על מדד log-liklihood.

באמצעות פונקציית predict, ניתן לספק תחזית.

zombie_prediction <- predict(zombie_ets, h = 36)
autoplot(zombie_prediction)

כדי לבחון את מידת הדיוק של המודלים, מדווחת הפונקציה ets על מספר סוגי שגיאה.

תרגיל

  1. בחרו אחד מקבצי הנתונים שמוצגים בתחילת הפרק, שאינם קשורים בזומבים, כלומר
    1. נתוני חיפוש על בריכות שחייה, או
    2. נתוני חיפוש על Data science, או Data mining, או Machine learning, או
    3. נתונים על מניות חברת Tesla.
  2. הציגו את הנתונים בתרשים מתאים - האם אתם מזהים דפוס עונתי או מגמה?
  3. התאימו מודל ets, אילו דפוסים מזהה המודל?
  4. בכל התרגילים שעסקו ב-Supervised learning השתמשנו בחלוקה ל-Train/test.
    1. חשבו - מה החלוקה המקבילה בבעית סדרות עיתיות? כיצד ניתן לחלק סדרה ל-Train/test ולבחון את שיטת החישוב?
    2. בצעו חלוקה ל-Train/test והעריכו את ביצועי המודלים שפיתחתם לעומת ביצועיהם של מודלים נאיביים. לצורך חישובים אלו, באפשרותכם להשתמש בפונקציה accuracy מתוך חבילת forecast.