ページ

2016年6月19日日曜日

【R】 成長(ロジスティック)曲線とマルコフ連鎖(2)

前回は、ロジスティック曲線ないしロトカ・ヴォルテラ・モデルとマルコフ連鎖との関係を整理してみた。要点は、モデルが非斉次(非線形性)という本質的な特徴を持っている、ということであった。

今回は斉次(線形)なマルコフ連鎖による銘柄転移モデルがいかにうまくないか、というあたりを説明することにしよう。



銘柄転移モデルへの応用例の振り返り

「マルコフ解析」(森村・高橋1995)で銘柄転移モデルが紹介されている。そのうちの一例 Theil and Rey(1966)は斉次マルコフ連鎖によってはうまく説明できていないように見える(斉次マルコフ連鎖について、Ehrenberg(1965)が痛烈な批判をし、それ以降、この手の論文がほとんど出なくなった、と言われている)。

Theil,H. and G.Rey , A quadratic programming approach to the estimation of transition probabilities, Management Science, 12(1966), pp.714-721 (あのTheil先生か?)
Ehrenberg, A.S.C. , An appraisal of Markov brand-switching models, Journal of Marketing Research, 2 (1965), pp.347-362

今回使うデータは、森村・高橋に掲載されている図から読み取りを行って作る。ちなみにグラフ数値の読み取りはフリーソフトGSYSを使わせていただいた。

GSYS 2.4

Theil and Reyの方法による遷移確率の推計

まずは、森村・高橋で紹介されているTheil and Reyの方法をまとめておこう。

t期について各銘柄の市場シェアデータを、

\[
\boldsymbol{y}_{t} = \left( \begin{array}{c} y_{1t} \\ \vdots \\ y_{nt} \end{array} \right) , \; t=0,1,\cdots,T
\]

と表現する。これが次のPを固定的な遷移確率行列とする斉次マルコフ連鎖に従っているものと考える。

\[
\boldsymbol{y}_{t} = \boldsymbol{P} \boldsymbol{y}_{t-1}
\]

ここでの問題はy0・・・yTからPを推定することである。そのために次のような制約付の最小自乗法を適用してみる。
\[
\begin{eqnarray} & & \frac{1}{2} \sum_{t=1}^{T} \sum_{i=1}^{n} \left( y_{it} - \sum_{k=1}^{n}p_{ik} y_{k,t-1} \right)^{2} \to \min \\ & \mbox{s.t.} & \sum_{i=1}^{n} p_{ij} = 1 , \; p_{ij \in J} = 0 \end{eqnarray}
\]

ここで最初の制約の意味は遷移確率行列であることを示している。2番目は遷移確率を強制的にゼロとするパスi → jの組み合わせの集合をJとしているのである。何故、強制的にゼロにする必要があるかというと、推定の結果、確率がマイナスで(あるいは小さすぎる値で)出てくるものをゼロとみなしたり、経験的に最初からそのパスの可能性をないものとすることがあるからである。(なお、制約条件の実際の取り扱いについて、最後の方に注釈を入れてあるので、参照されたい。)

制約付最小化問題なので、ラグランジェ乗数法で解く。
\[
\begin{eqnarray} {\cal L} & = & \frac{1}{2} \sum_{t=1}^{T} \sum_{i=1}^{n} \left( y_{it} - \sum_{k=1}^{n}p_{ik} y_{k,t-1} \right)^{2} \\ & & + \sum_{j=1}^{n} \lambda_{j} \left( \sum_{i=1}^{n} p_{ij} - 1 \right) + \sum_{(i,j) \in J} \mu_{ij} p_{ij} \end{eqnarray}
\]

Lをpijで偏微分することにより、次のような正規方程式を得る。

