4.3.1 時系列のスペクトル密度関数の推定 spectrum()
spectrum()
は時系列のスペクトル密度関数(spectral density)
を推定する．書式：
spectrum(x, ..., method = c("pgram", "ar"))
引数：x
一変量，もしくは多変量時系列method
文字列で，スペクトル密度を推定する手法を指定する．可能な手法は"pgram"
(
既定手法)
と"ar" (AR
モデルの当てはめに基づく)
...
密度推定法，もしくはプロットメソッドplot.spec()
に渡される追加引数 返り値： クラス"spec"
を持つオブジェクトで，少なくとも次の成分を持つリスト：freq
スペクトル密度が推定された周波数のベクトル(
可能性としてフーリエ周波数を 近似)
．単位は，単位時間当たり(
観測値間隔当たりではなく)
のサイクルの逆数spec freq
に対応する周波数位置に於けるスペクトル密度のベクトル(
一変量時系列の場合
)
，もしくは行列(
多変量時系列の場合)
coh
一変量時系列ではNULL
．多変量時系列の場合，coh
のi+(j-1)*(j-2)/2
列はx
のi, j
列(i<j)
間のコヒーレンシ(coherency
，干渉性)
の2
乗値を含むphase
一変量時系列ではNULL
．多変量時系列の場合，異なる時系列間のクロススペクトルのフェイズ
(cross-spectrum phase)
を含む．形式はcoh
と同様series
時系列の名前snames
多変量時系列の場合，成分時系列の名前method
スペクトル密度を計算した手法spectrum()
は二つのメソッドspec.gram(), spec.ar()
に対するラッパ関数 *3 で ある．スペクトル密度関数はS-PLUS
との互換性のために，スケール1/frequency(x)
を掛けてあり，関数の範囲は(-frequency(x)/2, +frequency(x)/2]
となる．より普 通には，スケールを2π
とし範囲を( − 0.5, 0.5]
としたり，1
として( − π, π]
とする．もし可能なら
plot.spec()
関数を用いて信頼区間がプロットされる．これは非対称であり，センターマークの幅はバンド幅と等値である．
注意：クラス
"spec"
のオブジェクトの既定のプロットは非常に複雑であり，エラーバー と既定のタイトル，副タイトルそして軸ラベルを持つ．これらの既定値は，適当な作図パ ラメータを補うことで，変更することができる．関連：
spec.ar(), spec.pgram(), plot.spec().
# 黄体形成ホルモン量データlhを使用(以下の図を参照)
> par(mfrow=c(2,2)) # 画面を2x2分割
> spectrum(lh); spectrum(lh, spans=3); spectrum(lh, spans=c(3,3))
> spectrum(lh, spans=c(3,5))
# 英国の呼吸器疾患死亡者数データUKLungDeathsを使用(以下の図を参照)
> spectrum(ldeaths, spans=c(3,3)); spectrum(ldeaths, spans=c(3,5))
> spectrum(ldeaths, spans=c(5,7))
> spectrum(ldeaths, spans=c(5,7), log="dB", ci=0.8)
# ARモデルフィットによるスペクトル密度推定例(以下の図を参照)
# 多変量時系列のプロットに付いてはspec.pgram()のヘルプを参照
> par(mfrow=c(1,2)) # 画面を1x2分割
> spectrum(lh, method="ar"); spectrum(ldeaths, method="ar")
*3wrapper function．(より複雑な)関数(群)を簡便に使用するための関数．
0.1 0.2 0.3 0.4 0.5
0.010.050.200.50
frequency
spectrum
Series: lh Raw Periodogram
bandwidth = 0.00601
0.1 0.2 0.3 0.4 0.5
0.050.100.200.501.00
frequency
spectrum
Series: lh Smoothed Periodogram
bandwidth = 0.0159
0.1 0.2 0.3 0.4 0.5
0.20.40.6
frequency
spectrum
Series: lh Smoothed Periodogram
bandwidth = 0.0217
0.1 0.2 0.3 0.4 0.5
0.20.40.6
frequency
spectrum
Series: lh Smoothed Periodogram
bandwidth = 0.0301
0 1 2 3 4 5 6
2e+031e+045e+042e+05
frequency
spectrum
Series: ldeaths Smoothed Periodogram
bandwidth = 0.173
0 1 2 3 4 5 6
2e+031e+045e+042e+05
frequency
spectrum
Series: ldeaths Smoothed Periodogram
bandwidth = 0.241
0 1 2 3 4 5 6
2e+035e+032e+045e+04
frequency
spectrum
Series: ldeaths Smoothed Periodogram
bandwidth = 0.363
0 1 2 3 4 5 6
35404550
frequency
spectrum (dB)
Series: ldeaths Smoothed Periodogram
bandwidth = 0.363, 80% C.I. is (−1.84, 2.41)dB
0.0 0.1 0.2 0.3 0.4 0.5
0.10.20.51.0
frequency
spectrum
Series: lh AR (3) spectrum
0.0 0.1 0.2 0.3 0.4 0.5
2e+035e+031e+042e+045e+041e+052e+055e+051e+062e+06
frequency
spectrum
Series: ldeaths AR (10) spectrum
スペクトル密度関数
(上左)黄体形成ホルモン量データ．原ス ペクトルと種々の平滑化
(上右) 英国肺疾患死亡者数データ．原ス ペクトルと種々の平滑化
(下左) AR モデルの当てはめによるス ペクトル推定．左は黄体形成ホルモン量 データ，右は英国肺疾患死亡者数データ
4.3.2 AR モデルの当てはめによりスペクトル密度を推定する， spec.ar()
spec.ar()
はx
にAR
モデルを当てはめ(
もしくは既に存在する当てはめ結果)
によ り，スペクトル密度を推定する．返り値はクラス
"spec"
のオブジェクトで，もしplot = TRUE
ならコンソールには出 力されない．書式：
spec.ar(x, n.freq, order = NULL, plot = TRUE, na.action = na.fail, method = "yule-walker", ...)
引数：
x
一変量(
多変量版はまだ未移植)
時系列．もしくはar()
による当てはめ結果n.freq
プロットされる点の数order
当てはめられるAR
モデルの次数．省略されると次数はAIC
で選ばれるplot
ペリオドグラムをプロットするか？na.action
欠損値に対する処理関数method ar()
当てはめに対する手法... plot.spec()
に引き渡される作図パラメータ返り値： 返り値はクラス
"spec"
のオブジェクト．もしplot=TRUE
なら値は不可視注意：
Thomson
等はAR
スペクトルは誤解を与えやすいと強く警告している．多変量ケースはまだ移植されていない．
4.3.3 ペリオドグラムの計算 spec.pgram()
spec.pgram()
は高速フーリエ変換によりペリオドグラムを計算する．オプションとして結果を修正
Daniel
平滑法(
終端の重みを半分にする移動平均)
による系列で平滑化する．もし
plot = TRUE
なら結果はコンソールに表示されない．書式：
spec.pgram(x, spans=NULL, kernel, taper=0.1, pad=0, fast=TRUE, demean=FALSE, detrend=TRUE, plot=TRUE,
na.action=na.fail, ...)
引数：
x
一変量もしくは多変量時系列spans
奇数からなるベクトルで，ペリオドグラムの平滑化に使われる修正Daniell
平滑法の幅を与える
kernel
もしくはクラス"tskernel"
の平滑法taper
テイパリングされるデータの割合．半cosine bell
テーパが，データの先頭と末尾のこの比率の部分に適用される
pad
埋められるデータの割合．データの長さを増すために，データの末尾にこの割合 分0
が加えられるfast
論理値．TRUE
なら，0
で埋めることにより，系列の長さを約数が極めて多い数 にするdemean
論理値．TRUE
なら，時系列から平均を引くdetrend
論理値．TRUE
なら，時系列から線形傾向を引く．同時に平均も引き去るplot
ペリオドグラムをプロットするか？na.action
欠損値処理関数... plot.spec()
に引き渡される作図パラメータ返り値： クラス
"spec" (spectrum()
を見よ)
のリストで，次の追加成分を持つ：kernel kernel
引き数．もしくはspans
から構成された核関数df
スペクトル密度関数推定値の分布を近似するカイ2
乗分布の自由度bandwidth Bloomfield
で定義されているような核関数平滑法の同値バンド幅taper taper
引き数の値pad pad
引き数の値detrend detrend
引き数の値demean demean
引き数の値生のペリオドグラムはスペクトル密度の一致推定量ではないが，隣接した値は漸近的に
独立である．したがって，スペクトル密度が滑らかという仮定の下で，生のペリオドグラ ムを平滑化することによりスペクトル密度の一致推定量を得ることができる．高速フーリ エ変換を助けるために，時系列はその長さが高度な合成数
(
約数が多い)
になるように0
で埋められる．これはpad
ではなくfast
で制御される．時系列の平均が引かれている ので，原点でのペリオドグラムは理論的に0
になる(
テーパリングにより影響を受けるか も知れない)
．これは平滑化の過程で隣接値の線形補間で置き換えられ，その周波数での 値は返されない．関連：
spectrum(), spec.taper(), plot.spec(), fft().
4.3.4 時系列の両端を削る spec.taper()
cosine-bell
関数*4 を，時系列x[,i]
の最初と最後のp[i]/2
分の観測値に適用する．返り値は新しい時系列である．
書式：
spec.taper(x, p=0.1)
引数：x
一変量もしくは多変量時系列p
削られる両端部分の割合．スカラー(
全ての時系列に適用される割合)
，もしくは 時系列の長さ(
各時系列に適用される割合)
のベクトル4.3.5 累積ペリオドグラムをプロットする cpgram()
cpgram()
は累積ペリオドグラムをプロットする．書式：
cpgram(ts, taper=0.1, main=paste("Series: ",deparse(substitute(ts))), ci.col="blue")
引数：
ts
一変量時系列taper
ペリオドグラム作製時にテーパリングする割合main
主タイトルci.col
信頼区間幅用の色返り値： 返り値は無く，副作用として正方領域に累積ペリオドグラムが作図される
# 黄体形成ホルモン量データlh使用(以下の図を参照)
> par(pty = "s", mfrow = c(2,2)) # 画面を2x2に分割．プロット点のタイプ指定
> cpgram(lh) # 原データのcpgram
> lh.ar <- ar(lh, order.max = 9) # AR(3)モデルを当てはめ
> cpgram(lh.arresid, main = "AR(3) fit to lh") # その残差系列の cpgram
# 英国の肺疾患死亡者数データldeathsを使用(以下の図を参照)
> cpgram(ldeaths) # 両性死亡者数数のcpgram
*4taperとは端を削って尖らすことを意味する．cosine bellとは関数(1−cos(x))/2を指す．この0か ら1まで値が変化する関数を時系列の両端部分に掛けることにより，時系列の値を両端で徐々に0に近 づける操作がcosine taperである．
0.0 0.1 0.2 0.3 0.4 0.5
0.00.20.40.60.81.0
frequency Series: lh
0.0 0.1 0.2 0.3 0.4 0.5
0.00.20.40.60.81.0
frequency AR(3) fit to lh
0 1 2 3 4 5 6
0.00.20.40.60.81.0
frequency Series: ldeaths
黄体形成ホルモン量データの累積 ペリオドグラム
原データ(上左)
AR(3) モデルを当てはめた残差
(上右)
英国の肺疾患死亡者数(両性)デー タの累積ペリオドグラ(下左)