S+NUOPT の例題をご紹介します。以下の例題は RNUOPTでも同様に解くことができますが、出力結果及びグラフに関しては RNUOPT のものとは異なる部分があります。
例題の動かし方
- サンプルモデル及びデータは全て S+NUOPT 及び RNUOPT に組み入れられており、それぞれの例題の最後に示してある「実行例」はコマンドウィンドウに直接コピー＆ペーストして実行することができます。
- ただし、以下の例題を実行する前に一度だけ S+NUOPT であれば module(nuopt)、RNUOPT であれば library(Rnuopt) を実行する必要があります。
- S+NUOPT で実行例をお試しいただくためには、お手許のマシンに S-PLUS 及び S+NUOPT がインストールされている必要があります。
- また、RNUOPT でお試しいただくためには、お手許のマシンに R 及び R パッケージ Rnuopt がインストールされている必要があります。
- S+NUOPT 及び RNUOPT の試用版のお申し込みについてはこちらからお問い合わせ下さい。（試用版は 30 日間限定で機能制限はありません）
目次
- 基本的なマルコビッツモデル
- さまざまなリスク尺度によるポートフォリオ最適化
- 絶対偏差をリスクとしたポートフォリオ最適化
- １次の下方部分積率（LPM）をリスクとしたポートフォリオ最適化
- CVaRをリスクとしたポートフォリオ最適化
- コンパクト分解（大規模ポートフォリオ最適化）
- 端株処理
- 銘柄グルーピング
- Maximum Drawdown
- Sharpe Ratio
- 半正定値行列の取得
- ロバスト最適化
- イールドカーブのフィッティング
- 格付け推移行列の推定
- ロジスティック回帰
基本的なマルコビッツモデル
4 銘柄の時系列データを用いて、マルコビッツモデル：
最小化
（収益率の分散）
制約
（最低期待収益率の確保）
, 
（空売りの禁止）
による最適ポートフォリオを求めます。ここでは
を1% とします。
元にするのは４銘柄の６０期間分の収益率データ R.60x4 です（S+NUOPTにサンプルとして含まれています）。S-PLUSでグラフ化してみましょう。
> matplot(rownames(R.60x4),R.60x4,pch=colnames(R.60x4),type="o")
|
次にS-PLUSの機能を使って最適化の材料となる分散
を求めます。
> Q <- var(R.60x4)
> Q
A B C D
A 0.017788838 0.005730923 0.004552905 0.003779794
B 0.005730923 0.026978759 0.005153147 0.008770707
C 0.004552905 0.005153147 0.005200091 0.004284923
D 0.003779794 0.008770707 0.004284923 0.014216287
平均
も求めましょう。
> rbar <- apply(R.60x4,2,mean)
> rbar
A B C D
-0.01483683 0.0147342 -0.01805572 0.01391076
NUOPTが認識するためには、S-PLUSのvector型オブジェクトはarray型として定義しておく必要があります。
> rbar <- as.array(rbar)
最適化モデルはS-PLUSの手続きとして定義しておきます（S+NUOPTにサンプルとして含まれています）。
> Marko
function(Q,rbar,rmin)
{
Asset <- Set()
Q <- Parameter(index=dprod(Asset,Asset),Q)
rbar <- Parameter(index=Asset,rbar))
i <- Element(set=Asset)
j <- Element(set=Asset)
x <- Variable(index=Asset)
risk <- Objective(type="minimize")
risk ~ Sum(x[j]*Q[i,j]*x[i],i,j)
Sum(rbar[j]*x[j],j) >= rmin
Sum(x[j],j) == 1
x[j] >= 0
}
モデルにデータを代入するには次のように Systemコマンドを使います。この命令で最適化問題が定義されます。
> s <- System(model=Marko,Q,rbar,0.01) Evaluating Marko(Q,rbar,0.01) ... ok! Expanding objective (1/4 name="risk") Expanding constraint (2/4) Expanding constraint (3/4) Expanding constraint (4/4)
展開した結果はsというオブジェクトになりますので、次のように内容を表示することができます。
> s 1-1 : -0.0148368*x[A]+0.0147342*x[B]-0.0180557*x[C]+0.0139108*x[D] >= 0 .01 2-1 : x[A]+x[B]+x[C]+x[D] == 1 3-1 : x[A] >= 0 3-2 : x[B] >= 0 3-3 : x[C] >= 0 3-4 : x[D] >= 0 risk<objective>: 0.0177888*x[A]*x[A]+0.00573092*x[A]*x[B]+0.00455291*x[A ]*x[C]+0.00377979*x[A]*x[D]+0.00573092*x[A]*x[B]+0.0269788*x[B]*x[B]+0.0 0515315*x[B]*x[C]+0.00877071*x[B]*x[D]+0.00455291*x[A]*x[C]+0.00515315*x [B]*x[C]+0.00520009*x[C]*x[C]+0.00428492*x[C]*x[D]+0.00377979*x[A]*x[D]+ 0.00877071*x[B]*x[D]+0.00428492*x[C]*x[D]+0.0142163*x[D]*x[D] (minimize)
では解きましょう。
> sol <- solve(s)
NUOPT 11.1.9 (NLP/LP/IP/SDP module)
<with META-HEURISTICS engine "wcsp"/"rcpsp">
, Copyright (C) 1991-2008 Mathematical Systems Inc.
NUMBER_OF_VARIABLES 4
NUMBER_OF_FUNCTIONS 3
PROBLEM_TYPE MINIMIZATION
METHOD TRUST_REGION_IPM
<preprocess begin>..........<preprocess end>
<iteration begin>
res=1.0e+001 .... 8.3e-005 5.9e-007
<iteration end>
STATUS OPTIMAL
VALUE_OF_OBJECTIVE 0.01087261487
ITERATION_COUNT 6
FUNC_EVAL_COUNT 9
FACTORIZATION_COUNT 8
RESIDUAL 5.879111059e-007
ELAPSED_TIME(sec.) 0.00
解ベクトル
（組み入れ比率）を取り出します。オブジェクトsに対して、current（現在の値の取得）という手続きを使います。
> xopt <- as.array(current(s,x))
> xopt
A B C D
0.07466462 0.1986515 0.06031033 0.6663735
attr(, "indexes"):
[1] "*"
棒グラフにしてみましょう。
> barplot(names=rownames(xopt),xopt)
 |
ポートフォリオの組み入れ比率に従った場合の収益率の時系列(pf)を出してみます。
> pf <- R.60x4 %*% as.vector(xopt)
分散は
>var(pf)
[,1]
[1,] 0.01087261
とあって、上で表示されている目的関数値（VALUE_OF_OBJECTIVE)とほぼ一致しています。
元の収益率グラフにこのポートフォリオの収益率の変動を重ねてみてみます。
> Raug <- cbind(R.60x4,pf)
> matplot(rownames(Raug),Raug,pch=c(".",".",".",".","O"),type="o")
 |
"O"のプロットがあるのがポートフォリオで、変動が抑えられていることがわかります。
# 表示
matplot(rownames(R.60x4),R.60x4,pch=colnames(R.60x4),type="o")
# 前処理
Q <- var(R.60x4)
rbar <- apply(R.60x4,2,mean)
rbar <- as.array(rbar)
# 展開
sys <- System(Marko, Q, rbar, 0.01)
# 実行
sol <- solve(sys)
# 解の取得
xopt <- as.array(current(sys,x))
# 図示
barplot(names=rownames(xopt),xopt)
# 確認
pf <- R.60x4 %*% as.vector(xopt)
var(pf)
Raug <- cbind(R.60x4,pf)
matplot(rownames(Raug),Raug,pch=c(".",".",".",".","O"),type="o")
|
さまざまなリスク尺度によるポートフォリオ最適化
ポートフォリオ最適化におけるリスク尺度（変動の測り方）には分散以外にも様々なものがあります。ここでは、分散以外のリスク尺度として、絶対偏差、1 次の下方部分積率（LPM）、Conditional Value at Risk（CVaR）をリスク尺度としたポートフォリオ最適化問題を扱います。これらのリスク尺度を用いたモデルはすべて線形計画問題として記述できます。
ここで元にするデータは8000ケースの５銘柄の収益率を含むR.8000x5です。
まずは分散最小化モデル（コンパクト分解表現）を解いてみましょう。S+NUOPTではモデルを次のように定義します（S+NUOPTにサンプルとして含まれています）。
> MinVar
function(r.d)
{
Asset <- Set()
j <- Element(set=Asset)
Sample <- Set()
t <- Element(set=Sample)
r <- Parameter(index=dprod(t,j),r.d)
rb <- Parameter(index=j)
rb[j] ~ Sum(r[t,j],t)/nrow(r.d)
x <- Variable(index=j)
s <- Variable(index=t)
f <- Objective(type=minimize)
f ~ Sum(s[t]*s[t],t)/nrow(r.d)
Sum(x[j],j) == 1
x[j] >= 0
Sum((r[t,j]-rb[j])*x[j],j) == s[t]
}
次のようにして最適化を行います。
> # 展開（データの代入） > sys <- System(MinVar, R.8000x5) > # 最適化実行 > sol <- solve(sys) > # 解の取得 > x <- as.array(current(sys,x))
さて、結果の解析に入りましょう。まず最適ポートフォリオxを見てみます。
> # 図示 > pie(1000*x,names=dimnames(x)[[1]],inner=1.2)
 |
