在 Python (SciPy) 中模擬相關股票收益

我希望在 Python 中生成具有股票間相關性的股票回報。但是，輸出行為不正常，並且可能具有意外的時間相關性，從而導致問題。

此程式碼旨在生成num_paths的相關股票收益，給定 Panda 的年化收益系列、恆定共變異數的 DataFrame 和日期的 DateIndex ( date_index )。

from pandas import DataFrame, concat
from scipy.stats import multivariate_normal

def correlated_returns(num_paths, returns, covariances, date_index, periods_per_year=1):
   period_returns = (1 + returns) ** (1 / periods_per_year) - 1 if periods_per_year != 1 else returns
   mn = multivariate_normal(period_returns, covariances / periods_per_year, allow_singular=True)

   digits = len(str(num_paths))
   paths = [DataFrame(mn.rvs(size=len(date_index)), index=date_index, columns=returns.index) for _ in range(num_paths)]
   keys = [f'Run {str(run_num).zfill(digits)}' for run_num in range(num_paths)]
   return concat(paths, axis='columns', keys=keys, names=['Run', 'Returns'])

我基於以下共變異數矩陣的測試程式碼在結果中顯示了一些可能相關的奇怪之處。

correlation = 0.2  # inter-stock correlation 0.18
annualized_return = 7 / 100  # Simulated return for each stock
stocks = [f'Stock {i}' for i in range(simulated_stocks)]
constituent_weights = DataFrame(1 / simulated_stocks, date_index, stocks)

returns = Series(annualized_return, stocks)
volatilities = Series(volatility, stocks)
correlations = DataFrame(correlation, stocks, stocks)
fill_diagonal(correlations.values, 1)
covariances = correlations.mul(volatilities, axis='index').mul(volatilities, axis='columns')

在使用相等的組成權重（上圖）的大量模擬中， index_returns似乎是負自相關的（在 Pandas 中使用 autocorr）。當相關性大於 0 時，該指數的年化回報率低於預期的 7%，但當相關性為 0 時等於 7%。非常奇怪*。*

for simulation in range(simulations):
   runs = correlated_returns(1, returns, covariances, date_index, frequency_scale)
   return_history = runs['Run 0']

   index_returns = return_history.mul(constituent_weights).sum(axis='columns') \
       .div(constituent_weights.sum(axis='columns'))

   annualized_return = (1 + index_returns[1:]).prod() ** (1 / simulation_years) - 1

我是否以某種方式濫用multivariate_normal？

我沒有執行你的程式碼，但這是 Jensen 不等式成立的一個例子。假設您查看單個資產在兩個時期內的回報， $ R_1 $ 和 $ R_2 $ . 此處的回報是指總回報，例如在您的情況下為 1.07。如果這些回報是獨立的，我們有 $ \mathrm{E}(R_1R_2) $ 等於 $ \mathrm{E}(R_1)\mathrm{E}(R_2) $ . 但是，這個產品不是你所看到的。您通過獲取年化回報 $ n $ -th 根，它是總回報的凹函式。定義 $ f $ 成為你的年化函式。然後，由 Jensen 不等式， $ \mathrm{E}(f(R_1R_2)) \le f(\mathrm{E}(R_1R_2)) $ . 如果您假設收益的特定分佈，您可以量化這種差距（例如，對於對數正態分佈的價格，幾何平均值將是算術平均值減去變異數的一半）。

一個數值範例（帶有 R 程式碼）：

ans <- NULL
for (n in c(1e3, 1e4, 1e5, 1e6, 1e7)) {

   R1 <- pmax(0, rnorm(n,mean = 1.07, sd = 0.2))
   R2 <- pmax(0, rnorm(n,mean = 1.07, sd = 0.2))

   ans <- rbind(ans, 
                c(n, 
                  mean(R1), mean(R2), mean(R1*R2), 
                  mean((R1*R2)^(1/2))))

}
colnames(ans) <- c("trials", "E[R1]", "[ER2]", "E[R1R2]", "E[f(R1R2)]")
ans
##        trials    E[R1]    [ER2]  E[R1R2] E[f(R1R2)]
## [1,]     1000 1.065672 1.076422 1.148000   1.061579
## [2,]    10000 1.066393 1.074114 1.145864   1.060501
## [3,]   100000 1.069893 1.069176 1.143758   1.059819
## [4,]  1000000 1.070105 1.069802 1.144723   1.060283
## [5,] 10000000 1.069981 1.070046 1.144939   1.060348

您會看到總回報 ( E[R1R2]) 符合預期（1.07 乘以 1.07 為 1.1449），但平均年化回報低於 0.07。

---

更新：Jensen 的不等式並不嚴格（即使是，差距仍然是微不足道的）。差距將取決於回報的分佈。你的例子應該接近對數正態的情況（它會隨著更多的句點變得更接近）。因此，差距是單期均值減去變異數的一半。在 500 種資產且沒有相關性的情況下，變異數將下降到幾乎為零，因此算術平均值和幾何平均值應該幾乎相同。隨著相關性的增加，變異數也會增加，因此幾何平均數（即年化回報）的估計值將減少。

這是一個包含 12 個句點和不斷增加的相關級別的程式碼範例

$$ I’d rather keep the first example, because I really simple examples :-) $$. 該功能randomReturns取自NMOF我維護的包。它返回一個大小為 12 乘以 500 的正態分佈變數的矩陣。

library("NMOF")
rhos <- c(0, 0.1, 0.3, 0.9)  ## levels of correlation to test
na <- 500                    ## number of assets
n <- 1e4                     ## number of trials
w <- rep(1/na, na)           ## equal weights

ans <- NULL
for (rho in rhos) {
   message(rho)
   
   results.geom <- numeric(n)
   results.TR <- numeric(n)
   for (i in seq_len(n)) {
       R <- randomReturns(na = na,
                          ns = 12,
                          mean = 0.07,
                          sd = 0.2,
                          rho = rho)
       results.geom[i] <- prod(R %*% w + 1)^(1/12) - 1
       results.TR[i]   <- prod(R %*% w + 1) - 1
   }

   ## compute E: expected annualized
   ##            return under lognormality
   C <- array(rho, dim = c(na, na))
   diag(C) <- 1
   E <- 0.07 - 0.5 * sum(0.2^2 * C)/na^2
   
   ans <- rbind(ans, 
                c(rho, mean(results.TR), mean(results.geom), E))

}
colnames(ans) <- c("rho", "E[ΣR]", "E[f(ΣR)]", "ln-expected")
ans
##      rho  E[R1R2] E[f(R1R2)] ln-expected
## [1,] 0.0 1.252762 0.06998791    0.069960
## [2,] 0.1 1.258257 0.06843583    0.067964
## [3,] 0.3 1.253640 0.06471616    0.063972
## [4,] 0.9 1.244869 0.05353727    0.051996

1.07^12 - 1 = 1.252在 10000 個樣本之後，對於所有相關級別，總回報大致符合預期 ( )。平均年化回報相當接近於對數正態下的預期。

引用自：https://quant.stackexchange.com/questions/68766