1. 문제 정의(탐색)
주제는 자유주제였고, 나는 개인적으로 금융 데이터를 다루어 보고 싶었다. 하지만 팀원들의 도메인 지식 분야가 달랐고, 나 혼자만 관심있는 주제를 끌고 가기엔 참여도 문제도 걱정되고 팀원들에게 크게 도움이 될 것 같지가 않았다. 해서 모두가 같이 겪고 있는 코로나 사태와, 코로나에 영향을 크게 받았을 것 같은 산업을 선정해 코로나가 진행됨에 따라 어떻게 영향을 받았을지 분석해보고, 코로나 상황이 심각/완화됨에 따라 산업별 매출액을 예측해 본 다음, 추후 투자 결정/사업 긴축 또는 확장 결정 등에 도움을 줄 수 있는 방향으로 진행해 보기로 했다.
2. 데이터 선정
처음엔 산업이 아니라 SSNIP 기준 시장에서 큰 파이를 차지하고 있는 회사의 주가 데이터 등도 같이 고려했다. 하지만 직관적으로 생각했을 때, 가장 먼저 염두해 둔 배달음식업체의 데이터를 찾기 힘들었고, 비상장 회사도 포함되어 있어 방향을 바꿔야 했다.
따라서 직관적으로 생각했을 때 코로나 이후 긍정적인 영향을 받았을 거라 생각한 배달 대행, 택배, 해외 운수업이 포함되어 있는 운송업의 매출액/총자산 데이터를 사용하기로 했다.
이와 대조적으로 부정적인 영향을 받았을 것이라 생각한 여가/스포츠산업을 같이 분석하기로 했다. 산업별 매출액 데이터는 한국은행경제포털에서 구할 수 있었고 코로나 확진자 수의 데이터를 WHO.int/data에서 구할 수 있었다.
3. 데이터 살펴보기
3-1 산업별 데이터
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from datetime import datetime
import pmdarima as pm
from statsmodels.tsa.seasonal import seasonal_decompose
from scipy import stats
import prophetpath = "c:/betatest/"
#### 운수업
df_trans = pd.read_csv(path + 'ECOS_TABLE_20220426_223049.csv',index_col = '항목명3') # 매출액/총자산 분기별 데이터
df_trans_year = pd.read_csv(path + 'ECOS_TABLE_20220427_002815.csv',index_col = '항목명3',thousands = ',') #연도별 매출액 데이터
df_trans_year2 = pd.read_csv(path + 'ECOS_TABLE_20220427_002516.csv',index_col = '항목명3',thousands = ',') # 연도별 총자산 데이터
df_trans = df_trans.iloc[0:2,5:] * 0.01
df_trans_year = df_trans_year.iloc[0,5:]
df_trans_year2 = df_trans_year2.iloc[0,5:]
df_trans_year = pd.concat([df_trans_year,df_trans_year2], axis = 1) ; del df_trans_year2
#### 예술, 스포츠 및 여가관련 서비스업
df_leisure = pd.read_csv(path + "ECOS_TABLE_20220505_145245.csv",index_col=('항목명3'))
df_leisure_year = pd.read_csv(path + "ECOS_TABLE_20220505_145120.csv",index_col=('항목명3'), thousands = ",")
df_leisure_year2 = pd.read_csv(path + "ECOS_TABLE_20220505_144425.csv",index_col=('항목명3'), thousands = ",")
df_leisure = df_leisure.iloc[0:2,5:] * 0.01
df_leisure_year = df_leisure_year.iloc[0,5:]
df_leisure_year2 = df_leisure_year2.iloc[0,5:]
df_leisure_year = pd.concat([df_leisure_year,df_leisure_year2], axis = 1) ; del df_leisure_year2필요없는 행과 열을 삭제해 주었다. 산업별 재무 데이터는 분기별 데이터와 연도별 데이터밖에 구하지 못했다.
더 나은 데이터를 찾지 못해 그대로 분석을 진행을 하긴 했지만, 시계의 수가 적었고 코로나에 영향을 받은 시기의 데이터는 더 적어 문제가 될 수 있을 것이라 예상했다.
분기별 데이터는 변화율(%) 단위기 때문에 100으로 나누어 주었다.
df_trans.T.info()
df_leisure.T.info()
<class 'pandas.core.frame.DataFrame'>
Index: 27 entries, 2015 1 to 2021 3
Data columns (total 2 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 총자산증가율 21 non-null float64
1 매출액증가율 21 non-null float64
dtypes: float64(2)
memory usage: 1.7+ KB
<class 'pandas.core.frame.DataFrame'>
Index: 27 entries, 2015 1 to 2021 3
Data columns (total 2 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 총자산증가율 21 non-null float64
1 매출액증가율 21 non-null float64
dtypes: float64(2)
memory usage: 648.0+ bytes
데이터는 2015년도 1분기 부터 2021년도 4분기까지 존재했다. 기간 자체도 길지 않지만, 분기별 데이터라 행의 갯수가 많지 않다. 학습할 데이터가 적어 과적합 문제 등이 걱정된다.
분기별 데이터를 살펴보았을 때, 결측치가 상당히 많았다. 전체 27개의 행 중 6개의 행이 Nan 값을 가지고 있었다.
추후 위에서 불러온 연도별 데이터를 이용해 보간하여 결측치를 채워넣어야한다.
인덱스 값이 str형식이라 date형식으로 변환해 주어야 시계열 분석을 진행할 수 있다.
3-2 코로나 관련 데이터
covid = list(filter(lambda x: 'WHO-COVID' in x ,os.listdir(path)))
covid = pd.read_csv(path + covid[0])
covid = covid[covid.iloc[:,1] == 'KR']
covid.index = pd.to_datetime(covid.iloc[:,0])
covid = covid.loc[idx,'Cumulative_cases']코로나 데이터의 경우 처음엔 확진자 수를 설명변수로 적용하려 했다. 하지만 추후 분석을 진행하며 코로나 초기의 확진자 수와 그 이후의 확진자 수가 양적으로 동일한 영향을 주지 않는다는 것을 알게되었다. 가령, 초기의 확진자 수가 이후의 확진자 수 보다 상대적으로 더 크게 영향을 주었다. 이후 추가할 데이터를 고려해 보았다.
20년도 초반부터 현재까지의 데이터가 있다. date행이 곧바로 date형식으로 바꿀 수 있는 형태여서 곧바로 바꾸어 주었다.
코로나 확진자 데이터는 수의 분포가 고르지 않아 log 변환 또는 복스콕스 변환을 해줄 필요를 느꼈다.
누적 확진자 데이터이기 때문에 해당 분기내 확진자 수로 바꾸어야 해당 기간 내 영향을 고려할 수 있다.
3-3 소비자 심리지수
#소비자 심리지수
csi = list(filter(lambda x: 'DT' in x ,os.listdir(path)))
csi = pd.read_csv(path + csi[0])
csi = csi.iloc[0,28:].T
csi = pd.DataFrame(csi.rolling(window = 3).mean()[2::3])csi 지수는 소비자가 느끼는 시장의 분위기를 반영하는 지수이다. 월별 데이터를 이동평균 처리해 분기별 데이터로 바꾸어 주었다.
4. 데이터 전처리
이미 위에서 간단하게 할 수 있는 데이터 처리는 해주었다. 이 단계에서는 데이터 전처리를 포함한 피쳐 엔지니어링(스케일링)을 같이 수행해 준다.
4-1 결측치 대체
위에서 보았던 결측치를 대체해 준다. 4분기의 데이터가 모두 NaN이었고, 연도별 데이터를 이용해 4분기의 변화율을 보간해 주었다.
def fill_na(df,df_year):
for z in range(df.shape[0]) :
for i in range(df.shape[1]) :
if np.isnan(df.iloc[z,i]) == True :
df.iloc[z,i] = df_year.iloc[(i+1)//4,z] / (df_year.iloc[(i+1)//4-1,z] * (df.iloc[z,(i-3)] + 1) * (df.iloc[z,(i-2)] + 1) * (df.iloc[z,(i-1)] + 1))-1
fill_na(df = df_trans,df_year = df_trans_year)
fill_na(df = df_leisure,df_year = df_leisure_year)
df_trans = df_trans.T
df_trans.columns = ['total_asset','total_take']
df_leisure = df_leisure.T
df_leisure.columns = ['total_asset','total_take']두 번 쓸꺼라 함수를 정의해 주고 4분기의 결측치를 연도별 데이터를 이용해 보간해 주었다.
지금 코드를 보니 동적 프로그래밍 알고리즘을 먼저 공부했으면 좀 더 편하게 채워 넣었겠다 싶다.
4-2 산업별 데이터 전처리
# datetime index로 변경
idx = pd.date_range("2015-1-1", "2021-9-30",freq = 'Q-DEC')
df_trans.index = idx
df_trans.index.name = 'date'
df_leisure.index = idx
df_leisure.index.name = 'date'
# 운수업 데이터 전처리
# 증가율 데이터에서 원단위 데이터로 변경 (필요하다면) 증가율 = 올해 매출액/전년도 매출액 -1 이 차분의 과정과 유사함.
def to_won(df,df_year):
for i in range(df.shape[0]):
if i == 0:
df.iloc[i,0] = (1+df.iloc[i,0]) * df_year.iloc[0,0]
df.iloc[i,1] = (1+df.iloc[i,1]) * df_year.iloc[0,1]
else :
df.iloc[i,0] = (1+ df.iloc[i,0]) * df.iloc[i-1,0]
df.iloc[i,1] = (1+ df.iloc[i,1]) * df.iloc[i-1,1]
to_won(df_trans,df_trans_year)
to_won(df_leisure,df_leisure_year)
del df_trans_year
del df_leisure_year인덱스를 datetime 인덱스로 바꾸어 주고 변화율(%)단위에서 당해 매출액으로 데이터를 바꾸어주었다.
변화율 단위를 유지한채 사용하면 회귀모형에서 평균이 0인 노이즈로 처리가 될 가능성이 있다.
4-3 COVID 확진자 데이터 전처리, scaling
#%% covid 데이터 분기내 확진자
idx = pd.date_range("2015-1-1", "2022-3-31",freq = 'Q-DEC')
covid_q = pd.DataFrame(0 for x in range(len(covid))).set_index(idx)
for i in range(len(covid)):
if i != 0:
covid_q.iloc[i,:] = covid.iloc[i,:]- covid.iloc[i-1,:]# 로그처리
covid_q_log = np.log10(covid_q +1)
covid_log = np.log10(covid +1)코로나 확진자 수 데이터를 n분기 - (n-1)분기 하여 해당 분기 내 확진자 수 데이터로 만들었다.
코로나 확진자 수의 경우 분산이 커 잔차의 분산이 불필요하게 커질 수 있기 때문에 로그처리를 해주었다. 0값을 많이 가지는 데이터이기 때문에 +1을 해주었다.
4-4 train, test Data 나누기
#%% train test 데이터 생성
def create_test(x):
return(x.loc[:'2021-9-30'],x.loc['2021-12-31':])
# 매출액 데이터
train_trans ,test_trans = create_test(df_trans)
train_leisure,test_leisure = create_test(df_leisure)
# covid 데이터
train_covid_q_log,test_covid_q_log = create_test(covid_q_log)
train_covid_log,test_covid_log = create_test(covid_log)
# 소비자 심리지수
train_csi,test_csi = create_test(csi)
train_csi_nor,test_csi_nor = create_test(csi_nor)원래는 3 ~ 4분기의 데이터를 test데이터로 분리해 예측이 잘 되었는지를 살펴볼 생각이었지만, 여기서 train데이터를 더 줄이면 코로나 유행 이후의 데이터가 너무 적어 모형이 제대로 적합하지 않을 것이라 생각해 포기했다.
대신 21년 3분기까지의 산업별 재무 데이터를 모두 사용하여 모델을 만들고 코로나 확진자 수 데이터나 csi지수 데이터를 사용해 2분기의 변화를 예측하고 추후 데이터가 업데이트 되었을 때 확인하기로 하였다.
5. EDA
5-1 간단한 시각화
데이터의 분포와 특징을 알기 위해 간단하게 그래프를 그렸다.
운수업 총매출, 여가업 총매출
운수업 총자산, 여가업 총자산
코로나 19 분기별 확진자, csi지수이다.
총 자산의 경우 매출 대비 코로나 유행에 한차례 늦은 반응을 보이는 것으로 보여 좀 더 즉각적인 반응을 보이는 매출액 데이터를 사용하기로 했다.
코로나 유행 직후인 2020년도 초에 운수업 매출액은 타격이 있었지만, 곧이어 가파른 성장세를 보인다.
여가업은 이와 대비되게 큰 폭으로 타격을 입었고, 이후 회복하긴했지만 전과 같은 수준엔 못미치고 있다.
하지만 총 자산 측면에서는 차이를 보이는데, 코로나 유행에서 성장의 가능성이 있는 운수업은 자산을 투자했을 것이고,
거리두기 규제로 심각하게 타격을 입은 여가업은 사내 유보를 택한 것으로 보인다.
소비자 심리지수 역시 2020년도초 가장 낮은 수치를 기록했다. 코로나 초기, 팬데믹에 의한 공포과 초기 규제로 과한 위축이 생겼다.
분기별 코로나 확진자 수의 경우에도 꾸준히 증가세를 보이고 있다.
5-2 설명변수간 상관관계
covid_q_log csi
covid_q_log 1.000000 0.732675
csi 0.732675 1.000000
위 그래프와 결과값은 코로나 분기별 확진자수와 csi점수간의 산포도와 상관관계이다.
코로나 유행 이후인 2020년도 이후의 데이터들만 대상으로 공분산성이 있는것으로 확인되어, 둘 중 하나의 데이터를 선택해야한다.
6. 모델 적합
6-1 모델 결정
적합할 모델을 결정하기 위해 몇가지 모형이 논의되었다. 머신러닝 계열로는 lstm을 고려해보았지만, 다변량 시계열모형을 적용하기엔 공부가 너무 부족하고, 시계 데이터가 너무 적어 오버피팅의 문제가 우려되었다. 단순 회귀 모형은 시계열의 특성을 반영하기 힘들었고, 벡터 자기 회귀모형의 경우 설명변수와 종속변수가 서로 영향을 미치는 모형인데, 직관적으로 Y(매출액)이 코로나 확진자에 직접적인 영향을 미쳤다고 생각하지 않다 적합하지 않다고 생각했다.
따라서 외부 충격에 의한 시계열 데이터를 분석하기 용이한 arimax 모형을 선택하여 분석을 진행했다.
패키지는 pmdarima에서 제공하는 auto-arima를 사용해 AIC를 기준으로 모델을 결정하기로 했다.
설명변수로는 코로나와 더 직접적으로 연관이 있는 코로나 확진자 데이터를 사용하기로 했다.
6-2 모델 적합(운송업)
result_trans = pm.auto_arima(train_trans.iloc[:,1],
exogenous = train_covid_q,
alpha = 0.05,
seasonal = True,
trace = True,
stepwise = True
)
result_trans.summary() SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 27
Model: SARIMAX(3, 0, 0) Log Likelihood -510.503
Date: Wed, 15 Jun 2022 AIC 1033.005
Time: 19:52:50 BIC 1040.780
Sample: 03-31-2015 HQIC 1035.317
- 09-30-2021
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
intercept 3.771e+07 5.93e-09 6.36e+15 0.000 3.77e+07 3.77e+07
0 4.619e+07 2e-08 2.3e+15 0.000 4.62e+07 4.62e+07
ar.L1 1.0180 0.396 2.569 0.010 0.241 1.795
ar.L2 0.0690 0.617 0.112 0.911 -1.140 1.278
ar.L3 -0.3524 0.413 -0.854 0.393 -1.161 0.457
sigma2 1.595e+15 3.94e-16 4.05e+30 0.000 1.6e+15 1.6e+15
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 246.35
Prob(Q): 0.95 Prob(JB): 0.00
Heteroskedasticity (H): 30.83 Skew: -3.47
Prob(H) (two-sided): 0.00 Kurtosis: 16.07
===================================================================================
AIC를 기준으로 (3,0,0)모형을 선택했다.
잔차에 대한 검정으로
융박스 검정은 Q = 0.95로 유의수준 0.05에서 잔차간 자기상관이 없다는 귀무가설을 기각하지 못했다.
쟈크베라 검정은 JB = 0.00으로 유의수준 0.05하에서 잔차가 정규분포를 이룬다는 귀무가설을 기각하였다.
이분산성 검정은 H=0.0으로 유의수준 0.05에서 분산이 동일하지 않다는 귀무가설을 기각하였다.
잔차를 살펴보았을 때, 코로나 유행 직후 큰 오차가 보인다. 급격한 매출 하락을 오차로 판단했다.
이를 특이치로 보고, 잔차의 분포와 QQ plot을 보았을때 잔차가 정규분포에 크게 벗어나지 않는것 같아 그대로 예측 곡선을 그려보았다.
빨간 점선은 예측선이고, 파란선이 실제값이다.
역시 코로나 초기의 급격한 변화를 예측하지 못한다. 해결 방법으로는 코로나 이전의 안정된 데이터와 코로나 이후의 데이터를 따로 모형을 만들거나, 코로나 확진자 데이터에 polynomal이나 exp같은 비선형적 수정을 가해 코로나 초기의 변화를 예측시키는 것이다. 하지만, 과적합 문제 때문에 특이치를 남겨두고 분석을 마무리 하기로 했다.
6-3 모델 적합(여가업)
result_leisure = pm.auto_arima(train_leisure.iloc[:,1],
exogenous = train_covid_q_log,
alpha = 0.05,
seasonal = True,
trace = True,
stepwise = True)
result_leisure.summary() SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 27
Model: SARIMAX(2, 0, 0) Log Likelihood -453.120
Date: Wed, 15 Jun 2022 AIC 916.240
Time: 20:45:05 BIC 922.719
Sample: 03-31-2015 HQIC 918.166
- 09-30-2021
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
intercept 1.762e+06 7.31e+04 24.123 0.000 1.62e+06 1.91e+06
0 4.527e+06 6.51e+05 6.954 0.000 3.25e+06 5.8e+06
ar.L1 1.2631 0.349 3.621 0.000 0.579 1.947
ar.L2 -0.3666 0.341 -1.074 0.283 -1.036 0.302
sigma2 2.132e+13 0.048 4.45e+14 0.000 2.13e+13 2.13e+13
===================================================================================
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 110.05
Prob(Q): 0.89 Prob(JB): 0.00
Heteroskedasticity (H): 9.00 Skew: -2.70
Prob(H) (two-sided): 0.00 Kurtosis: 11.28
===================================================================================
AIC를 기준으로 (2,0,0)모형을 선택했다.
운수업 모형과 마찬가지로 유의수준 0.05하에서 자기상관이 있다는 귀무가설을 기각하지는 못했지만, 정규분포를 따른다는 귀무가설과 동분산성 귀무가설은 기각되었다. 코로나 초기의 특이치가 크게 영향을 주는것 같다.
잔차의 분포를 보았을때, 자기상관은 없어 보이지만 QQ plot과 밀도그래프를 보았을때 정규분포를 따른다고 하지는 어려울 것 같다.
마찬가지로 코로나 유행 초기가 설명하기 어렵다. 코로나 유행이지나고 점차 안정되면서 초기 확진자수와 비례하지 않은 매출 감소가 모형에서는 설명되지 않는다.
7. 모형 평가
그래프상으로는 y축을 1e7, 1e8로 변환했기 때문에 비슷해보이지만, 두 산업의 총 매출액간 차이가 크기 때문에 MAPE 지수를 사용했다.
def MAPE(y_test, y_pred):
return np.mean(np.abs((y_test - y_pred) / y_test)) * 100
print(f"Trans model MAPE: {MAPE(train_trans.iloc[:,1], fit_values):.3f}")
print(f"Leisure model MAPE: {MAPE(train_leisure.iloc[:,1], fit_values):.3f}")Trans model MAPE: 11.826
Leisure model MAPE: 13.326
낮다고 할 수 없는 수치이지만, 특이치를 설명할 수 있는 방법이 있었다면 이보다 훨씬 낮았을거라 생각한다.
공부를 좀 더 해서 외생변수가 선형 관계가 아닌 비선형 관계, 일반화가법등으로 모형을 재설정 할 수 있게 만들고 싶다.
prophet등의 패키지도 시도해보았는데 arimax보다 오히려 MAPE가 떨어졌다. 외부 이벤트를 시계열 변수로 반영하는게 아니고 기간 설정으로 node를 조정하는 방식이라 큰 외부 충격에는 arimax가 더 적합하다.
후기
처음으로 주제 선정부터 분석까지 주도적으로 진행한 프로젝트였다. 하면할수록 아쉬운 점이 많았던 프로젝트였다.
나 자체가 변수관리부터 코드관리가 너무 안되어 실망했다. 파이썬을 배운지 한달정도부터 진행한 프로젝트라, 나중에 실력이 좀 나아질수록 전에 작성한 코드가 너무 조악해보였다.
사실 지금도 기껏해야 세달째라 실력이 좋은 사람이 보면 여전히 누더기처럼 보일지도...
특히 충분히 사용할 수 있었던 기법들을 데이터 때문에 포기한게 많았다. LSTM도 그렇고 자유도문제 때문에 코로나 확진자 수 이외의 변수를 사용하기 힘들었다. test데이터를 따로 나누는 것도 코로나 유행 이후 데이터가 너무 부족해 모형이 충분히 적합되지 않아 포기했다. acf와 pacf 그래프도 그려보았지만 node수가 적어 유의미한 해석이 힘들었다.
만약 데이터가 충분했다면 뉴스 등 신문에서 언급한 코로나 관련 단어들의 횟수를 사용하거나, 전국 코로나 거리두기 단계를 범주형 또는 이동평균처리하여 설명변수로 사용하는 시도를 해 볼 수 있었을 것 같다.
또는 사용한 모델을 직접 수정해 일반화 가법 모형/polynomal 회귀 모델을 적용해 비선형 관계를 고려했으면 더 좋은 결과가 나왔을 것 같다.
코랩을 사용해서 진행했는데, 서로 변수명이나 코드 진행 방식에 대한 합의가 선행되지 않아 기껏 작성한 코드들도 하나로 합치려고 한참 노력했다. 내가 팀장의 입장이라 원하는 코드를 좀 더 구체적으로 말했어야 했다고 생각한다.
다음에 프로젝트를 진행한다면 훨씬 나은 결과를 얻을 수 있을 것 같다. 아쉬웠던만큼 앞으로 주의할 점이 늘어났다.
'배우는 것들 > python' 카테고리의 다른 글
| 백준 11057번 - 오르막수 (0) | 2022.06.20 |
|---|---|
| 백준 1932번 정수삼각형 (DP) (0) | 2022.06.16 |
| return과 print의 차이 (0) | 2022.06.13 |
| [ERROR] window, anconda 환경에서 prophet 설치 (9) | 2022.05.29 |
| 데이터 과학을 위한 파이썬 머신러닝 (0) | 2022.04.04 |