次に、ポートフォリオに従った場合の収益率の分布のヒストグラムを書いてみます。
> rp <- as.matrix(R.8000x5) %*% as.vector(x) > graphsheet() > hist(rp)
 |
8000ケースのうち、わずかですが、収益率が非常に悪いケース（収益率 -0.6 〜 -0.4）が存在します。対象としているアセットが社債でデフォルトするリスクが含まれている場合にこのようなことが起きます。このようなアセットに対して分散最小化モデルを適用するのは危険であるといえるでしょう。
# 展開 sys <- System(MinVar, R.8000x5) # 実行 sol <- solve(sys) # 解の取得 x <- as.array(current(sys,x)) # 図示 pie(1000*x,names=dimnames(x)[[1]],inner=1.2) rp <- as.matrix(R.8000x5) %*% as.vector(x) hist(rp) |
絶対偏差をリスクとしたポートフォリオ最適化
リスク尺度として絶対偏差を用いた絶対偏差最小化モデルを解いてみましょう。このモデルは絶対値関数 abs() を用いずに、線形計画問題として定式化するテクニックを利用しています。
> MinMad
function(r.d)
{
Asset <- Set()
j <- Element(set=Asset)
Sample <- Set()
t <- Element(set=Sample)
r <- Parameter(index=dprod(t,j),r.d)
rb <- Parameter(index=j)
rb[j] ~ Sum(r[t,j],t)/nrow(r.d)
x <- Variable(index=j)
s <- Variable(index=t)
u <- Variable(index=t)
v <- Variable(index=t)
f <- Objective(type=minimize)
f ~ Sum(u[t]+v[t],t)/nrow(r.d)
Sum(x[j],j) == 1
x[j] >= 0
Sum((r[t,j]-rb[j])*x[j],j) == s[t]
u[t] >= 0
v[t] >= 0
s[t] == u[t] - v[t]
}
次は最適化の実行手順です。
> # 展開 > sys <- System(MinMad, R.8000x5) > # 実行 > sol <- solve(sys) > # 解の取得 > x <- as.array(current(sys,x))
ヒストグラムを見てみましょう。
> rp <- as.matrix(R.8000x5) %*% as.vector(x) > hist(rp)
 |
分散最小化モデルと比べて、線形計画問題として定式化できるという利点がありますが、大崩に対する感度が小さく、収益率 -0.6 以下のケースが見受けられます。
# 展開 sys <- System(MinMad, R.8000x5) # 実行 sol <- solve(sys) # 解の取得 x <- as.array(current(sys,x)) # 図示 pie(1000*x,names=dimnames(x)[[1]],inner=1.2) rp <- as.matrix(R.8000x5) %*% as.vector(x) hist(rp) |
１次の下方部分積率（LPM）をリスクとしたポートフォリオ最適化
下方リスクモデルの一つである1 次の下方部分積率（LPM）最小化モデルの結果を見てみましょう。このモデルは収益率 rG 以上のサンプル点については関知しないかわりに、目標収益率 rG を下回る部分の平均を最小化しようとします。
> MinLPM1
function(n.sample, r.d, rG)
{
Asset <- Set()
j <- Element(set=Asset)
Sample <- Set()
t <- Element(set=Sample)
r <- Parameter(index=dprod(t,j),r.d)
rb <- Parameter(index=j)
rb[j] ~ Sum(r[t,j],t)/nrow(r.d)
x <- Variable(index=j)
s <- Variable(index=t)
f <- Objective(type=minimize)
f ~ Sum(s[t],t)/nrow(r.d)
Sum(x[j],j) == 1
x[j] >= 0
s[t] >= 0
rG <= Sum(r[t,j]*x[j],j) + s[t]
}
rGは下回ってほしくない（リスクとみなせる）収益率を与えます。ここでは-0.4としましょう。
> # 展開 > sys <- System(MinLPM1, R.8000x5, -0.4) > # 実行 > sol <- solve(sys) > # 解の取得 > x <- as.array(current(sys,x)) > rp <- as.matrix(R.8000x5) %*% as.vector(x) > graphsheet() > hist(rp)
 |
