이 글은
2022.12.06 - [금융공학] - 기초자산 변동성: EWMA 변동성#1
에서 다루었던 EWMA 변동성을 직접 코딩하여 살펴보고, 250일 역사적 변동성과 비교해 보도록 하겠습니다.
Python Code
import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
year_scaler = np.sqrt(252)
index_name = '^KS11'
df_index = yf.download(index_name, start="2000-01-01", end="2022-11-01")
# print(df_index)
df = df_index[['Close']].copy()
df['rate'] = (df.Close - df.Close.shift(1)) / df.Close.shift(1)
vol_days = 250
vol_colname = 'vol_' + str(vol_days)
df[vol_colname] = year_scaler * df['rate'].rolling(window=vol_days).std()
ewma_lambda = 0.99
tolerance_dim = -8
cnt_weight = int(tolerance_dim / np.log10(ewma_lambda) + 1)
print('Number of weight: {}'.format(cnt_weight))
weight = (1 - ewma_lambda) * ewma_lambda ** np.arange(cnt_weight)
weight = weight[::-1]
df['rate_shift'] = df['rate'].copy().shift(1)
df['ewma_vol'] = np.sqrt(df['rate_shift'].rolling(cnt_weight).apply(lambda x: np.sum(weight * x ** 2))) * year_scaler
df = df.dropna()
df['ewma_vol_recursive'] = 0
init_ewma = df['ewma_vol'].iloc[0]
for i in range(len(df)):
if i == 0:
df['ewma_vol_recursive'].iloc[i] = init_ewma
else:
temp_value = df['rate'].iloc[i - 1] ** 2 * (1 - ewma_lambda) * year_scaler ** 2 + \
ewma_lambda * df['ewma_vol_recursive'].iloc[i - 1] ** 2
temp_value = np.sqrt(temp_value)
df['ewma_vol_recursive'].iloc[i] = temp_value
fig, ax = plt.subplots()
ax.plot(df.index, df[vol_colname], label=vol_colname)
ax.plot(df.index, df.ewma_vol, label='ewma vol')
ax.xaxis.set_major_locator(mdates.YearLocator(1))
ax.legend()
plt.show()
분석 대상은 KOSPI 지수입니다.
df 구성하는 부분까지는 [금융공학] - 기초자산 변동성: 역사적 변동성#1의 code에서 다루었으니, 넘어가고 EWMA 관련 부분만 살펴보겠습니다.
ewma_lambda = 0.99 #λ=0.99 로 세팅
tolerance_dim = -8 #tolerance를 10^(-8)로 세팅, 즉 1억분의 1
# 가중치의 합이 1과 차이나는 것이 1억분의 1이 되는 가중치 개수를 찾기 위함
cnt_weight = int(tolerance_dim / np.log10(ewma_lambda) + 1)
print('Number of weight: {}'.format(cnt_weight))
○ 허용오차를 $\epsilon$이라 할 때, $N =\log_{10}\epsilon / \log_{10}\lambda$ 를 만족하는 $N$을 찾는 과정입니다. $\log_{10}\epsilon = -8$ 이므로 이를 tolerance_dim이라는 변수로 관리합니다. 이렇게 찾은 $N$에 1을 더하여 int() 문법을 사용하여 자연수를 찾습니다.
weight = (1 - ewma_lambda) * ewma_lambda ** np.arange(cnt_weight)
# (1-λ)λ^(i-1) 이 i번째 가중치이므로 이런 가중치를 위에서 구한 N개 설정
weight = weight[::-1]
# 지수최근일자 등락률에 가장 큰 가중치가 곱해져야 하므로, 배열을 reverse시킴
df['rate_shift'] = df['rate'].copy().shift(1)
# 오늘의 ewma는 어제,그제,.. 의 등락률이 필요하므로, 한칸씩 shift함
df['ewma_vol'] = np.sqrt(df['rate_shift'].rolling(cnt_weight).apply(lambda x: np.sum(weight * x ** 2))) * year_scaler
# rate_shit 칼럼에서 N개만큼 rolling하여 weight_i * (r_{t-i)}^2 의 합을 구한 뒤 연변동성으로 변환
○ 두 번째 줄의 weight [::-1]은 순서를 거꾸로 하는 indexing 방법입니다. 일반적으로 numpy.array 가 있을 때, [::-1]을 붙여 reversed 할 수 있습니다.
○ 마지막 줄은 저번 글에서 살펴봤던 식
$$ \sigma_t ^2 = \sqrt {252} \sum_{i=1}^\infty (1-\lambda)\lambda^{i-1} r_{t-i}^2$$
을 구현한 코드입니다. lambda 라는 python 명령어가 나오는데요, 일회성으로 쓸 수 있는 간단한 함수를 작성할 때 유효합니다.
○ 마지막 줄의 apply(func)는 명령어 뜻 그대로, func라는 함수를 주어진 dataframe의 열이나, 행에 적용하라는 뜻입니다. apply 함수안에 인자를 func만 쓰면, func라는 함수를 열(column)에 적용하라는 의미입니다. 따라서 lambda 용법과 같이 쓰이는 경우가 많죠.
df = df.dropna() #N/A가 있는 행을 모두 지운다.
df['ewma_vol_recursive'] = 0 #EWMA recursive formula 값을 저장하기 위한 field 생성
init_ewma = df['ewma_vol'].iloc[0] # dropna로 na를 지운 df의 가장 오래전 ewma vol 추출
# 이값을 가지고 재귀적으로 모든 값을 만들어 낼 예정
for i in range(len(df)):
if i == 0:
df['ewma_vol_recursive'].iloc[i] = init_ewma # 0행의 초기값 대입
else:
# \sigma_t^2 = 252 (1-\lambda)r_{t-1}^2 + \lambda \sigma_{t-1}^2 을 구현
temp_value = df['rate'].iloc[i - 1] ** 2 * (1 - ewma_lambda) * year_scaler ** 2 + \
ewma_lambda * df['ewma_vol_recursive'].iloc[i - 1] ** 2
temp_value = np.sqrt(temp_value) # 제곱근을 취하여
df['ewma_vol_recursive'].iloc[i] = temp_value #다음 행 값 삽입
○ for문 안의 수식은
$$\sigma_t^2 = 252(1-\lambda)r_{t-1}^2 + \lambda \sigma_{t-1}^2 \tag{1}$$
로 재귀적으로 구하는 과정입니다.
결과
$\lambda=0.99$ 일 때, 가중치의 합이 tolerance 1억 분의 1 이하로 들어오는 가중치 개수는
Number of weight: 1833
입니다. 아주 많죠. 이를 이용하여 EWMA와 역사적 변동성 250일짜리를 같이 그려보면,
처럼 그려집니다.
250일 역사적 변동성과 마찬가지로 시장에 crash가 있었던 2008년 리먼사태, 2011년 미국 신용등급 강등, 2020년 초 COVID19에서 피크(peak)가 생깁니다.
따라서 시장 상황을 충실히 반영하고 있다고 볼 수 있는데요, 역사적 변동성의 최대 단점인 "미래의 변동성이 아닌 과거 움직임에서 추출한 값" 이라는 특징을 EWMA 변동성 역시 가지고 있습니다.
또한 $\lambda$을 어떻게 설정하는지에 따라 변동성의 움직임이 극명히 차이가 나는데요. 아래의 분석을 보실까요.
EWMA의 $\lambda$를 다르게 해보자.
import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
year_scaler = np.sqrt(252)
index_name = '^KS11'
df_index = yf.download(index_name, start="2000-01-01", end="2022-11-01")
df = df_index[['Close']].copy()
df['rate'] = (df.Close - df.Close.shift(1)) / df.Close.shift(1)
ewma_lambda = [0.99, 0.96, 0.94]
tolerance_dim = -8
df['rate_shift'] = df['rate'].copy().shift(1)
fig, ax = plt.subplots()
color_list = ['black', 'red','gray']
for i, elam in enumerate(ewma_lambda):
cnt_weight = int(tolerance_dim / np.log10(elam) + 1)
print('Number of weight(λ={}): {}'.format(elam, cnt_weight))
weight = (1 - elam) * elam ** np.arange(cnt_weight)
weight = weight[::-1]
col_name = 'ewma_vol(λ={})'.format(elam)
df[col_name] = np.sqrt(df['rate_shift'].rolling(cnt_weight).apply(lambda x: np.sum(weight * x ** 2))) * year_scaler
ax.plot(df.index, df[col_name], label=col_name, color=color_list[i])
ax.xaxis.set_major_locator(mdates.YearLocator(1))
ax.legend()
plt.show()
코드 자체는 위의 코드와 거의 비슷합니다. 결과를 보시면,
Number of weight(λ=0.99): 1833
Number of weight(λ=0.96): 452
Number of weight(λ=0.94): 298
당연히 $\lambda$값이 작아질수록, decay 속도가 빨라 보다 적은 개수의 가중치가 필요할 테고요, 변동성 그래프를 보면,
재귀적 함수 식인 식(1)을 보면, $\lambda$가 클수록 오늘의 변동성 $\sigma_t$는 어제의 변동성 $\sigma_{t-1}$에 훨씬 더 영향을 받는다는 것을 알 수 있습니다. 즉, 어제와 큰 차이가 없게 되므로 그래프 자체가 좀 완만하게 나오게 됩니다.(검은색 그래프). 반면 $\lambda$가 작으면, 최근 등락률 $r_{t-1}$이 말을 하게 되므로 이 값에 영향을 좀 더 받게 되겠죠. 따라서 피크가 더 심하게 생길 수 있는 것입니다. 특히 $\lambda =0.94$ 를 쓰면 EWMA 변동성이 80%을 넘는 시기가 존재하기도 하죠.
EWMA는 이렇게 들쑥날쑥한 변동성 시계열을 제공해주므로, 트레이더의 입장에서는 많은 검토를 해 보고 써야 할 것입니다. 현재에도 증권사의 Equity 헤지 하우스나 채권평가사 들도 이런 EWMA 변동성을 구해 금융상품을 평가하고 있습니다.
댓글