Rで経済時系列コミュのprecomp ver 0.0 仕様

xですが
| T |
|?T|
です。

虎の巻きってこれのこと？

http://www.ism.ac.jp/editsec/toukei/pdf/50-2-119.pdf

もう、この当たりは読めて実装できるレベルかもね？

＞教科書では、
＞G <- matrix(c(1,0,1,1),2,2);


> G <- matrix(c(1,0,1,1),2,2);
> Q <- matrix(c(1,0,0,1),2,2);
> G%*%Q%*%t(G)
　　　[,1] [,2]
[1,]　　2　　1
[2,]　　1　　1

> G <- matrix(c(1,0,0,1),2,2)
> Q <- matrix(c(1,0,0,1),2,2);
> G%*%Q%*%t(G)
　　　[,1] [,2]
[1,]　　1　　0
[2,]　　0　　1

＞温泉先生
ご指導有難うございました。
頭のたんこぶが、ビックリ仰天しています(笑)

T1AR2MSGモデルに『速度＆加速度』をプラスし、係数を加味しました。
冠の係数は、最近の暴落相場DATAで、経験則で当てはめ計算しています。

(0.7)速度 + (0.01)加速度

> AIC
[1] -18203.07

今後ともご指導よろしくお願いします。

>Gの違いは、システムノイズ

計算としては、どっちでも出来る事ですけど、Q1、Q2の意味としては

G
|1,0|
|0,1|

と書くと、

トレンド成分とトレンドの速度成分のノイズを合算して扱う事にします。となるし

|1,1|
|0,1|

と書くと、

これらは独立したノイズに繋がってるという事です。

>経験則で当てはめ計算しています。

http://www-b2.is.tokushima-u.ac.jp/theses/17/17kameyama_t.pdf

加速度も自動推定させた例です。

#2.2式が間違っていますが。。

>トレンドの変化そのものが、次へと影響を与えるということを重視すれば、

いや、そうだよねぇ。。と思っていたのもあって速度モデルを引っ張り出したと。。

まぁ、加速度もやる気になれば、付けれますが、相場ではどうかと。。

途中から無香料の説明になっているので、へぇ。。って感じで眺めて下さい。

関係ないが、徳島大学の卒論を読んで、

　エレベータ、エスカレータ、リフト等
　乗っている人数にかかわらず、一定の速度で、かつ安全に動作させる

をUKFで制御したら、どうか？ってあって、あの手の機械モノは詳しくないが
重量により、一定速度が、そもそも出てない機械なんだろうか？
また、こんな低速な機械の制御は、ＰＩＤじゃ、どうしていけないか？と
疑問に感じた。

まぁ、オマケ部分で、紙面潰しで書いたモノかもしれないが。。

>無香料の説明になっているので、へぇ。。って感じで眺めて下さい。

　　同窓の方が、『無香料の各種プログラム』を掲載されています。
　　言語はSciLab(Octave類似)だと思います。

　　　Unscented Kalman Filter を用いた連続時間非線形系のパラメータ推定
　　　http://www-b2.is.tokushima-u.ac.jp/theses/17/17sakao_t.pdf

うーん。

ちょい前までは、非線形な時間領域ではニューロが出たり、ニューロじゃ幾らなんでも学習が遅いので、カルマンで学習させるニューロなんかが幅を利かせていたし、今でもまだまだ、この辺りの研究をしてるところは結構あるんですよ。

無理やりっぽい感じはあるけれど、無香料なんかで、大きくこういう実用シーンは変わっていくんでしょうね。

早い話が実際の分布を真似して無理やり計算を合わせに来るわけだから、精度が高いのは寧ろ当たり前なんだろうけど、どこかスマートじゃ無いような気も、ちょっとします。

まぁ、私なんぞは、答えが合えば後はどうでも良いクチなんですが。。

トレンド方程式を作るときに

F
|2,-1|
|1, 0|

とすれば２階差ですが

F
|1,-1|
|1, 0|

とすれば、トレンドの打ち消しになります。
早い話が、トレンドの持つランダムノイズを打ち消すモデルです。

ここで

F
|1,-1,1|
|1, 0,0|
|0,0, A|

Aは1として

H
|1,0,0|
にします。

