[시계열 자료 분석] R에서 AirPassengers 데이터 선형계절추세모형 적합시키기

R에서 AirPassengers 데이터 선형계절추세모형(linear and seasonal trend model)에 적합시키는 방법

AirPassengers 데이터는 ts 타입의 데이터로, 1949년부터 1960년까지 매 월 한 포인트의 데이터를 가지고 있다.

> AirPassengers

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

1949 112 118 132 129 121 135 148 148 136 119 104 118

1950 115 126 141 135 125 149 170 170 158 133 114 140

1951 145 150 178 163 172 178 199 199 184 162 146 166

1952 171 180 193 181 183 218 230 242 209 191 172 194

1953 196 196 236 235 229 243 264 272 237 211 180 201

1954 204 188 235 227 234 264 302 293 259 229 203 229

1955 242 233 267 269 270 315 364 347 312 274 237 278

1956 284 277 317 313 318 374 413 405 355 306 271 306

1957 315 301 356 348 355 422 465 467 404 347 305 336

1958 340 318 362 348 363 435 491 505 404 359 310 337

1959 360 342 406 396 420 472 548 559 463 407 362 405

1960 417 391 419 461 472 535 622 606 508 461 390 432

---

1. 시계열 데이터 정상화

시간이 지날 수록 분산이 점차 커지는 것을 볼 수 있다. 시계열 데이터를 정상화시켜주기 위해 log변환을 한다.

* 분산의 상수화

변수 Zt에 대해 분산 \(\sigma^2_{t} = Var(Z_{t}) \)와 평균 \(\mu_{t} = E(Z_{t})\) 사이에는 보통 함수관계가 존재한다. 이 함수관계를 \(f\)라고 하면

$$\sigma^2_{t} = f(\mu_{t})$$

가 성립한다.

즉, 이 함수관계 때문에 분산이 일정하지 않고 변화하게 된다. 따라서 \(Z_{t}\)를 변환시켜 분산이 상수가 되도록 해야한다.

가장 많이 사용되는 함수 \(f\)의 형태로는 다음과 같은 두 가지가 있다.

① \(\sigma^2_{t} = f(\mu_{t}) c\mu^2_{t}\)

② \(\sigma^2_{t} = f(\mu_{t}) c\mu_{t}\)

①의 경우 \(Z_{t}\)를 자연로그 변환인 \(lnZ_{t}\) 해주고,

②의 경우 제곱근변화인 \(\sqrt(Z_{t})\)로 안정화시켜준다.

AirPassengers 데이터는 ①의 경우이기 때문에 \(lnZ_{t}\)로 변환해준다.

v <- c()
m <- c()
 
for (i in 1:length(AirPassengers)) {
  v <- c(v,var(AirPassengers[1:i]))
  m <- c(m,mean(AirPassengers[1:i]))
}
 
par(mfrow = c(1,2))
plot(v,m, main="분산과 평균의 관계")
slope <- (sort(m)[108]-sort(m)[36])/(sort(v)[108]-sort(v)[36])
inter <- sort(m)[36]-slope*sort(v)[36]
abline(a = inter, b = slope, col="red")
 
plot(v,m^2, main = "분산과 평균^2의 관계")
slope <- (sort(m^2)[108]-sort(m^2)[36])/(sort(v)[108]-sort(v)[36])
inter <- sort(m^2)[36]-slope*sort(v)[36]
abline(a = inter, b = slope, col="red")

[그림 1] \(\sigma^2_{t}와 \mu_{t}의 관계\) [그림 2] \(\sigma^2_{t}와 \mu^2_{t}의 관계\)

> AP <- log(AirPassengers)

---

2. 시계열 데이터의 체계적 성분 분해

시계열 데이터에서 추세성분과 계절성분을 따로 추정하기 위해 분해법으로 체계적 성분을 분리한다.

> AP.decomp <- decompose(AP, type = "multiplicative")

> plot(AP.decomp)

* decompose()

R에서 시계열 데이터는 stl() 또는 decompose() 함수로 trend와 seasonal 성분을 분해할 수 있다.

>> stl() in R 바로가기 : https://leedakyeong.tistory.com/entry/R-stl-%EC%9D%B4%EB%9E%80-stl-parameter-stl-swindow

decompose() 에서 type을 "multiplicative"로 설정해주면 \(Z_{t} = T_{t}*S_{t}*I_{t}\) 식에 부합하도록 분해할 수 있다.

> round(AP.decomp$trend*AP.decomp$seasonal*AP.decomp$random,2)

      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec

1949   NA   NA   NA   NA   NA   NA 5.00 5.00 4.91 4.78 4.64 4.77

1950 4.74 4.84 4.95 4.91 4.83 5.00 5.14 5.14 5.06 4.89 4.74 4.94

1951 4.98 5.01 5.18 5.09 5.15 5.18 5.29 5.29 5.21 5.09 4.98 5.11

