国立1防災科学技術センター研究報告第42号

1989年2月

519．281550．34

自己回帰モデルの係数の推定方法の比較
堀貞一喜＊

国立防災科学技術センター

ComparisonofAR

Coeffici6nts11lstimated
TimeSeries

from

Data

By

SadakiHORI

M｛㎝αlRθ∫ω舳α肋γ∫・γ1）づ3ατ〃P舳〃㎝，〃α〃Abs虹act
Th…t…g・…i・・（AR）m・d・1i・・h・…t・・i・・dby…t・f…ffi・i・・t・whi・h1i・…1ycombinetosuccessive

estimate

coefficients

present

signa1s

thetovalueARcontained

inancoefficients．The

observed

minimizeaクκ〃cκ㎝〃γoγwhich

and

value

the

predicted

timeseries，There

aretwotraditiona1method（method1）deteminesfmm

is

the

ca1cu1ated

past

asa

differenceobservations．In

themethods

theAR

between

the

other

method（m・th・d2）・whi・hi・・・…t1yp・・p…di・th・∫棚m・th・df・・…pPli・・ti・・t・・p・・t・・1
analyses，thedata

soAR

coefficients

that

they

are

estimateof

a1inear伽θsigna1．On

there1ationshipbiased1inearextract

detemined

autoregressiveto

re1ationshipbetweenother

hand，theAR

among

the

methods

becomes

c1ear

with

tme

va1ue

Gaussian

white

noise

by

method1are

tamination
ARof

coefficientsnoise

power．In

with

atime

seriesthree
theAR
makes

amountsnoto

alongercertainbetteraconsistent

negative

obtained

byof

two

methodsathepast〃〆砂data

time

to

the

estimatequite

estimate

effectleve1，thetheestimateincrease

betweentrue

seriesoftheofmethod2represent

presentκo｛∫2＿difference

of

the

the

mmber

of

the

the

data

＊第2研究部地震前兆解析研究室

一1一

ARof

of

two

estimatesareincrease

for

con−

of

the

the

data

is

examined

for

a

data

pointsusedinobtainedbycoefficients

thecomparesinusoidsthecoefficients．The

anlength．But

and

a

the

how

coefficientsARdata1ength

the

to

data1ength

the

two

the

sensitive

value．Fina11y，wecontaining

ofun−

byobtained

beseries

an

the

where

worse

time

andinvestigate

good

an

observed

numerica1experiments

sinusoids．When

correlation

afrom

coefficients，First，we

found

with

method2from
using

the

coefficientsratio．Then，then㎝decaying

estimate
theof

S／N

c㎝taining

method1becomesthrough

contrast，method2gives

re1atively1ower

calculation

nize

obtainedsignalsan〃伽κ㎝κ5伽α1．Method1gives

successive〃oづ∫θ＿伽2signals．The

comparedtheminimizewecandeviation

thewhose

recog−oftheresolvabi1ity
frequenciesare国立防災科学技術センター研究報告
very

c1ose

toeachmethod1．A11the
as

other

three

for

the

estimationpressionofthebyalso

findofAR

model

assumesthenoise．As

that

a

resultofmethod2isdemonstrate

the

how

method2is

coefficients．The

autoregressive

method1implicit1ytaminated

andexperiments

for

past

resu肚，the第42号

difference

a例ぬツtime

data

present

resu1tsmuch

thesな〃α〔s，inbetterthanthatsuperior

tomethod1

fromseries．The

predicting

1989年2月

the

least

present

thisconfusionsquaressignalinof

ex−condition

are

notin

con−

expression，independentofth・榊〃・・1・θbut・otind・pe・d・・t・fthe切・f・・1∫舳㎝tai・・di・the畑f〃α・Itsh・・1d
benoticedestimationthat

possib1yanyspectral

fai1tostudy

analysisanmethods

accurate

based

spectralonthetraditional

structureofdata．

autoregressive

