頭のリフレッシュ　〜　時系列分析３

今回はその時扱ったデータの差分について考えてみる。
Rでプロットするとこんな感じ。

なんだか周期的に、-50くらいへ下振れするようなので、周期性を調べてみる。
次の加法モデルを仮定して、この3つへ分解するということ。

観測値＝トレンド＋周期変動＋残差

Rのstatsパッケージの関数stlを使用する。

> basis.ts<-as.ts(basis)
> stl(basis.ts, s.window="per")
 stl(basis.ts, s.window = "per") でエラー:  一変量系列だけが許されます

あれ、出来ない。なんでや。「エラー: 一変量系列だけが許されます」って。
どうやら、numeric型の時系列データにしなければならないようだ。

> basis.ts<-ts(as.numeric(basis.ts))
> stl(basis.ts, s.window="per")
 stl(basis.ts, s.window = "per") でエラー: 
   系列は周期的でないか、2 未満の周期しか持ちません

うーん、今度は「系列は周期的でないか、2 未満の周期しか持ちません」だと？！
「周期を調べるのがお前（stl）の仕事じゃないのか？！意味わからん！」と思いつつ、このblogで「季節変動成分の目安となる周期インデックスを割り振ってやる」ことの必要性を知る。

Rで季節変動のある時系列データを扱う際に、何よりも最も重要なのはts(){stats}関数で生データをnumeric型で渡し、さらにfrequency引数を用いて季節変動成分の目安となる周期インデックスを割り振ってやることです。

まだあまり理解していないが、とりあえず指定してみる。

> basis.ts<-ts(as.numeric(basis.ts), frequency=150)
> plot(stl(basis.ts, s.window="per"))

とりあえず、エラーは出なくなった。

frequency=50としてみたらこんな感じ。さて、これからどうしようか。

> basis.ts<-ts(as.numeric(basis.ts), frequency=50)
> plot(stl(basis.ts, s.window="per"))
> summary(stl(basis.ts, s.window="per"))
 Call:
 stl(x = basis.ts, s.window = "per")

 Time.series components:
    seasonal              trend              remainder         
 Min.   :-15.616582   Min.   :-26.084161   Min.   :-198.06401  
 1st Qu.: -4.766603   1st Qu.: -8.817512   1st Qu.: -14.10818  
 Median :  0.754002   Median : -3.488457   Median :   2.43078  
 Mean   :  0.033760   Mean   : -3.957966   Mean   :  -0.07349  
 3rd Qu.:  4.448671   3rd Qu.:  1.416383   3rd Qu.:  18.45296  
 Max.   : 11.198481   Max.   : 11.624928   Max.   : 184.69391  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
      9.215       10.234    32.561        30.690
   %  30.0         33.3     106.1         100.0 

 Weights: all == 1

 Other components: List of 5
 $ win  : Named num [1:3] 12261 77 51
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 1227 8 6
 $ inner: int 2
 $ outer: int 0
>