大気海洋統計データ解析

-1-

第２章

第２章

時間領域における単変量解析 ..................................................................................................... 1

2.1.

変動の持続性................................................................................................................................. 1

2.2.

時系列モデリング ......................................................................................................................... 2

2.2.1.

赤色ノイズ............................................................................................................................. 2

2.2.2.

赤色ノイズの生成方法 ......................................................................................................... 4

2.2.3.

自己回帰過程(autoregressive process) .................................................................................. 5

2.2.4.

自己回帰モデルの推定 ......................................................................................................... 6

2.2.5.

移動平均(MA)過程................................................................................................................ 8

2.2.6.

自己回帰移動平均(ARMA)過程 .......................................................................................... 8

2.2.7.

自己回帰和分移動平均(ARIMA)過程 ................................................................................. 8

2.2.8.

コレログラム ......................................................................................................................... 82.3.

統計推定・検定の自由度 ............................................................................................................. 8

2.3.1.
2.4.

相関に関する無相関時間の導出 ....................................................................................... 11

フィルタリング........................................................................................................................... 13

参考文献....................................................................................................................................................... 14

第２章

時間領域における単変量解析

2.1.

変動の持続性

普通の統計の教科書では，N 個の独立なデータがあってなんとかは自由度 N-1 のかんとか分布に従う，と何気なく書いてある．
「独立な」は書いてない場合も多い．しかし，７月が猛暑であれ
ば多分８月も暑いと予想できるように，気温など大気・海洋の統計データ解析で扱うデータは，
時間・空間的に独立ではなくてある程度の持続性(persistency)を持つ．その持続時間よりも短いサンプリング間隔であれば，データを独立とみなすことができず，近くのデータ同士は遠くのデー
タ同士よりも強く相関する．この性質を，連続相関(serial correlation)があると言う．連続相関という用語はただし，具体的な相関係数を意味するのではなく，連続相関があるということの実際の計算はラグ相関などで行なわれている．
持続性をもたらすメカニズムの一つは海洋の熱的な慣性である．すなわち，2m の厚さの熱容量
が大気全層の熱容量に匹敵する海洋の温度が変化すると，その記憶は長く保たれる．海洋表面水温の経年変動が数ヶ月に渡って持続することは，Namias ら (1988)が研究し，この論文は 1970 年
代の気候レジームシフトについての先駆け研究としても有名である．しかも，熱的慣性はしかし数ヶ月後に失われるわけではない．大気からの熱的強制の偏差によって，冬季の海洋混合層内に温度偏差が作られると，浅い夏季の混合層の下で温度偏差が保存され，また次の年の冬季に表面水温に影響を与える．このメカニズムを，SST 偏差の再発現 (re-emergence)と呼ばれ，Alexanderら(1999)などによって研究されている．このような熱的慣性が存在すれば，それが大気にフィードバックすることで，大気もまたある程度の過去の状況に引きずられると考えられる．また，海洋の熱的な慣性によらなくても，大気自身のブロッキングも過去の記憶を未来に残すメカニズムとなっている．

大気海洋統計データ解析

3

-2-

第２章

Winter Arctic Oscillation Index

2
1
0
-1
-2
1950

1960

1970

1980

1990

2000

Sapporo Air-Temp
-1
-3

o

C

-2
-4
-5
-6
1950

1960

1970

1980

1990

2000

図 2.1 冬季の北極振動指数（上）と，札幌の気温（下）．点線が各年のデータで，実線が５年移動平均．各年が独立とは言い難い．

2.2.

時系列モデリング2.2.1.

赤色ノイズ

持続性を持つ現象について，何事かをモンテカルロ法によって調べようとすれば，Matlab のrandn で表現される乱数は，適切な時系列データの代理時系列とはなり得ない．この場合に最も広く使われかつ最も簡単な時系列の生成モデルは，次式で示される赤色ノイズモデル(red noisemodel)である．

y (t ) = ry (t − 1) + σ W ε (t )(2.1)

ここで時系列 y(t)を表現する．ここで，t は離散化された時間ステップ，r はラグ１自己相関係数，

