파이썬: PCA로 종합 경제 지표 만들기 (hard data index construction)

주성분분석(PCA)는 여러 차원의 데이터를 하나의 차원으로 압축하는데 종종 쓰이는 기법이다.

PCA를 통해 여러 경제 지표를 하나의 종합 인덱스로 압축하는 것도 가능하다.

방법론과 배경은 나중에 여유가 되면 작성해보려 하고, 이 글에서는 파이썬으로 일차적으로 돌려보는 걸 다뤄보려한다.

구체적으로는 아래 st.louis 논고에 나온 hard data 지표들을 모아 하나의 인덱스로 종합하는 것이다.

https://research.stlouisfed.org/publications/economic-synopses/2017/03/28/does-data-confusion-equal-forecast-confusion/

---

1.데이터 수집

-fred louis 사이트에서 데이터를 다운받는다. 데이터 목록은 다음과 같다.

DSPIC96	Real Disposable Personal Income, Billions of Chained 2012 Dollars, Monthly, Seasonally Adjusted Annual Rate
ALTSALES	Light Weight Vehicle Sales: Autos and Light Trucks, Millions of Units, Monthly, Seasonally Adjusted Annual Rate
NEWORDER	Manufacturers' New Orders: Nondefense Capital Goods Excluding Aircraft, Millions of Dollars, Monthly, Seasonally Adjusted
ANXAVS	Manufacturers' Value of Shipments: Nondefense Capital Goods Excluding Aircraft, Millions of Dollars, Monthly, Seasonally Adjusted
PAYEMS	All Employees, Total Nonfarm, Thousands of Persons, Monthly, Seasonally Adjusted
AWHI	Indexes of Aggregate Weekly Hours of Production and Nonsupervisory Employees, Total Private, Index 2002=100, Monthly, Seasonally Adjusted
UNRATE	Unemployment Rate, Percent, Monthly, Seasonally Adjusted
HOUST	New Privately-Owned Housing Units Started: Total Units, Thousands of Units, Monthly, Seasonally Adjusted Annual Rate
PRRESCONS	Total Private Construction Spending: Residential in the United States, Millions of Dollars, Monthly, Seasonally Adjusted Annual Rate
PNRESCONS	Total Private Construction Spending: Nonresidential in the United States, Millions of Dollars, Monthly, Seasonally Adjusted Annual Rate
INDPRO	Industrial Production: Total Index, Index 2017=100, Monthly, Seasonally Adjusted
BOPTEXP	Exports of Goods and Services, Balance of Payments Basis, Millions of Dollars, Monthly, Seasonally Adjusted

(*재고/판매비, value of state and local construction이 빠져있다. 결과에 크게 영향은 없겠지만 나중에 추가해보자.)

엑셀을 받아서 첫 시트에 공백 등을 제거해주고 아래와 같이 세팅한다. 파이썬 pandas read_excel로 데이터를 읽기 위한 작업이다. 

이후 구글드라이브에 업로드해서 코랩으로 돌릴 준비를 한다. 

 

2.데이터 코랩에 로드

-PCA 관련된 내용이 아니라 생략한다.

 

from google.colab import drive
drive.mount('/content/drive')

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
import os
warnings.filterwarnings('ignore')
%matplotlib inline

excel_direc1 = str('/content/drive/MyDrive/Colab Notebooks/')
excel_direc2= str('PCA ECON harddata.xlsx')
excel_direc = os.path.join(excel_direc1, excel_direc2)
print(excel_direc)
print(os.path.exists(excel_direc))

3.데이터 전처리 

1) 결측치를 빼고 날짜 열을 인덱스로 바꿔준 다음, 날짜 내림차순으로 정렬한다.

df = pd.read_excel(excel_direc, header =0, engine= 'openpyxl')
display(df)

df=df.dropna()
df.set_index('observation_date', drop= True, inplace =True )
df.sort_index(ascending = False, inplace= True )

display(df)

2) 데이터를 변환한다.

이를 위해 데이터와 변환방법을 대응시키는 딕셔너리를 형성한다. 

delta log는 로그 차분, log level은 로그 변환이다.

로그 차분은 %성장 추세를 보이는 데이터를 다룰 때 쓴다.

로그 변환은 이상치가 있어 skew가 있는 데이터의 분산을 줄이기 위해 쓴다. 

 데이터 변환 방법은 이하를  따라했다.

https://research.stlouisfed.org/publications/economic-synopses/2017/03/28/does-data-confusion-equal-forecast-confusion/

*Hard data 중에서 Light weight vehicle sales 등에만 log level을 적용한 건 그래프를 보면 알 수 있다. (다른 데이터와 달리) %성장 추세를 보이는 데이터가 아니고 이상치가 다수 있다. 

