모바일은 화면을 돌려 가로화면으로 보시는 게 읽으시기 편할 수 있습니다. 돌려서 보시는 걸 추천드릴게요!!
🕑 오차의 자기상관 해결 02. 회귀가변수 추가
🕑 1. 트렌드 (trend)
1차 자기상관이 심하게 있는 경우 고려해 줘야 하는 사안입니다.
다만 $X$변수에도 trend가 있기 때문에 회귀모형 적합 결과 오차항에는 trend가 남지 않을 수도 있습니다.
하지만 오차항에 trend 자기상관이 남아있다면, trend에 해당하는 가변수를 $X$변수로 추가 가능합니다.
트렌드 가변수로는 시간 가변수 $𝑡$ 또는 $𝑡^2$ 를 사용합니다.
🕑 트렌드 파악
우선 데이터에서 각 feature의 plot을 그려 트렌드 유무를 파악합니다.
df1 <- read.csv("Table6_1.csv")
## 사용할 부분만 ts데이터로 변경
## 년도별 데이터로 주기 설정
df1.ts <- ts(df1[, 2:8], start = 1947, frequency = 1)
plot(df1.ts, col="red")
전반적으로 점점 증가하는 트렌드를 확인 가능합니다. 따라서, 모델 피팅에 시간 가변수를 넣어 트렌드를 잡겠습니다!!
🕑 모델 피팅
시간 가변수 time을 추가해서 모델을 피팅하였습니다.
model.ts3 <- lm(lnconsump ~ lndpi + lnwealth + interest + time, data = df1)
summary(model.ts3)>>
Call:
lm(formula = lnconsump ~ lndpi + lnwealth + interest + time,
data = df1)
Residuals:
Min 1Q Median 3Q Max
-0.020318 -0.008144 0.001155 0.006197 0.025922
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6640849 0.3513251 1.890 0.06465 .
lndpi 0.7212201 0.0303831 23.738 < 2e-16 ***
lnwealth 0.1369181 0.0255755 5.353 2.28e-06 ***
interest -0.0024247 0.0007032 -3.448 0.00117 **
time 0.0051831 0.0015989 3.242 0.00214 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01094 on 49 degrees of freedom
Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
F-statistic: 3.377e+04 on 4 and 49 DF, p-value: < 2.2e-16
모델 피팅 결과, 시간 가변수 time이 유의미하다고 나왔습니다.
🕑 모델 성능 파악
☑️ 잔차 plot
plot(model.ts3$residuals,type="b")
☑️ acf plot
acf(model.ts3$residuals)
☑️ 더빈왓슨검정
dwtest(model.ts3)>>
Durbin-Watson test
data: model.ts3
DW = 1.3364, p-value = 0.001143
alternative hypothesis: true autocorrelation is greater than 0
durbinWatsonTest(model.ts,max.lag=3)
durbinWatsonTest(model.ts2,max.lag=3)
durbinWatsonTest(model.ts3,max.lag=3)>>
lag Autocorrelation D-W Statistic p-value
1 0.29775555 1.289232 0.000
2 -0.02840612 1.933710 0.678
3 0.01631809 1.678649 0.224
Alternative hypothesis: rho[lag] != 0
>>
lag Autocorrelation D-W Statistic p-value
1 -0.046323882 1.896780 0.680
2 -0.079590916 1.906867 0.932
3 0.008091062 1.718185 0.442
Alternative hypothesis: rho[lag] != 0
>>
lag Autocorrelation D-W Statistic p-value
1 0.26429515 1.336395 0.002
2 -0.11161897 2.086674 0.930
3 -0.09501358 1.926016 0.718
Alternative hypothesis: rho[lag] != 0
☑️ bg test
bgtest(model.ts3,order=1)
bgtest(model.ts3,order=3)>>
Breusch-Godfrey test for serial correlation of order up to 1
data: model.ts3
LM test = 4.4559, df = 1, p-value = 0.03478
>>
Breusch-Godfrey test for serial correlation of order up to 3
data: model.ts3
LM test = 6.6822, df = 3, p-value = 0.08275
트렌드에 의한 영향을 잡아주었음에도, 작기는 하지만 모델의 오차에 자기상관이 남아있습니다.
계절성은 없는 데이터이기 때문에 해결을 위해서는 차분 혹은 다른 방법을 사용해주면 해결은 가능할 것으로 보입니다.
🕑 2. 계절성 (seasonality)
12차 혹은 4차 자기상관이 심하게 있는 경우를 의미합니다.
X변수에도 계절성이 있다면 회귀모형 적합을 통해 오차항에는 계절성이 남지 않을 수도 있지만,
그렇지 않다면 계절성에 해당하는 가변수를 X변수로 추가 가능합니다.
계절가변수로는 월별 혹은 계절 더미 변수를 사용합니다.
plot(AirPassengers)
AirPassengers 데이터 plot 결과, 각 년도의 특정 시점에 올라갔다가 내려가는 주기적인 움직임을 보여주고 있습니다.
각 년도의 일정한 시기에 특정 움직임이 반복되고 있기 때문에, 계절성이 있다고 할 수 있습니다.
년도 별로 1~12 월의 계절가변수와 시간가변수를 넣어준 모델을 생성해 보겠습니다!!
## 계절더미변수 month
month<-rep(c(1:12), 12)
# 시간더미변수 time
time<-c(1:length(AirPassengers))
model<-lm(AirPassengers ~ as.factor(month) + time)
summary(model)>>
Call:
lm(formula = AirPassengers ~ as.factor(month) + time)
Residuals:
Min 1Q Median 3Q Max
-42.121 -18.564 -3.268 15.189 95.085
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 63.50794 8.38856 7.571 5.88e-12 ***
as.factor(month)2 -9.41033 10.74941 -0.875 0.382944
as.factor(month)3 23.09601 10.74980 2.149 0.033513 *
as.factor(month)4 17.35235 10.75046 1.614 0.108911
as.factor(month)5 19.44202 10.75137 1.808 0.072849 .
as.factor(month)6 56.61502 10.75254 5.265 5.58e-07 ***
as.factor(month)7 93.62136 10.75398 8.706 1.17e-14 ***
as.factor(month)8 90.71103 10.75567 8.434 5.32e-14 ***
as.factor(month)9 39.38403 10.75763 3.661 0.000363 ***
as.factor(month)10 0.89037 10.75985 0.083 0.934177
as.factor(month)11 -35.51996 10.76232 -3.300 0.001244 **
as.factor(month)12 -9.18029 10.76506 -0.853 0.395335
time 2.66033 0.05297 50.225 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 26.33 on 131 degrees of freedom
Multiple R-squared: 0.9559, Adjusted R-squared: 0.9518
F-statistic: 236.5 on 12 and 131 DF, p-value: < 2.2e-16
모델 summary 결과 6,7,8,9 월 계절 가변수가 유의미하다고 나왔습니다. 계절성을 잡은 것이라 생각해도 될 것 같습니다.
확실한 검증을 위해서, acf plot을 그려보겠습니다.
☑️ 원본데이터 acf plot
acf(AirPassengers)
☑️ 모델 acf plot
acf(model$residuals)
원본보다는 나아보이기는 하지만, 여전히 자기상관이 심하게 남아 있는 것 같습니다.
계절성을 어느 정도 잡았으니, 유의미한 변수만 남겨서 모델을 피팅한다면 좋은 모델이 될 것 같습니다.
🕑 3. 미국소비함수 데이터 분석
1973-2011년까지의 미국 신규주택건설에 대한 데이터입니다.
2007년 서브프라임모기지 사태로 인한 금융위기로 뚝 떨어지는 부분을 가진 데이터입니다.
• Hstart : 신규주택건설, 계절 조정 (단위:천)
• UN : 계절조정된 민간 실업률 (단위: %)
• M2: 계절조정된 M2 공급량 (단위:10억달러)
• Mgrate: 신규주택 모기지론 이자율 (단위: %)
• Primerate: 은행의 프라임레이트 (최대우대금리)(단위: %)
• RGDP: 실질 GDP (2005년 기준), 계절조정
📉 01. 테스트 모델 피팅
모델 피팅 전에, 데이터 Y변수의 분포를 알아보았습니다.
역시나 2007년과 2008년에 해당하는 부분에서 뚝 떨어지는 것을 확인할 수 있습니다.
data <- read.csv("Table6_10.csv")
head(data)
plot(data$hstart)
m <- lm(hstart ~ un + m2 + mgrate + primerate + rgdp, data = data)
summary(m)>>
Call:
lm(formula = hstart ~ un + m2 + mgrate + primerate + rgdp, data = data)
Residuals:
Min 1Q Median 3Q Max
-510.05 -148.54 -51.33 201.00 595.46
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1844.7354 740.3234 2.492 0.01791 *
un -169.7104 59.4540 -2.854 0.00739 **
m2 -0.1762 0.1254 -1.405 0.16924
mgrate 172.6015 67.5902 2.554 0.01547 *
primerate -129.6958 40.8864 -3.172 0.00326 **
rgdp 0.1144 0.1041 1.099 0.27975
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 268.6 on 33 degrees of freedom
Multiple R-squared: 0.5941, Adjusted R-squared: 0.5326
F-statistic: 9.66 on 5 and 33 DF, p-value: 9.443e-06
모델의 예측값에 따른 잔차의 plot을 통해 데이터의 어느 부분에서 잔차가 크게 나타났는지 확인해 보겠습니다.
plot(m$residuals~m$fitted.values)
abline(h = 0)
모델의 예측값이 600정도인 부분은 아마 금융 위기 시점의 데이터인 것으로 보입니다. 하지만 이러한 이상치에 대해서 크지 않은 잔차를 보여주었습니다. 예측값이 큰 부분에 있어서는 잔차가 어느 정도 있었지만, 전반적으로 비슷한 분포를 보여주고 있기 때문에 큰 문제는 되지 않을 것 같습니다.
모델의 잔차를 알아보았으니, 이번에는 acf plot을 통해 자기상관성을 알아보겠습니다.
acf(m$residuals)
아직은 1,2차 자기상관과 여러 자기상관이 남아있기 때문에, 이를 처리해주기 위해 몇 가지 변수 변환을 해보도록 하겠습니다.
📉 02. 변수 변환
## log
data$lnhstart <- log(data$hstart)
data$lnun <- log(data$un)
data$lnmgrate <- log(data$mgrate)
data$lnprimerate <- log(data$primerate)
data$lnrgdp <- log(data$rgdp)
data$lnm2 <- log(data$m2)
## 차분
newhstart <- diff(data$hstart, differences = 1)
newun <- diff(data$un, difference = 1)
newmgrate <- diff(data$mgrate, difference = 1)
newprimerate <- diff(data$primerate, differences = 1)
newm2 <- diff(data$m2, differences = 1)
newrgdp <- diff(data$rgdp, differences = 1)
newtime <- diff(data$time2, differences = 1)
## 로그변환 후 차분
newlnhstart <- diff(data$lnhstart, differences = 1)
newlnun <- diff(data$lnun, difference = 1)
newlnmgrate <- diff(data$lnmgrate, difference = 1)
newlnprimerate <- diff(data$lnprimerate, differences = 1)
newlnm2 <- diff(data$lnm2, differences = 1)
newlnrgdp <- diff(data$lnrgdp, differences = 1)
☑️ 차분 모델 해석 : $변수의 전 시점 대비 변화량으로 해석 (원래 $Y$ 변수의 단위 유지).
☑️ log 변환 모델 해석 : $변수의 증가율(%)로 해석.
☑️ log log 변환 모델 해석 : $변수 증가율의 증가율, 즉 %의 증가율로 해석. %p 단위.
☑️ log & 차분 모델 해석 : $변수의 전 시점 대비 변화율로 해석 (단위 유지 안됨, %로 변환가능(×100))
☑️ $Y$ 변수가 금융상품의 가격이라면 수익률 (로그 수익률이라고 부르기도 함)로 해석가능.
📉 03. 모델피팅
☑️ 차분한 데이터로 모델 피팅
m3 <- lm(newlnhstart ~ newm2 + newrgdp + newprimerate -1)
summary(m3)>>
Call:
lm(formula = newlnhstart ~ newm2 + newrgdp + newprimerate - 1)
Residuals:
Min 1Q Median 3Q Max
-0.24832 -0.12661 -0.04536 0.02798 0.36363
Coefficients:
Estimate Std. Error t value Pr(>|t|)
newm2 -0.0005333 0.0001075 -4.960 1.81e-05 ***
newrgdp 0.0005471 0.0001082 5.055 1.36e-05 ***
newprimerate -0.0491342 0.0126692 -3.878 0.000443 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1413 on 35 degrees of freedom
Multiple R-squared: 0.4947, Adjusted R-squared: 0.4514
F-statistic: 11.42 on 3 and 35 DF, p-value: 2.256e-05
☑️ 잔차 plot
plot(m3$residuals~m3$fitted.values)
abline(h = 0)
☑️ acf plot
acf(m3$residuals)
☑️ 더빈왓슨검정
dwtest(m3)>>
Durbin-Watson test
data: m3
DW = 1.1384, p-value = 0.001679
alternative hypothesis: true autocorrelation is greater than 0
☑️ 변수 간 상관계수 검정
cor(newm2,newrgdp)
cor(newm2,newprimerate)
cor(newrgdp,newprimerate)>> [1] -0.008227085
>> [1] -0.2242671
>> [1] 0.2920121
차분한 데이터로 모델을 피팅한 결과, acf로 알아본 자기상관 역시 크게 개선되었고 변수 간 상관관계 역시 적은 것으로 나타났습니다. 보다 좋은 모델을 피팅하기 위해서, 몇 가지 변수를 만들어보겠습니다.
📉 04. time 트렌드 변수 추가
먼저, 시간변수 $t$와 $t^2$을 만들어 모델에 추가하겠습니다.
## t
data$time<-c(1:length(data$hstart))
## t^2
data$time2<-c(1:length(data$hstart))^2
m4 <- lm(hstart ~ un + mgrate + primerate + time2 + time, data = data)
summary(m4)>>
Call:
lm(formula = hstart ~ un + mgrate + primerate + time2 + time,
data = data)
Residuals:
Min 1Q Median 3Q Max
-665.24 -122.79 9.96 120.15 586.46
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2669.9097 330.9011 8.069 2.60e-09 ***
un -275.7095 48.2457 -5.715 2.24e-06 ***
mgrate 270.5205 80.6590 3.354 0.002013 **
primerate -174.5266 44.8454 -3.892 0.000457 ***
time2 0.8167 0.5480 1.490 0.145621
time -36.2996 18.8691 -1.924 0.063046 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 262.3 on 33 degrees of freedom
Multiple R-squared: 0.613, Adjusted R-squared: 0.5543
F-statistic: 10.45 on 5 and 33 DF, p-value: 4.491e-06
☑️ acf plot
acf(m4$residuals)
다만 시간 변수의 추가만으로는 유의미하지도 않고 acf plot 결과도 좋지 않았습니다.
다른 방법을 사용해야 할 듯 합니다.
📉 05. 데이터 정보 범주형 변수 추가
언급한 대로, 해당 데이터는 금융위기가 포함된 데이터입니다. 따라서, 해당 부분에 대한 정보를 가지고 있는 가변수를 추가하여 이 부분의 영향을 고려한 모델을 만들어 보겠습니다. 금융위기에 의한 영향을 받은 2007년 이후의 데이터에 대해서는 1의 값을, 아닌 데이터는 0을 넣어주는 가변수 변환을 진행해주었습니다.
## structural change in HSTART
## 2007년 이전 / 이후에 대한 더미 변수 생성
## 2007년 -> 서브프라임모기지에 의한 금융위기
dd<-data$year
dd<-ifelse(data$year >= 2007, 1, 0)
## 데이터에 회귀가변수 GFC 추가
data$GFC = dd
m5 <- lm(hstart~ un + mgrate + primerate + GFC, data = data)
summary(m5)>>
Call:
lm(formula = hstart ~ un + mgrate + primerate + GFC, data = data)
Residuals:
Min 1Q Median 3Q Max
-545.76 -145.36 -17.84 202.15 469.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2256.95 197.61 11.421 3.49e-13 ***
un -164.29 45.61 -3.602 0.000997 ***
mgrate 145.41 61.72 2.356 0.024387 *
primerate -110.56 39.76 -2.781 0.008782 **
GFC -477.65 183.93 -2.597 0.013797 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 253 on 34 degrees of freedom
Multiple R-squared: 0.6291, Adjusted R-squared: 0.5855
F-statistic: 14.42 on 4 and 34 DF, p-value: 5.562e-07
☑️ acf plot
acf(m5$residuals)
새로 추가한 변수 GFC가 유의미하긴 했지만, 아직 자기상관은 남아있습니다. 해당 모델에 앞서 만든 시간가변수도 추가해서 모델을 피팅해보겠습니다.
📉 06. 타임트렌드 / 범주형 더미변수 추가
## 자기상관 해결
m6 <- lm(hstart ~ un + mgrate + primerate + time2 + time + GFC, data = data)
summary(m6)>>
Call:
lm(formula = hstart ~ un + mgrate + primerate + time2 + time +
GFC, data = data)
Residuals:
Min 1Q Median 3Q Max
-408.86 -117.74 23.91 109.14 386.04
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2361.7449 278.3356 8.485 1.07e-09 ***
un -222.3638 41.1425 -5.405 6.12e-06 ***
mgrate 265.0311 65.5132 4.045 0.000308 ***
primerate -161.1909 36.5525 -4.410 0.000109 ***
time2 2.0391 0.5299 3.848 0.000536 ***
time -68.0919 17.0533 -3.993 0.000357 ***
GFC -828.1515 194.9721 -4.248 0.000174 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 213 on 32 degrees of freedom
Multiple R-squared: 0.7525, Adjusted R-squared: 0.7061
F-statistic: 16.22 on 6 and 32 DF, p-value: 1.785e-08
☑️ acf plot
acf(m6$residuals)
시간 가변수와 데이터 세부사항 가변수를 모두 추가한 결과 정말 괜찮은 모델이 만들어졌습니다.
이렇게 여러 방식의 회귀가변수 추가를 통해서 자기상관을 해결하고, 모델의 성능을 정말 크게 향상시킬 수 있습니다.
지금까지 자기상관 해결을 위해 차분과 회귀가변수 추가의 방법을 사용해봤습니다.
다음 포스팅에서는 자기상관 해결을 위한 마지막 방법인 시계열 모형을 사용해보겠습니다.
'📈📉 비즈니스 어낼리틱스 > 🕑 시계열 분석' 카테고리의 다른 글
| 🕑시계열 데이터 분석 12 - 자기상관 해결 4. 자기회귀시차분포모형(ARDL Model) (2) | 2023.04.25 |
|---|---|
| 🕑시계열데이터 분석 11 - 자기상관 해결 3. 자기회귀모형 (AR Model) (0) | 2023.04.21 |
| 🕑시계열데이터 분석 09 - 자기상관 해결 1. 차분 (2) | 2023.04.20 |
| 🕑시계열데이터 분석 08 - 시계열 회귀분석 (0) | 2023.04.15 |
| 🕑시계열데이터 분석 07 - 평활기법 (0) | 2023.04.14 |