Pythonでモンテカルロ法を用いた建設コストリスクの定量化（多変量解析）

# やりたいこと
- 設備投資を行う際の建設コストのリスクを、モンテカルロ法を用いて定量化する
  - 建設コストを構成する各要素の確率分布と相関関係を設定し、乱数組を生成する

設備投資を行う際、全体費用の大きな割合を占める建設コストは、天候・物価・人件費などの様々な要因によって大きく変動するという不確実性を有しています。
ここでは、建設コストを各要素に分解し、その確率分布と相関関係を設定したうえで、モンテカルロ法により建設コストの変動リスクを定量的に評価することを試みます。

![image_s.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/02c5fc62-f650-81ca-4c3d-6b7692f76f16.png)

# モデル
簡単のため、建設工事は二つの洋上工事（**Offshore_1**,**Offshore_2**）および二つの陸上工事（**Onshore_1**,**Onshore_2**）によって構成されるものとし、建設コストはこれら個別のコストの合計と定義します。各工事コストのP10・最頻値・P90は以下の通り設定します。

<table>
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>P10</th>
      <th>mode</th>
      <th>P90</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Offshore_1</th>
      <td>100</td>
      <td>250</td>
      <td>400</td>
    </tr>
    <tr>
      <th>Offshore_2</th>
      <td>300</td>
      <td>350</td>
      <td>800</td>
    </tr>
    <tr>
      <th>Onshore_1</th>
      <td>900</td>
      <td>925</td>
      <td>950</td>
    </tr>
    <tr>
      <th>Onshore_2</th>
      <td>300</td>
      <td>550</td>
      <td>600</td>
    </tr>
  </tbody>
</table>

また、各工事コストの相関関係は以下の通りとします。
（洋上工事同士・陸上工事同士はそれぞれ強い相関[0.9]があり、洋上⇔陸上間では弱い相関[0.5]があることをイメージしてモデル化しています）

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Offshore_1</th>
      <th>Offshore_2</th>
      <th>Onshore_1</th>
      <th>Onshore_2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Offshore_1</th>
      <td>1.0</td>
      <td>0.9</td>
      <td>0.5</td>
      <td>0.5</td>
    </tr>
    <tr>
      <th>Offshore_2</th>
      <td>0.9</td>
      <td>1.0</td>
      <td>0.5</td>
      <td>0.5</td>
    </tr>
    <tr>
      <th>Onshore_1</th>
      <td>0.5</td>
      <td>0.5</td>
      <td>1.0</td>
      <td>0.9</td>
    </tr>
    <tr>
      <th>Onshore_2</th>
      <td>0.5</td>
      <td>0.5</td>
      <td>0.9</td>
      <td>1.0</td>
    </tr>
  </tbody>
</table>