ε は分散１のガウスノイズ．σ W はガウスノイズの標準偏差である．後に示すように，赤色ノイズ
は長周期成分のスペクトルが強いために，光のスペクトルとの類推から赤色ノイズと呼ばれる．
またガウスノイズは，白色ノイズ(white noise)とも呼ばれる．なお上の式で，r=0 とすれば，白色ノイズに他ならないので，赤色ノイズは白色ノイズを含んでいる．また，容易に見て取れるように，r が大きいほど持続性も大きい．
赤色ノイズの実現例を図 2-1.に示す．一見周期的に見える長期変動が赤色ノイズには存在し得
るが，後に示すように赤色ノイズは特定の周期性を持たない（スペクトルピークが無い）ので，視察では周期性があるように見えるとしても偶々そう見えるだけである．

大気海洋統計データ解析
r=0.20 (0.21), std=1.02 (0.91)

3

2

4

1

2

0

0

-1

-2

-2

-4

-3

0

50

100

r=0.80 (0.78) std=1.67 (1.53)

6

Amp

Amp

第２章

-3-

-6

0

Time

50

100
Time

図 2-2. ラグ１自己相関係数 0.2（左）と 0.8（右）に対する，赤色ノイズモデルの時系列．測定されたラグ１自己相関係数はパネル上部のカッコ内に示している．やはりパネル上部に理論的な標
準偏差と観測された標準偏差とを，観測された標準偏差にカッコを付けて示している．
赤色ノイズモデルを用いて，代理時系列を生成するには，通常(2.1)を用いる．この際，二つのパラメータ，r と σ W を知らなくてはならない．赤色ノイズの標準偏差を， σ R と表すことにすれば，

σ R2 ≡ y (t ) 2 = r 2 y (t − 1) 2 + 2rσ W y (t − 1)ε (t ) + σ W2 ε (t ) 2= r 2σ R2 + σ W2であるから，

σ W = (1 − r 2 )σ R

⇔ σ R = σ W / (1 − r 2 )(2.2)となる．ここで<>は期待値（アンサンブル平均）．観測された時系列を赤色ノイズと仮定すれば，観測された分散が σ R であるので，上式から σ W が得られる．
2

多くの場合，ラグ１相関係数は，時系列から通常の相関係数の計算と同じように求められる．
しかしこの方法で観測時系列からラグ１相関係数を求めると，求められたラグ１相関係数は若干小さく見積もられる，すなわちラグ１相関係数の推定が偏っていることが実は知られている．またラグ１相関係数の過小評価は(2.2)から白色ノイズ分散の過小評価につながり，これらのパラメータを使って行なう統計的なテストは有意となりやすくなる．この偏りはただしそれほど大きくなく，100 個のデータから求める場合では 0.02 程度である．ラグ１相関係数の不偏推定方法は，Allen and Smith (1996)に示されているので，万一気難しい査読者に攻撃されたりして興味が湧いた
読者は参照されたい．
なお，赤色ノイズモデルは，一次のマルコフ過程(first order Markov process)とも呼ばれる．マルコフ過程の定義は，未来の確率状態が過去によらずに決まる確立過程である．一次のマルコフ過大気海洋統計データ解析

-4-

第２章

程は多分，１時間ステップ前の状態には依存するが，それよりも前の状態には依存しないということであろう．

2.2.2.

赤色ノイズの生成方法

赤色ノイズを生成するには，生成するべき時系列の長さを tlng, 時系列の個数を nsmpl, ラグ１
自己相関係数を r1 とすると，例えば次のように生成することができる．

rtsrs=randn(tlng,nsmpl);% tlng×nsmpl の白色ノイズを生成

rtsrs(1,:)=rtsrs(1,:)/sqrt(1-r1^2);%最初の時間の赤色ノイズを規格化

for t=2:tlng
rtsrs(t,:)=rtsrs(t-1,:)*r1+rtsrs(t,:);end
各時系列の最初のデータは，白色ノイズを赤色ノイズと同じ標準偏差になるように規格化して与えている．
演習問題 2-1 r1=0.5, nsmpl=10000, tlng=10 で赤色ノイズを生成し，分散と標準偏差のヒストグラム
を描け．また，分散の平均が理論値に一致し，標準偏差の平均は一致しないことを示せ．m 番目nの時系列の分散を (1/ n)∑ y (t ) で求める．
2

i =1

回答

m

大気海洋統計データ解析

第２章

-5-

variance method 1, r=0.50

1400
1200
1000
800
600
400
200
0

0

1

2

3

4

5

6

2

2.5

3

std method 1, r=0.50

