01-intro-to-tsibble.Rmd

# 시계열 첫걸음, 관련 패키지의 이해 {#intro}

## `tsibble`, `feasts` 그리고 `fable` 패키지

시계열과 관련 된 R 패키지는 정말 많이 있다. 그 중에서 가장 유명한 패키지는 `forecast` 패키지이다. `fpp` 2판에서 메인 패키지로 활용되었던 `forecast` 패키지는 예측성능이 좋기로 유명하다. 하지만 `forecast`는 base 패키지를 기반으로 짜여져있는 패키지라서 `fpp` 3판에서 사용하게 될 패키지는 `tsibble`과 `feasts`, 그리고 `fable` 패키지를 주로 사용하게 될 예정이다. 세가지 패키지 모두 `tidyverse`스러운 코딩과 잘 맞는 함수들을 지원하는 것이 특징이다.

```{block, type='rmdwarning'}
### 주의하기

`forecast` 패키지는 다른 언어를 사용하는 커뮤니티에도 알려질 정도로 잘 짜여진 패키지이다. 꼭 한번은 사용법을 알아둘 것을 추천한다. [`fpp` 2판의 한글판](https://otexts.com/fppkr/)을 봐두는 것도 나쁘지 않는 선택이다.
```

* `tsibble` 패키지는 `tidyverse`의 대표 객체인 `tibble`을 시계열 데이터를 잘 다룰수 있도록 조정한 `tsibble` 객체를 지원한다. 

* `feasts` 패키지는 `Feature Extraction And Statistics for Time Series`의 줄임말로 시계열 데이터에서 특징을 뽑아내는데 사용된다. acf 함수라던지 lag 그래프를 그리는데에 사용된다.

* `fable`은 시계열 모델링에 관한 함수들이 모여져 있는 패키지이다. ARIMA, ETS와 같은 모델링 함수들을 가지고 있다.

## `fpp3` 패키지

`fpp3` 패키지는 Hyndman 교수와 Athanasopoulos 교수가 지은 Forecasting: principles and practice 3판에 딸림 R 패키지라고 생각하면 된다. `fpp3` 패키지를 로드시키면 앞에서 말했던 `tsibble`과 `feasts`, 그리고 `fable` 패키지가 같이 로드된다.

```{r message=TRUE}
library(fpp3)
```

## tsible 객체

```{r}
y <- tsibble(
  year = 2015:2019,
  observation = c(123, 39, 78, 52, 110),
  index = year
)
y
```

`tibble` 객체를 확장 시킨 `tidyverse`용 데이터 프레임. y변수를 선언할 때 어떤 열을 `index`로 할 것인지 선언해 줌. `[1Y]`라는 것을 보고 관찰값이 1년 단위 관측되었다는 것을 알 수 있음.

`tsibbledata` 패키지에 들어있는 `olympic_running` 데이터를 살펴보자.

```{r}
olympic_running
```

위의 결과값을 보면 key 변수가 2개 (Length, Sex) 로 설정이 되어있고, `[14]`라는 표시가 보인다. 이것은 두 개의 변수가 갖는 unique한 값들로 조합 할 수 있는 경우의 수가 14가지 경우가 된다는 것을 의미한다.

다음을 살펴보자.

```{r}
olympic_running %>% 
  distinct(Length)

olympic_running %>% 
  distinct(Sex)
```

## csv 파일로부터 `tsibble` 객체 만들기 연습

```{r}
prison <- readr::read_csv("https://OTexts.com/fpp3/extrafiles/prison_population.csv")
head(prison)
```

`prison` 변수에는 csv 파일로부터 불러온 시계열 정보가 들어있다. 날짜정보가 들어있는 첫번째 열을 인덱스로, count 정보를 제외한 나머지 열들을 key로 설정하여 tsibble로 바꿔보도록 하자.

```{r}
prison %<>% 
  janitor::clean_names() %>% 
  mutate(quarter = yearquarter(date),
         .keep = "unused") %>% 
  as_tsibble(key = state:indigenous,
             index = quarter)
prison
```

분기별 데이터를 `yearquarter` 사용해서 설정해줬던 것 처럼 다음 표과 같이 여러 시간 간격을 설정할 수 있는 함수를 `tsibble` 패키지에서는 제공한다. 자신의 데이터에 맞는 시간 간격을 설정해서 사용하자.

```{r echo=FALSE}
tibble::tribble(
   ~Frequency,                  ~Function,
     "Annual",                "start:end",
  "Quarterly",            "yearquarter()",
    "Monthly",              "yearmonth()",
     "Weekly",               "yearweek()",
      "Daily",         "as_date(), ymd()",
  "Sub-daily", "as_datetime(), ymd_hms()"
  ) %>% knitr::kable(caption = "tsibble 객체에 부여할 수 있는 시간 간격 종류들") %>%
  kableExtra::kable_styling(full_width = FALSE)
```



## Time plots