1952 5.14 5.19 5.26 5.20 5.21 5.38 5.44 5.49 5.34 5.25 5.15 5.27

1953 5.28 5.28 5.46 5.46 5.43 5.49 5.58 5.61 5.47 5.35 5.19 5.30

1954 5.32 5.24 5.46 5.42 5.46 5.58 5.71 5.68 5.56 5.43 5.31 5.43

1955 5.49 5.45 5.59 5.59 5.60 5.75 5.90 5.85 5.74 5.61 5.47 5.63

1956 5.65 5.62 5.76 5.75 5.76 5.92 6.02 6.00 5.87 5.72 5.60 5.72

1957 5.75 5.71 5.87 5.85 5.87 6.05 6.14 6.15 6.00 5.85 5.72 5.82

1958 5.83 5.76 5.89 5.85 5.89 6.08 6.20 6.22 6.00 5.88 5.74 5.82

1959 5.89 5.83 6.01 5.98 6.04 6.16 6.31 6.33 6.14 6.01 5.89 6.00

1960 6.03 5.97 6.04 6.13 6.16 6.28   NA   NA   NA   NA   NA   NA

> round(AP,2)

      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec

1949 4.72 4.77 4.88 4.86 4.80 4.91 5.00 5.00 4.91 4.78 4.64 4.77

1950 4.74 4.84 4.95 4.91 4.83 5.00 5.14 5.14 5.06 4.89 4.74 4.94

1951 4.98 5.01 5.18 5.09 5.15 5.18 5.29 5.29 5.21 5.09 4.98 5.11

1952 5.14 5.19 5.26 5.20 5.21 5.38 5.44 5.49 5.34 5.25 5.15 5.27

1953 5.28 5.28 5.46 5.46 5.43 5.49 5.58 5.61 5.47 5.35 5.19 5.30

1954 5.32 5.24 5.46 5.42 5.46 5.58 5.71 5.68 5.56 5.43 5.31 5.43

1955 5.49 5.45 5.59 5.59 5.60 5.75 5.90 5.85 5.74 5.61 5.47 5.63

1956 5.65 5.62 5.76 5.75 5.76 5.92 6.02 6.00 5.87 5.72 5.60 5.72

1957 5.75 5.71 5.87 5.85 5.87 6.05 6.14 6.15 6.00 5.85 5.72 5.82

1958 5.83 5.76 5.89 5.85 5.89 6.08 6.20 6.22 6.00 5.88 5.74 5.82

1959 5.89 5.83 6.01 5.98 6.04 6.16 6.31 6.33 6.14 6.01 5.89 6.00

1960 6.03 5.97 6.04 6.13 6.16 6.28 6.43 6.41 6.23 6.13 5.97 6.07

---

3. 추세성분 추정 :: 회귀모형 적합

분해법으로 분해된 성분 중 추세성분을 추정한다. x축을 시간(t), y축을 분해된 추세성분으로 plot을 그리면 다음과 같다.

> t <- 1:144

> plot(t, AP.decomp$trend, type = "l")

이를 x는 시간(t), y는 trend로 회귀모형에 적합했을 때 잔차에 대한 plot은 다음과 같다.

> t <- 1:144

> trend.fit1 <- lm(AP.decomp$trend~t)

> par(mfrow=c(2,2))

> plot(trend.fit1)

첫 번째 plot에서 곡선형태가 나타난다. 이를 해결하기 위해 회귀모형에 다항식을 추가한다.

> t <- 1:144

> t2 <- t^2

> trend.fit2 <- lm(AP.decomp$trend~t+t2)

> par(mfrow=c(2,2))

> plot(trend.fit2)

왼쪽 아래 plot을 보면 오차항들의 자기상관관계가 있어보인다. 이를 해결하기 위해 AR 모형에 적합시켜야 하나, 이는 나중에 따로 다루도록 하겠다.

* 잔차분석

오차(잔차)는 다음과 같은 가정을 만족해야 한다.

① 오차 \(\varepsilon_{t}\)의 평균은 0이고 분산은 \(\sigma^2_{\varepsilon}\)이다.

② 오차 \(\varepsilon_{t}\)들은 상관관계가 없다. 즉, 모든 \(t_{1}\)과 \(t_{2}(t_{1} \neq t_{2})\)에 대하여 

$$ Cov(\varepsilon_{t1},\varepsilon_{t2}) = E(\varepsilon_{t1},\varepsilon_{t2}) = 0 $$

③ 오차 \(\varepsilon_{t}\)들은 정규분포를 따른다.

실제값(○)과 모형에 적합한 값(-)으로 그림을 그리면 다음과 같다.

왼쪽 그림은 종속변수에 t하나만, 오른쪽은 종속변수에 t와 t^2을 추가한 모형이다.

첫 번째 모델에 대한 \(R^2\)는 0.9879이고, 두 번째 모델에 대한 \(R^2\)는 0.9957이다. 

