【Rで時系列解析】長時間相関時系列の生成 - ケィオスの時系列解析メモランダム

時系列の長時間相関解析とか，フラクタル解析とか，そう呼ばれる解析を試してみるためには，時系列のサンプルが必要です．そのために，Rを使って，時系列を生成する方法を紹介します．

ここでは，longmemoパッケージを使います．パッケージをインストールしていない場合は

install.packages("longmemo")

を実行してください．

データを保存するフォルダは，事前に各自作成して下さい．
例えば，

D:\Document\時系列解析

の場合は，バックスラッシュ("\"，日本語環境では￥マーク)を，"/"に書き直す必要がありますので，下記のRスクリプト3行目にある

DIR <- "...... ここに各自，データを保存したいフォルダへのパスを指定 ......"

の部分を

DIR <- "D:/Document/時系列解析"

に変更してください．

以上の点（パッケージのインストールとDIRの指定）に注意して，このスクリプトを実行してみてください．

############################
# ファイルの保存場所の指定
DIR <- "...ここにフォルダを指定..."
############################
# パッケージのインストールが必要な場合は以下を実行
# install.packages("longmemo")
####################################
# パッケージのロード
require(longmemo)
####################################
# データを保存するファイル名
FN <- "sample.csv"
####################################
# パラメタの設定
# ハースト指数(0 < H < 1)
H <- 0.8
# 時系列の長さ(長くすると，計算時間が長なり，ファイルサイズがでかくなる)
n <- 10000
####################################
# 非整数ガウスノイズの生成
#####
## 疑似乱数列を固定した場合はseedに整数値を設定すること　(設定しない場合は勝手にランダムっぽくなります)
# seed <- 1
# set.seed(seed)
#####
# ここで，時系列を生成しています．
x.fGn <- simFGN0(n,H)
####################################
# 時系列の可視化
par(mfrow=c(3,2)) # 2x2のアレンジでプロット
par(cex.lab=1.4,cex.axis=1.2) # ラベルなどのサイズ指定
par(pty="m")# プロットが長方形
# 時系列のプロット
plot(0:(n-1),x.fGn,"l",col=4,xlab="i",ylab="x[i]",main=paste("Sample time series (H=",H,")",sep=""),las=1,xaxs="i")
var.x <- var(x.fGn)
####################################
# パワースペクトルの推定
###
# スペクトル推定のために時系列を分割する部分区間の数
# (精度は高いのは１０以上，最低でも４推奨)
n.psd <- 10
####
n.sub <- floor(n/n.psd)
for(i in 1:n.psd){
 if(i == 1){
  f <- spectrum(x.fGn[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$freq
  psd <- data.frame(spectrum(x.fGn[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec) 
 }else{
  psd <- cbind(psd,spectrum(x.fGn[((i-1)*n.sub+1):(i*n.sub)],plot=FALSE)$spec)
 }
}
# fごとにPSDの平均を計算
psd.mean <- rowMeans(psd)
# psdの片側面積を分散と等しく設定
df <- f[2]-f[1]
psd.a <- sum(psd.mean)*df
psd.est <- psd.mean/psd.a*var.white
# log-logプロットに直線をフィット
fit <- lm(log10(psd.est)~log10(f))
# グラフを描画
plot(log10(f),log10(psd.est),xlab="log10 f",ylab="log10 S(f)","l",las=1,col=4,main=paste("Estimated PSD (slope=",format(fit$coefficients[[2]],digits=2),")",sep=""))
abline(fit,lwd=2,lty=2,col=2)
####################################
# 散布図
par(pty="s") # プロットが正方形
plot(x.fGn[1:(n-1)],x.fGn[2:n],pch=16,col=4,las=1,xlab="x[i]",ylab="x[i+1]",main="Scatter plot (lag = 1)")
plot(x.fGn[1:(n-2)],x.fGn[3:n],pch=16,col=4,las=1,xlab="x[i]",ylab="x[i+2]",main="Scatter plot (lag = 2)")
plot(x.fGn[1:(n-4)],x.fGn[5:n],pch=16,col=4,las=1,xlab="x[i]",ylab="x[i+4]",main="Scatter plot (lag = 4)")
plot(x.fGn[1:(n-8)],x.fGn[9:n],pch=16,col=4,las=1,xlab="x[i]",ylab="x[i+8]",main="Scatter plot (lag = 8)")
####################################
# 生成データの保存 
setwd(DIR)
write.table(x.fGn,FN, sep = ",", append=FALSE, quote=FALSE, col.names=FALSE,row.names=FALSE)

詳しい使い方は，Youtubeで解説します（公開予定）．