지난 글
2022.09.20 - [수학의 재미/행렬 이론] - n차원 상관계수행렬을 촐레스키 분해하면?
n차원 상관계수행렬을 촐레스키 분해하면?
지난 글 2022.09.19 - [수학의 재미/행렬 이론] - 촐레스키 분해 촐레스키 분해 이번 글은 2022.09.16 - [분류 전체보기] - 상관계수 행렬은 positive definite! 상관계수 행렬은 positive definite! 이번 글은 20..
sine-qua-none.tistory.com
에서 상관계수를 가지는 $n$개의 표준정규분포를 만들어 보았습니다. 복습을 간단히 해보겠습니다.
상관관계를 갖는 표준정규분포 만들기(복습)
1. 서로 독립인 표준정규분포 $n$개, $Z_1, Z_2, \cdots, Z_n$ 을 준비한다.
2. 상관계수 행렬 $\mathbf{R}$를 $\mathbf{R} = \mathbf{L}\mathbf{L}^t$로 촐레스키 분해한다.
3. 다음의 행렬을 만족하는 $X_1, X_2, \cdots, X_n$이 우리가 원하는 상관계수를 가지는 표준정규분포이다.
$$ \begin{pmatrix} X_1\\X_2\\ \vdots \\ X_n \end{pmatrix} = \mathbf{L} \begin{pmatrix} Z_1\\Z_2\\ \vdots \\ Z_n \end{pmatrix} \tag{1}$$
Python Code : 촐레스키 분해
상관관계를 갖는 표준정규분포를 만들기 위해서는 촐레스키 분해를 할 줄 알아야 합니다. Python에서는 촐레스키 분해 함수가 있습니다. 아래 예제를 볼까요?
import pprint
import scipy
import scipy.linalg
import numpy as np
CorrMat = np.array([[1, 0.8, 0.6],
[0.8, 1, 0.4],
[0.6, 0.4, 1]])
L_scipy = scipy.linalg.cholesky(CorrMat, lower=True)
U_scipy = scipy.linalg.cholesky(CorrMat, lower=False)
L_np = np.linalg.cholesky(CorrMat)
U_np = L_np.transpose()
print("Correlation Matrix:")
pprint.pprint(CorrMat)
print("L_scipy:")
pprint.pprint(L_scipy)
print("U_scipy")
pprint.pprint(U_scipy)
print("L_np:")
pprint.pprint(L_np)
print("U_np")
pprint.pprint(U_np)
Python으로 촐레스키 분해를 할 수 있는 방법은 2가지가 있습니다.
scipy.linalg.cholesky(A, lower=True)
scipy라는 수학/공학용 오픈소스 내, linalg(linear algebra) 라이브러리 내 Cholesky 함수입니다.
행렬 A를 input으로 받고, lower = True 이면 하삼각행렬을, lower = False 이면 상삼각행렬을 반환합니다.
numpy.linalg.cholesky(A)
numpy 라이브러리 내, linalg 내 Cholesky 함수입니다. 행렬 A를 받아 촐레스키 분해하여 상삼각행렬을 반환합니다.
코드를 간단하게 보겠습니다.
import pprint # array 출력을 예쁘게 하기위한 print 모듈
import scipy # scipy.linalg.cholesky 를 쓰기 위해 import
import scipy.linalg
import numpy as np # np.linalg.cholesky 를 쓰기 위해 import
CorrMat = np.array([[1, 0.8, 0.6],
[0.8, 1, 0.4],
[0.6, 0.4, 1]])
# 분석 대상이 되는 3*3 correlation matrix (확률변수 3개)
# 이 행렬이 positive semidefinite가 아니라면 Cholesky 분해는 오류가 남
L_scipy = scipy.linalg.cholesky(CorrMat, lower=True) # scipy 를 이용한 cholesky 분해(하삼각행렬 제공)
U_scipy = scipy.linalg.cholesky(CorrMat, lower=False)
L_np = np.linalg.cholesky(CorrMat) # numpy를 이용한 Cholesky 분해
U_np = L_np.transpose()
print("Correlation Matrix:")
pprint.pprint(CorrMat) # CorrMat를 아름답게 출력해줌 (pprint.pprint)
결과는 다음과 같습니다.
Correlation Matrix:
array([[1. , 0.8, 0.6],
[0.8, 1. , 0.4],
[0.6, 0.4, 1. ]])
L_scipy:
array([[ 1. , 0. , 0. ],
[ 0.8 , 0.6 , 0. ],
[ 0.6 , -0.13333333, 0.78881064]])
U_scipy
array([[ 1. , 0.8 , 0.6 ],
[ 0. , 0.6 , -0.13333333],
[ 0. , 0. , 0.78881064]])
L_np:
array([[ 1. , 0. , 0. ],
[ 0.8 , 0.6 , 0. ],
[ 0.6 , -0.13333333, 0.78881064]])
U_np
array([[ 1. , 0.8 , 0.6 ],
[ 0. , 0.6 , -0.13333333],
[ 0. , 0. , 0.78881064]])○ 기호 L은 하삼각행렬(Lower triangular matrix), U는 상삼각행렬(Upper triangular matrix)입니다. scipy와 numpy의 결과가 모두 똑같게 나옵니다.
Python Code : 상관관계가 있는 여러 개의 표준정규분포 난수 만들기
위의 예제에서 쓴
CorrMat = np.array([[1, 0.8, 0.6],
[0.8, 1, 0.4],
[0.6, 0.4, 1]])를 상관계수로 하는 3개의 확률변수를 만들어 보도록 합시다. Python Code는 아래와 같습니다.
import pprint
import scipy
import scipy.linalg
import numpy as np
import scipy.stats as stats
CorrMat = np.array([[1, 0.8, 0.6],
[0.8, 1, 0.4],
[0.6, 0.4, 1]])
L_scipy = scipy.linalg.cholesky(CorrMat, lower=True)
U_scipy = scipy.linalg.cholesky(CorrMat, lower=False)
L_np = np.linalg.cholesky(CorrMat)
U_np = L_np.transpose()
dataSet = np.random.normal(size=(10000, 3))
print('Dataset: ')
pprint.pprint(dataSet)
print('\n')
dataSet -= np.mean(dataSet, axis=0)
dataSet /= np.sqrt(np.var(dataSet, axis=0))
mean = np.mean(dataSet, axis=0)
var = np.var(dataSet, axis=0)
print('Mean of dataset:')
pprint.pprint(mean)
print('\n')
print('Variance of dataset: ')
pprint.pprint(var)
print('\n')
corr = np.corrcoef(dataSet.T)
print('correlation matrix of independent standard normal random')
pprint.pprint(corr)
corr_dataSet = dataSet.copy()
corr_dataSet[:, 0] = L_scipy[0][0] * dataSet[:, 0]
corr_dataSet[:, 1] = L_scipy[1][0] * dataSet[:, 0] + L_scipy[1][1] * dataSet[:, 1]
corr_dataSet[:, 2] = L_scipy[2][0] * dataSet[:, 0] + L_scipy[2][1] * dataSet[:, 1] + L_scipy[2][2] * dataSet[:, 2]
print('Mean of correlated dataset: ')
pprint.pprint(np.mean(corr_dataSet, axis=0))
print('\n')
print('Variance of correlated dataset: ')
pprint.pprint(np.var(corr_dataSet, axis=0))
print('\n')
print('Correlation matrix')
pprint.pprint(np.corrcoef(corr_dataSet.T))
print('\n')
print('Normal random test(1st column of dataset: ', stats.jarque_bera(corr_dataSet[:,0]))
print('Normal random test(2nd column of dataset: ', stats.jarque_bera(corr_dataSet[:,1]))
print('Normal random test(3rd column of dataset: ', stats.jarque_bera(corr_dataSet[:,2]))
간단하게 살펴보겠습니다.
dataSet = np.random.normal(size=(10000, 3))
# 변수 3개이고 각 변수당 샘플 갯수가 10000개인 데이터셋 준비○ 변수를 $X_1, X_2, X_3$이라 하고 각각이 10,000개의 샘플을 가지고 있는 상황입니다. 즉, 데이터셋은
$$ \begin{pmatrix} \mid & \mid &\mid // X_1 & X_2 & X_3\\ \mid & \mid &\mid $$
입니다.
○ numpy.random.normal 을 사용하여 표준정규분포 난수를 만듭니다.
dataSet -= np.mean(dataSet, axis=0) # 각 열의 평균을 구해 각 열에서 각각의 평균을 빼줌
dataSet /= np.sqrt(np.var(dataSet, axis=0)) # 각 열의 분산을 구해 각 열을 표준편차로 나눠줌○ 보다 정확한 Z-score을 구하는 작업입니다.
mean = np.mean(dataSet, axis=0)
var = np.var(dataSet, axis=0)
# Z- score 로 normalizing이 잘되었는지 확인
corr = np.corrcoef(dataSet.T) # correlation matrix를 구해봄
print('correlation matrix of independent standard normal random')
pprint.pprint(corr)○ numpy.corrcoef를 이용해 상관계수 행렬을 구할 수 있습니다. 주의할 점은 변수가 행의 위치, 변수의 샘플이 열의 위치로 와야 하므로 dataSet에 transpose를 취해야 합니다.
corr_dataSet = dataSet.copy() # dataSet의 복사본 만듬. 이 데이터셋을 상관계수 있는 데이터셋으로 만들 예정
# LX 를 하여 상관관계있는 데이터셋 만듬
corr_dataSet[:, 0] = L_scipy[0][0] * dataSet[:, 0]
corr_dataSet[:, 1] = L_scipy[1][0] * dataSet[:, 0] + L_scipy[1][1] * dataSet[:, 1]
corr_dataSet[:, 2] = L_scipy[2][0] * dataSet[:, 0] + L_scipy[2][1] * dataSet[:, 1] + L_scipy[2][2] * dataSet[:, 2]○ 상관관계있는 표준정규분포는
$$L\begin{pmatrix} X_1\\X_2\\X_3 \end{pmatrix}$$ 로 얻어집니다. 이 행렬 연산을 직접 한 것입니다.
print('Mean of correlated dataset: ')
pprint.pprint(np.mean(corr_dataSet, axis=0)) # corr_dataSet의 평균 구하기
print('\n')
print('Variance of correlated dataset: ')
pprint.pprint(np.var(corr_dataSet, axis=0)) # 분산구하기
print('\n')
print('Correlation matrix')
pprint.pprint(np.corrcoef(corr_dataSet.T)) # 상관계수행렬 구하기
print('\n')○ 새롭게 얻은 데이터셋 세 변수의 평균, 분산, 상관계수 행렬을 구해봅니다. 예상대로라면, 상관계수행렬은 초기에 주어진 행렬과 같게 나와야 하죠.
# Jarque_bera 검정으로 각 변수가 표준정규분포를 잘 따르는지 확인
print('Normal random test(1st column of dataset: ', stats.jarque_bera(corr_dataSet[:,0]))
print('Normal random test(2nd column of dataset: ', stats.jarque_bera(corr_dataSet[:,1]))
print('Normal random test(3rd column of dataset: ', stats.jarque_bera(corr_dataSet[:,2]))○ 마지막으로, 새롭게 구한 데이터셋의 각 열이 표준정규 분포를 따르는지 확인해 봅니다. 해당 사항은
2022.06.17 - [수학의 재미/확률분포] - 무엇이 더 정규분포인가? Jarque-Bera Test를 참고해 보시기 바랍니다.
○ jarque_bera()를 사용하기 위해 다음의 import가 필요합니다.
import scipy.stats as stats
이제 결과를 볼까요?
Dataset:
array([[-0.5261694 , -1.0196472 , 0.02098029],
[-0.00287988, 0.22709799, -0.64654555],
[ 0.39518587, 0.41786458, 0.83072939],
...,
[ 0.94961726, 0.34046586, 1.08219931],
[ 1.08533901, 2.5843262 , -0.8059873 ],
[-1.45422358, 1.2263708 , -0.02244433]])
Mean of dataset:
array([ 2.45581333e-17, 1.97397654e-17, -4.59774580e-17])
Variance of dataset:
array([1., 1., 1.])
correlation matrix of independent standard normal random
array([[1.00000000e+00, 2.93401986e-03, 2.40809137e-04],
[2.93401986e-03, 1.00000000e+00, 1.32450244e-02],
[2.40809137e-04, 1.32450244e-02, 1.00000000e+00]])
Mean of correlated dataset:
array([ 2.45581333e-17, 5.34239319e-17, -1.31450406e-17])
Variance of correlated dataset:
array([1. , 1.00281666, 0.99697242])
Correlation matrix
array([[1. , 0.80063365, 0.60070879],
[0.80063365, 1. , 0.40720862],
[0.60070879, 0.40720862, 1. ]])
Normal random test(1st column of dataset: Jarque_beraResult(statistic=0.8269825617967153, pvalue=0.6613373006224426)
Normal random test(2nd column of dataset: Jarque_beraResult(statistic=4.213901370438609, pvalue=0.12160822338450805)
Normal random test(3rd column of dataset: Jarque_beraResult(statistic=1.8961651732749543, pvalue=0.38748327723696296)○ 원래 데이터의 평균, 분산, 상관계수행렬이 정확히 나왔습니다. 원래 데이터는 3 확률변수가 서로 독립이므로 상관계수행렬은 단위행렬이 나옵니다.
○ 새롭게 얻은 데이터셋의 상관계수 행렬을 보세요. 처음에 우리가 가정했던 행렬과 차이가 거의 나지 않습니다.
○ Jarque-Bera 검정을 보면, 통계량이 모두 5.99보다 작죠( 여기 참고.) 따라서 표준정규분포라 가정할 수 있습니다.
결론
위와 같이 서로 독립인 표준정규분포에서 서로 상관관계가 있는(우리가 원하는 상관계수를 가지는) 새로운 표준정규분포 난수를 얻을 수 있습니다. 금융공학에서는 여러 개의 금융 자산이 한 데 모여 포트폴리오를 이루는데, 이 포트폴리오의 자산가치의 변동은 금융자산 각각의 움직임도 중요하지만, 서로 어떤 경향으로 움직이는 지도 아주 중요합니다.
위에서 살펴본 원리를 이용하여 여러 개의 금융자산이 움직이는 것을 어떻게 모델링하는지 차차 알아보기로 하겠습니다.
'수학의 재미 > 행렬 이론' 카테고리의 다른 글
| n차원 상관계수행렬을 촐레스키 분해하면? (0) | 2022.09.20 |
|---|---|
| 촐레스키 분해 (0) | 2022.09.19 |
| 3중 대각행렬의 풀이 (0) | 2022.07.29 |
| 상관계수와 상관계수 행렬 #2 (0) | 2022.07.06 |
| 상관계수와 상관계수 행렬 #1 (0) | 2022.07.05 |
댓글