2010年5月15日土曜日

Rで自己相関



前回に続いてRにNINO3のデータを渡し、今回は自己相関を求めてみる。

Rで自己相関を求める関数はacf である。
自己相関を求めるlagは72点(72ヶ月)までとした(lag.max=72)
自己相関の他に
自己共分散 type="covariance"
偏自己相関 type="partial"
も求めることができる。
以下、プログラム

#coding: utf-8
from read_NINO import read_NINO

#for R
import rpy2.robjects as robjects

data=read_NINO()
x=data["NINO3 ANOM"]

year=str(x.year[0])
month=str(x.month[0])

#data to R
nino3=robjects.FloatVector(x)
robjects.globalEnv["nino3"] = nino3

#r script
r=robjects.r
rscripts="""
nino3ts=ts(nino3,start=c(%(year)s,%(month)s),frequency=12)
png('NINO3_R_acf.png')

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

dev.off()
"""
r(rscripts %locals())

図の横の点線は95%の信頼区間をしめす(これより絶対値の小さい相関は95%の信頼区間で無相関ではないとは言えない)。 これを90%などにしたいときはci=0.9などとする。

参考
R の基本パッケージ base, stats 中の時系列オブジェクトの簡易解説
(RjpWiki)

 R言語による時系列分析
(hamadakoichi blog)

R による時系列分析入門
(書籍)

過去のRに関する記事

0 件のコメント: