遷移行列の冪乗根

状態遷移確率行列を２乗することで，２時点先の状態遷移確率行列を得ます．ですから，１年間の状態遷移確率行列の12乗根を求めれば，１ヶ月間の状態遷移確率行列になります．
ですが，SASの行列演算言語IMLでは，行列の冪乗演算子**は，-1以上の整数しか受付けません．

そこで，以下に，行列を固有値と固有ベクトル行列に分解して，冪乗根を求める方法を示します．
ここでは，平方根を求めていますが，##(1/2)を，##(1/12)にすれば12乗根が得られます．
（注：##は，IMLでの，行列の「各要素」の冪乗演算子です．）

警告！：すみません．固有値が複素数になるケースに対応していませんでした．
そのうち，直しますが，どなたか直してくれるとうれしいです（翔）
→たねやさんが回答してくれました．（感謝感謝）

*************************************************
rootmatrix.sas
翔
*************************************************;

proc iml;

A={
 0.5 0.3  0.2
,0.1 0.85 0.05
, 0    0    1
};
print A;


/*
固有値は実数であるとみなしてますが，
複素数の場合もあり，その計算ができていません！
申し訳ない．
*/
D=diag(eigval(A)[,1]);
U=eigvec(A);

RA=U*(D##(1/2))*inv(U);
print RA;

/*確認のため，2乗してみる*/
RA2=RA**(-1);
print RA2;

run;

/*
              A

      0.5       0.3       0.2
      0.1      0.85      0.05
        0         0         1


              RA

0.6989223 0.1858028 0.1152748
0.0619343 0.9156923 0.0223734
        0         0         1


             RA2

      0.5       0.3       0.2
      0.1      0.85      0.05
        0         0         1
 */

必ずしも解が見つかるとは限りませんが、その他にもIMLの非線形最適化サブルーチンを応用して利用することができます。最適化の結果で、目的関数が0に収束していることを確認してください。

*************************************************
rootmatrixnlp.sas
たねや
*************************************************;

proc iml;

A={
 0.5 0.3  0.2
,0.1 0.85 0.05
, 0    0    1
};



n=ncol(a);

start obj(x) global(a,n);
    f=(btran(x,1,n)**2-a)[##];
    return(f);
finish obj;

x0 = btran(a,1,n);
optn = {0 1};
call nlptr(rc,xres,"obj",x0,optn);
RA=btran(xres,1,n);
print RA;

/*
ABSGCONV convergence criterion satisfied.


                                   0.6989223 0.1858028 0.1152748
                                   0.0619343 0.9156923 0.0223734
                                   -4.881E-9 -6.064E-9         1
 */