1．はじめに地球物理学では，地震波形や，地殻変動の日値・時聞値変化など，時系列データを解析対象とすることが多い．こうした時系列データに対する数学的（パラメトリック）モデルとしては，正規性定常確率過程の見本過程と見なせる場合には，自己回帰（AR）モデル，移動平均（MA）モデル，自己回帰移動平均（ARMA）モデルが，非定常過程を含む場合，積分混合モデル（IAR，ARIMA）などがある（BoxandJenkins，1970）．これらのモデルのうち，積分混合モデルについては，地殻変動データの予測などに適用された例がある（大内・高橋，ユ981；橋爪・三雲，

1983）が，最も単純で，しかも扱い易いのは自己回帰モデルで，例えば地震波初動の自動検出（白井・徳弘，1979；浜口・鈴木，i981；森田・浜口，1981；1984），最大エントロピー法に基づくスペクトル解析（Burg，1967；Ur1ych
振動数や，Qの観測（Chao鼻ndandBishop，1975），そして地球自由振動の固有Gilbert，1980；堀・下平・深尾，1984；Horiθfα1、，1989）などに広く適用されている．

一般に，自己回帰モデルの係数を推定する場合，予測誤差を最小にする方法がとられるが，熊沢ら（1983）が提唱した存否法では，非Toep1itz型相関行列の最小固有値に対する固有べ一2一自己回帰モデルの係数の推定方法の比軟一一堀

クトルを用いて係数を推定している．両者の違いは，係数を推定するための過去のデータが，誤差を含んでいると見なすか，そうでないかという点である（大内・奥田，1986）が，その優劣について，詳細な議論はなされていない．本論文では，この2つの係数の推定方法の違いについて，数値実験を中心に解析をし，従来の推定方法の問題点を明らかにしたので報告する．

2．AR係数の推定方法

自己回帰過程は，1〃閉；一∞＜〃＜・・1を，無限に続く時系列信号とするとき，個々の信号間に以下の関係式が成り立つものとして定義される：〃（2−1）〃閉；、≧、w仇一ゴ・

ここで，（巧；ブ＝1，…，〃）は，AR係数であり，〃はAR次数である．この式では，ある時点での値が，それ以前の〃個の値の線形結合で表されることを意味している．熊沢ら（1983）の存否法では（2−1）を，〃（2−2）、≧榊一1＝0・の様に表しているが，（2−2）でαo＝一1とおけば（2−1）と同一の式になり，また，AR係数は相対値にのみ意味があるため，自己回帰過程の記述という意味では，両者に差はない．
ところで，我々が実際に解析対象とするデータ1π1は有限個であり，しかも様々な雑音1θ。1を含んでいる：（2−3）

伽：〃仇十θ仇；〃＝1，…，W．

一般に，12閉1は，信号1〃冊1とは独立であって，E［θ加］＝0，（2−4）
E［θ閉・閉1＝δ蜆榊σ、2，という性質を持つ，正規性白色雑音と仮定される．ただし，δ閉閉はクロネッカーのデルタ，σ、2は雑音の分散，そしてE［・］はアンサンブル平均を表している．

雑音を含んだ有限個の時系列データから，AR係数を推定する場合，AR次数〃を固定して

一3一

国立防災科学披術センター研究報告第42号

1989年2月

おいて，予測誤差∫力を最小にする方法が一般的である：
1〃M∫1一。一肩、、（κ1一，恥一j）2→舳・（2−5）（2−5）は，正規方程式，一÷一（2−6）

リα＝わ，

ただし，1〃

リ＝（｝＝（σ匂）＝（〃＿肌郭r灼一后）1〃＝1・・・ア＝（α1，…，α〃）T，→わ

1N一〃．｝（Σκ仇κr1，…，＝〃十1Σπ〃π閉＿后）T，＝〃十1村に帰着されるが，Nが十分大きいときは左辺の係数行列リの要素Σ。伽＿戸例＿を／（N一〃）は，相蜆＝〃十1関関数ρ〔：E［κ＿〆閉＿后］に置き換えられるとして，（2−7）

肋＝ρ，ただし，月＝（伽）＝（〜）＝（ρゴーセ），

ρ＝（ρ1，…， ρ〃）T，の様に表すことがある．ここで，（2−7）の係数行列Rは，Toep1it・行列になっている．（2−7）はYu1e＿Walker方程式とも呼ばれ，レビンソンーダービンのアルゴリズムによって，逐次的に解くことができる．また，AR係数〃はAlC（Akaike，1973；1977）等によって推定される．