1400
1200
1000
800
600
400
200
0

0

0.5

1

1.5

図 2-3 時系列の分散と標準偏差の分布図．各パネルの垂直線の黒は理論値，赤は時系列パラメー
タの平均値である．分散については両者がよく一致しているために，黒線が赤線に隠されている．

2.2.3.

自己回帰過程(autoregressive process)赤色ノイズでは，ラグ１の情報の伝達のみを考慮したが，同様の考えで多くのラグを含めることができる．このようなより一般的なモデルを，自己回帰モデル(autoregressive model)と呼ぶ．
p 次の自己回帰過程は

y (t ) + φ (1) y (t − 1) + .... + φ ( p) y (t − p) + σ W ε (t ) = 0(2.3)

で表現される．φ (i ) はモデルの係数である．時系列 y がこの過程に従うことを， yt∼AR(q)と書く．
赤色ノイズモデルは１次の自己回帰モデルである．係数 φ は AR1 ではラグ１相関係数に一致したが，２次以上の AR モデルでは一致しない．なお，赤色ノイズモデルとは，各項の係数は逆の符号としている．これは Matlab の慣用にあわせており，両方の流儀がある．
自己回帰モデルがどのような振る舞いをするかは，平均的な振る舞いについてしか議論できな
いので，まずアンサンブル平均， y (t ) = y (t ) ，を取って次式を得る．

y (t ) + φ (1) y (t − 1) + .... + φ ( p ) y (t − p ) = 0
ここで，大気海洋統計データ解析

-6-

第２章

y (t ) = λ t − p y (t − p )(2.4)

を仮定すると（つまり y (t ) = λ y (t − 1) ），λ p + φ (1)λ p −1... + φ ( p) = 0という特性方程式と呼ばれる式を得る．明らかに， λ > 1 では y が発散する解を持つし，λ = 1 でもノイズが加わる効果によって平均が徐々にゼロからずれて分散が増加する．このような，望ましくないな特徴をもたないためには，特性方程式の全ての根が λ < 1 である必要がある．この場
合，遠く離れたデータの影響は減衰するので，十分に初期値から離れれば初期値の影響は無視できる．この時 AR 過程は定常であると言う．AR1 では， φ (1) < 1 が定常の条件である．なお一般に，確率過程のアンサンブル平均で得られる平均・分散・自己相関関数とが，時間的に変化しな
いことを，確率過程が定常1であるという．

2.2.4.

自己回帰モデルの推定

自己回帰モデルの推定は，統計的2予測の実行，この章で説明する自由度の推定，次章で説明する Maximum Entropy Method (MEM) Spectrum を推定するなど，様々な局面で用いられる．この節では，自己回帰モデルの推定方法を説明しよう．
自己回帰モデルの係数 φ と，ラグ τ の自己共分散関数（ ρ (τ ) ≡ y (t ) y (t + τ ) ）と関係は，次の関係がある．

ρ (τ ) ≡ y (t ) y (t − τ )= − {φ (1) y (t − 1) + ... + φ ( p ) y (t − p ) + σ W ε (t )} y (t − τ )(2.5)= − φ (1) y (t − 1) y (t − τ ) + ... + φ ( p) y (t − p) y (t − τ ) + σ W ε (t ) y (t − τ )ここで，ラグがゼロでなければ，ρ (τ ) + φ (1) y (t − 1) y (t − τ ) + ... + φ ( p ) y (t − p ) y (t − τ ) + σ W ε (t ) y (t − τ ) = 0
ρ (τ ) + φ (1) ρ (τ − 1) + ... + φ ( p ) ρ (t − p ) = 0となる．上の式を分散 σ = ρ (0) で割って，自己共分散を自己相関 r に変えれば，
2r (τ ) + φ (1) r (τ − 1) + φ (2) r (τ − 2) + ... + φ ( p ) r (t − p ) = 0(2.6)

が得られる．この関係から AR1 については， r (1) = −φ (1) となるので，すでに述べたように AR1のモデル係数はラグ１自己相関係数に他ならないことが確認できる．ただしラグ 0 については，(2.5)から落とした白色ノイズの項も含めることが必要であり，この場合

1

厳密には，この定常性は，弱定常性あるいは共分散定常性であって，強定常とは区別される．
強定常性は， y (t ), t = 1,..., n の分布が， y (t + τ ) の分布と一致することを言う．正規分布に従う確率過程は，弱定常であれば強定常でもある．この資料では定常とは常に弱定常の意味で使われる．
2
予測には統計的な方法と，力学的なモデル計算による方法とがある．

