교차표에서 효과를 추정하는 방법은 아래 3가지가 있다.

교차표를 본다는 것은 범주형으로 이루어진 X,Y 의 변수간의 연관성을 파악하고 싶다는 의미이다. 

1) 비율의 차를 이용한 방법

2) 비율의 비를 이용한 방법

3) 오즈비를 이용한 방법 

흡연과 폐암의 관계

1) 비율의 차를 이용하는 방법

흡연자중 폐암 비율 = 90/1000 = 0.09

비흡연자중 폐암 비율 = 10/1000 = 0.01 

비율의 차 = 0.08 (risk difference)

2) 비율의 비를 이용하는 방법

0.09 / 0.01 = 9 (relative risk)

3) 오즈비를 이용하는 방법

(90/100) / (10/990) =  9.79 (odds ratio)

또는 p1 / (1-p1) / p2 / (1-p2) = ((90/1000) / (1-90/1000)) / ((10/1000) / (990/1000)) = 9.79 (odds ratio)

각 효과 추정 방법의 특징 및 장단점

비율의 비와 오즈비는 어떤 treatment의 효과(effect) 를 설명할 때 좋다. 비율의 비는 특히, 설명할 때 좋다. 

- 예를 들면, 어떤 위험인자가 질병에 미치는 영향이 있는지를 설명할 때는 비율의 비를 활용하는 것이 좋다. 

- 흡연의 reltavie risk 가 3이라는 말은 흡연을 하면 폐암 발생 위험을 3배 높인다고 해석할 수 있다. 

비율의 차는 전체 모수에서의 impact 를 설명할 때 좋다.

- 어떤 요인 A의 risk difference 는 10% 인데 relative risk 는 2라고 하자.

- 어떤 요인 B 의 risk difference 는 1%인데 relative risk 는 10이라고 하자.

- 이 때, 요인 B 의 effect size는 더 크지만, 실제 요인의 중요도는 A가 더 클 수 있다.  

오즈비는 y=1의 비율이 적을 때, 상대위험도와 값이 유사하다. 

- 만약의 y가 폐암과 같이 질병인 경우, P(Y=1) 은 유병률이다. 

- 즉, 유병률이 작은 질병의 경우 오즈비를 relative risk 처럼 해석할 수 있다.  

오즈비는 샘플이 불균형하게 추출한 경우에도 사용할 수 있는 지표이다. 

- 비율의 차 또는 비율의 비는 샘플링 바이어스의 영향을 받는다. 

- 만약, 흡연자 100명, 비흡연자100명을 선정해서 폐암여부를 비교할 때 비율의 비(relative risk) 에는 bias 가 생긴다. 이는 모집단에서 계산한 값과 차이가 생긴다는 의미이다. 

- 그러나, 오즈비의 경우 모집단에서의 값과 오즈비와 비교하여 bias 가 없게 된다. 

- 왜 오즈비는 샘플링 영향이 없는지 관련해서는 이 포스팅을 참고할 수 있다. 

ratio 에 로그를 취한 값은 유용하다.  

- 비율의 비 또는 오즈비는 매우 skew 된 값이다. 

- ratio는 0에서 무한대의 값을 갖는다.

- ratio 에 log 를 취해주면 -무한대~ +무한대의 값을 갖게 된다. 

- 만약, 어떤 A약의 효과가 B약의 효과보다 1.5배 있다 라는 것을 반대로 말하면 B 약의 효과가 비 약의 효과보다 1/1.5배 = 0.67배 있다라는 것이다. 그러나, 1.5배와 0.67배가 한눈에 역수 관계에 있다는 것을 알기 어렵다. 1.5배는 1로부터 0.5 떨어져 있고, 0.67은 1로부터 0.33 떨어져있다. 만약, 1.5와 0.67에 log를 취해주면, 각각 0.405, -0.405로 나오게 되어, 역관계에 있다는 것을 바로 확인할 수 있다. 

효과가 유의미한지 보려면 어떻게 할까? 

-> 신뢰구간을 보고, 이 값이 0을 포함하지 않으면 유의미하다고 판단할 수 있다.

비율의 차의 신뢰구간

$$ p_1 = n_{11}/n_1 $$

$$ p_2 = n_{21}/n_2 $$

p1,p2비율의 차의 standard error (s.e) 는 아래와 같다.  

$$ \sigma = \sqrt{\frac{p_1(1-p_1)}{n_1} + \frac{p_2(1-p_2)}{n_2}} $$

따라서 95% 신뢰구간은 1.96*±s.e 이다. 

