Pythonで長期積立投資シミュレーション

# 目的
- Pythonを用いて特定のリターン・リスク（標準偏差）を持つ資産に毎月積立投資を行った場合の長期資産推移を、モンテカルロ法によるシミュレーションで算定する

# 計算の概要
- `Pandas`、`Numpy`を用いて、対象期間における月次の対数収益率を、正規分布に従う乱数として生成する
- 得られた乱数を用いて、毎月積立投資を行った場合の資産推移をパーセンタイル値として出力する

<font color="red">(2024/1/29追記）</font>
- 実行例をWebアプリ化し、手軽に試せるようにしたものを投稿しました。

https://qiita.com/yasuyuki-s/items/4d22c9ab1f5005b64ae5

# 実行例
- 今回作成した解析用のクラスである`asset_model`の実行例を示します。
- 実装方法は後述します。

## 1. 資産のリターン・リスクを設定する
```python
# 期待対数収益率、標準偏差（%,年換算）
mu = 8
s = 20

#リターンmu、リスクsのモデルを定義
model = asset_model(mu,s)
```
- 対象資産の年次の対数収益率の`期待値（平均）`、`標準偏差`を入力します
- 今回の例では、全世界株式(`ACWI`)を想定した値（期待リターン8%、リスク20%)を設定しています。
- 任意のポートフォリオのリターン・リスクの算定には下記記事もご参考に。

https://qiita.com/yasuyuki-s/items/b0732596bf77f03c1a3d

## 2. モンテカルロシミュレーションを実行する
```python
# 投資期間（年）
dur_y = 10
# 試行回数
n = 20000

#投資期間dur_y、試行回数nのモンテカルロシミュレーションを実行
model.run_mc(dur_y,n)
```
- `投資期間`、モンテカルロ法の`試行回数`を入力しモンテカルロシミュレーションを実施します。
- 詳細は後述しますが、内部処理としては`投資期間(月換算)`×`試行回数`個の月次対数収益率を、前段で入力した期待値、標準偏差の正規分布に従い生成しています。
- 投資期間が長くなると、計算にはちょっと時間がかかります。

## 3. 積立投資を行った場合の資産推移を算定
```python
# 初期投資額
x_init = 0
# 毎月積立額
delta_m = 10

#初期投資x_init、毎月delta_mの積立投資を行った場合の資産推移を計算
model.exercise(x_init,delta_m)
```
- `初期投資額`、`毎月積立額`を入力し、前段で生成した対数収益率のデータを用いて、長期資産推移を算定します。
- 今回の例では、初期投資額を0、毎月の積立額を10としています。
### 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/099e1acc-f45b-4168-d80c-0edfa992ef62.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/1122c356-a622-f56d-2d2d-db358367a1ef.png)
- モンテカルロシミュレーションに基づき、年毎の総資産額の分布の平均値、P10～P90のパーセンタイル値、累計投資額（元本）をそれぞれ出力しています
- モデルと対象期間が同じであれば、`exercise`のみを再度実行することで、異なる初期投資額、積立額での計算が可能です。
- 例えば、上記の例における10年間の総投資額と同じ金額を、全て一括で初期投資し、積立は0とする場合は下記の通りになります。
```python
# 初期投資額
x_init = 1200
# 毎月積立額
delta_m = 0

#初期投資x_init、毎月delta_mの積立投資を行った場合の資産推移を計算
model.exercise(x_init,delta_m)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/f3061df3-4215-bc7f-f8b9-12f623e456e2.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/8347e991-841c-590e-491d-89f83bf684f9.png)

- 今回の例では、全額を一括投資するケースは積立投資のケースに比べ、運用期間の前半においてより大きな下振れリスクがある一方、10年後時点では下振れリスクに大きな違いはなく、より大きなアップサイドを期待できる結果となりました。

# 実装方法
- ソースコード全文は下記を参照してください

https://github.com/yasuyuki-s/Long_term_investment/

- `asset_model`クラスの構成は以下の通りです
```python
#ライブラリ読み込み
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from IPython.display import display

#対数収益率の正規分布を仮定した長期投資シミュレーションモデル
class asset_model:
    #初期化時の処理
    def __init__(self, mu_y, s_y):

    #対数収益率の長期モンテカルロシミュレーションを実施
    def run_mc(self, dur_y, n):
        
    #初期投資額x_init, 毎月積立額delta_mを投資した場合の資産推移を計算
    def exercise(self, x_init, delta_m):