---

4. 계절성분 추정 :: 계절추세모형 적합

분해법으로 분해된 계절성분은 다음과 같다.

> plot(1:144,AP.decomp$seasonal, type = "l")

> abline(v=seq(0,144,by=12),col="red")

데이터에서도 알 수 있듯이 12주기로 반복한다.

* 계절추세모형

계절추세모형은 주기가 \(s\)인 경우 다음 식으로 적합시킬 수 있다.

$$ Z_{t} = \beta_{0} + \sum_{i=1}^{m} A_{i}sin(\frac{2\pi i}{s}t + \phi_{i}) + \varepsilon_{t} $$

\(A_{i}\) : 진동폭

\(\phi_{i}\) : 편각

\(w_{i} = \frac{2\pi i}{s}\) : \(2\pi\) 단위시간 동안 관측되는 \(i\)번째 주기항의 빈도

\(\frac{s}{i}\) : \(i\)번째 주기항의 주기

$$ Z_{t} = \beta_{0} + \sum_{i=1}^{m} A_{i}sin(\frac{2\pi i}{s}t + \phi_{i}) + \varepsilon_{t} = \beta_{0} + \sum_{i=1}^{m} [\beta_{1i}sin(w_{i}t)+\beta_{2i}cos(w_{i}t)] + \varepsilon_{t} $$

이다.

ⅰ) m=1 일 때

 

> x11 <- sin(2*3.14*t/12)

> x21 <- cos(2*3.14*t/12)

> seasonal.fit1 <- lm(AP.decomp$seasonal~x11+x21)

> seasonal.pred1 <- predict(seasonal.fit1, newdata = data.frame(x11=x11, x21 = x21))

> plot(1:144,AP.decomp$seasonal, type="l")

> lines(seasonal.pred1, col="red", type="l")

실제값의 주기항이 3개처럼 보이기 때문에(한 주기동안 뾰족한 부분이 3개) m=3 일 때로 해보면 다음과 같다.

ⅱ) m=3 일 때

> x11 <- sin(2*3.14*t/12)

> x21 <- cos(2*3.14*t/12)

> x12 <- sin(2*3.14*t/12*2)

> x22 <- cos(2*3.14*t/12*2)

> x13 <- sin(2*3.14*t/12*3)

> x23 <- cos(2*3.14*t/12*3)

> seasonal.fit2 <- lm(AP.decomp$seasonal~x11+x21+x12+x22+x13+x23)

> seasonal.pred2 <- predict(seasonal.fit3, newdata = data.frame(x11=x11, x21 = x21, x12=x12, x22 = x22, x13=x13, x23 = x23))

> plot(1:144,AP.decomp$seasonal, type="l")

> lines(seasonal.pred2, col="red", type="l")

m=1 일 때와, m=3일 때 각각 \(R^2\)는 0.713과 0.944이다.

---

5. 선형계절추세모형

위에서 적합시킨 두 모델을 합쳐 최종적으로 하나의 모델을 만든다.

시간(t)에 대한 2차항을 포함시킨 trend를 추정하는 선형모형과 m=3일 때로 seasonal을 추정하는 계절추세모형으로 적합시킨 결과는 다음과 같다.

> plot(1:144,AP, type="l")

> lines(trend.pred2*seasonal.pred2, col="red", type="l")

실제값(-) 추정값(-)

처음에 분산을 안정화시켜주기 위해 \(ln(Z_{t})\) 변환을 했기 때문에 이를 다시 되돌려주면 다음과 같다.

> plot(1:144,AirPassengers, type="l")

> lines(exp(trend.pred2*seasonal.pred2), col="red", type="l")

---

6. 잔차분석

최종적으로 만들어진 모델의 적합성을 확인하기 위해 잔차분석을 시행한다.

> e <- AirPassengers-exp(trend.pred2*seasonal.pred2)

> par(mfrow=c(1,3))

> plot(e)

> plot(exp(trend.pred2*seasonal.pred2), e)

> qqnorm(e)

> qqline(e,col="red")

---

* trend 와 seasonal 을 따로 추정하지 않고 한 번에 추정하려면 다음과 같이 하면 된다.

t <- 1:144
t2 <- t^2
 
x11 <- sin(2*3.14*t/12)
x21 <- cos(2*3.14*t/12)
x12 <- sin(2*3.14*t/12*2)
x22 <- cos(2*3.14*t/12*2)
x13 <- sin(2*3.14*t/12*3)
x23 <- cos(2*3.14*t/12*3)
 
fit <- lm(AP ~ t+t2+x11+x21+x12+x22+x13+x23)

잔차 plot

> par(mfrow=c(2,2))

> plot(fit)

잔차와 적합값 사이에 어떠한 특징 없이 random한 형상을 보이며, 등분산성과 정규성 모두 만족한다.

실제값 vs 적합값

> plot(t,AirPassengers, type="l")

> lines(exp(pred), col="red")