\[
\begin{eqnarray} & & \sum_{t=1}^{T} y_{it} y_{j,t-1} - \sum_{k=1}^{n} \left( \sum_{t=1}^{T} y_{k,t-1} y_{j,t-1} \right) p_{ik} \\ & & = \left\{ \begin{array}{l l} \lambda_{j} & (i,j) \not \in J \\ \lambda_{j} + \mu_{ij} & (i,j) \in J \end{array} \right. \end{eqnarray}
\]

これを行列表記すると、次のようになる(eはすべての要素が1であるベクトル)。
\[
\boldsymbol{W} - \boldsymbol{P} \boldsymbol{Z} = \boldsymbol{e}\boldsymbol{\lambda}^{\prime} + \boldsymbol{M}
\]
ただし、
\[
\boldsymbol{W} = \boldsymbol{Y}_{1}^{\prime} \boldsymbol{Y}_{0} , \; \boldsymbol{Z} = \boldsymbol{Y}_{0}^{\prime} \boldsymbol{Y}_{0} , \; \boldsymbol{Y}_{0} = \left( \begin{array}{c c c} y_{10} & \cdots & y_{n0} \\ \vdots & & \vdots \\ y_{1,T-1} & \cdots & y_{n,T-1} \end{array} \right)
\]
\[
\boldsymbol{\lambda} = \left( \begin{array}{c} \lambda_{1} \\ \vdots \\ \lambda_{n} \end{array} \right) , \; \boldsymbol{M} = \left( m_{ij} \right), \; m_{ij} = \left\{ \begin{array}{c c} 0 & (i,j) \not \in J \\ \mu_{ij} & (i,j) \in J \end{array}\right.
\]

正規方程式からλを消去する。そのためにe'を左から掛ける。

\[
\boldsymbol{e}^{\prime}(\boldsymbol{W} - \boldsymbol{P} \boldsymbol{Z}) = \boldsymbol{e}^{\prime} \boldsymbol{e}\boldsymbol{\lambda}^{\prime} + \boldsymbol{e}^{\prime} \boldsymbol{M}
\]
変形すると(左辺はゼロベクトルになる)、
\[
\boldsymbol{0} = n \boldsymbol{\lambda}^{\prime} + \boldsymbol{e}^{\prime} \boldsymbol{M} \Longrightarrow \boldsymbol{\lambda}^{\prime} = - \frac{1}{n} \boldsymbol{e}^{\prime} \boldsymbol{M}
\]
λを正規方程式に代入して、Pについて解くと、
\[
\boldsymbol{P} = \left[\boldsymbol{W} + \left(\frac{\boldsymbol{e}\boldsymbol{e}^{\prime}}{n} - \boldsymbol{I} \right) \boldsymbol{M} \right] \boldsymbol{Z}^{-1}
\]

実際の推定作業の手順は次のとおりである。
  1. 最初に制約なしでPを計算してみる。
  2. 解を見て、ゼロ制約を置くべき集合Jを確定する。
  3. Jの要素の数だけの連立方程式を解き、Mを確定する。
  4. Mによって初期解を補正する。

実際にRで計算してみる

まずデータを読み込み、プロットしてみる。
filename <- "tabako.csv"
tabako <- read.table(filename,header=T,sep=",")
tabako[,2:4] <- tabako[,2:4]/rowSums(tabako[,2:4])

plot(tabako[,c(1,2)],type="l",ylim=c(0,0.7), ann=F, lty=1)
par(new=T)
plot(tabako[,c(1,3)],type="l",ylim=c(0,0.7),ann=F,lty=2) 
par(new=T)
plot(tabako[,c(1,4)],type="l",ylim=c(0,0.7),lty=3,ylab="prob") 
legend(1930, 0.15, c("Camel","Chesterfield","Lucky.Strike"),lty=c(1,2,3), merge = TRUE, bg='gray90') 
3行目は図読み取りソフトの補正を行っている。図は冒頭に示したものであり、前回の図がそのまま読み取れていることが確認できる。 次にデータ行列Yを用意する。
Y <- data.matrix(tabako[,2:4])
Y0 <- Y[1:(nrow(Y)-1),]
Y1 <- Y[2:nrow(Y),]
次に制約なしで正規方程式を解いて、遷移確率行列を求めてみる。
W <- t(Y1) %*% Y0
ZZ <- solve(Z)
P <- W %*% ZZ
結果は次のようになる。マイナスの値が出てしまっているので、補正が必要になる。
> P
                  Camel Chesterfield Lucky.Strike
Camel        0.56438232    0.5990459  -0.05429253
Chesterfield 0.08967302    0.7679806   0.09773020
Lucky.Strike 0.34594466   -0.3670265   0.95656233
ここでMとして次のような行列を想定する。 \[ \boldsymbol{M} = \left( \begin{array}{c c c} 0 & 0 & \mu_{13} \\ 0 & 0 & 0 \\ 0 & \mu_{32} & 0 \end{array} \right) \] すると、 \[ \left(\frac{\boldsymbol{e}\boldsymbol{e}^{\prime}}{n} - \boldsymbol{I} \right) \boldsymbol{M} = \frac{1}{3} \left( \begin{array}{c c c} 0 & \mu_{32} & -2 \mu_{13} \\ 0 & \mu_{32} & \mu_{13} \\ 0 & -2 \mu_{32} & \mu_{13} \end{array} \right) \] これから連立方程式を作ると(ZijはZの逆行列の要素)、 \[ \frac{1}{3} \left( \begin{array}{c c} -2z_{22} & z_{32} \\ z_{23} & -2z_{33} \end{array} \right) \left( \begin{array}{c} \mu_{32} \\ \mu_{13} \end{array} \right) = \left( \begin{array}{c} \Delta p_{32} \\ \Delta p_{13} \end{array} \right) \] ゆえに、 \[ \left( \begin{array}{c} \mu_{32} \\ \mu_{13} \end{array} \right) = 3 \left( \begin{array}{c c} -2z_{22} & z_{32} \\ z_{23} & -2z_{33} \end{array} \right)^{-1} \left( \begin{array}{c} \Delta p_{32} \\ \Delta p_{13} \end{array} \right) \] 具体的に計算してみる。
M <- P
M[P<0] <- - 3 * solve(matrix(c(-2*ZZ[2,2],ZZ[3,2],ZZ[2,3],-2*ZZ[3,3]),nrow=2)) %*% P[P<0]
M[P>0] <- 0
P2 <- (W + ((diag(0,3) +1)/3 - diag(3)) %*% M) %*% ZZ
結果は次のようになる。
> P2
                 Camel Chesterfield Lucky.Strike
Camel        0.6754543 3.910161e-01 1.776357e-15
Chesterfield 0.1857151 6.089839e-01 1.258319e-01
Lucky.Strike 0.1388306 3.552714e-15 8.741681e-01
> P
                  Camel Chesterfield Lucky.Strike
Camel        0.56438232    0.5990459  -0.05429253
Chesterfield 0.08967302    0.7679806   0.09773020
Lucky.Strike 0.34594466   -0.3670265   0.95656233
補正すると結構数字が変わる。 遷移確率が推定できたので、次は予測を行ってみる。
fct0 <- Y[1,]
fctbl <- fct0
for(i in 1:18){
 fct0 <- fct0 %*% t(P2)
 fctbl <- rbind(fctbl,fct0)
}
fctbl <- cbind(tabako[,1],fctbl)
結果をプロットすると次のような図が得られる。明らかに、データと比べて動きが静かすぎる。

 

 

この食い違いは遷移確率の補正のせいではなく、斉次であるという条件に由来するものと考えられる。どのように遷移確率の計算を精緻にしたとしても、この条件が変わらない限り、不自然な結果しか得られない。

 

まとめ
苦心して計算したわりには、適当な結果が得られていないのだった。このことは予想はできていたが、むしろ想定以上に計算に苦労させられた。

 

次回は、いよいよ非斉次モデルを使った分析に進む。

 

制約条件に関する注釈
制約について一言。
  1. 最初に行う推計では、まったく制約をかけなくて良い。その結果として不具合のある(i,j)について集合Jの中に入れて、再推計を行う。
  2. 最初の推計では遷移確率に関する第一の制約も不要である。なぜならば、用いるデータが市場シェアであるから。

2番目の事柄は次のことからわかる。まずすべての要素が1であるベクトルeをマルコフ連鎖の式の左から掛ける。
\[
\boldsymbol{e}^{\prime} \boldsymbol{y}_{t} = \boldsymbol{e}^{\prime} \boldsymbol{P} \boldsymbol{y}_{t-1}
\]

変形すると、
\[
1 = (\boldsymbol{e}^{\prime} + \boldsymbol{a}^{\prime}) \boldsymbol{y}_{t-1}
\]
(なお、ここでは次の仮定をおいている。)
\[
\boldsymbol{e}^{\prime} \boldsymbol{P} = \boldsymbol{e}^{\prime} + \boldsymbol{a}^{\prime}
\]

さらに変形すると、
\[
1 = 1 + \boldsymbol{a}^{\prime} \boldsymbol{y}_{t-1} \Longrightarrow \boldsymbol{a} = \boldsymbol{0} \Longrightarrow \boldsymbol{e}^{\prime} \boldsymbol{P} = \boldsymbol{e}^{\prime}
\]
だから一番目の制約は不要。