H_record

시계열 데이터 심화 강의의 마지막 단계인 실습프로젝트들 중 bitcoin을 포스팅하겠습니다.

MLP

from keras.models import Sequential, Model, load_model
from keras.layers import Input, Dense, Activation, Flatten, Dropout
from keras.layers import SimpleRNN, LSTM, GRU

###MLP
#Data Loading
location= './Data/Cry/Bitcoint.csv')
raw_all= pd.read_csv(location, index_col='Date')

#Parameters
criteria ='2020-01-01'
scaler= predprocessing.minMaxScaler()
sequence=60 
batch_size=32
epoch=10
verbose=1
dropout_ratio=0

# Feature Engineering
# train & test split
train = raw_all.loc[raw_all.index < criteria,:]
test = raw_all.loc[raw_all.inex >= criteria,:]
print('Train_size:', train.shape, 'Test_size:', test.shape)

## Scaling
train_scaled = scaler.fit_transform(train)
test_scaled = scaler.transform(test)

#X/Y split 
X_train, Y_train =[], []
for index in range(len(train_scaled)-sequence):
	X_train.append(train_scaled[index: index+sequence])
    Y_train.append(train_scaled[index + sequence])

X_test, Y_test =[], []
for index in range(len(test_scaled)-sequence):
	X_test.append(test_scaled[index: index + sequence])
    Y_test.append(test_scaled[index+sequence])
    
## Retype and Reshape
X_train, Y_train = np.array(X_train), np.array(Y_train) #리스트로 묶여 있는 걸 어레이로 변경
X_test, Y_test = np.array(X_test), np.array(Y_test)
print('X_train:', X_train.shape, 'Y_train:', Y_train.shape)
print('X_test:', X_test.shape, 'Y_test:', Y_test.shape)

X_train=X_train.reshape(X_train.shape[0], X_train.shpae[1])
X_test =X_test.reshape(X_test.shape[0], X_test.shape[1])
print('Reshaping for MLP')

#X_train, test shape is (row, sequence, feature)
print('X_train:', X_train.shape, 'Y_train:', Y_train.shape)
print('X_test:', X_test.shape, 'Y_test:', Y_test.shape)

model = Sequential()
model.add(Dense(128, input_shape=(X_train.shape[1],), activation='relu')) 
'''
input은 row를 빼고 다 넣으면된다. 왜? 미니배치(32)가 돌아가면서 끝까지 할테니까
핵심은 안에서 돌아가는 데이터가 무엇이냐가 중요하다
60개 feature가 히든레이어의 input 수
'''
#3차원으로 입력할 경우 (60,1)
# model.add(Dense(128, input_shape=(X_train.shape[1],X_train.shape[2]), activation='relu'))
'''
첫번째 outputshape는 (60,1) -> (none,60,128)로 바뀜
60개 feature가 128개 패턴으로 바뀐게 아니라 1 컬럼이 128개로 증식된것.
사실 이렇게 해도 무방하다 하지만 문제는 최종적으로 나오는 shape이다.
최종 shape는 (1035*60,1)가 나온다. Y_train 의 (1035,1)과 맞지가 않아서 
evaluate를 할수 없다
'''
model.add(Dropout(dropout_ratio)) #128*0.2(dropout) 수 만큼 '지우는것'
model.add(Dense(256, activation='relu'))
model.add(Dropout(dropout_ratio))
model.add(Dense(128, activation='relu'))
model.add(Dropout(dropout_ratio))
model.add(Dense(64, activation='relu'))
model.add(Dropout(dropout_ratio))
model.add(Dense(1))

model.complie(doptimaizer='adam', loss='mean_squared_error')
model.summary()

model_fit() =mmodel.fit(X_train, Y_train, batch_size=batch_size, epochs=epoch,verbose=verbose)