大気海洋統計データ解析

-7-

第２章

ρ (0) + φ (1) y (t − 1) y (t − τ ) + ... + φ ( p) y (t − p ) y (t − τ ) + σ W ε (t ) y (t ) = 0∴1 + φ (1)r (1) + φ (2)r (2) + ... + φ ( p )r ( p ) + σ W2 / σ 2 = 0(2.7)

である．自己相関は，対称関数（正負のラグで同じ値）であることを利用して，左辺のラグが１からｐまでの(2.6)および(2.7)を行列式であらわせば，r (2) r (1) 1r (1) r (1)
r (2)

#
 #
 r ( p − 1) r ( p − 2)r ( p) "
1 − σ W2 / σ 2  φ (1)  " r ( p − 1)  
r (1) 

φ (2) =  r (2) " r ( p − 2)   #  
%
# 
#
 

φ ( p)  


"
1 
 r ( p) 

(2.8)

これはちょっと書き換えて，行列が Toeplitz(テープリッツ行列：成分が対角線に沿って一定である行列)になるようにできる．

r (1)
r (2)
 1
 r (1)1r (1) r (2) r (1)
1

#
#
 #
 r ( p) r ( p − 1) r ( p − 2)"r ( p )   1   −σ W2 / σ 2 " r ( p − 1)   φ (1)  
0






=
" r ( p − 2) φ (2)
0


 
%
#  #  
#






"
1  φ ( p)  
0


(2.9)

この関係式が有名なユール・ウォーカー方程式(Yule-Walker equation)である．Matlab では lpc,
pyulear で利用できる．
ユール・ウォーカー方程式で自己回帰モデルを得るには，次数 p を決めてやらなくてはならな
い．次数を大きくすればするほど，モデルの当てはまりがよくなるのは当然なので，それだけを考えると，真の p は小さいのにもかかわらず大きな p を選んでしまい勝ちである．この落とし穴を回避し，より妥当に p を決めるためのよく使われる方法としては，赤池の情報量基準(Akaike’s
Information Criterion, AIC)が最小になる次数 p を用いる方法である．AIC は自己回帰モデルに限らず広く使われる概念であって，その詳しい説明は柴田(1996)の７章を参照されたい．自己回帰モ
デルについての AIC は，

AIC = N + N ln 2πσ W2 + 2( p + 1)

(2.10)

であり，ここで N は時系列データの個数である．上式のうち自己回帰モデルに依存する部分のみ取り出して，

AICの変化部分 = N ln σ W2 + ( p + 1)

(2.11)として AIC を用いることもよく行なわれる．AIC を用いて AR モデルを推定する手続きは次の通りである．
1) あらかじめ定めた次数の範囲で，各次数についてユール・ウォーカー方程式から AR モデルの

大気海洋統計データ解析

第２章

-8-

係数を求める．
2) 各次数の AR モデルの AIC を求める．
3) 最小となる AIC を与える次数の AR モデルを採用する．
なお，ユール・ウォーカー方程式の最大のラグ相関係数 φ ( p ) も未知数であるとして，同時に決定するバーグ(Burg)のアルゴリズムは，ユール・ウォーカーのアルゴリズムより若干優れていると考えられている．バーグのアルゴリズムは，Matlab では pburg で利用できる．
後に述べるように，自己回帰モデルにもとづいてラグ相関係数を推定する応用もある．その場
合，ラグが p 以下であればラグ相関係数は観測された相関係数であり，ラグが p よりも大きいのであれば(2.6)を用いてラグ相関係数が求められる．

2.2.5.

移動平均(MA)過程

q 次の移動平均過程 (moving-average process, MA process)はy (t ) + ε (t ) + θ (1)ε (t − 1) + θ (2)ε (t − 2) + .... + θ (q )ε (t − q),t = 1,..., nで表現され，時系列 y がこの過程に従うことを，yt∼MA(q)と書く．
無限次 MA 過程は適当な条件のもとでは AR(1)で表現される．任意の有限次の MA 過程は常に定常である．

2.2.6.

自己回帰移動平均(ARMA)過程