Aが1なので結局トレンドを取り出すわけですが、先ほどの打消しが効いているので、ループはこのまま計算を回しておいて、x[3]を本当のトレンドとすると、これはトレンドの持つランダム成分の抜けたトレンドが取り出せます。

更に速度項も付け足す事も出来るので、なかなか、おしゃれなトレンド推定になります。

トレンドの日変動（ランダム成分）が無いトレンド抽出法と言う事です。

まぁ、どんなだかは、試してみてください。

コードを上げればいいのですが、今日は現場に忘れてきました。。（ぼそ

> 来年には　precomp ver 0.0 仕様 のバァージョンアップ
> 準備がスマートに整うまで、暫し復習と休憩。。。。

> 来年の話には、少し早いですが・・・・・・・

　　米金融界が米国をつぶす
　　２００８年１１月２８日　　田中　宇
　　　http://tanakanews.com/081128bank.htm

　　　> 来年はドル崩壊と金高騰？
　　　> 金本位制復活論は、もはや嘲りの対象ではない？？？

　　　> 来年は「この世の終わり」に似た事態になりうる。
　　　> 中東での戦争もあり得る。


　　なにか明るい話でもないものですかねぇ〜〜

＞温泉先生
行列表現は苦手なので下記の質問をします。

F
|1,-1, 1, 1|
|1, 0, 0, 0|
|0, 0, A,-1|
|0, 0, 1, 0|

H
|1, 0, 0, 0|

　　Aは2として（ランダム成分が無い）２階差モデル


F
|1,-1, 1, 1|
|1, 0, 0, 0|
|0, 0, A, 1|
|0, 0, 0, 1|

H
|1, 0, 0, 0|

　　Aは1として（ランダム成分が無い）速度モデル


この行列表現でよろしいのでしょうか。

＞76　温泉先生

（ランダム成分が無い）２階差モデル 及び　速度モデル

　早速、試させて頂きました。有難う御座います。
　いづれの試行も、nkf,csv 及び　暴落data　で確認。

　両作図とも、ランダム成分消去の効果は抜群です。
　曲折点を明確にとらえていることには驚嘆しました。

　今後とも、休憩なしでスマートに(笑)　ご指導宜しくお願いします。

ランダム成分消去の効果確認できました。
一部分では先見性も見られます。
有難うございました。

＞行列表現は苦手なので

言われて見れば、これが重要なポイントなのかもね？

ここのカルマンフィルタを使えるかどうかの山場は

１）Rが動かせる
２）データを選択できる
３）作図できる

と言うあたりまでがRにまつわるスキルで

３）思った事を方程式にできる
４）方程式を行列にできる

と言うあたりが、落ちこぼれを作っているポイントなのかも。

どうやってやっているのかと聞かれても、それこそなんとなく
やっている事なので、たいして考えてもいない事なのですが、
中学生に教えるのだとすれば、線形代数的な要素が入っている
のが、難しい所かも。

見てはいるし、自分もやってみたいんだけど、なんだか難しくて
ついていけないし、線形代数と言われて、なんじゃそりゃ？って
人います？

カルマンフィルタを予測の公式だと捉えると、後は公式のマス目
に（行列）自分の思う予測の話を埋めていけば良い事なのだけど、
思う事はあっても、結局マス目が埋まらないなら使えないわけで。

今まで黙っていたけど、あたしもマス目を埋めれるようになりたい
って人がいるなら、このあたりをやりますけど？

＞「じゃあ、関係からマス目を作れる？」

工学としての、線形代数は、結局はこの能力を問われるわけで。

で、どうやったら身に付くか？と言われても自分の場合、それこそ
なんとなく身についたので、まぁ教科書的な例題っぽいのをやるのが
良いのか？と思い進めてきました。

で少しづつ例題から外れるような話も混ぜて行くのが一番現実的な
道筋なのかな？と。

＞それ以前の問題で、マス目を見ても連立方程式をイメージしたり、
＞理解できない人も多いのかも？

そうなんです。

マス目をなんとなく埋めていると言っても、連立条件なんかは
無意識に考慮しているわけで、全く適当に埋めているわけでも
ないわけで。。

機械的に、これはこうだ！と埋めていくのも最初は良い方法だとは
思うのですが、機械的に埋めるにしても、なんの素養も必要ないか？
と言われれば、そうでもないわけで。

