とあるデータのスペクトル解析
結論から言うと、今回は「やりたいことが実現できなーい。どうすればいいのか分からん!」という話。
なので、読んでもいいことないです。いつも以上に単なる備忘録です。
CSVを読み込んでtsに変換。
jniv<-read.csv("データ.csv") jniv.close<-jniv[2] jniv.close.ts<-ts(jniv.close) jniv.close.ts<-ts(as.numeric(jniv.close.ts))
スペクトル解析
spec.pgram(jniv.close.ts,log="no", xlim=c(0,0.002))
グラフ化せずに生データを見たいときは、「plot=FALSE」とする。
> spec.pgram(jniv10Year.close.ts,log="no",xlim=c(0,0.002), plot=FALSE) $freq [1] 0.0002222222 0.0004444444 0.0006666667 0.0008888889 0.0011111111 [6] 0.0013333333 0.0015555556 0.0017777778 0.0020000000 0.0022222222 [11] 0.0024444444 0.0026666667 0.0028888889 0.0031111111 0.0033333333 [16] 0.0035555556 0.0037777778 0.0040000000 0.0042222222 0.0044444444 ・・・・・・・・・・・・・・・・・ [2236] 0.4968888889 0.4971111111 0.4973333333 0.4975555556 0.4977777778 [2241] 0.4980000000 0.4982222222 0.4984444444 0.4986666667 0.4988888889 [2246] 0.4991111111 0.4993333333 0.4995555556 0.4997777778 0.5000000000 $spec [1] 3.063381e+03 2.232736e+04 3.555620e+04 4.828382e+03 1.111915e+04 [6] 3.021952e+03 9.457354e+03 1.109082e+04 8.995889e+02 2.844646e+02 [11] 9.599256e+02 5.158237e+03 4.985057e+03 2.780245e+03 5.226908e+03 [16] 8.107157e+02 1.879867e+03 1.869383e+03 2.184182e+02 1.788623e+03 [21] 5.462029e+03 1.010075e+03 7.125709e+02 3.086387e+03 2.357322e+03 [26] 8.773852e+02 2.078958e+02 8.100285e+02 9.715714e+01 4.294267e+02 [31] 2.259698e+03 1.571764e+03 1.130904e+03 2.759914e+02 8.386062e+01 ・・・・・・・・・・・・・・・・・・ [2236] 2.092352e+00 4.042008e-01 1.339844e+00 1.866943e+00 7.919000e-01 [2241] 2.731446e+00 9.306263e-01 1.085729e+00 3.989569e-01 1.309799e-01 [2246] 8.265553e-01 2.611334e-01 9.262818e-02 2.484483e+00 3.953657e+00 $coh NULL $phase NULL $kernel NULL $df [1] 1.75974 $bandwidth [1] 6.415003e-05 $n.used [1] 4500 $orig.n [1] 4420 $series [1] "jniv10Year.close.ts" $snames NULL $method [1] "Raw Periodogram" $taper [1] 0.1 $pad [1] 0 $detrend [1] TRUE $demean [1] FALSE attr(,"class") [1] "spec"
◎スペクトルの最大値を求める
max(spec.pgram(jniv10Year.close.ts,log="no",xlim=c(0,0.002), plot=FALSE)$spec)
ここで漸く本題。実は、この最大値に対応した周波数を知りたい。だけど求め方が分からないって話。
周波数は、$freqがnumeric型なので、which関数を使えば次のように(何番目の値かが)求められるかと思っていたが駄目だった。
> vec<-c(1,3,6,3,5) > class(vec) [1] "numeric" > which(vec==1) [1] 1 > which(vec==3) [1] 2 4 > which(vec==6) [1] 3
ところが、which(spec.pgram(jniv10Year.close.ts,log="no",xlim=c(0,0.002), plot=FALSE)$freq==スペクトルの最大値)としても戻り値はinteger(0)。
これは対応する添字が存在しない場合に返却される。
ここでお手上げ。
恐らく、数値計算の誤差の問題だろうから、自分で数値計算して求めないといけないのかなと思っているが、関数無いのかな????