오즈비의 신뢰구간

오즈비의 신뢰구간을 구하기 위해서는 로그 오즈비를 통해서 구하는 것이 좋다. 

로그 오즈비의 standard errer (s.e)는 아래와 같다. 

$$ \sigma = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} }  $$ 

따라서 로그 오즈비의 95% 신뢰구간은 1.96*±s.e이 이며, 이를 오즈비의 신뢰구간으로 변환하기 위해서는 exponential 을 취해주면 된다. 따라서 오즈비의 95% 신뢰구간은 exp(1.96*±s.e) 이다.

X-> Y 의 인과적 관계 파악을 위해, 간단한게 심플 회귀 분석을 진행할 수 있다. 만약 X 와 correlation 이 있고, Y 의 determinants 인 Z 라고 하는 변수가 보정되지 않는다면, omitted variable bias 가 발생한다. 이러한 상황에서 omitted variable bias 의 방향은 다음과 같이 알 수 있다. 

1) Z->X 에 영향을 주는 방향 

2) Z->Y 에 영향을 주는 방향

1) 2) 를 곱하면 이것이 bias 의 방향이 된다. 

예를 들어, 소득(X)이 의료비 지출(Y)에 주는 영향을 파악하려고 한다. 이 때, 건강 상태(이를 개인이 갖고 있는 질병의 갯수라고 하자) 를 보정하지 않으면, omitted variable bias 가 발생하게 된다. 

질병의 개수는 소득에 negative effect 이다. 질병의 개수가 증가할 수록 소득은 감소한다.

질병의 개수는 의료비 지출에 positive effect 이다. 질병의 개수가 증가할 수록 의료비 지출은 증가한다. 

1) 2) 를 곱하면 negative 가 되기 때문에 bias 의 방향은 negative 가 된다. 따라서 건강 상태를 변수로 포함하지 않고 소득과 의료비 지출의 관계를 파악하여 나온 회귀 계수는 underestimate 이 되었다고 볼 수 있다. 만약, 동일한 건강상태에 있는 사람들만을 대상으로 소득과 의료비 지출의 연관성은 더욱 강하게 측정될 것이다. 

충분 통계량 (sufficient statistics) 

충분 통계량 (sufficient statistics) 란 계수를 설명할 수 있는 충분한 정보를 가지고 있는 통계량입니다. 예를 들어, 평균이 theta 이고 분산이 1인 정규분포에서 X1,X2,X3...Xn을 관찰하였고, X의 표본 평균값을 구해서 3이 나왔습니다. 그러면 X의 평균값 이외에 X1,X2,X3,...Xn 각각에 대한 데이터를 추가적으로 봄으로써 theta의 추정에 도움이 되는 다른 정보를 얻을 수 있을까요? 만약 추가적인 정보가 없다면, X의 표본 평균은 theta 를 설명할 수 있는 충분한 통계량이라고 볼 수 있습니다. 좀더 포멀하게는 아래와 같이 충분 통계량을 정의할 수 있습니다.

$$ f(X|Y;\theta) \ne g(\theta) $$

위 식을 만족하면 Y는 theta 에 대한 충분통계량입니다. 위 식의 의미는 Y가 주어졌을 때, X의 pdf 가 theta 에 의존적이지 않다라는 의미입니다. 즉, Y가 theta 에 대한 모든 정보를 갖고 있다는 의미입니다.  

충분 통계량의 예시

예를 들어, 베르누이를 2회 시행해서 X1,X2 를 구했다고 합니다. X는 0 또는 1의 값을 갖습니다. 

이 때, X1 + X2 는 베르누이의 모수 p 에 대한 충분 통계량일까요? 

만약 X1+X2가 0이라면, X1,X2가 가질 수 있는 경우의 수는 (0,0) 입니다.

만약 X1+X2가 1이라면, (0,1), (1,0) 입니다.

만약 X1+X2가 2라면 (1,1) 입니다. 

즉, X1+X2 가 주어졌을 때, (X1,X2) 의 분포 (pmf) 는 theta 와 독립적입니다. 따라서 X1+X2 는 충분 통계량이라고 볼 수 있습니다. 이 의미는 만약에 X1+X2 가 1이라고 한다면, X1이 1이냐, X2가 1이냐는 적어도 계수를 추정함에 있어서 중요하지 않다는 것입니다. 위 예시에서는 베르누이 2회 시행에서 예시를 들었으나, N번의 베르누이 시행에서도 X1+...+Xn은 충분 통계량입니다. 만약 100번의 베르누이 시행에서 53번 성공했다라고 한다면, 몇 번째 시행에서 성공했는지는 계수 추정에 있어서 추가적인 정보를 주지는 않습니다. 