ところで，上述した係数の推定方法（以下，手法1という）は，厳密に云うと，雑音を含む時系列データの中の，完全に自己回帰的な信号部分の線形関係が得られるわけではない．

最小二乗条件（2−5）から得られる1αゴ1は，一4一白己回帰モデルの係数の推定方法の比較一堀〃κ・＝、恥rl＋・1；・＝1・…・w，（2−8）という線形関係を表す係数の不偏推定値になっている．この式では，現在のデータ北蜆は，雑音の現在値2刑とは独立であるが，データの過去値1κ刑＿ゴ；プ＝1，…，〃に含まれる雑音とは完全に独立でない。従って，手法1は過去の（雑音を含んだ）データと，現在のデータ中に存在する信号成分の問の線形関係を推定するためのものであると言える．

これに対して，存否法におけるAR係数は，（雑音を含んだ）データから自己回帰的な信号成分を除去した残差∫をを最小にする様に求められる：ユ」・ 〃∫・＝N一〃 身十、（、恥1r）2→伽・（2−9）

ただし，自明な解αo＝…α〃＝0以外の解を求めるために，AR係数を成分とするベクトル7＝（αo，…，α〃）の大きさは1に固定される：〃Σα2＝1．

（2−10）。一・ゴ

このような条件付き極値問題は，ラグランジュの未定乗数法によって，∂1w〃〃∂、冶〔ト〃斗細・）2一λ財一1）〕一・ll一…

仏（2−11）という方程式を解くことになり，結局，〃、昌（巧1一λδll）ザ0；后一0・…，〃，（2−12）

ただし，1NろFN＿〃蜆伽一戸・一ゲ〃＝0・・〃・という固有値問題に帰着する．行列P＝（㌦）の（〃十1）個の固有値及び固有ベクトルをそれぞれλ（γ），7（γ）（γ＝0，…，〃）とすると，7（γ）の成分を（2−9）に代入して得られる∫、＝∫ほ（γ）は，（2−12）から，λ（γ）に等しくなるので，結局，λ（γ〕のうちで最小の固有値λ（0）に対する固有ベクトルア（0）机AR係数の推定値となる．この手法（以下，手法2という）では，従来の方法とは異なって，連続する（〃十1）個のデータに含まれる信号成分の線形関係を記述する係数の不偏推定値が得られることになる．

一5一

■国立防災科学技術センター研究報告第42号

1989年2月

3．数値実験偶数次の自己回帰モデルで表現される信号は，次数の半数の（増大または減衰する）正弦波の重ね合わせと等価である（F．Oberg，1969）．今，信号〃冊が，L個の正弦波によって構成されているものとする．

止叫＝≧ん・・p（αμ△t）…（ωμ△t＋9）．

（3−1）

ト1

ここで，ん，α，ω1及び，9はそれぞれ，1番目の正弦波の初期振幅，増大率，角振動数及び初期位相を表しており，△fはサンプリング間隔である．この時のAR係数は，αo：1として残りは〃（：2L）元連立一次方程式，〃

止

占≧、、≧ん・・p／ト后）α1△1〕…／σ一后）ω

△叶9、〕α后（3−2）
L

rミん…（プα1△f）…（プω1△f＋91）・プー1，…，〃，を解くことによって厳密に求めることができる．

以下の数値実験では，（3−2）を解いて厳密解として得られたAR係数と，前述した2つの手法で得られたAR係数の推定値を比較することにする．厳密解と推定値の差は，AR係数を（〃十1）次元ベクトルと見なしたときの両者の方向のなす角δθで表すことにする：

Mδθr洲巧／1デ（o）llア1、（3−3）

ただし，7（O）：（α（Oも，…，α（％）及び，ア＝（α。，…，α〃）は，それぞれ厳密解より得られる

AR係数，及び，手法1または手法2から得られるAR係数の推定値である．

3・1

雑音のバワーとの関係

図1aと図1bはそれぞれ，δθと薙音のパワーの関係を直線スケールと対数スケールで示したものである．ここでは，時系列データとして2つの正弦波の重ね合わせを用いており，それぞれの初期振幅，初期位相，角振動数，増大率は，λ1＝λ2：1，？1＝92＝0，ω1＝一α1＝0．1，ω2＝一α2＝O．2とした．また，データ長N及びサンプリング問隔△サはそれぞれ1000，1とした．

