[시계열분석] 잔차진단 실습(Python) - 잔차진단 시각화 및 분석(bike-sharing-demand dataset)

[시계열분석] 잔차진단(1) - 백색잡음, 자기상관함수 : ysyblog.tistory.com/213

[시계열분석] 잔차진단(2) - 잔차진단 방향(정상성/정규분포/자기상관/등분산성) : ysyblog.tistory.com/214

해당 포스팅은 위 포스팅에 이어 진행됩니다.

데이터 코딩은 아래 포스팅에 이어 진행됩니다.

[시계열분석] 시계열 변수 추출 실습(Python)(1) - 시계열 분해 (bike-sharing-demand dataset) :ysyblog.tistory.com/209
[시계열분석] 시계열 변수 추출 실습(Python)(2) - 이동평균/지연값/증감폭/그룹화 (bike-sharing-demand dataset) :ysyblog.tistory.com/210
[시계열분석] 시계열 변수 추출 실습(Python)(3) - 종속변수들과 독립변수들과의 관계를 파악하기 위한 시각화 (bike-sharing-demand dataset) :ysyblog.tistory.com/211
[시계열분석] 시계열 변수 추출 실습(Python)(4) - 시계열 데이터 준비(train/test set 분리) (bike-sharing-demand dataset) :ysyblog.tistory.com/212
[시계열분석] 기본 모델링 실습(Python) - 모델링 및 분석 성능 평가(bike-sharing-demand dataset) : ysyblog.tistory.com/215

잔차진단 시각화

White Nose 시각화

• 첫번째 데이터는 굉장히 튀는 데이터(아웃라이어)이기 때문에 이를 제외하고 시각화합니다.

# 인덱스마다 번호 붙이기
Resid_tr_reg1['RowNum'] = Resid_tr_reg1.reset_index().index
Resid_tr_reg1

sns.set(palette="muted", color_codes=True, font_scale=2)
sns.lmplot(data=Resid_tr_reg1.iloc[1:], x='RowNum', y='Error', 
           fit_reg=True, line_kws={'color': 'red'}, size=5.2, aspect=2, ci=99, sharey=True)

 

• White Noise가 아니다.(아래쪽으로 감소하는 추세가보임) -> 정상시계열이 아님

정규분포 시각화

sns.distplot(Resid_tr_reg1['Error'].iloc[1:], norm_hist='True', fit=stats.norm) 
# norm_hist='True'는 히스토그램으로 그린다.
# fit=stats.norm은 noum,정규분포그래프를 그려주므로써 정규분포와 가까운지 비교할 수 있게 한다.

 

자기상관관계 시각화

figure, axes = plt.subplots(1, 4, figsize=(30,5))
pd.plotting.lag_plot(Resid_tr_reg1['Error'].iloc[1:], lag=1, ax=axes[0]) #lag를 볼 수 있음
pd.plotting.lag_plot(Resid_tr_reg1['Error'].iloc[1:], lag=5, ax=axes[1])
pd.plotting.lag_plot(Resid_tr_reg1['Error'].iloc[1:], lag=10, ax=axes[2])
pd.plotting.lag_plot(Resid_tr_reg1['Error'].iloc[1:], lag=50, ax=axes[3])

 

• lag가 커질수록 모여있는 형태로 변한다.
• 1일때는 상관관계가 있는것처럼 보였지만 50일때는 그렇지 않는다(줄어들었다).
• 잔차는 시간이 흘러도 상관성이 일정해야하는데 그렇지 않기 때문에 잔차가 백색잡음이 아닌 것이다.

#100개의 lag를 그리고 수치로 표현
figure, axes = plt.subplots(2,1,figsize=(12,5))
# 자기상관함수(ACF)
figure = sm.graphics.tsa.plot_acf(Resid_tr_reg1['Error'].iloc[1:], lags=100, use_vlines=True, ax=axes[0]) 

 

• 0번째 값은 1인데, 그 이유는 et와 et(자기자신)과의 상관관계기 때문이다.
• 1번째 값이 위의 첫번째 그림의 선형성 수치이다.
• leg5,lag10,lag50도 위와 같다.
• 파란색 음영(use_vlines)이 있는데 파란색음영을 벗어나면 자기상관이 있다는 것이다.
• 즉, 자기상관이 있다.
• 약한상관성

# 편자기상관함수(PACF)
figure = sm.graphics.tsa.plot_pacf(Resid_tr_reg1['Error'].iloc[1:], lags=100, use_vlines=True, ax=axes[1])

 

• 시점간의 영향을 제거한 자기상관함수이다.
• 강한상관성 : 제약조건이나 상황이 acf보다 엄격하기 때문.

 

잔차진단 수치화

정상성 테스트(ADF/KPSS)

정상성을 테스트하는 방법인 ADF test의 함수는 다음과 같다.