Maximum likelihood estimation

MLE는 관측한 데이터를 가장 잘 설명할 수 있는 계수(parameter) 를 찾는 것이 목적이라고 할 수 있습니다. 데이터 분석을 하다보면 데이터를 관측하고, 이를 통해 모델의 계수를 추정해야할 때가 있습니다. 우선 likelihood 를 정의해보면, 특정 계수에서 데이터를 관찰할 가능성을 의미합니다. 

예를 들어, 계수가 0.4인 베르누이분포를 100번 시행해서 45번의 성공이 발생했다고 합시다. likelihood 란 계수의 함수이며, 이 경우에서는 아래와 같이 나타낼 수 있습니다. 

$$ L(\theta) = {100 \choose 45}\theta^{45}(1-\theta)^{55} $$

이 그래프는 계수가 0.1일 때 데이터를 관찰할 확률, 0.2 일때 관찰할 확률 등등.. 그 값을 모두 구해 계수와 likelihood 의 관계를 나타내는 그래프라고 볼 수 있습니다. 이 때, likelihood 는 계수가 0.45 일 때 최댓값을 갖습니다. 따라서 MLE 는 0.45 입니다. 이 경우 MLE 와 실제 계수에 차이가 있습니다. 하지만, 관측하는 데이터의 갯수가 증가할 수록, MLE 는 실제 계수인 0.4로 수렴해갈 것입니다 (하지만 모든 분포에서 실제 계수에 수렴하게 되는 것은 아닙니다).  

Probability 와 Likelihood 는 다른 개념입니다. Likelihood 는 데이터가 주어져있고, 이를 통해 계수를 추정하기 위한 것이라고 볼 수 있습니다. Probability 란, 어떤 모델의 계수가 주어졌을 때, 특정 outcome 을 예측하기 위한 것입니다. 

Maximum likelihood 를 만드는 지점을 찾기

그러면 likelihood function 의 값을 최대화하는 계수값을 어떻게 구할 수 있을까요? 일반적으로는 likelihood function 을 미분해서 값이 0인 지점을 찾으면 됩니다. 또는 계산상의 이점을 위해 log likelihood 를 미분해서 값이 0이 되는 지점을 찾습니다. likelihood 를 미분해서 0이 되는 지점이나, log likelihood 를 미분해서 0이 되는 지점이나 값은 같기 때문입니다. 

예를 들어, 앞서 살펴본 베르누이 분포에서 MLE 를 통해 계수 추정을 해봅시다. 베르누이 분포의 pmf 는 아래와 같이 정의 됩니다. 

$$ \theta^{x}(1-\theta)^{1-x} $$

베르누이 분포의 likelihood 는 아래와 같이 쓸 수 있습니다.

$$ L(\theta) =  \theta^{\sum_{i=1}^{n}{x_i}} 1-\theta^{n- \sum_{i=1}^{n}{x_i} }  $$  

그러면 log likelihood는 아래와 같습니다. 

$$ l(\theta) = \log \theta \sum_{i=1}^{n}{x_i} + \log (1-\theta)(n-\sum_{i=1}^{n}{x_i} ) $$

log likelihood 를 미분하면 아래와 같습니다. (이를 score 라고 부르기도 합니다.)

$$ l'(\theta) = \frac{1}{\theta} \sum_{i=1}^{n}{x_i} - \frac{1}{1-\theta}( n-\sum_{i=1}^{n}{x_i} )$$

위 score를 0으로 만드는 theta (계수) 값은 표본 평균입니다. (이와 같이 많은 경우, MLE 는 표본 평균인 경우가 많습니다.)

$$ \hat{\theta} = \bar{X} $$ 

또한 이 경우 표본 평균의 기댓값이 모평균이기 때문에 unbiased estimation 이라고 볼 수 있습니다.

MLE 는 항상 unbiased 는 아니다. 

하지만 MLE 의 경우 모든 경우에 unbiased 인 것은 아닙니다. 예를 들어, uniform distribution(0,theta) 의 계수를 MLE 로 추정값은 관측한 데이터중 최댓값이 됩니다(관련된 자세한 설명은 생략하겠습니다). 그리고 이값의 기댓값은 아래와 같이 n에 의존하는 값이 되기 때문에 biased estimation 이라고 할 수 있습니다.(참고링크) 

MLE 의 분산과 Efficiency