けど、これくらいの線形代数は、学校で、どのくらい習う話かと
言えば、せいぜい１時間ぶんくらい教養過程で、ちょろっと習う
程度の話だものね。

∇やdivなんかと一緒に習ったような記憶があるが、あまり定かでも
ない話。

たぶんきちんとした線形代数の知識なんかなくても間に合う話だと
思うんだけどね？

>「じゃあ、関係からマス目を作れる？」となると、立ち往生しています。

私は、もちろん落第塾生です。

まずは、
１）関係から連立方程式を組み立てる。
２）それから、マス目を作る

３）Rが動かせる
４）データを選択できる
５）作図できる

できた図で『関係とマス目』を再確認にし、
さらに『事象と関係式』の認識を深めていくことになります。

　　この繰り返しが、温泉塾の裾野を広げ、かつ、
　　塾生の理解レベルも向上することに繋がれば、結構なことですね。

以前、Rを動かすお勉強と併行し
下記のWebも参考に暇つぶししていました。

我楽多頓陳館 http://www.snap-tck.com/tckinf.html
　　統計学入門
　　　http://www.snap-tck.com/room04/c01/stat/stat.html
　　ベクトルと行列
　　　http://www.snap-tck.com/room04/c01/matrix/matrix.html

ベクトルと行列を書き始めたら、急にWebへのアクセスが減少したそうです。(筆者の嘆き)
やはり、工学系以外の多数の人には馴染みが薄いものなのでしょうね。
そういう私も工学系以外ですが。。。

＞ベクトルと行列を書き始めたら、急にWebへのアクセスが減少

諸工学の礎は、このあたりにあるわけで、どこの学校でも、このあたりを
教えはするんですが、なんと言うか、こういうの使えると、どういいの？
と言う部分は旨く伝えていないわけでね。。

なんというか、大工ごっこしたいなら、トンカチの使い方は習いたいかも
しれないが、大工ごっこに興味ないのに、トンカチの使い方だけ教えられ
ても、やっぱり退屈するというか。。

音楽なんかの習い事は３才くらいから始めるなんて良くあるけれど、これは
右も左もわからんうちになら、トンカチの使い方から仕込む事はできても
大人になってからトンカチの使い方から始めても、だるいだけと言うか。。

けれど一方で、大人の音楽教室なんてのもあって、基礎は無視しても、楽しい
から音楽をやろうと言うコースもあると思うわけで。。

目の前で、1960年生まれのカルマンフィルタが動き出し、どうやら自在に
予測が出来る公式らしいことはわかってきたが、自在に使えるかどうかは、
結局はベクタやマトリクスを工学的な見地から扱えるかどうかに掛かって
いるものらしい。

と言う情況になったわけだけど、だから線形代数をしっかりやれば良いかも？
と思うのは、どうかな？とも思う事なんです。

まぁ、確かにやったほうが良い事とは思いますけど、中学で習うような連立
方程式が解ければ、解く能力そのものは充分なような気もするんです。

工学的な視点とは、連立方程式を解く能力ではなくて、問題を連立方程式に
してしまう能力の事で、解くのは計算機にでも任せれば良い事。

まぁこんな考えでは、３才から英才教育を受けたピアニストには、なれない
のかも知れないけれど、大人の音楽教室に行った程度には、なれるかも
知れないわけで。。

と思っているんですけどね？

左　Trend成分除去Random成分抽出(nkfDATA)画面2x2
中　Trend成分除去Random成分抽出(nkfDATA)期近50拡大作図
右　Random成分除去Trend抽出法(nkfDATA).j 期近50拡大作図

大変興味深いデータが作図されました。

中図に注目です。
他の左、右図が右肩上がり、中図のみ下降傾向

感想
強烈な マイナス下降傾向のRandom成分()が見られます
年末特有の動きとみるか、その後の波乱を予知しているのか？？

F
|2,-1|
|1,0 |

と書けば２階差だが、この2,-1は、どこから出てくるのかと言う話。

１階差とは、たとえば

?Tt = Tt - Tt-1

ここで?Tt=0 つまりランダムウオークとすると

0 = Tt - Tt-1

なので、

Tt = Tt-1

となる。これが１階差なので

F
1

という事。

２階差は、１階差の１階差の事なので

?^2Tt = ?Tt - ?Tt-1

ここで

?Tt = Tt - Tt-1

?Tt-1 = Tt-1 - Tt-2

