【Rで時系列解析】相乗対数正規過程に従う非ガウスゆらぎ

自然界で観測される時系列には正規分布に従わないものがたくさんあります．正規分布はガウス分布とも呼ばれるので，正規分布に従わない確率過程は「非ガウス過程」と呼ばれます．時系列であれば，「非ガウス時系列」ですが，私はもっとやわらかい響きが好きなので，ここでは「非ガウスゆらぎ」と呼ぶことにします．

　自然界に広く見られる非ガウスゆらぎのモデルとして，私が約20年間追求し続けているものがあります (下記の文献参照)．それが「相乗対数正規過程」(multiplicative log-normal process)です．

この過程は，正規分布に従う確率過程 $\{W[n]\}$ と，それと独立な正規分布に従う過程 $\{Y[n]\}$ を使って，

　 $X[n] = W[n] e^{Y[n]}$

の形で記述できます．

$e^{Y[n]}$ が対数正規分布になっていて，それが掛かっているので相乗対数正規過程と呼んでいます． $e^{Y[n]}$ の部分は，正規分布の標準偏差が確率的に変動することを表していて，結果として $X$ の従う分布は非ガウスになります．そんなの教科書に載ってないし，聞いたことない，という人がほとんどだと思います．この過程の分布は，いわゆる「裾の厚い分布」(fat tailed distribution)になります．

　が，しかし，この過程で記述できる時系列は，自然界にあふれていると断言します．とはいえ，間違っているかもしれないので，私の意見は信じないで，皆さん自身で確かめていってください．

　これから少しずつ，相乗対数正規過程について説明していきます．今日は，相乗対数正規過程の確率密度関数を描くRスクリプトを作ったのでそれを載せておきます．以下では，確率密度関数を描くときに，縦軸が対数軸になっていることに注意してください．この描き方に違和感を持つかもしれませんが，相乗対数正規過程の話をするときは，縦軸は対数目盛にするので慣れてください．この場合，正規分布は2次関数になるので，正規分布と非正規分布 (非ガウス分布)を見分けるのには便利です．

　相乗対数正規過程には，非ガウスパラメタ $\lambda$ があります． $\lambda$ が0のとき，正規分布で， $\lambda$ が増加していくと，中心部が鋭くなり，裾が厚くなっていきます． $\lambda$ を変化したときの，確率密度関数をいくつか載せておきます．以下のRスクリプトを実行すると，相乗対数正規過程 (独立同分布過程になってます)の乱数から推定した分布と理論的分布が重ねて描かれます．図中の灰色破線が，正規分布です．

Rスクリプト

##############################
# non-Gaussian parameter
lambda <- sqrt(0.6)
##############################
# Define integrand
multiLogNorm <- function(sig){
 return(1/(2*pi)/lambda/sig^2*exp(-(log(sig)+lambda^2)^2/(2*lambda^2))*exp(-xi^2/(2*sig^2)))
}
##############################
# Plotting range [-x.max, x.max]
x.max <- 10
###
# Number of data points for numerical test
N <- 10000
##############################
# Numerical integration of the multiplicative log-normal model
x.list <- seq(0,x.max,length.out=500)
n.x <- length(x.list)
pdf <- c()
for(i in 1:n.x){
 xi <- x.list[i]
 pdf[i] <- integrate(multiLogNorm, lower=0,upper=Inf)$value
}
x <- c(-rev(x.list),x.list)
pdf <- c(rev(pdf),pdf)
##############################
# Random numbers obeying the multiplicative log-normal model
y <- rnorm(N,-lambda^2,lambda)
w <- rnorm(N)
z <- w*exp(y)
# PDF estimation
pdf.est <- hist(z,breaks="FD",plot=FALSE)
##############################
# plot
y.min <- min(pdf.est$density[pdf.est$density>0])
y.max <- max(pdf.est$density)
x.max <- 10
plot(pdf.est$mids[pdf.est$density>0],pdf.est$density[pdf.est$density>0],pch=16,col="lightpink",cex.axis=1.2,cex.lab=1.5,
   log="y",xlim=c(-x.max,x.max),ylim=c(y.min/4,y.max*2),yaxt="n",xlab="x",ylab="Probability Density")
curve(dnorm(x),col=8,lwd=3,lty=2,add=TRUE)
par(new=TRUE)
plot(x,pdf,"l",log="y",col=2,lwd=3,xlim=c(-x.max,x.max),ylim=c(y.min/4,y.max*2),xaxt="n",yaxt="n",xlab="",ylab="")
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.01)
a123<-paste(collapse=",",paste(sep="","expression(10^",-5:5,")"))
a123b<-paste(sep="","axis(2,las=1,at=10^(-5:5),label=c(",a123,"),cex.axis=1.2,tck=-0.02)")
eval(parse(text=a123b))