그럼에도 불구하고 MLE 는 관측 데이터가 많으면 많을 수록 가능한 unbiased estimator 집합들이 가질 수 있는 분산의 최솟값으로 수렴한다는 큰 장점을 가지고 있습니다. 즉, MLE 는 이론적으로 나올 수 있는 가장 작은 분산을 가진다는 의미이며, 이를 "Asymptotically efficient 하다" 라고 부르기도 합니다. 

이 때, 이론적으로 나올 수 있는 분산의 최솟값을 Rao-Cramer Lower Bound 라고 부릅니다. 어떤 unbiased estimator Y가 있을 때, Y의 분산의 최솟값은 아래과 같습니다. 

$$ var(Y) \ge \frac{1}{nI_1{\theta}} $$

여기서 I는 fisher information 이고 아래와 같습니다. I_1 은 데이터 하나로 구한 fisher information 입니다. 

$$ I(\theta) = E(-l''(\theta))   $$

$$  I(\theta) = nI_1(\theta)$$

베르누이 분포에서 fisher information 과 rao-cramer lower bound 를 구해봅시다. 데이터 1개에서 log likelihood 를 구하면 fisher information 을 구하고 이를 통해 rao-cramer lower bound 를 구할 수 있습니다.

$$l(\theta) = xlog\theta + (1-x)log(1-\theta)$$

$$l'(\theta) = \frac{x}{\theta} - \frac{1-x}{1-\theta}$$

$$l''(\theta) = \frac{x}{\theta^2}-\frac{1-x}{(1-\theta)^2}$$

$$I_1(\theta) = \frac{1}{\theta(1-\theta)}$$

rao-cramer lower bound 는 아래와 같습니다. 

$$ var(\hat{\theta}) \ge \frac{\theta(1-\theta)}{n} $$ 

실제 계수 추정값의 분산을 구해봅시다. 

$$ var(\hat{\theta}) = var(\bar{X}) = var(\frac{  \sum_{i=1}^{n}{x_i}  }{n}) = \frac{1}{n^2} \sum_{i=1}^{n}{var(x_i)} = \frac{\theta(1-\theta)}{n}  $$

이 값은 rao-cramer lower bound 의 값과 같습니다. 따라서 베르누이 분포에서 MLE 는 최소 분산을 가지며, efficient estimator 라고 할 수 있습니다. 또한 unbiased 이기 때문에, minimum variance unbiased  estimator (MVUE) 라고 부르기도 합니다. 

Autoregressive Processes

Autoregressive process 란 history 가 현재 값에 직접적인 영향을 주는 time series 를 말한다. 식으로 표현하면 아래와 같다. 

현재시점 (t) 의 값은 과거 시점들 (t-1 ~ t-p) 의 값을 가중치를 두고 합한 값에 error term (Z) 을 더한 값이다. 

$$ X_t = Z_t + \phi_1(X_{t-1})  + \phi_2(X_{t-2}) ... + \phi_p(X_{t-p}) $$

Example

p=2 이고, 가중치가 0.7, 0.2 인 autoregressive process 는 아래와 같다.

$$ X_t = Z_t + 0.7X_{t-1} + 0.2X_{t-2} $$

이를 R 코드로 구현하면 아래와 같다. 선차트를 통해 보면 현재 값이 과거 값과 높은 상관성이 있다는것을 확인할 수 있다. correlogram 을 통해 가까운 시간에 측정된 값이 현재값과 더 높은 상관성이 있다는것을 확인할 수 있다. (가중치가 0.7, 0.2 이므로)

set.seed(2017)
X.ts <- arima.sim(list(ar=c(0.7,0.2)), n=1000)
par(mfrow=c(2,1))
plot(X.ts, main="AR(2) Time series, phi1=0.7, phi2=0.2")


X.acf <- acf(X.ts, main="Autocorrelation of AR(2) Time series")
X.acf

Moving average process 와의 관계 

Autoregressive process 는 moving average process 의 무한 수열로 나타낼 수 있다. 

차수 (p) 가 1인 AR 을 생각해보자.

아래와 같이 식을 쓸 수 있다. (Z는 평균이 0, 분산이 sigma^2 을 따른다고 가정하고, phi 를 theta 로 치환하자.) 

$$ X_t = Z_t + \phi X_{t-1} = Z_t + \phi Z_{t-1} + \phi^2 X_{t-2} ... = Z_t + \theta_1 Z_{t-1} + \theta_2 X_{t-2} ... $$

AR (1) 의 통계량

