第４回R勉強会＠東京 で話してきた -「R言語による時系列分析」

「R言語による時系列分析」を、第４回R勉強会＠東京 (2010/04/24) で話してきました。
第4回R勉強会＠東京（Tokyo.R#04） : ATND

双方向形式の進行で質疑応答も含め80分で行いました。 私の双方向形式の講義進行で発言・質問して下さった方々、懇親会で話せた方々、どうもありがとうございました。資料、ソースコード、実行結果を以下に記載します。

講師資料

[Ｒ勉強会][データマイニング] R言語による時系列分析

View more presentations from Koichi Hamada.

AGENDA

• 自己紹介
• 時系列分析とは
• データ操作
• モデル
  • 確率密度
  • 期待値と分散
  • ARモデル
  • ARMA/ARIMAモデル
  • ARFIMAモデル
  • ARCHモデル
• 最後に

資料中リンク

今回の資料アップではPowerpointの資料をキレイにSlideShareにアップする方法を行っています。(今まではPPTをそのままアップしていたのですが、形が崩れてしまっていました。)

ただこの方法では、資料内のリンクが外れてしまうようなので、以下に、資料内の各リンクを記載します。

内容	リンク
データマイニング+WEB勉強会＠東京 Google Group	http://groups.google.com/group/webmining-tokyo
Twitter (#TokyoWebmining アナウンスも)	http://twitter.com/hamadakoichi
自己紹介	http://iddy.jp/profile/hamadakoichi
毎週３時間ダンスコーチ	http://www.youtube.com/hamadakoichi
博士論文	http://hosi.phys.s.u-tokyo.ac.jp/~koichi/PhD-thesis.pdf
PRML Revenge #1	http://atnd.org/events/4115
LBIビジネス講演会	http://www.lbi.gr.jp/modules/eguide/event.php?eid=31
Tokyo.R #7	http://wiki.livedoor.jp/syou6162/d/Tsukuba.R%237

R言語ソースコードと実行結果(抜粋)

資料中に記載されている「R言語の実行可能なソースコード」と「実行結果」の抜粋を以下に記載します。

時系列データの形式と属性

ソース

class(UKgas) #型の表示
start(UKgas) #開始時間
start(UKgas) #終了時間
frequency(UKgas) #時間ごとの観測回数
window(UKgas,c(1975,2),c(1979,3)) #時系列の切り出し

実行結果

> class(UKgas) #型の表示
[1] "ts"
> start(UKgas) #開始時間
[1] 1960    1
> start(UKgas) #終了時間
[1] 1960    1
> frequency(UKgas) #時間ごとの観測回数
[1] 4
> window(UKgas,c(1975,2),c(1979,3)) #時系列の切り出し
      Qtr1  Qtr2  Qtr3  Qtr4
1975       321.8 177.7 409.8
1976 593.9 329.8 176.1 483.5
1977 584.3 395.4 187.3 485.1
1978 669.2 421.0 216.1 509.1
1979 827.7 467.5 209.7    

時系列表示

ソース

ts.plot(UKgas) #時系列表示

実行結果

ラグ処理

ソース

tmp <- ts(1:120,start=c(1995,6),frequency=12)#時系列オブジェクト作成
class(tmp)#型の表示
tmp #データ表示

実行結果

> tmp <- ts(1:120,start=c(1995,6),frequency=12)#時系列オブジェクト作成
> class(tmp)#型の表示
[1] "ts"
> tmp #データ表示
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1995                       1   2   3   4   5   6   7
1996   8   9  10  11  12  13  14  15  16  17  18  19
1997  20  21  22  23  24  25  26  27  28  29  30  31
1998  32  33  34  35  36  37  38  39  40  41  42  43
1999  44  45  46  47  48  49  50  51  52  53  54  55
2000  56  57  58  59  60  61  62  63  64  65  66  67
2001  68  69  70  71  72  73  74  75  76  77  78  79
2002  80  81  82  83  84  85  86  87  88  89  90  91
2003  92  93  94  95  96  97  98  99 100 101 102 103
2004 104 105 106 107 108 109 110 111 112 113 114 115
2005 116 117 118 119 120    

差分

ソース

ts.plot(UKgas) #時系列表示
ts.plot(diff(UKgas))#差分の時系列表示

実行結果

時系列データの統計量

ソース

op<- par(mfrow=c(1,3)) #描画設定：1行3列画面表示
acf(UKgas) # UKgasの自己相関プロット
acf(UKgas,type="covariance") #UKgasの自己共分散プロット
acf(UKgas,type="partial")    #UKgasの偏相関プロット
par(op)#作業前の描画設定に戻す

実行結果

スペクトル分析

ソース

op <- par(mfrow=c(1,4))#描画設定：1行3列画面表示
ts.plot(UKgas) #UKgasの時系列表示
spec.pgram(UKgas) #スペクトル解析：Fast Fourier Transformation(FFT)でpreiodogram算出
spec.pgram(UKgas,spans=c(3,3)) #Daniell平滑化(終端の重み半分での移動平均）を用いたスペクトル解析
spectrum(UKgas,method="ar") #自己回帰(AutoRegression)スペクトル解析
par(op)#作業前の描画設定に戻す

実行結果

ARモデル

ソース

result <- ar(UKgas) #自己回帰モデルを求める
result #結果の表示
summary(result) #結果概要の表示

result$order #次数
result$ar #係数
result$aic #情報量基準 AIC
result$resid #残差

実行結果

> result <- ar(UKgas) #自己回帰モデルを求める
> result #結果の表示

Call:
ar(x = UKgas)

Coefficients:
      1        2        3        4        5        6  
 0.5284  -0.2951   0.5844   0.3489   0.0436  -0.2709  

Order selected 6  sigma^2 estimated as  8496 
> summary(result) #結果概要の表示
             Length Class  Mode     
order          1    -none- numeric  
ar             6    -none- numeric  
var.pred       1    -none- numeric  
x.mean         1    -none- numeric  
aic           21    -none- numeric  
n.used         1    -none- numeric  
order.max      1    -none- numeric  
partialacf    20    -none- numeric  
resid        108    ts     numeric  
method         1    -none- character
series         1    -none- character
frequency      1    -none- numeric  
call           2    -none- call     
asy.var.coef  36    -none- numeric  
> 
> result$order #次数
[1] 6
> result$ar #係数
[1]  0.5283543 -0.2950795  0.5843581  0.3489334  0.0436247 -0.2709347
> result$aic #情報量基準 AIC
         0          1          2          3          4          5          6 
210.901838 172.327892 169.210792  23.118177   5.487060   6.233848   0.000000 
         7          8          9         10         11         12         13 
  1.870267   3.864514   4.145583   6.041680   8.012734   9.437802  11.338009 
        14         15         16         17         18         19         20 
 13.054909  15.024158  16.923270  18.821075  20.767077  22.746861  23.546878 
> result$resid #残差
            Qtr1        Qtr2        Qtr3        Qtr4
1960          NA          NA          NA          NA
1961          NA          NA -16.5105970 -15.5797229
1962   1.3345706  -2.3562659 -15.1521582 -12.8411247
1963   4.4391020 -16.2988690 -11.9961810 -25.1842864
1964 -11.0887011 -10.1459008 -14.8095946 -10.9605092
1965  -1.0789915  -5.8988827 -10.5310923 -11.3919090
1966   2.8957223 -12.6996593 -12.2360652 -15.8516772
1967  -0.4086809  -1.9832162  -9.1272299 -14.7601502
1968  12.8559373  -2.7453608 -14.9442404 -20.4671974
1969  13.9106165  -0.6184421 -13.5154616 -11.0730546
1970  -7.6264080  -4.6173176  52.3086300 -57.8427761
1971  74.6081756 -95.4255172   8.1170932  57.6766859
1972  20.2962762  -1.3814550 -27.3815187  72.2425142
1973  -1.9531400  11.8478073 -52.9774380  43.4686310
1974  45.0283670  31.7420844 -32.7149745  43.9502016
1975  13.7374828  32.5237181 -56.7049192  32.2153362
1976  80.7403115 -11.8486568 -25.7619222  56.6565968
1977  -9.3898598  75.9812196 -67.1689121  53.8269680
1978  39.3368973  48.1479112 -37.1483146  37.2788211
1979 152.0353431 -11.0467834 -23.5350079  -7.3341375
1980  69.4714686 -73.6926404  39.7052873  94.6030401
1981  35.3963450   4.6089229 -49.6412402  71.7849861
1982  74.4718045  10.6182195 -39.3433397   4.4690531
1983  43.1555798  81.5995997 -41.0599132  44.2127564
1984  74.8983047 -26.7772388  -5.1032682  36.0848241
1985 156.1266286  -2.6794890  38.6536935   5.1290691
1986 150.1372066   9.1595882  59.0100486 -62.7249513

予測

ソース

arresult <- ar(UKgas) #自己回帰モデルを求める
result <- predict(arresult,n.ahead=20) #予測する
result$pred #予測値 時系列オブジェクト
result$se #標準誤差 時系列オブジェクト
ts.plot(UKgas,result$pred) #表示：元データ、予測値

実行結果

> arresult <- ar(UKgas) #自己回帰モデルを求める
> result <- predict(arresult,n.ahead=20) #予測する
> result$pred #予測値 時系列オブジェクト
          Qtr1      Qtr2      Qtr3      Qtr4
1987 1053.9970  600.7282  316.9538  748.8121
1988  981.4623  546.8309  308.6325  708.0975
1989  912.3671  504.7132  297.3734  674.2334
1990  849.5541  468.2595  289.8132  643.3941
1991  793.8445  437.2220  284.6219  616.0413
> result$se #標準誤差 時系列オブジェクト
          Qtr1      Qtr2      Qtr3      Qtr4
1987  92.17482 104.24961 104.25994 111.21667
1988 137.88835 147.29676 147.31935 151.18307
1989 170.10928 175.66459 175.66783 178.40438
1990 191.54920 195.23883 195.24297 197.23555
1991 206.93975 209.38176 209.40298 210.89961

参考文献

Rによる時系列分析入門

Rによるデータサイエンス データ解析の基礎から最新手法まで

関連リンクまとめ

第４回R勉強会＠東京の関連リンク

• 第4回R勉強会＠東京（Tokyo.R#04） : ATND
• 第4回R勉強会＠東京（Tokyo.R#04）を開催しました - yokkunsの日記
• http://d.hatena.ne.jp/nokuno/20100424/1272121385

• R勉強会： 第３回データマイニング+WEB 勉強会＠東京 (Tokyo.Webmining#3) を開催しました−「R言語による クラスター分析 -活用編-」 - hamadakoichi blog

• R勉強会： 第２回データマイニング+WEB 勉強会＠東京 (Tokyo.Webmining#2) を開催しました - 「はじめてでもわかる R言語によるクラスター分析」- - hamadakoichi blog
• R言語プログラミング： クラスター分析 - 階層的クラスタリング - hamadakoichi blog
• 第２回 データマイニング+WEB 勉強会＠東京 (Tokyo.Webmining#2) を開催します - hamadakoichi blog
• R勉強会 講師をしてきた： 「はじめてでもわかる 統計解析・データマイニング R言語入門」 - hamadakoichi blog
• 第３回R勉強会＠東京 で話してきた -「プロセス・リソース・グラフと数理統計解析」 - hamadakoichi blog

• [データマイニング+WEB勉強会][R勉強会] はじめてでもわかる 統計解析・データマイニング R言語入門
• [データマイニング+WEB勉強会][R勉強会] はじめてでもわかる R言語によるクラスター分析 - 似ているものをグループ化する-
• [R勉強会][データマイニング] プロセス・リソース・グラフと数理統計解析