p,q 次 の 自 己 回 帰 移 動 平 均 過 程 (autoregressive-moving average process, ARMA process) を ，yt∼ARMA(p,q)と書く(ハーベイ,1985,２章(1.25))y (t ) + φ (1) y (t − 1) + .... + φ ( p) y (t − p)+ε (t ) + θ (1)ε (t − 1) + .... + θ (q)ε (t − q) = 0,

t = 1,..., T

p,q が共にゼロで無い場合に混合過程と呼ぶ． ARMA 過程の自己回帰部分が，AR 過程の定常性を満たす場合には，ARMA 過程もまた定常となる．

2.2.7.

自己回帰和分移動平均(ARIMA)過程

自己回帰和分移動平均(autoregressive-integrated-moving moving average)過程，略して ARIMA 過程と呼ぶ．ARIMA 過程は，y(t)-y(t-1)が ARMA 過程となる過程である．ARIMA 過程は，局所的なトレンドを表現することができ，そのために一定の非定常性を表現することができる．

2.2.8.

コレログラム

非負のラグに対する標本自己相関関数のプロットをコレログラム(correlogram)と呼ぶ．なお，
自己相関関数とは，ラグの関数として複数の自己相関係数を扱う場合に用いる呼び名である．コ
レログラムはしばしば大気海洋の研究でも使われ，その場合主としてコレログラムに振動現象が見られるかどうかが議論されることが多い．

2.3.

統計推定・検定の自由度

連続相関がある場合に，どのようにして推定・検定に用いる自由度を推定すればよいのかには，

大気海洋統計データ解析

-9-

第２章

自己相関関数を用いて推定する標準的な手法があるので，それを紹介しよう．基本的な考え方は，実効的に独立な標本間の時間(time between effectively independent samples, Trenberth (1984)およびMetz (1991)の呼び名)，もしくは有効無相関時間(effective decorrelation time)と呼ばれる Te で，デー
タのサンプル数を割って，有効自由度(effective number of degrees of freedom)を得る．平均・分散・相関および回帰に関する有効無相関時間を下にまとめて示そう．
平均について

N
N
 τ 
 τ 
Te = ∑ 1 −  rxx (τ ) = 1 + 2∑ 1 −  rxx (τ )NN
τ =− N τ =1 

(2.12)

分散についてNN τ  2 τ  2
Te = ∑ 1 −  rxx (τ ) = 1 + 2∑ 1 −  rxx (τ )NN
τ =− N τ =1 

(2.13)N τ Te = ∑ 1 −   rxx (τ ) ryy (τ ) + rxy (τ )ryx (τ ) N
τ =− N 
相関・回帰についてN τ = 1 + rxy (0)ryx (0) + 2∑  1 −   rxx (τ )ryy (τ ) + rxy (τ )ryx (τ ) N
τ =1 

(2.14)

ここで r は相関関数，N は標本の数（時系列の長さ）， τは ラグである．相関関数は，共分散関数
を全分散（つまり分母の和が 1…N である）で規格化して

1 N −τ
1 N −τ
+
τ
x
(
t
)
x
(
t
)
∑
∑ y(t ) y(t + τ )N − τ t =1
N − τ t =1
rx x (τ ) ≡, ryy (τ ) ≡
,
1 N
1 N
2
2
∑ x(t )∑ y(t )
N t =1
N t =1rxy (τ ) ≡1 N −τ
∑ x(t ) y(t + τ )N − τ t =1
1 N
1 N
x(t ) 2 ∑ y (t ) 2∑
N t =1
N t =1

定義される．注意するべきことは，無相関時間は（従って有効自由度は）現象に内在する量ではなく，データの長さにも，またどの統計解析を行うかにも依存することである．なお，上の式を使う場合の参考文献を上げておこう．どの式についても，Lieth (1982)で上げられているようだが，残念ながら元論文を入手できていない．平均については統計分野の Bayley and Hammersley (1946)
や大気海洋では，Lieth (1982)や Trenberth (1984)が妥当なとこである．相関は Bartlett が 1930 年代に定式化したらしくその内容は Bartlett(1955)にも記されている．相関について，大気海洋に導入したのは Davis (1976)である．分散は誰が出したのかよく分からない．Lieth は 1982 年の論文の前にもいくつか関連する話題について論文やレポートを書いているようだが(Trenberth 1984 のreference による)，どれも入手しづらい文献で彼の実際の貢献がどのようなものか判断できかねている．
なお，実際の計算に当たっては，大きなラグについては上式の分子の推定が使えるデータの数が少なく，良い推定が得られないという問題がある．そこで，Katz (1982)や Trenberth (1984)では，
AR モデルを推定して，その AR モデルに基づいて自己相関係数を計算する方法を勧めている．こ大気海洋統計データ解析

