상관관계를 가지는 표준정규분포 변수 여러개 만들기(python)

 

지난 글

2022.09.20 - [수학의 재미/행렬 이론] - n차원 상관계수행렬을 촐레스키 분해하면?

n차원 상관계수행렬을 촐레스키 분해하면?

에서 상관계수를 가지는 $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보다 작죠( 여기 참고.) 따라서 표준정규분포라 가정할 수 있습니다.

 

 

 

결론

위와 같이 서로 독립인 표준정규분포에서 서로 상관관계가 있는(우리가 원하는 상관계수를 가지는) 새로운 표준정규분포 난수를 얻을 수 있습니다. 금융공학에서는 여러 개의 금융 자산이 한 데 모여 포트폴리오를 이루는데, 이 포트폴리오의 자산가치의 변동은 금융자산 각각의 움직임도 중요하지만, 서로 어떤 경향으로 움직이는 지도 아주 중요합니다.

위에서 살펴본 원리를 이용하여 여러 개의 금융자산이 움직이는 것을 어떻게 모델링하는지 차차 알아보기로 하겠습니다.