なので

?^2Tt = Tt -2Tt-1 + Tt-2

となる。

同様に２階差がランダムウオークだとすると

Tt = 2Tt-1 - Tt-2

となり

F
|2,-1|
|1, 0|

を得る。

同様にn階差は、これを繰り返して行く事により求められる。

またトレンド打消しモデルは左辺への移行を行わないだけの話。

練習として、３階差、４階差のFを求めてみましょう。

はい、あってますよ。

階差処理の事を和分I(n)と呼びます。

経済時系列はI(1)程度で定常になる為、あまり大きな次数は必要無いですが
まぁ、中身はこういう処理と言う事で。

たとえばAR(2)なら

F
|1,0|
|0,1|

H
|yn-1,yn-2|

ですが、

F
|1,1,0,0|
|0,1,0,0|
|0,0,1,1|
|0,0,0,1|

H
|yn-1,0,Yn-2,0|

とすれば、速度付きのAR(2)これでも、２階差でも次数は２だけど、これも次数は２。
性能は時と場合による？

そういえば、(a-b)^nの係数ですね。

いつもコードをあげないので、レヴィさんのコードに行列だけ書き込みました。

＃図はデタラメです。

トレンドのAR(2)の係数の両方に速度を持たせたため、かなりイケイケな予測と言う感じ。

トレンドが明らかな時は、かなりヨサゲですが。。
--------------------------------------------------------

nkf<-read.csv("nkf.csv");
nkflog<-log(nkf[,2]);
yobs<-matrix(t(nkflog),1);
N <-length(yobs); time<-1:N; N2<-6;

#データ確認
plot(time,yobs,type="l");

# State Space Representation
#
# x(t+1) = F x(t) + G w(t);# cov{w(t) w(t)} = {Q 0}
# y(t) = H x(t) + v(t);# {v(t),v(t)} = {0 R}

F <- matrix(c(1,0,0,0,0,0, 1,1,0,0,0,0, 0,0,1,0,0,0, 0,0,1,1,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),6,6);
G <- matrix(c(1,0,0,0,0,0, 1,1,0,0,0,0, 0,0,1,0,0,0, 0,0,1,1,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),6,6);
H <- matrix(c(1,0,1,0,1,0),1,6);
Q <- matrix(c(1,0,0,0,0,0, 0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),6,6);
R <- matrix(c(1),1,1);
# ----------------------------------
# Kalman Filter
# ----------------------------------
#
# initial value
x <- matrix(0,N2,1); # x(0|0)
P <- diag(N2); # P(0|0)
ypre <- matrix(numeric(N),1);#予測値
Spre <- matrix(numeric(N),1);#誤差共分散
arpre <- matrix(numeric(N*2),2);#AR(2)
Tpre <- matrix(numeric(N),1);#推定トレンド
MApre<- matrix(numeric(N),1);#MA移動平均
K <- matrix(0,N2,1);
#
x[1,1] = yobs[1,5];
#
for(t in 5:N){
#
# prediction
x <- F %*% x; # x(t+1|t) <- x(t|t)
P <- F %*% P %*% t(F) + G %*% Q %*% t(G); # P(t+1|t) <- P(t|t)
H[1,1] = 1; #推定トレンド
H[1,2] = 0;
H[1,3] = yobs[t-1] - yobs[t-2];
H[1,4] = 0;
H[1,5] = yobs[t-2] - yobs[t-3];
H[1,6] = 0;
ypre[t] <- H %*% x;#予測値
Spre[t] <- sqrt(H %*% P %*% t(H)+R);#誤差共分散
Tpre[t] = x[1,1]; #推定トレンド
MApre[,t] = x[2,1]; #MA移動平均
arpre[,t] = x[5:6,1];#AR(2)係数行列

#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R); # Kalman Gain
x <- x + K * (yobs[,t] - ypre[t]); # x(t|t) <- x(t|t-1)
P <- P - K %*% H %*% P; # P(t|t) <- P(t|t-1)
}
par(ask=TRUE);#作図確認
par(mfrow=c(3,2));#画面を3×2 に分割

#画面3-1-1　原データ(黒)
Min<-min(ypre[,5:N])-0.5;Max<-max(ypre[,5:N])+0.5;
plot(time(5:N),yobs[,5:N],type="l",ylim=c(Min,Max));