第２章

-10-

の方法は，ラグが大きなラグ相関の値が大きく，無相関時間の推定に影響が強い時に使うとよいだろう．もし，ラグが大きなラグ相関が小さいのであれば，AR モデルの推定を経由する必要はそれほどない．また Trenberth(1984)では，ラグが大きいと上式の分子の推定が使えるデータの数がN −τ

少なくて悪くなるために，あえて偏りのある推定値， (1/ N )∑ x(t ) x(t + τ ) を使うことも提案し
t =1

ている．
結局，有意性の検定を行うには，検定の種類に応じて，上のいずれかの無相関時間を求め，それから得られる有効自由度で各々の検定の自由度を置き換える．また無相関時間の推定には，標本相関係数をそのまま使うと危ない場合には AR モデルを経由して，相関係数を推定する．AR モ
デルの推定は，ユール・ウォーカー方程式で行い，AIC で次数を決める．というように，データが独立である統計の教科書では簡単な有意性検定も，独立ではない時系列解析では一気に大変になる．

2.3.1.

平均に関する無相関時間の導出

真の平均を x ，標本平均を { x} とする．標本の２乗から真の平均の２乗を引くと，v ≡ ({ x} − x )2= { x} − 2 { x} x + x
2

2= { x} − x22(2.15)

である．標本平均を定義に従って書き直して，

1 N
v =  ∑ x(t )  N t =1


2

− x

2

(2.16)であるので，ここで x を偏差と真の平均の和で表して， x = x '+ µ , µ = x ，として

1 N
v =  ∑ x '(t ) + µ  N t =1


2

− x

2

=

