統計処理言語Rのすすめ：向日葵工房２：SSブログ

	ブログをはじめるログイン

統計処理言語Rのすすめ　[プログラミング] [編集]

こんにちは。時系列データを統計処理できるアプリに関心がありますがが，Mac備え付けのNumbersはかなり貧弱で使い物になりません（2014/12現在の個人的感想）。そこで，有償のExcel for Macを試用してみようといろいろ調べたのですが，Microsoftに問い合わせても試用版がなぜかうまくインストールできずに断念。そうこうしているうちに，Rという無償の統計処理アプリの存在を知り，2014/10にインストールして使い始めました。結論としては，大変よいものなので皆さんにも大いにおすすめしたいと思います。

　R言語とは何か
ぐぐればいくらでも出てきますが，R言語とは http://ja.wikipedia.org/wiki/R言語にあるとおり，学術分野で発展を遂げてきたオープンソースのフリーソフトウェアのことで，正確にはプログラミング言語というよりはコンソール環境のことです。

R言語は有償のS言語から派生したもののようで，ぐぐった印象ではS言語にお金がかけられない大学の研究室の多くが使っているようです。無償だけあってユーザーの広がりは半端なく，使い方に困ったときなどは公式のヘルプサイトだけでなくユーザーの私的サイトが大いに参考になるため，導入時のハードルは極めて低いです。
　

　R言語のインストール
インストールについてのアドバイスもぐぐればいくらでも出てきますが，ワタクシは http://www.iu.a.u-tokyo.ac.jp/~kadota/R_install_mac.pdf を参考にしました。すなわち，バージョン選択においては最新のものではなくMac OS X 10.6（Snow Leopard）用を使うのが無難だということなので従いました。

　コンソール環境とプログラムコードの関係
ある程度プログラミング経験のある人がR言語を使い始める際にまず知っておくべきことは，R言語のコンソール環境とプログラムコードの関係だと思います。

R言語を起動させると，Windowsのコマンドプロンプトと似たようなコンソール画面が現れます。ワタクシの世代ではMS-DOSやUNIXライクで，むしろ懐かしささえ覚えます。コンソール画面では「>」という表示によってコマンド入力を促していますが，ここで入力すべきコマンドの文法がR言語のプログラムコードの文法そのものです。

こう書くと「ぶつ切りのコマンドはともかく，条件文とかループ文とか関数定義といった複雑なコードには対応していないのではないか」とつい早とちりしてしまいそうですが，他の言語と同様にもちろん対応しています。どういう仕組みかというと，複雑なコードの場合には他の言語と同様に（）や｛｝を用いた入れ子構造が発生しますが，コンソール画面ではすべての入れ子構造が完結しない限り，いくつ改行しても「>」の表示が出ず，コマンドがまだ終わっていないぞというメッセージが発せられるというわけです。

一連のプログラムコードをどうやって一度にコンソール環境に伝えるかというと，source関数を呼び出して引数にプログラムコードを書いたテキストファイルを書くだけです。あるいは，コンソール画面とは別にプログラムコード用のエディタをR上で立ち上げ，コンソール環境に伝えたいコード部分を選択してCommand+Enterしてもよい。とても簡単です。

　プログラム例
ここでは，ドル円の日次四本値時系列データ

Date,Open,High,Low,Close
1998/12/31,114.95,114.96,112.56,113.60
1999/01/01,113.24,114.67,113.20,113.60
1999/01/04,113.30,114.46,111.58,112.02
1999/01/05,112.15,112.42,110.49,111.58
1999/01/06,111.72,113.50,111.56,112.88
1999/01/07,112.88,113.02,109.75,111.31
1999/01/08,………

を読み込んで，適当な時系列処理をした上でグラフ表示するというものです。解説は省略しますが，雰囲気だけ味わってもらえればと思います。

#R，H，L，MaH，MaLのデータフレームから2000-01-01以降のリターン時系列を作る関数
f<-function(m2,trig1,trig2,trig3){
	m1<-ifelse(m2$H>=m2$MaH*trig1,m2$MaH*trig1,m2$R)+ifelse(m2$L<=m2$MaL*trig1,-m2$MaL*trig1,-m2$R)+trig3*(ifelse(m2$H>=m2$MaH*trig2,m2$R-m2$MaH*trig2,0)+ifelse(m2$L<=m2$MaL*trig2,-m2$R+m2$MaL*trig2,0))
	m0<-subset(m1,index(m1)>="2000-01-01")
	return(m0)
}

setwd("/fx_research_r") #参照フォルダを変更する
library(zoo) #時系列ライブラリを使う
m1<-read.zoo("usdjpy.csv",header=TRUE,sep=",") #ヘッダ付き時系列のcsvファイルをデータフレームとして読み込む
m2<-merge(R=log(m1$Close/lag(m1$Close,k=-1)),H=ifelse(m1$High<=lag(m1$Close,k=-1),0,log(m1$High/lag(m1$Close,k=-1))),L=ifelse(m1$Low>=lag(m1$Close,k=-1),0,log(m1$Low/lag(m1$Close,k=-1)))) #前日引値対比の引値（R），高値（H），安値（L）を作成する
library(TTR) #移動平均などが含まれたライブラリを使う
rho<-14
m2<-merge(m2,MaH=lag(EMA(m2$H,rho,wilder=FALSE),k=-1),MaL=lag(EMA(m2$L,rho,wilder=FALSE),k=-1)) #前日のR，H，Lの移動平均を作成する

mm<-f(m2,0.325,0.75,1)

#グラフ作成
plot(as.vector(cumsum(mm)),type="l",xaxt="n")
mmm<-axTicksByTime(mm,"years",format.labels="%Y-%m-%d")
axis(side=1, at=mmm, labels=names(mmm))
abline(v=mmm,h=seq(-1,3,by=0.5),lty=2)

実行結果