먼저 R Studio를 실행하여 과제를 저장할 파일부터 만듭니다.
새 스크립트를 열고
과제를 제출할 이름으로 저장합니다.
저장버튼 누르고 -> 폴더경로 지정해주고(영문으로만)-> 과제명.R로 정해서 save 버튼을 누르면 됩니다.
작성은 Script 창에서 합니다. Console 창에서 하면 저장할 수 없습니다.
이제 문제를 풀어봅니다.
1. Holt-Winters 계절지수평활법을 이용한 예측
(1) quantmod package의 getSymbols 함수를 이용하여 FRED에서 다음의 data를 download 받는다.
quantmod 패키지의 getsymbols 함수를 이용한다고 하였으므로 라이브러리 명령어를 사용하여
퀀트모드를 불러옵니다.
> library("quantmod")
FRED에서 TLOFCON의 데이터를 가져옵니다.
> getSymbols("TLOFCON",src="FRED")
(2)이렇게 받은 data 중 2013년 1월부터 2017년 12월까지의 data를 tcs라는 변수에 저장한다.
> tcs<-TLOFCON["2013/2017"]
(3) tcs를 ts함수를 써서 ts class 변수로 변환하여 tcs.ts라는 변수에 저장한다.
> tcs.ts<-ts(tcs,frequency=12,start=c(2013,1))
(4) tcs.ts에 ets함수를 사용하여 Holt-Winters 계절지수평활법을 적용하고 결과를 tcs.fit에 저장한다. ets함수 적용시 다음에 유의한다.
Holt-Winters 계절지수평활법 적용
수평성(Level) 평활지수는 0.5로 고정
모형을 추정하는데 있어 MSE(mean square error)를 최소화하는 방법으로 실행
모형선택에 있어 BIC 선정기준 채택
이건 ets 함수를 사용할 줄 알아야 합니다. 정석은 ?ets 실행해서 보는 것이고..
강의 노트의 4장에 자세히 기록되어 있으나 문제풀이에 필요한 부분위주로 살펴보자면
자, ets 함수에 들어가는 다양한 변수 중에 여기서 우리에게 주어진 조건을 빨간색으로 칠했습니다.
그리고 강의자료 예제로 사용한 ets 함수를 보면
ets(
+ 넣을 변수/여기서는 tcs.ts
+ 모델/Holt-Winters 계절지수평활법 적용이니 위의 설명대로 "AAA",
+ 수평성 평활지수 0.5는 뭔말이죠. 평활상수를 말씀하시는 거 같은데...
평활지수는 처음 들었어요. 평활상수 alpha=0.5
+ 그리고 bic를 입력해야 하니 ic=c(bic")
+)
여기까지를 코드로 입력하면
> ?ets
> tcs.fit<-ets(tcs.ts,model="AAA",alpha=0.5,ic=("bic"))
(5) summary함수를 이용하여 tcs.fit 내용을 출력한다.
> summary(tcs.fit)
실행화면은 다음과 같이 도출됩니다.
(6)원시계열자료인 tcs.ts는 검은색 실선으로 ets함수로 적합한 값들은 붉은색 점선 및 o표식으로 나타내어 다음의 그래프를 그린다.
강의자료의 4장 평활법을 보면, 유사한 문제가 나옵니다.
이를 참조해서 코드를 입력하면(문제에서 요구하지 않았으므로 legend는 넣을 필요 없습니다)
> plot(tcs.fit$x)
> lines(tcs.fit$fitted,col="red",type="o",lty=2)
(7) forecast함수를 사용하여 2018년 1월부터 4월까지 예측하고 그 결과를 tcs_f에 저장한다. 그리고, tcs_f 출력하고 plot(tcs_f) 출력한다.
먼저 forecast함수를 사용하여 2018년 1월부터 4월까지 예측합니다.
역시 강의자료에서 다음의 내용을 참조할 수 있습니다.
> forecast(tcs.fit,h=4)
> tcs_f<-forecast(tcs.fit,h=4)
> tcs_f
> plot(tcs_f)
(8) 모형의 잔차(tcs.fit의 잔차)의 ACF와 PACF를 조사한다. 아울러 해당 잔차의 Ljung-Box test를 lag=12로 검정한다.
역시 강의자료4장에서 다음의 자료를 참조할 수 있습니다.
> plot(tcs.fit$residuals)
> acf(tcs.fit$residuals)
> pacf(tcs.fit$residuals)
> Box.test(tcs.fit$residuals,lag=12,type="Ljung-Box")
유의수준 5%에서 p-value가 너무 낮으므로 귀무가설을 기각합니다. 즉 whitenoise가 아닙니다.
(9) 처음 시작시에 받은 “Total Construction Spending: Office”의 2018년 1월부터 4월의 실제값과 tcs.fit에서 예측한 동일기간 값들을 다음과 같은 그래프로 나타낸다. 검은색 실선은 실제값이고 붉은색 점선 및 o표식은 예측값이다. y축 값들의 범위를 5300에서 6200으로 지정한다. 이는 plot함수의 ylim 옵션을 이용하여 지정할 수 있다.
먼저 2018년 1월부터 4월의 실제값
> ts(TLOFCON["2018-01/2018-04"],frequency=12,start=c(2018,1))
tcs.fit에서 예측한 동일기간 값
> tcs_f$mean
실행하면
이 둘을 한 화면에 그래프로 나타내려면
> plot(ts(TLOFCON["2018-01/2018-04"],frequency=12,start=c(2018,1)),ylim=range(5300:6200))
> lines(tcs_f$mean,ylim=range(5300:6200),col="red",type="b",lty=2)
(참고로 type="o" 인 경우는 o표식 안에 선이 간섭되고, type="b"인 경우는 o표식에 선이 간섭되지 않음
1번의 Code 정리
1.
(1)
library("quantmod")
getSymbols("TLOFCON",src="FRED")
(2)
tcs<-TLOFCON["2013/2017"]
(3)
tcs.ts<-ts(tcs,frequency=12,start=c(2013,1))
(4)
library("forecast")
?ets
tcs.fit<-ets(tcs.ts,model="AAA",alpha=0.5,ic=("bic"))
(5)
summary(tcs.fit)
(6)
plot(tcs.fit$x)
lines(tcs.fit$fitted,col="red",type="o",lty=2)
(7)
forecast(tcs.fit,h=4)
tcs_f<-forecast(tcs.fit,h=4)
tcs_f
plot(tcs_f)
(8)
plot(tcs.fit$residuals)
acf(tcs.fit$residuals)
pacf(tcs.fit$residuals)
Box.test(tcs.fit$residuals,lag=12,type="Ljung-Box")
(9)
ts(TLOFCON["2018-01/2018-04"],frequency=12,start=c(2018,1))
tcs_f$mean
plot(ts(TLOFCON["2018-01/2018-04"],frequency=12,start=c(2018,1)),ylim=range(5300:6200))
lines(tcs_f$mean,ylim=range(5300:6200),col="red",type="b",lty=2)
2.요소분해법을 이용한 DJIA 분석
(1) quantmod package의 getSymbols 함수를 이용하여 Yahoo! Finance에서 Dow Jones Industrial Average 지수를 내려받는다. 내려받은 OHLC data 중 2008년부터 2017년까지의 data를 dji_s에 저장한다.
야후 파이낸스에서 검색하면
이렇게 Dow Jones Industrial Average의 symbol이 '^DJI'임을 알 수 있습니다.
퀀트모드 불러와서 겟심볼로 ^DJI를 내려받고, 문제대로 데이터를 추출해서 dji_s에 저장합니다.
> library("quantmod")
> getSymbols("^DJI")
^DJI는 DJI로 받아졌음을 알 수 있습니다. 내용파악을 위해 head로 살펴보고 dji_s로 저장합니다.
> head(DJI)
> dji_s<-DJI["2008/2017"]
(2) dji_s를 월별 data로 변환하여 dji_mth에 저장하고 dji_mth 중 수정종가만을 dji_p에 저장한다.
그리고 head(dji_p)를 출력한다.
먼저 월별 data로 변환하여 저장하고
> dji_mth<-to.monthly(dji_s)
dji_mth의 데이터에서 수정종가의 명칭을 보기 위해 head(dji_mth)를 실행합니다.
> head(dji_mth)
dji_s.Adjusted임을 확인하고
> dji_p<-dji_mth$dji_s.Adjusted
> head(dji_p)
가 도출됩니다.
(3) ts함수를 이용하여 dji_p를 ts class 변수로 변환하여 dji.ts에 저장하고 dji.ts를 출력한다.
> dji.ts<-ts(dji_p,frequency=12,start=c(2008,1))
> dji.ts
(4) decompose함수를 이용하여 가법모형의 요소분해법을 dji.ts에 적용하고 그 결과를 dji.d에 저장한다. 또한, summary함수를 이용하여 dji.d의 summary를 출력한다.
> dji.d<-decompose(dji_ts,type="additive")
> summary(dji.d)
(5) par(mfrow=c(2,1))을 이용하여 dji.d의 원시계열과 추세를 한 그래프에 그린다.
(4)번문제 summary에서 dji.d가 x, seasonal, trend, random, figure로 구분됨을 알 수 있습니다.
> par(mfrow=c(2,1))
> plot(dji.d$x)
> plot(dji.d$trend)
실행하면
의 그래프가 그려집니다.
(6) par(mfrow=c(2,1))을 이용하여 dji.d의 계절변동과 잔차를 한 그래프에 그린다.
> par(mfrow=c(2,1))
> plot(dji.d$seasonal)
> plot(dji.d$random)
실행하면
의 그래프가 그려집니다.
2번의 Code 정리
2.
(1)
library("quantmod")
getSymbols("^DJI")
head(DJI)
dji_s<-DJI["2008/2017"]
(2)
dji_mth<-to.monthly(dji_s)
head(dji_mth)
dji_p<-dji_mth$dji_s.Adjusted
head(dji_p)
(3)
dji.ts<-ts(dji_p,frequency=12,start=c(2008,1))
dji.ts
(4)
dji.d<-decompose(dji_ts,type="additive")
summary(dji.d)
(5)
head(dji.d)
par(mfrow=c(2,1))
plot(dji.d$x)
plot(dji.d$trend)
(6)
par(mfrow=c(2,1))
plot(dji.d$seasonal)
plot(dji.d$random)
3. ARMA 모델링
(1) quantmod package의 getSymbols 함수를 이용하여 FRED에서 다음의 data를 download 받는다.
> getSymbols("EXUSEU",src="FRED")
EXUSEU가 어떤 데이터인지 봐야겠죠. 그냥 EXESEU를 실행합니다.
1999-01-01부터 2018-06-01까지의 환율 데이터(U.S. Dollars to One Euro)를 알 수 있습니다.
(2) 이렇게 받은 data 중 1999년 1월부터 2017년 12월까지의 data를 ts함수를 이용하여 ts class 변수로 변환한 다음 usdeur라는 변수에 저장한다. 또한, 2018년 1월부터 2018년 6월까지 data는 역시 ts class 변수로 변환한 후 usdeur.os라는 변수에 저장한다.
> usdeur<-ts(EXUSEU["1999-01/2017-12"],frequency=12,start=c(1999,1))
> usdeur.os<-ts(EXUSEU["2018-01/2018-6"],frequency=12,start=c(2018,1))
(3) 앞서 저장한 usdeur에 대하여 auto.arima함수를이용하여 ARMA모형 추정을 하여 usdeur.fit에 결과를 저장한다. 다음에 유의한다.
stationary 옵션을 이용하여 usdeur data를 차분하지 않고 (즉, d=0) ARMA 모형을 추정한다.
seasonal 옵션을 이용하여 ARMA모형 추정에 있어 계절변동을 고려하지 않는다.
auto.arima 함수를 이용하는 법을 먼저 봅니다.
강의노트 6장을 보면 다음과 같이 forecast package에 오토아리마 함수가 있음을 알 수 있습니다.
먼저 forecast 패키지를 불러온 후 오토 아리마를 사용합니다. (필요시 ?auto.arima로 검색해도 되겠죠)
그런데,전 이해가 안됩니다. 오토아리마가 뭐고 알마를 왜 알마함수가 아닌 오토아리마함수로..검색해보니 어렴풋하게 이해를 도와주는 자료(http://www.dodomira.com/2016/04/21/arima_in_r/)를 보자면
ARIMA는 AR 모형과 MA 모형을 결합한 것(I) 인데, ARIMA 모델에서는 세가지 모형을 위한 세가지 파라미터 (p, d, q)가 필요합니다. 파라미터를 구하기 위해서는 AR(p)모형의 p차수 MA(q)의 q차수 그리고 트랜드를 제거하여 안정시계열로 만들기 위한 I(d)의 차분 차수 d를 결정하기 위해 KPSS test4 , ACF, PACF를 그려보아야 합니다.
하지만 R에서는 auto.arima를 사용해 자동으로 p, d, q를 결정할 수가 있습니다.
정리하자면 'auto.arima 함수를 사용하면 p,d,q가 자동으로 결정되니 arma에서의 p,q 또한 자동으로 결정되어 편리하게 사용할 수 있다' 입니다.(다른 참고 교재에서는 auto.arima : ARIMA 모형의 차수결정 프로그램으로 정리되어있습니다)
문제의 조건에 해당하는 빨간줄을 보시면 자료는 usdeur, 차분차수d=0, 계절변동을 고려하지 않으니 D=0
> auto.arima(usdeur,d=0,D=0)
> usdeur.fit<-auto.arima(usdeur,d=0,D=0)
(4) summary함수를 이용하여 usdeur.fit의 summary를 출력한다
> summary(usdeur.fit)
(5) 2017년의 usdeur와 추정된 ARMA 모형의 적합값을 다음과 같은 그래프로 그린다.
(window함수를 이용하면 2017년 해당값들만 추려낼 수 있다. 4장 강의슬라이드 참조)
의 내용으로 윈도우 함수를 써봅니다.
> plot(window(usdeur.fit$x,start=c(2017,1),end=c(2017,12)))
> lines(window(usdeur.fit$fitted,start=c(2017,1),end=c(2017,12)),type="b",col="red",lty=2)
(6) forecast함수를 이용하여 usdeur의 2018년 1월-6월 예측을 하고 그 결과를 usdeur_f에 저장한
다. 아울러 usdeur_f를 출력한다
usdeur은 이미 2017까지만 저장되어 있으므로, usdeur자료 끝 후의 6개월간을 추정합니다
> usdeur_f<-forecast(usdeur,h=6)
> usdeur_f
(7) usdeur.os (2018년 1월-6월 실제값)과 (6)에서 예측한 값을 다음과 같은 그래프로 나타낸다. 검은색 실선이 실제값이고 붉은색 점선과 o표식이 예측값을 나타낸다.
> plot(usdeur.os)
> lines(usdeur_f$mean,col="red",type="b",lty=2)
실행하면 다음과 같이 나옵니다.
3번의 Code 정리
3.
(1)
getSymbols("EXUSEU",src="FRED")
EXUSEU
(2)
usdeur<-ts(EXUSEU["1999-01/2017-12"],frequency=12,start=c(1999,1))
usdeur.os<-ts(EXUSEU["2018-01/2018-6"],frequency=12,start=c(2018,1))
(3)
plot(usdeur)
library("forecast")
?auto.arima
auto.arima(usdeur,d=0,D=0)
usdeur.fit<-auto.arima(usdeur,d=0,D=0)
(4)
summary(usdeur.fit)
(5)
plot(window(usdeur.fit$x,start=c(2017,1),end=c(2017,12)))
lines(window(usdeur.fit$fitted,start=c(2017,1),end=c(2017,12)),col="red", type="b", lty=2)
(6)
usdeur_f<-forecast(usdeur,h=6)
usdeur_f
(7)
plot(usdeur.os)
lines(usdeur_f$mean,col="red",type="b",lty=2)
여기까지 오시느라, 수고하셨습니다.
궁금한 점은 언제든 물어봐주세요. 제가 대답할 수 있는 한 답변드리겠습니다.