そして，このデータに標準偏差σ、の正規性白色雑音を加え，σ、を10■8から1r2まで，増大させていった。図ユ中白丸と黒丸は，それぞれ手法1と手法2の結果である．予想されるように，一6一自已回帰モデルの係数の推定方法の比較

堀m〆OO．
Oo○グOO①ω旧Oo●●OO●OO●。OOd．・一○o．●●．．●．●．．○●●一．．．｝10■810一？

10−5

10−6

10，410−3

lO−2σ。

（・）

Oo88●．

。

Oo●．・○．・1．
■

O8。！

①ωO。O

O．．○● ●●■、●

●

・1
■●．
●

もF，408「O■8107101010一10■10−2σ。

（b）

図1．

雑音のパワーの増大に伴うAR係数の推定値の変化を直線スケール（a）と対数スケール（b）で示したもの．白丸と黒丸はそれぞれ手法1と手法2の結果である．σ、：正規性
白色雑音の標準偏一差．δθ：AR係数の推定値及び真の値からなる2つのベクトルの方向の差（ラジアン）．

Fig．1

Variation

ofof

theAR

shown

the1inearsca1e（a）andin

the1ogoflandmethod2respective1y。σ、：Standard

deviationofGaussian

whitenoise．δθ：Differenceindirections

and

true

vectorconsti−inresults

tutedbyARthemethodestimate

ofcoefficients

with副nsca1e（b）一〇pen

theestimatedcoefficients（radian）一一7一

andvectorincreaseoffi11edcirc1es

the

the

noise

power

representthe国立防災科学披術センター研究報告第42号

1989年2月

雑音のパワーが小さい場合（σ、＜10■6），両者ともδθ＜10山2で良好な結果が得られている．

しかし，手法1の結果は，直線スケールで見るとσ、＝2×1016からδθが急激に増大し，σ、：10−5では約0．5（約3ぴ）にも達しており，真の値からかなりずれている．これに対し，手法2の結果は，対数スケールで見ると，ほぼδθが直線的に増大しているカ㍉σ、＝ユO■3でもδ

θは約O．1（約6。）と十分小さい．以上から，手法1で得られるAR係数の推定値は，手法2によるものと比べると，雑音の混入に対してはるかに敏感であることが分かる．

3．2

データ長との関係

次に，雑音のパワーを固定しておいて，データ長〃を変えたときのAR係数の推定値の変化を調べる．図2aと図2bはWを101から105まで増大させたときのδθの値を，それぞれ直線スケールと対数スケールで示したものである．ここで用いたデータも2つの正弦波と正規性白色雑音の重ね合わせであるが，データ長の影響だけを見るため，増大率は2つとも0にして，全データを通して∫／N比を一定にした．振幅，位相など他のパラメータは3．1と同一である．
また雑音の標準偏差σ、は10■4に固定した．

データ長が有限の場合，雑音は程度の差こそあれ，「有色」となる．しかし，データ長が大きくなれば（2−4）で表される「正規性」と「白色性」の傾向が強くなってくると考えられる．従って，減衰しない信号ならば，データを長くとればとるほど雑音と分離されてAR係数の推定値については良好な結果が得られるはずである．ところが図2をみると，手法1ではデー
タ長が200ポイントになるまではδθは減少していくが，それ以上はいくらデータを増やしていってもδθはほぼ一定で，よい結果は得られてこない．これに対して手法2では，直線スケー

ルで見ると，データ長を100ポイントまで増やした時点でδθの値はほぼ0となって，真の値に近いAR係数を推定していることが分かる．また，対数スケールで見ると，手法1の結果は，データ長を大きくしても全く向上していないのに対して，手法2の結果では，多少のばらつきはあるものの，データ長を増やした分，結果が向上するという，wとδθの間の一貫した負の相関が認められる．

一8一自己回帰モデルの係数の推定方法の比較

ω堀●o●

O

OOO oCO・。。。O帖

O●● ●●●●m8O

OOOO00①OOOωO●●、o90●●●｝103l0210●

●●●