pd.Series(sm.tsa.stattools.adfuller(Resid_tr_reg1['Error'])[0:4]

통계량, p-value, 체크한 lag수, p-value기준들(1%,5%,10%)을 보기 위해 4개만 추출한다.

Stationarity = pd.Series(sm.tsa.stattools.adfuller(Resid_tr_reg1['Error'])[0:4], index=['Test Statistics', 'p-value', 'Used Lag', 'Used Observations'])
for key, value in sm.tsa.stattools.adfuller(Resid_tr_reg1['Error'])[4].items(): #p-value값들을 순서대로 넣음(1%,5%,10%)
    Stationarity['Critical Value(%s)'%key] = value
Stationarity['Maximum Information Criteria'] = sm.tsa.stattools.adfuller(Resid_tr_reg1['Error'])[5]
Stationarity = pd.DataFrame(Stationarity, columns=['Stationarity'])
Stationarity

 

• p-value가 0보다 작으므로 정상성데이터이다.

다른 정상성 테스트읜 kpss_test의 함수는 다음과 같다.

Stationarity_kpss = stationarity_kpss_test(Y_Data, Target_name)

정규분포 테스트(shapiro)

Normality = pd.DataFrame([stats.shapiro(Resid_tr_reg1['Error'])], index=['Normality'], columns=['Test Statistics', 'p-value']).T
Normality

 

• p-value가 0보다 작으므로 정규분포가 아니라는 것이다.

자기상관 테스트(Ljung–Box test)

Autocorrelation = pd.concat([pd.DataFrame(sm.stats.diagnostic.acorr_ljungbox(Resid_tr_reg1['Error'], lags=[1,5,10,50])[0], columns=['Test Statistics']),
                             pd.DataFrame(sm.stats.diagnostic.acorr_ljungbox(Resid_tr_reg1['Error'], lags=[1,5,10,50])[1], columns=['p-value'])], axis=1).T
Autocorrelation.columns = ['Autocorr(lag1)', 'Autocorr(lag5)', 'Autocorr(lag10)', 'Autocorr(lag50)']
Autocorrelation

등분산성 테스트(Goldfeld–Quandt test)

Heteroscedasticity = pd.DataFrame([sm.stats.diagnostic.het_goldfeldquandt(Resid_tr_reg1['Error'], X_train.values, alternative='two-sided')],
                                  index=['Heteroscedasticity'], columns=['Test Statistics', 'p-value', 'Alternative']).T #array 형식으로 넣어야함
Heteroscedasticity

 

한번에 보기

    Score = pd.concat([Stationarity_adf, Stationarity_kpss, Normality, Autocorrelation, Heteroscedasticity], join='outer', axis=1)
    index_new = ['Test Statistics', 'p-value', 'Alternative', 'Used Lag', 'Used Observations',
                 'Critical Value(1%)', 'Critical Value(5%)', 'Critical Value(10%)', 'Maximum Information Criteria']
    Score.reindex(index_new)

 

• 결과 : ADF(정상성), KPSS(비정상성), 정규분포가 아니며, 자기상관이 있으며, 등분산이 아니다..
• 데이터를 다시 재구성할 필요성이 있음

 

아웃라이어라고 빼도 될까

위의 결과는 맨 앞의 아웃라이어를 빼고 했었는데, 모든 데이터를 넣고 분석해보면 다음과 같다.

• 결과 : ADF(정상성), KPSS(정상성), 정규분포가 아니며, 자기상관이 없으며, 등분산이 아니다.
• 데이터가 하나 더 들어갔을 뿐인데 분석결과는 상당히 달라졌다.
• 즉, 아웃라이어라고 해서 무조건 빼는 것이 아니라, 신중히 생각해 보아야 한다.
  
  

코드 요약

### Functionalize
### Error analysis
def stationarity_adf_test(Y_Data, Target_name): #정상성테스트 코드
    if len(Target_name) == 0:
        Stationarity_adf = pd.Series(sm.tsa.stattools.adfuller(Y_Data)[0:4],
                                     index=['Test Statistics', 'p-value', 'Used Lag', 'Used Observations'])
        for key, value in sm.tsa.stattools.adfuller(Y_Data)[4].items():
            Stationarity_adf['Critical Value(%s)'%key] = value
            Stationarity_adf['Maximum Information Criteria'] = sm.tsa.stattools.adfuller(Y_Data)[5]
            Stationarity_adf = pd.DataFrame(Stationarity_adf, columns=['Stationarity_adf'])
    else:
        Stationarity_adf = pd.Series(sm.tsa.stattools.adfuller(Y_Data[Target_name])[0:4],
                                     index=['Test Statistics', 'p-value', 'Used Lag', 'Used Observations'])
        for key, value in sm.tsa.stattools.adfuller(Y_Data[Target_name])[4].items():
            Stationarity_adf['Critical Value(%s)'%key] = value
            Stationarity_adf['Maximum Information Criteria'] = sm.tsa.stattools.adfuller(Y_Data[Target_name])[5]
            Stationarity_adf = pd.DataFrame(Stationarity_adf, columns=['Stationarity_adf'])
    return Stationarity_adf

def stationarity_kpss_test(Y_Data, Target_name): #정상성테스트 코드
    if len(Target_name) == 0:
        Stationarity_kpss = pd.Series(sm.tsa.stattools.kpss(Y_Data)[0:3],
                                      index=['Test Statistics', 'p-value', 'Used Lag'])
        for key, value in sm.tsa.stattools.kpss(Y_Data)[3].items():
            Stationarity_kpss['Critical Value(%s)'%key] = value
            Stationarity_kpss = pd.DataFrame(Stationarity_kpss, columns=['Stationarity_kpss'])
    else:
        Stationarity_kpss = pd.Series(sm.tsa.stattools.kpss(Y_Data[Target_name])[0:3],
                                      index=['Test Statistics', 'p-value', 'Used Lag'])
        for key, value in sm.tsa.stattools.kpss(Y_Data[Target_name])[3].items():
            Stationarity_kpss['Critical Value(%s)'%key] = value
            Stationarity_kpss = pd.DataFrame(Stationarity_kpss, columns=['Stationarity_kpss'])
    return Stationarity_kpss

def error_analysis(Y_Data, Target_name, X_Data, graph_on=False):
    for x in Target_name:
        Target_name = x
    X_Data = X_Data.loc[Y_Data.index]

    if graph_on == True:
        ##### Error Analysis(Plot)
        Y_Data['RowNum'] = Y_Data.reset_index().index

        # Stationarity(Trend) Analysis
        sns.set(palette="muted", color_codes=True, font_scale=2)
        sns.lmplot(x='RowNum', y=Target_name, data=Y_Data, fit_reg='True', size=5.2, aspect=2, ci=99, sharey=True)
        del Y_Data['RowNum']

        # Normal Distribution Analysis
        figure, axes = plt.subplots(figsize=(12,8))
        sns.distplot(Y_Data[Target_name], norm_hist='True', fit=stats.norm, ax=axes)

        # Lag Analysis
        length = int(len(Y_Data[Target_name])/10)
        figure, axes = plt.subplots(1, 4, figsize=(12,3))
        pd.plotting.lag_plot(Y_Data[Target_name], lag=1, ax=axes[0])
        pd.plotting.lag_plot(Y_Data[Target_name], lag=5, ax=axes[1])
        pd.plotting.lag_plot(Y_Data[Target_name], lag=10, ax=axes[2])
        pd.plotting.lag_plot(Y_Data[Target_name], lag=50, ax=axes[3])

        # Autocorrelation Analysis
        figure, axes = plt.subplots(2,1,figsize=(12,5))
        sm.tsa.graphics.plot_acf(Y_Data[Target_name], lags=100, use_vlines=True, ax=axes[0])
        sm.tsa.graphics.plot_pacf(Y_Data[Target_name], lags=100, use_vlines=True, ax=axes[1])

    ##### Error Analysis(Statistics)
    # Checking Stationarity
    # Null Hypothesis: The Time-series is non-stationalry
    Stationarity_adf = stationarity_adf_test(Y_Data, Target_name)
    Stationarity_kpss = stationarity_kpss_test(Y_Data, Target_name)

    # Checking of Normality
    # Null Hypothesis: The residuals are normally distributed
    Normality = pd.DataFrame([stats.shapiro(Y_Data[Target_name])],
                             index=['Normality'], columns=['Test Statistics', 'p-value']).T

    # Checking for Autocorrelation
    # Null Hypothesis: Autocorrelation is absent
    Autocorrelation = pd.concat([pd.DataFrame(sm.stats.diagnostic.acorr_ljungbox(Y_Data[Target_name], lags=[1,5,10,50])[0], columns=['Test Statistics']),
                                 pd.DataFrame(sm.stats.diagnostic.acorr_ljungbox(Y_Data[Target_name], lags=[1,5,10,50])[1], columns=['p-value'])], axis=1).T
    Autocorrelation.columns = ['Autocorr(lag1)', 'Autocorr(lag5)', 'Autocorr(lag10)', 'Autocorr(lag50)']

    # Checking Heteroscedasticity
    # Null Hypothesis: Error terms are homoscedastic
    Heteroscedasticity = pd.DataFrame([sm.stats.diagnostic.het_goldfeldquandt(Y_Data[Target_name], X_Data.values, alternative='two-sided')],
                                      index=['Heteroscedasticity'], columns=['Test Statistics', 'p-value', 'Alternative']).T
    Score = pd.concat([Stationarity_adf, Stationarity_kpss, Normality, Autocorrelation, Heteroscedasticity], join='outer', axis=1)
    index_new = ['Test Statistics', 'p-value', 'Alternative', 'Used Lag', 'Used Observations',
                 'Critical Value(1%)', 'Critical Value(5%)', 'Critical Value(10%)', 'Maximum Information Criteria']
    Score.reindex(index_new)
    return Score
# error_analysis(Resid_tr_reg1[1:], ['Error'], X_train, graph_on=True)

 

해당 포스팅은 패스트캠퍼스의 파이썬을 활용한 <시계열 데이터 분석 A-Z 올인원 패키지> 강의를 듣고 정리한 내용입니다