시계열에서 제일 기본적인 그래프이다. 시간에 따라서 우리가 보고싶은 반응 변수값을 그려준다. 앞에서 정의한 `prison` 시계열 데이터에서 특정 그룹의 시간에 따른 죄수 수의 변화를 그려보자.

```{r}
prison %>% 
  filter(state == "ACT", gender == "Female",
         legal =="Remanded", indigenous == "ATSI") %>% 
  autoplot(count) +
  labs(title = "시간에 따른 죄수 인구의 변화",
       y = "죄수 (명)",
       x = "시간 (분기)")
```

## 시계열 자료 구성 요소

* 추세 (Trend): 전반적인 방향성의 존재 유무
* 계절성 (Seasonal): 계절에 따른 반복유무 (고정 빈도)
* 주기성(cycle): 고정된 빈도가 아닌 증가, 감소 형태

## Seasonal plot

`vic_elec` 데이터는 호주의 빅토리아 주의 전력 수요량을 30분 단위로 기록한 데이터이다. 

```{r}
head(vic_elec)
```

이 자료에서 하루를 기준으로 전력 수요가 어떠한 패턴으로 움직이는지 보고싶을 때 계절성 그래프 (Seasonal plot)가 유용함.

```{r fig.cap="하루 기준 전력 수요량의 변화 패턴"}
vic_elec %>% 
  janitor::clean_names() %>% 
  filter(year(time) == 2012) %>% 
  gg_season(demand, period = "day") +
  theme(legend.position = "none") +
  labs(y="MWh", title="Electricity demand: Victoria")
```

한 주를 기준으로 그려보면 새로운 패턴을 발견할 수 있다.

```{r fig.cap="한주 기준 전력 수요량의 변화 패턴"}
vic_elec %>% 
  janitor::clean_names() %>% 
  filter(year(time) == 2012) %>% 
  gg_season(demand, period = "week") +
  theme(legend.position = "none") +
  labs(y="MWh", title="Electricity demand: Victoria")
```

```{block note-dataloader, type='rmdnote'}
### 혼자 학습하기

연도를 기준으로 계절성 그래프를 그려보고, 계절에 따른 전력 수요량 패턴을 해석해보자.
```

## Scatter plots

전력 수요량 데이터의 경우 전력 수요량을 나타내는 `demand` 변수가 존재하고, 이에 대응하는 `temperature` 변수가 존재한다. 두 개의 변수를 독립적으로 시계열 데이터로 보고 시간 그래프를 그려보면 다음과 같을 것이다.

```{r}
vic_elec %>% 
  janitor::clean_names() %>% 
  filter(year(time) == 2012) %>% 
   autoplot(demand) +
  labs(y = "GW",
       title = "Half-hourly electricity demand: Victoria")
```


```{r}
vic_elec %>% 
  janitor::clean_names() %>% 
  filter(year(time) == 2012) %>% 
   autoplot(temperature) +
  labs(y = "Degrees Celsius",
       title = "Half-hourly electricity demand: Victoria")
```

하지만, 위의 두 그래프를 시간 변수를 무시한 채 산점도를 그려보면 어떻게 될까? 다음과 같이 산점도를 그릴 경우 온도와 전력 수요량 간의 새로운 패턴을 발견 할 수 있다.

```{r}
vic_elec %>% 
  janitor::clean_names() %>% 
  filter(year(time) == 2012) %>% 
  ggplot(aes(x = temperature, y = demand)) +
  geom_point() +
  labs(x = "Temperature (degrees Celsius)",
       y = "Electricity demand (GW)")
```

전력수요량의 경우 온도가 낮은 경우와 높은 경우 증가하는 2차 곡선 형태의 패턴을 보인다. 이것은 겨울철 난방과 여름철 냉방을 위한 전력 수요에 기인한 것으로 보인다.

```{r}
prison %>%
  group_by(state) %>%
  summarise(count = sum(count)) %>% 
  ggplot(aes(x = quarter, y = count)) +
  geom_line() +
  facet_grid(vars(state), scales = "free_y") +
  labs(title = "Australian prison population by states",
       y= "population")
```


## Lag plots

다음은 Lag plots을 공부해보자. 이것을 이해하기 위해서는 먼저 Lag라는 것은 무엇인가를 이해해야 하는데, 사실 lag라는 단어는 대한민국 모든 사람이 알고 있는 단어이다. 왜냐하면 컴퓨터 렉 걸렸어! 라고 할 때 렉이 바로 영어로 lag이기 때문.

Lag는 뒤쳐지다라는 뜻인데, 어떤 시계열 벡터가 있을 때, 특정 시점을 기준으로 이전 시점을 가리킨다. 쉽게 생각해서 시차라고 생각해 볼 수 있다.

예를 들어, 다음과 같은 숫자들이 있다고 하자.

```{r echo=FALSE, results='asis'}
x <- cumsum(1:10)
glue::glue(str_c(x, collapse = ", "))
```

위의 숫자들을 시계열 벡터 `x`라고 할 때, `x` 벡터에 대응하는 `lag 1`, `lag 2` 벡터는 다음과 같다.