```

## 計算の流れ
1. 年次リターン・リスクを月換算
2. 投資期間、試行回数に応じた数の対数収益率に相当する乱数を生成
2. 投入資金の、各時点における推移を計算
2. 計算結果を出力

- 積立投資の長期資産推移の計算の考え方は以下の通りです。
- ある時点$i$における資産価格を$P_i$としたとき、対数収益率$lr$は以下の通り表されます。
```math
lr_{i+1}=\ln \frac{P_{i+1}}{P_{i}}  
```
- 上記の対数収益率の定義および対数の性質より、($t>i$)となる時点$i$で投資した金額$X_i$の、時点$t$における金額$X_{i,t}$は以下の通り求められます。
```math
\begin{align}
X_{i,t} &=  X_i * (\frac{X_{i+1}}{X_{i}}*\frac{X_{i+2}}{X_{i+1}}* \dots *\frac{X_{t}}{X_{t-1}})  \\\
&=X_i * \exp (lr_{i+1}+lr_{i+2}+\dots+lr_t)
\end{align}
```
- 以上より、時点$t$における積立投資の累計投資額$Capital_t$および総資産額$TotalAsset_t$は以下の通り求められます。
```math
Capital_t = \sum_{i=0}^{t-1}X_i 
```
```math
TotalAsset_t = \sum_{i=0}^{t-1}X_{i,t} 
```
- 各時点の投資額$X_i$は入力値として与えられます。
- モンテカルロシミュレーションでは、対数収益率$lr$を乱数により生成します。


## 1. 年次リターン・リスクを月換算
```python:class asset_model
def __init__(self, mu_y, s_y):
    #パーセント表記から少数へ変換
    mu_y = mu_y / 100
    s_y = s_y / 100
    #入力された対数収益率の年次リターン・リスク（標準偏差）を月次へ変換
    self.mu_m = mu_y /12
    self.s_m = s_y / np.sqrt(12)
```
- 引数として入力された年あたりの対数収益率の平均・標準偏差を月あたりへ変換し、格納します。
## 2. 投資期間、試行回数に応じた数の対数収益率に相当する乱数を生成

```python:class asset_model
    def run_mc(self, dur_y, n):
        #入力された投資期間を年→月へ変換
        dur_m = dur_y * 12
        # [(投資期間_月),(試行回数)]の正規分布N(μ,σ)に従う乱数行列として毎月の対数収益率を生成
        A1 = np.random.normal(self.mu_m, self.s_m,size=[dur_m,n])        
        #次元を拡張し、[(投資期間_月),(投資期間_月),(試行回数)]の行列とする
        A2 = np.tile(A1,[dur_m,1,1])
        #行列を転置し、[(試行回数),(投資期間_月),(投資期間_月)]の行列とする
        A3 = np.transpose(A2,(2,0,1))
        #各試行の上三角行列を取得
        A4 = np.triu(A3)
        #計算結果を格納
        self.A = A4
        self.dur_m = dur_m
        print("Monte Carlo Calculation of "+ str(dur_y) + " years finished")
```
- `A1 = np.random.normal(self.mu_m, self.s_m,size=[dur_m,n]) `
 前段で計算した平均`mu_m`、標準偏差`s_m`に従う正規分布の乱数行列を作成します
 行列のサイズは[`投資期間_月`,`試行回数`]となります
- 以下、投資期間1年(`dur_m=12`）、試行回数`n=5`の例を示します
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/f6a30253-65cc-7364-7f63-0119b72bfb04.png)
- `A2 = np.tile(A1,[dur_m,1,1])`
- `A3 = np.transpose(A2,(2,0,1))`
- `A4 = np.triu(A3)`
 それぞれ、行列に以下のような変形操作を加え、最終的に[`試行回数`,`投資期間
 _月`,`投資期間_月`]の形状の上三角行列を得ます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/8580ad8d-5fbb-b601-0ef8-6912415d225d.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/5918c5a7-0470-23a2-e902-1a70be5ea49a.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/92ae32ee-8050-315d-d4a5-d7920f58b4d0.png)



## 3. 投入資金の、各時点における推移を計算
```python:class asset_model
def exercise(self, x_init, delta_m):
    #月別の投資額のリスト
    x = delta_m * np.ones(self.dur_m)
    #初期投資額を追加
    x[0] = x[0] + x_init
    
    #総資産額のリスト
    ta_list = []
    #累計投資額のリスト
    capital = [] 
    
    #総資産額の推移を1年毎(12ヶ月毎)に集計する
    for i in range(11, self.dur_m,12):
        #集計時点までの対数収益率行列を切り出し
        A_tmp = self.A[:,:i+1, :i+1] 
        #集計時点までの投資額のリストを切り出し
        x_tmp = x[:i+1]
        #対数収益率行列の各行の和を計算
        A_tmp = A_tmp.sum(axis=2)
        #対数収益率から、価格変動比へ変換
        A_tmp = np.exp(A_tmp)
        #各試行の、投資期間経過後の資産額を計算
        ta = np.dot(A_tmp,x_tmp)
        #計算結果を格納
        ta_list.append(ta)
        
        #累計投資額を格納
        capital.append(x_tmp.sum())