目標収益率 -0.4 を下回るケースが少なく、結果的に分散最小化モデルや絶対偏差最小化モデルと比べて大崩れするケースが少なくなっていることがわかります。
# 展開 sys <- System(MinLPM1, R.8000x5, -0.4) # 実行 sol <- solve(sys) # 解の取得 x <- as.array(current(sys,x)) # 図示 pie(1000*x,names=dimnames(x)[[1]],inner=1.2) rp <- as.matrix(R.8000x5) %*% as.vector(x) hist(rp) |
CVaRをリスクとしたポートフォリオ最適化
Conditional Value at Risk（CVaR）最小化モデルを解いてみましょう。CVaR は収益率の損失がある確率水準（パーセント点）
を上回るときの平均損失のことです。つまり、収益率の下位
（
は銘柄数）の平均損失を最小化します。この例では
、
ですので、収益率の下位 240 個のサンプル点の平均損失を最小化します。
> MinCVaR
function(r.d, beta)
{
Asset <- Set()
j <- Element(set=Asset)
Sample <- Set()
t <- Element(set=Sample)
r <- Parameter(index=dprod(t,j),r.d)
rb <- Parameter(index=j)
rb[j] ~ Sum(r[t,j],t)/nrow(r.d)
x <- Variable(index=j)
s <- Variable(index=t)
VaR <- Variable()
f <- Objective(type=minimize)
f ~ Sum(s[t],t)/(nrow(r.d)*(1-beta)) + VaR
Sum(x[j],j) == 1
x[j] >= 0
s[t] >= 0
Sum(r[t,j]*x[j],j) + VaR + s[t] >= 0
}
> # 展開
> sys <- System(MinCVaR, R.8000x5, 0.97)
> # 実行
> sol <- solve(sys)
> # 解の取得
> x <- as.array(current(sys,x))
> # 図示
> graphsheet()
> rp <- as.matrix(R.8000x5) %*% as.vector(x)
> hist(rp)
|
1 次の下方部分積率（LPM）最小化モデル同様、分散最小化モデルや絶対偏差最小化モデルと比べて大崩れするケースが少なくなっていることがわかります。
# 展開 sys <- System(MinCVaR, R.8000x5, 0.97) # 実行 sol <- solve(sys) # 解の取得 x <- as.array(current(sys,x)) # 図示 pie(1000*x,names=dimnames(x)[[1]],inner=1.2) rp <- as.matrix(R.8000x5) %*% as.vector(x) hist(rp) |
コンパクト分解（大規模ポートフォリオ最適化）
基本的なマルコビッツモデルの弱点として
（収益率の分散）
という式に表れる分散・共分散行列を取得しなければならないという点があります。昨今実務家が解かねばならない1000銘柄〜3000銘柄のポートフォリオではこのような分散・共分散行列を陽に生成すると計算負荷が増大してしまいます。
そのような場合には
：各銘柄の収益率のヒストリカルデータ（各列ｊは、期間
にわたる銘柄ｊの収益率の観測値を並べたもの）
を直接データとして入力、
と表されることを利用して、補助変数
を定義して次のような制約を導入し、
リスクを次のように書き換えます。
こうすれば、分散・共分散行列を直接生成しなくともよいので、効率的に計算を行うことができます。 1000銘柄のポートフォリオを求めてみます。データとして月次５年分（６０期間）のデータR.60x1000が得られているとします。
次のモデルでは分散共分散行列を作らず、R.60x1000を直接入力データとして最適化を行うことができます。
> MinVar
function(r.d)
{
Asset <- Set()
j <- Element(set=Asset)
Sample <- Set()
t <- Element(set=Sample)
r <- Parameter(index=dprod(t,j),r.d)
rb <- Parameter(index=j)
rb[j] ~ Sum(r[t,j],t)/nrow(r.d)
x <- Variable(index=j)
s <- Variable(index=t)
f <- Objective(type=minimize)
f ~ Sum(s[t]*s[t],t)/nrow(r.d)
Sum(x[j],j) == 1
x[j] >= 0
Sum((r[t,j]-rb[j])*x[j],j) == s[t]
}
> # 展開
> sys <- System(MinVar, R.60x1000)
> # 実行
> sol <- solve(sys)
> # 解の取得
> x <- as.array(current(sys,x))
同様のテクニックは収益率が
のようにファクターリターン
によって表現されているファクターモデルの場合にも利用することができます。このとき収益率の分散は次のように表現されますので
中間変数
なる線形制約式を導入し分散は
と表現することができます。
コンパクト分解は効率的であるのみならず、下方リスクモデルなど様々なリスク尺度によるポートフォリオモデルを記述する場合の基本となります。
# 展開 sys <- System(MinVar, R.60x1000) # 実行 sol <- solve(sys) # 解の取得 x <- as.array(current(sys,x)) # 図示 eps <- 1e-4 barplot(x[x>eps],names=names(x[x>eps])) |
端株処理
一般にポートフォリオ最適化問題においては、取引額を連続量として求めることが多いのですが、実際には各銘柄ごとの最小取引額の整数倍を取引しなければなりません。連続量として求められた取引額を何らかの形で最小取引額の整数倍に補正する作業のことを端株処理と呼びます。
紹介するモデルでは、各銘柄において端株を切り上げるなら 1 切り下げるなら 0 を取るような 0-1 整数 d を用いて、取引の総額が与えられた値の±1% の範囲に収まるという制約の下で、各サンプル点における収益額の絶対偏差を最小化しています。
ここでは、前項のコンパクト分解モデルによって得られた解を端株処理してみましょう。コンパクト分解モデルの解 x は既に得られおり、取引の総額を 50億円 とします。また、以下のような最小取引額データ RoundLot.unit が与えられているとします。
> RoundLot.unit
A0001 A0002 A0003 A0004 A0005 A0006 A0007 A0008 A0009
10880000 5400000 873000 1875000 7040000 4980000 3360000 15100000 5390000
A0010 A0011 A0012 A0013 A0014 A0015 A0016 A0017 A0018
2360000 4330000 6160000 994000 1385000 15150000 2200000 4390000 18000000
A0019 A0020 A0021 A0022 A0023 A0024 A0025 A0026 A0027
2760000 3790000 2395000 2730000 4430000 4160000 2960000 2265000 5410000
…
次が切り上げ、切り下げを決定する最適化モデルです。比較的規模の大きな０−１計画問題ですのでメタヒューリスティクス解法 wcsp を用いて解きます。
> RoundLot
function(r.d, unit.d, x.d, total)
{
Asset <- Set()
j <- Element(set=Asset)
Sample <- Set()
t <- Element(set=Sample)
r <- Parameter(index=dprod(t,j),r.d)
rb <- Parameter(index=j)
rb[j] ~ Sum(r[t,j],t)/nrow(r.d)
unit <- Parameter(index=j,unit.d)
base <- Parameter(index=j,(total*x.d)%/%unit.d*unit.d)
d <- IntegerVariable(type="binary",index=j)
fund <- Expression(index=j)
fund[j] ~ base[j] + unit[j]*d[j]
0.99*total <= Sum(fund[j],j) <= 1.01*total
s <- Expression(index=t)
s[t] ~ 0.01*Sum((r[t,j]-rb[j])*fund[j],j)
softConstraint(1,0,1)
s[t] == 0
}
> # 総額（50億円）
> total <- 5e9
> # 展開
> sys <- System(RoundLot, R.60x1000, RoundLot.unit, RoundLot.x, total)
次で解を求めます。wcsp を用いる場合には時間制限が必要です。
> # 実行（wcspを指定して解く） > nuopt.options(method="wcsp") > nuopt.options(maxtim=20) > sol <- solve(sys)
解が得られたら、次のようにして、解 dとそれを採用した場合の取引額を取得しましょう。
> # 解の取得 > d <- as.array(current(sys,d)) > fund <- as.array(current(sys,fund))
次のようにして、連続量から計算した理想的な投資量と具体的な投資量を定義します。
> n <- sum(fund > 0) > ideal <- (total* RoundLot.x)[fund > 0] # 理想的な投資量 > up <- d[fund > 0] > real <- fund[fund > 0] # 具体的な投資量
理想的な取引額(ideal)、切り上げるか否か(up)、実際の取引額(real)を表示します。
> cbind(ideal, up, real)
ideal up real
A0072 2.268026e+007 0 22275000
A0113 1.116033e-006 1 3740000
A0150 4.873316e+007 1 50140000
A0461 1.049896e+007 0 10296000
A0508 9.945964e+007 0 97500000
A0647 1.501497e+007 1 15030000
A0664 1.513937e+007 1 16450000
A0712 1.759865e-006 1 3880000
A0713 1.312477e+007 0 12705000
A0717 1.757595e+007 0 17492400
A0725 2.784842e-007 1 609000
A0746 1.929582e+008 0 188160000
A0751 1.009987e+008 0 99960000
A0870 3.162644e+007 1 32890000
A0875 7.202727e+006 1 9060000
A0907 2.333477e+007 0 17970000
A0919 2.447410e+007 0 21700000
A0935 7.131422e+006 0 6769000
A0945 4.145754e+007 0 38940000
A0962 1.749186e+007 1 18480000
A0975 1.265938e+009 0 1265850000
A0977 1.487018e+008 1 150960000
A0987 2.505445e+009 0 2505170000
ideal up real
A0997 385938255 0 385640000
グラフに描きます。
> barplot(rbind(ideal,real),beside=T,col=rbind(rep(2,n),rep(3,n)))
> legend(1,2.5e9,c("ideal","real"),fill=2:3)
 |
