모바일은 화면을 돌려 가로화면으로 보시는 게 읽으시기 편할 수 있습니다. 돌려서 보시는 걸 추천드릴게요!!
🕑 오차의 자기상관 해결 03. AR 모델 : 자기회귀모형
🕑 오차의 자기상관 해결방법
1. 변수변환 & 회귀분석
2. 회귀 가변수를 이용한 회귀모형
3. 시계열 모형 : 선형회귀모형의 형태 : 시간이 다른 변수값(과거값)을 사용한다는 점에서 다름
• 자기회귀모형 (Autoregressive, AR model)
- y 변수를 과거의 y값으로 적합. .
• 자기회귀시차분포모형 (Autoregressive Distributed Lag –ARDL)
- y변수를 x변수와 과거의 x값, 과거의 y값으로 적합.
• 자기상관오차회귀모형 (Regression model with autoregressive error)
- y변수를 X변수로 적합. 오차항의 자기상관은 오차항의 AR로 적합.
• 그 외 시계열 모형 (Moving average (MA), ARMA, ARIMA)
4. 회귀모형에서 회귀계수의 추정방법을 바꾸는 방법
☑️ Y 변수만 분석할 경우 사용하는 방법
• 시계열분해기법(time series decomposition method)
• 지수평활법(exponential smoothing)
• ARIMA(autoregressive integrated moving average)
☑️ X-Y 변수의 관계를 분석할 경우 -> X의 변화의 따른 Y의 변화 파악
• 회귀모형
• 시계열회귀모형
• 자기회귀모형
• 자기회귀시차분포모형
• 자기상관오차회귀모형
• 전이함수잡음모형(transfer function noise model)
• 개입분석(intervention analysis)
이 중에서 자기회귀모형, 자기회귀시차분포모형, 자기상관오차회귀모형을 시계열 모형이라고 합니다.
🚩 시계열 모형을 위한 가정 - 안정적 시계열 : Stationary Condition (정상성 조건)
- a) 평균이 해당기간 동안 일정해야 하며, 트렌드가 없어야 합니다.
- b) 분산이 해당기간 동안 일정해야 합니다.
- c) 공분산 혹은 자기상관은 시차에만 의존해야 합니다. 즉, 일정 시차를 두고 반복되는 추세를 보줘야 합니다.
- d) 위 가정을 만족하지 못하면 변수변환 (로그, 차분)을 통해 정상시계열로 만들어 다음의 모델을 사용합니다.
• AR
• DL
• Regression model with autoregressive error
• 혹은 좀더 복잡한 시계열 모형을 사용하여 해결
🕑 AR 모델 : 자기회귀모형 (Autoregressive model, AR model)
ARIMA 모형의 한 종류로써, 종속변수에 대한 회귀모형의 독립변수로 종속변수의 과거값만을 이용한 모형입니다.
$AR(p) : Y_t = β_0 + β_1Y_{t-1}+β_2Y_{t-2}...,+β_pY_{t-p}+u_t$
p-차 자기회귀모형으로, 종속변수의 p 시점 이전까지 모든 과거값을 독립변수로 사용하는 AR모형입니다. 일반적인 선형회귀에 사용하는 추정법인 OLS가 아닌 MLE(Maximum Likelihood Estimation)니 CSS-ML(conditional sum of squaredMLE)을 사용합니다.
🕑 AR모형에서 잔차의 자기상관 검정
1차 자기상관만을 보여주는 Dw test 의 한계를 극복하는 BoxPierce, LjungBox test를 사용합니다.
여러 개의 자기상관계수를 한꺼번에 검정하는 방법입니다.
자기상관이 없다는 귀무가설을 가지고 있기 때문에, 검정의 p-value가 0.05보다 커야 좋은 모델입니다.
🕑 예제: 미국의 인플레이션율
1983년 12월 부터 2006년 5월까지의 미국의 소비자 물가지수를 이용해 구한 1984년 1월 부터 2006년 5월까지의 월간 인플레이션율 (monthly inflationrate) 데이터입니다. 분석을 위해 자기회귀 모형 (AR model) 사용할 것이며, 만약 독립변수인 임금변화율이 동기간에 대해 조사되었다면 자기회귀시차분포모형도 가능합니다.
data<-read.csv(file="inflation2.csv", header = FALSE)
## 년도, 월, 임금, 임금상승률, 소비자물가지수, 인플레이션율
names(data)<-c("YEAR", "MONTH", "WAGE", "PCWAGE", "CPI", "INFLN")
☑️ 독립변수 plot
plot(data$CPI)
원본 데이터의 경우 단순 CPI를 가지고 있는 데이터입니다.
하지만 plot 결과 트렌드가 명확해서 stationary condition을 만족하지 못하였기에, 이를 로그 변환하고 차분해서 조건을 만족시키는 것이 좋습니다. 즉, $ln(CPI_{t}) - ln(CPI_{t-1})$ 의 과정을 거쳐야 했는데, 이는 $(CPI_{t} - CPI_{t-1})/(CPI_{t-1})$, 즉 증가율과 유사합니다. 따라서 증가율에 대한 column을 계산해 독립변수의 로그변환과 차분을 대신할 것입니다. 그 열이 아래 보이는 INFLN 이고, 이 변수를 종속변수로 사용해서 AR 모델을 만들겁니다!!
summary(data)## YEAR MONTH WAGE PCWAGE
## Min. :1983 Min. : 1.000 Min. : 8.320 Min. :-0.2260
## 1st Qu.:1989 1st Qu.: 3.000 1st Qu.: 9.812 1st Qu.: 0.1688
## Median :1995 Median : 6.000 Median :11.530 Median : 0.2442
## Mean :1995 Mean : 6.456 Mean :11.958 Mean : 0.2572
## 3rd Qu.:2000 3rd Qu.: 9.000 3rd Qu.:14.155 3rd Qu.: 0.3635
## Max. :2006 Max. :12.000 Max. :16.620 Max. : 0.6565
## NA's :1
## CPI INFLN
## Min. :101.4 Min. :-0.6551
## 1st Qu.:124.5 1st Qu.: 0.1470
## Median :151.1 Median : 0.2454
## Mean :149.6 Mean : 0.2560
## 3rd Qu.:173.8 3rd Qu.: 0.3694
## Max. :201.9 Max. : 1.2158
## NA's :1☑️ 종속변수 plot
plot(data$INFLN,type='l')
## 시계열 모델 사용을 위해 ts데이터로 변환
inf <- ts(data$INFLN[2:270],start=c(1984,1), frequency=12)
plot(inf)
INFLN plot 결과 분산이 커지는 몇 부분이 발견되었습니다. 시장에 shock이 있었을 것이라 가정하고 전반적인 분산이 일정하다는 가정 하에 모델링을 진행할 것입니다. 우선 Y변수만 사용해서 모델을 구축할 것이기에, AR모델을 사용합니다. 차수를 먼저 결정해 주겠습니다!!
📉 01. AR모델의 차수 결정 : Pacf( ) 함수
데이터를 정상 시계열로 만든 후에는, 시계열 모형을 통해 분석 혹은 예측을 진행하면 됩니다!!
모든 시계열 모형은 자기상관성을 잡기 위해 몇 번째 차수까지 보정할지를 인수로 넣어줘야 합니다.
이 인수를 알아보기 위해서, acf 함수를 사용했던 선형회귀 모형과는 달리 AR 모델의 경우 Pacf 함수를 사용합니다.
partial autocorrelation을 계산해줍니다. 다른 시간의 $Y$ 값이 고정된 상태에서 $Y_t$와 $Y_{t-n}$ 간의 자기상관을 나타냅니다.
acf가 튀는 차수까지 잡아서 AR모델의 차수를 정해줍니다. 또한 시계열 데이터의 경우 일반 acf가 아닌 forecast 패키지의 Pacf를 사용해야 x축이 차수로 변경되므로, 용이하게 사용할 수 있습니다.
☑️ acf plot
Acf(inf)
☑️ pacf plot
library(forecast)
Pacf(inf)
plot 결과, 1차, 2차, 6차, 12차에서 자기상관이 높게 나타났습니다. 어느 차수까지 잡아줘야 할지 정해보겠습니다.
📉 02. 자동으로 적절한 ARMA 차수를 잡아주는 함수 : auto.arima(Y) 함수
R 패키지에는 자동으로 적절한 차수를 잡아주는 auto.arima() 함수가 존재합니다.
AR 모델의 경우 종속변수를 사용해서 종속변수를 피팅하기 때문에, 인수로 종속변수 $Y$를 넣어주면 됩니다.
이 함수는 log likelihood가 최소가 되게끔 차수를 결정해줍니다.
## AR(2) = Y_t = β0 + β1Y_t-1 + β2Y_t-2 + ε_t
library(forecast)
auto.arima(inf)## Series: inf
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## 0.2858 -0.2533 -0.9237
## s.e. 0.0653 0.0639 0.0362
##
## sigma^2 = 0.03862: log likelihood = 56.26
## AIC=-104.51 AICc=-104.36 BIC=-90.15함수 사용 결과, ARIMA(2,1,1)로 차수를 결정해주었습니다.
위 결과처럼, AR + I + MA 개념에서 함수의 결과 차수가 (2, 1, 1) 이면
• AR = 2 : AR모델의 차수로 2를 넣어줘야 합니다.
• I = 1 : 1차 차분 한번 진행. 만약 이 차수가 0이라면, 데이터가 정상성을 만족함을 의미합니다.
• MA = 1 : moving average 모델의 차수로 1을 넣어줘야 합니다.
차수로 (2,0,0) 을 넣어주면 본 모델은 AR(2) 모델이 됩니다.
만약 (2,1,1) 을 넣어주면 ARIMA(2,1,1) 모델이 되고,
(0,1,1) 을 넣어주면 IMA(1,1) 모델이 되는 겁니다.
생각보다 단순하죠 😗??
이러한 방법으로 ARIMA 모델의 차수를 정해줍니다. 하지만 MA 모형은 해석이 어려워 잘 사용하지는 않습니다.
위 결과를 보면, AR 모델의 경우 차수 2를 주고 차분은 한 번 진행하며, MA 차수로는 1을 넣어줘야 합니다.
다만, INFLN을 구하면서 이미 차분을 한번 진행하였으므로 차분은 하지 않습니다.
이번 포스팅에서 다루는 내용은 AR 모델이기 때문에 모델에 차수로 (2,0,0) 만 넣어주도록 하겠습니다!!
하지만 함수가 AR차수를 잡고 MA차수를 잡았기에 AR모델만으로는 자기상관이 완전히 잡히지 않을 것 같기는 합니다.
📉 03. ARIMA model fit : arima( )함수
$AR(2) = Y_t = β0 + β1Y_{t-1} + β2Y_{t-2} + ε_t$
arima(Y, order = c(AR차수, I차수, MA차수)) 형태로 함수를 사용합니다.
예를 들어, 차수로 (2,0,0)을 넣어주었다면 단순 AR(2) 모델을 사용하게 됩니다.
모델을 summary 해도 p-value는 주지 않아 유의미 정도를 파악하기 쉽지는 않지만, coefficient의 절댓값을 standard error의 절댓값으로 나눈 값이 2 이상이면 유의미하다고 할 수 있습니다.
이 모델이 완전히 적합했다면, 더 이상 오차항에 자기상관이 남아있으면 안됩니다.
Acf 함수로 자기상관 확인 후 만약 자기상관이 남아있다면, order를 바꿔 ARIMA 모델의 차수를 변경해줘야 합니다.
☑️ AR(2) 모형
model1<-arima(inf,order=c(2,0,0))
summary(model1)##
## Call:
## arima(x = inf, order = c(2, 0, 0))
##
## Coefficients:
## ar1 ar2 intercept
## 0.3619 -0.1821 0.2561
## s.e. 0.0602 0.0606 0.0147
##
## sigma^2 estimated as 0.03894: log likelihood = 54.76, aic = -101.53
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0001505561 0.1973391 0.1425864 -Inf Inf 0.8001155 0.01561765☑️ acf plot
Acf(model1$residuals)
모델의 결과식은 아래와 같습니다.
$AR(2) = Y_t = 0.25 + 0.36Y_{t-1} - 0.18Y_{t-2} + ε_t$
하지만 6차 자기상관이 남아있기에 AR 모델의 차수를 늘려 자기상관을 잡아줘야 합니다.
📉 04. AR 모델의 차수 증가
model2<-arima(inf,order=c(6,0,0))
summary(model2)##
## Call:
## arima(x = inf, order = c(6, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 intercept
## 0.3660 -0.2215 0.1021 0.0093 0.0351 0.1332 0.2581
## s.e. 0.0608 0.0647 0.0660 0.0667 0.0662 0.0623 0.0204
##
## sigma^2 estimated as 0.03759: log likelihood = 59.44, aic = -102.88
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0009482493 0.1938749 0.1392514 -Inf Inf 0.7814011 -0.01027995
☑️ acf plot
Acf(model2$residuals)
어느 정도 해결되기는 했지만, 차수를 6차까지 늘리면 유의미하지 않은 3,4,5차수를 사용해야 한다는 단점이 생깁니다.
자기상관 차수 검정을 통해 정말 보정이 필요한 차수를 알아보겠습니다.
📉 05. AR 모델 자기상관 검정 : Box.test(MODEL$residual, lag, type) 함수
p-value 형태로 자기상관의 유무를 알려주는 함수입니다.
lag 인수는 보려는 자기상관 계수의 범위를 지정해줍니다.
귀무가설이 자기상관이 없다는 것이기 때문에, 귀무가설을 채택하기 위해서 p-value가 0.05보다 크기를 원합니다.
☑️ 2차 Ljung Test
Box.test(model1$residual,lag=2,type="Ljung")##
## Box-Ljung test
##
## data: model1$residual
## X-squared = 0.28491, df = 2, p-value = 0.8672
☑️ 2차 Box-Pierce Test
Box.test(model1$residual,lag=2,type="Box-Pierce")##
## Box-Pierce test
##
## data: model1$residual
## X-squared = 0.28095, df = 2, p-value = 0.8689
☑️ 3차 Box-Pierce Test
Box.test(model1$residual,lag=3,type="Box-Pierce")##
## Box-Pierce test
##
## data: model1$residual
## X-squared = 2.9031, df = 3, p-value = 0.4068
☑️ 4차 Box-Pierce Test
Box.test(model1$residual,lag=4,type="Box-Pierce")##
## Box-Pierce test
##
## data: model1$residual
## X-squared = 3.4413, df = 4, p-value = 0.4869
☑️ 6차 Box-Pierce Test
Box.test(model1$residual,lag=6,type="Box-Pierce")##
## Box-Pierce test
##
## data: model1$residual
## X-squared = 8.4221, df = 6, p-value = 0.2088
☑️ 11차 Box-Pierce Test : 자기상관 존재
Box.test(model1$residual,lag=11,type="Box-Pierce")##
## Box-Pierce test
##
## data: model1$residual
## X-squared = 20.352, df = 11, p-value = 0.04075
☑️ 12차 Box-Pierce Test
Box.test(model1$residual,lag=12,type="Box-Pierce")##
## Box-Pierce test
##
## data: model1$residual
## X-squared = 20.458, df = 12, p-value = 0.0589
☑️ 13차 Box-Pierce Test : 자기상관
Box.test(model1$residual,lag=13,type="Box-Pierce")##
## Box-Pierce test
##
## data: model1$residual
## X-squared = 24.304, df = 13, p-value = 0.02845
위에서 11차, 13차 AR 모델의 자기상관이 남아 있으므로, seasonality를 파악해야 합니다.
따라서, 모델 피팅 시 seasonal 옵션을 사용해줍니다.
seasonal = list(c(1, 0, 0)) 옵션 : AR 모델에 12차 항을 추가함. 즉, 모델에 $Y_{t-12}$ 항을 추가해줍니다.
model2<-arima(inf,order=c(2,0,0), seasonal = list(order = c(1, 0, 0), period =12))
summary(model2)##
## Call:
## arima(x = inf, order = c(2, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12))
##
## Coefficients:
## ar1 ar2 sar1 intercept
## 0.3624 -0.1766 -0.0264 0.2559
## s.e. 0.0603 0.0624 0.0697 0.0144
##
## sigma^2 estimated as 0.03892: log likelihood = 54.84, aic = -99.67
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -8.709875e-05 0.1972847 0.1423641 -Inf Inf 0.7988682 0.01537737
하지만 모델 summary 결과 계절 항이 유의미하지는 않은 것으로 보입니다.
☑️ acf plot
Acf(model2$residual)
☑️ 2차 Ljung Test
##Box Pierce or Ljung-Box tests
Box.test(model2$residual,lag=2,type="Ljung")##
## Box-Ljung test
##
## data: model2$residual
## X-squared = 0.29489, df = 2, p-value = 0.8629
☑️ 2차 Box-Pierce Test
Box.test(model2$residual,lag=2,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 0.29078, df = 2, p-value = 0.8647
☑️ 3차 Box-Pierce Test
Box.test(model2$residual,lag=3,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 3.0014, df = 3, p-value = 0.3914
☑️ 4차 Box-Pierce Test
Box.test(model2$residual,lag=4,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 3.5437, df = 4, p-value = 0.4713
☑️ 5차 Box-Pierce Test
Box.test(model2$residual,lag=5,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 4.165, df = 5, p-value = 0.5259
☑️ 6차 Box-Pierce Test
Box.test(model2$residual,lag=6,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 8.8773, df = 6, p-value = 0.1806
☑️ 11차 Box-Pierce Test : 자기상관 존재
Box.test(model2$residual,lag=11,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 21.142, df = 11, p-value = 0.03193
☑️ 12차 Box-Pierce Test : 자기상관 존재
Box.test(model2$residual,lag=12,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 21.142, df = 12, p-value = 0.04834
☑️ 12차 Box-Pierce Test : 자기상관 존재
Box.test(model2$residual,lag=13,type="Box-Pierce")##
## Box-Pierce test
##
## data: model2$residual
## X-squared = 24.936, df = 13, p-value = 0.02353Box Test 결과 11차 ACF가 해결이 안 됐고, 그 결과도 model1과 별반 다르지 않았습니다. seasonality를 잡아도 해당 변수가 유의미하지 않았고 잡지 않은 것과 별반 차이가 없기에, AR(2) 모형인 model1을 최종 모델로 채택해였습니다. 물론 완벽하지는 않은 모델입니다. 다른 시계열 함수를 사용한다면, 조금 더 나아질 수 있을 것입니다.
📉 06. prediction : predict(model, n.ahead) 함수
n.ahead 옵션 : 데이터의 마지막 날짜부터 몇 단위까지 예측할 것인지 지정합니다.
☑️ model1 predict
predict(model1, n.ahead=6)## $pred
## Jun Jul Aug Sep Oct Nov
## 2006 0.2626406 0.2237040 0.2431482 0.2572761 0.2588472 0.2568426
##
## $se
## Jun Jul Aug Sep Oct Nov
## 2006 0.1973391 0.2098627 0.2101056 0.2107652 0.2108068 0.2108123
☑️ model2 predict
predict(model2, n.ahead=6)## $pred
## Jun Jul Aug Sep Oct Nov
## 2006 0.2687100 0.2158133 0.2359122 0.2316717 0.2585063 0.2806217
##
## $se
## Jun Jul Aug Sep Oct Nov
## 2006 0.1972847 0.2098389 0.2100287 0.2106265 0.2106678 0.2106718
각 모델의 예측 결과는 위와 같습니다. 원본 데이터 이후 자료가 있다면 모델의 성능도 파악할 수 있을 것입니다.
다음 포스팅에서는 자기회귀시차분포모형, ARDL Model에 대해 알아보도록 하겠습니다.
'📈📉 비즈니스 어낼리틱스 > 🕑 시계열 분석' 카테고리의 다른 글
| 🕑시계열 데이터 분석 13 - 자기상관 해결 5. 자기상관오차회귀모형 (0) | 2023.04.25 |
|---|---|
| 🕑시계열 데이터 분석 12 - 자기상관 해결 4. 자기회귀시차분포모형(ARDL Model) (2) | 2023.04.25 |
| 🕑시계열데이터 분석 10 - 자기상관 해결 2. 회귀가변수 (2) | 2023.04.20 |
| 🕑시계열데이터 분석 09 - 자기상관 해결 1. 차분 (2) | 2023.04.20 |
| 🕑시계열데이터 분석 08 - 시계열 회귀분석 (0) | 2023.04.15 |