トップ > 製品概要 > 派生製品 > SNUOPT > SNUOPT の例題

SNUOPTの例題

SNUOPT の例題をご紹介します。以下の例題は RNUOPTでも同様に解くことができますが、出力結果及びグラフに関しては RNUOPT のものとは異なる部分があります。

例題の動かし方

目次

基本的なマルコビッツモデル

4 銘柄の時系列データを用いて、マルコビッツモデル:

最小化(収益率の分散)

制約

(最低期待収益率の確保)

, (空売りの禁止)

による最適ポートフォリオを求めます。ここではを1% とします。

元にするのは4銘柄の60期間分の収益率データ R.60x4 です(SNUOPTにサンプルとして含まれています)。S-PLUSでグラフ化してみましょう。

> matplot(rownames(R.60x4),R.60x4,pch=colnames(R.60x4),type="o")

4銘柄の60期間分の収益率データのグラフ

次に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

Numerical Optimizer が認識するためには、S-PLUSのvector型オブジェクトはarray型として定義しておく必要があります。

> rbar <- as.array(rbar)

最適化モデルはS-PLUSの手続きとして定義しておきます(SNUOPTにサンプルとして含まれています)。

> 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ケースの5銘柄の収益率を含むR.8000x5です。

まずは分散最小化モデル(コンパクト分解表現)を解いてみましょう。SNUOPTではモデルを次のように定義します(SNUOPTにサンプルとして含まれています)。

> 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)

1次の下方部分積率(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)

1次の下方部分積率(LPM)をリスクとしたポートフォリオ最適解のヒストグラム

目標収益率 -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)

CVaRをリスクとしたポートフォリオ最適解のヒストグラム

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銘柄のポートフォリオではこのような分散・共分散行列を陽に生成すると計算負荷が増大してしまいます。

そのような場合には

:各銘柄の収益率のヒストリカルデータ(各列jは、期間にわたる銘柄jの収益率の観測値を並べたもの)

を直接データとして入力、

と表されることを利用して、補助変数を定義して次のような制約を導入し、

リスクを次のように書き換えます。

こうすれば、分散・共分散行列を直接生成しなくともよいので、効率的に計算を行うことができます。 1000銘柄のポートフォリオを求めてみます。データとして月次5年分(60期間)のデータ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
…

次が切り上げ、切り下げを決定する最適化モデルです。比較的規模の大きな0−1計画問題ですのでメタヒューリスティクス解法 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)

理想的な取引額(ideal)と実際の取引額(real)の棒グラフ

<< 実行例 >>
# 展開(総額 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(MinMaxDD, 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]))

Maximum Drawdown 最小化モデルの最適解

先ほど計算した 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)

各期における資産額 W,累積最大値 U,時間を逆方向に見た累積最小値 V

<< 実行例 >>
# 価格データの作成
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 最大化問題は非線形計画問題(目的関数が非線形)となりますが SNUOPT でそのまま解くことができます。

> 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))

半正定値行列の取得

相関を持ついくつかのアセットの値動きをモンテカルロ法でシミュレーションする場合などには、アセット間の相関行列を与える必要があります。その際、相関行列は原理的にコレスキー分解が可能であるように正定値である必要があります。しかし、正定値であるという性質は直観的に明らかではありません。SNUOPT に備わった半正定値計画アルゴリズムを用いれば、与えられた行列(正定値であるとは限らない)に最も近い正定値な行列を算出させることができます。

例えば次のすべて正の要素から成る 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の要素ごとの積の和を表すものとします。

次の SNUOPT のスクリプトは、分散共分散行列の各要素の下限を表す行列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)

ロジスティック回帰

SNUOPT を用いてロジスティック回帰を行うモデルを作成することもできます。

> 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

SNUOPT には、数理計画モデルを自在に設計できるというメリットがあります。例えば、推定したいパラメータに何らかの制約を課して回帰を行ったり、ロジット関数を 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