1 N N
x '(t ) + µ )( x '(s ) + µ ) − µ 22 ∑∑ (
N s =1 t =1

=

1 N N
∑∑ ( x '(t ) x '(s ) + ( x '(s ) + x '(t )) µ + µ 2 ) − µ 2
N 2 s =1 t =1

{

1 N N
∑∑ x '(t ) x '(s ) + x '(s ) + x '(t ) x + µ 2
N 2 s =1 t =1
1 N N
= 2 ∑∑ x '(t ) x '(s )N s =1 t =1=}− µ(2.17)
2である．ここで， τ = j − i が同じ値を持つものについての和という形にまとめれば，v=

1 N
σ2 N  τ 
N
(
)
−
τ
ρ
τ
=
(
)
∑
∑ 1 −  rxx (τ )xx
N 2 τ =− NN τ =− N  N (2.18)

大気海洋統計データ解析

第２章

-11-

互いに無相関な正規分布する標準偏差が σの 分布を， N e 点で平均すれば分散は

σ 2 / N e = σ 2Te / N(2.19)

である．ただし Te は上で述べた有効無相関時間で Te = N / N e で定義される．(2.18)と(2.19)を比較
すると，
N
 τ 
Te = ∑ 1 −  rxx (τ )Nτ =− N (2.20)となって，これは，平均に関する有効無相関時間と自己相関関数との関係式(2.12)に他ならない．
なお，(2.17)の偏差の定義を標本平均からのずれとすると以下のようになって，簡単な形は得ら
れない．

1 N

v =  ∑ x '(t ) + x  N t =1


2.3.2.

2

− x

2

=

1 N N
2
x '(t ) + x )( x '(s ) + x ) − x
2 ∑∑ (
N s =1 t =1=

1
2
x '(t ) x '(s ) + ( x '(s ) + x '(t ) ) x + x 2 ) − x
2 ∑∑ (
N s =1 t =1

=

1 N N
x '(t ) x '(s ) + x '(s ) + x '(t ) x + x 2
2 ∑∑N s =1 t =1N

(2.21)

N

{

}− x

2

相関に関する無相関時間の導出(2.14)式を導出してみよう．<>を期待値もしくは真の平均とし，{}を標本平均(sample mean) とする．すなわち，ある観測量 q(i), i=1,…,N について {q} = (1/ N )

∑

N
n =1

q (i ) である．二つの時系列

x, y について，真の共分散を < xy >= V ，sample variance を {xy} = V + v としよう．問題なのは，
真の共分散と標本分散(sample variance) が一致せず，標本分散が通常大きいことである．
定義から，v = {xy}− < xy > であるから，v = {xy} − 2{xy} < xy > + < xy > となる．ここで，2

2

2

v2 の真の平均は，

v 2 = {xy}2 − 2 {xy} xy + xy

2

= {xy}2 − xy2で与えられる．最後の等号では， {q} = q を用いた．さらに，この式を書きなおして，

大気海洋統計データ解析

第２章

-12-

1 N
 1 N2v 2 =  ∑ x(t ) y (t )   ∑ x( s ) y ( s )  − xy
 N t =1
  N s =1

=

1 N N
2
x(t ) y(t ) x( s ) y( s ) − xy
2 ∑∑N t =1 s =1

=

1 N N
2
x(t ) y(t ) x( s ) y ( s ) − xy
2 ∑∑N t =1 s =1となる．ここで第一項の<>内が結合正規分布するとすれば，結合正規分布する４つの確率変数の積の期待値は，次のように，二つの変数間の共分散の全ての組み合わせの和で表される(証明は例
えば，川嶋・酒井，p.33)．すなわち， xm , xn , x p , xq を４つの結合正規分布する確率変数であるとすれば，xm xn x p x q ≈ xm xnx p x q + xm x p

xn x q + xm x q

xn x pであることを用いて，

1 N N
2
 x(t ) x( s ) y(t ) y( s ) + x(t ) y( s) y(t ) x( s ) + x(t ) x(t ) y ( s ) y ( s )  − xy
2 ∑∑ 
N t =1 s =1
1 N N
= 2 ∑∑  x(t ) x( s ) y(t ) y( s ) + x(t ) y( s) y(t ) x( s ) N t =1 s =1となる．ここで， τ = j − i が同じ値を持つものについての和という形にまとめれば，
v2 =

v2 =

1 N  τ 
∑ 1 −   ρ xx (τ ) ρ yy (τ ) + ρ xy (τ ) ρ yx (τ ) N τ =− N  N  (2.22)となる．ただし， ρは 次式で定義される自己共分散関数および相互共分散関数である．

ρ xx (τ ) = x(t ) x(t + τ ) , ρ yy (τ ) = y (t ) y (t + τ ) ,ρ xy (τ ) = x(t ) y (t + τ ) , ρ yx (τ ) = y (t ) x(t + τ )である．さらに，意味を掴みやすい表現を得るために，x, y 各々の分散を Qx = x

2, Qy = y 2 と書いて，rxx (τ ) = ρ xx (τ ) / Qx , ryy (τ ) = ρ yy (τ ) / Q y , rxy (τ ) = ρ xy (τ ) / Qx Q yと定義される相関係数を用いれば，(2.22)は，

v2 =


QxQy  N  τ  ∑ 1 −   rxx (τ )ryy (τ ) + rxy (τ )ryx (τ )  N τ =− N  N 
(2.23)と表わされる．
一方，互いに無相関かつ自己相関を持たない正規分布する確率過程が与える，標本時系列

大気海洋統計データ解析

-13-

第２章

x(t ) ', y (t ) ', t = 1,..., N の共分散の期待値は，
Qx Qy
1 N
x(t ) y (t ) =∑N t =1Nである．そこで，上の共分散の増加分， v

2，に含まれる独立なデータの数を N e (有効自由度，

the effective number of degrees of freedom (effective sample size), Metz (1991))と定義すれば，

v2 = Qx Q y
Qx Q y  N 
Qx Q yτ = ∑  1 −   rxx (τ ) ryy (τ ) + rxy (τ ) ryx (τ )   =N τ =− N NNeN ∆t / Teと書くことができる．ここで，Te は有効無相関時間である．明らかに，Nτ Te = ∆t ∑  1 −   rxx (τ ) ryy (τ ) + rxy (τ ) ryx (τ ) N
τ =− N 

であり，これは(2.14)に他ならない．積分時間スケールから自由度は DOF = N ∆t / Te と求めることができる．上式の[]内第二項は通常小さいので，無視されることが多い．
なお，得られているデータ長 N の範囲で，相関係数が十分に速くゼロになるならば，(2.14)の近似表現として，∞Te = ∆t ∑  rxx (τ ) ryy (τ ) + rxy (τ ) ryx (τ ) 
τ =−∞も用いられる．ここで， τの和を取る範囲 ( −∞, ∞ ) は実際的には，相関係数がゼロとみなせる範囲とすれば良い．この場合，積分時間スケールは，データを生成する物理過程によって定まる定数であると考えることができる．

2.4.

フィルタリング

時系列解析の様々な局面で，フィルタリングが用いられる．大きく分けると，フーリエ変換・逆変換を使う方法と，時間領域で重み付き移動平均を適用する方法（非再帰的フィルター，Finite
Impulse Reponse (FIR) フィルターとも呼ばれる），時間領域で重み付き移動平均を再帰的に使う方
法（再帰的フィルター，Infinite Impulse Response (IIR)フィルターとも呼ばれる）の３つである．
フーリエ変換・逆変換は通常のデータ解析であれば，まず考えるべき第一の方法であり，一般にFIR・IIR フィルターよりも良い特性をもつ（詳しくは Press ら(1993)を参照のこと）．時間領域の（非再帰的な）重み付き移動平均の代表例は，重みを等しくしたｎ点移動平均である．この移動
平均はフィルターの性能はあまりよくないが，その単純さと広く使われているので大気海洋コミュニティーで理解されやすいというなじみ深さの点から，これで済む程度であればこれを使うのは良い方法である．バターワース・フィルターなどの再帰的フィルターは，時間領域の非再帰的なフィルターに比べると，一般に高い性能を実現しやすい．一方で，位相がずれるなどの問題もある．また，フィルターの性能はフーリエ変換・逆変換に比べて高いわけではない．リアルタイ
ムにデータを処理する必要があれば FIR か IIR を使わねばならず，その場合には通常 IIR が用いられる．

大気海洋統計データ解析

-14-

第２章

Matlab ではフーリエ変換が fft，逆変換が ifft である．時間領域のフィルターは filter で，特に位相のずれをなくしたい時は，時間領域に順方向と逆方向の２回フィルターを通す filtfilt も用意されている．様々なフィルターの設計もサポートされており，例えば butter でバターワースフィ
ルターが得られる．

参考文献
Allen, Myles R., and Leonard A. Smith, 1996: Monte Carlo SSA: Detecting irregular oscillations in the
Presence of Colored Noise. J. Climate, 9 (12), 3373–3404.
Alexander , M. A., C. Deser, and M. S. Timlin, 1999: The Reemergence of SST anomalies in the North
Pacific Ocean. J. Climate, 12 (8), 2419-2433.
Bartlett, M. S., 1955: An introduction to stochastic processes with special reference to methods andapplications. Cambridge university press, pp. 388. (third edition が 1978)
Bayley G. V. and J. M. Hammersley, 1946: The 'effective' number of independent observations in
autocorrelated time series. J. Roy. Statist. Soc. Suppl., 8(2), 184-197.
Davis, R. E., 1976: Predictability of sea surface temperature and sea level pressure anomalies over the
North Pacific Ocean. J. Phys. Oceanogr. 6 (3), 249-266.
Hanan, E. J. 1970: Multiple Time Seies. Wiley.
Katz, R. W., 1982: Statistical evaluation of climate experiments with genal circulation models: A
parametric time series modeling approach. J. Atmos. Sci., 39, 446-1455. 1982.
川嶋

弘尚・酒井

英昭：現代スペクトル解析．pp. 372, 森北出版株式会社，1989.

Lieth, C. E.,1982: Statistical methods for the verification of long and short range forecasts. ECMWF
seminar on Problems and prospects in long and medium range weather forecasting. 313-334.
[available from ECMWF]
Metz, W., 1991 : Optimal relationship of large-scale flow patterns and the barotropic feedback due tohigh-frequency eddies, J. Atmos. Sci., 48, 1141-1159.
Namias, J., X.,

Yuan, and D. R. Cayan, 1988: Persistence of North Pacific Sea Surface Temperature andAtmospheric Flow Patterns. J. Climate, 1, 682-703.
Press, W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, 1993: Numerical Recipes in C [日本語版],
技術評論社, pp. 685.
Trenberth K. E., 1984: Some effects of finite sample size and persistence on meteorological statistics. Part
I: autocorrelations. Mon. Wea. Rev., 112, 2359-2368.