- 참고 : Statistics and Data Analysis for Financial Engineering
1. Hstarts 데이터셋
Package "Ecdat"에서 제공하는 Hstarts 데이터셋은 1960년 1분기부터 2001년 4분기까지 캐나다의 분기별 도시 주택 착공 건수의 로그값이다.
hs hssa
[1,] 7.98933 NA
[2,] 8.83961 NA
[3,] 8.94841 NA
[4,] 8.98907 NA
[5,] 8.39087 NA
[6,] 9.09204 NAx <- ts(Hstarts[,1],
start = 1960, # 시계열 시작 날짜
frequency = 4) # 분기별 시계열로 1년에 4번 관측
par(mfrow=c(3,1)) # 3개의 그래프를 한 화면에 출력
# 시계열 그림
plot(x,
xlab = "year", ylab = "log(starts)",
type = "l", # 선으로 표시
main="(a)")
# 자기상관계수 그림
acf(x,
main = "(b)",
xlab = "lag")
# 분기별 상자그림
quart = rep(1:4, 42) # 분기
boxplot(x ~ quart,
xlab = "quarter", ylab = "log(starts)",
main = "(c)")
Result! (a) 그래프를 통해 해당 시계열은 계절성이 강하다는 것을 알 수 있다.
(b) 그래프를 통해 4의 배수인 시차에서는 자기상관이 크며, 다른 시차에서는 천천히 0으로 감소한다는 것을 알 수 있다.
(c) 그래프를 통해 해당 시계열의 관측값은 1분기에 가장 낮으며, 2분기에는 높다는 것을 알 수 있다.
par(mfrow=c(3,2)) # 한 행에 2개의 그래프를 출력 -> 총 3개 행으로 6개의 그래프가 출력됨
plot(diff(x), # 1번 차분한 시계열 그림
xlab = "year",
type = "l",
main = "(a) nonseasonal differencing")
acf(diff(x), # 1번 차분한 시계열의 자기상관계수 그림
main = "(b) nonseasonal differencing",
xlab = "lag")
plot(diff(x, lag = 4), # 주기 4만큼 계절 차분한 시계열 그림
type = "l",
xlab = "year",
main = "(c) seasonal differencing")
acf(diff(x, lag = 4), # 계절 차분한 시계열의 자기상관계수 그림
main = "(d) seasonal differencing",
xlab = "lag")
plot(diff(diff(x, diff = 1), lag = 4), # 1번 차분한 후 주기 4만큼 계절 차분한 시계열 그림
type = "l",
xlab = "year",
main = "(e) seasonal & nonseasonal differencing")
acf(diff(diff(x, diff = 1), lag = 4), # 1번 차분한 후 주기 4만큼 계절 차분한 시계열의 자기상관계수 그림
main = "(f) seasonal & nonseasonal differencing",
xlab = "lag")
Caution! 함수 diff()를 이용하여 차분을 수행할 수 있다. \(k\)번 차분 \(\Delta^kY_{t}\)을 수행하고 싶으면 옵션 diff = k를 입력하면 되고, 주기 \(s\)만큼 계절차분 \(\Delta_sY_{t}\)을 수행하고 싶으면 옵션 lag = s를 입력하면 된다.
Result! (a)와 (b) 그래프를 통해 1번 차분한 시계열은 계절성을 보이며, 높은 자기상관을 가진다는 것을 알 수 있다.
(c)와 (d) 그래프를 통해 원 시계열에 계절차분을 수행함으로써 계절성이 제거된 것을 알 수 있다.
(e)와 (f) 그래프를 통해 원 시계열에 계절차분과 차분을 동시에 수행함으로써 계절성이 제거되고 시간의 흐름에 따라 평균이 변하지 않는다는 것을 알 수 있다.
1-1. SARIMA 모형
pacman::p_load("forecast")
# 1. ARIMA(1,1,1)(1,1,1)[4] 모형 구축
pacman::p_load("forecast")
fit1 <- arima(x,
c(1,1,1), # ARIMA(p,d,q)
seasonal = list(order = c(1,1,1), # 계절성(P,D,Q)
period = 4)) # 주기 s
fit1
Call:
arima(x = x, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 4))
Coefficients:
ar1 ma1 sar1 sma1
0.6883 -0.8819 -0.1550 -0.7659
s.e. 0.2386 0.1875 0.0974 0.0812
sigma^2 estimated as 0.02573: log likelihood = 64.13, aic = -118.27Caution! 함수 arima()의 옵션 seasonal을 이용하여 SARIMA 모형을 구축할 수 있다.
Result! 모형 추정 결과에 의하면, 구축된 ARIMA(1,1,1)(1,1,1)[4] 모형은 \((1-0.6883B)(1+0.1550B^4)(1-B)(1-B^4)Y_t=(1-0.8819B)(1-0.7659B^4)\epsilon_t\)이다.
# 잔차를 이용한 모형 진단
Box.test(residuals(fit1), lag = 10,
type = "Ljung-Box",
fitdf = 4) # 추정한 phi, Phi, theta, Theta 개수
Box-Ljung test
data: residuals(fit1)
X-squared = 5.8219, df = 6, p-value = 0.4434Result! 귀무가설 \(H_0 : \rho(1)=\rho(2)=\cdots=\rho(10)=0\)에 대한 검정 결과에 따르면, \(p\)값이 0.4434이므로 유의수준 0.05에서 \(p\)값이 0.05보다 크기 때문에 귀무가설을 기각하지 못한다. 즉, 잔차에 대해 시차 10까지의 자기상관계수 \(\rho(1), \rho(2), \cdots, \rho(10)\) 중 유의한 자기상관계수가 적어도 1개 존재한다는 증거가 부족하며, 해당 시계열에 대해 ARIMA(1,1,1)(1,1,1)[4] 모형을 가정하는 것이 적절하다.
# 2. ARIMA(1,1,1)(0,1,1)[4] 모형 구축
fit2 <- arima(x,
c(1,1,1), # ARIMA(p,d,q)
seasonal = list(order = c(0,1,1), # 계절성(P,D,Q)
period = 4)) # 주기 s
fit2
Call:
arima(x = x, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))
Coefficients:
ar1 ma1 sma1
0.6748 -0.8901 -0.8220
s.e. 0.1419 0.1046 0.0512
sigma^2 estimated as 0.02609: log likelihood = 62.91, aic = -117.82Result! 모형 추정 결과에 의하면, 구축된 ARIMA(1,1,1)(0,1,1)[4] 모형은 \((1-0.6748B)(1-B)(1-B^4)Y_t=(1-0.8901B)(1-0.8220B^4)\epsilon_t\)이다.
# 잔차를 이용한 모형 진단
Box.test(residuals(fit2), lag = 10,
type = "Ljung-Box",
fitdf = 3) # 추정한 phi, Phi, theta, Theta 개수
Box-Ljung test
data: residuals(fit2)
X-squared = 8.45, df = 7, p-value = 0.2946Result! 귀무가설 \(H_0 : \rho(1)=\rho(2)=\cdots=\rho(10)=0\)에 대한 검정 결과에 따르면, \(p\)값이 0.2946이므로 유의수준 0.05에서 \(p\)값이 0.05보다 크기 때문에 귀무가설을 기각하지 못한다. 즉, 잔차에 대해 시차 10까지의 자기상관계수 \(\rho(1), \rho(2), \cdots, \rho(10)\) 중 유의한 자기상관계수가 적어도 1개 존재한다는 증거가 부족하며, 해당 시계열에 대해 ARIMA(1,1,1)(0,1,1)[4] 모형을 가정하는 것이 적절하다.
1-2. 예측
1-1. SARIMA 모형에서 구축한 ARIMA(1,1,1)(0,1,1)[4] 모형과 ARIMA(1,1,1)(0,1,1)[4] 모형은 적절한 것으로 판단되나, 두 번째 모형이 첫 번째 모형보다 모수 개수가 적기 때문에 두 번째 모형이 더 합리적으로 보인다(모형의 간결성). 두 번째 모형을 이용한 예측은 다음 코드를 통해 수행된다.
# 예측
pred <- forecast(fit2, h = 16) # 미래 16시점까지 예측
pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2002 Q1 9.036139 8.829131 9.243147 8.719547 9.352730
2002 Q2 9.581109 9.317976 9.844241 9.178683 9.983535
2002 Q3 9.475528 9.180979 9.770076 9.025055 9.926000
2002 Q4 9.395209 9.080061 9.710357 8.913232 9.877186
2003 Q1 9.004073 8.661140 9.347006 8.479602 9.528544
2003 Q2 9.563662 9.200980 9.926344 9.008988 10.118336
2003 Q3 9.467946 9.089962 9.845930 8.889869 10.046022
2003 Q4 9.394284 9.003641 9.784926 8.796847 9.991720
2004 Q1 9.007639 8.595867 9.419411 8.377888 9.637390
2004 Q2 9.570259 9.141715 9.998803 8.914858 10.225660
2004 Q3 9.476588 9.033912 9.919264 8.799574 10.153602
2004 Q4 9.404306 8.949180 9.859432 8.708250 10.100361
2005 Q1 9.018593 8.542702 9.494484 8.290780 9.746405
2005 Q2 9.581841 9.088848 10.074834 8.827873 10.335809
2005 Q3 9.488594 8.980786 9.996402 8.711969 10.265219
2005 Q4 9.416598 8.895471 9.937725 8.619603 10.213593
2. AirPassengers 데이터셋
R에서 제공하는 AirPassengers 데이터셋은 1949년 1월부터 1960년 12월까지 월별 비행기 탑승객 수가 기록되어져 있다.
# 데이터 불러오기
data("AirPassengers")
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# 시계열 그림
plot(AirPassengers)
Result! 해당 시계열은 추세와 계절성을 모두 가지고 있는 비정상시계열임을 알 수 있다.
# 자기상관계수 그림
acf(AirPassengers)
Result! 자기상관계수 ACF는 천천히 감소하므로 비정상시계열임을 알 수 있다.
par(mfrow = c(2,1)) # 2개의 그래프를 한 화면에 출력
plot(diff(diff(AirPassengers, diff = 1), lag = 12), # 1번 차분한 후 주기 12만큼 계절 차분한 시계열 그림
type = "l")
# 자기상관계수 그림
acf(diff(diff(AirPassengers, diff = 1), lag = 12))
Result! 계절차분과 차분을 수행한 AirPassengers 시계열은 정상시계열로 보인다.
2-1. SARIMA 모형
# 1. ARIMA(1,1,1)(1,1,1)[12] 모형 구축
fit <- arima(AirPassengers,
c(1,1,1), # ARIMA(p,d,q)
seasonal = list(order = c(1,1,1), # 계절성(P,D,Q)
period = 12)) # 주기 s
fit
Call:
arima(x = AirPassengers, order = c(1, 1, 1), seasonal = list(order = c(1, 1,
1), period = 12))
Coefficients:
ar1 ma1 sar1 sma1
-0.1387 -0.2027 -0.9228 0.8329
s.e. 0.5860 0.6123 0.2386 0.3519
sigma^2 estimated as 130.8: log likelihood = -506.15, aic = 1022.3Result! 모형 추정 결과에 의하면, 구축된 ARIMA(1,1,1)(0,1,1)[12] 모형은 \((1+0.1387B)(1+0.9228B^{12})(1-B)(1-B^{12})Y_t=(1-0.2027B)(1+0.8329B^{12})\epsilon_t\)이다.
# 잔차를 이용한 모형 진단
Box.test(residuals(fit), lag = 10,
type = "Ljung-Box",
fitdf = 4) # 추정한 phi, Phi, theta, Theta 개수
Box-Ljung test
data: residuals(fit)
X-squared = 10.476, df = 6, p-value = 0.106Result! 귀무가설 \(H_0 : \rho(1)=\rho(2)=\cdots=\rho(10)=0\)에 대한 검정 결과에 따르면, \(p\)값이 0.106이므로 유의수준 0.05에서 \(p\)값이 0.05보다 크기 때문에 귀무가설을 기각하지 못한다. 즉, 잔차에 대해 시차 10까지의 자기상관계수 \(\rho(1), \rho(2), \cdots, \rho(10)\) 중 유의한 자기상관계수가 적어도 1개 존재한다는 증거가 부족하며, 해당 시계열에 대해 ARIMA(1,1,1)(0,1,1)[12] 모형을 가정하는 것이 적절하다.
2-2. 예측
# 예측
pred <- forecast(fit, h = 48) # 미래 48시점까지 예측
pred Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 1961 449.2223 434.5545 463.8901 426.7899 471.6547
Feb 1961 424.3072 406.7452 441.8691 397.4485 451.1658
Mar 1961 459.0023 438.6160 479.3886 427.8241 490.1805
Apr 1961 497.7220 474.9011 520.5430 462.8204 532.6237
May 1961 509.7462 484.7209 534.7715 471.4733 548.0190
Jun 1961 568.1989 541.1491 595.2488 526.8297 609.5681
Jul 1961 655.7364 626.8032 684.6696 611.4868 699.9859
Aug 1961 641.1725 610.4712 671.8737 594.2190 688.1259
Sep 1961 546.3741 514.0013 578.7470 496.8642 595.8841
Oct 1961 496.7500 462.7878 530.7123 444.8092 548.6908
Nov 1961 427.6964 392.2158 463.1769 373.4335 481.9592
Dec 1961 471.2915 434.3550 508.2281 414.8019 527.7811
Jan 1962 484.9196 441.0846 528.7546 417.8797 551.9595
Feb 1962 458.8334 411.0374 506.6295 385.7356 531.9312
Mar 1962 487.3741 435.6837 539.0644 408.3205 566.4277
Apr 1962 529.1175 473.8375 584.3975 444.5741 613.6609
May 1962 540.1970 481.5427 598.8514 450.4929 629.9012
Jun 1962 602.8458 541.0014 664.6902 508.2629 697.4286
Jul 1962 689.8873 625.0095 754.7651 590.6652 789.1094
Aug 1962 673.9982 606.2226 741.7738 570.3443 777.6521
Sep 1962 576.2455 505.6910 646.8000 468.3417 684.1493
Oct 1962 529.0428 455.8148 602.2708 417.0502 641.0354
Nov 1962 458.1931 382.3859 534.0004 342.2559 574.1303
Dec 1962 500.3163 422.0147 578.6180 380.5643 620.0683
Jan 1963 517.2610 431.8276 602.6945 386.6018 647.9203
Feb 1963 492.2555 401.9707 582.5402 354.1768 630.3341
Mar 1963 526.4753 431.3724 621.5782 381.0279 671.9226
Apr 1963 565.4285 465.7688 665.0883 413.0121 717.8450
May 1963 577.3797 473.3587 681.4007 418.2933 736.4661
Jun 1963 636.1565 527.9505 744.3626 470.6696 801.6435
Jul 1963 723.6557 611.4204 835.8910 552.0066 895.3047
Aug 1963 708.9894 592.8647 825.1142 531.3919 886.5869
Sep 1963 613.9629 494.0748 733.8511 430.6099 797.3160
Oct 1963 564.5258 440.9889 688.0627 375.5924 753.4592
Nov 1963 495.3334 368.2525 622.4144 300.9799 689.6870
Dec 1963 538.8149 408.2861 669.3437 339.1883 738.4415
Jan 1964 552.6991 415.1965 690.2018 342.4069 762.9914
Feb 1964 526.6965 383.8656 669.5273 308.2555 745.1374
Mar 1964 555.6757 407.5409 703.8105 329.1231 782.2284
Apr 1964 597.2036 443.9709 750.4363 362.8545 831.5528
May 1964 608.3505 450.1811 766.5198 366.4513 850.2496
Jun 1964 670.7002 507.7441 833.6563 421.4803 919.9201
Jul 1964 757.7771 590.1707 925.3834 501.4453 1014.1088
Aug 1964 741.9824 569.8515 914.1133 478.7309 1005.2339
Sep 1964 644.4402 467.9006 820.9798 374.4462 914.4342
Oct 1964 597.0650 416.2242 777.9058 320.4928 873.6371
Nov 1964 526.3433 341.3012 711.3854 243.3459 809.3408
Dec 1964 568.5714 379.4213 757.7215 279.2913 857.8515