위와 같이 식을 써서 moving average 처럼 표현하면 AR process 에서의 X(t) 의 기댓값과 분산을 쉽게 구해볼 수 있다. 

$$ E(X_t) = 0 $$

$$ V(X_t) = \sigma^2 \sum^{\infty}_{i=0}\theta_i^2  $$ 

time series 가 stationarity 를 만족하기 위해서는 분산이 t 에따라 바뀌지 않고, 일정해야한다. 그렇기 때문에 phi 의 절댓값이 1보다 작은것은 stationarity 에 대한 필요 조건이라고 할 수 있다. 

AR(1) 의 auto covariance function 

MA process 의 acf 와 유사한 형태라는것을 확인할 수 있다. (link)

$$ \gamma(k) = \sigma^2 \sum^{\infty}_{i=0} \theta_i \theta_{i+k} $$ 

AR(1) 의 auto correlation coefficient 

$$ \rho(k) = \frac{\sum^{\infty}_{i=0} \theta_i \theta_{i+k}}{\sum^{\infty}_{i=0} \theta_i \theta_i} $$ 

theta 를 phi 로 치환하여 AR(1) 의 auto covariance 와 auto covariance coefficient 를 구해보자. (무한등비수열 공식 사용하여 정리) 

$$ \gamma(k) = \sigma^2 \frac{\phi^k}{1-\phi^2} $$

$$ \rho(k) = \phi^k $$

AR Process 의 Stationarity 를 확인하는 방법

Example: AR(1) process

$$ X_t = Z_t + \phi X_{t-1} $$ 

위 AR(1) process 에서 Z 를 제외한 나머지 텀들을 한쪽으로 옮겨서 아래와 같은 식을 만들 수 있다.

$$ \phi(B) = 1-\phi B $$ 

이 때, 우변을 0으로 만드는 B 의 해를 찾는다. 해는 B = 1/phi 이다. B 의 해가 단위원 (unit circle) 바깥에 있는 것이  stationarity 를 만족하기 위한 조건이 된다. 따라서 AR(1) 모델에서는 phi 의 절댓값이 1 미만이어야 stationarity 를 만족한다. 

Example: AR(2) process 

$$ X_t = \frac{1}{3} X_{t-1} + \frac{1}{2} X_{t-2} + Z_t $$

$$ \phi(B) = 1-\frac{1}{3} B-\frac{1}{2} B^2 $$

위 식에서 B 의 해를 찾으면 (-2+sqrt(76))/6, (-2-sqrt(76))/6 이 된다. 위 두 값이 모두 단위원 바깥에 있기 때문에 위 AR(2) process 는 stationarity 를 만족한다. 

Week stationarity

이전 포스팅에서 개념을 직관적으로 소개하였다. 이번엔 좀 더 포멀한 정의와 함께 예시를 통해 설명해보려고한다.

아래 조건을 만족할 때 weekly stationary 라고 부른다. 

1) 시간에 따른 평균이 같다. 

$$ \mu(t) = \mu  $$ 

2) Auto covariance function 이 time spacing 에만 의존한다. (t2=t1+tau 라고 생각하면 이해가 쉽다.)

$$ \gamma(t_1, t_2) = \gamma(t_2-t_1) = \gamma(\tau) $$

: 이는 시간에 따른 분산이 같다는 조건을 포함하는 조건이다. 

Examples

실제 시계열 데이터의 예시를 통해 stationarity 에 대해 더 이해해보려고한다. 

1) White noise 는 stationarity 를 만족한다.

White noise model

$$ X_t \sim N(0, \sigma) $$

White noise model은 시간에 따른 평균이 같다. 

$$ \mu = 0 $$

White noise model은 Auto covariance function 이 time spacing 에만 의존한다.

$$ \gamma(t_1, t_2) = 0, \ when \ t_1 = t_2 $$

$$ \gamma(t_1, t_2) = \sigma^2, \ when \ t_1 \neq t_2 $$

2) Random walk 는 stationarity 를 만족하지 않는다.

Random walk model

아래 식에서 X(t) 를 random walk model 이라고 한다. 

$$ Z_t \sim iid(\mu, \sigma^2) $$

$$ X_t = X_{t-1} + Z_t = \sum^{t}_{i=1}Z_i $$

따라서 Random walk model 의 시간에 따른 평균은 t*mu 이고, 분산은 t*sigma^2 이다. 만약 Z의 평균이 0이라고 가정하더라도 분산이 시간에 따라 점점 커진다는 것을 알 수 있다. 따라서 Random walk model 은 stationarity 를 만족하지 않는다. 

