Я имею дело с некоторыми симуляциями в Python, и у меня есть следующий код:
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.sandbox.distributions import multivariate as mt
# Simulation parameters
n = 100 # number of assets
T = 360 # time-series size
p = 100 # number of controls
b = 4 # number of relevant controls
# Risk price of newly proposed factor
lambdag1 = 0.5
# Other risk prices
lambdah1 = 0.5 * np.random.multivariate_normal(np.zeros(b), cov = 0.1 * np.identity(b)).reshape((b, 1))
# Zero-beta rate
gamma_0 = 0.5
# Simulating moments (covariances, expected returns and betas)
# Computing cross-sectional residuals (nx1)
Ce1 = 0.5 * np.random.multivariate_normal(mean = np.zeros(n), cov = np.identity(n)).reshape((n,1))
# Computing Ch1 (nxb)
Ch1 = 0.1 * np.transpose(np.random.multivariate_normal(mean = np.zeros(n), cov = np.identity(n), size = b))
# Computing Cg1 (nx1)
xi1 = 0.5 * np.random.normal()
chi1 = 0.5 * np.random.multivariate_normal(mean = np.zeros((b)), cov = np.identity(b)).reshape((1, b))
Cg1 = np.ones((n,1))*np.transpose(xi1) + Ch1.dot(np.transpose(chi1)) + Ce1
# Computing Cz1
eta1 = np.random.multivariate_normal(np.zeros(b), cov = np.identity(b)).reshape((1,b))
Cz1 = Cg1 - Ch1.dot(np.transpose(eta1))
# Computing the cross-section of expected returns
Er = np.ones((n,1)) * gamma_0 + Cg1.dot(lambdag1) + Ch1.dot(lambdah1)
# Computing the betas
# Computing the correlation structure for shocks z1
rhoz1 = 0.5
# Defines the covariance matrix for the important controls. Some correlation is desired to challenge the LASSO.
# Too much "independence" would be unreal. Using a non-random co-variance matrix
rhoh1 = 0.5
Sigmah1 = np.zeros((b, b))
for i in range(0, b):
for j in range(0, b):
Sigmah1[i, j] = 0.5 * rhoh1**(np.abs(i - j))
# Use orthogonality conditions to derive the betas
# betag1 = Cz1 * 1/Sigmaz1[0,0]
betag1 = Cz1 * 1/rhoz1
betah1 = Ch1.dot(np.linalg.inv(Sigmah1)) - betag1.dot(eta1)
# Creating Sigmau
rhou = 0.5
Sigmau = np.zeros((n,n))
for i in range(0, n):
for j in range(0, n):
Sigmau[i, j] = rhou**(np.abs(i - j))
# Defining covariances for h2
theta0 = np.zeros((1, p -b))
red_fac = int(np.floor(0.5*(p-b))) # indices of redundant factors in h2
theta1 = 0.5 * np.hstack((np.transpose(np.random.multivariate_normal(mean = np.zeros(b + 1),
cov = np.identity(b+1),
size = red_fac)),
np.zeros((b+1, p - b - red_fac))))
Cepsilon = np.random.multivariate_normal(mean = np.zeros(p-b), cov = 0.5 * np.identity(p-b), size = n)
# The factor shrinking covariances is a way magnitudes agree. This is important for LASSO
Ch2 = 0.1 * (np.ones((n,1)).dot(theta0) + np.hstack((Cg1, Ch1)).dot(theta1) + Cepsilon)
# Implied covariance matrix for returns
Sigmar = (Cz1 * 1/rhoz1).dot(np.transpose(Cz1)) + Ch1.dot(np.linalg.inv(Sigmah1)).dot(np.transpose(Ch1)) + Sigmau
phi = np.linalg.inv(Sigmar).dot(Ch2)
После этого я определяю следующую функцию:
def gen_data(n, p, T, b):
"""This function simulates one sample dataset of returns and factors from previously computed covariances and betas.
Outputs:
R (n x T): matrix of returns for n assets and T periods
H (p X T): matrix of p factors and T periods
G (1 X T): matrix of the evolution of the newly proposed factor (T periods)
"""
# Computing the important controls
H1 = np.transpose(np.random.multivariate_normal(mean = np.zeros(b), cov = Sigmah1, size = T))
# Computing the shocks z1t and the factor g_t stacked in G
Z1 = np.random.multivariate_normal(mean = np.zeros(T), cov = np.identity(T)).reshape(1, T)
G = eta1.dot(H1) + Z1
# Simulating the shocks
# Reference for multivariate Student-t: https://github.com/statsmodels/statsmodels/blob/master/statsmodels/sandbox/distributions/multivariate.py
U = mt.multivariate_t_rvs(m = np.zeros(n), S = Sigmau, df = 5, n = T).T
# Finally simulating the returns
R = Er.dot(np.ones((1, T))) + betag1.dot(G) + betah1.dot(H1) + U
# Simulating other factors
H2 = Ch2.T.dot(np.linalg.inv(Sigmar)).dot(R - Er.dot(np.ones((1,T)))) + np.random.multivariate_normal(mean = np.zeros(p-b), cov = np.identity(p-b), size = T).T
H = np.vstack((H1, H2))
# Outputs
return(R, G, H)
И теперь я вызываюэто просто, чтобы проверить, что все в порядке
a, b, c = gen_data(n, p, T, b)
Пока все хорошо.Все имеет желаемые размеры и так далее.Ошибка начинает появляться, когда я запускаю вышеупомянутый вызов еще раз.Если я это сделаю, у меня будет следующая ошибка:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-11-cb3db1328026> in <module>()
----> 1 a, b, c = gen_data(n, p, T, b)
2 #print(time.time() - tic)
<ipython-input-10-6f0fd2e2ab70> in gen_data(n, p, T, b)
10
11 # Computing the important controls
---> 12 H1 = np.transpose(np.random.multivariate_normal(mean = np.zeros(b), cov = Sigmah1, size = T))
13
14 # Computing the shocks z1t and the factor g_t stacked in G
TypeError: only integer scalar arrays can be converted to a scalar index
Сначала я не могу понять сообщение об ошибке.Во-вторых, я не могу понять, как он работает нормально в первый раз и не работает во второй раз.Это не имеет смысла для меня.Извините, если это глупо, я новичок в Python.Есть идеи по этому поводу?Заранее большое спасибо!