Rでwavelet

前回までのようにRにNINO3のデータを渡し、今回は連続wavelet解析を行ってみる。

Wavelet解析にはwmtsaパッケージを用いる。

以下がプログラムである。
プロットする時、series=TRUEでもとの時系列もプロットする。

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_wavelet.png')

#import library wmtsa
library("wmtsa")
## calculate the CWT of the sunspots series using
## a Mexican hat wavelet (gaussian2)
a<-wavCWT(nino3ts)    
## print the result
print(a)
## plot an image of the modulus of the CWT and the
## time series
plot(a,series=TRUE) 

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

結果

Continuous Wavelet Transform of nino3ts
---------------------------------------
Wavelet           : Mexican Hat (Gaussian, second derivative)
Wavelet variance  : 1
Length of series  : 724
Sampling interval : 0.08333333
Number of scales  : 73
Range of scales   : 0.0833333333333333 to 60.3333333333333

x軸は年。Y軸は0が1年（2^0)、２が4年(2^2)に対応する。

参考
A Practical Guide to Wavelet Analysis
http://paos.colorado.edu/research/wavelets/


書籍
Wavelet Methods for Time Series Analysis (Cambridge Series in Statistical and Probabilistic Mathematics)