10105

N（・）

O

O Oα刃0D OO●OO◎Oα〕Q＝□〕OO

OOOOCl＝0O

OO

●●●

■
■、●

o■

●

●

Φω

．●、●

●

O

●

●
■

●●．、■

O

●

●●
●．

103

l02

10l04

105

N（b）

図2．

データ長の増大に伴うAR係数の推定値の変化．N：計算に用いたデータのポイント数、

Fig．2

Variationof

theestimateshownintheresu1tsof

method1mdpoints

thelinearoftheAR

coefficientsscale（a）andin

thewithlog

anmethod2respectively、〃：Numbercu1ation、

一9一

increase

thenumberofandfi11edcircles

represent

data

pointsusedscale．0pen

ofof

in

data

cal−

■国立防災科学技術センター研究報告

3．3

第42号

1989年2月

分解能

最後の雑音パワーとデータ長を一定にしたときの分解能について，手法1と手法2を比較してみた．時系列データとしては，2つの正弦波と正規性白色雑音を合成したものを用いた．振

幅，位相，増大率は2つとも同一で，それぞれ，ユ，0，0とした．また，データ長は10000
ポイント，サンプリング間隔は1，雑音の標準偏差は10■3にした．そして，角振動数はω、＝0．25／2π（ナイキスト周波数の1／2），ω2＝ω1＋△ω／2πとして，△ωとδθの関係を調べた．その結果を図3に示す．周波数の接近した2つの正弦波を合成すると，「うなり」が発生し，その周期は，2つの周波数の差の逆数になることが知られているが，今回の実験に用いたデータ長は，△ω／2π＝10■4でちょうどうなり1周期分に相当することになる．図3を見るとここでも手法2の方が手法1よりも優れていることが明らかになっているが，特に，手法2の方は△ω／2π＝10■5，つまり，うなりの周期が全データ長の10倍にも達する時でさえ，δθが10−2まで減少しているのは注目に値する．

一1

一■一1一「■「一τ1［88888㎜888・％
●．

0●．
．．

O

■

8O ％■

■

OO
●OO

■o O

●

（〜）

○●

■

●

O

●
●

○●

、．

lO

図3．

7

10一lO1
△ω／2兀

1O．4

10

2つの正弦波の振動数の差の増大に伴うAR係数の推定値の変化．△ω：角振動数の差（サイクル毎サンプリング間隔）

Fig・3V・・i・ti㎝・fth…tim・t・・fth・AR…ffi・i・・t・with・・i・・・・・…fth・f・・q・…ydiff。。一
enceofthetwosinusoids．Openandfi11edcircles

represent

theresu1tsof

method1andm・th・d2…p・・ti・・ly・△ω：F・・q・…ydiff・・・…b・tw…tw・w・…1・m・・t・（・y・1・p・・

one

samp1ing

interva1）．

一10一自己回帰モデルの係数の推定方法の比較一堀

4．議論と結論
今回行った単純な3つの数値実験では，AR係数を時系列データから推定する場合，いずれも，従来の予測誤差を最小にする方法より，信号を除去した残差を最小にする存否法の推定方法の方が優れていることが分かった．従来の方法が，雑音の混入に対して敏感で，データ長を大きくとることの利点が生かされにくいのは，データ中に含まれる信号成分と雑音を完全に独立なものとして扱っていないためであることは明らかである．つまり，自己回帰モデルを記述する場合，本来分けて記述すべき（2−1）ないし（2−2）と（2−3）を，暗黙のうちにユつにまとめて（2−8）で表現することの，本質的な問題点が今回の研究で判明した．

近年，自己回帰モデルはスペクトル解析法にも適用されるようになってきて，例えばChaoand

Gi1bert（1980）は，地球の自由振動の固有周波数を複素領域で求める手法を提唱している．

同じように自己回帰モデルを自由振動の解析に適用する手法は，例えば，存否法を応用したものとして堀・下平・深尾（1984）やHoriθ川．（1989）も考案している．地球自由振動の場合，そのスペクトル構造は比較的単純で，自己回帰モデルを適用するのにふさわしいデータであるといえる．しかし，Chaoand

