728x90
반응형
Forecasting: Principles and Practice 2nd Edition을 공부한 내용을 기록, 정리하고 있습니다.
5.10 연습문제
1번
a. elecdaily 테이블을 그래프로 나타내고, 온도를 설명변수로 사용해서 수요값에 대한 회귀모델 찾기.
일반적으로 30도까지는 온도가 올라갈수록 추울때보다 난방에 에너지를 덜 사용하게 되어 음의 관계가 나온다.
elecdaily 데이터의 모델의 계수, 절편값 구하기
> tslm(Demand ~ Temperature, data=daily20)
Call:
tslm(formula = Demand ~ Temperature, data = daily20)
Coefficients:
(Intercept) Temperature
39.212 6.757
autoplot(daily20, facets=TRUE)
daily20 %>%
as.data.frame %>%
ggplot(aes(x=Temperature, y=Demand)) +
geom_point() +
geom_smooth(method="lm", se=FALSE)
fit <- tslm(Demand ~ Temperature, data=daily20)
checkresiduals(fit)
forecast(fit, newdata=data.frame(Temperature=c(15,35)))
b. 잔차 그래프. 특별히 outlier는 보이지 않는다?
> df <- as.data.frame(daily20)
> df[, "residuals"] <- as.numeric(residuals(fit.daily20))
> ggplot(df, aes(x=Temperature, y=Demand)) + geom_point() + xlab("temperature") + ylab("demand")
> ggplot(df, aes(x=Temperature, y=residuals)) + geom_point() + xlab("temperature") + ylab("Residuals")
c. 최고기온이 15도일때, 다음날의 전기 소비량을 예측하는 모델을 사용하고, 최고기온이 35도 일때 얻은 예측값과 비교. 신뢰할 만 한가?
기온이 15도 → 39.2117 + 6.7572*15 = [1] 140.5697
기온이 35도 → 39.2117 + 6.7572*35 = [1] 275.7137
신뢰할만한가... 아마도...? 근데 15도는 데이터가 없어서 잘 모르겠음..?
> fit.daily20 <- tslm(Demand ~ Temperature, data=daily20)
> summary(fit.daily20)
Call:
tslm(formula = Demand ~ Temperature, data = daily20)
Residuals:
Min 1Q Median 3Q Max
-46.060 -7.117 -1.437 17.484 27.102
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.2117 17.9915 2.179 0.0428 *
Temperature 6.7572 0.6114 11.052 1.88e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 22 on 18 degrees of freedom
Multiple R-squared: 0.8716, Adjusted R-squared: 0.8644
F-statistic: 122.1 on 1 and 18 DF, p-value: 1.876e-09d. 예측 구간 만들어보기.
> forecast(fit, newdata=data.frame(Temperature=c(15,35)))
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
4.285714 140.5701 108.6810 172.4591 90.21166 190.9285
4.428571 275.7146 245.2278 306.2014 227.57056 323.8586전제조건 확인. 히스토그램이 좀 애매하긴 하다. 그외 다른 부분은 자기상관관계가 없음으로 뜸
> checkresiduals(fit)
Breusch-Godfrey test for serial correlation of order up to 5
data: Residuals from Linear regression model
LM test = 3.8079, df = 5, p-value = 0.5774
e. elecdaily에 있는 모든 데이터를 가지고 수요 vs 기온 그래프를 그려봅시다. 이것이 여러분의 모델과 관련하여 어떤 의미가 있습니까?
적은 데이터보다 전체 데이터를 사용했을 때, 선형회귀를 적용하기 어려운 데이터가 될 수도 있다..?
2번
a. plot 그리고 설명하기.
시간이 지날수록 감소하는 트렌드, 중간중간 전쟁기간동안에는 없는 데이터. 계절성은 없음.
b. regression line 그리기
> time_mens400 <- time(mens400)
> fit.mens400 <- tslm(mens400 ~ time_mens400, data=mens400)
> autoplot(mens400) + geom_abline(slope = fit.mens400$coefficients[2], intercept = fit.mens400$coefficients[1], color='red')
c. residuals against year
ACF에서는 자기상관관계 있음. 브로이쉬-갓프레이에서는 유의미하지 않음. → 자기상관이 막 엄청 크지는 않다. fitted line에는 잘 맞는거 같다.
> checkresiduals(fit.mens400)
Breusch-Godfrey test for serial correlation of order up to 6
data: Residuals from Linear regression model
LM test = 3.6082, df = 6, p-value = 0.7295
d. 2020년의 prediction interval 구하기
> forecast(fit.mens400, newdata = data.frame(time_mens400=2020))
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2020 42.04231 40.44975 43.63487 39.55286 44.53176
3번
easter(ausbeer) 쳐보고 해석하기. 해당 코드를 돌렸을 때, 이스터가 들어가는 분기는 1, 아닌 분기는 0으로 표시된다. 이스터는 주로 봄에 일어나므로 1분기나 2분기에 포함된다. 이스터가 1분기에 시작해서 2분기에 끝나는 식이면 1956년처럼 쪼개져서 기록된다.
## 너무 길어서 잘랐음
> easter(ausbeer)
Qtr1 Qtr2 Qtr3 Qtr4
1956 0.67 0.33 0.00 0.00
1957 0.00 1.00 0.00 0.00
1958 0.00 1.00 0.00 0.00
1959 1.00 0.00 0.00 0.00
1960 0.00 1.00 0.00 0.00
1961 0.33 0.67 0.00 0.00
1962 0.00 1.00 0.00 0.00
1963 0.00 1.00 0.00 0.00
1964 1.00 0.00 0.00 0.00
1965 0.00 1.00 0.00 0.00
1966 0.00 1.00 0.00 0.00
1967 1.00 0.00 0.00 0.00
1968 0.00 1.00 0.00 0.00
1969 0.00 1.00 0.00 0.00
5번
a. 시간 그래프 그리고 설명하기 증가하는 트렌드, 그리고 반기마다 치솟는 계절성이 있다.
b. 반기마다 치솟는 계절성이 있는데, 치솟는 정도가 일정하지 않고 시간이 지날수록 더 많이 치솟고 있어, 선형적인 모양은 아니다. 따라서 로그를 씌워 이 값을 좀더 선형적으로 만들어 볼 수 있다.
c. 데이터에 로그를 취한 것을 선형추세, 계절성 가변수, 서핑축제 가변수를 넣어 맞춰 보기
> fit.fancy <- tslm(log(fancy) ~ trend + season)
> summary(fit.fancy)
Call:
tslm(formula = log(fancy) ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-0.41644 -0.12619 0.00608 0.11389 0.38567
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.6058604 0.0768740 98.939 < 2e-16 ***
trend 0.0223930 0.0008448 26.508 < 2e-16 ***
season2 0.2510437 0.0993278 2.527 0.013718 *
season3 0.6952066 0.0993386 6.998 1.18e-09 ***
season4 0.3829341 0.0993565 3.854 0.000252 ***
season5 0.4079944 0.0993817 4.105 0.000106 ***
season6 0.4469625 0.0994140 4.496 2.63e-05 ***
season7 0.6082156 0.0994534 6.116 4.69e-08 ***
season8 0.5853524 0.0995001 5.883 1.21e-07 ***
season9 0.6663446 0.0995538 6.693 4.27e-09 ***
season10 0.7440336 0.0996148 7.469 1.61e-10 ***
season11 1.2030164 0.0996828 12.068 < 2e-16 ***
season12 1.9581366 0.0997579 19.629 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1858 on 71 degrees of freedom
Multiple R-squared: 0.9527, Adjusted R-squared: 0.9447
F-statistic: 119.1 on 12 and 71 DF, p-value: < 2.2e-16
더보기
서핑축제 어떻게 넣나 했는데 인터넷의 누군가는 이렇게 해결했구나.
# Create Festival Dummy Variable - Note: Skip first year.
d_festival <- rep(0, length(fancy))
d_festival[seq_along(d_festival)%%12 == 3] <- 1
d_festival[3] <- 0
# Make Dummy Variable a Timeseries
d_festival <- ts(d_festival, freq = 12, start=c(1987,1))
fancydata <- data.frame(fancy.log, d_festival)
# Create Model
fancyfit <- tslm(fancy.log ~ trend + season + d_festival, data = fancydata)
d. 시간에 따른 잔차, 적합값에 따른 잔차 그려보기.
시간에 따른 잔차는 아래처럼 나온다. 일반적으로 0을 중심으로 왔다갔다 하고 있다.
적합값에 따른 잔차는 아래처럼 나온다. 특별한 패턴은 안보이므로, 잘 맞춰져 있는 것 같다.
> res <- residuals(fit.fancy)
> plot(res, xlab="year", ylab="residuals")
>
> fit_values <- as.vector(fitted(fit.fancy))
> fit_res <- as.vector(residuals(fit.fancy))
> plot(fit_values, fit_res, xlab="fitted values", ylab="residuals")
e. 월별 잔차에 대해 박스플롯 그려보기. 모델에 문제가 있는가?
> res.month <- data.frame(res, month=rep(1:12, 7))
> boxplot(res ~month, data=res.month)
f. 계수값들이 각 변수에 대해 말하고 있는 것은?
데이터에 로그를 씌웠기 때문에 1:1 스케일로 해석하기는 어렵다. 하지만 계수값을 통해 해당 변수가 예측하고자 하는 값에 어느정도로 큰 영향을 미치는지, 그 영향이 양인지 음인지 정도는 알 수 있다.
g. 브로쉬-갓프리?
p-value가 5%보다 작아서 귀무가설이 기각된다. 잔차에 자기상관관계가 있다.
> checkresiduals(fit.fancy)
Breusch-Godfrey test for serial correlation of order up to 17
data: Residuals from Linear regression model
LM test = 37.152, df = 17, p-value = 0.003209
h. 모델을 사용해 94~96년의 월별 매출을 예측해보기. 예측구간도 만들어보기
> future_data <- data.frame(dummy_festival = rep(0, 36))
> fc_fit_fancy <- forecast(fit.fancy, newdata=future_data)
>
> autoplot(fc_fit_fancy)
> fc_fit_fancy
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 1994 9.509264 9.246995 9.771532 9.105002 9.913525
Feb 1994 9.782700 9.520432 10.044969 9.378439 10.186962
Mar 1994 10.249256 9.986988 10.511525 9.844995 10.653518
Apr 1994 9.959377 9.697108 10.221645 9.555115 10.363638
May 1994 10.006830 9.744562 10.269098 9.602569 10.411091
Jun 1994 10.068191 9.805923 10.330459 9.663930 10.472452
Jul 1994 10.251837 9.989569 10.514105 9.847576 10.656099
Aug 1994 10.251367 9.989099 10.513635 9.847106 10.655628
Sep 1994 10.354752 10.092484 10.617020 9.950491 10.759014
Oct 1994 10.454834 10.192566 10.717102 10.050573 10.859095
Nov 1994 10.936210 10.673942 11.198478 10.531949 11.340471
Dec 1994 11.713723 11.451455 11.975991 11.309462 12.117984
Jan 1995 9.777979 9.512777 10.043182 9.369195 10.186763
Feb 1995 10.051416 9.786214 10.316618 9.642632 10.460200
Mar 1995 10.517972 10.252770 10.783174 10.109188 10.926756
Apr 1995 10.228092 9.962890 10.493295 9.819308 10.636876
May 1995 10.275546 10.010343 10.540748 9.866762 10.684330
Jun 1995 10.336907 10.071704 10.602109 9.928123 10.745691
Jul 1995 10.520553 10.255351 10.785755 10.111769 10.929337
Aug 1995 10.520083 10.254880 10.785285 10.111299 10.928867
Sep 1995 10.623468 10.358266 10.888670 10.214684 11.032252
Oct 1995 10.723550 10.458347 10.988752 10.314766 11.132334
Nov 1995 11.204926 10.939723 11.470128 10.796142 11.613710
Dec 1995 11.982439 11.717236 12.247641 11.573655 12.391223
Jan 1996 10.046695 9.777950 10.315440 9.632451 10.460940
Feb 1996 10.320132 10.051387 10.588877 9.905887 10.734376
Mar 1996 10.786688 10.517943 11.055433 10.372443 11.200932
Apr 1996 10.496808 10.228063 10.765553 10.082564 10.911053
May 1996 10.544261 10.275516 10.813007 10.130017 10.958506
Jun 1996 10.605623 10.336878 10.874368 10.191378 11.019867
Jul 1996 10.789269 10.520524 11.058014 10.375024 11.203513
Aug 1996 10.788798 10.520053 11.057543 10.374554 11.203043
Sep 1996 10.892184 10.623439 11.160929 10.477939 11.306428
Oct 1996 10.992266 10.723521 11.261011 10.578021 11.406510
Nov 1996 11.473641 11.204896 11.742386 11.059397 11.887886
Dec 1996 12.251155 11.982409 12.519900 11.836910 12.665399
i. 예측값과 구간을 변형해 원데이터와 스케일 맞추기
> df <- as.data.frame(fc_fit_fancy)
> df <- exp(df)
> head(df, 10)
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 1994 13484.06 10373.35 17527.60 9000.202 20201.76
Feb 1994 17724.45 13635.50 23039.57 11830.533 26554.69
Mar 1994 28261.51 21741.71 36736.44 18863.703 42341.27
Apr 1994 21149.61 16270.49 27491.85 14116.722 31686.25
May 1994 22177.42 17061.19 28827.88 14802.755 33226.11
Jun 1994 23580.87 18140.87 30652.19 15739.516 35328.75
Jul 1994 28334.55 21797.90 36831.38 18912.452 42450.69
Aug 1994 28321.23 21787.65 36814.06 18903.560 42430.74
Sep 1994 31405.93 24160.73 40823.80 20962.508 47052.23
Oct 1994 34711.77 26703.92 45120.97 23169.053 52005.01
j. 어떻게 예측을 개선할 수 있는가?
위에 자기상관관계가 조금 있다고 나왔으니 그 부분을 좀더 파 보면 되지 않을까..!
6번
gasoline 시계열은 미국 완성품 모터 가솔린 제품 공급에 대한 주별 데이터임. 기간은 91년 2.2 ~ 17년 1.20. 2004년 말까지의 데이터만 가지고 살펴보기.
a. 추세를 고려하는 조화 회귀로 데이터를 맞춰봅시다. 푸리에 항의 수를 바꾸면서 실험해봅시다. 관측된 가솔린과 적합값을 그래프로 나타내보고 관찰한 것을 설명해봅시다.
푸리에 항의 갯수가 커질수록 예측하는 파란색 데이터 선의 진폭이 자잘해지는 모습. 5, 10, 20으로 비교했을 때는 20이 가장 원 데이터의 모습을 잘 따라가고 있습니다.
K=5
K=10
K=20
> fourier.gasoline <- tslm(gasoline2 ~ trend + fourier(gasoline2, K=5))
> fourier.gasoline2 <- tslm(gasoline2 ~ trend + fourier(gasoline2, K=10))
> fourier.gasoline3 <- tslm(gasoline2 ~ trend + fourier(gasoline2, K=15))
> fourier.gasoline4 <- tslm(gasoline2 ~ trend + fourier(gasoline2, K=20))
> fourier.gasoline5 <- tslm(gasoline2 ~ trend + fourier(gasoline2, K=25))
>
> autoplot(gasoline2, series="data") + autolayer(fitted(fourier.gasoline), series="fitted")
> autoplot(gasoline2, series="data") + autolayer(fitted(fourier.gasoline2), series="fitted")
> autoplot(gasoline2, series="data") + autolayer(fitted(fourier.gasoline4), series="fitted")
b. AICc 또는 CV 값을 최소화하여 넣어야 할 푸리에 항의 적절한 개수를 골라봅시다.
CV를 기준으로는 K = 10이 가장 적합하며, AICc를 보았을 때는 15가 가장 적합하다.
> CV(fourier.gasoline) # k=5
CV AIC AICc BIC AdjR2
0.0768856 -1796.2437141 -1795.7138742 -1737.0611115 0.8283319
> CV(fourier.gasoline2) # k=10
CV AIC AICc BIC AdjR2
7.281821e-02 -1.834743e+03 -1.833112e+03 -1.730035e+03 8.397505e-01
> CV(fourier.gasoline3) # k=15
CV AIC AICc BIC AdjR2
7.330875e-02 -1.830854e+03 -1.827489e+03 -1.680621e+03 8.410505e-01
> CV(fourier.gasoline4) # k=20
CV AIC AICc BIC AdjR2
7.343301e-02 -1.830709e+03 -1.824950e+03 -1.634951e+03 8.431449e-01
> CV(fourier.gasoline5) # k=25
CV AIC AICc BIC AdjR2
7.499553e-02 -1.817360e+03 -1.808513e+03 -1.576077e+03 8.422319e-01
c. checkresiduals() 함수를 사용하여 마지막 모델의 잔차를 확인해봅시다. 잔차가 상관 관계 테스트를 통과하지 못하더라도, 아마도 예측값과 예측 구간에 큰 차이를 만들어내는 정도로 심각하지는 않을 것입니다. (상관관계가 통계적으로 유의미하더라도 상대적으로 작다는 것에 주목합시다.)
아래 그림을 보았을 때, 잔차의 등분산성을 확인 가능하다. ACF를 보았을때는 자기상관관계가 있지만 크게 두드러지지 않는다. 잔차의 평균은 0으로 봐도 될 것이다.
> checkresiduals(fourier.gasoline3)
Breusch-Godfrey test for serial correlation of order up to 104
data: Residuals from Linear regression model
LM test = 155.47, df = 104, p-value = 0.0008107
d. 조화 회귀로 예측하기 위해서, 푸리에 항의 미래 값을 예측할 필요가 있을 것입니다. 이러한 작업을 다음과 같이 할 수 있습니다.여기에서 fit은 tslm()으로 맞춘 모델이고, K는 fit을 만들 때 사용하는 푸리에 항의 개수, h는 필요한 예측 범위입니다. 데이터의 다음 연도를 예측해봅시다.
fc <- forecast(fit, newdata=data.frame(fourier(x,K,h)))e. 2005년의 실제 데이터에 따라 예측값 그려보기.
7번
huron 데이터는 휴론 강의 1875년부터 1972년까지의 수위기록이 피트 단위로 있음.
a. 데이터 그려보고 특징 말해보기
계절성..은 잘 모르겠고, 추세는 장기간 감소하다가 1960년대 이후 일시적 상승하는 느낌..?
b. 선형 회귀로 맞춰보고 1925년에 매듭이 있는 조각별 선형 추세모델과 비교
c. 위 두가지 모델에서 1980년까지의 예측값을 얻어보고 얻은 예측값 설명하기
> huron
Time Series:
Start = 1875
End = 1972
Frequency = 1
[1] 10.38 11.86 10.97 10.80 9.79 10.39 10.42 10.82 11.40 11.32 11.44 11.68 11.17 10.53 10.01 9.91
[17] 9.14 9.16 9.55 9.67 8.44 8.24 9.10 9.09 9.35 8.82 9.32 9.01 9.00 9.80 9.83 9.72
[33] 9.89 10.01 9.37 8.69 8.19 8.67 9.55 8.92 8.09 9.37 10.13 10.14 9.51 9.24 8.66 8.86
[49] 8.05 7.79 6.75 6.75 7.82 8.64 10.58 9.48 7.38 6.90 6.94 6.24 6.84 6.85 6.90 7.79
[65] 8.18 7.51 7.23 8.42 9.61 9.05 9.26 9.22 9.38 9.10 7.95 8.12 9.75 10.85 10.41 9.96
[81] 9.61 8.76 8.18 7.21 7.13 9.10 8.25 7.91 6.89 5.96 6.80 7.68 8.38 8.52 9.74 9.31
[97] 9.89 9.96
> ## linear regression model
> huron.fit <- tslm(huron ~ trend)
> ## create piecewise linear regression
> huron.t <- time(huron) ## huron의 시간축 불러온다
> huron.t.break <- 1925 ## huron의 시간축 중 1925년에 knot를 지정해준다.
> ## knot 앞뒤로 시계열 데이터 생성하고 piecewise linear regression fit
> huron.pt <- ts(pmax(0, huron.t - huron.t.break), start = 1925)
> fit.pw <- tslm(huron ~ huron.t + huron.pt)
>
> ## 예측값을 제한 값 피팅
> autoplot(huron)
+ autolayer(fitted(huron.fit), series="linear")
+ autolayer(fitted(fit.pw), series ="piecewise")
>
> h <- 10
> fcsts.linear <- forecast(huron.fit, h=h) ## linear regression forecast
> ## piecewise의 예측값은 전체 데이터와 knot로 구한 데이터 두개를 합쳐서 구한다.
> huron.t.new <- huron.t[length(huron.t)] + seq(h)
> huron.pt.new <- huron.pt[length(huron.pt)] + seq(h)
> huron.new <- cbind(huron.t = huron.t.new, huron.pt = huron.pt.new) %>% as.data.frame()
> fcsts.piecewise <- forecast(fit.pw, newdata=huron.new)
> ## plot the models & forecasts
> autoplot(huron) + autolayer(fitted(huron.fit), series="linear")
+ autolayer(fitted(fit.pw), series ="piecewise")
+ autolayer(fcsts.linear, series="linear", PI=FALSE)
+ autolayer(fcsts.piecewise, series="piecewise", PI=FALSE)
+ xlab("year") + ylab("water level")
+ ggtitle("Huron Lake Data")
+ guides(color = guide_legend(title =" "))반응형
'Math & Statistics > Forecasting: Principles and Practice' 카테고리의 다른 글
| 챕터 6 시계열 분해 - 파트 2 (0) | 2022.03.19 |
|---|---|
| 챕터 6 시계열 분해 - 파트 1 (0) | 2022.03.12 |
| 챕터 5 - 상관관계, 인과관계, 예측 (0) | 2022.03.05 |
| 챕터 5 - 회귀분석, 비선형 회귀 (0) | 2022.03.05 |
| 챕터 5. 예측변수 선택 (0) | 2022.02.26 |