3) Moving average model 는 stationarity 를 만족한다. 

moving average model 

$$ X_t = Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} ... + \theta_q Z_{t-q} $$

$$ Z_t \sim Normal(\mu, \sigma) $$

moving average 의 parameter q 와 가중치 theta 를 고정해놓고 계산을 하면, 평균과 분산은 t 와는 관계 없이 고정된다는 것을 알 수 있다. 따라서 moving average model 은 stationarity 를 만족한다. 

추가적으로 Moving average model 의 auto covariance function 을 구해보자. 

moving average model 은 stationarity 를 만족하기 때문에 auto covariance function 은 time spacing 에만 의존한다. 또한 이전 포스팅에서 time spacing 이 최대 q 인 경우에만 자기상관성이 존재한다는 것을 correlogram 을 통해 확인할 수 있었다. moving average model 의 노이즈의 평균이 0일 때를 가정하고 covaraicne 를 구해보자.

$$ Z_t \sim Normal(0, \sigma) $$

$$ Cov(X_t, X_{t+k}) = E(X_t X_{t+k}) - E(X_t)E(X_{t+k}) =  E(X_t X_{t+k}) $$

위 기댓값을 정리하면 아래와 같은 식이 된다.

$$ Cov(X_t, X_{t+k}) = E(X_t X_{t+k}) = \sigma^2 \sum_{i=0}^{q-k} \beta_i \beta_{i+k}, \ when \ k \leq q $$

$$ Cov(X_t, X_{t+k}), \ when \ k \gt q $$ 

참고) 위 식의 정리에는 아래 기댓값과 분산의 성질을 이용하면 된다.

$$ V(X) = E(X^2) - E(X)^2 $$

$$ E(XY) = E(X)E(Y), \ when \ X, \  Y \ is \ independent $$ 

Moving average process

특정 시점 t에서의 주가를 X_t 라고하자. 또한 특정 시점 t 에서의 회사의 공지 Z_t (noise) 가 주가에 영향을 미친다고 하자. 그런데 과거 시점 (t-1, t-2...) 에 회사의 공지도 주가에 영향을 미친다. 이런 경우에 X_t 를 다음과 같이 모델링할 수 있다.

$$ X_t = Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} ... + \theta_q Z_{t-q} $$

$$ Z_t \sim Normal(\mu, \sigma) $$

이 때, q는 어떤 시점의 noise 까지 현재값에 영향을 미치는지를 의미하며, MA(2) 는 이와 같이 정의된다. 

MA(2) Process

$$ X_t = Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} $$

$$ Z_t \sim Normal(\mu, \sigma) $$

MA(2) process 의 예를 들면 아래와 같다. 

$$ X_t = Z_t + 0.7 Z_{t-1} + 0.2 Z_{t-2} $$

MA(2) simulation

위에 예시로 든 MA(2) process 를 R 을 통해 simulation 해보자.

# noise 생성
noise <-rnorm(10000)

ma_2 = NULL

# ma(2) 생성을 위한 loop
for (i in 3:10000) {
  ma_2[i] = noise[i] + 0.7*noise[i-1]+0.2*noise[i-2]
}

# shift 
moving_average_process <- ma_2[3:10000]
moving_average_process <- ts(moving_average_process)

par(mfrow=c(2,1))

plot(moving_average_process, main = "A moving average process of order 2", ylab = "")
acf(moving_average_process, main = "Correlogram of ma (2)", ylab = "")

correlogram 을 보면 time step 이 0,1,2 인 경우에만 상관성이 있는 것을 확인할 수 있다. 우선, time step 이 0 인 경우는 항상 auto correlation coef 1이다. 또한 현재값에는 최대 2 time step 전의 noise 까지 반영이 되기 때문에, 최대 2 time step 의 값과 상관성이 있다는 것을 확인할 수 있따. 

Random walk model

아래와 같이 정의되는 X_t 를 random walk 이라고 한다. X_t 는 이전 time step 에서의 값 X_t-1 에 noise Z가 더해진 값이다. random sampling 과 다른점은 현재값이 이전값에 더해진다는것이다. 이는 랜덤하게 어떤 한 방향으로 걷는것과 비슷하다. 매번 시작점에서 한발짝 걷는 것이 아니라 한발짝 걸어서 도착한 곳에서 다시 한발짝을 간다. 

$$ X_t = X_{t-1} + Z_t $$

$$ Z_t \sim Normal(\mu, \sigma) $$

이러한 random walk 모델에서 X_t 는 이전 time step 에서의 값 X_t-1 과 매우 큰 연관성을 갖는다. 따라서 non-stationary time series 데이터이다. 