# 展開（総額 50 億）
total <- 5e9
sys <- System(RoundLot, R.60x1000, RoundLot.unit, RoundLot.x, total)
# 実行
nuopt.options(method="wcsp", maxtim=20)
sol <- solve(sys)
# 解の取得
fund <- as.array(current(sys,fund))
d <- as.array(current(sys,d))
# 比較
n <- sum(fund > 0)
ideal <- (total* RoundLot.x)[fund > 0]
up <- d[fund > 0]
real <- fund[fund > 0]
cbind(ideal,up,real)
barplot(rbind(ideal,real),beside=T,col=rbind(rep(2,n),rep(3,n)))
legend(1,2.5e9,c("ideal","real"),fill=2:3)
nuopt.options(nuopt.options(NA))
|
| ※S-PLUS や R でバッファに出力する設定をしている場合、ハングアップしているように見える場合がございますのでご注意下さい。 |
銘柄グルーピング
ポートフォリオ最適化後に適当な端株処理を行い、各銘柄に対する投資額が決まれば、あとは実際に買い／売り注文を出すだけです。しかし、一度に巨額の注文を出すのは、マーケットインパクト等を考慮すれば、あまり好ましいことではありません。そこで、銘柄を複数のグループに分けてグループ毎に注文を出すことがあります。このとき、各グループの性質がそれほど大きく異ならないようにしたい、というのは自然に発生する要求です。ここでは、各銘柄にはいくつかのファクター値（例えば、投資額）が与えられており、グループごとのファクター値が均等になるように銘柄をグルーピングすることを目的とします。
ここでは、50 個の銘柄を 5 個のグループに分けることを考えます。各銘柄は次のように、F0 から F9 までの 10 個のファクター値で特徴づけられるとします。
> Basket.fcoef
F0 F1 F2 F3 F4 F5 F6 F7 F8 F9
1 0.81 0.29 0.42 0.00 -0.14 0.03 0.04 1.83 -14.89 -0.20
2 1.29 0.39 0.41 0.04 0.10 0.18 -0.19 4.43 -15.35 0.23
3 1.00 -0.16 0.32 0.02 -0.02 0.25 0.14 2.12 -26.86 0.19
4 -2.75 -0.24 -0.40 0.01 0.29 -0.15 -0.40 -2.09 28.49 -0.18
5 -0.64 -0.12 -0.19 -0.14 0.00 0.22 -0.05 -1.84 24.50 0.72
6 -5.94 0.05 -0.20 0.03 -0.08 -0.24 -0.06 -1.57 49.74 -0.17
7 1.85 0.03 0.59 0.04 0.31 0.28 0.12 2.12 -42.61 -0.88
8 5.20 0.07 0.22 0.03 -0.07 0.01 0.19 -2.98 -2.71 -0.28
9 12.82 0.12 0.24 0.03 -0.10 0.13 0.47 3.40 -42.22 0.21
10 4.15 0.05 0.12 0.01 -0.05 0.09 0.10 -0.35 23.63 0.22
11 -32.92 0.12 -0.11 0.03 -0.12 -0.12 0.17 -1.81 16.64 -0.17
12 4.72 -0.20 -0.08 -0.03 -0.17 -0.04 0.20 3.84 -6.97 0.34
13 4.75 0.31 0.67 0.00 -0.01 0.48 -0.45 7.06 -31.22 0.18
14 -12.12 -0.07 -0.04 -0.07 0.36 -0.10 -0.21 -4.90 4.29 0.06
15 -7.92 -0.22 -0.26 -0.01 0.06 -0.01 -0.21 0.84 8.11 0.83
16 6.30 0.03 0.23 0.02 -0.07 0.19 0.00 2.46 -16.53 0.14
17 -30.70 -0.52 -0.04 -0.02 0.33 0.17 -0.22 -2.99 81.08 -0.18
18 -5.14 -0.09 -0.46 0.01 -0.17 -0.29 0.04 -1.05 64.62 -0.35
19 -5.92 0.23 -0.02 -0.01 0.08 -0.20 -0.23 -0.04 -11.09 -0.04
20 -12.54 -0.16 -0.06 0.00 0.18 0.06 -0.07 1.78 -2.60 -0.23
21 2.46 0.29 0.56 0.06 0.06 0.41 0.14 0.99 24.71 0.17
22 2.91 -0.26 0.16 0.01 -0.02 0.00 0.03 -1.26 12.49 0.16
23 3.78 0.11 0.16 0.00 -0.06 0.18 0.10 -0.66 10.73 0.14
F0 F1 F2 F3 F4 F5 F6 F7 F8 F9
24 1.83 0.16 0.35 0.01 0.03 -0.05 0.11 -2.04 15.40 0.17
25 3.02 0.15 0.51 0.00 -0.10 0.07 0.34 3.21 -17.92 0.35
26 -17.63 0.00 0.12 0.03 -0.05 0.03 -0.03 1.55 -39.72 -0.19
27 -1.54 -0.60 -0.54 0.01 0.13 -0.16 -0.27 -2.41 20.25 -0.18
28 1.56 -0.15 -0.22 0.03 -0.29 -0.20 0.19 -3.81 33.04 -0.01
29 1.56 -0.15 -0.22 0.03 -0.29 -0.20 0.19 -3.81 33.04 -0.01
30 5.65 -0.18 -0.28 0.02 -0.15 -0.19 0.12 -4.15 33.42 0.15
31 1.36 0.40 0.57 0.06 -0.24 0.11 0.18 2.04 0.87 0.17
32 -1.69 -1.01 -0.48 -0.01 0.01 0.01 -0.37 -1.61 -6.85 -0.20
33 -1.32 0.00 -0.27 0.01 0.07 0.02 -0.44 3.43 -20.62 -0.22
34 -1.60 -0.14 -0.44 0.02 0.09 -0.33 0.16 -0.34 15.63 -0.25
35 2.46 0.00 -0.06 0.10 0.10 -0.11 0.10 3.93 -16.98 -0.88
36 -3.28 -0.03 -0.27 -0.02 0.03 -0.25 -0.03 -1.54 -19.38 -0.30
37 -1.48 -0.17 -0.28 0.03 0.21 0.03 -0.22 -2.21 1.19 -0.12
38 -2.38 -0.24 -0.02 -0.01 0.17 0.09 -0.01 -2.26 -14.43 0.25
39 -0.50 -0.25 -0.24 0.02 0.04 0.19 0.10 2.10 -1.32 -0.21
40 -0.35 -0.44 -0.27 0.01 0.11 0.28 -0.26 -2.88 -3.57 0.66
41 0.78 0.32 0.22 0.02 -0.13 0.00 -0.71 0.42 -26.84 -0.88
42 -14.75 0.00 -0.11 0.04 0.17 -0.23 -0.10 0.30 -5.03 0.00
43 0.46 0.86 0.45 0.07 0.12 -0.09 0.37 5.16 -71.47 -0.12
44 13.48 -0.47 -0.17 -0.01 0.06 0.06 0.09 1.86 -1.98 -0.31
45 -1.83 -0.54 -0.44 0.01 -0.49 -0.21 -0.25 2.32 0.31 0.40
46 -2.56 0.25 0.31 0.00 -0.06 0.19 0.02 3.26 -6.25 -0.28
F0 F1 F2 F3 F4 F5 F6 F7 F8 F9
47 3.83 -0.48 -0.06 -0.02 -0.09 0.10 0.20 2.18 5.77 -0.40
48 -3.67 0.41 -0.03 0.03 0.14 -0.18 -0.04 -4.86 0.09 0.56
49 6.17 -0.26 0.09 -0.01 -0.05 0.13 0.16 1.81 44.54 0.00
50 4.11 0.45 0.35 0.00 0.24 -0.03 0.52 -0.47 -15.80 0.19
それぞれのグループに割り当てられた銘柄での F0 〜 F9 の合計が、与えられた範囲に収まること、また、特に重要な F2,F3,F4,F5,F6,F9 の属性については与えられた値にできるだけ近いことを要請します。
そのための数理計画モデルは以下のようになります。
> Basket
function(ng.d, fcoef.d, flow.d, fbar.d, fhigh.d, W.d)
{
M <- Set()
m <- Element(set=M)
J <- Set()
j <- Element(set=J)
K <- Set(1:ng.d)
k <- Element(set=K)
f <- Parameter(index=dprod(j,m), fcoef.d)
flow <- Parameter(index=m, flow.d)
fbar <- Parameter(index=m, fbar.d)
fhigh <- Parameter(index=m, fhigh.d)
W <- Parameter(index=m, W.d)
u <- IntegerVariable(index=dprod(j,k),type="binary")
.F <- Expression(index=dprod(m,k))
.F[m,k] ~ Sum(f[j,m]*u[j,k], j)
selection(u[j,k],k)
scf <- 1000
scf*fhigh[m] >= scf*.F[m,k] >= scf*flow[m]
softConstraint(scf,1)
W[m]*(.F[m,k] - fbar[m]) == 0
g <- Expression(index=j)
g[j] ~ Sum(u[j,k]*k, k)
}
各ファクターのグループごとの下限、
> Basket.flow F0 F1 F2 F3 F4 F5 F6 F7 F8 F9 -15 -15 -15 -15 -15 -15 -15 -15 -15 -15
上限、
> Basket.fhigh F0 F1 F2 F3 F4 F5 F6 F7 F8 F9 15 15 15 15 15 15 15 15 15 15
平均、
> Basket.fbar F0 F1 F2 F3 F4 F5 F6 F7 F8 F9 0 0 0 0 0 0 0 0 0 0
特に大事なものを定義する重みベクトル
> Basket.W F0 F1 F2 F3 F4 F5 F6 F7 F8 F9 0 0 1 1 1 1 1 0 0 1
が与えられているものとします。モデルを展開して求解します。グループの数は "5" としています。
> # 展開 > sys <- System(Basket, 5, Basket.fcoef, Basket.flow, Basket.fbar, Basket.fhigh, Basket.W) > # 実行 > nuopt.options(method = "wcsp",maxtim = 10) > sol <- solve(sys) > # 解の取得 > gnum <- as.array(current(sys,g))
得られた解について、グループごとのファクター値がどのようになっているか、適当にグルーピングした場合と比較します。fopt は上記のモデルを使ってグルーピングした場合、frand は適当にグルーピングした場合の結果です。
> # 検証
> fopt <- rbind(apply(Basket.fcoef[gnum == 1, ], 2, sum), apply(Basket.fcoef[
> gnum == 2, ], 2, sum), apply(Basket.fcoef[gnum == 3, ], 2, sum),
> apply(Basket.fcoef[gnum == 4, ], 2, sum), apply(Basket.fcoef[gnum ==
> 5, ], 2, sum))
> frand <- rbind(apply(Basket.fcoef[1:10, ], 2, sum), apply(Basket.fcoef[11:
> 20, ], 2, sum), apply(Basket.fcoef[21:30, ], 2, sum), apply(
> Basket.fcoef[31:40, ], 2, sum), apply(Basket.fcoef[41:50, ], 2,
> sum))
> fopt
F0 F1 F2 F3 F4 F5 F6 F7 F8 F9
[1,] -14.65 0.76 0.36 0.16 0.09 -0.23 -0.23 -4.52 13.94 0.59
[2,] -14.13 -0.13 -0.41 0.24 -0.56 -0.32 1.01 -0.57 14.82 -0.77
[3,] -14.72 -0.24 0.06 0.16 1.02 0.35 -0.07 -2.66 14.22 -0.38
[4,] -14.60 -0.30 0.69 0.03 0.41 0.62 -0.37 8.34 14.89 0.42
[5,] -14.76 -2.15 0.11 -0.06 -0.49 0.19 -0.57 9.92 13.50 -0.61
> frand
F0 F1 F2 F3 F4 F5 F6 F7 F8 F9
[1,] 17.79 0.48 1.53 0.07 0.24 0.80 0.36 5.07 -18.28 -0.14
[2,] -91.49 -0.57 -0.17 -0.08 0.47 0.14 -0.98 5.19 106.33 0.58
[3,] 3.60 -0.63 0.60 0.20 -0.74 -0.11 0.92 -12.39 125.44 0.75
[4,] -8.78 -1.88 -1.76 0.21 0.59 0.04 -0.79 0.66 -65.46 -1.10
[5,] 6.02 0.54 0.61 0.13 -0.09 -0.26 0.26 11.98 -76.66 -0.84
適当にグルーピングした場合と比べて、均質なグループが得られていることが分かります。
# 展開 sys <- System(Basket, 5, Basket.fcoef, Basket.flow, Basket.fbar, Basket.fhigh, Basket.W) # 実行 nuopt.options(method = "wcsp",maxtim = 10) sol <- solve(sys) # 解の取得 gnum <- as.array(current(sys,g)) # 検証 fopt <- rbind( apply(Basket.fcoef[gnum==1,],2,sum), apply(Basket.fcoef[gnum==2,],2,sum), apply(Basket.fcoef[gnum==3,],2,sum), apply(Basket.fcoef[gnum==4,],2,sum), apply(Basket.fcoef[gnum==5,],2,sum) ) frand <- rbind( apply(Basket.fcoef[1:10,],2,sum), apply(Basket.fcoef[11:20,],2,sum), apply(Basket.fcoef[21:30,],2,sum), apply(Basket.fcoef[31:40,],2,sum), apply(Basket.fcoef[41:50,],2,sum) ) fopt frand nuopt.options(nuopt.options(NA)) |
| ※S-PLUS や R でバッファに出力する設定をしている場合、ハングアップしているように見える場合がございますのでご注意下さい。 |
Maximum Drawdown
Maximum Drawdown はポートフォリオ最適化問題におけるリスク尺度の一つであり、ある特定の期間において、資産額の最も大きな減少量のことを指します。特徴として、分散や下方部分積率と異なり、入力データをサンプル点ではなく時系列として扱う、という点が挙げられます。ここでは所与の期間において Maximum Drawdown が最も小さくなるような各投資対象に対する投資量を決定します。なお、投資量は初期時点でのみ決定され、その後リバランスは行われないものとします。
まずは、収益率データを初期価格を 1 とした価格データに変換します。
> tmp <- rbind(rep(0,ncol(R.521x95)),R.521x95)+1 > P.521x95 <- apply(tmp,2,cumprod)
以下は Maximum Drawdown 最小化モデルです。変数 x は初期時点における各投資対象に対する投資量です。変数 W は時点 t における資産額、U は時点 t における W の累積最大値、V は時点 t における W の時間を逆方向に見た累積最小値です。Maximum Drawdown は U と V の差の最大値として表現することができます。
> MinMaxDD
function(P.d)
{
Asset <- Set()
j <- Element(set=Asset)
Period <- Set()
t <- Element(set=Period)
P <- Parameter(index=dprod(t,j),P.d)
x <- Variable(index=j)
U <- Variable(index=t)
V <- Variable(index=t)
W <- Variable(index=t)
maxdd <- Variable()
risk <- Objective()
risk ~ maxdd
W[t] == Sum(P[t,j]*x[j],j)
U[t,t>1] >= U[t-1,t>1]
U[t] >= W[t]
V[t-1,t>1] <= V[t,t>1]
V[t] <= W[t]
U[t] - V[t] <= maxdd
Sum(x[j],j) == 1
x[j] >= 0
}
システムを作成して解きます。
> sys <- System(MaxDD, P.521x95) > sol <- solve(sys)
解を取得します。
> x <- as.array(current(sys,x))
ここで、確認のために価格データと組入比率から Maximum Drawdown の値を求める関数を定義します。各期における資産額 W、累積最大値 U、時間を逆方向に見た累積最小値 V も求めます。
> calc.MaxDD
function(P,x)
{
W <- as.vector(P %*% as.vector(x))
U <- cummax(W)
V <- rev(cummin(rev(W)))
return(list(maxdd=max(U-V),W=W,U=U,V=V))
}
Maximum Drawdown を計算します。
> res <- calc.MaxDD(P.521x95,x) > res$maxdd
得られた解を棒グラフとして表示します。
> eps <- 1e-4 > barplot(x[x>=eps],names=names(x[x>=eps]))
 |