Gi1bert（ユ980）のAR係数の推定方法は基本的には手法1に基づいており，地震記象のような，かなり雑音を含んだ時系列データに対する効果には疑問が持たれる．さらに，自由振動を解析する場合，時には，1週問を超えるようなデータを用いる場

合があるが，このようなデータに対しては，長いデータほど雑音と信号の分離がうまくいく手
法2に基づくスペクトル解析法，つまり存否法の方がより精度の高い結果が得られると考えられる．

今回明らかにした，従来のAR係数の推定方法における問題点は，予測フィルターや，ARMAモデルの推定にも関係があると思われるが，それについては今後さらに検討すべきであろう．

参考文献
1）Ak・i…H・（・973）11・f・…ti・・t・・・・…d・…t…i…ft・・…im・li・。1i・。。d。。i．ci．1、．2〃サ、ザ舳㎝α1∫榊・∫｛・閉㎝〃・舳αf1㎝τんθ・び，・d・・B・N−P・t・・…dF．C・・ki，PP．267−281．Ak．d．mi．Ki．d。，

Budapest．
2）Ak・ik・・H・（1977）1O…t・・・・…i・i・・ti・…i・・ip1・、P伽・伽1パ〃。∫1舳㎝肋、τ、、d、閉，N．r。．

Hol1and，27−41．

3）B…G・E・…1・・ki…G・M・（・970）1・i・・…i・・…1・・i・．F・舳物α舳舳。1，H．1d。。Da。，Sa，Francisco．
4）Burg，J・P・（1967）：Maximum

entropyspectra1ana1ysis．Pα枇γ仰召∫6舳〃、㍑加37肋λ舳〃αけ〃2㎜、κ㎝α1〃2〃伽g，Soc．Explor．Geophysics，Oklahoma

City

Ok1a．，Oct．31．

一11一

国立防災科学技術センター研究報告
5）Chao，B．F．and
quencyseismic

Gi1bert，F、（1980）1Autoregressive

第42号

estimationof

ユ989年2月

comp1exeigenfrequencies

inlow

fre−

spectra．G20肋ツ3I∫R．λs伽㎝、∫oc。，63，641−657．

6）FrOberg，C．E．（1969）：Introduction

toNumerical

Ana1ysis，2ndedn一，pp．347−349，Addis㎝一Wesley，Reading，Massachusetts．

7）浜口博之・鈴木次郎（1979）：微小地震P波初動の自動検測とその信頼性について．自然災害特別研
究成果，A−54−2，62−83．

8）橋爪道郎・三雲健（1983）：地殻変動観則データの多変量解析一AR1MAモデルによる予測の試み一．

地震2，36，515−530．

9）堀貞喜・下平

1984，

No1，

慶・深尾良夫（1984）：存否法によるマルチプレットの解析．地震学会講演予稿集，

B59．

10）Hori，S．，Fukao，Y．，Kumazawa，M．，Furumoto，M．，Yamamoto，A．（1989）：A
analysis

andits

app1ication

tothe

Earthsfreenew

methodofosciuations：the∫㎝一切method。∫．Gω助ツ∫。地s，in11）熊沢峰夫・古本宗充・深尾良夫・山本明彦・水谷

spectralpress．

仁（1983）：存否スペクトル法一原理と特長，地震学会講演予稿集，1983，No．1，B33．

12）森田裕一・浜口博之（1981）：2次元自己回帰過程による∫波初動の自動検測．地震2，34，223−240．

13）森田裕一・浜口博之（1984）：自己回帰過程による地震波初動の自動検測とその信頼区間．地震2，37，281−293．

14）大内徹・奥田暁（1986）：地震観測データの統計解析（1）一最小二乗法と主成分分析法一．地震学会講演予稿集，1986，No．2，A95．

15）大内徹・高橋享（1981）：地震学で扱う非定常時系列データの解析一予測と異常検出の試み一、
地震2，34，455−464．
16）白井克彦．徳弘一路（1979）：地震波初動の自動検出．地震2，32，・141−148．
17）Ur1ych，T．J．andBishop，T．N．（1975）：Maximum

entropyspectra1ana1ysisand

autoregressive

decom−

position肋．Gω伽s．∫似2、物s．，13，183−200．

（1988年6月20日原稿受理）

一12一