#画面3-1-2　原データ(黒)と予測値(赤)
plot(time(5:N),yobs[,5:N],type="l",ylim=c(Min,Max));
par(new=T)
plot(time(5:N),ypre[,5:N],type="l",col="red",xlab="",ylab="",main="",ylim=c(Min,Max));

#画面3-2-1 誤差共分散
plot(time(5:N),Spre[,5:N],type="l",ylim=c(1,2));

#画面3-2-2　AR係数 a1(cyan),a2(magenta)
Min2<-min(arpre[2,]);Max2<-max(arpre[1,]);
plot(time,arpre[1,],type="l",col="cyan",ylab="ar",ylim=c(Min2,Max2));
par(new=T)
plot(time,arpre[2,],type="l",col="magenta",ylab="",ylim=c(Min2,Max2));

#画面3-3-1 MA移動平均
plot(time,MApre,type="l",col="cyan",xlab="",ylab="");

#画面3-3-2 推定トレンド
plot(time,Tpre,type="l",ylim=c(Min,Max));

#以下、期近50を拡大作図
par(ask=TRUE);#作図確認→OKなら図をクリック
par(mfrow=c(1,1));#画面を1×1
N50<-N-49;time50<-N50:N;length(time50)
yobs50<-yobs[N50:N];length(yobs50)
ypre50<-ypre[N50:N];length(ypre50);
plot(time50,yobs50,type="l",ylim=range(ypre50));
par(new=T)
plot(time50,ypre50,type="l",col="red",xlab="",ylab="",ylim=range(ypre50));

外す時は、思いっきり外しているけど、いかにも相場予測って感じ。

こういう風に予測してみたいもんだね。って感じする。

まぁ、最終的にはシグナルだと思えば、結構いいのかも。

>１階差とは、たとえば
>?Tt = Tt - Tt-1
>ここで?Tt=0 つまりランダムウオークとすると
>0 = Tt - Tt-1

これは強制的にランダムウオークを仮定しているので0に
しているわけですが、実際の時系列がランダムウオークに
近くなる場合、Tt - Tt-1は0に近くなると言う事でもある
わけ。

つまり、トレンドそのものを打ち消しあう関係。

この値が、どちらかに振られると言う事は、時系列は
打ち消ししても、更に勾配を持つと考えたほうが妥当
だよと言う事ですが、和分次数が不足気味なのでは？
と言う事。

２階差の打ち消しにすれば、和分の次数は充分なので、
こういう事は起きなくなるんじゃないかな？

どういう風に打ち消しモデルを使うかと言えば例えば

F
|1,-1,1|
|1, 0,0|
|0, 0,1|

これもトレンド方程式なのだけど、お判りですか？

トレンドとは、ランダムノイズの乗った、だいたいの値の
事である。と定義しているわけです。

＞どういう風に打ち消しモデルを使うかと言えば例えば

トレンドを打ち消しても残る成分つまり、トレンドの
速度の存在です。

経済時系列では、絶対に明示的に打ち消しを書かないと
いけない事は、いつもいつも必要になることでも無い
でしょうね。

ほとんどが１回で定常になるし、株で言えば倒産でも
して整理ポストにでも行かない限り、そんなに大きく
打ち消しそこない成分そのものがが値を持つことも
ないだろうし。

さきほど紹介した例は、たとえばほとんど時系列は
I(1)で定常になるのだけれど、時々そうも言って
いられない時も発生してるような時に使うんです。

株で言えば、整理ポストから、ふたたび復活するような
事を繰り返す銘柄でもあれば、有効なのかも知れない
けれど、まぁあまり考え難いかも。

為替で言えば、政府介入なんかをモデリングするので
あれは、有効なのかも？

他に思いつくのは、テロなんかおきた場合の情報伝播速度
とか？

まぁ?Ttがあると考えたほうが妥当だよと言う根拠でも
あるわけなので、速度項が必要か？と言うのを確認する
方法でもあると考えられると言う事です。

2階差だと、これは加速度項の存在ですよね？

檜木さん、どうもです。

パスカルの三角形って言葉が浮かんできました…

トレンドモデルで

F
|1,0, 0|
|0,1,-1|
|0,1, 0|

H
|1,1,0|

とすれば、これも速度付きモデル。
この表現からだと、加速度も付けやすいかもね？