先ほど計算した W, U, V を図示します。
> m <- cbind(res$W,res$U,res$V)
> matplot(1:nrow(m),m,name=colnames(m),type="l")
> legend(0,1.08,c("W","U","V"),lty=1:3)
 |
# 価格データの作成
tmp <- rbind(rep(0,ncol(R.521x95)),R.521x95)+1
P.521x95 <- apply(tmp,2,cumprod)
# 展開
sys <- System(MinMaxDD, P.521x95)
# 実行
sol <- solve(sys)
# 解の取得
x <- as.array(current(sys,x))
# 確認
calc.MaxDD <- function(P,x)
{
W <- as.vector(P %*% as.vector(x))
U <- cummax(W)
V <- rev(cummin(rev(W)))
return(list(maxdd=max(U-V),W=W,U=U,V=V))
}
res <- calc.MaxDD(P.521x95,x)
res$maxdd
# 図示
eps <- 1e-4
barplot(x[x>=eps],names=names(x[x>=eps]))
m <- cbind(res$W,res$U,res$V)
matplot(1:nrow(m),m,name=colnames(m),type="l")
legend(0,1.08,c("W","U","V"),lty=1:3)
|
Sharpe Ratio
ポートフォリオ最適化問題の目的関数として Sharpe Ratio と呼ばれる指標を用いることがあります。Sharpe Ratio は
と表すことができます。ここで、
は投資対象の集合、
は組入比率、
は平均収益率、
は収益率の分散共分散行列、
は安全資産の収益率を表します。Sharpe Ratio 最大化問題は非線形計画問題（目的関数が非線形）となりますが S+NUOPT でそのまま解くことができます。
> Sharpe
function(Q.d,rb.d)
{
Asset <- Set()
j <- Element(set=Asset)
k <- Element(set=Asset)
Q <- Parameter(index=dprod(j,k),Q.d)
rb <- Parameter(index=j,rb.d)
x <- Variable(index=j)
x[j] ~ 0.1
f <- Objective(type="maximize")
f ~ (Sum(rb[j]*x[j],j)-0.002)/Sum(Q[j,k]*x[j]*x[k],j,k)^0.5
Sum(x[j],j) == 1
x[j] >= 0
}
早速解いてみましょう。まずは収益率データから収益率の分散共分散行列と平均収益率を求めます。
> Q <- var(R.60x200) > rb <- apply(R.60x200,2,mean) > rb <- as.array(rb)
システムを作成、求解を行い、解を得ます。
> sys <- System(Sharpe, Q, rb) > sol <- solve(sys) > x <- as.array(current(sys,x))
# 前処理 Q <- var(R.60x200) rb <- apply(R.60x200,2,mean) rb <- as.array(rb) # 展開 sys <- System(Sharpe, Q, rb) # 実行 sol <- solve(sys) # 解の取得 x <- as.array(current(sys,x)) # 図示 eps <- 1e-4 barplot(x[x>eps],names=names(x[x>eps])) |
実は Sharpe Ration 最大化問題は目的関数の分子を
とおいて、変数を
と変換することによって、等価な二次計画問題（目的関数が二次式）に置き換えることができます。目的関数を一般の非線形ではなく、二次式として表現できるということは、数理計画問題を解く上で大きなメリットとなります。
> Sharpe.qp
function(Q.d,rb.d)
{
Asset <- Set()
j <- Element(set=Asset)
k <- Element(set=Asset)
Q <- Parameter(index=dprod(j,k),Q.d)
rb <- Parameter(index=j,rb.d)
w <- Variable(index=j)
f <- Objective(type="minimize")
f ~ Sum(Q[j,k]*w[j]*w[k],j,k)
Sum((rb[j]-0.002)*w[j],j) == 1
w[j] >= 0
}
解いてみましょう。
> sys <- System(Sharpe.qp, Q, rb) > nuopt.options(eps=1e-10) > sol <- solve(sys) > w <- as.array(current(sys,w))
得られた解 w は元の問題の解 x に変換することができます。
> x.qp <- w/sum(w)
得られた解を比較すると同等の解が得られていることがわかります。
> eps <- 1e-4
> x[x>=eps]
A071 A080 A081 A087 A103 A113 A139
0.1317304 0.06223414 0.1040298 0.03068688 0.02559427 0.03757759 0.003409419
A189 A190 A195 A197 A200
0.235764 0.1985228 0.005258932 0.164983 0.0002087379
> x.qp[x.qp>=eps]
A071 A080 A081 A087 A103 A113 A139
0.1317538 0.06221308 0.1040643 0.03067362 0.02548632 0.03759601 0.003395785
A189 A190 A195 A197 A200
0.2357851 0.1985498 0.005192057 0.1649768 0.0003129249
# 前処理 Q <- var(R.60x200) rb <- apply(R.60x200,2,mean) rb <- as.array(rb) # 展開 sys <- System(Sharpe.qp, Q, rb) # 実行 nuopt.options(eps=1e-10) sol <- solve(sys) # 解の取得 w <- as.array(current(sys,w)) x <- w/sum(w) # 図示 eps <- 1e-4 barplot(x[x>=eps],names=names(x[x>=eps])) nuopt.options(nuopt.options(NA)) |
半正定値行列の取得
相関を持ついくつかのアセットの値動きをモンテカルロ法でシミュレーションする場合などには、アセット間の相関行列を与える必要があります。その際、相関行列は原理的にコレスキー分解が可能であるように正定値である必要があります。しかし、正定値であるという性質は直観的に明らかではありません。S+NUOPT に備わった半正定値計画アルゴリズムを用いれば、与えられた行列（正定値であるとは限らない）に最も近い正定値な行列を算出させることができます。
例えば次のすべて正の要素から成る 10×10の対称行列は正定値ではありませんのでコレスキー分解することができません。
> Cormat.A
X01 X02 X03 X04 X05 X06 X07 X08 X09 X10
X01 1.00 0.17 0.5 0.2 0.1 0.0 0.2 0.0 0.8 0.7
X02 0.17 1.00 0.1 0.9 0.6 0.4 0.0 0.0 0.0 0.0
X03 0.50 0.10 1.0 0.2 0.0 0.0 0.0 0.2 0.0 0.2
X04 0.20 0.90 0.2 1.0 0.2 0.2 0.2 0.2 0.0 0.0
X05 0.10 0.60 0.0 0.2 1.0 0.4 0.0 0.0 0.0 0.0
X06 0.00 0.40 0.0 0.2 0.4 1.0 0.0 0.0 0.0 0.0
X07 0.20 0.00 0.0 0.2 0.0 0.0 1.0 0.0 0.0 0.0
X08 0.00 0.00 0.2 0.2 0.0 0.0 0.0 1.0 0.0 0.0
X09 0.80 0.00 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
X10 0.70 0.00 0.2 0.0 0.0 0.0 0.0 0.0 0.0 1.0
これは、S-PLUS を使って求めた固有値：
> eigen(Cormat.A)$values [1] 2.6192649 2.0903778 1.2534419 1.0818527 1.0000000 0.836797 [7] 0.7030020 0.6042132 -0.0285787 -0.1603709
となり、最小固有値（-0.1603709）が負になることからわかります。次は与えられた行列とフロベニウスノルムの意味で近い、正定値行列を求めるための数理計画モデルです。
> Cormat
function(A, minEig = 0.001)
{
N <- Set()
i <- Element(set = N)
j <- Element(set = N)
A <- Parameter(A, index = dprod(i, j))
minEig <- Parameter(minEig)
X <- Variable(index = dprod(i, j))
M <- SymmetricMatrix(dprod(i, j))
M[i, j, i >= j] ~ X[i, j]
diffnrm <- Objective(type = minimize)
diffnrm ~ Sum((X[i, j] - A[i, j]) * (X[i, j] - A[i, j]), i, j, i >= j)
# 最小固有値が minEig 以上
M >= minEig
# 初期値
X[i, i] == 1
# 対称行列であることを要請
X[i,j,i<j] == X[j,i]
}
> s <- System(model=Cormat,Cormat.A)
> sol <- solve(s)
> Aopt <- as.array(current(s,X))
このようにして算出された Aopt は次のような行列です。各要素は元の行列とそれなりに近い値になっていることにご注意ください。
> round(Aopt,2)
X01 X02 X03 X04 X05 X06 X07 X08 X09 X10
X01 1.00 0.16 0.46 0.19 0.09 0.01 0.18 0.01 0.71 0.63
X02 0.16 1.00 0.10 0.88 0.59 0.40 0.01 0.00 0.01 0.00
X03 0.46 0.10 1.00 0.20 0.00 0.00 0.01 0.20 0.03 0.22
X04 0.19 0.88 0.20 1.00 0.21 0.20 0.20 0.20 0.00 0.00
X05 0.09 0.59 0.00 0.21 1.00 0.40 0.00 0.00 0.01 0.00
X06 0.01 0.40 0.00 0.20 0.40 1.00 0.00 0.00 0.00 0.00
X07 0.18 0.01 0.01 0.20 0.00 0.00 1.00 0.00 0.01 0.01
X08 0.01 0.00 0.20 0.20 0.00 0.00 0.00 1.00 -0.01 0.00
X09 0.71 0.01 0.03 0.00 0.01 0.00 0.01 -0.01 1.00 0.05
X10 0.63 0.00 0.22 0.00 0.00 0.00 0.01 0.00 0.05 1.00
これは、Cormat.A にフロベニウスノルムの意味で最も近い、最小固有値が 0.001 以上の行列です。
> eigen(Aopt)$values [1] 2.580796251 2.028567070 1.236807350 1.067534416 0.968329999 [6] 0.830140682 0.685764222 0.600045416 0.001012264 0.001002330
となることから、Aopt の最小固有値が 0.001 以上であることが確認できます。
# 確認 Cormat.A eigen(Cormat.A)$values # 展開 sys <- System(Cormat, Cormat.A) # 実行 sol <- solve(sys) # 解の取得 X <- as.array(current(sys,X)) # 確認 X eigen(X)$values |
ロバスト最適化
ポートフォリオのリスクを投資対象の収益率の分散共分散行列を用いて計測するマルコビッツモデルにおいて、与えられた分散共分散行列は、しばしば不確実性を伴います。この不確実性に対してロバストな解を得る以下の問題を考えます。
以下の典型的な平均・分散モデルを考えます。ここでは期待収益率とリスクの線形結合（効用関数）を目的関数としています。
ここで、
を期待リターン、
をリターンの分散共分散行列、
をポートフォリオの重み、
をリスク回避係数とします。
分散共分散行列
に不確実性が伴うとして、以下のロバストポートフォリオ最適化問題を考えます。なお
は、不確実性が伴うことによって取りうる
に関する行列の集合とします。
として、分散共分散行列
の各要素が
のような区間を持つとします。このとき(2)は、新たに行列
を用いて以下のような問題に置き換えることができます。なお行列A、Bに対し、
は、AとBの要素ごとの積の和を表すものとします。
次の S+NUOPT のスクリプトは、分散共分散行列
の各要素の下限を表す行列sigL、上限を表す行列sigU、収益率の期待値mu、およびlambdaが与えられているとき、上記の問題を解くものです。
> Robust
function(sigL, sigU, mu, lambda)
{
Asset <- Set()
i <- Element(set = Asset)
j <- Element(set = Asset)
sigL <- Parameter(sigL, index = dprod(i, j))
sigU <- Parameter(sigU, index = dprod(i, j))
lambda <- Parameter(lambda)
mu <- Parameter(mu, index = i)
U <- Variable(index = dprod(i, j))
L <- Variable(index = dprod(i, j))
x <- Variable(index = i)
V <- Set(c(as.list(Asset),"dummy")) # Asset に "dummy" を加えた集合 V を作る
v <- Element(set = V)
w <- Element(set = V)
M <- SymmetricMatrix(dprod(v, w))
M[i, j] ~ U[i, j] - L[i, j]
M[j, "dummy"] ~ x[j]
M["dummy", "dummy"] ~ 1
f <- Objective(type = maximize)
f ~ Sum(mu[i] * x[i], i) - lambda * Sum(sigU[i, j] * U[i, j]
- sigL[i, j] * L[i, j], i, j)
M >= 0
Sum(x[i], i) == 1
U[i, j, i > j] == U[j, i]
L[i, j, i > j] == L[j, i]
U[i, j] >= 0
L[i, j] >= 0
}
データは以下のように与えられているとします。
> Robust.sigU # 分散共分散の要素の値の上限
X1 X2 X3 X4 X5 X6 X7 X8
X1 3.0 1.0 2.5 -0.9 1.00 -0.4 2.50 2.0
X2 1.0 4.5 -1.4 1.2 0.50 -1.1 0.50 2.6
X3 2.5 -1.4 6.0 -0.5 1.20 1.0 0.40 -0.1
X4 -1.0 1.2 -0.5 6.5 1.60 0.5 1.30 -0.8
X5 1.0 0.5 1.2 1.6 5.50 -1.3 0.16 -1.9
X6 -0.4 -1.1 1.0 0.5 -1.30 11.0 -0.90 0.6
X7 2.5 0.5 0.4 1.3 0.16 -0.9 6.50 -1.2
X8 2.0 2.6 -0.1 -0.8 -1.90 0.6 -1.20 13.5
> Robust.sigL # 分散共分散の要素の値の下限
X1 X2 X3 X4 X5 X6 X7 X8
X1 1.9 -1.4 -1.000 -1.000 1.00 -0.5 -1.00 0.10
X2 -1.4 3.5 -1.500 0.500 -1.20 -1.2 0.40 0.60
X3 -1.0 -1.5 5.000 -1.125 1.10 0.9 0.30 -0.20
X4 -1.0 0.5 -1.125 6.000 1.50 0.4 1.20 -0.90
X5 1.0 -1.2 1.100 1.500 4.50 -1.4 0.15 -2.00
X6 -0.5 -1.2 0.900 0.400 -1.40 9.0 -1.00 0.50
X7 -1.0 0.4 0.300 1.200 0.15 -1.0 5.50 -1.25
X8 0.1 0.6 -0.200 -0.900 -2.00 0.5 -1.25 11.50
> Robust.mu # 収益率の期待値
X1 X2 X3 X4 X5 X6 X7 X8
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
次のようにして解きます。
> s <- System(model=Robust,Robust.sigL,Robust.sigU,Robust.mu,1.0) > # 求解 > sol <- solve(s) > # 図示 > xopt <- as.array(current(s,x)) > barplot(names=rownames(xopt),xopt)
次のようなポートフォリオが得られました。
 |