#데이터 열에 대응되는 전처리 방법 딕셔너리 형성
transform_method_dic = {
'DSPIC96': 'deltalog',
'ALTSALES': 'loglevel',
'NEWORDER':	'deltalog',
'ANXAVS'	:'deltalog',
'PAYEMS':	'deltalog',
'AWHI'	:'deltalog',
'UNRATE':'loglevel',
'HOUST'	: 'deltalog',
'PRRESCONS':	'deltalog',
'PNRESCONS':	'deltalog',
'INDPRO'	:'deltalog',
'BOPTEXP':'deltalog',
 }

display(transform_method_dic )

이후 읽은 데이터프레임 열별로 데이터 변환을 적용시킨다. 

데이터프레임 열 이름들을 columns 리스트로 바꾸고,

for문을 사용, 위에서 만든 딕셔너러 key로 사용하여  columns 리스트 원소별로 (데이터프레임 열 이름)변환 방법 value를 조회하고,

그 value에 따라 데이터프레임을 변환해준다. 

columns = list(df)
print(columns)

for i in range(0,len(columns)) :

  if  transform_method_dic[ columns[i] ] == 'loglevel' :
        df.iloc[:, i]  = np.log10( df.iloc[:, i] )

  elif  transform_method_dic[ columns[i] ] == 'deltalog' :
        df.iloc[:, i]  = np.log10( df.iloc[:, i] ) -  np.log10( df.iloc[:, i].shift(periods=-1))

  i=i+1

display(df)

3)데이터 열별로 정규화를 해준다. ( 분포 mean 0, sd =1) 

-PCA는 데이터 스케일에 크게 영향을 받아 일반적으로 정규화를 해줘야 한다.

#데이터 열별로 정규화
norm_df=(df-df.mean())/df.std()

norm_df=norm_df.dropna()

#norm_df['UNRATE'] = norm_df['UNRATE'].apply(lambda x: x*-1)

display(norm_df)

#https://stackoverflow.com/questions/26414913/normalize-columns-of-a-dataframe

(실업률에 -1을 곱해줘야 할 것 같았는데, 생각해보니 할 필요없었다.)

저렇게 단순하게 짜도 되나 싶었지만, 판다스는 자동적으로 열단위 계산을 수행해주기 때문에 괜찮다.

결측치 제거를 한 번 더 해주는건 앞서 로그 차분하는 과정에서 결측치가 다시 생겼기 때문이다. 

 

4.PCA 수행 

- scikit-learn 라이브러리에서 PCA 모듈을 제공하기 때문에 그대로 쓰면 된다.

#PCA 수행
from sklearn.decomposition import PCA

pca = PCA(n_components=3)
PC = pca.fit_transform(norm_df)

PCdf = pd.DataFrame(data=PC, columns = ['PC1','PC2','PC3'])

PCdf['PC1'] = PCdf['PC1'].apply(lambda x: x*-1)
PCdf['PC2'] = PCdf['PC2'].apply(lambda x: x*-1)
PCdf['PC3'] = PCdf['PC3'].apply(lambda x: x*-1)

PCdf.index = index_date = norm_df.index.tolist()

display(PCdf)
display(pca.explained_variance_ratio_)

loading_matrix = pca.components_.T * np.sqrt(pca.explained_variance_)
display(loading_matrix )

 

주성분(PC1)은 전체 데이터 분산의 26%를 설명한다. 기존의 레퍼런스를 보면 주성분1로 원하는 인덱스를 형성한다. 

PCA의 적재계수 부호는 임의로 정해지기에 마음대로 바꿔도 결과에 영향을 미치지 않는다.

상식에 부합하게 부호를 바꿔주는게 일반적이다. 

5.그래프화

-주성분1을 matlotlib으로 그래프로 표현한다.

#주성분1을 그래프화

import matplotlib.pyplot as plt
#import matplotlib.dates as dates
#import matplotlib.ticker as ticker

plt.figure(figsize=(15,6))
plt.plot(index_date, PCdf.iloc[:, 0], label= 'PC1 index')
#plt.plot(index_date, PCdf.iloc[:, 1], label= 'PC2 index')
#plt.plot(index_date, PCdf.iloc[:, 2], label= 'PC3 index')
plt.legend()
plt.grid(True)
plt.title('Hard data Index')

다시 원래 복제하려했던 인덱스와 눈으로 비교해보면 꽤 비슷해보인다. 축 단위가 차이가 있긴한데, 데이터의 일부 차이 때문일 수 있다. 

*추가로 시간이 지나며 데이터가 추가되면 기존의 주성분 index값도 변화한다. 왜냐하면 정규화 과정의 표본분산 및 평균이 데이터 추가에 따라 변화하기 때문이다.