# 乱数生成
各要素の確率分布は三角分布に従うものとして、設定した入力条件を満足する乱数を生成します。
三角分布のクラス**triangular_fit**は[こちらのコード](https://qiita.com/yasuyuki-s/items/6c43dd5541978ce7778f)を使用しています。
また、乱数の生成方法について[詳しくはこちら](https://qiita.com/yasuyuki-s/items/6e54abe7fee69b677c1e)を参照してください。

```Python:multi_correlated_triangular.py
import numpy as np
from triangular_dist import triangular_fit
from scipy.stats import norm,spearmanr,triang
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

#試行回数
n = 50000

#各要素の名前
name = ["Offshore_1","Offshore_2","Onshore_1","Onshore_2"]
#P10,最頻値、P90のリスト
elements = np.array([[100,250,400],
                     [300,350,800],
                     [900,925,950],
                     [300,550,600]])
#各要素の相関関係
cor_coef = np.array([[1.0,0.9,0.5,0.5],
                     [0.9,1.0,0.5,0.5],
                     [0.5,0.5,1.0,0.9],
                     [0.5,0.5,0.9,1.0]])

#三角分布の準備
tri_dists = []
for ele in elements:
    #与えられたP10,mode,P90にフィットする三角形分布を生成
    tri_dists.append(triangular_fit(ele[0],ele[1],ele[2]))

#相関のある多変量正規分布の乱数z_listを生成
mu = np.zeros(elements.shape[0])
sigma = cor_coef
z_list = np.random.multivariate_normal(mu,sigma,n)
#正規累積分布関数から、一様分布の乱数u_listを生成
u_list= norm.cdf(z_list)
#逆関数法を用いてu_listから任意の分布を持つ乱数を生成
a_list = []
for i in range(u_list.shape[1]):
    #三角分布に従う乱数を生成
    a_list.append(tri_dists[i].ppf(u_list[:,i]))
a_list = np.transpose(a_list)

#DataFrameを生成
df = pd.DataFrame(a_list,columns=name)
#生成した乱数のリストを出力
df
```

><table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Offshore_1</th>
      <th>Offshore_2</th>
      <th>Onshore_1</th>
      <th>Onshore_2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>271.674899</td>
      <td>423.619032</td>
      <td>922.862146</td>
      <td>436.283326</td>
    </tr>
    <tr>
      <th>1</th>
      <td>178.707851</td>
      <td>354.904091</td>
      <td>882.976324</td>
      <td>217.474613</td>
    </tr>
    <tr>
      <th>2</th>
      <td>128.821805</td>
      <td>403.111938</td>
      <td>915.485525</td>
      <td>448.719875</td>
    </tr>
    <tr>
      <th>3</th>
      <td>351.455486</td>
      <td>734.054223</td>
      <td>920.970889</td>
      <td>521.520472</td>
    </tr>
    <tr>
      <th>4</th>
      <td>442.248604</td>
      <td>903.573617</td>
      <td>948.630335</td>
      <td>551.485066</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>49995</th>
      <td>147.568593</td>
      <td>395.741033</td>
      <td>899.264791</td>
      <td>331.173835</td>
    </tr>
    <tr>
      <th>49996</th>
      <td>382.085926</td>
      <td>649.281076</td>
      <td>912.453528</td>
      <td>315.281751</td>
    </tr>
    <tr>
      <th>49997</th>
      <td>340.035343</td>
      <td>627.800319</td>
      <td>926.774642</td>
      <td>458.619111</td>
    </tr>
    <tr>
      <th>49998</th>
      <td>187.235699</td>
      <td>463.307672</td>
      <td>917.633978</td>
      <td>490.614796</td>
    </tr>
    <tr>
      <th>49999</th>
      <td>131.474132</td>
      <td>387.047341</td>
      <td>922.874209</td>
      <td>449.360773</td>
    </tr>
  </tbody>
</table>
<p>50000 rows × 4 columns</p>

# 生成した乱数の可視化
生成した乱数の散布図行列を出力してみます。

```python:
#生成した乱数の散布図行列を出力
pd.plotting.scatter_matrix(df,figsize=(10, 10),hist_kwds={"bins":30},s=5)
plt.show()
```

>![scatter_matrix.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/3960cf54-952c-8cb6-2c62-e3a1a07d9f8b.png)

対角成分は各変数のヒストグラム、それ以外は散布図を示しています。
設定した入力条件を満足してそうですね。

念のため相関係数も確認してみます。

```python:
#スピアマンの順位相関係数を算定、ヒートマップで出力
plt.figure(figsize=(10, 10)) 
sns.heatmap(df.corr(method='spearman'),vmax=1, vmin=-1,square=True, annot=True, cbar=False)
plt.show()
```

>![coeff.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/7803a439-a543-1913-d907-ac4221ff65a3.png)

概ね入力条件で設定した相関関係の通りですね。

# 建設コストの算定
得られた乱数列の各行を足し合わせて、合計の建設コストを算定します。

```python:
#データフレームの各行の合計値を計算
sum_column = df.sum(axis=1)

#合計値のパーセンタイル値を計算
percentile = [10,25,50,75,90]
df_percentile = pd.DataFrame(np.percentile(sum_column,percentile), 
                      index=["P10","P25","P50","P75","P90"], columns=["Overall Cost"]).T

#ヒストグラムを出力
plt.hist(sum_column, bins=30, density=1)
plt.title("Overall Cost")
plt.show()

#合計値のパーセンタイル値を出力
df_percentile
```

>![overall.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/11edd9e9-580f-39e5-b0c9-2531f12f1c93.png)

><table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>P10</th>
      <th>P25</th>
      <th>P50</th>
      <th>P75</th>
      <th>P90</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Overall Cost</th>
      <td>1679.657667</td>
      <td>1885.347016</td>
      <td>2140.773801</td>
      <td>2429.011128</td>
      <td>2676.415322</td>
    </tr>
  </tbody>
</table>

これで合計の建設コストの分布を得る事ができました。
あとはリスク許容度に応じて、適切なパーセンタイル値を選択すればよさそうです。

#重回帰分析
得られた建設コストに対して、個別の要素がどの程度影響を及ぼしているのかを評価するため、
重回帰分析を行い標準偏回帰係数を比較します。重回帰分析にはscikit-learnを使用します。

```python:
from sklearn import linear_model
from scipy.stats import zscore

#重回帰分析モデルを作成
clf = linear_model.LinearRegression()

#データフレームの各列を標準化
df_normalized=zscore(df, ddof=1)
sum_normalized = zscore(sum_column, ddof=1)
#説明変数に標準化した各行のデータを使用
X = df_normalized
#目的変数に標準化した各行の合計値を使用
Y = sum_normalized

#モデルにフィッティング
clf.fit(X,Y)

#偏回帰係数を計算
df_reg_coef = pd.DataFrame({"Elements":df.columns,
                    "Regression Coefficients":np.abs(clf.coef_)}).sort_values(by='Regression Coefficients',
                    ascending=False) 
#降順でグラフに描画
sns.barplot(x="Regression Coefficients",y="Elements",data=df_reg_coef,palette="Blues_r")
plt.show()

#偏回帰係数を出力
df_reg_coef
```

>![regression.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1350685/47ef506d-b27a-2243-17d0-26aec6cb8991.png)

><table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Elements</th>
      <th>Regression Coefficients</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>1</th>
      <td>Offshore_2</td>
      <td>0.501799</td>
    </tr>
    <tr>
      <th>3</th>
      <td>Onshore_2</td>
      <td>0.301040</td>
    </tr>
    <tr>
      <th>0</th>
      <td>Offshore_1</td>
      <td>0.296749</td>
    </tr>
    <tr>
      <th>2</th>
      <td>Onshore_1</td>
      <td>0.049565</td>
    </tr>
  </tbody>
</table>

今回のケースでは、全体の建設コストは、

- P10～P90の変動幅が最も大きい**Offshore_2**の影響を大きく受けている
- 一方で、絶対値は大きいものの変動幅が小さい**Onshore_1**の影響は相対的に小さい

ということが見て取れます。

# 最後に
今回のケースでは、建設コストは各要素の単純足し合わせとしましたが、例えば**[単価][数量]**といった要素を設定し、

>[建設コスト]=[単価]×[数量]

とすれば、単価と数量の変動を考慮することも可能ですし、それらを組み合わせることで柔軟にモデルを構築可能です。

# 参考にしたサイト
- [scikit-learn で線形回帰 (単回帰分析・重回帰分析)](https://pythondatascience.plavox.info/scikit-learn/%E7%B7%9A%E5%BD%A2%E5%9B%9E%E5%B8%B0)
- [scatter_matrix（スカッター行列）について](https://qiita.com/mitsumizo/items/efddc80f166e900f97ef)