# 展開 sys <- System(Robust, Robust.sigL, Robust.sigU, Robust.mu, 1) # 実行 sol <- solve(sys) # 解の取得 xopt <- as.array(current(sys,x)) barplot(names=rownames(xopt),xopt) |
イールドカーブのフィッティング
非線形なモデル関数をデータに合わせこむ問題は次の形の非線形な数理計画問題として一般的に記述することができます。
変数：
（モデルパラメータ）
最小化：
（誤差の二乗和）
制約：
,
｛観測点｝（誤差の定義）
（パラメータについての制約）
ここで、
はパラメータに対して線形あるいは非線形な式として定義されているモデル式（
：パラメータ、
：変数）で、
個の観測点
（
｛観測点｝）において、観測値
が与えられているとします。
はパラメータ値に関する制限を記述した制約式です。
重回帰はこの最も基本的な形で、
が線形関数である場合に対応します。
が線形関数でもパラメータの値に制限（制約）がある場合にもこの数理計画問題を解く必要があります。
金融工学に表れるフィッティングで最も多く現れるのは、スポットレート（割引債の利回り）を、観測された利付債の現在価格から推定するイールドカーブフィッティングです。
償還期間
における額面価格 100、クーポンレート 1% の利付債の理論現在価格
は、スポットーレート
を用いて以下のように表すことができます。
償還期間
の利付債の現在価格
が以下のように観測されたとします。
> Yield.term
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
5 5 6 6 6 6 6 6 7 7 7 7 7 7 8 8 8 8 8 8 9 9 9 9 9
54 55 56 57 58 59 60
9 10 10 10 10 10 10
> Yield.price
1 2 3 4 5 6 7 8 9 10 11
100.39 102.15 99.24 101.32 99.72 101.51 100.62 99.34 98.27 98.31 100.97
12 13 14 15 16 17 18 19 20 21 22 23 24
101.73 99.78 100.47 98.19 99.55 99.79 98.4 96.6 96.93 96.8 94.91 96.28 95.33
25 26 27 28 29 30 31 32 33 34 35 36 37
95.13 93.5 93.42 91.56 92.67 96.28 89.66 89.46 91.21 91.42 93.32 89.76 90.18
38 39 40 41 42 43 44 45 46 47 48 49 50
88.58 87.41 90.33 87.5 87.66 86.06 85.55 84.74 88.78 86.79 88.76 84.95 84.31
51 52 53 54 55 56 57 58 59 60
87.24 84.73 83.76 84.02 81.75 84.46 82.51 85.42 81.45 85.36
次の数理計画モデルは各観測点
における理論現在価格
と観測された現在価格
の差の二乗和：
を最小にするスポットーレート
を求めます。
> Yield
function(telem.d, tvalue.d, Svalue.d)
{
Term <- Set()
Point <- Set()
t <- Element(set=Term)
i <- Element(set=Point)
r <- Variable(index=t)
telem <- Parameter(index=t,telem.d)
tvalue <- Parameter(index=i,tvalue.d)
Svalue <- Parameter(index=i,Svalue.d)
d <- Expression(index=t)
d[t] ~ 1 / (1+0.01*r[t])^telem[t]
S <- Expression(index=i)
S[i] ~ Sum(d[t],t,t<=tvalue[i]) + Sum(100/(1+0.01*r[t])^telem[t],
t,t==tvalue[i])
diff <- Expression(index=i)
diff[i] ~ S[i] - Svalue[i]
err <- Objective()
err ~ Sum(diff[i]*diff[i], i)
0 <= r[t]
}
次はデータを与えて解きます。データは各債券の償還期間（tvalue.d）および市場で観測された価格（Svalue）となります。telem.dは
の各項のべきの値のベクトルです。
解いてみます。
> # 展開 > sys <- System(Yield, Yield.telem, Yield.term, Yield.price) > # 実行 > sol <- solve(sys) > # 解の取得 > r <- as.array(current(sys,r)) > # 図示 > plot(r,type="b")
以下のイールドカーブが得られます。
 |