Random walk model simulation in R

R 로 random walk 모델을 만들어보자. 아래는 1000개의 random walk 데이터를 생성하는 예제이다. 시계열 그래프를 그려보면, 이 데이터는 non-stationary time series 데이터라는 것을 확인할 수 있다. 구간을 나눠서보면 트렌드를 보이기 때문이다. 

x <- NULL
x[1] <- 0
for(i in 2:1000){
  x[i] <- x[i-1]+rnorm(1)
}

random_walk <- ts(x)
plot(random_walk, main="A random walk", ylab="", xlab=" Days", col="black")

random walk 데이터에서 correlogram 을 그려보자. 인접한 time step 에서 auto correlation coefficient 가 큰 패턴을 보이기 때문에 non-stationary time series 라는 것을 다시 확인할 수 있다. 

acf_result <- acf(random_walk)

random walk 모델에서 noise Z는 stationary time series 라고 볼 수 있다. 

$$ Z_t \sim Normal(\mu, \sigma) $$

noise 가 stationary time series 라는 것을 데이터로 실제로 확인해보자. 

random_walk_diff <- diff(random_walk)
plot(random_walk_diff, main="A random walk diff", ylab="", xlab=" Days", col="black")

참고자료

https://people.duke.edu/~rnau/411rand.htm

Auto correlation coefficient 

앞선 포스팅에서 auto covariance coefficient 에 대해 설명하였다. auto covariance coefficient 은 time series 데이터에서의 각각의 time point 간 연관성을 의미하는데, stationary time series 에서는 k 라고하는 parameter 에 의해 달라진다. auto covariance coefficient 의 추정값 c_k 는 아래와 같이 계산된다. 

$$ c_k = \frac{\sum^{N-k}_{t=1}(x_t - \bar{x})(x_{t+k}-\bar{x})} {N-k} $$

이번에는 auto correlation coefficient 에 대해 정리해보려고 한다. auto correlation coefficient 도 auto covariance coefficient 와 마찬가지로 time series 데이터에서 time step 별 값의 연관성을 의미하는데 범위를 -1~1로 조정한 것으로 이해할 수 있다. 마치 공분산과 상관계수의 관계와 같다.

auto correlation coefficient 의 계산식은 아래와 같다. 

$$ -1 \leq \rho_k = \frac{\gamma_k}{\gamma_0} \leq 1 $$

이 때, rho_k 의 추정값 r_k 는 아래와 같다. c0 는 time step (=lag) 이 0일 때의 auto covariance coefficient 로 이는 분산과 같다.  

$$ \rho_k \sim r_k = \frac{c_k}{c_0} $$

$$ c_0 = \frac{\sum^{N}_{t=1}(x_t-\bar{x})^2}{N}  $$

Correlogram

rho_k 의 추정값을 k 에 따라 구한 뒤, 이를 시각화해서 표현한 것을 correlogram 이라고 한다. 이 때, rho_0 은 항상 1이다. (자기 자신과의 연관성이기 때문이다.) 

$$ r_0 = \frac{c_0}{c_0} , r_1 = \frac{c_1}{c_0} ... $$ 

R 에서 correlogram 은 acf 함수를 통해 쉽게 그려볼 수 있다.

예제 1)

아래 R 코드는 100개의 표준정규분포를 따르는 데이터를 만든 후, correlogram 을 그리는 코드이다. 파란선은 연관성이 유의한지에 대한 임계치를 의미한다. 유의한 데이터 포인트가 하나 밖에 없고, lag 에 따른 패턴이 보이지 않으므로, 전체적으로 시계열 데이터가 자기상관성이 없다고 결론 내릴 수 있다.  

purely_random_process <- ts(rnorm(100))
print(purely_random_process)
plot(purely_random_process)

auto_correlation_coef_by_lags <- acf(purely_random_process)
print(auto_correlation_coef_by_lags)

예제 2)

실제 데이터를 correlogram 을 그려보자. 다음은 모 어플리케이션의 월간 활성 이용자수 (MAU, monthly active user) 추이이다. 이 서비스는 점점 성장하는 추이를 보여주고 있다. 시계열 데이터의 관점에서는 시간에 따른 평균의 변화 (trend) 를 보이는 non-stationary time series 이다. 

위 데이터에서 correlogram 을 그리면 아래와 같이 나타난다. lag 에 따른 auto correlation coef 의 패턴이 보이며 (점점 감소), 인접한 데이터 포인트에서는 유의한 상관성을 보이고 있는 것을 확인할 수 있다. 