plt.plot(pd.DataFrame(model_fit.history())
plt.grid(True)
plt.show()

#prediction
Y_train_pred = model.predict(X_train)
Y_test_pred = model.predict(x_test)

#evaluation
result = model.evaluate(X_test, Y_test_pred)
if scaler !=[]:
	Y_train =scaler.inverse_transform(Y_train)
    Y_train_pred = scaler.inverse_transform(Y_train_pred)
    Y_test - scaler.inverse_tainsform(Y_test)
    Y_test_pred - scaler.inverse_transform(y_test_pred)
Score_MLP, Residual_tr, Residual_te = evaluation_trte(pd.DataFrame(Y_train), Y_train_pred.flatten(),
													  pd.DataFrame(Y_test), Y_test_pred.flatten(),
                                                      graph_on=True)
display(Score_MLP)                      

위 예시에서는 input_shape= (X_train.shape[1],) 였는데 이렇게 실행하기 위해서는 사전에 3차원을 2차원으로 축소해야 한다. 

그러면 3차원을 그대로 사용할려면 어떻게 해야 할까?

input_shape= (X_train.shape[1],X_train.shape[2]) 로 변경 가능하다. 즉 (60,1)로 넣을 수가 있다.

하지만 이렇게 진행할 경우 최종 shape가 Y_train shape와 맞지 않아서 Y값끼리 비교가 불가능하다. evaluation을 진행하지 못한다

 input_shape= (X_train.shape[1],) 왼쪽처럼 넣게 되면 feature의 수가 60개이다. 이 값이 히든레이어의 수만큼 변형이 되면서 특징을 뽑을 수가 있는데 

 input_shape= (X_train.shape[1],X_train.shape[2])로 넣게 되면 (60,1)로 입력이 되고 Dense(128)를 만나게 되면 feature를 60개로 인식을 안 하고 1로 인식한다. 즉 1->128로 늘어난다. 이 방법이 틀린 게 아니다. 아래 이미지로 추가 설명을 한다면 

dense_49 => (none, 60, 64)이다. 이를 Dense(1) 통과시키면 (none,60,1)이 된다. 대체 어떻게 변한 걸까?

행이 None으로 지정되는 이유는 데이터의 갯수는 계속해서 추가될 수 있기 때문에 딥러닝 모델에서는 주로 행을 무시하는데 이해를 돋기위해서 여기서의 none를 1035라고 하겠습니다(1035는 전체데이터)

None은 총 row인 1035이다. 즉 60*1 이 1035개가 있다는 뜻이 된다. 즉 각 row마다 1035개가 있다.

결국 (60*1035,1)이라는 소리인데 이 모양은 Y_train (1035,1)와 비교는 할 수 없음을 딱 봐도 알 수 있다. 

따라서 문제 해결을 위해서는 60*1035를 1035로 만들어야 한다. 아래 그림을 보면서 설명해본다

dense_49의 형태를 직접 그려 봤다. (60,64) 매트릭스가 1035개가 있는 형태이다. (위 그림의 1034 오타 -> 1035 ), 60*64사이즈의 A4용지가 1035장 쌓여 있는 것이다. 60*64=3840이다. 이 둘을 곱해서 가로로 나열하여 1035를 쌓아버린다. 이때 사용되는 것이 flatten()이다. 최종적으로 세로가 1035*3840의 매트릭스가 되고 이를 Dense(1)에 넣게 되면 3840+1(bias)가 1로 변한다

아래 Param을 보면 3841로 출력되어 있다. (3841 * 1 = 3841)

최종적으로 y_train과 비교가 가능한 (1035,1) 형태가 되었다.

 RNN

# Data Loading
location = './Data/Cryptocurrency/Bitcoin.csv'
raw_all = pd.read_csv(location, index_col='Date')
raw_all.index = pd.to_datetime(raw_all.index)

# Parameters
criteria = '2020-01-01'
scaler = preprocessing.MinMaxScaler()
sequence = 60
batch_size = 32
epoch = 10
verbose = 1
dropout_ratio = 0

# Feature Engineering
## Train & Test Split
train = raw_all.loc[raw_all.index < criteria,:]
test = raw_all.loc[raw_all.index >= criteria,:]
print('Train_size:', train.shape, 'Test_size:', test.shape)

## Scaling
train_scaled = scaler.fit_transform(train)
test_scaled = scaler.transform(test)

## X / Y Split
X_train, Y_train = [], []
for index in range(len(train_scaled) - sequence):
    X_train.append(train_scaled[index: index + sequence])
    Y_train.append(train_scaled[index + sequence])
X_test, Y_test = [], []
for index in range(len(test_scaled) - sequence):
    X_test.append(test_scaled[index: index + sequence])
    Y_test.append(test_scaled[index + sequence])

## Retype and Reshape
X_train, Y_train = np.array(X_train), np.array(Y_train)
X_test, Y_test = np.array(X_test), np.array(Y_test)
print('X_train:', X_train.shape, 'Y_train:', Y_train.shape)
print('X_test:', X_test.shape, 'Y_test:', Y_test.shape)

# RNN
model = Sequential()
#return_sequences: 출력된 아웃풋을 그대로 받아 올건지(이전 히든스테이트가 학습에 사용되도록) 
#X_train.shape[1]: 시간 정보 와 X_train.shape[2]: 컬럼을 입력으로 사용
model.add(SimpleRNN(128, input_shape=(X_train.shape[1], X_train.shape[2]), return_sequences=True, activation='relu'))
model.add(Dropout(dropout_ratio)) 
model.add(SimpleRNN(256, return_sequences=True, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(SimpleRNN(128, return_sequences=True, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(SimpleRNN(64, return_sequences=True, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(Flatten())
model.add(Dense(1))
model.compile(optimizer='adam', loss='mean_squared_error')
model.summary()
'''
컬럼이 1 -> 128개로 증식 x
60개 시퀀스를 128개 시퀀스 패턴으로 바뀐다.
즉, 60*1 이라는 패턴을 60*128개로 바꾸는 느낌으로 해석하는게 맞다
시퀀스 자체를 늘린다. 
'''
model_fit = model.fit(X_train, Y_train, 
                      batch_size=batch_size, epochs=epoch,
                      verbose=verbose)

plt.plot(pd.DataFrame(model_fit.history))
plt.grid(True)
plt.show()

# prediction
Y_train_pred = model.predict(X_train)
Y_test_pred = model.predict(X_test)

# evaluation
result = model.evaluate(X_test, Y_test_pred)
if scaler != []:
    Y_train = scaler.inverse_transform(Y_train)
    Y_train_pred = scaler.inverse_transform(Y_train_pred)
    Y_test = scaler.inverse_transform(Y_test)
    Y_test_pred = scaler.inverse_transform(Y_test_pred)
Score_RNN, Residual_tr, Residual_te = evaluation_trte(pd.DataFrame(Y_train), Y_train_pred.flatten(), 
                                                      pd.DataFrame(Y_test), Y_test_pred.flatten(), graph_on=True)
display(Score_RNN)

# error analysis
# error_analysis(Residual_te, ['Error'], pd.DataFrame(X_train.reshape(X_train.shape[0], X_train.shape[1])), graph_on=True)

LSTM

# Data Loading
location = './Data/Cryptocurrency/Bitcoin.csv'
raw_all = pd.read_csv(location, index_col='Date')
raw_all.index = pd.to_datetime(raw_all.index)

# Parameters
criteria = '2020-01-01'
scaler = preprocessing.MinMaxScaler()
sequence = 60
batch_size = 32
epoch = 10
verbose = 1
dropout_ratio = 0

# Feature Engineering
## Train & Test Split
train = raw_all.loc[raw_all.index < criteria,:]
test = raw_all.loc[raw_all.index >= criteria,:]
print('Train_size:', train.shape, 'Test_size:', test.shape)

## Scaling
train_scaled = scaler.fit_transform(train)
test_scaled = scaler.transform(test)

## X / Y Split
X_train, Y_train = [], []
for index in range(len(train_scaled) - sequence):
    X_train.append(train_scaled[index: index + sequence])
    Y_train.append(train_scaled[index + sequence])
X_test, Y_test = [], []
for index in range(len(test_scaled) - sequence):
    X_test.append(test_scaled[index: index + sequence])
    Y_test.append(test_scaled[index + sequence])

## Retype and Reshape
X_train, Y_train = np.array(X_train), np.array(Y_train)
X_test, Y_test = np.array(X_test), np.array(Y_test)
print('X_train:', X_train.shape, 'Y_train:', Y_train.shape)
print('X_test:', X_test.shape, 'Y_test:', Y_test.shape)

# LSTM
model = Sequential()
model.add(LSTM(128, input_shape=(X_train.shape[1], X_train.shape[2]), return_sequences=True, activation='relu'))
model.add(Dropout(dropout_ratio)) 
model.add(LSTM(256, return_sequences=True, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(LSTM(128, return_sequences=True, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(LSTM(64, return_sequences=False, activation="relu"))
'''
마지막 부분에서 return을 False로 놓으면 시퀀스 입력이 중지가 되고
3차원이 시간정보가 사라지면서 2차원으로 줄어든다 즉, flatten기능을 대신 할수있다.
'''
model.add(Dropout(dropout_ratio)) 
model.add(Dense(1))
model.compile(optimizer='adam', loss='mean_squared_error')
model.summary()
model_fit = model.fit(X_train, Y_train, 
                      batch_size=batch_size, epochs=epoch,
                      verbose=verbose)

plt.plot(pd.DataFrame(model_fit.history))
plt.grid(True)
plt.show()

# prediction
Y_train_pred = model.predict(X_train)
Y_test_pred = model.predict(X_test)

# evaluation
result = model.evaluate(X_test, Y_test_pred)
if scaler != []:
    Y_train = scaler.inverse_transform(Y_train)
    Y_train_pred = scaler.inverse_transform(Y_train_pred)
    Y_test = scaler.inverse_transform(Y_test)
    Y_test_pred = scaler.inverse_transform(Y_test_pred)
Score_LSTM, Residual_tr, Residual_te = evaluation_trte(pd.DataFrame(Y_train), Y_train_pred.flatten(), 
                                                      pd.DataFrame(Y_test), Y_test_pred.flatten(), graph_on=True)
display(Score_LSTM)

# error analysis
# error_analysis(Residual_te, ['Error'], pd.DataFrame(X_train.reshape(X_train.shape[0], X_train.shape[1])), graph_on=True)

 GRU

# Data Loading
location = './Data/Cryptocurrency/Bitcoin.csv'
raw_all = pd.read_csv(location, index_col='Date')
raw_all.index = pd.to_datetime(raw_all.index)

# Parameters
criteria = '2020-01-01'
scaler = preprocessing.MinMaxScaler()
sequence = 60
batch_size = 32
epoch = 10
verbose = 1
dropout_ratio = 0

# Feature Engineering
## Train & Test Split
train = raw_all.loc[raw_all.index < criteria,:]
test = raw_all.loc[raw_all.index >= criteria,:]
print('Train_size:', train.shape, 'Test_size:', test.shape)

## Scaling
train_scaled = scaler.fit_transform(train)
test_scaled = scaler.transform(test)

## X / Y Split
X_train, Y_train = [], []
for index in range(len(train_scaled) - sequence):
    X_train.append(train_scaled[index: index + sequence])
    Y_train.append(train_scaled[index + sequence])
X_test, Y_test = [], []
for index in range(len(test_scaled) - sequence):
    X_test.append(test_scaled[index: index + sequence])
    Y_test.append(test_scaled[index + sequence])

## Retype and Reshape
X_train, Y_train = np.array(X_train), np.array(Y_train)
X_test, Y_test = np.array(X_test), np.array(Y_test)
print('X_train:', X_train.shape, 'Y_train:', Y_train.shape)
print('X_test:', X_test.shape, 'Y_test:', Y_test.shape)

# GRU
model = Sequential()
model.add(GRU(128, input_shape=(X_train.shape[1], X_train.shape[2]), return_sequences=True, activation='relu'))
model.add(Dropout(dropout_ratio)) 
model.add(GRU(256, return_sequences=True, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(GRU(128, return_sequences=True, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(GRU(64, return_sequences=False, activation="relu"))
model.add(Dropout(dropout_ratio)) 
model.add(Dense(1))
model.compile(optimizer='adam', loss='mean_squared_error')
model.summary()
model_fit = model.fit(X_train, Y_train, 
                      batch_size=batch_size, epochs=epoch,
                      verbose=verbose)

plt.plot(pd.DataFrame(model_fit.history))
plt.grid(True)
plt.show()

# prediction
Y_train_pred = model.predict(X_train)
Y_test_pred = model.predict(X_test)

# evaluation
result = model.evaluate(X_test, Y_test_pred)
if scaler != []:
    Y_train = scaler.inverse_transform(Y_train)
    Y_train_pred = scaler.inverse_transform(Y_train_pred)
    Y_test = scaler.inverse_transform(Y_test)
    Y_test_pred = scaler.inverse_transform(Y_test_pred)
Score_GRU, Residual_tr, Residual_te = evaluation_trte(pd.DataFrame(Y_train), Y_train_pred.flatten(), 
                                                      pd.DataFrame(Y_test), Y_test_pred.flatten(), graph_on=True)
display(Score_GRU)

# error analysis
error_analysis(Residual_te, ['Error'], pd.DataFrame(X_train.reshape(X_train.shape[0], X_train.shape[1])), graph_on=True)

*본 내용은 시계열데이터 강의 내용 중 일부분을 요약한 내용입니다

개요

오늘은 시계열 딥러닝에 대한 내용을 포스팅합니다. 

시계열 딥러닝은 구간추정값 없이 오직 점추정값을 출력하고 연산과정에서의 설명력이 떨어지므로 신뢰성이 낮다는 단점이 있습니다. 그런데도 불구하고 시계열 딥러닝을 이용하는 이유는 뭘까요? 대략 세가지의 장점이 있습니다.

1. Feature Engineering을 자동으로 구해줍니다.

히든레이어를 통해서 데이터 안에 숨겨진 패턴을 찾아낼 수 있습니다. 히든레이어의 층이 깊어 질수록 1차,2차,3차 쭉죽 늘어나게 됩니다. 그래서 3차식에서 알수 있는 특징이 있고, 4차식에서 알수 있는 특징을 찾아 낼수 있습니다.

2. 입력& 출력의 자유도 높습니다.

입력을 어떻게 줄건지, 출력은 어떤 형태로 할건지 사용자가 정할수 있는 자유도가 있습니다. 

3. 시퀀스 길이를 제한없이 다양한 패턴으로 적용할수있다.

시퀀스길이를 사용자 마음대로 조절이 가능합니다. 4일치 데이터로 미래를 예측할건지 60일치 데이터로 예측할건지 정할수가 있다는 거죠 게다가 높은 정확도도 딥러닝을 사용 할수 밖에 없는 이유 중 하나입니다.

본격적으로 시계열 딥러닝의 알고리즘에 대한 소개를 하겠습니다.

Recurrent Neural Network: RNN

입력층 > 은닉층 > 출력층으로 연결된 단방향 신경망 외에 이전 출력값이 다시 입력으로 연결되는 순환신경망

RNN은 AR을 적용한 딥러닝 버전입니다. AR은 자기 자신의 이전 값을 입력으로 사용합니다. 이와 마찬가지로 RNN도 자기 자신의 이전의 값을 가져오는데, 여기서 차이점은 AR은 데이터를 가져오고  RNN은 히든레이어 정보를 가져옵니다.

왜 차이가 나는걸까요? RNN의 입력값은 batch사이즈 만큼 들어옵니다. 따라서 매번 층을 돌릴때마다 데이터를 누적시키는건 상식적으로 생각해도 비효율적입니다. 이 batch사이즈 만큼 들어온 데이터를 대표할 만한 값이 있어야 하는데 이를 히든레이어가 합니다. 히든레이어에 데이터의 특징 정보가 있기 때문에 다음 히든레이어를 계산 시 함께 연산하면 됩니다. 위 이미지를 보면 녹색 동그라미가 이에 해당됩니다. 주황색의 f() 는 비선형으로 만들어 주는 활성함수 입니다.

역전파(Backpropagation Through Time: BPTT)

앞으로 전진하는 Feed Forward Network(FFN)과 달리 각 시점마다 Loss Function에서 추정된 Error를 이전 시점으로 전파함

공통 파라미터: 각 시점마다 출력 희망값(Y^)과 실제 출력값(Y)의 Error로 가중치(w)들을 업데이트 함

역전파를 사용해서 가중치를 업데이트 할수 있습니다. 하지만 이는 또 다른 문제를 야기 시킵니다. 저번 포스팅에서 언급했었듯이 미분은 곱셈의 연산이 있기 때문에 0.5 *0.5 = 0.25 가 되는 것처럼 점점 값이 작아져서 잔차가 0이 됩니다.

그렇기 때문에 전파가 되는 양이 점점 작아지는 문제가 생깁니다. 

LeRu를 사용을 지양하는 이유

-> RNN계열은 이전 값들이 사용되므로 한번 0으로 빠지면 이전 정보가 손실이 되기때문에 tanh를 여전히 사용한다

Long Short-Term Memory: LSTM

곱하기로 인한 정보 손실을 막고자 덧셈을 이용한 알고리즘입니다. 

C라는 기억셀을 추가하여 이를 바탕으로 외부 계층에 은닉상태 H를 출력합니다. C는 LSTM 계층 내에서만 전달되며 과거부터 t까지 모든 정보 저장됩니다.

Ft = sigmoid(새로운 X데이터* 가중치 W + 이전 히든레이어의 정보값*가중치 U)

이 Ft에다가 이전의 모~든 정보를 기억하고있는 Ct-1를 곱해 줌으로써 정보를 잊게 합니다. 

정리: RNN은 시간이 지나면 이전 값을 잊어버리는 반면, LSTM은 각 게이트를 통해 이전 입력값이 계속 저장되어 필요한 시점에 출력 및 반영됨

장기기억셀(C)는 덧셈으로 되어있기 때문에 곱셈에 의해 작아지는 전파력이 감소하지 않게 됩니다.

Gated Recurrent Unit: GRU

LSTM보다 구조를 간결하게 만들어서 빠른속도와 유사한 성능을 보입니다.

변경된 점

1. LSTM의 Gate 3개 -> GRU Gate 2개

: Reset Gate & Update Gate(LSTM의 Forget+Input과 유사)

2. Output Gate가 없는 LSTM이기에 메모리셀에 담기는 정보 양 증가

3. LSTM의 Ct와 Ht가 하나의 벡터 Ht로 통합

4. GRU가 LSTM보다 학습할 가중치가 적어짐

업데이트를 관장하는 U는 새로운 데이터의 정보를 얼만큼 반영하느냐를 정합니다.

1-U는 위 설명처럼 얼마나 잊을지 결정하는데, U가 1이면 (새로운 입력정보를 모두 사용하겠다) ht 게이트는 작동을 할 것이고 forget을 담당하는 1-U는 0 이 되므로 기능을 발휘 하지 못합니다. 즉 잊을 부분이 없다는 겁니다. U을 구성요소를 보면 이전 히든레이어값이 들어 있습니다. 이 히든레이어에 t-1정보가 있을것입니다. 따라서 과거 정보와 현재 입력으로 들어가는 정보를 모두 사용한다면 1이 될것이고, 버릴 부분이 없다는 의미로 1-U =0이 되어서 forget할 부분이 없다는 의미 입니다.

만일 U의 비율이 0.5라면 어떤 해석이 가능할까요

U에서는 sigmoid(새로운 정보+이전 정보)의 결과가 0.5 이므로 새 정보+ 이전 정보에서 절반은 반영하고 절반은 버린다는 의미로 해석될것입니다. 이때 ht는 절반만 들어간 U * reset gate + 새로운 정보 

절반이 남은 1-U는 ht-1과 곱해집니다. 만일 1-U상태로 출력이 된다면 vanishing Gradients 문제는 해결되지 않습니다.

왜냐면 다음 t+1 시점에서 또 forget이 작동하기 때입니다. 따라서 새로운 정보를 '더하기' 해줘야 하는데 이때 ht로 더해줍니다. 

*본 내용은 시계열데이터 강의 내용 중 일부분을 요약한 내용입니다

개요

오늘은 딥러닝에 대해서 포스팅 하겠습니다. 발전 개요는 생략했고 중요하다고 생각하는 것들만 정리해 보았습니다.

단층퍼셉트론(Single-Layered Perceptron: SLP)

위 이미지는 단층 구조입니다. H0와 이를 가리키는 X들만 따로 구별해서 본다면 우리가 알고 있는 회귀분석 y=f(x) 구조와 같습니다. 즉 ho=f(x) 구조 형식을 띕니다. 즉 각 노드가 하나의 회귀분석이며 딥러닝은 이러한 회귀분석을 여러개 계산하는 것입니다.  H0 -> Y으로 넘어 갈때 활성함수가 들어가게됩니다. 활성함수는 비선형으로 만들어 주는 장치입니다. 따라서 정리하면 뉴런(노드)의 수만큼 회귀분석 한 다음에 이를 결합 하고 또다시 회귀 분석을 합니다 

Y= g(H0, H1, ... , Hk-2) 이때 g는 활성함수입니다. 

히든레이어가 1개밖에 없다는 의미는 오직 선형으로만 구분 할수 있다는 것입니다. 아래 이미지를 보겠습니다

위 예시로는 한 직선으로 구분이 가능합니다.  하지만 ...

XOR 같은 문제는 한 직선으로 구분을 못하게 됩니다. 이 문제를 구분 할려면 비선형이 필요하게 되고 이를 만들기 위해서는 히든레이어를 2개 이상 쌓아야 합니다. 아래 1번 문항이 이에 해당됩니다.

2번의 의미는 가중치 업데이트가 안된다는 것입니다. 어찌저찌 해서 회귀분석 수십개를 돌린 후 Y라는 결과를 한번 딱 냈다고 했을때, 실제 Y와 예측치 Y 간의 잔차를 토대로 잔차가 줄어드는 방향으로 가중치를 조절해야 하는데, 단층에서는 이 조절할 장치가 없습니다. 초기에 정해진 가중치가 업데이트가 되지않고 유지가 되어서 최적의 값을 찾지 못합니다. 

다층퍼셉트론(Multi-Layered Perceptron: MLP)

XOR문제를 해결하기 위해서 2개 이상의 히든 레이어를 쌓은 형태입니다. 하지만 여전히 가중치 조절은 해결하지 못한 상태입니다. 

히든 레이어가 1개면 1차식, 2개면 2차식, 3개면 3차식 형태로 증가합니다. 따라서 다양한 비선형 특징을 추출할수 있게 됩니다. 이로써 XOR 문제는 해결이 가능해 졌습니다. 이제 가중치 업데이트를 어떻게 하면 될건가 문제인데 

이는 역전파를 통해서 해결이 가능합니다.

역전파(Back Propagation: BP)

실제 Y와 예측치 Y의 차이인 error를 학습 시 사용했던 파라미터의 수만큼 잘게 쪼갭니다. 그러면 위 이미지의 노란색 화살표에 그려진 e값이 분산이 됩니다. 정리를 하자면, 단층 퍼셉트론 설명할때 노드의 수만큼 회귀분석이 실행이 된다고 했는데, 이 각각의 회귀분석에게 잔차를 나눠주고 각자 할당 받은 잔차를 가장 낮게 가중치를 조절하는 것입니다. 마치 회사 대표가 이번년도 매출은 100억원이다 라고 정해 놓으면 각 부서의 차장이나 책임자들이 목표를 소분해서 가장효율적인 방법을 찾아 일하는 것과 비슷해 보입니다. 역전파를 통해 학습 업데이트를 해결 할수 있었습니다.

하지만 이는 또다른 문제를 야기 했습니다. 이때 사용되는 활성함수는 sigmoid인데, 바로 vanishing Gradients가 발생하는 것입니다. 

히든레이어가 많아지면 잔차 값이 0이 됩니다. 그러면 가장 처음부분에서 계산된 파라미터, 여러개의 회귀분석은 잔차값을 받지 못하여 업데이트를 못하게 됩니다. 역전파 과정에서 미분끼리 연산과정이 발생합니다. 예를들어 0.5, 0.5 이 둘을 곱하면 0.25로 줄어 들듯 뒤로 갈수록 잔차가 줄어들고 결국 0밖에 남지 않는다는 것이죠

이 문제는 LeRu라는 활성함수를 통해 해결할수 있습니다.

딥러닝(Deep Learning: DL)

인공신경망의 다층구조가 심화(Deep)된 네트워크 구조의 알고리즘

 (시계열)회귀분석과 (딥러닝)인공신경망의 차이

y=f(x)의 형태를 취하고 있는 반면 인공신경망은 중간에 히든레이어가 추가된 형태입니다.

Y를 구하기 위한 연관 과정을 위 식에서 확인 할 수 있습니다.  Vi *  h(x) 에서 h(x)는 단층에서 언급했었던 내용처럼 히든레이어를 회귀식으로 취했을때 형태입니다. h(x) 안에는 X들과 가중치 값들을 포함하고 있습니다.

위 내용의 이유로 현재 딥러닝을 이용한 비지니스 수익창출하기가 아직은 어렵다고 합니다. 딥러닝은 사람의 손을 많이 탑니다. 자유도가 높아서 값이 매번 바뀌고 가중치의 설명력이 부족하기때문에 정확성이 높아도 신뢰하기가 어렵습니다.

오늘 비선형 알고리즘을 복습 했습니다

워낙 복잡해서 강사가 내용만 훑어주는 정도로 넘어갔는데

그것만으로도 복잡했습니다 

정말 멘탈 털털 털렸습니다

마코프 영역은 뭔가 뜬 구름 잡는 것처럼 이해가 안되서 상당히 힘이 들었어요

지금도 이해 못하고 있습니다 ;;

통계학과가 아닌게 이렇게 비통할수가 없네요ㅠ 

다음주 부터는 딥러닝 파트로 넘어가고, 실습 프로젝트까지 해서 다음주 금요일까지 마무리 하는걸로

계획을 짜놨습니다 그 이후에는 데이콘에 시계열 데이터로 진행하는 대회가 있더라구요 

연습삼아 그 경연에 참가 할려고 생각중입니다. 

아마 5월에나 끝날것같아요 

*본 내용은 시계열데이터 강의 내용 중 일부분을 요약한 내용입니다

전 포스팅까지는 선형알고리즘을 다뤘다면 오늘 포스팅은 비선형 알고리즘에 대한 간략한 개념정도만

올리겠습니다. 비선형이 상당히 많은 내용과 복잡한 내용이 있다고 하셔서 깊은 내용까지는 강의에서 다루지 않았어요

확실히 선형보다는 복잡한 수식이 많았습니다. 저도 완전히 100% 숙지 한게 아니라서 부족한 설명이 있더라도

양해 부탁드리겠습니다

상태 공간 모형(State Space Models)

상태 공간 모형에 앞서서 위 이미지를 통해서 비선형이 뭔지를 먼저 설명드리겠습니다.

지구의 온도를 측정하기 위해서 태평양 온도 수치를 녹색이라 하겠습니다. 파란색도 태평양 온도이긴 하지만 

실제로 측정한 값이 아니고 다른 우회적인 방법을 통해서 얻어낸 근사값이라고 가정 합니다.

이때 Y는 지구의 온도, X를 근사값으로 얻어낸 파란색이라고 설정하겠습니다.

여기서 우리는 이 X가 정말 태평양 온도를 잘 대변하고 있는지를 의심해봐야 합니다. 위 예에서는 실제값과 근사값을 나눠서 설명드렸지만 실제에서는 근사값인지 실제값인지, 근사값이 실제값과 같은지 등 아무 정보가 없습니다.

여지껏 올린 포스팅에서는 X에 대해서 어떠한 의심없이 바로 X로 적용 시켰습니다. 

비선형은 '의심' 이라는 관점에서 출발 합니다. 이 X를 그냥 넘기지 말고 정확한 값으로 '추정' 해보겠다는 것입니다.

따라서 X는 정확한 데이터가 아니니까 이를 별도로 추정해 보는 행위가 비선형의 목적입니다.

기존 방정식은 Y=f(x)+잡음 입니다. 이 X를 녹색값이 되도록 추정할려면 아래와 같은 식이 있어야 합니다.

X= g(x')+잡음 

상태 전이식이 X를 추정하기 위한 부분입니다. 상태변수 라는 것은 이 X처럼 추정하는 것들을 일컫는 말입니다.

위 사실을 통해 알수 있는 점은 X를 추정할수 있다는 개념이고 이 변수 이름을 상태변수라고 일컫는 것입니다.

지수평활법(Simple Exponential Smoothing)

정리하면 추세나 계절성 패턴이 없는 데이터를 예측할때 사용가능합니다. 

왜일까요? 생각을 해봤는데 이유를 아래 식과 같이 설명해보겠습니다.

제 나름 이유를 말씀드리자면 1번 식 or 2번 식 모두 아주 먼 과거 데이터에도 똑같은 가중치를 매긴다는 점에서

추세와 계절성이 성립이 안되는 게 아닐까 합니다.

추세라면 증감의 경향을 나타내는데 이는 현재 데이터와 과거로 갈수록 데이터간의 차이가 나게됩니다. 

예를 들어 증가하는 경향의 데이터라면 현 시점의 데이터는 값이 100이고 맨 처음 데이터 값은 5라고 했을때 차이가 많이 납니다. 다음 시점을 예측하기 위해서 5가 나온 시점과 100이 나온 시점을 같은 가중치로 계산하겠다는 것입니다.

상식적으로 납득이 안됩니다. 가장 최근 데이터 일수록 가중치를 높게 주는게 경향 파악에 유리 할 것입니다.

그래서 위 식을 이용할 경우에는 아무 패턴이 없는 데이터에만 효과가 있는게 아닌가 생각합니다.

선형 추세 알고리즘 by Holt

추세패턴이 잡을때 사용하는 예측알고리즘입니다. 추세패턴이라면 ARIMA를 생각나게 하는데요

이 알고리즘은 Level이라는 상태변수와 b라는 상태변수로 구성되여 있습니다. 상태변수니까 각각을 구성하는 방정식이 존재합니다. 추정한다는 것입니다.

여기에서 ℓt는 시간 tt에서 시계열의 수준 추정 값, bt는 시간 t에서의 시계열의 추세(기울기) 추정 값, 0≤α≤1은 수준에 대한 매개변수, 0≤β∗≤1은 추세에 대한 매개변수를 나타냅니다.

강의에서 훑고 지나가는 정도라서 자세한 원리는 잘 모르겠지만 일단 여기서 알고 넘어 가야 하는건 추세를 반영한 상태 변수가 있다는 것 정도입니다.

추세알고리즘이 있다면 당연히 계절성 알고리즘도 있습니다.

계절 알고리즘 by Holt-Winter

선형 추세 알고리즘에 계절성을 반영한 예측 알고리즘입니다. 위에서 본 추세 알고리즘에서 계절성이 추가됩니다.

시계열 데이터는 분해가 가능한데 이때 덧셈으로 분해가 되는지 곱셈으로 분해가 되는지에 따라 addive, Multiplicative 로 구별되였습니다. 이와 마찬가지로 각 경우에 따라 분해 방식이 조금씩 달라지는데 자세한 것은 아래를 참고해주세요

로컬레벨 모형(Local Level Model)

랜덤워크 모형에 관측잡음이 추가된 것으로 랜덤워크 과정을 따르는 단변수 상태변수 𝜇𝑡μt를 가집니다.

맨 위에서 지구온도를 Y, 태평양 온도를 X라고 정했는데 여기서 X가 램덤워크 일 경우입니다.

따라서 Xt= Xt-1+잡음이 되고 Yt=Xt+ 잡음이 됩니다. 랜덤워크는 비 정상성이므로 차분을 해줘야하는데 

이는 위에서 설명한 것처럼 ARIMA를 연상케 합니다.

그러면 이 둘의 차이점은 뭘까? 

예시를 보면 배를 찾기 위해서 바다의 진동에 의해 이동된 잡음을 이용해서 거리를 구할때 로컬 레벨 모형에 해당됨을 이해할수가 있습니다. 

실습: 로컬레벨 모형 데이터생성 및 이해

# 데이터 생성
np.random.seed(123)
model_generator = KalmanFilter(k_endog=1, k_states=1,
                               transition=[[1]], selection=[[1]], state_cov=[[10]],
                               design=[[1]], obs_cov=[[1]])
y_gener, x_gener = model_generator.simulate(100)

# 데이터 생성
plt.figure(figsize=(10,6))
plt.plot(y_gener, 'r:', label="Observation Values")
plt.plot(x_gener, 'g-', label="State Values")
plt.legend()
plt.title("Simulation of Local Level Model ($\sigma_w^2 = 10$, $\sigma_v^2 = 1$)")
plt.show()

# 데이터 생성
np.random.seed(123)
model_generator = KalmanFilter(k_endog=1, k_states=1,
                               transition=[[1]], selection=[[1]], state_cov=[[1]],
                               design=[[1]], obs_cov=[[10]])
y_gener_target, x_gener_target = model_generator.simulate(100)

# 데이터 생성
plt.figure(figsize=(10,6))
plt.plot(y_gener_target, 'r:', label="Observation Values")
plt.plot(x_gener_target, 'g-', label="State Values")
plt.legend()
plt.title("Simulation of Local Level Model ($\sigma_w^2 = 10$, $\sigma_v^2 = 1$)")
plt.show()

# 로컬레벨모형 추정
fit = sm.tsa.UnobservedComponents(y_gener_target, level='local level').fit()
fit.summary()

이 둘의 그래프 모양이 다릅니다. state values 는 랜덤워크에 해당되는 값입니다. 빨간 실선은 실제값입니다.

실제값에서 랜덤워크가 차지하는 비중을 볼수 있습니다.

왼쪽의 녹색은 실제값과 거의 유사합니다. 즉 랜덤워크의 비중이 크다는 걸 알수 있습니다(state_cov=10)

반대로 오른쪽은 랜덤워크가 차지 하는 비중이 작습니다. Yt= Xt+ 잡음 에서  잡음의 크기가 크고 Xt는 작다는 얘기입니다. 

상태공간모형 기반의 시계열 "구조화모형"

여러개의 보이지않은 성분들로 분리하는 것이 구조화모형의 특성입니다.

아래 실습은 주기특성과 AR 특성을 각각 추출하는 내용입니다. 사이클 먼저 보도록 하겠습니다

# 구조화모형 추정
# rwalk: 랜덤워크
# damped_cycle: 안정적으로 사이클을 추정
model = sm.tsa.UnobservedComponents(unemployment_US, level='rwalk', cycle=True,
                                    stochastic_cycle=True, damped_cycle=True)
result = model.fit(method='powell')
display(result.summary())
'''
unemployment_US를 랜덤워크와 stochastic_cycle로 분해한다
'''
# 추정 시각화
fig, axes = plt.subplots(2, figsize=(10, 6))
axes[0].plot(unemployment_US.index, unemployment_US.UNRATE,
             label='Unemployment Rate')

# result.level.smoothed 추정된 값
axes[0].plot(unemployment_US.index, result.level.smoothed, label='Random Walk')
axes[0].legend(loc='upper left')
axes[0].set(title='Level/Trend Component')
axes[1].set(title='Cycle Component')
axes[1].plot(unemployment_US.index, result.cycle.smoothed, label='Cyclic')
axes[1].legend(loc='upper left')
fig.tight_layout()
plt.show()

# 잔차진단
result.plot_diagnostics(figsize=(10, 8))
plt.show()

데이터에서 싸이클 패턴만 뽑힌 시각화입니다. 밑에 있는 AR 시각화 비교하면 상당히 비슷한 모양입니다.

# 구조화모형 추정
model = sm.tsa.UnobservedComponents(unemployment_US, level='rwalk', autoregressive=4)
result = model.fit(method='powell')
display(result.summary())

# 추정 시각화
fig, axes = plt.subplots(2, figsize=(10,6))
axes[0].plot(unemployment_US.index, unemployment_US.UNRATE, label='Unemployment Rate')
axes[0].plot(unemployment_US.index, result.level.smoothed, label='Random Walk')
axes[0].legend(loc='upper left')
axes[0].set(title='Level/Trend Component')
axes[1].set(title='AR Component')
axes[1].plot(unemployment_US.index, result.autoregressive.smoothed, label='AR')
axes[1].legend(loc='upper left')
fig.tight_layout()
plt.show()

# 잔차진단
result.plot_diagnostics(figsize=(10,8))
plt.show()

확률적 변동성 모형 (Stochastic Volatility Model)

식이 복잡해지고 고차원으로 들어가면  잔차는 SARIMAX를 시용해도 WN이 되지 않습니다.
이때 이 잔차를 변동성이라고 부르고 이를 상태변수로 취급을해서 g(잔차)로 변동성을 추정한다는 개념

따라서 변동성 패턴을 잡기 위한 모형입니다.

오늘은 강의 진도보다 더 중요한 자소서를 작성했습니다.

자소서는 할때마다 늘 짜릿하고 신선하네요 너무 힘이 듭니다 ㅎㅎㅎ

공부하는 것보다 더 힘든 것같아요

오전과 오후 13시까지는 유튜버 이형님을 보면서 어떻게 써야 할지 나름 생각하고 고민하는 시간을 가졌어요

말만 들으면 뭔가 알것같고 면접관이 원할만한 내용이 뭔질 딱 촉이 오는데 

막상 적을려고 하면 꽉 막힙니다

운동 갔다 와서 쉬지도 않고 내리 적고만 있었습니다 

사실 기존에 작성된 자소서가 있었지만 개선해야될 부분이 너무 많아서 전면 재 수정을 해야 했습니다 

아무튼 .. 작업이 이제 막 끝냈고 12시 마감이였는데 5분전에 제출은 할수있었습니다. 

미리 자소서를 확인해 봤더라면 더 완성도 있게 만들었을 텐데  

많이 아쉽네요 .. 점 찍어둔 기업이라서 신중히 해야 했었는데 ㅠ

아무튼 오늘은 이만 자고 내일 다시 포스팅 하겠습니다

*본 내용은 시계열데이터 강의 내용 중 일부분을 요약한 내용입니다

공적분 모형(Cointegration Model)

공적분이란

두 비정상성 시계열을 선형조합(Linear Combination)하여 생성한 시계열의 적분(Integration) 차수가 낮아지거나 정상상태가 되는 경우입니다.

위 정의가 어떤 의미인가를 생각해 보자

일단 두 비정상성 시계열 데이터가 필요해 보이고 이 둘 사이에서 뭔가 ''공통''인 부분이 있습니다. 결론부터 말씀드리자면 이 공통이어야 하는 것은 방향입니다. 보통 술에 취한 사람과 개로 예를 듭니다.

술 취한 사람은 비틀비틀 걸어 다녀서 어디로 튈지 모릅니다. 가다가 갑자기 넘어질수도 있고 중심을 잃어서 갑자기 왼쪽으로 휙 방향을 틀 수도 있습니다. 만일 술 취한 사람이 두 명이라면 이 두 명은 각자 가는 방향이 다를 것입니다. 여기서 술 취한 사람이 목줄 하고 있는 개가 있다면 얘기가 달라집니다. 개도 자기 페이스대로 이리저리 왔다 갔다 하겠지만 어찌 됐든 술 취한 주인을 따라가게 되어 있습니다. 이때 애완견의 위치를 나타내는 시계열과 취한 사람의 위치를 나타내는

시계열은 공적분이라고 합니다.

공적분이 되기위한 조건은 아래와 같습니다

• 각각의 시계열들이 모두 같은 차분의 수(order of integraion)를 가진다. (조건1)
• 시계열들의 선형결합으로 만들어진 새로운 시계열은 기존의 시계열들보다 더 낮은 차분의 수를 가져야 한다.(조건 2)

추가적인 설명을 하자면

Y1(술 취한 사람 위치 시계열)과 Y2(애완견 위치 시계열)을 각각 봤을 때 비정상 시계열 데이터입니다. 정상성으로 만들기 위해서 차분이 각각 2씩 필요해야 합니다. (조건 1: 차분의 수가 같아야 한다)

또한 Y1과 Y2를 선형 결합을 했을 때 만들어지는 비정상 시계열을 차분하게 될 때 차분의 수가 2보다 미만이어야 합니다(조건 2: 기존 시계열 차분의 수 보다 작아야 한다)

만일 애초에 차분의 수가 하나씩 필요했다면 선형 결합으로 인해 차분의 수가 0이어야 한다. 즉 선형 결합을 하면 정상성 데이터가 된다는 뜻입니다.

차분의 수가 각자 1일 경우에 대해서 좀 더 자세히 말씀드리겠습니다

둘 다 비 정상성인 랜덤워크가 있을 때 각자 방정식은 위와 같습니다

이 두 개를 선형 결합한 대문자 Y의 수식 풀이를 보면 결국 잔차만 남게 될 때 이 잔차가 WN(화이트 노이즈)라면 차분이 필요가 없게 되므로 Y1, Y2는 공적분 관계입니다. 

이를 실습으로 확인해 보겠습니다

실습: 공적분의 이해

# 데이터 생성
np.random.seed(1)
Y1= np.random.randn(500).cumsum() #차분을 하면 정상성 -> diff, 누적을 하면 비정상성 -> cumsum
Y2= 0.6*Y1 + np.random.randn(500)  #랜덤워크의 60% + 추가적인 랜덤을 추가(np.random.randn(500))  => Y1과 스케일 차이가 발생
Y = 0.6 * Y1 - Y2 #선형결합 (왼쪽 식 계산결과 e만 남는다. 즉 잔차는 WN)

#시각화 
plt.figure(figsize=(10,4))
plt.plot(Y1, 'g', label ='Y1')
plt.plot(Y2, 'b', label ='Y2')
plt.legend()
plt.show()

plt.figure(figsize=(10,4))
plt.plot(Y, label='Y * 0.6*Y1 - Y2')
plt.legend()
plt.show()

#VECM 공적분 확인
coint_result = sm.tsa.coint(Y1,Y2)
pd.DataFrame([coint_result[0], coint_result[1]], index=['statistics', 'p-value'], columns=['output'])

• 공적분 계산 결과: p-value= 0

실습 2

#데이터 로딩
import pandas_datareader.dataas web
import datetime

start = datetime.datetime(2014 ,1, 1)
end = datetime.datetime(2018, 12, 31)
raw1 = web.DataReader("005930.KS", 'yahoo', start, end)
raw2 = web.DataReader("005935.KS", 'yahoo', start, end)

# 데이터 시각화
## 시계열 움직임
raw = pd.concat([raw1.Closem raw2..Close], axis=1).dropna() #종가 비교
raw.columns = ["SE", "SE_PS"]
raw.plot(figsize=(10,5))
plt.show()

#비교 움직임
plt.figure(figsize=(8,6))
plt.plot(raw.SE, raw.SE_PS, alpha = 0.7)
plt.scatter(raw.SE.values[-1:], raw.SE_PS.values[-1:], c='r', s=300)
plt.show()

#각 비정상성 차수 추론
target = raw.SE.copy()
integ_result = pd.Series(sm.stattools.adfuller(target)[0:4],
					index=['Test Statistics', 'p-value', 'Used Lag', 'Used Observations'])
Y1_integ_order = 0
if integ_result[1] > 0.1:
    Y1_integ_order = Y1_integ_order + 1
    
    
target = raw.SE_PS.copy()
integ_result = pd.Series(sm.tsa.stattools.adfuller(target)[0:4], 
                         index=['Test Statistics', 'p-value', 'Used Lag', 'Used Observations'])
Y2_integ_order = 0
if integ_result[1] > 0.1:
    Y2_integ_order = Y2_integ_order + 1
print('Y1_order: ', Y1_integ_order, 'Y2_order: ', Y2_integ_order)

#회귀분석 적합_모델링
Y=raw.SE
X=raw.SE_Ps
X=sm.add_constant(X)
fit=sm.OLS(Y,X).fit()
display(fit.summary()) #회귀분석을 통해 coint(공적분)를 확인하기 위함

# 회귀분석 시각화
import seaborn as sns
plt.figure(figsize=(10,7))
sns.regplot(x='SE_PS', y='SE', data=raw) #regplot:데이터를 점으로 나타내면서 선형성을 함께 확인한다.
plt.tight_layout()
plt.show()

#공적분 시계열 시각화
y_integ = raw.SE - fit.params[1] 
plt.figure(figsize=(10,7))
Y_integ.plot()
plt.show()

# 공적분 시계열 비정상성 테스트
target = fit.resid 
#회귀분석 잔차결과 정상성으로 판정되었으므로 삼성전자와 삼성전자우량주는 서로 유사하게 움직인다.(공적분이다)
display(pd.Series(sm.tsa.stattools.adfuller(target)[0:4],
                  index=['Test Statistics', 'p-value', 'Used Lag', 'Used Observations']))
                  
# VECM 공적분 테스트
coint_result = sm.tsa.coint(raw.SE, raw.SE_PS)
pd.DataFrame([coint_result[0], coint_result[1], index=['statistics', 'p-value'], columns=['output'])

이 두 데이터 간에 선형성이 있다는 걸 알 수 있다.(방향이 같다) 따라서 공적분의 존재 가능성을 생각해 볼 수 있습니다.

언뜻 보기에는 하향 추세를 보이는 패턴이 보이는 것 같은데 자세한 것은 아래 통계를 통해 판별해 보겠습니다.

회귀분석의 잔차를 정상성 테스트를 한 결과입니다. p-value가 0이므로 정상성임을 알 수 있습니다.

따라서 잔차가 차분이 필요 없는 WN이므로 회귀분석을 통해 공적분이라는 사실을 알수 있습니다.

VECM 공적분 테스트 결과 5% 미만이므로 공적분입니다.

공적분을 공부하다 보니까 강아지와 술 취한 사람을 예로는 부족한 것 같아서 제가 이해하기 편하게끔 생각해 봤습니다

주로 팀 게임을 할 때 구성원 간의 실력이 다 다릅니다. 분명히 못하는 사람은 여기저기서 패작을 할 테고 에이스들은 캐리를 할 것입니다. 그리고 실력은 달라도 승리를 한다 라는 방향은 같습니다. 이때 못하는 사람과 에이 스은 서로 공적분 관계입니다. 

하지만 악의적으로 패작을 하는 못된 사람들이 있습니다. 이 사람들은 팀이 패배하는 게 목적이므로 에이스와 패작하는 사람의 방향은 다르므로 서로 공적분이 아닙니다.

위 예시가 말이 되는지 모르겠지만 저는 이렇게 생각하니까 좀 더 확 와 닿았습니다.

*본 내용은 시계열데이터 강의 내용 중 일부분을 요약한 내용입니다

오늘의 포스팅 내용은 다변량 선형확률입니다. 지금껏 진행 한 것은 Y가 1개인 경우였습니다. 다변량은 Y가 2개 이상인 것을 말합니다. 매출과 광고비가 있을때 단변량에서는 매출이 Y라고 정해놓고 시작했습니다. 하지만 반대로 광고비가 종속변수로써 영향을 더 많이 받는다면 말이 달라지게 됩니다. 따라서 어떤 변수가 가장 영향을 많이 받는 종속변수 인지를 알고자 할때 필요한게 단변량입니다. 

단변량 알고리즘 중에서 대표적으로 벡터자기회귀 모형을 먼저 보겠습니다.

벡터자기회귀 모형(Vector Autoregressive Model)

교재에서는 위 내용을 VAR의 정의로 적혀 있는데 직관적으로 이해하기 힘들어서 구글링을 해봤습니다.

• 나의 과거에만 영향을 받는 모델: 자기회귀
• 나의 과거 뿐만 아니라 다른 변수의 과거에서도 영향을 받는 모델: 벡터자기회귀

다른 변수의 과거까지 수식으로 함께 넣다보니까 벡터 형식으로 표현이 가능합니다.

VAR(1)인 경우: 각 변수 한 타임 과거

Yt-1안에 A라는 변수 At-1과 B라는 변수 Bt-1이 함계 들어 있어서 이를 간단하게 표현한게 위 식입니다.

At-1,Bt-1 계수를 벡터형식으로 모아 놓은게 A1입니다. 참고로 꺽쇠가 있다는 것은 그 안에 여러개가 들어 있다는 뜻입니다.

 실습: VAR 데이터생성 및 이해

Y1t, Y2t를 묶어서 표현 한것이 Y[t]입니다.  A1 벡터의 구성요소는 0.2, 0.3, -0.6, 1.1 입니다.

적용하는 방법은 실습을 통해 확인해 봅니다.

import numpy as np
import matplotlib.pyplot as plt
import statsmodels
import statsmodels.api as sm

#차수 입력
intercept = np.array([5,3])
matrix_A = np.array([[0.2,0.3], [-0.6,1.1]])
residual_covarinace = np.array([[1,0.8], [0.8, 2]]) 

# VAR 데이터 생성
fit = statsmodels.tsa.vector_ar.var_model.VARProcess(matrix_A, intercept, residual_covariance)
#시뮬레이션 시각화1
simul_num = 100
fit.plotsim(steps=sumul_sum, sseed=123) #100개 만큼 임의로 생성 후 시각화
plt.tight_layout()
plt.show() #첫번째, 두번째 그래프

# 시뮬레이션 시각화2_fit의 계수를 기반으로 시뮬레이션
# 임의로 데이터 두개를 시각화한 사실이 중요하다
# 예를들어 하나는 한국주식시장, 다른 하나는 미국 주식시장
simul_num =100
#sigma_u = corr과 같음
simul_values = statsmodels.tsa.vector_ar.util.varsim(fit.coefs, fit.intercept, fit.sigma_u, steps=simul_num)
plt.figure(figsize=(10,5))
plt.plot(simul_values)
plt.tight_layout()
plt.show() #두개가 겹친 그래프

#ACF 시각화
fit.plot_acorr()
plt.tight_layout()
plt.show()
'''
대각: Y1(1,1)에 대한 자기상관, Y2(2,2)에 대한 자기상관 (둘다 자기상관이 있음)
반 대각: Y1,Y2 서로에 대한 corr
'''
#VAR 모형적합
fit = sm.tsa.VAR(simul_values).fit()
display(fit.summary())

# 예측 및 시각화
# endog: 각 시점 마다의 계수가 출력된다
forecast_num = 20
pred_var =fit.forecast(fit.model.endog[-1;], steps = forecast_num)
pred_var_ci = fit.forecast_interval(fit.model.endog[-1:], steps=forecast_num)
fit.plot_forecast(forecast_num)
plt.tight_layout()
plt.show()

#잔차진단
fit.plot_acorr()
plt.tight_layout()
plt.show()
'''
각각의 autocorr이 작기때문에 화이트노이즈에 가깝다는 추론 할수있다.
'''

Results for equation이 y1과 y2 두개가 있습니다. y1을 한국시장, y2를 미국 시장이라고 했을때 방정식이 두개로 나뉘어서 어떤 변수가 종속변수로써 더 적합한지를 분석하게 됩니다. (equation=방정식) 

한국시장은 0.085의 유의하지 않는게 하나 발견이 되었고 미국시장은 모두 유의하다고 나왔습니다. 따라서 미국시장이 다른 변수에 의해 많은 영향을 받고 있기 때문에 한국시장보다 미국시장이 종속변수로써 더 적합하다고 할수 있습니다.

실선+도트는 신뢰구간이고 실선만 있는 것은 점추정 값입니다. 마지막 값이후로 steps 만큼 예측을 하게 됩니다.

임펄스 응답 함수(Impulse Response Function)

VAR 모형은 여러개의 시계열 상호상관관계를 기반으로 각각의 변수가 다른 변수에 어떤 영향을 주는지 임펄스 반응 함수로 알 수 있습니다. 즉 위 실습에서 y1은 한국주식시장, y2는 미국주식시장 이라는 가정에서 y1이 변화를 했을때 y2에 미치는 영향을 표시한 것입니다. 실습코드는 위 실습 내용과 이어짐을 알려 드립니다

# 임펄스반응함수 추정
#irf: impulse responsible function
#20라는 건 20lag을 뜻함
fit.irf(forecast_num).plot()
plt.tight_layout()
plt.show()

y1이 변화할때 y2에 끼치는 영향은 마이너스를 보이고 반대로 y2가 변했을때 y1에 끼치는 영향은 플러스입니다.

만일 주식시장아니라 소비(y1), 매출(y2)라고 했을때 소비가 줄면 매출도 감소하는 효과가 있다 라고 해석할수 있습니다. 

실습2: 거시경제 VAR 모형화 및 해석

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels
import statsmodels.api as sm

#데이터 로딩
raw = sm.datasets.macrodata.load_pandas().data
dates_info =raw[['year','quarter']].astype(int).astype(str)
raw.index = pd.DateTimeIndex(sm.tsa.datatools.dates_from_str(dates_info['year']+'Q'+date_info['quarter']))
raw_use = raw.iloc[:,2:5]

#데이터 시각화
raw_use.plot(subplots=True, figsize(12,5))
plt.tught_layout()
plt.show()

raw_use.diff(1).dropna().plot(subplots=True, figsize(12,5))
plt.tight_layout()
plt.show()

# VAR 모형적합
#추세를 제거한게 정상화가 될 확률이 높기때문에 이걸로 VAR모형적합을 한다
raw_use_return = raw_use.diff(1).dropna()
fit = sm.tsa.VAR(raw_use_return).fit(maxlags=2) #p=2
display(fit.summary())
'''
realcons(소비)이 가장 유의한 변수가 많다. 6가지 변수가 realcons에 영향을 준다는 뜻이며
가장 종속성이 짙은걸 알수있다.
'''
#예측 및 시각화
forecast = 20
fit.plot_forecast(forecast_num)
plt.tight_layout()
plt.show()

#임펄스반응함수 추정
fit.irf(forecast_num).plot()
plt.tight_layout()
plt.show()

# 잔차진단
fit.plot_acorr()
plt.tight_layout()
plt.show()
'''
튀는 값이 있긴하지만 전반적으로 WN에 가깝다는 것만 알고 넘어간다다
'''

1.realcons, realinv는 특정 시차동안 꾸준하게 양수의 효과를 보인다.
2.realgdp->realcons:  gdp가 1만큼증가시 realcons가 오히려 감소한다. 즉 gdp가 증가한다고 해서 바로 소비로 이어지지 않는다 오히려 줄인다 
3. realinv->realcons: 투자가 늘면 소비도 증가한다.
4.realcons -> realgdp: 소비를 늘리면 gdp가 1이상 늘어나는 긍정적인 효과를 보인다.
그래서 gdp를 늘리려면 realinv를 늘리는 것 보다는 소비를 늘리는게 효과가 좋다

각각의 autocorr이 작기때문에 화이트노이즈에 가깝다는 추론 할수있다.

 그래인저 인과관계 모형(Granger Causality Model)

Spurious Regression: X와 Y가 관련(인과관계)이 없거나 논리적인 스토리가 없음에도 단순하게 상관성이 높다라는 공격

X,Y 둘간의 관련성이 없는데 단순히 상관성 있다는 것이 마치 인과성이 있다고 해석하는게 Spurious(거짓된, 겉으로만 그럴싸한) Regression 용어이다.

그래프상에서 아이스크림과 상어공격이 언뜻 상관성이 있어보인다. 아이스크림 매출 증가시 상어 공격 건수도 증가하므로 이 둘은 인과관계가 있다고 착각하는게  Spurious Regression이다.

단지 여름이라서 아이스크림 판매와 상어 공격이 많아진것이 올바른 인과관계이다. 

이를 더 확장해서 생각해 보면 실생화의 징크스도 해당이 된다.

• 징크스: 머리를 염색했더니 시험점수가 100점이 나오더라

머리염색을 한 후에 시험을 치뤘더니 100점이 나왔다는 징크스얘기는 머리염색이 원인이라는 인과관계 형태를 보인다. 

머리염색과 시험성적의 인과관계를 우리는 규명할수가 없다. 하지만 정말로 머리염색을 한 후 (선행) 의 시험성적은 영향을 받는지에 대한 인과관계는 확인해 볼수 있다. 다시 말하면 인과관계는 파악 못하지만 시간 흐름에 따른 인과관계를 파악하고자 할때를 Granger Causality라는 개념을 사용해서 확인할수 있다. 아래 예를들어본다.

• 추론불가한 문제: "닭이 먼저인가 달걀이 먼저인가?" (인과관계)
• 추론가능한 문제: "닭과 달걀의 생성순서 별 서로의 영향력은 어떤가?" (Granger 인과관계)

Granger Causality를 사용하기 위한 조건들이 있다. 

just use Y식을 보면 X(다른 변수)가 없고 오직 자기자신의 과거만 해당이 된다. 현재는 과거의 나에게 얼만큼 영향을 주는지를 알수있고 

use X and Y에서 베타X 가 있는데 위 이미지처럼 Y가 달걀 X가 닭이라고 했을때 달걀의 생산량은 닭의 생상량에 얼마큼 영향을 받고 있는지를 알수가 있다.

이해를 돕기위해서 아래 실습을 실행한다.

실습: 닭과 달걀의 생산량을 통한 Granger 인과관계 이해

# 닭과 달걀을 시간적 순서로 봤을때 어떤게 선행성이 있는가
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

#데이터 로딩
location = './Data/ChickenEggProduction/Chicken_Egg_Production.txt'
raw_all = pd.read_csv(location, sep='\t')
raw_all.head()

#데이터 시각화 
plt.figure(figsize=(10,4))
plt.plot(raw_all['Chicken'], 'g', label='Chicken')
plt.legend()
plt.show()

plt.figure(figsize=(10,4))
plt.plot(raw_all['Egg'], 'b', label='Egg')
plt.legend()
plt.show()

# 정상성 변환시각화
# 정상성을 가져야 실행이 되므로 diff를 적용
plt.figure(figsize=(10,4))
plt.plot(raw_all['Chicken'].diff(1), 'g', label='Diff(Chicken)')
plt.legend()
plt.show()

plt.figure(figsize=(10,4))
plt.plot(raw_all['Egg'].diff(1),'b',label='Diff(Egg)')
plt.legend()
plt.show()

# Granger Causality 테스트
print('\n[Egg -> Chicken]')
#maxlag 몇개의 파라미터를 적용할것이가(lag값)
#addcons 상수항 적용 여부
#맨 처음 x의 입력은 array 형식이여야 하므로 values
'''
알고리즘 설계상 두번째 컬럼이 첫번째 컬럼에 영향을 주는지를 분석하므로 달걀이 두번째 컬럼으로 와야한다. 
iloc[:,1:] -> 컬럼이 닭,달걀순이므로 두번째 컬럼인 달걀이 닭에 Granger Causality가 있는지를 본다
'''
granger_result1 = sm.tsa.stattools.grangercausalitytests(raw_all.diff(1).dropna().iloc[:,1:].values, maxlag=3, verbose=True)
'''
lag1~3 통계량을 봤을때 p-value가 5%보다 작으므로 X가 Y에 영향을 준다
'''
print('\n[Chicken -> Egg]')
#순서 바꾸기 닭->달걀
granger_result2= sm.tsa.stattools.grangercausalitytests(raw_all.diff(1).iloc[:,[2,1]].values, maxlag=3, verbose=True)
'''
전혀 도움이 되고 있지 않다.

즉, 달걀->닭이다. 
'''
# 의사결정
# 닭이 달걀을 낳으면 그 수는 약 3년후까지 닭의 수에 영향을 준다
# 왜 3년? lag3까지 유의하게 나왔기 때문
# 닭의 수가 많아진다고해서 달걀을 많이 낳지는 않는다
# 달걀 -> 닭 (Granger Causality) 달걀이 닭에 Granger Causality 경향을 보인다.

달걀->닭 인 경우 lag3까지 모두 유의하며, 닭-> 달걀 인 경우 lag1까지만 유의하다

(좀더 자세한 설명은 코드 맨 밑에 있습니다)