```{r}
x <- c(1,3,6,10,15,21,28,36,45,55)
lag_x_1 <- lag(x, n = 1)
lag_x_2 <- lag(x, n = 2)
```

```{r echo=FALSE}
lab_tab <- tibble(x = x, lag1 = lag_x_1, lag2 = lag_x_2)
kable(lab_tab, caption = "x에 대한 lag 1, lag 2 벡터") %>%
  kableExtra::kable_styling(full_width = TRUE)
```

Lag 플랏은 시계열 벡터와 각 lag 값에 대응하는 벡터를 사용하여 산점도를 그려준다. 이 그래프를 사용하면 좋은 점은 데이터의 주기성을 다시 확인 할 수 있다는 점이다. 다음은 호주의 분기별 맥주 생산량 정보를 담고 있는 `aus_production` 데이터이다.

```{r echo=FALSE}
library(DT)
aus_production %>% 
  mutate(Quarter = as.Date(Quarter)) %>% 
  datatable(caption = '호주의 맥주 분기별 생산량')
```

lag plot은 다음과 같이 맥주 생산량 데이터의 계절성을 체크하는데 유용하게 사용할 수 있다. lag 4와 lag 8을 살펴보면 대각선을 따라서 점들이 분포하는 것을 알 수 있다. 반대로 lag 2와 lag 4의 경우에는 좌에서 우로 뻗어 내리는 대각선 형태로 점들이 분포한 것을 알 수 있다.

```{r}
aus_production %>%
  janitor::clean_names() %>% 
  filter(year(quarter) >= 2000) %>% 
  gg_lag(beer, geom = "point") +
  labs(x = "lag(Beer, k)")
```

슬기로운 통계생활에서 진행하는 [R을 사용한 기초 통계 수업](https://www.youtube.com/playlist?list=PLKtLBdGREmMnLbQnqGEfpCBtkGj2g_d-B) 을 충실히 이수한 분이라면 특정 lag에 따라서 상관계수의 절대값이 크게 나올 것이라는 생각을 할 수 있을 것이다. 이러한 현상과 관련이 있는 개념이 바로 시계열에서의 자기 상관성이다.

## Autocorrelation

기초 통계 시간에 배웠던 상관계수를 구하는 식을 살펴보자. 다음과 같이 $\underline{x}=\{x_1, x_2, ..., x_n\}$와 $\underline{y}=\{y_1, y_2, ..., y_n\}$의 $n$개의 표본이 주어졌을 때, 두 벡터의 표본 상관계수를 구하는 식은 다음과 같다.

$$
{\displaystyle r_{xy}={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{{\sqrt {\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}{\sqrt {\sum _{i=1}^{n}(y_{i}-{\bar {y}})^{2}}}}}}	
$$

자기 상관성은 특정 lag 값에 대응하는 벡터와 자기 자신과의 상관계수를 측정한 값으로 이해할 수 있다. lag $k$에 대한 (표본) 자기 상관계수 $r_k$는 다음과 같이 계산 할 수 있다.

$$
r_{k}=\frac{\sum_{t=k+1}^{T}\left(x_{t}-\overline{x}\right)\left(x_{t-k}-\overline{x}\right)}{\sum_{t=1}^{T}\left(x_{t}-\overline{x}\right)^{2}}
$$

분자에 위치한 summation기호는 $k+1$ 부터 시작하고, 분모는 $1$부터 시작하는 것에 주의하자. 앞에서 정의한 `x` 벡터의 $r_2$를 계산하면 다음과 같다.

```{r}
# ACF 값 직접 구하기: r_2
sum((x - mean(x)) * (lag(x, n = 2) - mean(x)),
    na.rm = TRUE) / (var(x) * (length(x) - 1))
```

위와 같은 자기 상관계수 값은 `k`값이 주어지면 값이 나오는 함수 형태로 볼 수 있다. 우리는 이것을 자기상관함수(autocorrelation function)라고 부른다. 자기상관함수 값은 `ACF` 함수를 사용해 다음과 같이 구할 수 있다.

```{r}
my_x <- tsibble(value = x,
               time = 1:length(x),
               index = time)
my_x %>% 
  ACF(value, lag_max = 9)
```

이론상으로 자기상관계수는 주어진 벡터의 길이보다 하나 작은 $n-1$ 시점까지 자기 상관계수를 구할 수 있는데, `autoplot()`은 이 값들을 각 시점마다 바차트 형태로 표현해 준다.

```{r}
my_x %>% 
  ACF(value, lag_max = 9) %>% 
  autoplot() +
  labs(title = "x 벡터의 자기상관함수",
       x = "Time lags",
       y = "Autocorrelation function")
```

이러한 자기상관함수 그래프는 시계열 데이터의 특징을 파악하는데에 중요한 역할을 하게 된다. 또한 , 이러한 특성을 잘 잡아내기 때문에 ACF 함수의 형태를 기준으로 시계열 데이터를 분류하는데에도 사용된다.