Time series data (시계열 데이터) 

어떤 종류의 데이터이든 상관 없으며, 그저 시간에 따라 수집된 데이터를 시계열 데이터 (timeseries data) 라고 한다. 

예를 들어, 일별 코로나 확진자수는 1일이라고 하는 time step 으로 수집된 시계열 데이터의 한 종류이다. 

Week stationary time series 

week stationary time series 란 다음의 조건을 만족한다. 

1) 시간에 따른 평균 (mean) 에 변화가 없다. 

2) 시간에 따른 분산 (variance) 의 변화가 없다. 

3) 주기적인 등락 (flucation) 이 없다.

이러한 조건을 만족하긴 위해서는 time series 의 한 섹션 (A 섹션) 고른 후, 다른 섹션 (B 섹션) 을 골랐을 때, A, B 섹션이 비슷하면 된다. 

Stochastic process

random variable 의 collection - X1,X2,X3 .. 가 있다고 하자. 이들이 각각 다른 모수를 가진 분포를 따를 때, 이를 stochastic process 라고 한다. stochastic process 의 반대개념은 deterministic process 이다. deterministic process 는 모든 step (t) 에 대해서 예측 가능하다. 예를 들어, 어떤 함수에 대한 미분함수는 특정 X 에서의 Y 값을 정확하게 알 수 있다. 이와 반대로 stochastic process 는 매 step 이 random 이기 때문에 어떤 확률 분포에서 왔다는 것만을 알 수 있을 뿐, 값을 정확하게 예측할 수 없다.

$$ X_t \sim distribution(\mu_t, \sigma_t) $$

예를 들어, 다음과 같은 시계열 데이터가 있다고 해보자. 

$$ X_1 = 30, X_2 = 29, X_3 = 57 ... $$ 

시계열 데이터를 바라보는 한 가지 관점은 stochastic process 의 실현 (realization) 으로 보는것이다. 매 timestep 별로 어떤 확률 변수가 정해지고 우리는 그 확률변수에서 나온 하나의 샘플값을 관찰하는 것이다. 

Autocovariance function

stationary time series 라는 가정을 하자. 두 가지 timestep s,t 에서의 covariance 를 정의해볼 수 있다. (확률 변수이기 때문)

$$ \gamma(s,t) = Cov(X_s, X_t) $$

$$ \gamma(s,s) = Var(X_s)   $$  

또한 아래처럼 covariance function 을 정의할 수 있는데, 이 함수는 stationary time series 라는 가정 하에 t 에 따라서는 값이 바뀌지 않으며, k가 결정하는 함수가 된다. (아래 식에서 c 는 추정값이다.) 이러한 time step (k) 에 따른 공분산의 식을  autocovariance function  이라고 한다. 

$$ \gamma_k = \gamma(t, t+k) \sim c_k $$ 

즉, stationary time series 에서는 Cov(X1, X2) 나 Cov(X10,X11) 이나 기댓값은 같다고 할 수 있다. 그 이유는 데이터에서 두 가지 섹션을 선택했을 때, 그 모습이 똑같다고 기대하는것이 stationary time series 이기 때문이다.  

또한 gamma(t_k, t) 는 autocovariance coefficient 라고 하며, stochastic process 에서의 실제 autocovariance 값이다. 데이터를 통해 구한 c_k 를 통해 autocovariance coefficient 를 추정한다. 

Autocovariance coefficient

그러면 Autocovariance coefficient 의 추정값은 어떻게 구할까? timestep 을 k 라고 할 때, 추정값은 아래와 같다. 

$$ c_k = \frac{\sum^{N-k}_{t=1}(x_t - \bar{x})(x_{t+k}-\bar{x})} {N} $$

이는 time series 데이터에서 k time step 만큼 차이나는 점들의 묶음 (x_t, x) 을 확률 변수의 관찰값으로 놓고 추정한 covariance 와 같다.

참고로 확률 변수 X,Y 의 covariance 의 추정값은 아래와 같이 구할 수 있다. 

$$ s_{xy} =  \frac{\sum^{N}_{t=1}(x_t - \bar{x})(y_t-\bar{y})} {N-1} $$

R 에서는 acf 함수를 통해 auto covariance coefficient 추정값을 계산할 수 있다. 

purely_random_process <- ts(rnorm(100))
print(purely_random_process)
plot(purely_random_process)

auto_covariance_coef_by_lags <- acf(purely_random_process, type = "covariance")
print(auto_covariance_coef_by_lags)