Numpy線形代数その2
台湾から帰還。。
ｼｭｯｼｬｼﾀｸﾅｲ。。。
今日の朝はこれで終了。
（平均分散アプローチについて学んでたら終わった）
import numpy as np import numpy.linalg as lin import matplotlib.pyplot as plt #最小分散ポートフォリオの計算 Mu=np.array([1.0,3.0,1.5,6.0,4.5]) Stdev=np.array([5.0,10.0,7.5,15.0,11.0]) CorrMatrix = np.array([[1.00, 0.25, 0.18, 0.10, 0.25], [0.25, 1.00, 0.36, 0.20, 0.20], [0.18, 0.36, 1.00, 0.25, 0.36], [0.10, 0.20, 0.25, 1.00, 0.45], [0.25, 0.20, 0.36, 0.45, 1.00]]) Sigma=np.diag(Stdev).dot(CorrMatrix).dot(np.diag(Stdev))
- np.diag
np.diagは初めて。
In [1]: import numpy as np
In [2]: help(np.diag)
Help on function diag in module numpy.lib.twodim_base:
diag(v, k=0)
Extract a diagonal or construct a diagonal array.
See the more detailed documentation for ``numpy.diagonal`` if you use this
function to extract a diagonal and wish to write to the resulting array;
whether it returns a copy or a view depends on what version of numpy you
are using.
Parameters
----------
v : array_like
If `v` is a 2-D array, return a copy of its `k`-th diagonal.
If `v` is a 1-D array, return a 2-D array with `v` on the `k`-th
diagonal.
k : int, optional
Diagonal in question. The default is 0. Use `k>0` for diagonals
above the main diagonal, and `k<0` for diagonals below the main
diagonal.
Returns
-------
out : ndarray
The extracted diagonal or constructed diagonal array.
See Also
--------
diagonal : Return specified diagonals.
diagflat : Create a 2-D array with the flattened input as a diagonal.
trace : Sum along diagonals.
triu : Upper triangle of an array.
tril : Lower triangle of an array.
Examples
--------
>>> x = np.arange(9).reshape((3,3))
>>> x
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> np.diag(x)
array([0, 4, 8])
>>> np.diag(x, k=1)
array([1, 5])
>>> np.diag(x, k=-1)
array([3, 7])
>>> np.diag(np.diag(x))
array([[0, 0, 0],
[0, 4, 0],
[0, 0, 8]])
ぶっちゃけよくわからなかったので、試して見た
Mu=np.array([1.0,3.0,1.5,6.0,4.5]) Stdev=np.array([5.0,10.0,7.5,15.0,11.0]) CorrMatrix = np.array([[1.00, 0.25, 0.18, 0.10, 0.25], [0.25, 1.00, 0.36, 0.20, 0.20], [0.18, 0.36, 1.00, 0.25, 0.36], [0.10, 0.20, 0.25, 1.00, 0.45], [0.25, 0.20, 0.36, 0.45, 1.00]]) print(Stdev) print (np.diag(Stdev))
結果
[ 5. 10. 7.5 15. 11. ] [[ 5. 0. 0. 0. 0. ] [ 0. 10. 0. 0. 0. ] [ 0. 0. 7.5 0. 0. ] [ 0. 0. 0. 15. 0. ] [ 0. 0. 0. 0. 11. ]]
こんなんがあるのか
引数であるStdevを対角要素にもつ対角行列を作るんすね
続けていく
Sigma=np.diag(Std).dot(CorrMatrix).dot(np.diag(Std)) iota=np.ones(Mu.shape) inv_Sigma=lin.inv(Sigma) print (Sigma) print (Mu.shape) print (iota) print (inv_Sigma)
結果
[[ 25. 12.5 6.75 7.5 13.75 ] [ 12.5 100. 27. 30. 22. ] [ 6.75 27. 56.25 28.125 29.7 ] [ 7.5 30. 28.125 225. 74.25 ] [ 13.75 22. 29.7 74.25 121. ]] (5,) [1. 1. 1. 1. 1.] [[ 4.47987105e-02 -4.51004983e-03 -1.24347718e-03 7.17380912e-04 -4.40574739e-03] [-4.51004983e-03 1.21306899e-02 -4.86483457e-03 -8.70624550e-04 3.52683222e-05] [-1.24347718e-03 -4.86483457e-03 2.28222730e-02 -8.18315526e-04 -4.07386288e-03] [ 7.17380912e-04 -8.70624550e-04 -8.18315526e-04 5.70192503e-03 -3.22127446e-03] [-4.40574739e-03 3.52683222e-05 -4.07386288e-03 -3.22127446e-03 1.17353428e-02]]
- Numpy配列ブロードキャスト
Mu=np.array([1.0,3.0,1.5,6.0,4.5])#期待値 Std=np.array([5.0,10.0,7.5,15.0,11.0]) #分散 CorrMatrix = np.array([[1.00, 0.25, 0.18, 0.10, 0.25], [0.25, 1.00, 0.36, 0.20, 0.20], [0.18, 0.36, 1.00, 0.25, 0.36], [0.10, 0.20, 0.25, 1.00, 0.45], [0.25, 0.20, 0.36, 0.45, 1.00]]) Sigma=np.diag(Std).dot(CorrMatrix).dot(np.diag(Std))#分散共分散行列 iota=np.ones(Mu.shape) inv_Sigma=lin.inv(Sigma) A=Mu.dot(inv_Sigma).dot(iota) B=Mu.dot(inv_Sigma).dot(Mu)
このあと、このように続いていって、13行目からNumpy行列の掛け算が始まるんだけど、ここでつまづいた
1dNumpy配列（ベクトル）と2dNumpy配列の掛け算って、どう考えるんだ？？？
いや、ブロードキャストと言われても、いまいちピンときていない、、
明日簡単な配列の計算でもしてみながら慣れてみよう、
現状、足し算の時は理解できる
a=np.array([[1,2], [3,4]]) b=np.array([3,4]) print (a+b)
結果
[[4 6] [6 8]]
b.shapeが(2,)→(1,2)にまず変換（次元を揃える）
さらに(1,2)→(2,2)に変換（要素数を揃える）
だから、bは最終的に
b=np.array([3,4], [3,4])
として、aに足される
じゃあ、掛け算はどうかというと
掛け算って次元とか要素数を全く同じにする必要なくないかと思ってしまう
うーん、続く
冷静に、そりゃあ行列の積で要素数合わせる必要ないよな
a=np.array([[1,2], [3,4]]) b=np.array([3,4]) print (a.dot(b)) print (b.dot(a))
結果
[11 25] [15 22]
np.dot
スカラーが絡むと、成分ごとの積となります。 ベクトルとベクトルの積は、内積となります。 その他はすべて行列積となります。その際、左がベクトルの場合は行ベクトル、右がベクトルの場合は列ベクトルと解釈されます。
そうらしい
結局、ブロードキャストは関係なくて、Numpy行列の積についての話だったなあ
ってか、行列の詳しい話すっとばしてとりあえずコード写すみたいなことしてたら何もわからなくなった、、、
行列をちょっと勉強した、続き
minSigma=1/np.sqrt(C) F_Sigma=np.linspace(minSigma,max(Std)*1.05,250) F_PosiReturn=(A+np.sqrt(D*(C*F_Sigma**2-1)))/C F_NegaReturn=(A-np.sqrt(D*(C*F_Sigma**2-1)))/C
結果
pyfin_mvf_example1.py:34: RuntimeWarning: invalid value encountered in sqrt F_PosiReturn=(A+np.sqrt(D*(C*F_Sigma**2-1)))/C pyfin_mvf_example1.py:35: RuntimeWarning: invalid value encountered in sqrt F_NegaReturn=(A-np.sqrt(D*(C*F_Sigma**2-1)))/C
計算結果は出たんだけど、なんかエラーというほどでもRuntimeWarningとやらが出てきた
ルートの中に負の値が入ってるみたいだ
minSigma=1/np.sqrt(C) F_Sigma=np.linspace(minSigma,max(Std)*1.05,250) F_PosiReturn=(A+np.sqrt(D*(C*F_Sigma**2-1)))/C F_NegaReturn=(A-np.sqrt(D*(C*F_Sigma**2-1)))/C print (F_Sigma) print (D*(C*F_Sigma**2-1))
結果
[ 4.44212563 4.48753878 4.53295193 4.57836508 4.62377823 4.66919138 4.71460453 4.76001768 4.80543083 4.85084398 4.89625713 4.94167028 4.98708343 5.03249658 5.07790973 5.12332288 5.16873603 5.21414918 5.25956233 5.30497548 5.35038863 5.39580178 5.44121493 5.48662808 5.53204123 5.57745438 5.62286753 5.66828068 5.71369383 5.75910698 5.80452013 5.84993328 5.89534643 5.94075958 5.98617273 6.03158588 6.07699903 6.12241218 6.16782533 6.21323848 6.25865163 6.30406478 6.34947793 6.39489108 6.44030423 6.48571738 6.53113053 6.57654368 6.62195683 6.66736998 6.71278313 6.75819628 6.80360943 6.84902258 6.89443573 6.93984888 6.98526203 7.03067518 7.07608833 7.12150148 7.16691463 7.21232778 7.25774093 7.30315408 7.34856723 7.39398038 7.43939353 7.48480668 7.53021983 7.57563298 7.62104613 7.66645928 7.71187243 7.75728558 7.80269873 7.84811188 7.89352503 7.93893818 7.98435133 8.02976448 8.07517763 8.12059078 8.16600393 8.21141708 8.25683023 8.30224338 8.34765653 8.39306968 8.43848283 8.48389598 8.52930913 8.57472228 8.62013543 8.66554858 8.71096173 8.75637488 8.80178803 8.84720118 8.89261433 8.93802749 8.98344064 9.02885379 9.07426694 9.11968009 9.16509324 9.21050639 9.25591954 9.30133269 9.34674584 9.39215899 9.43757214 9.48298529 9.52839844 9.57381159 9.61922474 9.66463789 9.71005104 9.75546419 9.80087734 9.84629049 9.89170364 9.93711679 9.98252994 10.02794309 10.07335624 10.11876939 10.16418254 10.20959569 10.25500884 10.30042199 10.34583514 10.39124829 10.43666144 10.48207459 10.52748774 10.57290089 10.61831404 10.66372719 10.70914034 10.75455349 10.79996664 10.84537979 10.89079294 10.93620609 10.98161924 11.02703239 11.07244554 11.11785869 11.16327184 11.20868499 11.25409814 11.29951129 11.34492444 11.39033759 11.43575074 11.48116389 11.52657704 11.57199019 11.61740334 11.66281649 11.70822964 11.75364279 11.79905594 11.84446909 11.88988224 11.93529539 11.98070854 12.02612169 12.07153484 12.11694799 12.16236114 12.20777429 12.25318744 12.29860059 12.34401374 12.38942689 12.43484004 12.48025319 12.52566634 12.57107949 12.61649264 12.66190579 12.70731894 12.75273209 12.79814524 12.84355839 12.88897154 12.93438469 12.97979784 13.02521099 13.07062414 13.11603729 13.16145044 13.20686359 13.25227674 13.29768989 13.34310304 13.38851619 13.43392934 13.4793425 13.52475565 13.5701688 13.61558195 13.6609951 13.70640825 13.7518214 13.79723455 13.8426477 13.88806085 13.933474 13.97888715 14.0243003 14.06971345 14.1151266 14.16053975 14.2059529 14.25136605 14.2967792 14.34219235 14.3876055 14.43301865 14.4784318 14.52384495 14.5692581 14.61467125 14.6600844 14.70549755 14.7509107 14.79632385 14.841737 14.88715015 14.9325633 14.97797645 15.0233896 15.06880275 15.1142159 15.15962905 15.2050422 15.25045535 15.2958685 15.34128165 15.3866948 15.43210795 15.4775211 15.52293425 15.5683474 15.61376055 15.6591737 15.70458685 15.75 ] [-9.94290886e-19 1.84051059e-04 3.69974157e-04 5.57769293e-04 7.47436468e-04 9.38975682e-04 1.13238693e-03 1.32767023e-03 1.52482556e-03 1.72385292e-03 1.92475233e-03 2.12752378e-03 2.33216726e-03 2.53868278e-03 2.74707034e-03 2.95732994e-03 3.16946158e-03 3.38346526e-03 3.59934098e-03 3.81708873e-03 4.03670852e-03 4.25820036e-03 4.48156423e-03 4.70680014e-03 4.93390808e-03 5.16288807e-03 5.39374009e-03 5.62646416e-03 5.86106026e-03 6.09752840e-03 6.33586858e-03 6.57608080e-03 6.81816506e-03 7.06212135e-03 7.30794969e-03 7.55565006e-03 7.80522247e-03 8.05666692e-03 8.30998341e-03 8.56517194e-03 8.82223250e-03 9.08116511e-03 9.34196975e-03 9.60464643e-03 9.86919515e-03 1.01356159e-02 1.04039087e-02 1.06740735e-02 1.09461104e-02 1.12200193e-02 1.14958003e-02 1.17734533e-02 1.20529783e-02 1.23343754e-02 1.26176445e-02 1.29027856e-02 1.31897988e-02 1.34786840e-02 1.37694413e-02 1.40620706e-02 1.43565719e-02 1.46529453e-02 1.49511907e-02 1.52513082e-02 1.55532977e-02 1.58571592e-02 1.61628928e-02 1.64704984e-02 1.67799760e-02 1.70913257e-02 1.74045474e-02 1.77196412e-02 1.80366070e-02 1.83554449e-02 1.86761547e-02 1.89987366e-02 1.93231906e-02 1.96495166e-02 1.99777146e-02 2.03077847e-02 2.06397268e-02 2.09735410e-02 2.13092272e-02 2.16467854e-02 2.19862157e-02 2.23275180e-02 2.26706923e-02 2.30157387e-02 2.33626571e-02 2.37114476e-02 2.40621101e-02 2.44146446e-02 2.47690512e-02 2.51253298e-02 2.54834804e-02 2.58435031e-02 2.62053979e-02 2.65691646e-02 2.69348034e-02 2.73023143e-02 2.76716972e-02 2.80429521e-02 2.84160790e-02 2.87910780e-02 2.91679491e-02 2.95466922e-02 2.99273073e-02 3.03097944e-02 3.06941536e-02 3.10803849e-02 3.14684881e-02 3.18584634e-02 3.22503108e-02 3.26440302e-02 3.30396216e-02 3.34370851e-02 3.38364206e-02 3.42376281e-02 3.46407077e-02 3.50456593e-02 3.54524830e-02 3.58611786e-02 3.62717464e-02 3.66841862e-02 3.70984980e-02 3.75146818e-02 3.79327377e-02 3.83526656e-02 3.87744656e-02 3.91981376e-02 3.96236816e-02 4.00510977e-02 4.04803858e-02 4.09115460e-02 4.13445782e-02 4.17794824e-02 4.22162587e-02 4.26549070e-02 4.30954274e-02 4.35378198e-02 4.39820842e-02 4.44282207e-02 4.48762292e-02 4.53261097e-02 4.57778623e-02 4.62314869e-02 4.66869836e-02 4.71443523e-02 4.76035930e-02 4.80647058e-02 4.85276906e-02 4.89925475e-02 4.94592764e-02 4.99278773e-02 5.03983503e-02 5.08706953e-02 5.13449123e-02 5.18210014e-02 5.22989625e-02 5.27787957e-02 5.32605009e-02 5.37440781e-02 5.42295274e-02 5.47168487e-02 5.52060421e-02 5.56971075e-02 5.61900449e-02 5.66848544e-02 5.71815359e-02 5.76800895e-02 5.81805150e-02 5.86828127e-02 5.91869823e-02 5.96930240e-02 6.02009378e-02 6.07107236e-02 6.12223814e-02 6.17359112e-02 6.22513131e-02 6.27685871e-02 6.32877331e-02 6.38087511e-02 6.43316411e-02 6.48564032e-02 6.53830373e-02 6.59115435e-02 6.64419217e-02 6.69741720e-02 6.75082942e-02 6.80442886e-02 6.85821549e-02 6.91218933e-02 6.96635038e-02 7.02069862e-02 7.07523408e-02 7.12995673e-02 7.18486659e-02 7.23996365e-02 7.29524792e-02 7.35071939e-02 7.40637807e-02 7.46222395e-02 7.51825703e-02 7.57447731e-02 7.63088480e-02 7.68747950e-02 7.74426140e-02 7.80123050e-02 7.85838680e-02 7.91573031e-02 7.97326103e-02 8.03097894e-02 8.08888406e-02 8.14697639e-02 8.20525592e-02 8.26372265e-02 8.32237659e-02 8.38121773e-02 8.44024607e-02 8.49946162e-02 8.55886437e-02 8.61845433e-02 8.67823149e-02 8.73819585e-02 8.79834742e-02 8.85868619e-02 8.91921217e-02 8.97992535e-02 9.04082573e-02 9.10191331e-02 9.16318811e-02 9.22465010e-02 9.28629930e-02 9.34813570e-02 9.41015931e-02 9.47237012e-02 9.53476813e-02 9.59735335e-02 9.66012577e-02 9.72308540e-02 9.78623222e-02 9.84956626e-02 9.91308749e-02 9.97679594e-02 1.00406916e-01 1.01047744e-01 1.01690445e-01 1.02335017e-01 1.02981462e-01 1.03629779e-01]
なぜ2つ目の配列の最初の項が負になっているんだ
print (F_Sigma[0]) print (C) print (C*F_Sigma[0]**2) print (C*F_Sigma[0]**2-1)
4.442125625150906 0.05067786697220465 0.9999999999999999 -1.1102230246251565e-16
配列の初項が負となっているのは、普通の結果っぽい
なので、とりあえずルートの中身を絶対値にした
F_PosiReturn=(A+np.sqrt(abs(D*(C*F_Sigma**2-1))))/C F_NegaReturn=(A-np.sqrt(abs(D*(C*F_Sigma**2-1))))/C print (F_PosiReturn)
[1.34613316 1.61383466 1.72568172 1.81215788 1.88560503 1.95078981 2.01014999 2.06512965 2.11666705 2.16541206 2.2118365 2.25629553 2.29906429 2.34036092 2.38036179 2.41921182 2.45703174 2.4939234 2.52997357 2.56525692 2.59983817 2.63377391 2.66711388 2.6999021 2.73217773 2.76397577 2.79532765 2.82626175 2.8568037 2.88697685 2.91680244 2.9462999 2.97548706 3.0043803 3.03299472 3.06134429 3.08944192 3.1172996 3.14492847 3.17233893 3.19954067 3.22654274 3.25335362 3.27998127 3.30643315 3.33271628 3.35883728 3.38480237 3.41061745 3.43628806 3.46181946 3.48721665 3.51248432 3.53762698 3.56264886 3.58755403 3.61234631 3.63702938 3.66160674 3.6860817 3.71045746 3.73473704 3.75892336 3.78301919 3.80702719 3.83094989 3.85478976 3.87854911 3.9022302 3.92583519 3.94936614 3.97282504 3.99621381 4.01953429 4.04278825 4.06597739 4.08910336 4.11216774 4.13517207 4.15811781 4.18100639 4.20383918 4.22661751 4.24934266 4.27201587 4.29463833 4.31721122 4.33973564 4.36221268 4.3846434 4.40702882 4.42936991 4.45166764 4.47392294 4.4961367 4.51830979 4.54044307 4.56253736 4.58459345 4.60661211 4.62859411 4.65054018 4.67245101 4.69432732 4.71616976 4.73797899 4.75975566 4.78150037 4.80321374 4.82489634 4.84654875 4.86817153 4.88976523 4.91133035 4.93286744 4.95437697 4.97585946 4.99731536 5.01874516 5.0401493 5.06152822 5.08288236 5.10421215 5.12551799 5.14680028 5.16805943 5.18929581 5.2105098 5.23170177 5.25287207 5.27402105 5.29514906 5.31625644 5.33734351 5.35841058 5.37945798 5.40048601 5.42149497 5.44248515 5.46345684 5.48441032 5.50534587 5.52626375 5.54716423 5.56804758 5.58891403 5.60976384 5.63059726 5.65141451 5.67221584 5.69300147 5.71377163 5.73452653 5.75526639 5.77599143 5.79670185 5.81739784 5.83807963 5.85874739 5.87940131 5.9000416 5.92066843 5.94128198 5.96188243 5.98246996 6.00304474 6.02360693 6.04415671 6.06469422 6.08521964 6.10573312 6.12623481 6.14672486 6.16720343 6.18767066 6.20812668 6.22857165 6.2490057 6.26942896 6.28984158 6.31024367 6.33063538 6.35101682 6.37138813 6.39174942 6.41210081 6.43244243 6.4527744 6.47309682 6.4934098 6.51371347 6.53400793 6.55429329 6.57456965 6.59483711 6.61509579 6.63534578 6.65558718 6.67582009 6.69604461 6.71626082 6.73646883 6.75666872 6.77686059 6.79704453 6.81722061 6.83738894 6.85754959 6.87770265 6.8978482 6.91798632 6.9381171 6.9582406 6.97835692 6.99846612 7.01856829 7.03866349 7.05875181 7.0788333 7.09890805 7.11897612 7.13903759 7.15909252 7.17914098 7.19918304 7.21921876 7.2392482 7.25927143 7.27928852 7.29929952 7.3193045 7.33930351 7.35929662 7.37928389 7.39926537 7.41924111 7.43921119 7.45917564 7.47913453 7.49908792 7.51903584 7.53897837 7.55891554 7.57884741 7.59877403 7.61869546 7.63861173 7.6585229 7.67842902 7.69833014]
最後に結果載せておく
import numpy as np import numpy.linalg as lin import matplotlib.pyplot as plt Mu=np.array([1.0,3.0,1.5,6.0,4.5])#期待値 Std=np.array([5.0,10.0,7.5,15.0,11.0]) #分散 CorrMatrix = np.array([[1.00, 0.25, 0.18, 0.10, 0.25], [0.25, 1.00, 0.36, 0.20, 0.20], [0.18, 0.36, 1.00, 0.25, 0.36], [0.10, 0.20, 0.25, 1.00, 0.45], [0.25, 0.20, 0.36, 0.45, 1.00]]) Sigma=np.diag(Std).dot(CorrMatrix).dot(np.diag(Std))#分散共分散行列 iota=np.ones(Mu.shape) inv_Sigma=lin.inv(Sigma) A=Mu.dot(inv_Sigma).dot(iota) B=Mu.dot(inv_Sigma).dot(Mu) C=iota.dot(inv_Sigma).dot(iota) D=B*C-A**2 V_Target=np.linspace(Mu.min(),Mu.max(),num=5) V_Risk=np.zeros(V_Target.shape) V_Weight=np.zeros((V_Target.shape[0],Mu.shape[0])) for idx,Target_return in enumerate(V_Target): V_Weight[idx,:]=(C*Target_return-A)/D*lin.solve(Sigma,Mu)\ +(B-A*Target_return)/D*lin.solve(Sigma,iota) V_Risk[idx]=(C/D)*(Target_return-A/C)**2+1.0/C #可視化 minSigma=1/np.sqrt(C) F_Sigma=np.linspace(minSigma,max(Std)*1.05,250) F_PosiReturn=(A+np.sqrt(abs(D*(C*F_Sigma**2-1))))/C F_NegaReturn=(A-np.sqrt(abs(D*(C*F_Sigma**2-1))))/C fig=plt.figure(1) plt.plot(F_Sigma,F_PosiReturn,'k-') plt.plot(F_Sigma,F_NegaReturn,'k--') plt.show()
np.diag
引数であるStdevを対角要素にもつ対角行列を作る
対角行列は行列の各行や各列へ同じ値をかけるときに利用する。
例
Std=np.array([5.0,10.0,7.5,15.0,11.0]) test=np.diag(Std) test0=np.ones((5,5)) print (test) print (test0) print (test0.dot(test)) print (test.dot(test0))
結果
[[ 5. 0. 0. 0. 0. ] [ 0. 10. 0. 0. 0. ] [ 0. 0. 7.5 0. 0. ] [ 0. 0. 0. 15. 0. ] [ 0. 0. 0. 0. 11. ]] [[1. 1. 1. 1. 1.] [1. 1. 1. 1. 1.] [1. 1. 1. 1. 1.] [1. 1. 1. 1. 1.] [1. 1. 1. 1. 1.]] [[ 5. 10. 7.5 15. 11. ] [ 5. 10. 7.5 15. 11. ] [ 5. 10. 7.5 15. 11. ] [ 5. 10. 7.5 15. 11. ] [ 5. 10. 7.5 15. 11. ]] [[ 5. 5. 5. 5. 5. ] [10. 10. 10. 10. 10. ] [ 7.5 7.5 7.5 7.5 7.5] [15. 15. 15. 15. 15. ] [11. 11. 11. 11. 11. ]]