# 展開 sys <- System(Yield, Yield.telem, Yield.term, Yield.price) # 実行 sol <- solve(sys) # 解の取得 r <- as.array(current(sys,r)) S <- as.array(current(sys,S)) # 図示 plot(r,type="b") plot(Yield.price,S,xlab="Observed",ylab="Theoretical") |
格付け推移行列の推定
格付け（Rating）とは、企業の発行する社債の元本、利息の支払い能力をランク形式で表示したものです。格付け会社は独自の調査結果のもと、A、B、C などの記号を用いて対象社債のリスク度合いを示します。格付け推移行列とは、ある格付け評価を受けている企業が、一定期間後にどのような格付けとなるかについての確率を表す行列のことをいいます。
ここでは、1 年分の格付け推移行列から 1 ヶ月分の格付け推移行列を推定する問題を扱います。
> Rating
function(q0.d)
{
Rating <- Set()
i <- Element(set=Rating)
j <- Element(set=Rating)
k <- Element(set=Rating)
q0 <- Parameter(index=dprod(i,j),q0.d)
q <- Variable(index=dprod(i,j))
q2 <- Expression(index=dprod(i,j))
q4 <- Expression(index=dprod(i,j))
q8 <- Expression(index=dprod(i,j))
q12 <- Expression(index=dprod(i,j))
q2[i,j] ~ Sum(q[i,k]*q[k,j], k)
q4[i,j] ~ Sum(q2[i,k]*q2[k,j], k)
q8[i,j] ~ Sum(q4[i,k]*q4[k,j], k)
q12[i,j] ~ Sum(q8[i,k]*q4[k,j], k)
diff <- Expression(index=dprod(i,j))
diff[i,j] ~ q0[i,j] - q12[i,j]
diffnrm <- Objective()
diffnrm ~ Sum(diff[i,j]^2, i, j)
Sum(q[i,j], j) == 1
0 <= q[i,i] <= 1
0 <= q[i,j, i!=j] <= 0.05
q[i,i] ~ 0.9
}
> # 展開
> sys <- System(Rating, Rating.Q0)
> # 実行
> sol <- solve(sys)
> # 解の取得
> q <- as.array(current(sys,q))
> # 確認
> q12 <- diag(1,nrow(q))
> dimnames(q12) <- dimnames(q)
> for(i in 1:12) q12 <- q12 %*% q
> # 表示
> r.order <- order(rownames(Rating.Q0)) # 格付けの順番を定義
得られた 1 ヶ月分の格付け推移行列を表示します。
> round(q[r.order,r.order],digits=4)
AAA AA A BBB BB B CCC CC C
AAA 0.9970 0.0030 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
AA 0.0031 0.9946 0.0023 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
A 0.0001 0.0038 0.9945 0.0016 0.0000 0.0000 0.0000 0.0000 0.0000
BBB 0.0000 0.0000 0.0031 0.9952 0.0014 0.0003 0.0000 0.0000 0.0000
BB 0.0000 0.0000 0.0000 0.0098 0.9885 0.0004 0.0013 0.0000 0.0000
B 0.0000 0.0000 0.0000 0.0000 0.0069 0.9889 0.0036 0.0006 0.0000
CCC 0.0000 0.0000 0.0000 0.0000 0.0000 0.0025 0.9970 0.0005 0.0000
CC 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0007 0.9972 0.0021
C 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0023 0.9977
1 ヶ月分の格付け推移行列を 12 乗、つまり 1 年分に変換した値を示します。
> round(q12[r.order,r.order],digits=4)
AAA AA A BBB BB B CCC CC C
AAA 0.9649 0.0346 0.0004 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
AA 0.0355 0.9381 0.0261 0.0002 0.0000 0.0000 0.0000 0.0000 0.0000
A 0.0014 0.0432 0.9364 0.0186 0.0003 0.0001 0.0000 0.0000 0.0000
BBB 0.0000 0.0007 0.0348 0.9454 0.0152 0.0036 0.0002 0.0000 0.0000
BB 0.0000 0.0000 0.0018 0.1071 0.8716 0.0045 0.0149 0.0001 0.0000
B 0.0000 0.0000 0.0000 0.0041 0.0734 0.8752 0.0400 0.0071 0.0001
CCC 0.0000 0.0000 0.0000 0.0000 0.0010 0.0274 0.9650 0.0064 0.0001
CC 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0082 0.9674 0.0241
C 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0265 0.9733
与えられた 1 年分の格付け推移行列を表示します。
> round(Rating.Q0,digits=4)
AAA AA A BBB BB B CCC CC C
AAA 0.9651 0.0349 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
AA 0.0356 0.9382 0.0262 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
A 0.0014 0.0433 0.9364 0.0186 0.0003 0.0000 0.0000 0.0000 0.0000
BBB 0.0000 0.0000 0.0352 0.9456 0.0154 0.0038 0.0000 0.0000 0.0000
BB 0.0000 0.0000 0.0000 0.1078 0.8720 0.0049 0.0153 0.0000 0.0000
B 0.0000 0.0000 0.0000 0.0000 0.0747 0.8762 0.0410 0.0081 0.0000
CCC 0.0000 0.0000 0.0000 0.0000 0.0000 0.0278 0.9654 0.0068 0.0000
CC 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0083 0.9675 0.0242
C 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0266 0.9734
与えられた1 年分の格付け推移行列と推定した1 ヶ月分の格付け推移行列を 1 年分に変換した行列を見比べ、1 ヶ月分の格付け推移行列が適切に推定できていることを確かめることができます。
# 展開 sys <- System(Rating, Rating.Q0) # 実行 sol <- solve(sys) # 解の取得 q <- as.array(current(sys,q)) # 確認 q12 <- diag(1,nrow(q)) dimnames(q12) <- dimnames(q) for(i in 1:12) q12 <- q12 %*% q r.order <- order(rownames(Rating.Q0)) round(q[r.order,r.order],digits=4) round(q12[r.order,r.order],digits=4) round(Rating.Q0,digits=4) |
ロジスティック回帰
S+NUOPT を用いてロジスティック回帰を行うモデルを作成することもできます。
> LogReg
function(X.d,t.d)
{
M <- Set()
i <- Element(set=M)
L <- Set()
l <- Element(set = L)
X <- Parameter(index=dprod(l,i),X.d)
t <- Parameter(index=l,t.d)
a <- Variable(index=i)
a0 <- Variable()
y <- Expression(index=l)
y[l] ~ exp(Sum(a[i]*X[l,i],i)+a0) / (1+exp(Sum(a[i]*X[l,i],i)+a0))
mle <- Objective(type="maximize")
mle ~ Sum(t[l]*log(y[l])+(1-t[l])*log(1-y[l]),l)
}
ここでは、例として Iris virsinica と Iris versicolor の判別分析を扱います。
> # 展開 > sys <- System(LogReg,LogReg.X,LogReg.t) > # 実行 > sol <- solve(sys) > # 解の取得 > a <- as.array(current(sys,a)) > a0 <- as.array(current(sys,a0))
得られた解について、学習データとは別のテストデータで予測力を評価します。
> # 予測
> eta <- a0 + LogReg.test.X %*% as.vector(a)
> p <- exp(eta)/(1+exp(eta))
> res <- cbind(round(LogReg.test.t,0),round(p,0),eta)
> dimnames(res)[[2]] <- c("answer","predict","eta")
> res
S+NUOPT には、数理計画モデルを自在に設計できるというメリットがあります。例えば、推定したいパラメータに何らかの制約を課して回帰を行ったり、ロジット関数を 2 次式で表して予測力向上のために 2 次の係数行列に半正定値制約を課したりすることができます。
# 展開
sys <- System(LogReg,LogReg.X,LogReg.t)
# 実行
sol <- solve(sys)
# 解の取得
a <- as.array(current(sys,a))
a0 <- as.array(current(sys,a0))
# 予測
eta <- a0 + LogReg.test.X %*% as.vector(a)
p <- exp(eta)/(1+exp(eta))
res <- cbind(round(LogReg.test.t,0),round(p,0),eta)
dimnames(res)[[2]] <- c("answer","predict","eta")
res
|