```
- 以下、投資期間1年(`dur_m=12`）、試行回数`n=5`の例の続きです。
- 前段で得られた上三角行列は、$i$ヶ月目の月次の対数収益率を$lr_i$とすると、次のように表されます。
```math
\boldsymbol{A}=
\begin{bmatrix}
lr_1 & lr_2&  \dots & lr_{12} \\\  
0 &lr_2&\dots & lr_{12} \\\ 
\vdots & \vdots & \ddots & \vdots \\\ 
0 & 0 & \dots & lr_{12} \\\ 
\end{bmatrix}
```
- 行列$\boldsymbol{A}$の各行を足して、さらに指数関数を掛けて得られる列ベクトル$\boldsymbol{a}$と、各月の積立投資額の列ベクトル$\boldsymbol{x}$の内積を計算することで、$t=12$の総資産額$TotalAsset_t$を求めています。
```math
\boldsymbol{a}=
\begin{bmatrix}
\exp (lr_1+\dots+lr_{12}) \\\  
\exp (lr_2+\dots+lr_{12}) \\\ 
\vdots \\\ 
\exp (lr_{12}) \\\ 
\end{bmatrix} \quad
\boldsymbol{x}=
\begin{bmatrix}
X_0\\\  
X_1\\\ 
\vdots \\\ 
X_{11}\\\ 
\end{bmatrix}
```
```math
\begin{align}
TotalAsset_{12} &= \boldsymbol{a} \cdot \boldsymbol{x} \\\
&= X_0 *\exp (lr_1+\dots+lr_{12})+\dots+X_{11}*\exp(lr_{12})
\end{align}
```
- 実際の計算処理上は、これにモンテカルロ法の試行回数の次元（例では`n=5`)が加わります。
- 以下の例では、各月の積立額は簡単のため`X=1`としています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/1d8f66cd-b09a-3483-9526-437bf86c4657.png)



## 4. 計算結果を出力
```python:def exercise(続き）
    #1年毎の各種統計値を計算
    #平均
    mean = np.mean(ta_list, axis=1)
    #パーセンタイル値（np.percetileは下位○%を入力するが、出力としては上位○%の表記とする）
    p10 = np.percentile(ta_list,90,axis=1)
    p25 = np.percentile(ta_list,75,axis=1)
    p50 = np.percentile(ta_list,50,axis=1)
    p75 = np.percentile(ta_list,25,axis=1)
    p90 = np.percentile(ta_list,10,axis=1)
    
    year = np.arange(self.dur_m/12)+1
    
    #計算結果をDataFrameに格納
    df_result = pd.DataFrame(np.transpose([year,mean,p10,p25,p50,p75,p90,capital]),columns=["Year","Mean","P10","P25","P50","P75","P90","Capital"])
    pd.options.display.float_format = '{:.1f}'.format
    display(df_result)
        
    #グラフ出力
    left, width = 0.1, 0.65 
    bottom, height = 0.1, 0.65 
    fig_main = [left, bottom, width, height] #散布図の定義
    fig =plt.figure(figsize=(10,7)) 
    ax = fig.add_axes(fig_main)
    
    ax.fill_between(df_result["Year"],df_result["P10"],df_result["P90"], alpha=0.3, label="Range P10-P90")
    ax.fill_between(df_result["Year"],df_result["P25"],df_result["P75"], alpha=0.3, label="Range P25-P75")
    ax.plot(df_result["Year"],df_result["P50"], label="P50")
    ax.plot(df_result["Year"], df_result["Capital"], label="Capital")
    ax.set_ylim(0, max(df_result["P10"])*1.1)
    ax.set_xlim(1,max(df_result["Year"]))
    ax.set_xlabel("Year")
    ax.set_ylabel("Total assets")
    ax.grid(True)
    ax.legend()    
```
- 計算結果の平均や各種パーセンタイル値を整理した上で、`matplotlib`を用いて図化しています